Two Efficient Second Order Stabilized Semi-Implicit Schemes for the Cahn-Hilliard Phase-Field Equation

# Two Efficient Second Order Stabilized Semi-Implicit Schemes for the Cahn-Hilliard Phase-Field Equation

Lin Wang Haijun Yu NCMIS & LSEC, Institute of Computational Mathematics and Scientific/Engineering Computing, Academy of Mathematics and Systems Science, Beijing 100190, China School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
###### 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 semi-implicit schemes for the Cahn-Hilliard phase-field equation. The first one uses Crank-Nicolson method and the second one uses backward differentiation formula to discretize linear terms. In both schemes, the nonlinear bulk forces are treated explicitly with second-order 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 step-size 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, Cahn-Hilliard equation, energy stable, stabilized semi-implicit scheme, second order time marching

## 1 Introduction

In this paper, we consider numerical approximation for the Cahn-Hilliard equation

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩ϕt=γΔμ,(x,t)∈Ω×(0,T],μ=−εΔϕ+1εf(ϕ),(x,t)∈Ω×(0,T],ϕ|t=0=ϕ0(x),x∈Ω, (1.1)

with Neumann boundary condition

 ∂nϕ=0,∂nμ=0,x∈∂Ω. (1.2)

Here is a bounded domain with a locally Lipschitz boundary, is the outward normal, is a given time, is the phase-field 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 Cahn-Hilliard equation was originally introduced by Cahn and Hilliard cahn_free_1958 () to describe the phase separation and coarsening phenomena in non-uniform systems such as alloys, glasses and polymer mixtures. If the term in equation (1.1) is replaced with , one get the Allen-Cahn equation, which was introduced by Allen and Cahn allen_microscopic_1979 () to describe the motion of anti-phase boundaries in crystalline solids. The Cahn-Hilliard equation and the Allen-Cahn equation are two widely used phase-field models. In a phase-field 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 phase-field models due to the good mathematical properties that phase-field 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 phase-field 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 bi-harmonic operator makes the equations very stiff and the stiffness is not purely diffusive. In the Allen-Cahn equation, the stiffness is mainly due to the small parameter , while in the Cahn-Hilliard 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 Cahn-Hilliard equations based on the the mathematical properties of the equation. One of the basic properties of the Cahn-Hilliard equation is that they can be viewed as the gradient flow of the Ginzburg-Landau energy functional

 Eε(ϕ):=∫Ω(ε2|∇ϕ|2+1εF(ϕ))dx (1.3)

in . More precisely, by taking the inner product of (1.1) with , and taking integration, we immediately find the following energy law:

 Eε(ϕ(t))−Eε(ϕ0)=−γ∫t0∫Ω|∇μ|2dx=−γ∫t0∥ϕt∥2−1dx,∀t>0, (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 step-size 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 step-size 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 secant-line 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 secant-type method for Cahn-Hilliard equation need a small time step-size to guarantee the semi-discretized nonlinear system has a unique solution (cf. for instance du_numerical_1991 (); barrett_finite_1999 ()). To remove the restriction on time step-size, a diffusive three-step Crank-Nicolson scheme was introduced in guo_h2_2016 () and diegel_stability_2016 () coupled with a second order convex splitting. After time-discretization, 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 Cahn-Hilliard equation, it first appeared in guillen-gonzalez_linear_2013 (); guillen-gonzalez_second_2014 () as a Lagrange multiplier method. It then generalized by Yang et al. and successfully extended to handle several very complicated nonlinear phase-field 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 semi-implicit 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 well-structured linear system which is easier to solve. The modified energy in IEQ is an order-consistent approximation to the original system energy. At each time step, it needs to solve a linear system with time-varying coefficients.

Another trend of improving numerical schemes for phase-field models focuses on algorithm efficiency. Chen and Shen chen_applications_1998 (); zhu_coarsening_1999 () studied stabilized semi-implicit Fourier-spectral method for Cahn-Hilliard equation. The space variables are discretized by using a Fourier-spectral method whose convergence rate is exponential in contrast to the second order convergence of a usual finite-difference method. The time variable is discretized by using semi-implicit 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 time-stepping stabilized semi-implicit method for a 2-dimensional epitaxial growth model. He et al he_large_2007 () proposed a similar large time-stepping methods for the Cahn-Hilliard 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 Allen-Cahn equation and Cahn-Hilliard equation in mixed formulation shen_numerical_2010 (), which leads to unconditionally energy stable first-order linear schemes and second-order linear schemes with reasonable stability conditions. This idea was followed up in feng_stabilized_2013 () for the stabilized Crank-Nicolson schemes for phase field models. Another stabilized second-order Crank-Nicolson scheme with a new convex-concave splitting of the energy is proposed for a tumor-growth 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 convex-splitting 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 second-order time marching schemes based on Crank-Nicolson and second-order 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 second-order stabilized schemes for the Cahn-Hilliard 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 2-dimensional 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

 −Δv1=v  in Ω, ∂v1∂n=0  on ∂Ω,

and .

For any given function of , we use to denote an approximation of , where is the step-size. We will frequently use the shorthand notations: , , , and . Following identities will be used frequently as well

 2(a−b,a)=∥a∥2−∥b∥2+∥a−b∥2, (2.5)
 (Dτhn+1,hn+1)=14τ(∥hn+1∥2+∥2hn+1−hn∥2−∥hn∥2−∥2hn−hn−1∥2+∥δtthn+1∥2), (2.6)
 (δtthn+1,hn+1)=12(∥hn+1∥2−∥hn∥2)−12(∥hn∥2−∥hn−1∥2)+12∥δtthn+1∥2−∥δthn∥2. (2.7)

To prove energy stability of the numerical schemes, we assume that the derivative of in equation (1.1) is uniformly bounded:

 maxϕ∈R|f′(ϕ)|≤L, (2.8)

where is a non-negative constant.

### 2.1 The stabilized linear Crank-Nicolson scheme

Suppose and are given, our stabilized liner Crank-Nicolson scheme (abbr. SL-CN) calculates iteratively, using

 ϕn+1−ϕnτ=γΔμn+12, (2.9) μn+12=−εΔ(ϕn+1+ϕn2)+1εf(32ϕn−12ϕn−1)−AτΔδtϕn+1+Bδttϕn+1, (2.10)

where and are two non-negative constants to stabilize the scheme.

###### Theorem 2.1.

Under the condition

 A≥L216ε2γ,B≥L2ε, (2.11)

the following energy dissipation law

 En+1C≤EnC−(2√Aγ−L2ε)∥δtϕn+1∥2−(B2−L4ε)∥δttϕn+1∥2,∀n≥1, (2.12)

holds for the scheme (2.9)-(2.10), where

 En+1C=Eε(ϕn+1)+(L4ε+B2)∥δtϕn+1∥2. (2.13)
###### Proof.

Pairing (2.9) with , (2.10) with , and combining the results, we get

 ε2(∥∇ϕn+1∥2−∥∇ϕn∥2)+1ε(f(^ϕn+12),δtϕn+1)=−γτ∥∇μn+1∥2−Aτ∥∇δtϕn+1∥2−B(δttϕn+1,δtϕn+1). (2.14)

Pairing (2.9) with , using Cauchy-Schwartz inequality, we get

 2√Aγ∥δtϕn+1∥2=−2√Aγτ(∇μn+1,∇δtϕn+1)≤γτ∥∇μn+1∥2+Aτ∥∇δtϕn+1∥2. (2.15)

To handle the term involving , we expand and at as

 F(ϕn+1) =F(^ϕn+12)+f(^ϕn+12)(ϕn+1−^ϕn+12)+12f′(ξn1)(ϕn+1−^ϕn+12)2, F(ϕn) =F(^ϕn+12)+f(^ϕn+12)(ϕn−^ϕn+12)+12f′(ξn2)(ϕn−^ϕn+12)2,

where is a number between and , is a number between and . Taking the difference of above two equations, we have

 F(ϕn+1)−F(ϕn)−f(^ϕn+12)(ϕn+1−ϕn)=12f′(ξn1)[(ϕn+1−^ϕn+12)2−(ϕn−^ϕn+12)2]−12(f′(ξn2)−f′(ξn1))(ϕn−^ϕn+12)2=12f′(ξn1)δtϕn+1δttϕn+1−18(f′(ξn2)−f′(ξn1))(δtϕn)2≤L4(|δtϕn+1|2+|δttϕn+1|2)+L4|δtϕn|2.

Multiplying the above equation with , then taking integration leads to

 1ε(F(ϕn+1)−F(ϕn)−f(^ϕn+12)δtϕn+1,1)≤L4ε(∥δtϕn+1∥2+∥δttϕn+1∥2+∥δtϕn∥2). (2.16)

For the term involving , by using identity (2.5) with , one get

 −B(δttϕn+1,δtϕn+1)=−B2∥δtϕn+1∥2+B2∥δtϕn∥2−B2∥δttϕn+1∥2. (2.17)

Summing up (2.14)-(2.17), we obtain

 ε2(∥∇ϕn+1∥2−∥∇ϕn∥2)+1ε(F(ϕn+1)−F(ϕn),1)+B2∥δtϕn+1∥2−B2∥δtϕn∥2≤−2√Aγ∥δtϕn+1∥2+L4ε∥δtϕn+1∥2+L4ε∥δtϕn∥2−B2∥δttϕn+1∥2+L4ε∥δttϕn+1∥2, (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. SL-BDF2) calculate iteratively, using

 3ϕn+1−4ϕn+ϕn−12τ=γΔμn+1, (2.19) μn+1=−εΔϕn+1+1εf(2ϕn−ϕn−1)−AτΔδtϕn+1+Bδttϕn+1, (2.20)

where and are two non-negative constants.

###### Theorem 2.2.

Under the condition

 B≥Lε;A≥γα2L216ε2−α1ε2τ,0≤α1≤1,0<α2≤1, (2.21)

the following energy dissipation law

 En+1B≤EnB−14τγ∥δttϕn+1∥2−1−(1−α1)ε2∥∇δtϕn+1∥2−(1−α2)1τγ∥δtϕn+1∥2−1−(2√α2γ(A+α1ε2τ)−L2ε)∥δtϕn+1∥2−(B2−L2ε)∥δttϕn+1∥2,∀n≥1, (2.22)

holds for the scheme (2.19)-(2.20), where

 En+1B=Eε(ϕn+1)+14τγ∥δtϕn+1∥2−1+(L2ε+B2)∥δtϕn+1∥2. (2.23)
###### Proof.

Pairing (2.19) with , plus (2.20) paired with , we get

 (1γDτϕn+1,(−Δ)−1δtϕn+1)=ε(Δϕn+1,δtϕn+1)−1ε(f(^ϕn+1),δtϕn+1)−Aτ∥∇δtϕn+1∥2−B(δttϕn+1,δtϕn+1). (2.24)

By applying (2.5) and integration by parts, following identities hold

 −(1γDτϕn+1,(−Δ)−1δtϕn+1)=−1τγ∥δtϕn+1∥2−1−14τγ(∥δtϕn+1∥2−1−∥δtϕn∥2−1+∥δttϕn+1∥2−1), (2.25)
 ε(Δϕn+1,δtϕn+1)=−ε2(∥∇ϕn+1∥2−∥∇ϕn∥2+∥∇δtϕn+1∥2), (2.26)
 −B(δttϕn+1,δtϕn+1)=−B2∥δtϕn+1∥2+B2∥δtϕn∥2−B2∥δttϕn+1∥2. (2.27)

To handle the term involves , we expand and at as

 F(ϕn+1) =F(^ϕn+1)+f(^ϕn+1)(ϕn+1−^ϕn+1)+12f′(ξn1)(ϕn+1−^ϕn+1)2, F(ϕn) =F(^ϕn+1)+f(^ϕn+1)(ϕn−^ϕn+1)+12f′(ξn2)(ϕn−^ϕn+1)2,

where is a number between and , is a number between and . Taking the difference of above two equations, using the fact and , we obtain

 F(ϕn+1)−F(ϕn)−f(^ϕn+1)δtϕn+1=12f′(ξn1)(δttϕn+1)2−12f′(ξn2)(δtϕn)2≤L2|δttϕn+1|2+L2|δtϕn|2. (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

 χ∥∇δtϕn+1∥2+α2τγ∥δtϕn+1∥2−1≥2√χα2τγ∥δtϕn+1∥2, (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 Cahn-Hilliard equation coupled with the Navier-Stokes equations have a sharp-interface 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 SL-BDF2 scheme is stable with any , if

 τ≤8ε3L2γ.

If small enough, the combination of the term and term controls the term as well, since

 ∥δttϕn+1∥2≤2∥δtϕn+1∥2+2∥δtϕn∥2.

A direct calculation shows that the SL-BDF2 scheme is stable with any , , if

 τ≤8ε325L2γ. (2.31)
###### Remark 2.3.

Similar stabilization skills work for Allen-Cahn 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 semi-discretized scheme SL-BDF2 for the Cahn-Hilliard equation in the norm of (Refer WangYu2017b () for the error estimate of the SL-CN 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 :

1. , , and elsewhere. There exist two non-negative constants , such that

 ϕ2≤B0+B1F(ϕ),∀ϕ∈R. (3.32)
2. . and are uniformly bounded, or, satisfies (2.8) and

 maxϕ∈R|f′′(ϕ)|≤L2, (3.33)

where is a non-negative constant.

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 non-negative constants such that

 m0:=1|Ω|∫Ωϕ0(x)dx ∈(−1,1), (3.34) Eε(ϕ0):=ε2∥∇ϕ0∥2+1ε∥F(ϕ0)∥L1 ≲ε−2ρ1. (3.35)
###### Assumption 3.3.

By Assumption 3.2 and property (1.4), we know that the exact solution of (1.1) satisfies

1. , or .

We further assume that:

1. , or ,

2. , or ,

3. , or ,

4. , or .

Here are non-negative constants.

###### Assumption 3.4.

We also assume that an appropriate scheme is used to calculate the numerical solution at first step, such that

 m1:=1|Ω|∫Ωϕ1(x)dx=m0, (3.36)
 Eε(ϕ1)≤Eε(ϕ0)≲ε−2ρ1, (3.37)
 1τ∥ϕ1−ϕ0∥2−1≲ε−2ρ1, (3.38)
 1τ∥ϕ1−ϕ0∥2≲ε−2ρ1−2, (3.39)

and exist a constant , such that

 ∥e1∥2−1+τε∥∇e1∥2≲ε−σ1τ4. (3.40)

Note that the solution obtained by the first order stabilized linear scheme

 ϕ1−ϕ0τ=Δμ1,μ1=−εΔϕ1+1εf(ϕ0)+Bδtϕ1, (3.41)

with satisfies the assumption (3.36)-(3.40) with .

###### Lemma 3.1.

Suppose that satisfy Assumption 3.1, satisfy Assumption 3.2, and satisfy (2.21) with , . Then, the following estimates holds for the numerical solution of (2.19)-(2.20)

 1|Ω|∫Ωϕn(x)dx=m0,n=1,…,N+1, (3.42)
 max1≤n≤N+1{ε∥∇ϕn∥2+1εF(ϕn)+14τ∥δtϕn∥2−1+L2ε∥δtϕn∥2}≲ε−2σ2, (3.43)

where for , for .

###### Proof.

(i) Equation (3.42) is obtained by integrating (2.19).

(ii) Equation (3.43) is a direct result of the energy estimate (2.22). Actually, when condition (2.21) with is satisfied, we have

 En+1B+n+1∑k=2{14τ∥δttϕk∥2−1+12τ∥δtϕk∥2−1+ε4∥∇δtϕk∥2}≤E1B,0≤n≤N. (3.44)

By using assumptions (3.37), (3.38), (3.32) and

 ∥ϕi∥2≤B0|Ω|+B1∥F(ϕi)∥L1≲1+ε1−2ρ1,i=0,1, (3.45)

we get

 E1B=Eε(ϕ1)+14τ∥δtϕ1∥2−1+(L2ε+B2)∥δtϕ1∥2≲ε−2ρ1+ε−1(1+ε1−2ρ1)+≲ε−2max{1/2,ρ1}.

Note that if , we can get using (3.39). Taking the maximum of the left hand side of (3.44), we get

 max0≤n≤NEn+1B≲ε−2σ2. (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

 ∥en+1∥2−1+∥2en+1−en∥2−1+2Aτ2∥∇en+1∥2+2Aτ2∥δt∇en+1∥2+τε∥∇en+1∥2+∥δtten+1∥2−1+4Bτ∥en+1∥2≤∥en∥2−1+∥2en−en−1∥2−1+2Aτ2∥∇en∥2+C1τε−3∥2en−en−1∥2−1+C2τ4ε−max{ρ3+1,ρ4+3,ρ5+5},n≥1, (3.47)

and

 max1≤n≤N{∥en+1∥2−1+∥2en+1−en∥2−1+2Aτ2∥∇en+1∥2}+N∑n=1(2Aτ