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.
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]:where is a bounded regular domain in , and denotes a finite time. For the above equation, we impose the periodic or following homogeneous Neumann boundary conditionson . Here is the unit outward normal vector of the domain . The function 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 and is a given energy potential, taking its global minimum value at and . In other words, the values -1 and 1 of the function prevail in whereas the transition between them forms a thin interface layer.
The fractional CH equation (1) can be viewed as an -gradient flow of the following fractional analogue version of Ginzburg-Landau free energy functional [10]:The same as the standard CH equation [1], one intrinsic property of the fractional-in-space CH equation is energy-decay property in sense ofAnother intrinsic property of the fractional-in-space CH equation, as well as the standard CH equation, is the law of conservation of mass
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 Section 2, 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 , , where . The new second-order scheme can be written aswhere is a stabilizing term with , and is approximated by the following unconditionally stable Crank-Nicolson scheme [11]:where and .
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. For , definewhere and suppose . Then the semidiscrete scheme (6) is unconditionally energy stable, i.e., for all Here is inner product of .
Proof. We use to denote the norm and to the inner product in . Taking the inner product for (6) with yieldsNow we estimate the terms -, respectively. For the time difference term, by the definition of the norm , it is easy to obtain For the highest-order diffusion, using the identitywe can obtain By virtue of the identity the nonlinear term becomes By means of the identitythe concave diffusive term turns out to be Finally, the stabilizing term becomes For any and the formula from the last row in , in virtue of the Cauchy-Schwarz inequality, we havePutting everything together, we can getConsequently, dropping some unnecessary terms, we can obtainThis completes the proof.
Remark. The function consists of three parts. The formula comes from a typical free energy functional The formula and comes from the time difference term and the concave diffusive term , respectively.
For the initialization scheme (7), we have the following stability result.
Theorem 2 (see [11]). The initialization scheme (7) has the energy-decay propertyMoreover, we have the following mass conservation law.
Theorem 3. For (1), the semidiscrete scheme (6) satisfies the mass conservation law:
Proof. For m=0, it is easy to have based on Theorem 3.2. in [11].
For , rewrite the scheme (6) as Then we haveHenceThen combined with the result (please see Theorem 3.2. in [11]), it is easy to verify that 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 can be defined in various ways. A review of these existing definitions can be seen in [14]. One important and effective definition for is defined using Fourier transform; see [15–17].
Without loss of generality, we set and multi-indexThen, for , suppose the -dimensional Laplacian has a complete set of orthonormal eigenfunctions corresponding to eigenvalues on the bounded region , i.e., . Let Then for any , it is easy to havewhere and will depend on the specified boundary conditions.
For periodic boundary condition, Here, the imaginary unit is denoted by .
For homogeneous Neumann boundary condition, In this work, we use (31) to handle the fractional Laplacian .
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 , we have the following nonlinear scheme for the discrete -th Fourier coefficients
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: for find such thatwith .
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 , , , and the initial condition
To test the accuracy of the numerical solutions, we choose a small grid size and calculate the -norm error using . We further estimate the rate of convergence using the formula as
The computational results corresponding to = are shown in Tables 1-2 for different time steps and numbers of fixed-point iterations . 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 in (35) is sufficient to get accurate results.
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:
Figure 1 shows the profiles of the phase-field at 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.
(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 with both kinds of the boundary conditions and the following initial condition:
The simulation parameters are chosen as follows:
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 , the solution with periodic boundary condition has arrived the steady state at , whereas the solutions for are still unstable, which further support the results of Section 3.2.
We also study three important physical quantities of the space fractional CH equation (1): the energy , the mass , and the surface roughness SR using the formula as
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 for (see Figure 3(d)) and for (see Figure 3(j)), respectively whereas the solution with Neumann boundary condition is stable at around for (see Figure 3(a)) and for (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.
(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):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 is considered on with the following initial condition
In this test, we fix the spatial step . In order to compare the performance of our method (denoted by ) with the method (40) (denoted by ), the numerical results and the computational times (in seconds) for different are presented in Table 3. It is observed that the scheme is not convergence when for and for , whereas our scheme can work well even with =0.02 for and =0.0075 for . 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.
The numerical solutions at 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 (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.
(a) : τ = 0.004, = 2
(b) : τ = 0.001, = 1.6
(c) : τ = 0.02, = 2
(d) : τ = 0.0075, = 1.6
(e) = 2
(f) = 1.6
Figures 4(e) and 4(f) show the evolution of the free energy and surface roughness 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 and for scheme (35) is almost the same for (40).
Remark. From Table 1 we find that the critical time step for is much smaller than that of . 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.
Data Availability
The data of numerical experiment in the manuscript used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work of Zhifeng Weng is in part supported by the Natural Science Foundation of China (no. 11701197), and the research of Langyang Huang is in part supported by the Natural Science Foundation of China (no. 11701196) and Progmotion Program for Young and Middle-Aged Teacher in Science and Technology Research of Huaqiao University (ZQN-YX502).