# Numerical Approximation of the Space Fractional Cahn-Hilliard Equation

Authors: Zhifeng Weng; Langyang Huang; Rong Wu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In this paper, a second-order accurate (in time) energy stable Fourier spectral scheme for the fractional-in-space Cahn-Hilliard (CH) equation is considered. The time is discretized by the implicit backward differentiation formula (BDF), along with a linear stabilized term which represents a second-order Douglas-Dupont-type regularization. The semidiscrete schemes are shown to be energy stable and to be mass conservative. Then we further use Fourier-spectral methods to discretize the space. Some numerical examples are included to testify the effectiveness of our proposed method. In addition, it shows that the fractional order controls the thickness and the lifetime of the interface, which is typically diffusive in integer order case.


---

## Body

## 1. Introduction

The Cahn-Hilliard model [1, 2] was originated from a phase-separation model in a binary alloy. The Cahn-Hilliard model has been widely used in many complicated moving interface problems in materials science and fluid dynamics through a phase-field approach [3–5]. Hence, it is important to study the physical properties of Cahn-Hilliard equations by accurate and efficient numerical algorithms. The numerical analysis of Cahn-Hilliard equations and the references therein have been studied in the past [3, 6, 7].However, in the original formulation of the physical model [1], nonlocal interactions were parts of CH model, but in the subsequent mathematical literatures CH model was approximated by local one by assuming slow spatial variations. In recent years, fractional phase-field models have attracted the attention of many researchers, and a relevant number of works [8–12] have been devoted to the mathematical/numerical analysis of the space (time) fractional phase-field models. The space fractional phase-field models can provide a more accurate description of the formation and phase change mechanism by considering the long-range interactions among particles. So it is necessary to obtain high accuracy numerical solutions and design numerical methods for the space fractional phase-field models that satisfy the energy decays with respect to time.In this paper, we focus on the following fractional-in-space Cahn-Hilliard (CH) equation [9]:(1)∂u∂t=Δϵ2-Δα/2u+fu,x,t∈Ω×0,T,ux,0=u0x,x∈Ω,where Ω is a bounded regular domain in Rd(d=1,2,3), and T>0 denotes a finite time. For the above equation, we impose the periodic or following homogeneous Neumann boundary conditions(2)n·∇u=0,n·∇fu+ϵ2-Δα/2u=0,on Ω. Here n is the unit outward normal vector of the domain Ω. The function u denotes the concentration for one of the two metallic components of the alloy, and ϵ is the interfacial energy in a phase-separation phenomenon. The nonlinear function f(u)=F′(u) and F(u)=1/4(u2-1)2 is a given energy potential, taking its global minimum value at u=-1 and u=1. In other words, the values -1 and 1 of the function u prevail in Ω whereas the transition between them forms a thin interface layer.The fractional CH equation (1) can be viewed as an H-1-gradient flow of the following fractional analogue version of Ginzburg-Landau free energy functional [10]:(3)Eαu=∫Ωϵ22-Δα/4u2+Fudx.The same as the standard CH equation [1], one intrinsic property of the fractional-in-space CH equation is energy-decay property in sense of(4)ddtEαu=-∫Ω∇ϵ-Δα/2u+fu2dx≤0.Another intrinsic property of the fractional-in-space CH equation, as well as the standard CH equation, is the law of conservation of mass (5)ddtMu=ddt∫Ωux,tdx=0.Comparing with the well-known classic CH equation [1], the fractional-in-space CH equation is equipped with a nonlocal diffusion operator and can describe more practical phenomena for modeling phase transitions of microstructures in materials. On the other hand, it evidently brings more computational costs in numerical simulations; thus efficient and accurate numerical methods are highly desired. Recently, Bosch and Stoll [8] used Fourier spectral method for fractional vector-valued CH model combined with an implicit first-order convex splitting scheme. Numerical tests showed the superiority of the proposed fractional CH inpainting approach over its nonfractional version. Ainsworth and Mao [9] proposed a first-order semi-implicit Fourier-Galerkin scheme for solving the 2D model (1) and detailed theoretical research and simulation studies were given. Weng [11] et al. proposed an unconditionally energy stable nonlinear Fourier spectral scheme for model (1) by introducing a stabilizing term.In our previous work [11], we have proposed second-order unconditionally energy stable Fourier spectral scheme for fractional-in-space Cahn-Hilliard equation using the Crank-Nicolson methodology and nonlinear stabilized term. In this paper, new schemes are proposed combining the implicit backward differentiation formula and stabilized term. Mainly, the new schemes are three-level schemes, and the stabilized term is linear. The newly proposed BDF schemes are unconditionally energy stable second-order accurate in time. Then we develop a Fourier spectral scheme to effectively treat the significantly increased memory requirement and computational complexity for space fractional CH models. The energy stability and mass conservation will be proven in the semidiscretized case in time. Compared with the literature [11], our method is more efficient for solving fractional-in-space Cahn-Hilliard in numerical experiments.The rest of the paper is organized as follows. In Section2, we develop the numerical schemes and prove their unconditional energy stability and mass conservation in the semidiscretized case in time. In Section 3, we present some numerical simulations to demonstrate the accuracy and efficiency of the proposed method. Finally, some concluding remarks are presented in Section 4.


## 2. Numerical Schemes

In this section, an unconditionally energy stable and mass conservation of the second order BDF-type scheme in time will be proven. Then we will propose fully discrete scheme based on the implicit backward difference scheme for the time and Fourier spectral method for the space to solve the fractional equation (1).


### 2.1. An Unconditionally Energy Stable and Mass Conservation Semidiscrete Scheme

First, we use a finite difference method to discretize the time direction of the fractional CH equation (1). Let tm=mτ, m=0,1,…,M, where τ=T/M. The new second-order scheme can be written as(6)3um+1-4um+um-12τ=Δϵ2-Δα/2um+1+um+13-2um+um-1-AτΔ2um+1-um,where AτΔ2(um+1-um) is a stabilizing term with A>0, and u1 is approximated by the following unconditionally stable Crank-Nicolson scheme [11]:(7)u1-u0τ=Δϵ2-Δα/2u1/2+f1/2-12Δu1-u02u1/2,where u1/2=u0+u1/2 and f1/2=f0+f1/2.For this semidiscretized problem, we have the following energy stability estimate, which is a slight modification of the Theorem 4.1 for the standard CH equation in [13].Theorem 1.
Form≥1, define(8)Eαum+1,um=Eαum+1+14τum+1-um-12+12um+1-um02,where Eα(u)=ϵ2/2(-Δ)α/4u2+1/4uL44-1/2u2 and suppose A≥1/16. Then the semidiscrete scheme (6) is unconditionally energy stable, i.e., for all τ>0(9)Eαum+1,um≤Eαum,um-1.Here um-12≔(um,(-Δ)-1um) is inner product of H-1(Ω).Proof.
We use·2 to denote the L2 norm and (·,·) to the inner product in L2. Taking the L2 inner product for (6) with (-Δ)-1(um+1-um) yields(10)
3
u
m
+
1
-
4
u
m
+
u
m
-
1
2
τ
,
-
Δ
-
1
u
m
+
1
-
u
m
+
ϵ
2
-
Δ
1
+
α
/
2
u
m
+
1
,
-
Δ
-
1
u
m
+
1
-
u
m
+
-
Δ
u
m
+
1
3
,
-
Δ
-
1
u
m
+
1
-
u
m
+
Δ
2
u
m
-
u
m
-
1
,
-
Δ
-
1
u
m
+
1
-
u
m
+
A
τ
-
Δ
2
u
m
+
1
-
u
m
,
-
Δ
-
1
u
m
+
1
-
u
m
=
J
1
+
J
2
+
J
3
+
J
4
+
J
5
=
0
.Now we estimate the terms J1-J5, respectively. For the time difference term, by the definition of the norm um-12, it is easy to obtain (11)J1=12τ3um+1-4um+um-1,-Δ-1um+1-um=14τ6um+1-8um+2um-1,-Δ-1um+1-um=τ54um+1-umτ-12-14um-um-1τ-12+τ34um+1-2um+um-1τ2-12.
For the highest-order diffusion, using the identity(12)b,b-a=12b2-12a2+12b-a2,we can obtain (13)J2=ϵ2-Δ1+α/2um+1,-Δ-1um+1-um=ϵ22-Δα/4um+12--Δα/4um2+ϵ22-Δα/4um+1-um2.By virtue of the identity (14)b3,b-a=14b4-a4+14b2-a22+12bb-a2,the nonlinear term becomes (15)J3=um+13,um+1-um=14um+1L44-umL44+14um+12-um22+12um+1um+1-um2.By means of the identity(16)b,b-a=12b2-12a2+12b-a2,the concave diffusive term turns out to be (17)J4=-2um-um-1,um+1-um=-um,um+1-um-um-um-1,um+1-um=-12um+12-um2+12um+1-um2+um+1-2um+um-1,um+1-um-um+1-um2=-12um+12+12um+1-um2--12um2+12um-um-12+12um+1-2um+um-12-12um+1-um2.Finally, the stabilizing term becomes (18)J5=Aτ∇um+1-um2.
For anyσ>0 and the formula from the last row in J4, in virtue of the Cauchy-Schwarz inequality, we have(19)12um+1-um2≤12∇um+1-∇umum+1-um≤τ4σ∇um+1-∇um2+στ4um+1-umτ-12
Putting everything together, we can get(20)Eαum+1,um+τ1-σ4um+1-umτ-12+τA-14σ∇um+1-um2≤Eum,um-1.
Consequently, dropping some unnecessary terms, we can obtain(21)Eαum+1,um≤Eαum,um-1.This completes the proof.Remark. The function (22)Eαum+1,um=Eαum+1+14τum+1-um-12+12um+1-um02consists of three parts. The formula Eα(um+1) comes from a typical free energy functional Eα(u). The formula 1/4τum+1-um-12 and 1/2um+1-um02 comes from the time difference term J1 and the concave diffusive term J4, respectively.For the initialization scheme (7), we have the following stability result.Theorem 2 (see [11]).
The initialization scheme (7) has the energy-decay property(23)Eαu1≤Eαu0.Moreover, we have the following mass conservation law.Theorem 3.
For (1), the semidiscrete scheme (6) satisfies the mass conservation law: (24)Mum+1=Mum,m=0,1,2,…,M-1Proof.
For m=0, it is easy to haveM(u1)=M(u0) based on Theorem 3.2. in [11].
Form>0, rewrite the scheme (6) as (25)3um+1-4um+um-12τ=Δω,ω=ϵ2-Δα/2um+1+um+13-2um+um-1-AτΔum+1-um.
Then we have(26)3um+1-4um+um-12τ,1=Δω,1=-∇ω,∇1=0.Hence(27)3um+1-um,1=um-um-1,1.Then combined with the result (u0,1)=(u1,1) (please see Theorem 3.2. in [11]), it is easy to verify that (28)Mum+1=Mum.This completes the proof.


### 2.2. Spatial Discretization Using Spectral Methods

In this subsection, we will give the fully discrete scheme for model (1). We know that the fractional Laplacian (-Δ)α/2 can be defined in various ways. A review of these existing definitions can be seen in [14]. One important and effective definition for (-Δ)α/2 is defined using Fourier transform; see [15–17].Without loss of generality, we setΩ=[a,b]d and multi-index(29)μ=p1,p2,…,pd∣pi∈Zμ+=p=p1,p2,…,pd∣pi∈N.Then, for p∈μ, suppose the d-dimensional Laplacian (-Δ) has a complete set of orthonormal eigenfunctions φp corresponding to eigenvalues λp on the bounded region [a,b]d, i.e., (-Δ)φp=λpφp. Let (30)Uα=u=∑l=1d∑pl=-∞∞u^pφp,u^p=u,φp,∑l=1d∑pl=-∞∞u^pλpα/2<∞,.Then for any u∈Uα, it is easy to have(31)-Δα/2u=∑l=1d∑pl=-∞∞λpα/2u^pφp,where λp and φp will depend on the specified boundary conditions.For periodic boundary condition,(32)λp=∑l=1d2plπb-a2,φp=1b-ad/2∏l=1dexp⁡i2plπx-ab-a,p∈μ.Here, the imaginary unit is denoted by i.For homogeneous Neumann boundary condition,(33)λp=∑l=1nplπb-a2,φp=2b-ad/2∏l=1dcos⁡plπx-ab-a,p∈μ+.In this work, we use (31) to handle the fractional Laplacian (-Δ)α/2.After discretizing governing equation in time as described in the above subsection, now we will further give the fully discrete scheme based on Fourier spectral method.By applying the Fourier transform to both sides of (6), together with f(u)=u3-u, we have the following nonlinear scheme for the discrete p-th Fourier coefficients (34)3upm+1^-4upm^+upm-1^2τ=-ϵ2λp1+α/2upm+1^-λpupm+13^+2λpqupm^-λpqupm-1^-Aτλp2upm+1^-upm^.For the nonlinear case, we need use an iterative method to solve the resulting algebraic system. Thus the above nonlinear term is treated using the following fixed point iteration: forl=1,2,…,L find upm+1,l^ such that(35)upm+1,l^=13+2τϵ2λp1+α/2+2Aτ2λp24+4τλp+2Aτ2λp2upm^-1+2τλpupm-1^-2τλpupm+1,l-13^,with upm+1,0^=upm^.It also preserves the mass conservation and energy dissipation properties in the full discrete sense. The proofs are similar. We thus omit them for simplicity.


## 2.1. An Unconditionally Energy Stable and Mass Conservation Semidiscrete Scheme

First, we use a finite difference method to discretize the time direction of the fractional CH equation (1). Let tm=mτ, m=0,1,…,M, where τ=T/M. The new second-order scheme can be written as(6)3um+1-4um+um-12τ=Δϵ2-Δα/2um+1+um+13-2um+um-1-AτΔ2um+1-um,where AτΔ2(um+1-um) is a stabilizing term with A>0, and u1 is approximated by the following unconditionally stable Crank-Nicolson scheme [11]:(7)u1-u0τ=Δϵ2-Δα/2u1/2+f1/2-12Δu1-u02u1/2,where u1/2=u0+u1/2 and f1/2=f0+f1/2.For this semidiscretized problem, we have the following energy stability estimate, which is a slight modification of the Theorem 4.1 for the standard CH equation in [13].Theorem 1.
Form≥1, define(8)Eαum+1,um=Eαum+1+14τum+1-um-12+12um+1-um02,where Eα(u)=ϵ2/2(-Δ)α/4u2+1/4uL44-1/2u2 and suppose A≥1/16. Then the semidiscrete scheme (6) is unconditionally energy stable, i.e., for all τ>0(9)Eαum+1,um≤Eαum,um-1.Here um-12≔(um,(-Δ)-1um) is inner product of H-1(Ω).Proof.
We use·2 to denote the L2 norm and (·,·) to the inner product in L2. Taking the L2 inner product for (6) with (-Δ)-1(um+1-um) yields(10)
3
u
m
+
1
-
4
u
m
+
u
m
-
1
2
τ
,
-
Δ
-
1
u
m
+
1
-
u
m
+
ϵ
2
-
Δ
1
+
α
/
2
u
m
+
1
,
-
Δ
-
1
u
m
+
1
-
u
m
+
-
Δ
u
m
+
1
3
,
-
Δ
-
1
u
m
+
1
-
u
m
+
Δ
2
u
m
-
u
m
-
1
,
-
Δ
-
1
u
m
+
1
-
u
m
+
A
τ
-
Δ
2
u
m
+
1
-
u
m
,
-
Δ
-
1
u
m
+
1
-
u
m
=
J
1
+
J
2
+
J
3
+
J
4
+
J
5
=
0
.Now we estimate the terms J1-J5, respectively. For the time difference term, by the definition of the norm um-12, it is easy to obtain (11)J1=12τ3um+1-4um+um-1,-Δ-1um+1-um=14τ6um+1-8um+2um-1,-Δ-1um+1-um=τ54um+1-umτ-12-14um-um-1τ-12+τ34um+1-2um+um-1τ2-12.
For the highest-order diffusion, using the identity(12)b,b-a=12b2-12a2+12b-a2,we can obtain (13)J2=ϵ2-Δ1+α/2um+1,-Δ-1um+1-um=ϵ22-Δα/4um+12--Δα/4um2+ϵ22-Δα/4um+1-um2.By virtue of the identity (14)b3,b-a=14b4-a4+14b2-a22+12bb-a2,the nonlinear term becomes (15)J3=um+13,um+1-um=14um+1L44-umL44+14um+12-um22+12um+1um+1-um2.By means of the identity(16)b,b-a=12b2-12a2+12b-a2,the concave diffusive term turns out to be (17)J4=-2um-um-1,um+1-um=-um,um+1-um-um-um-1,um+1-um=-12um+12-um2+12um+1-um2+um+1-2um+um-1,um+1-um-um+1-um2=-12um+12+12um+1-um2--12um2+12um-um-12+12um+1-2um+um-12-12um+1-um2.Finally, the stabilizing term becomes (18)J5=Aτ∇um+1-um2.
For anyσ>0 and the formula from the last row in J4, in virtue of the Cauchy-Schwarz inequality, we have(19)12um+1-um2≤12∇um+1-∇umum+1-um≤τ4σ∇um+1-∇um2+στ4um+1-umτ-12
Putting everything together, we can get(20)Eαum+1,um+τ1-σ4um+1-umτ-12+τA-14σ∇um+1-um2≤Eum,um-1.
Consequently, dropping some unnecessary terms, we can obtain(21)Eαum+1,um≤Eαum,um-1.This completes the proof.Remark. The function (22)Eαum+1,um=Eαum+1+14τum+1-um-12+12um+1-um02consists of three parts. The formula Eα(um+1) comes from a typical free energy functional Eα(u). The formula 1/4τum+1-um-12 and 1/2um+1-um02 comes from the time difference term J1 and the concave diffusive term J4, respectively.For the initialization scheme (7), we have the following stability result.Theorem 2 (see [11]).
The initialization scheme (7) has the energy-decay property(23)Eαu1≤Eαu0.Moreover, we have the following mass conservation law.Theorem 3.
For (1), the semidiscrete scheme (6) satisfies the mass conservation law: (24)Mum+1=Mum,m=0,1,2,…,M-1Proof.
For m=0, it is easy to haveM(u1)=M(u0) based on Theorem 3.2. in [11].
Form>0, rewrite the scheme (6) as (25)3um+1-4um+um-12τ=Δω,ω=ϵ2-Δα/2um+1+um+13-2um+um-1-AτΔum+1-um.
Then we have(26)3um+1-4um+um-12τ,1=Δω,1=-∇ω,∇1=0.Hence(27)3um+1-um,1=um-um-1,1.Then combined with the result (u0,1)=(u1,1) (please see Theorem 3.2. in [11]), it is easy to verify that (28)Mum+1=Mum.This completes the proof.


## 2.2. Spatial Discretization Using Spectral Methods

In this subsection, we will give the fully discrete scheme for model (1). We know that the fractional Laplacian (-Δ)α/2 can be defined in various ways. A review of these existing definitions can be seen in [14]. One important and effective definition for (-Δ)α/2 is defined using Fourier transform; see [15–17].Without loss of generality, we setΩ=[a,b]d and multi-index(29)μ=p1,p2,…,pd∣pi∈Zμ+=p=p1,p2,…,pd∣pi∈N.Then, for p∈μ, suppose the d-dimensional Laplacian (-Δ) has a complete set of orthonormal eigenfunctions φp corresponding to eigenvalues λp on the bounded region [a,b]d, i.e., (-Δ)φp=λpφp. Let (30)Uα=u=∑l=1d∑pl=-∞∞u^pφp,u^p=u,φp,∑l=1d∑pl=-∞∞u^pλpα/2<∞,.Then for any u∈Uα, it is easy to have(31)-Δα/2u=∑l=1d∑pl=-∞∞λpα/2u^pφp,where λp and φp will depend on the specified boundary conditions.For periodic boundary condition,(32)λp=∑l=1d2plπb-a2,φp=1b-ad/2∏l=1dexp⁡i2plπx-ab-a,p∈μ.Here, the imaginary unit is denoted by i.For homogeneous Neumann boundary condition,(33)λp=∑l=1nplπb-a2,φp=2b-ad/2∏l=1dcos⁡plπx-ab-a,p∈μ+.In this work, we use (31) to handle the fractional Laplacian (-Δ)α/2.After discretizing governing equation in time as described in the above subsection, now we will further give the fully discrete scheme based on Fourier spectral method.By applying the Fourier transform to both sides of (6), together with f(u)=u3-u, we have the following nonlinear scheme for the discrete p-th Fourier coefficients (34)3upm+1^-4upm^+upm-1^2τ=-ϵ2λp1+α/2upm+1^-λpupm+13^+2λpqupm^-λpqupm-1^-Aτλp2upm+1^-upm^.For the nonlinear case, we need use an iterative method to solve the resulting algebraic system. Thus the above nonlinear term is treated using the following fixed point iteration: forl=1,2,…,L find upm+1,l^ such that(35)upm+1,l^=13+2τϵ2λp1+α/2+2Aτ2λp24+4τλp+2Aτ2λp2upm^-1+2τλpupm-1^-2τλpupm+1,l-13^,with upm+1,0^=upm^.It also preserves the mass conservation and energy dissipation properties in the full discrete sense. The proofs are similar. We thus omit them for simplicity.


## 3. Numerical Experiments

In this section, we investigate the performance of the obtained scheme (35). Four numerical experiments in 1D and 2D are carried out, some of the results are compared with the theoretical analysis from previous section, and all of them are consistent with the theoretical results.


### 3.1. Convergence Test

In order to verify the convergence rate of the proposed method, a 1D fractional CH problem (1) with homogeneous Neumann boundary condition is first considered. We take Ω=[-π,π], T=1, ϵ2=0.1, and the initial condition u0x=0.1cos⁡(x).To test the accuracy of the numerical solutions, we choose a small grid sizeh=2π/100 and calculate the L2-norm error using Errττ/2=Uτ/2-Uτ2. We further estimate the rate of convergence using the formula as Rate=log⁡2(Errττ/2/Errτ/2τ/4).The computational results corresponding toα=1.8,1.4 are shown in Tables 1-2 for different time steps and numbers of fixed-point iterations L. One may see that the convergence rate in time is close to 2 for all cases, which was expected from the theoretical analysis in the previous section. Moreover, it is clear that the iteration number L=2 in (35) is sufficient to get accurate results.Table 1
The discreteL2-norm errors and corresponding temporal convergence rates with α=1.8, ϵ2=0.1, h=π/50 and different numbers of fixed-point iterations L for Section 3.1.
L τ 1 20 1 40 1 80 1 160 1 320 2 E r r τ τ / 2 2.1916e-3 8.2837e-4 2.3878e-4 6.2590e-5 1.5914e-5 Rate - 1.40 1.79 1.93 1.98  4 E r r τ τ / 2 2.1959e-3 8.2959e-4 2.3905e-4 6.2649e-5 1.5927e-5 Rate - 1.40 1.80 1.93 1.98  10 E r r τ τ / 2 2.1959e-3 8.2959e-4 2.3905e-4 6.2649e-5 1.5927e-5 Rate - 1.40 1.80 1.93 1.80Table 2
The same as in Table1, but with α=1.4.
L τ 1 20 1 40 1 80 1 160 1 320 2 E r r τ τ / 2 1.1681e-2 5.6884e-3 1.8195e-3 4.9437e-4 1.2724e-4 Rate - 1.04 1.64 1.88 1.96  4 E r r τ τ / 2 1.1687e-2 5.6882e-3 1.8190e-3 4.9414e-4 1.2717e-4 Rate - 1.04 1.64 1.88 1.96  10 E r r τ τ / 2 1.1687e-2 5.6882e-3 1.8190e-3 4.9414e-4 1.2717e-4 Rate - 1.04 1.64 1.88 1.96


### 3.2. The Relationship amongα, the Width of the Transition Layer, and the Equilibrium Solution

In this numerical experiment, we investigate the relationship amongα, the width of the transition layer, and the equilibrium solution. We still use the same initial condition as in Section 3.1, and the simulation parameters are chosen as follows:(36)ϵ2=0.04,T=20,h=π50,τ=1100,L=2.Figure1 shows the profiles of the phase-field at t=20 and time evolution of solution for 1D fractional CH equation with α=2, 1.5, 1.2. We observe that small α can generate sharper interface than large α does; a similar phenomenon has also been found in [18] for the space fractional Allen-Cahn equation. Furthermore, due to the long-range dependence and the heavy-tailed influence of the fractional diffusion process, the lifetime of the unstable interface is largely prolonged as α decreases.Figure 1
Comparison of the interface widths and numerical solutions at T=20 with differentα for Section 3.2.
(a)
u(x,T) with differentα (b)
u(x,t) withα = 2 (c)
u(x,t) withα = 1.5 (d)
u(x,t) withα = 1.2


### 3.3. Effects on Phase Behavior ofα and Boundary Conditions

In order to investigate effects on phase behavior of periodic and homogeneous Neumann boundary conditions, the third problem is considered on(-1,1)2×(0,T] with both kinds of the boundary conditions and the following initial condition: (37)u0x,y=1,ifx-0.52+y-0.52<164,1,ifx+0.52+y+0.52<164,0,else.The simulation parameters are chosen as follows:(38)ϵ2=0.01,T=20,h=140,τ=1100,L=2.For the 2D fractional CH equation, we study the coalescence dynamics of two separated bubbles by varying the fractional orderα and boundary conditions. The snapshots of the solutions for both kinds of the boundary conditions and various values of α are shown in Figure 2. We observe that the solutions for different boundary conditions have significant differences and converge to different steady states. As α decreases, the interface thickness becomes thinner and the lifetime of the unstable interface is largely prolonged. For example, for α=2, the solution with periodic boundary condition has arrived the steady state at t=2, whereas the solutions for α=1.8 are still unstable, which further support the results of Section 3.2.Figure 2
Section3.3 for α=2 and α=1.8 with both boundary conditions: the first two rows for α=2 with periodic and Neumann boundary conditions, respectively; the next two rows for α=1.8 with periodic and Neumann boundary conditions, respectively.We also study three important physical quantities of the space fractional CH equation (1): the energy Eα, the mass M, and the surface roughnessSR using the formula as(39)SRt=1Ω∫Ωux,y,t-u¯t2dΩ,whereu¯t=1Ω∫Ωux,y,tdΩ.The corresponding quantities for both kinds of the boundary conditions and various values ofα are shown in Figure 3, which represents the nonincreasing energy and conservation of mass of numerical solutions. Moreover, combined with the snapshots in Figure 2, we observe that for the space fractional CH equation (1) with periodic boundary condition, the time to reach steady state is much shorter than problem (1) with periodic boundary condition. For example, the energy curves in Figure 3 indicate that the solutions with periodic boundary condition arrive the steady state at around t=1.8 for α=2 (see Figure 3(d)) and t=2.5 for α=1.8 (see Figure 3(j)), respectively whereas the solution with Neumann boundary condition is stable at around t=6 for α=2 (see Figure 3(a)) and t=17 for α=1.8 (see Figure 3(g)). This phenomenon is in good agreement with that reported in [11]. Meanwhile, the above results show again that the lifetime of the unstable interface is largely prolonged as the fractional order α decreases.Figure 3
Section3.3 for the curves of energy, mass, and surface roughness computed with h=1/40 and τ=1/100: the first two rows, respectively, for Neumann and periodic boundary conditions with α=2; the next two rows, respectively, for Neumann and periodic boundary conditions with α=1.8.
(a)
Energy withα = 2 (b)
Mass withα = 2 (c)
Surface roughness withα = 2 (d)
Energy withα = 2 (e)
Mass withα = 2 (f)
Surface roughness withα = 2 (g)
Energy withα = 1.8 (h)
Mass withα = 1.8 (i)
Surface roughness withα = 1.8 (j)
Energy withα = 1.8 (k)
Mass withα = 1.8 (l)
Surface roughness withα = 1.8


### 3.4. Comparison of Different Methods

Moreover, for comparison purpose we review an alternative numerical method [11], which is an unconditionally energy stable Crank-Nicolson (CN) scheme for space fractional CH equation (1):(40)upqm+1^-upqm^τ=12-λpqupqm3^+λpqupqm^-ϵλpq1+α/2upqm^-λpqupqm+13^+λpqupqm+1^-ϵ2λpq1+α/2upqm+1^+14λpqum+1-um2um+1+um^.For this scheme, the spatial and temporal discretization errors are also of spectral and second-order accuracy, respectively.Since the time to reach steady state for model (1) with periodic boundary condition is much shorter than the Neumann case (as shown in Section 3.3), we only consider the periodic boundary condition for convenient comparison. The fourth problem from [11, 19] with ϵ2=0.01 is considered on (0,2π)2×(0,20] with the following initial condition (41)u0x,y=0.05sin⁡xsin⁡y+0.001.In this test, we fix the spatial steph=2π/80. In order to compare the performance of our method (denoted by Mnew) with the method (40) (denoted by Mold), the numerical results and the computational times (in seconds) for different α are presented in Table 3. It is observed that the scheme Mold is not convergence when τ≥0.0013 for α=2 and τ≥0.005 for α=1.6, whereas our scheme Mnew can work well even with τ=0.02 for α=2 and τ=0.0075 for α=1.6. This further suggests that the obtained scheme in this paper can lead to a substantial reduction of the CPU time, while keeping high accuracy of the computed solution.Table 3
A comparison of the CPU time (in seconds) used by our scheme (35) (denoted by Mnew) and the method in [11] (denoted by Mold).
α=2 α=1.6 M o l d M n e w M o l d Mnew  τ ≥0.005 0.004 >0.02 0.02 0.01 ≥ 0.0013 0.001 > 0.0075 0.0075 0.004  Convergence No Yes No Yes Yes No Yes No Yes Yes  CPU - 89.8 - 10.5 20.0 - 332.4 - 28.6 52.0The numerical solutions att=20 are presented in Figures 4(a)–4(d), it is clear that the interface thickness becomes thinner as the fractional power is decreased. Besides that fact that both methods (35) and (40) reach the same equilibrium solution, only the scheme (40) arrives to an inverse equilibrium solution for α=2 (of course, the surface area of all numerical solutions reached the minimum value.). As a conclusion we claim that from the computational point of view, the scheme (35) is better than (40) because we can take bigger time steps (reducing the computational cost) and still obtain a good behavior of the scheme.Figure 4
Comparison of the numerical solutions at T=20 andEα, SR with different α and τ for Section 3.4.
(a)
M o l d: τ = 0.004,α = 2 (b)
M o l d: τ = 0.001,α = 1.6 (c)
M n e w: τ = 0.02,α = 2 (d)
M n e w: τ = 0.0075,α = 1.6 (e)
α = 2 (f)
α = 1.6Figures4(e) and 4(f) show the evolution of the free energy Eα and surface roughness SR in time for both schemes. The energies and roughnesses do not evolve in the same way but all of them reach the same equilibrium states. Furthermore, with the decrease of the time steps, the dynamic of Eα and SR for scheme (35) is almost the same for (40).Remark. From Table 1 we find that the critical time step for α=1.6 is much smaller than that of α=2. We would like to comment with the fact that since the thickness of the interface is inversely proportional to α (see Figures 4(a)–4(d)), the case for small α requires a finer time step. In addition, since phase-field models tend to sharp-interface models as the thickness of the interfaces tends to zero; this case is more interesting from the physical point of view.


## 3.1. Convergence Test

In order to verify the convergence rate of the proposed method, a 1D fractional CH problem (1) with homogeneous Neumann boundary condition is first considered. We take Ω=[-π,π], T=1, ϵ2=0.1, and the initial condition u0x=0.1cos⁡(x).To test the accuracy of the numerical solutions, we choose a small grid sizeh=2π/100 and calculate the L2-norm error using Errττ/2=Uτ/2-Uτ2. We further estimate the rate of convergence using the formula as Rate=log⁡2(Errττ/2/Errτ/2τ/4).The computational results corresponding toα=1.8,1.4 are shown in Tables 1-2 for different time steps and numbers of fixed-point iterations L. One may see that the convergence rate in time is close to 2 for all cases, which was expected from the theoretical analysis in the previous section. Moreover, it is clear that the iteration number L=2 in (35) is sufficient to get accurate results.Table 1
The discreteL2-norm errors and corresponding temporal convergence rates with α=1.8, ϵ2=0.1, h=π/50 and different numbers of fixed-point iterations L for Section 3.1.
L τ 1 20 1 40 1 80 1 160 1 320 2 E r r τ τ / 2 2.1916e-3 8.2837e-4 2.3878e-4 6.2590e-5 1.5914e-5 Rate - 1.40 1.79 1.93 1.98  4 E r r τ τ / 2 2.1959e-3 8.2959e-4 2.3905e-4 6.2649e-5 1.5927e-5 Rate - 1.40 1.80 1.93 1.98  10 E r r τ τ / 2 2.1959e-3 8.2959e-4 2.3905e-4 6.2649e-5 1.5927e-5 Rate - 1.40 1.80 1.93 1.80Table 2
The same as in Table1, but with α=1.4.
L τ 1 20 1 40 1 80 1 160 1 320 2 E r r τ τ / 2 1.1681e-2 5.6884e-3 1.8195e-3 4.9437e-4 1.2724e-4 Rate - 1.04 1.64 1.88 1.96  4 E r r τ τ / 2 1.1687e-2 5.6882e-3 1.8190e-3 4.9414e-4 1.2717e-4 Rate - 1.04 1.64 1.88 1.96  10 E r r τ τ / 2 1.1687e-2 5.6882e-3 1.8190e-3 4.9414e-4 1.2717e-4 Rate - 1.04 1.64 1.88 1.96


## 3.2. The Relationship amongα, the Width of the Transition Layer, and the Equilibrium Solution

In this numerical experiment, we investigate the relationship amongα, the width of the transition layer, and the equilibrium solution. We still use the same initial condition as in Section 3.1, and the simulation parameters are chosen as follows:(36)ϵ2=0.04,T=20,h=π50,τ=1100,L=2.Figure1 shows the profiles of the phase-field at t=20 and time evolution of solution for 1D fractional CH equation with α=2, 1.5, 1.2. We observe that small α can generate sharper interface than large α does; a similar phenomenon has also been found in [18] for the space fractional Allen-Cahn equation. Furthermore, due to the long-range dependence and the heavy-tailed influence of the fractional diffusion process, the lifetime of the unstable interface is largely prolonged as α decreases.Figure 1
Comparison of the interface widths and numerical solutions at T=20 with differentα for Section 3.2.
(a)
u(x,T) with differentα (b)
u(x,t) withα = 2 (c)
u(x,t) withα = 1.5 (d)
u(x,t) withα = 1.2


## 3.3. Effects on Phase Behavior ofα and Boundary Conditions

In order to investigate effects on phase behavior of periodic and homogeneous Neumann boundary conditions, the third problem is considered on(-1,1)2×(0,T] with both kinds of the boundary conditions and the following initial condition: (37)u0x,y=1,ifx-0.52+y-0.52<164,1,ifx+0.52+y+0.52<164,0,else.The simulation parameters are chosen as follows:(38)ϵ2=0.01,T=20,h=140,τ=1100,L=2.For the 2D fractional CH equation, we study the coalescence dynamics of two separated bubbles by varying the fractional orderα and boundary conditions. The snapshots of the solutions for both kinds of the boundary conditions and various values of α are shown in Figure 2. We observe that the solutions for different boundary conditions have significant differences and converge to different steady states. As α decreases, the interface thickness becomes thinner and the lifetime of the unstable interface is largely prolonged. For example, for α=2, the solution with periodic boundary condition has arrived the steady state at t=2, whereas the solutions for α=1.8 are still unstable, which further support the results of Section 3.2.Figure 2
Section3.3 for α=2 and α=1.8 with both boundary conditions: the first two rows for α=2 with periodic and Neumann boundary conditions, respectively; the next two rows for α=1.8 with periodic and Neumann boundary conditions, respectively.We also study three important physical quantities of the space fractional CH equation (1): the energy Eα, the mass M, and the surface roughnessSR using the formula as(39)SRt=1Ω∫Ωux,y,t-u¯t2dΩ,whereu¯t=1Ω∫Ωux,y,tdΩ.The corresponding quantities for both kinds of the boundary conditions and various values ofα are shown in Figure 3, which represents the nonincreasing energy and conservation of mass of numerical solutions. Moreover, combined with the snapshots in Figure 2, we observe that for the space fractional CH equation (1) with periodic boundary condition, the time to reach steady state is much shorter than problem (1) with periodic boundary condition. For example, the energy curves in Figure 3 indicate that the solutions with periodic boundary condition arrive the steady state at around t=1.8 for α=2 (see Figure 3(d)) and t=2.5 for α=1.8 (see Figure 3(j)), respectively whereas the solution with Neumann boundary condition is stable at around t=6 for α=2 (see Figure 3(a)) and t=17 for α=1.8 (see Figure 3(g)). This phenomenon is in good agreement with that reported in [11]. Meanwhile, the above results show again that the lifetime of the unstable interface is largely prolonged as the fractional order α decreases.Figure 3
Section3.3 for the curves of energy, mass, and surface roughness computed with h=1/40 and τ=1/100: the first two rows, respectively, for Neumann and periodic boundary conditions with α=2; the next two rows, respectively, for Neumann and periodic boundary conditions with α=1.8.
(a)
Energy withα = 2 (b)
Mass withα = 2 (c)
Surface roughness withα = 2 (d)
Energy withα = 2 (e)
Mass withα = 2 (f)
Surface roughness withα = 2 (g)
Energy withα = 1.8 (h)
Mass withα = 1.8 (i)
Surface roughness withα = 1.8 (j)
Energy withα = 1.8 (k)
Mass withα = 1.8 (l)
Surface roughness withα = 1.8


## 3.4. Comparison of Different Methods

Moreover, for comparison purpose we review an alternative numerical method [11], which is an unconditionally energy stable Crank-Nicolson (CN) scheme for space fractional CH equation (1):(40)upqm+1^-upqm^τ=12-λpqupqm3^+λpqupqm^-ϵλpq1+α/2upqm^-λpqupqm+13^+λpqupqm+1^-ϵ2λpq1+α/2upqm+1^+14λpqum+1-um2um+1+um^.For this scheme, the spatial and temporal discretization errors are also of spectral and second-order accuracy, respectively.Since the time to reach steady state for model (1) with periodic boundary condition is much shorter than the Neumann case (as shown in Section 3.3), we only consider the periodic boundary condition for convenient comparison. The fourth problem from [11, 19] with ϵ2=0.01 is considered on (0,2π)2×(0,20] with the following initial condition (41)u0x,y=0.05sin⁡xsin⁡y+0.001.In this test, we fix the spatial steph=2π/80. In order to compare the performance of our method (denoted by Mnew) with the method (40) (denoted by Mold), the numerical results and the computational times (in seconds) for different α are presented in Table 3. It is observed that the scheme Mold is not convergence when τ≥0.0013 for α=2 and τ≥0.005 for α=1.6, whereas our scheme Mnew can work well even with τ=0.02 for α=2 and τ=0.0075 for α=1.6. This further suggests that the obtained scheme in this paper can lead to a substantial reduction of the CPU time, while keeping high accuracy of the computed solution.Table 3
A comparison of the CPU time (in seconds) used by our scheme (35) (denoted by Mnew) and the method in [11] (denoted by Mold).
α=2 α=1.6 M o l d M n e w M o l d Mnew  τ ≥0.005 0.004 >0.02 0.02 0.01 ≥ 0.0013 0.001 > 0.0075 0.0075 0.004  Convergence No Yes No Yes Yes No Yes No Yes Yes  CPU - 89.8 - 10.5 20.0 - 332.4 - 28.6 52.0The numerical solutions att=20 are presented in Figures 4(a)–4(d), it is clear that the interface thickness becomes thinner as the fractional power is decreased. Besides that fact that both methods (35) and (40) reach the same equilibrium solution, only the scheme (40) arrives to an inverse equilibrium solution for α=2 (of course, the surface area of all numerical solutions reached the minimum value.). As a conclusion we claim that from the computational point of view, the scheme (35) is better than (40) because we can take bigger time steps (reducing the computational cost) and still obtain a good behavior of the scheme.Figure 4
Comparison of the numerical solutions at T=20 andEα, SR with different α and τ for Section 3.4.
(a)
M o l d: τ = 0.004,α = 2 (b)
M o l d: τ = 0.001,α = 1.6 (c)
M n e w: τ = 0.02,α = 2 (d)
M n e w: τ = 0.0075,α = 1.6 (e)
α = 2 (f)
α = 1.6Figures4(e) and 4(f) show the evolution of the free energy Eα and surface roughness SR in time for both schemes. The energies and roughnesses do not evolve in the same way but all of them reach the same equilibrium states. Furthermore, with the decrease of the time steps, the dynamic of Eα and SR for scheme (35) is almost the same for (40).Remark. From Table 1 we find that the critical time step for α=1.6 is much smaller than that of α=2. We would like to comment with the fact that since the thickness of the interface is inversely proportional to α (see Figures 4(a)–4(d)), the case for small α requires a finer time step. In addition, since phase-field models tend to sharp-interface models as the thickness of the interfaces tends to zero; this case is more interesting from the physical point of view.


## 4. Conclusions

In this work, we have presented a second order BDF-type scheme in time to solve fractional-in-space CH model. The energy stability and mass conservation of the semi-discrete scheme have been proved. Moreover, the space has been discretized by Fourier spectral method. Several numerical examples have verified our theoretical results. Via numerical comparison, we have shown that, by introducing a second-order Douglas-Dupont-type regularization based on the implicit backward difference scheme, larger time step can also guarantee stability and accuracy. We observe from the numerical results the lifetime of the unstable interface is largely prolonged as the fractional orderα decreases and small α can generate sharper interface.



---

*Source: 3163702-2019-04-01.xml*