Energy Stable Second Order Linear Schemes for the AllenCahn PhaseField Equation^{†}^{†}thanks: Feb 14, 2018, and accepted date (The correct dates will be entered by the editor).
Abstract
Phasefield model is a powerful mathematical tool to study the dynamics of interface and morphology changes in fluid mechanics and material sciences. However, numerically solving a phase field model for a real problem is a challenge task due to the nonconvexity of the bulk energy and the small interface thickness parameter in the equation. In this paper, we propose two stabilized second order semiimplicit linear schemes for the AllenCahn phasefield equation based on backward differentiation formula and CrankNicolson method, respectively. In both schemes, the nonlinear bulk force is treated explicitly with two secondorder stabilization terms, which make the schemes unconditional energy stable and numerically efficient. By using a known result of the spectrum estimate of the linearized AllenCahn operator and some regularity estimate of the exact solution, we obtain an optimal second order convergence in time with a prefactor depending on the inverse of the characteristic interface thickness only in some lower polynomial order. Both 2dimensional and 3dimensional numerical results are presented to verify the accuracy and efficiency of proposed schemes.
Keywords. AllenCahn equation; energy stable; stabilized semiimplicit scheme; second order scheme; error estimate
AMS subject classifications. 65M12; 65M15; 65P40
1 Introduction
In this paper, we consider numerical approximation for the AllenCahn equation with Neumann boundary condition
\hb@xt@.01(1.1)  
\hb@xt@.01(1.2) 
Here is a bounded domain with a locally Lipschitz boundary, is the outward normal, is a given time, is the phasefield variable. , the bulk force, is the derivative of a given energy function , which is usually nonconvex with two or more than two local minima. One commonly used energy function for twophase problem is the doublewell potential . is the thickness of the interface between two phases. , called mobility, is related to the characteristic relaxation time of the system. The homogeneous Neumann boundary condition implies that no mass loss occurs across the boundary walls. The equation (LABEL:eq:AC) is introduced by Allen and Cahn [1] to describe the process of phase separation in multicomponent alloy systems. It can be regarded as the gradient flow with respect to the GinzburgLandau energy functional
\hb@xt@.01(1.3) 
The corresponding energy dissipation is given as
\hb@xt@.01(1.4) 
Another popular phase field model is the CahnHilliard equation, which is the the gradient flow with respect to the GinzburgLandau energy functional. It was originally introduced by Cahn and Hilliard [4] to describe the phase separation and coarsening phenomena in nonuniform systems such as alloys, glasses and polymer mixtures.
The AllenCahn equation and the CahnHilliard equation are widely used in modeling many interface problems due to their good mathematical properties (cf. e.g. [16, 7, 15, 14, 53, 44] ). However, the small parameter and the nonconvexity of energy function make the numerical approximation of a phase field equation a challenging task, especially the design of time marching schemes. It is wellknown that if a fully explicit or implicit time marching scheme is used, a tiny time stepsize is required for the semidiscretized scheme to be stable or uniquely solvable since the nonlinear function is neither convex nor concave. A very popular approach to obtain unconditional stable time marching schemes is the so called convex splitting method which appears to be introduced by Elliott and Stuart [17], and popularized by Eyre [18], in which, the convex part of is treated implicitly and the concave part of is treated explicitly. This method has been applied to various gradient flows (see e.g. [2, 40, 21, 22, 19]). Traditional convex splitting schemes are first order accurate. Recently, several extensions to second order schemes were proposed based on either the CrankNicolson scheme (see e.g.[2, 6, 29, 12, 8, 36]), or second order backward differentiation formula (BDF2) [47, 35]. In all convex splitting schemes, no matter first order or second order, one usually obtains an uniquely solvable nonlinear convex problem at each time step.
There are another types of second order unconditional stable schemes for the phase field equations. In [13], Du and Nicolaides proposed a secantline method which is energy stable and second order accurate. It is used and extended in several other works, e.g. [25, 20, 9, 26, 2, 54, 3]. Similar to the convex splitting method, the secantline method leads to nonlinear semidiscretized system, which need special efforts to solve. Recently, an augmented Lagrange multiplier(ALM) method was proposed in [27, 28] to get second order linear energy stable schemes. The idea is generalized as invariant energy quadratization (IEQ) by Yang et al. and successfully applied to handle several very complicated nonlinear phasefield models (see e.g. [49, 30, 50, 51]). Based on similar methodology, a new variant called scalar auxiliary variable (SAV) method is developed by Shen et al. [38, 39]. In the ALM and IEQ approach, nonlinear semidiscretized systems are avoided, but one has to solve variablecoefficient systems, while in the SAV scheme, one only need to solve some linear systems with constant coefficients. Different to other methods, the energy in ALM, IEQ and SAV approach is a modified one which also depends on the auxiliary variable.
Our methods are based on semiimplicit discretization and stabilization skill. To improve the numerical stability of solving phasefield equations, semiimplicit schemes were proposed by Chen and Shen[5] and Zhu et al.[55]. Although not unconditionally stable, semiimplicit schemes allow much larger time stepsizes than explicit schemes. To further improve the stability, Xu and Tang proposed stabilized semiimplicit methods for epitaxial growth model in [45]. The proposed schemes have extraordinary numerical stability even though the mathematical proof of the stability is not complete. Similar schemes was developed for phase field equation by He et al.[31] and Shen and Yang[40], where the latter one adopted a mixed form for the CahnHilliard equation, by using a truncated double well potential such that the assumption is satisfied, the unconditional energy stability was proved for the first order stabilized scheme. It is worth to mention that with no truncation made to , Li et al [34, 33] proved that the energy stable property can be obtained as well, but a much larger stability constant need be used.
In this paper, we develop two secondorder unconditionally energy stable linear schemes for the AllenCahn equation based on the schemes proposed in [45] and [40]. The energy dissipation is guaranteed by including two second order stabilization terms, the first one is directly from [45], the other one is inspired by the work [44]. We also carry out an optimal error estimate for the time semidiscretized schemes. For the phase field equations, the error bounds will depend on the factor of exponentially if one uses a standard procedure. By using a spectrum estimate result of de Mottoni and Schatzman [10, 11] and Chen[7] for the linearized AllenCahn operator, we are able to get an optimal error estimate with a prefactor depend on only in some lower polynomial order for small . This spectrum estimate argument was first used by Feng and Prohl [23, 24] for an implicit first order scheme for phase field equations. It was also applied by Kessler et al. [32] to derive a posteriori error estimate for adaptive time marching. Similar analysis for a firstorder stabilized semiimplicit scheme of the AllenCahn equation in given by Yang [48]. Recently, Feng and Li [21] , Feng et al. [22] extended this spectrum estimate argument to first order convex splitting scheme coupled with interior penalty discontinuous Galerkin spatial discretization for AllenCahn and CahnHilliard equation, respectively. To our best knowledge, our analysis is the first such result for second order linear schemes. In summary, the proposed methods have several merits: 1) They are second order accurate; 2) They lead to linear systems with constant coefficients after time discretization; 3) The stability and error analysis bases on weak formulations, so both finite element method and spectral method can be used for spatial discretization to satisfy discretized energy dissipation law. 4) The methods can be easily used in more complicated systems. Note that, similar approach can be extended to the CahnHilliard equation [42, 43], where Lipschitz condition of is assumed based on physical intuition and the analyses are more tedious.
The remain parts of the paper is organized as follows. In Section 2, we present the two secondorder stabilized schemes for the AllenCahn equation and prove they are energy stable. The error estimate to derive a convergence rate that does not depend on exponentially is then constructed in Section 3. Detailed implementation and numerical experiments for problems in both 2dimensional and 3dimensional tensorproduct domain are presented in Section 4 to verify our theoretical results. We end the paper with some conclusions in Section 5.
2 The two second order stabilized linear schemes
We first introduce some notations which will be used throughout the paper. We use to denote the standard norm of the Sobolev space . In particular, we use to denote the norm of ; to denote the norm of ; and to denote the norm of . Let represent the inner product. For , we define , and denote .
For any given function of , we use to denote an approximation of , where is the stepsize. We will frequently use the shorthand notations: , , , and . Following identities will be used frequently as well
\hb@xt@.01(2.1)  
\hb@xt@.01(2.2) 
To prove energy stability of the numerical schemes, we assume that the derivative of in equation (LABEL:eq:AC) is uniformly bounded, i.e.
\hb@xt@.01(2.3) 
where is a nonnegative constant.
Remark 2.1
Note that the commonly used double well potential does not satisfy the above assumption. But, thanks to the maximum principle that the AllenCahn equation has (cf. e.g. [1, 7, 23, 48]), the solution to equation (LABEL:eq:AC) will be bounded by value and if the initial condition is bounded by and . So it is safe to modify the doublewell energy for larger than to be quadratic growth without affecting the exact solution if the initial condition is bounded by and , such that assumption (LABEL:eq:Lip) is satisfied. This argument also applies to the assumption (LABEL:eq:Lip2) in next section.
2.1 The stabilized linear BDF2 scheme.
Suppose and are given, our stabilized linear BDF2 scheme (SLBDF2) calculate iteratively, using
\hb@xt@.01(2.4) 
where and are two nonnegative constants to stabilize the scheme.
Theorem 2.1
Assume that (LABEL:eq:Lip) is satisfied. Under the condition
\hb@xt@.01(2.5) 
the following energy dissipation law
\hb@xt@.01(2.6) 
holds for the scheme (LABEL:BDF1), where
\hb@xt@.01(2.7) 
Proof. Pairing (LABEL:BDF1) with , we get
\hb@xt@.01(2.8) 
By integration by parts, following identities hold
\hb@xt@.01(2.9)  
\hb@xt@.01(2.10)  
\hb@xt@.01(2.11) 
To handle the term involves in (LABEL:cs7), we expand and at as
where is a number between and , is a number between and . Taking the difference of above two equations, using the fact and , we obtain
\hb@xt@.01(2.12) 
Taking inner product of the above equation with constant , then combining the result with (LABEL:cs7), (LABEL:cs8), (LABEL:cs9) and (LABEL:cs111), we obtain
\hb@xt@.01(2.13) 
Combining the above equation and the inequality , we get energy dissipation law (LABEL:eq:BDF:Edis).
Remark 2.2
From equation (LABEL:eq:BDF:ABcond), we see that the SLBDF2 scheme is stable with any nonnegative including , if
\hb@xt@.01(2.14) 
If one takes time step size even smaller,
\hb@xt@.01(2.15) 
then the SLBDF2 scheme is stable with any any nonnegative and , including the case .
On the other hand side, if we take
\hb@xt@.01(2.16) 
then the SLBDF2 scheme is unconditional stable for any .
2.2 The stabilized linear CrankNicolson scheme.
Suppose and are given, our stabilized linear CrankNicolson scheme (SLCN) calculate iteratively, using
\hb@xt@.01(2.17) 
where and are two nonnegative constants.
Theorem 2.2
Assume that (LABEL:eq:Lip) is satisfied. Under the condition
\hb@xt@.01(2.18) 
the following energy law holds
\hb@xt@.01(2.19) 
for the scheme (LABEL:accn1), where we define
\hb@xt@.01(2.20) 
Proof. Pairing the equation (LABEL:accn1) with , we get
\hb@xt@.01(2.21) 
We use Taylor expansion at ,
\hb@xt@.01(2.22)  
\hb@xt@.01(2.23) 
Subtracting (LABEL:acta2) from (LABEL:acta1) and the definition of , we have
\hb@xt@.01(2.24) 
which give us
\hb@xt@.01(2.25) 
Plugging (LABEL:accn6) into (LABEL:accn5), we obtain
\hb@xt@.01(2.26) 
By the definition of and , , we get the desired results.
Remark 2.3
If we take
\hb@xt@.01(2.27) 
then the SLCN scheme is unconditional stable for any .
On the other hand, by using the inequality , it is easy to prove that when , the SLCN scheme (LABEL:accn1) is stable for
\hb@xt@.01(2.28) 
Remark 2.4
To make SLBDF2 and SLCN scheme be unconditionally stable, i.e. stable for any time step size , we need take . This seems that need to be very large in a real simulation. But actually, it is not necessary. It was showed by Magaletti et al. [37] and Xu et al. [46] that the phasefield CahnHilliard–NavierStokes model for binary fluids has a fast convergence with respect to when the phenomenological mobility . We expect a similar scaling of in the AllenCahn model, which means is not necessary as large as . Furthermore, it is proved that the numerical interface for the AllenCahn equatoin converges with the rate before the singularities appear[23, 21], which suggests that we don’t need to take as small as a physical interface.
Remark 2.5
Recently, Li, Qiao and Tang [34], Li and Qiao [33] studied several first order and second order stabilized semiimplicit Fourier schemes, respectively, for the CahnHilliard equation with doublewell potential
\hb@xt@.01(2.29) 
Without a Lipschitz condition on , they proved that those schemes are unconditionally stable when very large stability constant used. For example, according to Theorem 1.3 in [33], for a classical second order semiimplicit stabilized scheme proposed by Xu and Tang [45] applied to the CahnHilliard equation, the stabilization constant need to be as large as to make the scheme unconditionally stable (Note that the in [33] corresponds to in this paper). However, the constants in our schemes are only of order , respectively. The reasons are in two aspects. Firstly, for the AllenCahn equation, since its solution satisfies a maximum principle, it is reasonable to modify the defined in (LABEL:eq:doublewell0) for , such that the Lipschitz condition (LABEL:eq:Lip) is satisfied. Secondly, we use two stabilization terms instead of only one stabilization term. Our approach can be extended to the CahnHilliard equation with quadratic growth energy as well[43, 42].
3 Convergence analysis
In this section, we shall establish the error estimate of the two proposed schemes for the AllenCahn equation in the norm of . We will shown that, if the interface is well developed in the initial condition, the error bounds depend on only in some lower polynomial order for small . Let be the exact solution at time to the AllenCahn equation (LABEL:eq:AC) and be the solution at time to the time discrete numerical scheme (LABEL:BDF1) (or (LABEL:accn1)), we define error function . Obviously .
Before presenting the detailed error analysis, we first make some assumptions. For simplicity, we take in this section, and assume . We use notation in the way that means that with positive constant independent of .
Assumption 3.1
We make following assumptions on : , for , such that and are uniformly bounded, i.e. satisfies (LABEL:eq:Lip) and
\hb@xt@.01(3.1) 
where is a nonnegative constant.
Since the solution of AllenCahn equation satisfies maximum principle (see Remark LABEL:rmk:2.1), one can always modify for large such that Assumption LABEL:ap:1 hold without affecting the exact solution.
Assumption 3.2

We assume that there exist nonnegative constants such that
\hb@xt@.01(3.2) \hb@xt@.01(3.3) \hb@xt@.01(3.4) \hb@xt@.01(3.5) 
Assume that an appropriate scheme is used to calculate the numerical solution at first step, such that
\hb@xt@.01(3.6) \hb@xt@.01(3.7) Then it is easy to get
\hb@xt@.01(3.8) \hb@xt@.01(3.9) 
There exist a constant ,
\hb@xt@.01(3.10)
Given Assumption LABEL:ap:1 LABEL:ap:02 (i), we have following estimates for the exact solution to the AllenCahn equation.
Lemma 3.1
Let be the exact solution of
(LABEL:eq:AC), under the condition of Assumption
LABEL:ap:1 and LABEL:ap:02 (i),
the following regularities holds:

;

;

;

;

.
Proof. Take in equation (LABEL:eq:AC), we have
\hb@xt@.01(3.11) 

Pairing (LABEL:eq:AC11) with and taking integration by parts on the second term, we get
\hb@xt@.01(3.12) After integration over and using the inequality (LABEL:ap:022), we obtain (i).

We differentiate (LABEL:eq:AC11) in time to obtain
\hb@xt@.01(3.13) Pairing (LABEL:eq:AC21) with yields
\hb@xt@.01(3.14) Integrate (LABEL:eq:AC22) over , yields
\hb@xt@.01(3.15) The assertion then follows from (i) and the inequality (LABEL:init1).

Testing (LABEL:eq:AC21) with , we get
\hb@xt@.01(3.16) Integrating (LABEL:eq:AC31) over , we get
\hb@xt@.01(3.17) and by using (i) and the inequality (LABEL:init2) of Assumption LABEL:ap:02, we obtain (iii).

We differentiate (LABEL:eq:AC21) in time to derive
\hb@xt@.01(3.18) Testing (LABEL:eq:AC41) with and using for , we have
\hb@xt@.01(3.19) Integrating (LABEL:eq:AC42) over , we obtain
\hb@xt@.01(3.20) The assertion then follows from (i) (ii) (iii) and the inequality (LABEL:init3).

Testing (LABEL:eq:AC41) with , we have