Energy Stable Second Order Linear Schemes for the Allen-Cahn Phase-Field EquationFeb 14, 2018, and accepted date (The correct dates will be entered by the editor).

Energy Stable Second Order Linear Schemes for the Allen-Cahn Phase-Field Equationthanks: Feb 14, 2018, and accepted date (The correct dates will be entered by the editor).

Lin Wang 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 (wanglin@lsec.cc.ac.cn).    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 (hyu@lsec.cc.ac.cn).
Abstract

Phase-field 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 non-convexity of the bulk energy and the small interface thickness parameter in the equation. In this paper, we propose two stabilized second order semi-implicit linear schemes for the Allen-Cahn phase-field equation based on backward differentiation formula and Crank-Nicolson method, respectively. In both schemes, the nonlinear bulk force is treated explicitly with two second-order stabilization terms, which make the schemes unconditional energy stable and numerically efficient. By using a known result of the spectrum estimate of the linearized Allen-Cahn 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 2-dimensional and 3-dimensional numerical results are presented to verify the accuracy and efficiency of proposed schemes.

Keywords. Allen-Cahn equation; energy stable; stabilized semi-implicit scheme; second order scheme; error estimate

AMS subject classifications. 65M12; 65M15; 65P40

1 Introduction

In this paper, we consider numerical approximation for the Allen-Cahn 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 phase-field variable. , the bulk force, is the derivative of a given energy function , which is usually non-convex with two or more than two local minima. One commonly used energy function for two-phase problem is the double-well 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 multi-component alloy systems. It can be regarded as the gradient flow with respect to the Ginzburg-Landau 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 Cahn-Hilliard equation, which is the the gradient flow with respect to the Ginzburg-Landau energy functional. It was originally introduced by Cahn and Hilliard [4] to describe the phase separation and coarsening phenomena in non-uniform systems such as alloys, glasses and polymer mixtures.

The Allen-Cahn equation and the Cahn-Hilliard 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 non-convexity of energy function make the numerical approximation of a phase field equation a challenging task, especially the design of time marching schemes. It is well-known that if a fully explicit or implicit time marching scheme is used, a tiny time step-size is required for the semi-discretized 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 Crank-Nicolson 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 secant-line 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 secant-line method leads to nonlinear semi-discretized 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 phase-field 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 semi-discretized systems are avoided, but one has to solve variable-coefficient 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 semi-implicit discretization and stabilization skill. To improve the numerical stability of solving phase-field equations, semi-implicit schemes were proposed by Chen and Shen[5] and Zhu et al.[55]. Although not unconditionally stable, semi-implicit schemes allow much larger time step-sizes than explicit schemes. To further improve the stability, Xu and Tang proposed stabilized semi-implicit 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 Cahn-Hilliard 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 second-order unconditionally energy stable linear schemes for the Allen-Cahn 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 semi-discretized 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 Allen-Cahn 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 first-order stabilized semi-implicit scheme of the Allen-Cahn 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 Allen-Cahn and Cahn-Hilliard 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 Cahn-Hilliard 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 second-order stabilized schemes for the Allen-Cahn 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 2-dimensional and 3-dimensional tensor-product 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 step-size. 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 non-negative 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 Allen-Cahn 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 double-well 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 (SL-BDF2) calculate iteratively, using

\hb@xt@.01(2.4)

where and are two non-negative 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:cs11-1), 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 SL-BDF2 scheme is stable with any non-negative including , if

\hb@xt@.01(2.14)

If one takes time step size even smaller,

\hb@xt@.01(2.15)

then the SL-BDF2 scheme is stable with any any non-negative and , including the case .

On the other hand side, if we take

\hb@xt@.01(2.16)

then the SL-BDF2 scheme is unconditional stable for any .

2.2 The stabilized linear Crank-Nicolson scheme.

Suppose and are given, our stabilized linear Crank-Nicolson scheme (SL-CN) calculate iteratively, using

\hb@xt@.01(2.17)

where and are two non-negative 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 SL-CN scheme is unconditional stable for any .

On the other hand, by using the inequality , it is easy to prove that when , the SL-CN scheme (LABEL:accn1) is stable for

\hb@xt@.01(2.28)
Remark 2.4

To make SL-BDF2 and SL-CN 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 phase-field Cahn-Hilliard–Navier-Stokes model for binary fluids has a fast convergence with respect to when the phenomenological mobility . We expect a similar scaling of in the Allen-Cahn model, which means is not necessary as large as . Furthermore, it is proved that the numerical interface for the Allen-Cahn 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 semi-implicit Fourier schemes, respectively, for the Cahn-Hilliard equation with double-well 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 semi-implicit stabilized scheme proposed by Xu and Tang [45] applied to the Cahn-Hilliard 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 Allen-Cahn 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 Cahn-Hilliard 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 Allen-Cahn 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 Allen-Cahn 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 non-negative constant.

Since the solution of Allen-Cahn 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 non-negative 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 Allen-Cahn 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:

  1. ;

  2. ;

  3. ;

  4. ;

  5. .

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