Linear, Second order and Unconditionally Energy stable schemes for The Viscous Cahn-Hilliard Equation with hyperbolic relaxation

5.1 MB·

by  Xiaofeng Yang

LINEAR, SECOND ORDER AND UNCONDITIONALLY ENERGY STABLE
SCHEMES FOR THE VISCOUS CAHN-HILLIARD EQUATION WITH
HYPERBOLIC RELAXATION
XIAOFENGYANG† ANDJIAZHAO‡
7
1
0 Abstract. In this paper, we consider the numerical approximations for the fourth order viscous
2 Cahn-Hilliard equation with the hyperbolic relaxation. The main challenge in solving such a dif-
fusive system numerically is how to develop high order temporal discretization for the hyperbolic
n
andnonlineartermsthatallowslargetimestepwhilepreservingtheunconditionalenergystability,
a
J i.e., the energy dissipative structure at the time-discrete level. We resolve this issue by developing
twosecondordertimemarchingschemesusingtherecentlydeveloped“InvariantEnergyQuadrati-
9
zation”approachwhereallnonlineartermsaretreatedsemi-explicitly. Ineachtimestep,oneonly

needstosolvealinearsystemthatissymmetricpositivedefinite. Werigorouslyproveallproposed
]
A schemes are unconditionally energy stable, which is further verified by time step refinement test
N numerically. Moreover,anumberof2Dand3Dnumericalsimulationsarepresentedtodemonstrate
thestability,accuracyandefficiencyoftheproposedschemes.
.
h
t
a
m
1. Introduction
[

The viscous Cahn-Hilliard (VCH) system and/or its perturbed verison with a hyperbolic re-
1
v laxation (HR) effect, had been well-studied recently, where the topics are mainly focused on the
6 well-posedness, sharp interface limit or global attractor, etc., see [2,3,6,9–11,14,15,17–19,23–30,
6 33,35,40,43–45,54] and the references therein. Formally, the governing equation of the VCH-HR
0
system is slightly different from the most classical Cahn-Hilliard (CH) equation, that is proposed
2
0 by Cahn and Hillard in [5], by incorporating a strong damping term (or called “viscosity”) and a
. hyperbolic relaxation term (or called “inertia”). The viscosity term is firstly proposed by Novick-
1
0 Cohen in [43] to introduce an additional regularity and some parabolic smoothing effects. Its
7 system can be viewed as a singular limit of the phase field equations for phase transitions (cf. [18]).
1 The hyperbolic relaxation term was proposed by Galenko et.al.in [23–28,35], in order to describe
:
v strongly non-equilibrium decomposition generated by rapid solidification under supercooling into
Xi the spinodal region occurring in certain materials (e.g., glasses). Since the VCH-HR system com-
bines the hyperbolic relaxation and the viscosity together, it is mathematically more complicated
r
a and tractable than the CH or VCH system (cf. [33,45,67]).
Before developing efficient numerical schemes to solve the VCH-HR system, we notice that its
reduced version, the classical CH equation is now widely applied to model the interfacial dynamics
in various scientific fields (cf. [5,7,16,34,36,39,41,55] and the references therein). The CH equation
and its analogous counterpart, the Allen-Cahn equation, are both categorized as representative
equations of phase field type models. From the numerical point of view, although all unknown
variables in phase field models are continuous and smooth, the induced systems are still very stiff
Keywordsandphrases. Phase-field,Linear,Cahn-Hilliard,EnergyStability,Secondorder,InvariantEnergyQuadra-
tization.
†Corresponding author, Department of Mathematics, University of South Carolina, Columbia, SC 29208; Email:
[email protected].
‡DepartmentofMathematics,UniversityofNorthCarolinaatChapelHill,27599;Email: [email protected].
1
2 XIAOFENGYANGANDJIAZHAO
where the stiffness is induced by the order parameter representing the thickness of the interface.
Therefore, when solving such stiff systems, it is desirable to establish efficient numerical schemes
that can verify the so called “energy stable” property at the discrete level irrespectively of the
coarseness of the discretization. In what follows, those algorithms will be called unconditionally
energy stable or thermodynamically consistent. Schemes with this property is specially preferred
sinceitisnotonlycriticalforthenumericalschemetocapturethecorrectlongtimedynamicsofthe
system,butalsoprovidessufficientflexibilityfordealingwiththestiffnessissue. Inspiteofthis,itis
necessary to point out a basic fact that larger time step will definitely induce larger computational
errors. In other words, the schemes with unconditional energy stability can allow unconditionally
large time step only for the sake of the stability concern. In practice, the controllable accuracy is
one of the most important factors to measure whether a scheme is reliable or not. Therefore, if
one attempts to use the time step as large as possible while maintaining desirable accuracy, the
only possible solution is to develop more accurate schemes, e.g., the second order schemes with
unconditionally energy stability, that is the main focus of this paper.
It is remarkable that, unlike the enormous numerical studies on the classical CH system, almost
all research related to the VCH or VCH-HR system had been focused on the theoretical PDE
analysis with very few numerical analysis or algorithm design. More precisely, to the best of
the authors’ knowledge, no schemes have been claimed to posses the following three properties,
namely, easy-to-implement, unconditionally energy stability and second order accuracy for the
VCH-HR model. This is because two additional numerical difficulties, the discretization of the
viscous effect as well as the hyperbolic inertia, emerge except the regular stiffness issue induced
by the nonlinear double well potential. At the very least, even for the reduced version, the CH
system, the algorithm design is still challenging. It can be seen clearly from the following fact
that some severe stability conditions on the time step will occur if the nonlinear term is discretized
in some normal ways like fully implicit or explicit type approaches. Such a time step constraint
can cause very high the computational cost in practice [4,21,49]. Many efforts (primarily for CH
system) had been done in order to remove this constraint and two commonly used techniques were
developed, namely, the nonlinear convex splitting approach [20,31,54,56], and the linear stabilized
approach [8,37,38,42,46–53,58,59,61,63,65,66]. The convex splitting approach is unconditionally
energystable,butitproducesthenonlinearscheme,thustheimplementationiscomplicatedandthe
computational cost might be high. The linear stabilized approach is linear so it is efficient and very
easy to implement. But, its stability requests a special property (generalized maximum principle)
satisfied by the classical PDE solution and the numerical solution, that is very hard to prove in
general. Moreover,itisdifficulttoextendtosecond-orderwhilepreservingtheunconditionalenergy
stability (cf. [49]).
Therefore, in order to develop some more efficient and accurate time marching schemes for
solving the VCH-HR equation, we use the Invariant Energy Quadratization (IEQ) approach, which
has been successfully applied to solve various phase field type models, see the authors’ recent work
in [57,60,62,64]. Its idea, that is very simple but quite different from those traditional methods
like implicit, explicit, nonlinear splitting, or other various tricky Taylor expansions to discretize
the nonlinear potentials, is to make the free energy quadratic. To be more specific, the free energy
potential is transformed into the quadratic form forcefully via the change of variables. Then, upon
a simple reformulation, all nonlinear terms are treated by the semi-explicit way, which in turn
yields a linear system. We develop two second order schemes, in which, one is based on the Crank-
Nicolson, and the other is based on the Adam-Bashforth (BDF2). These schemes are not only
easy-to-implement (linear system), but also unconditionally energy stable (with a discrete energy
dissipation law). Moreover, we show that the linear operator of all schemes are symmetric positive
NUMERICAL SCHEMES FOR THE VCH-HR MODEL 3
definite,sothatonecansolveitusingthewell-developedfastmatrixsolversefficiently(CGorother
Krylov subspace methods). Through various 2D and 3D numerical simulations, we demonstrate
stability and accuracy of the proposed schemes.
Therestofthepaperisorganizedasfollows. InSection2, wepresentthewholesystemandshow
the energy law in the continuous level. In Section 3, we develop the numerical schemes and prove
their unconditional stabilities. In Section 4, we present various 2D and 3D numerical experiments
tovalidatetheaccuracyandefficiencyoftheproposednumericalschemes. Finally,someconcluding
remarks are presented in Section 5.
2. Model Equations
Nowwegiveabriefdescriptionforthemodelequations. Weconsiderabinaryalloyinabounded
domain Ω Rd,d=2,3 with ∂Ω Lipschitz continuous, and define φ(x,t) as relative concentration
∈
of one material component, J the diffusion flux, then the balance law for concentration gives
(2.1) φ + J =0.
t
∇·
In order to describe the evolution for φ, one need to introduce a constitutive assumption on J, i.e.,
δE
(2.2) αJ +J = ( +βφ ).
t t
−∇ δφ
Here E(φ) is the total free energy that takes the form as
(cid:15)2
(2.3) E(φ)= φ2+F(φ) dx,
2|∇ |
ZΩ(cid:16) (cid:17)
where α 0 is the relaxation parameter, β 0 is the viscosity parameter, and (cid:15) is the positive
≥ ≥
constant that measures the interfacial width.
For the choice of nonlinear potential F(φ), we can choose either (i) double well (Ginzburg-
Landau) potential where
(2.4) F(φ)=φ2(φ 1)2;
−
or (ii) Flory-Huggins potential (cf. [45]) where
(2.5) F(φ)=(1 φ)ln(1 φ)+φlnφ+θφ(1 φ),θ >0.
− − −
By combining (2.1) and (2.2), the governing PDE reads as follows,
(2.6) αφ +φ =λ∆µ,
tt t
(2.7) µ= (cid:15)2∆φ+f(φ)+βφ ,
t
−
where f(φ)=F (φ), i.e., f(φ)=2φ(φ 1)(2φ 1) for double well potential and f(φ)=ln( φ )+
0 − − 1 φ
θ(1 2φ) for Flory-Huggins potential. The boundary conditions can be −
−
(2.8) (i) all variables are periodic;or (ii)∂ φ =∂ µ =0,
n ∂Ω n ∂Ω
| |
where n is the unit outward normal on ∂Ω.
When α = β = 0, the system degenerates to the standard CH system that conserves the local
mass density. When α=0, the volume conservation will only hold provided φ (0,x)dx=0. To
6 Ω t
see this, by taking the L2 inner product of (2.6) with 1, one can obtain directly
R
d
(2.9) α φ (t,x)dx+ φ (t,x)dx=0.
t t
dt
ZΩ ZΩ
4 XIAOFENGYANGANDJIAZHAO
This actually is an ODE system for time, and its solution is
t
(2.10) φ (t,x)dx=exp(− ) φ (0,x)dx.
t t
α
ZΩ ZΩ
Therefore, by setting φ (0,x)dx=0, we obtain
Ω t
R
(2.11) φ (t,x)dx= φ (t,x)dx=0.
t tt
ZΩ ZΩ
Hereandafter, foranyfunctionsg (x),g (x) L2(Ω),wedenotetheinnerproductandL2 norm
1 2
∈
as
(2.12) (g ,g )= g (x)g (x)dx, g = g (x)2dx.
1 2 1 2 1 1
k k | |
ZΩ ZΩ
We define the inverse Laplace operator u ( L2(Ω)) with udx=0 v :=∆ 1u (v H1(Ω)) by
∈ Ω 7→ − ∈
R
∆v =u, udx=0,
(2.13)
 ZΩ
with the boundary conditions either (i) v is periodic,or (ii) ∂ v =0.
 n ∂Ω
|
Lemma 2.1. The model system (2.6)-(2.7) satisfies an energy dissipation law as,
d (cid:15)2 α
(2.14) φ2+F(φ)+ ∆ 1φ 2 dx= ∆ 1φ 2 β φ 2 0.
− t − t t
dt 2|∇ | 2|∇ | −k∇ k − k k ≤
ZΩ(cid:16) (cid:17)
Proof. We introduce a new variable ψ =φ . Since ψdx= ψ dx=0, using the operator ∆ 1,
t Ω Ω t −
we rewrite the system (2.6)-(2.7) as follows,
R R
(2.15) α∆ 1ψ +∆ 1ψ = (cid:15)2∆φ+f(φ)+βφ .
− t − t
−
By taking the L2 inner product of (2.15) with φ , we have
t
d (cid:15)2
(2.16) α(∆ 1ψ ,ψ)+(∆ 1ψ,ψ) β φ 2 = φ2+F(φ) dx.
− t − t
− k k dt 2|∇ |
ZΩ(cid:16) (cid:17)
Since ψdx=0, we can find another auxillary variable p such that p=∆ 1ψ, i.e.,
Ω −
R
(2.17) ∆p=ψ, ψdx=0,
ZΩ
with the boundary condition specified in (2.13). By taking the L2 inner product of (2.17) with p,
that is
(2.18) (p,ψ)= p 2 =(∆ 1ψ,ψ).
−
−k∇ k
We differentiate (2.17) with time t to obtain
(2.19) ∆p =ψ .
t t
BytakingtheL2 innerproductof (2.19)withp, weobtain(ψ ,p)=(∆p ,p)= 1d p 2. Hence,
t t −2 tk∇ k
we derive
α
(2.20) α(∆ 1ψ ,ψ)=α(ψ ,∆ 1ψ)=α(ψ ,p)= d p 2.
− t t − t t
−2 k∇ k
By combining (2.16)-(2.18)-(2.20), we obtain
d (cid:15)2 α
(2.21) φ2+F(φ)+ p2 dx= p 2 β φ 2 0.
t
dt 2|∇ | 2|∇ | −k∇ k − k k ≤
ZΩ(cid:16) (cid:17)
This will ensure the modified energy of the VCH-HR system (2.6)-(2.7) nonincreasing in time. (cid:3)
NUMERICAL SCHEMES FOR THE VCH-HR MODEL 5
3. Numerical Schemes
We now construct two semi-discrete time marching numerical schemes for solving the model sys-
tem(2.6)-(2.7)-(2.8)andprovetheirenergystabilitiesbasedontheInvariantEnergyQuadratization
(IEQ) approach. The key point of the IEQ method is to make the nonlinear potential quadratic via
the change of auxiliary variables. It is feasible since we notice that the nonlinear potential F(φ)
is always bounded from below, in either the double well form (for the Ginzberg-Landau potential)
or logarithmic form (for the Flory-Huggins potential). Thus, in general, we could rewrite the free
energy functional F(φ) into the following equivalent form
(3.1) F(φ)=(F(φ)+B) B,
−
where B is some constant to ensure F(x)+B 0, x R, and define an auxilliary function U as
≥ ∀ ∈
(3.2) U = F(φ)+B.
Thus the total energy of (2.3) turns to a nepw form
(cid:15)2
(3.3) E(φ,U)= φ2+U2 B dx.
2|∇ | −
ZΩ
Then we obtain an equivalent PDE syst(cid:0)em by taking the tim(cid:1) e derivative for the new variable U:
(3.4) αψ +ψ =∆µ,
t
(3.5) µ= (cid:15)2∆φ+UH +βφ ,
t
−
1
(3.6) U = Hφ ,
t t
2
(3.7) ψ =φ ,
t
where
f(φ)
(3.8) H(φ)= , f(φ)=F (φ).
0
F(φ)+B
The boundary conditions for the new system are still (2.8) since the equation (3.6) for the new
p
variable U is simply an ODE with time. The initial conditions read as
(3.9) φ =φ ,ψ =0,
(t=0) 0 (t=0)
| |
(3.10) U = F(φ )+B,
(t=0) 0
|
where we simply set the initial profile of ψ to bepzero pointwise.
It is clear that the new transformed system (3.4)-(3.7) still retains a similar energy dissipative
law. By applying the inverse Laplace operator ∆ 1 to (3.4), taking the L2 inner product of it with
−
φ , of (3.6) with U, using (2.18) and (2.20), and summing them up, we can obtain the energy
t
−
dissipation law of the new system (3.4)-(3.5) as
d (cid:15)2 α
(3.11) φ2+U2+ ∆ 1ψ 2 dx= ∆ 1φ 2 β ψ 2 0.
− − t
dt 2|∇ | 2|∇ | −k∇ k − k k ≤
ZΩ(cid:16) (cid:17)
Remark 3.1. We emphasize that the new transformed system (3.4)-(3.7) is exactly equivalent to
the original system (2.6)-(2.7), since (3.2) can be easily obtained by integrating (3.6) with respect to
the time. For the time-continuous case, the potentials in the new free energy (3.3) are the same as
the Lyapunov functional in the original free energy of (2.3), and the new energy law (3.11) for the
transformed system is also the same as the energy law (2.14) for the original system as well. We
will develop unconditionally energy stable numerical schemes for time stepping of the transformed
6 XIAOFENGYANGANDJIAZHAO
system (3.4)-(3.7), and the proposed schemes should formally follow the new energy dissipation law
(3.11) in the discrete sense, instead of the energy law for the originated system (2.14).
Remark3.2. IfF(φ)=φ2(φ 1)2, weletB =0, thusH(φ)=2φ 1. Atthistime, theIEQmethod
− −
is exactly the same as the so-called Lagrange multiplier method developed in [32]. We remark that
the Lagrange multiplier method in [32] only works for the fourth order polynomial potential (φ4).
This is because the term φ3 (the first order derivative of φ4) can be decomposed into a multiplication
of two factors: λ(φ)φ, where λ(φ)=φ2. In [32], this Lagrange multiplier term λ(φ) is then defined
as the new auxiliary variable U. However, for other type potentials, e.g., the F-H potential, if
following the Lagrange method in [32], a new variable U will be defined as λ(φ)= 1ln( φ ), which
φ 1 φ
is unworkable for algorithms design. −
Remark3.3. IfF(φ)=(1 φ)log(1 φ)+φlogφ+θφ(1 φ),followingtheworkin[12],weregularize
− − −
the logarithmic bulk potential by a C2 piecewise function. More precisely, for any 0 < σ 1, the
(cid:28)
regularized free energy is
φlnφ+ (1−φ)2 +(1 φ)lnσ σ +θφ(1 φ), if φ 1 σ,
2σ − − 2 − ≥ −
(3.12)F(φ)= φlnφ+(1 φ)ln(1 φ)+θφ(1 φ), if σ φ 1 σ,

− − − ≤ ≤ −
(1 φ)ln(1 φ)+ φ2 +φlnσ σ +θφ(1 φ), if φ σ.
b − − 2σ − 2 − ≤
For convenience, we consider the problem formulated with the substitute F(φ), but omit the
in the notation. Now the domain for the regularized functional F(φ) is now ( , ). Hence, we
−∞ ∞
do not need to worry about that any small fluctuation near the domain bounbdary ( 1,1) of theb
−
numerical solution can cause the overflow. In [12], the authors proved the error bound between the
regularized PDE and the original PDE is controlled by σ up to a constant. For this case, we simply
take B =1 that can ensure F(x)+B >0, x R.
∀ ∈
The time marching numerical schemes are developed to solve the new transformed system (3.4)-
b
(3.7). The proof of the unconditional stability of the schemes follows the similar lines as in the
derivation of the energy law (3.11). Let δt > 0 denote the time step size and set tn = nδt for
0 n N with the ending time T =Nδt.
≤ ≤
3.1. Crank-NicolsonScheme. WefirstdevelopasecondorderversonthatisbasedontheCrank-
Nicolson, that reads as follows.
Scheme 1. Compute U1 and φ1 by assuming U 1 = U0 and φ 1 = φ0 for the initial step. Having
− −
computed (φn, Un) and (φn 1,Un 1), with n 1, we update φn+1 and Un+1 as follows:
− −
≥
ψn+1 ψn ψn+1+ψn
(3.13) α − + =∆µn+1,
δt 2
φn+1+φn Un+1+Un φn+1 φn
(3.14) µn+1 = (cid:15)2∆ + H?+β − ,
− 2 2 δt
1
(3.15) Un+1 Un = H?(φn+1 φn),
− 2 −
ψn+1+ψn φn+1 φn
(3.16) = − ,
2 δt
where
f(φ?) 3 1
(3.17) H? = , φ? = φn φn 1.
−
F(φ?)+B 2 − 2
p
NUMERICAL SCHEMES FOR THE VCH-HR MODEL 7
The boundary conditions are either
(3.18) (i) φn+1,µn+1 are periodic;or (ii)∂ φn+1 = µn+1 n =0.
n ∂Ω ∂Ω
| ∇ · |
Since the nonlinear coefficient H of the new variables U are treated explicitly, we can rewrite
the equations (3.15) and (3.16) as follows:
H?
Un+1 = φn+1+gn,
2 1
(3.19)
 2
ψn+1 = δtφn+1+g2n,
where g1n =(Un−H2?φn), g2n =(−δ2tφn−ψn). Thus (3.13)-(3.14) can be rewritten as the following
linear system
(3.20) αφn+1 = ∆µn+1+gn,
3
(3.21) µn+1 = P (φn+1)+gn,
1 4
b
where
α 1 2
α=( + ) ,
δt 2 δt

(cid:15)2 1 β
(3.22) Pb1(φn+1)=−2∆φn+1+ 4H?H?φn+1+ δtφn+1,
 α 1 α 1
gn = ( + )gn+( )ψn,
 3 − δt 2 2 δt − 2
gn = −(cid:15)2∆φn+ 1H?(gn+Un) βφn.
 4 2 2 1 − δt


Therefore, we can solve φn+1 and µn+1 directly from (3.20) and (3.21). Once we obtain φn+1,

the ψn+1,Un+1 are automatically given in (3.19). Furthermore, for any ψ that enjoys the same
boundary condition as φ in (3.18), we have
(cid:15)2 1 β
(3.23) (P (φ),ψ)= ( φ, ψ)+ (H?φ,H?ψ)+ (φ,ψ).
1
2 ∇ ∇ 4 δt
Therefore,thelinearoperatorP (φ)issymmetric(self-adjoint). Moreover,foranyφwith φdx=
1 Ω
0, we have
R
(cid:15)2 1 β
(3.24) (P (φ),φ)= φ 2+ Hnφ 2+ φ 2 0,
1
2k∇ k 4k k δtk k ≥
where “=” is valid if and only if φ 0.
≡
We first show the well-posedness of the linear system (3.13)-(3.16) (or (3.20)-(3.21)) as follows.
Theorem 3.1. The linear system (3.20)-(3.21) admits a unique solution in H1(Ω), and the linear
operator is symmetric positive definite.
Proof. From (3.13), by taking the L2 inner product with 1 and notice ψ0 =0, we derive
α 1 α 1
(3.25) ( + ) ψn+1dx=( ) ψndx=0.
δt 2 δt − 2
ZΩ ZΩ
From (3.16), we have
(3.26) φn+1dx= φndx= = φ0dx.
···
ZΩ ZΩ ZΩ
8 XIAOFENGYANGANDJIAZHAO
Let V = 1 φ0dx, V = 1 µn+1dx, and we define
φ Ω Ω µ Ω Ω
| | | |
(3.27) R φn+R1 =φn+1 V , µn+1 =µn+1 V .
φ µ
− −
Thus, from (3.20)-(3.21), (φn+1,µn+1) are the solutions for the following equations with unknowns
(φ,w), b b
(3.28) b b αφ ∆w =fn,
−
(3.29) w+V P (φ)=gn,
µ 1
−
where fn =gn αV , fndx=0, gn =bgn+ 1H?H?α + βα , φdx=0 and wdx=0.
3 − φ Ω 4 4 0 δt 0 Ω Ω
Applying ∆ 1 to (3.28) and using (3.29), we obtain
− − R R R
(3.30) b α∆ 1φ+P (φ) V = ∆ 1fn gn.
− 1 µ −
− − − −
Let us express the above linear system (3.30) as Aφ=b,
(i). Foranyφ andφ inHb1(Ω)satisfytheboundaryconditions(2.8)and φ dx= φ dx=0,
1 2 Ω 1 Ω 2
using integration by parts, we derive
R R
(A(φ1),φ2)= α(∆−1φ1,φ2)+(P1(φ1),φ2)
−
(3.31) C ( ∆ 1φ ∆ 1φ + φ φ + φ φ )
1 − 1 − 2 1 2 1 2
≤ k∇ kk∇ k k∇ kk∇ k k kk k
b
C φ φ .
2 1 H1 2 H1
≤ k k k k
Therefore, the bilinear form (A(φ1),φ2) is bounded φ1,φ2 H1(Ω).
∀ ∈
(ii). For any φ H1(Ω), it is easy to derive that, ,
∈
(cid:15)2 1 β
(3.32) (A(φ),φ)=αk∇∆−1φk2+ 2k∇φk2+ 4kH?φk2+ δtkφk2 ≥C3kφk2H1,
for φdx=0 from Poincare inequality. Thus the bilinear form (A(φ),ψ) is coercive.
Ω b
Then from the Lax-Milgram theorem, we conclude the linear system (3.30) admits a unique
R
solution in H1(Ω).
For any φ ,φ with φ dx=0 and φ dx=0, we can easily derive
1 2 Ω 1 Ω 2
(3.33) R (AφR1,φ2)=(φ1,Aφ2).
Thus A is self-adjoint. Meanwhile, from (3.32), we derive (Aφ,φ) 0, where “=” is valid if only if
≥
φ=0. This concludes the linear operator A is positive definite. (cid:3)
The energy stability of the scheme (3.13)-(3.16) (or (3.20)-(3.21)) is presented as follows.
Theorem 3.2. The scheme (3.13)-(3.16) (or (3.20)-(3.21)) is unconditionally energy stable satis-
fying the following discrete energy dissipation law,
1 (pn+1+pn) 2 φn+1 φn 2
(3.34) (En+1 En )= ∇ β − 0,
δt cn2 − cn2 − 2 − δt ≤
(cid:13) (cid:13) (cid:13) (cid:13)
where (cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:15)2 α
(3.35) E = φ 2+ U 2+ p 2 B Ω.
cn2
2k∇ k k k 2k∇ k − | |
Proof. First, we combine (3.13) and (3.14) together and apply the ∆ 1 to obtain
−
α ψn+1+ψn
∆ 1(ψn+1 ψn)+∆ 1
− −
δt − 2
(3.36)
φn+1+φn Un+1+Un φn+1 φn
= (cid:15)2∆ + H?+β − .
− 2 2 δt
NUMERICAL SCHEMES FOR THE VCH-HR MODEL 9
Secondly, by taking the L2 inner product of (3.36) with φn+1 φn, we obtain
−
α 1
(∆ 1(ψn+1 ψn),φn+1 φn)+ (∆ 1(ψn+1+ψn),φn+1 φn)
− −
δt − − 2 −
(3.37)
(cid:15)2 Un+1+Un β
= ( φn+1 2 φn 2)+( H?,φn+1 φn)+ φn+1 φn 2.
2 k∇ k −k k 2 − δtk − k
Thirdly, by taking the L2 inner product of (3.15) with (Un+1+Un), we obtain
−
1
(3.38) ( Un+1 2 Un 2)= ( H?(φn+1 φn),Un+1+Un).
− k k −k k − 2 −
Fourthly, define pn+1 =∆ 1ψn+1. By subtracting with the n-step, we obtain
−
(3.39) ∆(pn+1 pn)=ψn+1 ψn.
− −
From (3.16) and (3.39), we derive
α α
(∆ 1(ψn+1 ψn),φn+1 φn)= (pn+1 pn,ψn+1+ψn)
−
δt − − 2 −
α
(3.40) = (pn+1 pn,∆(pn+1+pn))
2 −
α α
= ( pn+1 2 pn 2),
− 2k∇ k − 2k∇ k
and
1 δt
(∆ 1(ψn+1+ψn),φn+1 φn)= (pn+1+pn,ψn+1+ψn)
−
2 − 4
δt
(3.41) = (pn+1+pn,∆(pn+1+pn))
4
δt
= (pn+1+pn) 2.
−4k∇ k
Finally, by combining (3.37), (3.38), (3.40) and (3.41), we obtain
(cid:15)2 α
( φn+1 2 φn 2)+ Un+1 2 Un 2+ ( pn+1 2 pn 2)
2 k∇ k −k∇ k k k −k k 2 k∇ k −k∇ k
(3.42)
δt β
= (pn+1+pn) 2 φn+1 φn 2,
−4k∇ k − δtk − k
that concludes the theorem. (cid:3)
Remark 3.4. The proposed schemes follow the new energy dissipation law (3.11) formally in-
stead of the energy law for the originated system (2.14). In the time-discrete case, the energy
E(φn+1,Un+1) (defined in (3.35)) can be rewritten as a first order approximation to the Lyapunov
functionals in E(φn+1) (defined in (2.14)), that can be observed from the following facts, heuristi-
cally. Assuming the case for double well potential, from (3.15), we have
(3.43) Un+1 ( F(φn+1)+B)=Un ( F(φn)+B)+R ,
n+1
− −
where Rn+1 = O((φn+1 φpn)(φn+1 2φn +φn−1)).pSince Rk = O(δt3) for 0 k n+1 and
− − ≤ ≤
U0 =( F(φ0)+B), by mathematical induction we can easily get
(3.44) p Un+1 = F(φn+1)+B+O(δt2).
p
3.2. Adam-Bashforth Scheme. We further develop another second order scheme based on the
Adam-Bashforth backward differentiation formula (BDF2), as follows.
10 XIAOFENGYANGANDJIAZHAO
Scheme 2. Compute U1 and φ1 by assuming U 1 = U0 and φ 1 = φ0 for the initial step. Having
− −
computed (φn, Un) and (φn 1,Un 1), with n 1, we solve φn+1, Un+1 as follows:
− −
≥
3ψn+1 4ψn+ψn 1
(3.45) α − − +ψn+1 =∆µn+1,
2δt
3φn+1 4φn+φn 1
(3.46) µn+1 = (cid:15)2∆φn+1+Un+1H†+β − − ,
− 2δt
1
(3.47) 3Un+1 4Un+Un 1 = H (3φn+1 4φn+φn 1),
− † −
− 2 −
3φn+1 4φn+φn 1
(3.48) ψn+1 = − −
2δt
where
f(φ )
(3.49) H = † ,φ =2φn φn 1.
† † −
F(φ )+B −
†
The boundary conditions are
p
(3.50) (i) φn+1,µn+1 are periodic;or (ii)∂ φn+1 = µn+1 n =0.
n ∂Ω ∂Ω
| ∇ · |
Similar to the Crank-Nicolson scheme, we can rewrite the equations (3.47) and (3.48) as follows:
H
Un+1 = †φn+1+hn,
2 1
(3.51)

3
ψn+1 = 2δtφn+1+hn2,
where hn1 =(U±− H2†φ±), hn2 = 23δtφ± with S± = 4Sn−3Sn−1 for any variable S. Thus (3.45)-(3.46)
can be rewritten as the following linear system
(3.52) αφn+1 =∆µn+1+hn,
3
(3.53) µn+1 =P (φn+1)+hn,
2 4
e
where
1 3β
P (φn+1)= (cid:15)2∆φn+1+ H H φn+1+ φn+1,
2 † †
− 2 2δt
 3α 3α
hn = ( +1)hn+ ψ ,
(3.54)  3 − 2δt 2 2δt ±
 1 3β
hn4 = 2H†hn1 − 2δtφ±,
3α 3

α=( +1) .
 2δt 2δt


Actually, we can solve φn+1 and µn+1 directly from (3.52) and (3.53). Once we obtain φn+1, the

e
ψn+1,Un+1 is automatically given in (3.51). Furthermore, we notice
1 3β
(3.55) (P (φ),ψ)=(cid:15)2( φ, ψ)+ (H φ,H ψ)+ (φ,ψ),
2 † †
∇ ∇ 2 2δt
if ψ enjoys the same boundary condition as φ in (3.50). Therefore, the linear operator P (φ) is
2
symmetric (self-adjoint). Moreover, for any φ with φdx=0, we have
Ω
1
(3.56) (P (φ),φ)=(cid:15)2 φ 2+R H φ 2 0,
2 †
k∇ k 2k k ≥
where “=” is valid if and only if φ 0.
≡