Two Efficient Second Order Stabilized SemiImplicit Schemes for the CahnHilliard PhaseField Equation
Abstract
Efficient and energy stable high order time marching schemes are very important but not easy to construct for the study of nonlinear phase dynamics. In this paper, we propose two linearly stabilized second order semiimplicit schemes for the CahnHilliard phasefield equation. The first one uses CrankNicolson method and the second one uses backward differentiation formula to discretize linear terms. In both schemes, the nonlinear bulk forces are treated explicitly with secondorder linear stabilization terms. This treatment leads to linear equations with constant coefficients, thus robust and efficient solution procedures are guaranteed. Moreover, both schemes are unconditionally stable. Rigorous error analysis shows that, when the time stepsize is small enough, second order accurate in time is obtained with a prefactor controlled by a fixed power of , where is the characteristic interface thickness. Numerical results are presented to verify the accuracy and efficiency of proposed schemes.
keywords:
phase field model, CahnHilliard equation, energy stable, stabilized semiimplicit scheme, second order time marching1 Introduction
In this paper, we consider numerical approximation for the CahnHilliard equation
(1.1) 
with Neumann boundary condition
(1.2) 
Here is a bounded domain with a locally Lipschitz boundary, is the outward normal, is a given time, is the phasefield variable. Function , with is a given energy potential with two local minima, e.g. the double well potential . The two minima of produce two phases, with the typical thickness of the interface between two phases given by . , called mobility, is related to the characteristic relaxation time of the system.
The CahnHilliard equation was originally introduced by Cahn and Hilliard cahn_free_1958 () to describe the phase separation and coarsening phenomena in nonuniform systems such as alloys, glasses and polymer mixtures. If the term in equation (1.1) is replaced with , one get the AllenCahn equation, which was introduced by Allen and Cahn allen_microscopic_1979 () to describe the motion of antiphase boundaries in crystalline solids. The CahnHilliard equation and the AllenCahn equation are two widely used phasefield models. In a phasefield model, the information of interface is encoded in a smooth phase function . In most parts of the domain , the value of is close to local minima of . The interface is a thin layer of thickness connecting regions of different local minima. It is easier to deal with dynamical process involving morphology changes of interfaces using phasefield models due to the good mathematical properties that phasefield equations have. For this reason, phase field models have been the subject of many theoretical and numerical investigations for several decades(cf., for instance, du_numerical_1991 (), elliott_error_1992 (), chen_spectrum_1994 (), caffarelli_l_1995 (), elliott_cahnhilliard_1996 (), eyre_unconditionally_1998 (), furihata_stable_2001 (), liu_phase_2003 (), feng_error_2004 (), kessler_posteriori_2004 (), shen_numerical_2010 (), condette_spectral_2011 ()).
But, numerically solving the phasefield equations is not an easy task. Firstly, the small parameter requires a very high spatial and temporal grid resolution. Secondly, the small parameter in the nonlinear bulk force and the biharmonic operator makes the equations very stiff and the stiffness is not purely diffusive. In the AllenCahn equation, the stiffness is mainly due to the small parameter , while in the CahnHilliard equation, stiffness due to is mixed with stiff Laplace operator, which makes it much harder to solve. Nevertheless, a lot of numerical schemes have been proposed to solving the CahnHilliard equations based on the the mathematical properties of the equation. One of the basic properties of the CahnHilliard equation is that they can be viewed as the gradient flow of the GinzburgLandau energy functional
(1.3) 
in . More precisely, by taking the inner product of (1.1) with , and taking integration, we immediately find the following energy law:
(1.4) 
where the is defined in Section 2. Since the nonlinear energy is neither a convex nor a concave function, treating it fully explicit or implicit in a time discretization will not lead to an efficient scheme. In fact, if the nonlinear force is treated fully explicitly, the resulting scheme will require a very tiny stepsize to be stable(cf. for instance shen_numerical_2010 ()). On the other hand, treating it fully implicitly will lead to a nonlinear system, for which the solution existence and uniqueness requires a restriction on stepsize as well (cf. e.g. feng_error_2004 ()). One popular approach to solve this dilemma is the convex splitting method introduced by Eyre eyre_unconditionally_1998 (), in which the convex part of is treated implicitly and the concave part treated explicitly. The scheme is of first order accurate and unconditional stable. In each time step, one need solve a nonlinear system. The solution existence and uniqueness is guaranteed since the nonlinear system corresponds to a convex optimization problem. The convex splitting method was used widely, and several second order extensions were derived in different situations condette_spectral_2011 (); baskaran_energy_2013 (); chen_linear_2014 (); guo_h2_2016 (), etc. Another type unconditional stable scheme is the secantline method proposed by Du and Nicolaidesdu_numerical_1991 (). It is also used and extended in several other works, e.g. furihata_stable_2001 (); kim_conservative_2004 (); feng_fully_2006 (); condette_spectral_2011 (); gomez_provably_2011 (); baskaran_energy_2013 (); zhang_adaptive_2013 (); benesova_implicit_2014 (). Like the fully implicit method, the usual second order convex splitting method and the secanttype method for CahnHilliard equation need a small time stepsize to guarantee the semidiscretized nonlinear system has a unique solution (cf. for instance du_numerical_1991 (); barrett_finite_1999 ()). To remove the restriction on time stepsize, a diffusive threestep CrankNicolson scheme was introduced in guo_h2_2016 () and diegel_stability_2016 () coupled with a second order convex splitting. After timediscretization, one get a nonlinear but unique solvable problem at each time step.
Recently, a new approach termed as invariant energy quadratization (IEQ) was introduced to handle the nonlinear energy. When applying to CahnHilliard equation, it first appeared in guillengonzalez_linear_2013 (); guillengonzalez_second_2014 () as a Lagrange multiplier method. It then generalized by Yang et al. and successfully extended to handle several very complicated nonlinear phasefield models yang_linear_2016 (); han_numerical_2017 (); yang_efficient_2017 (); yang_yu_efficient_2017 (). In the IEQ approach, a new variable which equals to the square root of is introduced, so the energy is written into a quadratic form in terms of the new variable. By using semiimplicit treatments to all the nonlinear terms in the equations, one get a linear and energy stable scheme. It is straightforward to prove the unconditional stability for both first order and second order IEQ schemes. Comparing to the convex splitting approach, IEQ leads to wellstructured linear system which is easier to solve. The modified energy in IEQ is an orderconsistent approximation to the original system energy. At each time step, it needs to solve a linear system with timevarying coefficients.
Another trend of improving numerical schemes for phasefield models focuses on algorithm efficiency. Chen and Shen chen_applications_1998 (); zhu_coarsening_1999 () studied stabilized semiimplicit Fourierspectral method for CahnHilliard equation. The space variables are discretized by using a Fourierspectral method whose convergence rate is exponential in contrast to the second order convergence of a usual finitedifference method. The time variable is discretized by using semiimplicit schemes which allow much larger time step sizes than explicit schemes. Xu and Tang xu_stability_2006 () introduced a different stabilized term to build large timestepping stabilized semiimplicit method for a 2dimensional epitaxial growth model. He et al he_large_2007 () proposed a similar large timestepping methods for the CahnHilliard equation, in which a stabilized term (resp. ) is added to the nonlinear bulk force for the first order(resp. second order) scheme. Shen and Yang applied similar stabilization skill to AllenCahn equation and CahnHilliard equation in mixed formulation shen_numerical_2010 (), which leads to unconditionally energy stable firstorder linear schemes and secondorder linear schemes with reasonable stability conditions. This idea was followed up in feng_stabilized_2013 () for the stabilized CrankNicolson schemes for phase field models. Another stabilized secondorder CrankNicolson scheme with a new convexconcave splitting of the energy is proposed for a tumorgrowth system by Wu et al.wu_stabilized_2014 (). Those time marching schemes all lead to linear systems, which are easier to solve than nonlinear systems resulting from traditional convexsplitting schemes, in which the nonlinear convex force is treated implicitly. On the other hand, when the nonlinear force is treated explicitly, one need to introduce a proper stabilization term and a suitably truncated nonlinear function instead of to prove the unconditionally energy stable property with a reasonable stabilization constant. It is worth to mention that with no truncation made to , Li et al li_characterizing_2016 (); li_second_2017 () proved that the energy stable property can be obtained as well, but a much larger stability constant need be used.
In this paper, we propose two secondorder time marching schemes based on CrankNicolson and secondorder backward differentiation formula (BDF2) approximation respectively. In both schemes, explicit extrapolation are used for the nonlinear force with small extra stabilization terms which consist to the order of the schemes added to guarantee energy dissipation. We also give an error analysis in norm. The new 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 analysis bases on Galerkin formulation, so both finite element method and spectral method can be used for spatial discretization to conserve volume fraction and satisfy discretized energy dissipation law.
The remain part of this paper is organized as follows. In Section 2, we present the two secondorder stabilized schemes for the CahnHilliard equation and prove they are energy stable. In Section 3, we present an error estimate to derive a convergence rate that does not depend on exponentially. Detailed implementation and numerical results for problems in a 2dimensional square domain are presented in Section 4 to verify our theoretical results. We end the paper with some concluding remarks 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. In addition, define for
where stands for the dual product between and . We denote . For , let , where is the solution to
and .
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
(2.5) 
(2.6) 
(2.7) 
To prove energy stability of the numerical schemes, we assume that the derivative of in equation (1.1) is uniformly bounded:
(2.8) 
where is a nonnegative constant.
2.1 The stabilized linear CrankNicolson scheme
Suppose and are given, our stabilized liner CrankNicolson scheme (abbr. SLCN) calculates iteratively, using
(2.9)  
(2.10) 
where and are two nonnegative constants to stabilize the scheme.
Theorem 2.1.
Proof.
Pairing (2.9) with , (2.10) with , and combining the results, we get
(2.14) 
Pairing (2.9) with , using CauchySchwartz inequality, we get
(2.15) 
To handle the term involving , we expand and at as
where is a number between and , is a number between and . Taking the difference of above two equations, we have
Multiplying the above equation with , then taking integration leads to
(2.16) 
For the term involving , by using identity (2.5) with , one get
(2.17) 
Summing up (2.14)(2.17), we obtain
(2.18) 
which is the energy estimate (2.12). ∎
2.2 The stabilized linear BDF2 scheme
Suppose and are given, our stabilized linear BDF2 scheme (abbr. SLBDF2) calculate iteratively, using
(2.19)  
(2.20) 
where and are two nonnegative constants.
Theorem 2.2.
Proof.
Pairing (2.19) with , plus (2.20) paired with , we get
(2.24) 
By applying (2.5) and integration by parts, following identities hold
(2.25) 
(2.26) 
(2.27) 
To handle the term involves , 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
(2.28) 
Taking inner product of the above equation with constant , then combining the result with (2.24), (2.25), (2.26) and (2.27), we obtain
(2.29) 
Combining the above equation and the inequality
(2.30) 
with , , we get the energy estimate (2.22). ∎
Remark 2.1.
The constant defined in equation (2.11) and (2.21) seems to be quite large when is small, but it is not necessarily true. Since usually is a small constant related to . For example, it was pointed out in magaletti2013sharp () that, the CahnHilliard equation coupled with the NavierStokes equations have a sharpinterface limit when , while gives the fastest convergence. On the other hand, the numerical results in Section 4 shows can take much smaller values than those defined in (2.11) and (2.21) when nonzero values are used.
Remark 2.2.
From equation (2.21), we see that the SLBDF2 scheme is stable with any , if
If small enough, the combination of the term and term controls the term as well, since
A direct calculation shows that the SLBDF2 scheme is stable with any , , if
(2.31) 
Remark 2.3.
Similar stabilization skills work for AllenCahn equation as well, but the stabilization term should be changed to . Please refer WangYu2017b () for more details.
3 Convergence analysis
In this section, we shall establish the error estimate of the semidiscretized scheme SLBDF2 for the CahnHilliard equation in the norm of (Refer WangYu2017b () for the error estimate of the SLCN scheme). 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 equation (1.1) and be the solution at time to the time discrete numerical scheme (2.19)(2.20), 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 sense that means that with positive constant independent of .
Assumption 3.1.
We make following assumptions on :

, , and elsewhere. There exist two nonnegative constants , such that
(3.32)
Some regularities of the initial condition (i.e. ) and exact solution are necessary for the error estimates. To save space, we give these regularities results as assumptions.
Assumption 3.2.
We assume , and there exist constant and nonnegative constants such that
(3.34)  
(3.35) 
Assumption 3.3.
Assumption 3.4.
We also assume that an appropriate scheme is used to calculate the numerical solution at first step, such that
(3.36) 
(3.37) 
(3.38) 
(3.39) 
and exist a constant , such that
(3.40) 
Note that the solution obtained by the first order stabilized linear scheme
(3.41) 
Lemma 3.1.
Proof.
(ii) Equation (3.43) is a direct result of the energy estimate (2.22). Actually, when condition (2.21) with is satisfied, we have
(3.44) 
By using assumptions (3.37), (3.38), (3.32) and
(3.45) 
we get
Note that if , we can get using (3.39). Taking the maximum of the left hand side of (3.44), we get
(3.46) 
∎
The volume conservation property (3.42) is the basis for the error estimate. Because of (3.42), and belong to such that we can define norm and use Poincare inequality for those quantities.
Proposition 3.0.
(Coarse error estimate) , following error estimate holds
(3.47) 
and