Stabilized energy factorization approach for Allen-Cahn equation with logarithmic Flory-Huggins potential

TL;DR

This paper introduces a stabilized energy factorization method for the Allen-Cahn equation with logarithmic Flory-Huggins potential, ensuring linearity, stability, and ease of implementation while preserving key physical properties.

Contribution

The paper develops a novel stabilized energy factorization scheme that handles the nonlinear Flory-Huggins potential efficiently and guarantees maximum principle and energy stability.

Findings

The scheme is unconditionally energy stable.

All nonlinear terms are treated semi-implicitly.

Numerical results confirm stability and effectiveness.

Abstract

The Allen--Cahn equation is one of fundamental equations of phase-field models, while the logarithmic Flory--Huggins potential is one of the most useful energy potentials in various phase-field models. In this paper, we consider numerical schemes for solving the Allen--Cahn equation with logarithmic Flory--Huggins potential. The main challenge is how to design efficient numerical schemes that preserve the maximum principle and energy dissipation law due to the strong nonlinearity of the energy potential function. We propose a novel energy factorization approach with the stability technique, which is called stabilized energy factorization approach, to deal with the Flory--Huggins potential. One advantage of the proposed approach is that all nonlinear terms can be treated semi-implicitly and the resultant numerical scheme is purely linear and easy to implement. Moreover, the discrete maximum principle and unconditional energy stability of the proposed scheme are rigorously proved using the discrete variational principle. Numerical results are presented to demonstrate the stability and effectiveness of the proposed scheme.

Equations

$$\frac{\partial \varphi (\mathbf{x},t)}{\partial t}-{ϵ}^2\Delta \varphi (\mathbf{x},t)+f(\varphi (\mathbf{x},t))=0,\mathbf{x}\in \Omega ,$$

$$F(\varphi )=\varphi \ln (\varphi )+(1-\varphi )\ln (1-\varphi )+\theta \left(\varphi -{\varphi }^2\right),0<\varphi <1,$$

$$E(\varphi )={\int }_{\Omega }\left(F(\varphi )+\frac{1}{2}{ϵ}^2|\nabla \varphi {|}^2\right)d\mathbf{x},$$

$$\frac{\partial }{\partial t}E(\varphi )=-{\int }_{\Omega }{\left|\frac{\partial \varphi }{\partial t}\right|}^{2}d\mathbf{x}\le 0.$$

$$\mathcal{C}_{h}=\big{\{}\phi:(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}})\mapsto\mathbb{R},~{}~{}0\leq i\leq N-1,~{}0\leq j\leq M-1\big{\}},$$

$$\mathcal{U}_{h}=\big{\{}u:(x_{i},y_{j+\frac{1}{2}})\mapsto\mathbb{R},~{}~{}0\leq i\leq N,~{}0\leq j\leq M-1\big{\}},$$

$$\mathcal{V}_{h}=\big{\{}v:(x_{i+\frac{1}{2}},y_{j})\mapsto\mathbb{R},~{}~{}0\leq i\leq N-1,~{}0\leq j\leq M\big{\}}.$$

$$\mathcal{U}^{0}_{h}=\big{\{}u\in\mathcal{U}_{h}|~{}u_{0,j+\frac{1}{2}}=u_{N,j+\frac{1}{2}}=0,~{}~{}0\leq j\leq M-1\big{\}},$$

$$\mathcal{V}^{0}_{h}=\big{\{}v\in\mathcal{V}_{h}|~{}v_{i+\frac{1}{2},0}=v_{i+\frac{1}{2},M}=0,~{}~{}0\leq i\leq N-1\big{\}}.$$

$${\nabla }_{h,x}{\varphi }_{i,j+\frac{1}{2}}=\frac{{\varphi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\varphi }_{i-\frac{1}{2},j+\frac{1}{2}}}{h},1\le i\le N-1,0\le j\le M-1,$$

$${\nabla }_{h,y}{\varphi }_{i+\frac{1}{2},j}=\frac{{\varphi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\varphi }_{i+\frac{1}{2},j-\frac{1}{2}}}{h},0\le i\le N-1,1\le j\le M-1.$$

$$D_x{u}_{i+\frac{1}{2},j+\frac{1}{2}}=\frac{u_{i+1,j+\frac{1}{2}}-u_{i,j+\frac{1}{2}}}{h},D_y{v}_{i+\frac{1}{2},j+\frac{1}{2}}=\frac{v_{i+\frac{1}{2},j+1}-v_{i+\frac{1}{2},j}}{h},$$

$${\left(\varphi ,{\varphi }^{\prime }\right)}_h=h^2\sum _{i=0}^{N-1}\sum _{j=0}^{M-1}{\varphi }_{i+\frac{1}{2},j+\frac{1}{2}}{\varphi }_{i+\frac{1}{2},j+\frac{1}{2}}^{\prime },\varphi ,{\varphi }^{\prime }\in {\mathcal{C}}_h,$$

$${\left(u,u^{\prime }\right)}_h=h^2\sum _{i=1}^{N-1}\sum _{j=0}^{M-1}{u}_{i,j+\frac{1}{2}}{u}_{i,j+\frac{1}{2}}^{\prime },u,u^{\prime }\in {\mathcal{U}}_h^0,$$

$${\left(v,v^{\prime }\right)}_h=h^2\sum _{i=0}^{N-1}\sum _{j=1}^{M-1}{v}_{i+\frac{1}{2},j}{v}_{i+\frac{1}{2},j}^{\prime },v,v^{\prime }\in {\mathcal{V}}_h^0.$$

$$\Vert \varphi {\Vert }_h={\left(\varphi ,\varphi \right)}_h^{1/2},\Vert u{\Vert }_h={\left(u,u\right)}_h^{1/2},\Vert v{\Vert }_h={\left(v,v\right)}_h^{1/2},$$

$$\Vert {\nabla }_h\varphi {\Vert }_h^2=\Vert {\nabla }_{h,x}\varphi {\Vert }_h^2+\Vert {\nabla }_{h,y}\varphi {\Vert }_h^2.$$

$${\left(u,{\nabla }_{h,x}\varphi \right)}_h=-{\left(D_{x}u,\varphi \right)}_h,u\in {\mathcal{U}}_h^0,\varphi \in {\mathcal{C}}_h,$$

$${\left(v,{\nabla }_{h,y}\varphi \right)}_h=-{\left(D_{y}v,\varphi \right)}_h,v\in {\mathcal{V}}_h^0,\varphi \in {\mathcal{C}}_h,$$

$$-{\left({\Delta }_h\varphi ,\phi \right)}_h={\left({\nabla }_h\varphi ,{\nabla }_h\phi \right)}_h,\varphi ,\phi \in {\mathcal{C}}_h.$$

$$\Vert {\nabla }_h{\varphi }_{-}{\Vert }_h^2\le -{\left({\Delta }_h\varphi ,{\varphi }_{-}\right)}_h,$$

$$\Vert {\nabla }_h{\varphi }_{+}{\Vert }_h^2\le -{\left({\Delta }_h\varphi ,{\varphi }_{+}\right)}_h.$$

$$E(\varphi )={\int }_{\Omega }\sum _{i=1}^M{\Phi }_i(\varphi )\cdot {\Psi }_i(\varphi )d\mathbf{x},$$

$$\begin{gathered} E({\varphi }^{n+1})-E({\varphi }^n)\le {\int }_{\Omega }\mu \left({\varphi }^{n+1},{\varphi }^n,\cdots ,{\varphi }^0) \\ ({\varphi }^{n+1}-{\varphi }^n\right)d\mathbf{x}. \end{gathered}$$

$$E_{\nabla }(\varphi )=\frac{1}{2}{\int }_{\Omega }{ϵ}^2|\nabla \varphi {|}^{2}d\mathbf{x},$$

$$E_{\nabla }({\varphi }^{n+1})-E_{\nabla }({\varphi }^n)$$

$$={\int }_{\Omega }{ϵ}^2\nabla {\varphi }^{n+1}\cdot \nabla \left({\varphi }^{n+1}-{\varphi }^n\right)d\mathbf{x}-\frac{1}{2}{\int }_{\Omega }{ϵ}^2|\nabla ({\varphi }^{n+1}-{\varphi }^n){|}^{2}d\mathbf{x}$$

$$\le {\int }_{\Omega }{ϵ}^2\nabla {\varphi }^{n+1}\cdot \nabla \left({\varphi }^{n+1}-{\varphi }^n\right)d\mathbf{x}$$

$$=-{\int }_{\Omega }\left({\varphi }^{n+1}-{\varphi }^n\right){ϵ}^2\Delta {\varphi }^{n+1}d\mathbf{x},$$

$$E_{\nabla }({\varphi }^{n+1})-E_{\nabla }({\varphi }^n)$$

$$=-\frac{1}{2}{\int }_{\Omega }\left({\varphi }^{n+1}-{\varphi }^n\right){ϵ}^2\Delta \left({\varphi }^{n+1}+{\varphi }^n\right)d\mathbf{x},$$

$$H_a(\varphi )=\varphi \ln \left(\varphi \right),H_b(\varphi )=\left(1-\varphi \right)\ln \left(1-\varphi \right),H_c(\varphi )=\theta \varphi \left(1-\varphi \right),$$

$$F(\varphi )=\left(\lambda +1\right)\left(H_a(\varphi )+H_b(\varphi )\right)+H_c(\varphi )-\lambda \left(H_a(\varphi )+H_b(\varphi )\right),$$

Full text

1 Introduction

It is known that the logarithmic Flory–Huggins potential [5, 55, 49, 8, 19, 10, 20, 23, 52, 49] and double well potential [24, 26, 11, 60, 37, 21, 7, 38, 53, 42] are the most useful energy potentials in phase field models. Moreover, the Flory–Huggins energy potential is more realistic than the double well potential from the viewpoint of physics[5, 8]. However, compared to the double well potential, there are more scarce works of numerical algorithm developments and numerical analysis devoted to the Flory–Huggins potential [55]. There are two possible reasons: one is the difficulty resulting from the strong nonlinearity of this potential function, and the other is the singularity in the situations that the phase variable approaches the limits 0 and 1. So it becomes a very challenging issue to design efficient algorithms for phase-field models with the Flory–Huggins potential. The Allen–Cahn equation, as well as the Cahn–Hilliard equation [6, 5], is the fundamental equation of phase-field models, so it has become an attractive research topic in mathematical analysis [13] and numerical simulation [16, 17]. The maximum principle and energy dissipation law are two intrinsic properties of the Allen–Cahn equation [13]. Numerical schemes preserving these key properties are highly preferred for solving this equation. However, it is still challenging to design numerical schemes that mimic these properties. A brief review on the developments of such schemes in the literature will be stated as below. The main goal of this work is to propose linear, unconditionally energy stable and maximum principle preserving numerical schemes for the Allen–Cahn equation with the Flory–Huggins potential.

The traditional explicit and fully implicit time marching schemes for bulk free energy functions have been shown to suffer from very severe constraints on the time step size [4, 16, 10]. Therefore, it becomes a major issue to design the semi-implicit energy stable numerical schemes that remove the constraints on the time step size. A popular approach used for designing such schemes is the convex splitting approach [12, 14, 43, 51, 3, 11, 60, 37, 21, 44, 36, 59], which produces unconditionally energy stable schemes by the use of implicit treatment for the convex terms of the energy functions and explicit treatment for the concave terms. Besides various phase-field models, it has also been successfully applied to deal with the Helmholtz free energy based on a realistic equation of state in recent years [41, 15, 32, 29, 40, 33, 34]. However, it usually results in the nonlinear schemes; for instance, the convex splitting approach for the logarithmic Flory–Huggins potential leads to the scheme involving the implicit logarithmic function. As a consequence, practical implementation of such schemes requires nonlinear iterative solvers and computational cost may be expensive accordingly. The other commonly used conventional approach for constructing linear numerical schemes is the linear stabilization approach [7, 38, 53, 54, 42, 26, 27], which simply treats all nonlinear terms by the fully explicit way and introduces a linear stabilization term to remove the time step constraint. We observe that the stabilization approach is effective for the double well potential, but it works not well for the logarithmic potential probably because of more complicate nonlinearity of it. In [28, 32], a nonlinear stabilization approach has been proposed for the Peng-Robinson equation of state (PR-EOS) [39] that is one of the most useful tools in petroleum industry and chemical engineering. For the logarithmic Flory–Huggins potential, the instability will occur and numerical results may be out of normal range when the energy parameter takes a large value. Inspired by the approach in [28, 32], we will incorporate the nonlinear stabilization term for the logarithmic Flory–Huggins potential to ensure the symmetric positive definiteness and discrete maximum principle of the resultant scheme in the case with a large energy parameter.

In recent years, the novel auxiliary variable approaches have been developed and successfully applied to design linear numerical schemes for various diffuse interface models [55, 56, 57, 61, 35, 30, 58, 45, 46, 62, 63]. The first approach is the so-called invariant energy quadratization (IEQ) approach [55, 56, 57, 61] that has been successfully applied to devise efficient, linear schemes for various phase-field models intensively in recent years. The basic idea of IEQ is to define a set of auxiliary variables and then transform the original free energy function into a quadratic form. The second approach is the scalar auxiliary variable (SAV) approach proposed in [45], which differs from IEQ that it uses a scalar auxiliary variable instead of the space-dependent auxiliary variable. The convergence and error estimates for the SAV schemes are studied in [46]. Very recently, a generalized positive auxiliary variable, termed gPAV, is proposed in [58] inspired by SAV and IEQ. Numerical schemes developed by the auxiliary variable approaches are linear and easy to implement. As a result, such approaches have rapidly become the useful and successful tools for simulating a variety of diffuse interface models. The auxiliary variable approaches use the transformed energies that are equivalent to the original energies at the continuous level, but there are the discrete errors between the original energies and transformed energies at the time-discrete level [55], thus the resultant schemes may not preserve the original energy dissipation law although the transformed energy dissipation can be proved.

More recently, a novel energy factorization (EF) approach is proposed in [31] to design linear, efficient numerical schemes for the diffuse interface model with PR-EOS. The basic idea of EF is to factorize an energy function/term into a product of several factors, which can be treated in the energy difference by the use of individual properties of each factor. Compared with the convex splitting approach, the EF approach can produce the linear semi-implicit schemes. It is different from the IEQ/SAV approach that the EF approach never introduces any new independent energy variable, and thus the resultant schemes can preserve the original energy dissipation law. The EF approach has been successfully applied to design the semi-implicit linear schemes for the PR-EOS model [31]. In this paper, we will propose the stabilized EF approach to deal with the logarithmic Flory–Huggins energy potential, and from this, we will propose a linear semi-implicit discrete scheme inheriting the original energy dissipation law. The proposed scheme is efficient and very easy to implement.

For the phase-field models, the phase variables usually comply with the specific maximum principles in terms of their physical meanings. It is known that the Allen-Cahn equation satisfies the maximum principle [13]. A discrete scheme preserving the maximum principle is highly preferred for solving such equations since it can eliminate spurious numerical solutions such that it can not only ensure the physical reasonability of numerical results but also improve the long-time stability of numerical simulation substantially. However, the efforts regarding the maximum principle preserving discrete schemes and numerical analysis are even more scarce partially because such schemes are actually very challenging due to particularities of involved spatial and temporal discretization [8]. In [47], the commonly used numerical method for the Allen-Cahn equation with the double well potential was proved to preserve the maximum principle by the matrix analysis approach. For the logarithmic Flory–Huggins energy potential, the implicit Euler scheme was analyzed in [10] under the restricted constraint on the time step size due to the fully implicit treatment for all energy terms, and recently, the discrete maximum principle of the semi-implicit convex splitting scheme was analyzed in [8]. More recently, the maximum principle of a linear numerical scheme for the PR-EOS diffuse interface model has been proved in [31]. However, there are no results regarding the linear semi-implicit scheme with the maximum principle for the phase-field models with the logarithmic Flory–Huggins energy potential. This gap will be filled in this paper. The discrete maximum principle of the proposed scheme will be rigorously proved using the discrete variational principle [31]. Appropriate stability strategy plays a crucial role in preserving the discrete maximum principle.

The cell-centered finite difference (CCFD) method [48] is applied as the spatial discretization method. The CCFD method can be equivalent to a special mixed finite element method with quadrature rules [2] and has been applied for the phase field models [18, 22, 25, 51, 50, 9]. The discrete variational principle will be used to carry out theoretical analysis of the proposed scheme.

The new aspects of the current work are listed as follows: (1) the stabilized energy factorization approach to treat the logarithmic Flory–Huggins potential; (2) the purely linear numerical scheme inheriting the original energy dissipation law; (3) the discrete maximum principle of the proposed scheme.

The rest of this paper is structured as follows. In Section 2, we describe the Allen–Cahn equation with Flory–Huggins potential. In Section 3, we introduce the discrete function spaces and discrete operators based on CCFD. In Sections 4, we propose the stabilized energy factorization approach to deal with the logarithmic Flory–Huggins potential, and then present the fully discrete scheme. In Section 5, we carry out theoretical analysis, including the well posedness of the discrete solutions, discrete maximum principle and unconditional energy stability. In Section 6, numerical results are presented to demonstrate the effectiveness and stability of the proposed scheme. Finally, the concluding remarks are given in Section 7.

2 Allen–Cahn equation with Flory–Huggins potential

This paper is concerned with the development of efficient numerical schemes for the following Allen–Cahn equation [1]

[TABLE]

where $t$ is the time, $\Omega$ denotes the bounded domain with smooth boundaries and $\epsilon$ is a positive constant measuring the interfacial width. Here, $\phi$ represents the phase variable and $f(\phi)=F^{\prime}(\phi)$ is the chemical potential, where $F(\phi)$ is the Helmholtz free energy. Meanwhile, the equation (2.1) is subject to specified initial condition and homogeneous Neumann/Dirichlet or periodic boundary conditions.

In this paper, we consider the logarithmic Flory–Huggins energy potential [5, 55, 49]

[TABLE]

where $\theta>2$ is the energy parameter. It is noticed in [49] and also shown in Figure 1 that the choice of $\theta>2$ is necessary for $F(\phi)$ since in this case it has two wells and admits two phases; otherwise, for $\theta\leq 2$, only a single well exists as well as a single phase. It is clear that the maximum principle, i.e. $0<\phi<1$, is crucial to ensure that $F(\phi)$ and $f(\phi)$ are well defined in both mathematics and physics.

The equation (2.1) follows the energy dissipation law; in fact, we define the energy function

[TABLE]

then $E(\phi)$ is decreasing with time as follows

[TABLE]

3 Discrete function spaces and discrete operators

In this section, we describe the discrete function spaces, discrete operators and discrete variational principle [48, 18, 22, 25, 51, 50, 9, 31, 30], which are based on the cell-centered finite difference (CCFD) method [48, 2]. Here, we consider the two-dimensional case only, but the forms of the three-dimensional case are similar. Let $\Omega=[l_{x}^{a},l_{x}^{b}]\times[l_{y}^{a},l_{y}^{b}]$, where $l_{x}^{a}<l_{x}^{b}$ and $l_{y}^{a}<l_{y}^{b}$. For simplicity, the uniform division is used and the mesh size is denoted by $h=\frac{l_{x}^{b}-l_{x}^{a}}{N}=\frac{l_{y}^{b}-l_{y}^{a}}{M}$, where $N$ and $M$ are positive integers. The grid points are denoted by $x_{i}=l_{x}^{a}+ih$ and $y_{j}=l_{y}^{a}+jh$, and furthermore, the intermediate points are denoted by $x_{i+\frac{1}{2}}=x_{i}+\frac{1}{2}h$ and $y_{j+\frac{1}{2}}=y_{j}+\frac{1}{2}h$.

The discrete function spaces are defined as follows

[TABLE]

[TABLE]

[TABLE]

Here, we use the identification $\phi_{i+\frac{1}{2},j+\frac{1}{2}}=\phi(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}})$, etc. The homogeneous Neumann boundary condition is considered in this paper, and the case with the homogeneous Dirichlet or periodic boundary condition can be formulated and analyzed analogously. We introduce the subsets of $\mathcal{V}_{u}$ and $\mathcal{V}_{v}$ involving the homogeneous Neumann boundary condition as

[TABLE]

[TABLE]

For $\phi\in\mathcal{C}_{h}$, the discrete gradient operator is defined as $\nabla_{h}=[\nabla_{h,x},\nabla_{h,y}]^{T}$, where $\nabla_{h,x}\phi\in\mathcal{U}_{h}^{0}$ and $\nabla_{h,y}\phi\in\mathcal{V}_{h}^{0}$ are defined as follows

[TABLE]

For $u\in\mathcal{U}_{h}$ and $v\in\mathcal{V}_{h}$, the discrete divergence operator is defined as $\nabla_{h}\cdot[u,v]^{T}=D_{x}u+D_{y}v$, where $D_{x}$ and $D_{y}$ are defined as follows

[TABLE]

where $0\leq i\leq N-1,~{}0\leq j\leq M-1$. The discrete Laplace operator is defined as $\Delta_{h}=\nabla_{h}\cdot\nabla_{h}$.

We define the following discrete inner-products:

[TABLE]

[TABLE]

[TABLE]

The discrete norms for $\phi\in\mathcal{C}_{h}$, $u\in\mathcal{U}_{h}^{0}$ and $v\in\mathcal{V}_{h}^{0}$ are denoted as

[TABLE]

[TABLE]

The following discrete variational formulas are obtained by direct calculations [51, 50, 9, 31, 30]

[TABLE]

[TABLE]

[TABLE]

The following lemma is a direct consequence of Lemma 5.1 in [31].

Lemma 3.1**.**

Let $\phi_{-}=\min(\phi-a,0)$ and $\phi_{+}=\max(\phi-b,0)$, where $\phi\in\mathcal{C}_{h}$ and $a<b$. Then we have

[TABLE]

[TABLE]

We note that for the homogeneous Dirichlet and periodic boundary conditions, we have the similar results to (3.8) and (3.9).

4 Stabilized energy factorization approach and discrete scheme

Let $\phi^{n}$ denote the discrete function at the time level $t=t_{n}$, where $n\geq 0$. The key idea of the energy factorization (EF) approach is that we first factorize the energy function as follows

[TABLE]

where $\Phi_{i}(\phi)$ and $\Psi_{i}(\phi)$ are the energy factors, $M\geq 1$, and then we derive the following energy inequality using the properties of the factors

[TABLE]

Thus, $\mu$ is the discrete general chemical potential, which is usually a functional of $\phi^{n+1}$ and $\phi^{n}$, and more generally may rely on $\phi^{n-1},\cdots,\phi^{0}$ as well. For a specific energy $E(\phi)$, there may be several different factorization forms, but not all of the resultant discrete chemical potentials $\mu$ are linear with respect to $\phi^{n+1}$. So an ingenious factorization approach is required for the sake of obtaining a linear efficient scheme.

As a simple example, we consider the gradient free energy

[TABLE]

and we can derive that

[TABLE]

which leads to the classical implicit chemical potential $\mu_{\nabla}^{n+1}=-\epsilon^{2}\Delta\phi^{n+1}$. The other treatment for $E_{\nabla}(\phi)$ is described as follows

[TABLE]

which yields the classical semi-implicit chemical potential $\mu_{\nabla}^{n+1}=-\frac{1}{2}\epsilon^{2}\Delta\left(\phi^{n+1}+\phi^{n}\right)$.

4.1 Stabilized energy factorization approach for Flory–Huggins potential

In order to treat the logarithmic Flory–Huggins potential, we introduce a stabilized energy factorization approach, which combines the stability technique with the energy factorization approach proposed in [31] for the PR-EOS based free energy. To apply the energy factorization approach, we define

[TABLE]

and from this, we rewrite $F(\phi)$ as

[TABLE]

where $\lambda\geq 0$ is a constant. The stability term is crucial to ensure the symmetric positive definiteness and discrete maximum principle of the resultant scheme in the case with a large $\theta$, which will be demonstrated in Subsection 5.1 and Subsection 5.2 respectively.

It is clearly observed that $H_{a}(\phi)$ as well as $H_{b}(\phi)$ can be factorized into the product of a linear function and a logarithm function, while $H_{c}(\phi)$ is a product of two linear functions. For $\phi^{n}>0$ and $\phi^{n+1}>0$, we can deduce that

[TABLE]

where the last inequality is obtained using the concavity of $\ln(\phi)$ as

[TABLE]

Noting that $\ln(1-\phi)$ is also concave, we have

[TABLE]

Thus, for $\phi^{n}<1$ and $\phi^{n+1}<1$, we deduce that

[TABLE]

For $H_{c}(\phi)$, we apply the factorization approach to deal with $H_{c}(\phi)$ as follows

[TABLE]

Due to the concavity of $H_{a}(\phi)$, we have $\ln(\phi^{n})-\ln(\phi^{n+1})\leq\frac{1}{\phi^{n+1}}\left(\phi^{n}-\phi^{n+1}\right)$, thus we deduce

[TABLE]

It is similar to deduce that

[TABLE]

Using (4.8)-(4.12), we define the discrete chemical potential as

[TABLE]

From (4.7)-(4.13), we obtain the following energy inequality

[TABLE]

There is the other alternative approach for $H_{c}(\phi)$ on the basis of the concavity of $H_{c}(\phi)$, which is expressed as

[TABLE]

In this case, the chemical potential $f(\phi^{n},\phi^{n+1})$ is defined as below

[TABLE]

It is observed from numerical tests that the treatment given in (4.10) performs better than that of (4.15).

4.2 Fully discrete scheme

The semi-implicit discrete scheme for the Allen-Cahn equation with the logarithmic Flory–Huggins potential is stated as: given $\phi^{n}\in\mathcal{C}_{h}$, find $\phi^{n+1}\in\mathcal{C}_{h}$ such that

[TABLE]

where $\tau$ denotes the time step size, $f(\phi^{n},\phi^{n+1})$ is defined in (4.13) or (4.16) and $\phi^{0}$ is provided by the initial condition.

In what follows, we will consider only $f(\phi^{n},\phi^{n+1})$ defined in (4.13) in theoretical analysis, but the scheme with (4.16) can be analyzed in a very similar routine.

For convenience of theoretical analysis, using (4.13), we rewrite (4.17) as the following equivalent from

[TABLE]

where $\nu_{\theta,\lambda}(\phi)$ and $r_{\theta,\lambda}(\phi)$ are defined as follows

[TABLE]

[TABLE]

Remark 4.1**.**

Since the linear discrete chemical potentials are obtained using the energy factorization approach, the proposed discrete scheme is linear, easy to implement, and preserve the original energy dissipation law. However, the commonly used convex splitting approach for the logarithmic Flory–Huggins potential results in the nonlinear discrete chemical potential as well as the nonlinear schemes [8].

5 Theoretical analysis

In this section, we will prove the well-posedness and maximum principle of the discrete solution. In particular, we will demonstrate that $\lambda$ is essential to ensure the maximum principle. The unconditional energy stability of the proposed scheme will be proved as well.

5.1 Existence and uniqueness

For given $\phi^{n}$, we define the linear operator as follows

[TABLE]

Lemma 5.1**.**

Assume that $0<\phi^{n}<1$. For given $\theta>2$, if $\lambda$ is chosen such that

[TABLE]

then ${\bf Q}^{n}_{\theta,\lambda}$ is symmetric and positive definite.

Proof.

${\bf Q}^{n}_{\theta,\lambda}$ is symmetric since for $\phi,\varphi\in\mathcal{C}_{h}$, we have

[TABLE]

We now prove that ${\bf Q}^{n}_{\theta,\lambda}$ is positive definite. Let $\xi(s)=\frac{1}{s}+\frac{1}{1-s}$ for a scalar number $0<s<1$, then the derivative of $\xi(s)$ is expressed as

[TABLE]

We observe that $\xi^{\prime}(s)\leq 0$ for $s\in(0,\frac{1}{2}]$ and $\xi^{\prime}(s)\geq 0$ for $s\in[\frac{1}{2},1)$, thus

[TABLE]

For $\phi\in\mathcal{C}_{h}$ and $\phi\neq 0$, we deduce that

[TABLE]

which yields the positive definiteness of ${\bf Q}^{n}_{\theta,\lambda}$. ∎

Remark 5.1**.**

For $2<\theta<4$, the condition (5.2) holds if $\lambda=0$. But for $\theta\geq 4$, a positive stability constant satisfying $\lambda>\frac{1}{4}\theta-1$ is required to ensure the positive definiteness of ${\bf Q}^{n}_{\theta,\lambda}$. We also note that the linear system of the scheme with (4.16) is symmetric and positive definite for any $\theta>2$ and $\lambda\geq 0$.

We now prove the existence and uniqueness of the discrete solution of (4.18).

Theorem 5.1**.**

Assume that $0<\phi^{n}<1$. For any time step size $\tau>0$, there exists a unique $\phi^{n+1}$ to solve (4.18) in $\mathcal{C}_{h}$ provided that $\lambda$ is chosen to satisfy (5.2).

Proof.

It suffices to prove that the following homogeneous equation has a unique zero solution in $\mathcal{C}_{h}$

[TABLE]

which can be rewritten as

[TABLE]

As a direct consequence of the symmetric positive definiteness of ${\bf Q}^{n}_{\theta,\lambda}$, there is a unique solution $\phi\equiv 0$. ∎

5.2 Discrete maximum principle

The maximum principle of the discrete solution essentially relies on the stability constant for given energy parameter $\theta$, thus we need to employ the proper stability constant $\lambda$ to ensure this key property. We denote

[TABLE]

[TABLE]

In order to ensure the maximum principle of the proposed scheme, we introduce the following condition: for given $\theta>2$, $\lambda$ shall be taken to satisfy (5.2) and

[TABLE]

Remark 5.2**.**

We illustrate the reasonability of the condition (5.10) in Figure 2. Here, the values of $\lambda$ are taken as follows

[TABLE]

It can be observed from Figure 2 that the condition (5.10) naturally holds for $2<\theta\leq 3$ without any stability term, while it is still true for $3<\theta\leq 6$ when we take the small values (1 and 2) for the stability constant $\lambda$ as in (5.11). We have also checked some cases of $\theta>6$ and found that the condition (5.10) is true if $\lambda$ is properly chosen. We also emphasize that (5.10) is a good guide to choose $\lambda$ since it is clearly independent of discrete solutions. For the scheme with (4.16), we have the similar condition for the choice of the stability constant.

The following theorem demonstrates that the scheme (4.18) with the appropriate stability constant preserves the discrete maximum principle.

Theorem 5.2**.**

Assume that $0<\phi^{0}<1,$ and $\lambda$ is chosen to satisfy (5.2) and (5.10) for given $\theta>2$. For any time step size, the solutions of the scheme (4.18) satisfy $0<\phi^{n}<1$, $n\geq 1$.

Proof.

By induction, assuming $0<\phi^{n}<1$, we prove $0<\phi^{n+1}<1$. There exists a unique $\phi^{n+1}\in\mathcal{C}_{h}$ in terms of Theorem 5.1. Let us define $\phi_{-}^{n+1}=\min(\phi^{n+1},0)$. It is obtained from (4.18) that

[TABLE]

The definition of $\phi_{-}^{n+1}$ yields

[TABLE]

Moreover, it is apparent that $\phi_{-}^{n+1}\leq 0$. Due to $\phi^{n}>0$, we have

[TABLE]

Thus the first term on the left-hand side of (5.2) is bounded as

[TABLE]

Applying (3.8), the second term on the left-hand side of (5.2) is bounded as

[TABLE]

The third term on the left-hand side of (5.2) is bounded below

[TABLE]

Using the condition (5.10) and taking into account $\phi_{-}^{n+1}\leq 0$, the right-hand side of (5.2) is estimated as

[TABLE]

We combine (5.13)-(5.18) and obtain

[TABLE]

It follows from (5.19) that $\|\phi_{-}^{n+1}\|_{h}^{2}=0$ and consequently, we get $\phi^{n+1}\geq 0$. We further prove that $\phi^{n+1}>0$ by contradiction. Suppose that there exists at least a grid point $(i_{0}+\frac{1}{2},j_{0}+\frac{1}{2})$ such that $\phi^{n+1}_{i_{0}+\frac{1}{2},j_{0}+\frac{1}{2}}=0$, then the equation (4.18) at $(i_{0}+\frac{1}{2},j_{0}+\frac{1}{2})$ has the form

[TABLE]

The right-hand side of (5.20) can be estimated using (5.10) as

[TABLE]

For the two terms on the left-hand side of (5.20), we obtain

[TABLE]

which contradict (5.21). Therefore, $\phi^{n+1}>0$.

To prove $\phi^{n+1}<1$, we define $\phi_{+}^{n+1}=\max(\phi^{n+1}-1,0)$. We get from (4.18) that

[TABLE]

The definition of $\phi_{+}^{n+1}$ implies that

[TABLE]

Taking into account $\phi^{n}<1$, we have

[TABLE]

The first term of (5.2) is bounded below using (5.24) and (5.25)

[TABLE]

Similar to (5.16) and (5.17), we have

[TABLE]

[TABLE]

where (5.27) results from (3.9). Applying the condition (5.10) and $\phi_{+}^{n+1}\geq 0$, we derive

[TABLE]

Finally, we reach

[TABLE]

which implies that $\|\phi_{+}^{n+1}\|_{h}^{2}=0$, and thus, $\phi^{n+1}\leq 1$. We further prove that $\phi^{n+1}\neq 1$ in any grid cell by contradiction. Suppose that $\phi^{n+1}_{i_{0}+\frac{1}{2},j_{0}+\frac{1}{2}}=1$, and then the equation (4.18) at $(i_{0}+\frac{1}{2},j_{0}+\frac{1}{2})$ becomes

[TABLE]

The right-hand side of (5.31) can be estimated using (5.10) as

[TABLE]

while the two terms on the left-hand side of (5.31) are estimated as

[TABLE]

Substituting (5.32) and (5.33) into (5.31) yields a contradiction, and consequently, we obtain $\phi^{n+1}<1$. ∎

5.3 Unconditional energy stability

For $\phi\in\mathcal{C}_{h}$, the discrete total free energy is defined as

[TABLE]

Theorem 5.3**.**

Assume that $0<\phi^{0}<1$ and $\lambda$ is chosen to satisfy (5.2) and (5.10) for given $\theta>2$. For any time step size $\tau$, the total free energy of the scheme (4.18) is dissipated with time steps, i.e.,

[TABLE]

Proof.

Theorems 5.1 and 5.2 argue that there exists a unique $\phi^{n}\in\mathcal{C}_{h}$ and $0<\phi^{n}<1$ for $n\geq 1$. The inequality (4.14) yields

[TABLE]

On the other hand, using (3.7), we can derive that

[TABLE]

Using (4.17), (5.36) and (5.3), we deduce

[TABLE]

which yields (5.35). ∎

Remark 5.3**.**

We note that (5.10) is just a sufficient condition to ensure the discrete maximum principle. The energy dissipation law always holds as long as $0<\phi^{n}<1$, $n\geq 1$. We apply the discrete chemical potential (4.13) for solving the Cahn–Hilliard equation and find that $\phi^{n}$ always falls in $(0,1)$, thus the energy dissipation is still validated.

6 Numerical results

In this section, we present numerical results to show the effectiveness of the proposed scheme. In all numerical experiments, the computational domain is $\Omega=[-1,1]^{2}$, which is divided using a uniform mesh with $100\times 100$ elements. Here, we take $\epsilon=0.05$. The proposed scheme admits a very time step size in theory, so we take the time step size $\tau=10^{10}$ for the purpose of verifying this feature.

In this example, we compare the discrete chemical potentials proposed in (4.13) and (4.16). We take $\theta=3$ and $\lambda=0$. The initial condition of the phase variable is given as

[TABLE]

Figure 3 depicts that the maximum and minimum values of the phase variable computed by the proposed scheme with (4.13) or (4.16) at each time point. Although the initial maximum and minimum values of $\phi^{n}$ are set to approach the singular points, all of $\phi^{n}$ never go beyond $(0,1)$ in the subsequent time, and consequently the maximum principle is always preserved by both (4.13) and (4.16).

Figure 4 demonstrates that the total energies are always decreasing with time steps. The proposed scheme with (4.13) has more rapid dissipation rate than that with (4.16), thus the square-shape phase simulated by the former can be shrinking into a circle more rapidly. This phenomenon is observed in Figure 5, which illustrates the dynamical evolution process of the phase variable at different time steps. The superiority of (4.13) results from the semi-implicit treatment for $H_{c}(\phi)$, while it is treated explicitly in (4.16).

7 Conclusions

The stabilized energy factorization approach has been developed to treat the logarithmic Flory–Huggins potential semi-implicitly. The stability term can eliminate the instability caused by the large energy parameter. Compared to the prevalent convex-splitting approach and auxiliary variable approaches, such approach leads to the simply linear numerical scheme inheriting the original energy dissipation law. To our best knowledge, the proposed scheme is the first such linear, unconditionally original energy stable scheme for the logarithmic Flory–Huggins potential. Moreover, the proposed scheme is rigorously proved to satisfy the discrete maximum principle under the appropriate choices of the stability constant. Numerical results are in good agreement with theoretical analysis.

In addition, the proposed scheme can be extended to the Cahn–Hilliard equation and the corresponding numerical results (not presented here) demonstrate that the maximum principle and unconditional energy stability are still validated as long as the proper stability constant is taken although theoretical proof may be not available for the maximum principle of the Cahn–Hilliard equation.