A decoupled finite element method for a modified Cahn-Hilliard-Hele-Shaw system

In this paper, a decoupled finite element method for a modified Cahn-Hilliard-Hele-Shaw system with double well potential is constructed. The proposed scheme is based on the convex splitting of the energy functional in time and uses the mixed finite element discretzation in space. The computation of the velocity u is separated from the computation of the pressure p by using an operatorsplitting strategy. A Possion equation is solved to update the pressure at each time step. Unconditional stability and error estimates are analyzed in detail theoretically. Furthermore, the theoretical part is verified by several numerical examples, whose results show that the numerical examples are consistent with the results of the theoretical part.


Introduction
The Cahn-Hilliard-Hele-Shaw system is a very important mathematical model which describes the motion of a viscous incompressible fluid between two closely spaced parallel plates and can be viewed as the simplification of the Cahn-Hilliard-Navier-Stokes system [1][2][3]. The model are widely applied in different fields, such as simulations of nonlinear tumor growth and neovascularization [4][5][6][7], spinodal decomposition in a Hele-Shaw cell [8], and two-phase flow in porous medium [9,10], etc.
The Cahn-Hilliard-Hele-Shaw system is a gradient system coupled with fluid motion, which is difficult to solve because of its complex form. For this model, purely explicit methods are limited by strict time step constraints for stability, and completely implicit numerical methods must contend with potentially large systems of nonlinear algebraic equations [11]. There have been many effective numerical schemes for the Cahn-Hilliard-Hele-Shaw system. Guo et al. proposed a semi-implicit time integration scheme based on convex splitting technique, and proved the unconditional stability of the fully discrete scheme of the Cahn-Hilliard-Hele-Shaw system [12]. S.M. Wise put forward an unconditionally stable finite difference scheme for the Cahn-Hilliard-Hele-Shaw [13]. Chen et al. established a finite difference simulation of Gagliardo-Nirenberg-type inequalities to analyze stability and convergence [14]. Liu et al. developed a mixed finite element numerical scheme for the Cahn-Hilliard-Hele-Shaw system and proved its unconditional stability [15]. Guo carried out a numerical analysis for the Cahn-Hilliard-Hele-Shaw system with variable mobility and logarithmic Flory-Huggins potential [16]. The above mentioned works are numerical methods to solve the Cahn-Hilliard-Hele-Shaw system.
However, there are few researches on the modified Cahn-Hilliard-Hele-Shaw system.
The modified Cahn-Hilliard equation (also named Cahn-Hilliard-Oono equation) used to suppress phase coarsening in [17] is as follows where φ 0 := 1 |Ω| Ω φ 0 (x)dx. More works on the modified Cahn-Hilliard equation can be found in [18][19][20][21]. For the modified Cahn-Hilliard equation, when θ = 0, the equation becomes the classical Cahn-Hilliard equation [22][23][24]. When the modified Cahn-Hilliard is coupled with the Darcy equation, the modified Cahn-Hilliard-Hele-Shaw equation can be obtained. Jia et al. introduced a novel finite element method for the modified Cahn-Hilliard-Hele-Shaw system [25], in which the time discretization was based on the convex splitting of the energy functional in the modified Cahn-Hilliard equation. Of course, the above numerical methods are directly solved based on the coupling equation, and the solving process is complicated. To solve this kind of problem, many decoupled methods have been proposed to solve the Cahn-Hilliard-Hele-Shaw system in recent years. Han [26] presented a decoupled unconditionally stable numerical scheme for the Cahn-Hilliard-Hele-Shaw system with variable viscosity, in which the operator-splitting strategy and the pressure-stabilization technique were used to completely decouple the nonlinear Cahn-Hilliard equation from pressure. Similar strategies were also adopted in [27]. Then, Gao [28] studied the fully decoupled numerical scheme of the Cahn-Hilliard-Hele-Shaw model, in which the scalar auxiliary variable method was used to deal with the nonlinear term in the free energy. Similarly, decoupled schemes are also effectively used in other systems and models recently. Zhao et al. [29] developed an energy-stable scheme for a binary hydrodynamic phase field model of mixtures of nematic liquid crystals and viscous fluids. A second-order decoupled energy-stable schemes for Cahn-Hilliard-Navier-Stokes equations was suggested in [30]. For thermodynamically consistent models, Zhao [31] investigated a general numerical framework for designing linear, energy stable, and decoupled numerical algorithms. However, to the best of our knowledge, there are few researches on decoupling methods of the modified Cahn-Hilliard-Hele-Shaw system, it will be the purpose of our paper.
Based on Eqs (1.1)-(1.3), the modified Cahn-Hilliard-Hele-Shaw system with double well potential is given by in Ω × (0, T ), (1.10) where Ω ∈ R d (d = 2, 3). φ is the concentration field, u is the advective velocity, ε > 0 is the constant to measure the thickness of the transition layer between the two phases, and µ is the chemical potential. f (φ) is the derivative of the double well potential F(φ), ξ is an auxiliary variable. p and γ represent the pressure and the dimensionless surface tension parameter, respectively. n is the unit outer normal of the boundary ∂Ω. when θ = 0, the equation becomes the classical Cahn-Hilliard-Hele-Shaws equation.
With regard to the double well potential corresponding to f (φ) in Eq (1.2), the followingF(φ) can be taken [32][33][34] (1.11) Correspondingly, the derivatives ofF(φ) can be split as followš (1.12) F(φ) and f (φ) are replaced byF(φ) and its derivativef (φ), which are still recorded as F and f for simplicity. Typically, the free energy functional of a modified Cahn-Hilliard-Hele-Shaw system with double well potential is given by In this paper, a decoupled finite element scheme for a modified Cahn-Hilliard-Hele-Shaw system with double well potential is proposed. The temporal discretization is based on the convex splitting of the energy functional in the modified Cahn-Hilliard equation, and the spacial discretization is carried out by the mixed finite element method. The computation of the velocity u is separated from the computation of the pressure p by using an operator-splitting strategy. A Possion equation is solved to update the pressure at each time step. We prove that the the proposed scheme is unconditionally stable in energy, and the error analyses are obtained. Finally, the numerical results verify the theoretical analysis. The rest of this article is structured as follows.The finite element discrete scheme of the Cahn-Hilliard-Hele-Shaw system combing with the convex splitting is given in Section 2; The theoretical preparations and stability of the proposed numerical scheme are proved in Section 3; The error analyses of the proposed scheme are addressed in Section 4; Some numerical examples are given to verify the previous theory in Section 5, and the conclusion is given in Section 6.

Semi-discrete scheme
Let L 2 (Ω) is a space of square integrable function and H k (Ω), H k 0 (Ω) denote the usual Sobolev spaces. L 2 (Ω) inner product and its norm are denoted by (u . The weak formulation of the modified Cahn-Hilliard-Hele-Shaw system with double well potential can be written as where, Let N be a positive integer and 0 = t 0 < t 1 < · · · < t N = T be a uniform partition of [0, T ], where The semi-discrete scheme of the modified Cahn-Hilliard-Hele-Shaw system with double well potential is as follows. For n ≥ 0, find {φ n+1 , µ n+1 , ξ n+1 , p n+1 } such that where the velocity is given by Combing with the idea of the literatures [26,35], the computation of the modified Cahn-Hilliard equations (2.2)-(2.4) are decoupled from Eq (2.5) after substituting u n+1 into Eq (2.2), since the pressure is explicit in Eq (2.6). The velocity u n+1 is regarded as an intermediate velocity by using the incremental projection method similar to the Navier-Stokes equation. The real velocityũ n+1 is obtained from the intermediate velocity and satisfies Then Eq (2.6) and Eq (2.7) are added together to obtain the original Eq (1.7). If the divergence operator is applied to both side of Eq (2.7), the real velocityũ n+1 will vanished. We have (2.8)

Fully discrete scheme
Let T h = {K} be a regular partition of the domain Ω that is divided into triangles with the size h = max 0≤i≤N h i . S h is a piecewise polynomial space, which is defined as where P k (x, y) is a polynomial of degree at most r.
Let us denote . The corresponding fully discrete scheme have the following expression, find where the velocity is given by 3. The analysis of stability and have the following estimates, Define the operator T h :Ĥ −1 →Ĥ 1 through the following variational problems, Lemma 3.1. [12,15] Let ζ, ϕ ∈Ĥ −1 and set where (·, ·) −1,h defines an inner product on theĤ −1 and its corresponding H −1 norm is written as Consequently, for ∀χ ∈Ĥ 1 , ζ ∈Ĥ −1 , Furthermore, the following Poincaré inequalities holds, Lemma 3.2. [12,15] Projection operator P is linear and satisfies the following properties

11)
and where p h ∈Ŝ h is the unique solution to (3.14) Lemma 3.3. [14,15] Projection operator P h is linear and satisfies the following properties

15)
and Lemma 3.4. [12,15] Suppose that w ∈ H q (Ω) with the compatible boundary condition w · n = 0 on ∂Ω and q ∈ H q+1 (Ω), then Then for any h, τ, ε > 0, n ≥ 0, scheme (2.9)-(2.12) satisfies the following property, where η is a number between φ n h and φ n+1 h , we have Then, choosing w h = −(φ n+1 h − φ n h ) and using the fact that (a, Next, we take inner product of Eq (2.13) with Now, taking q h = τ γ p n h and using the fact that (a The proof is completed.
there is a constant C > 0 independent of τ and h, such that the following estimates hold for any τ, h > 0, Proof. Summing the Eq (3.19) from i = 0 to N, we get The proof is completed.

Numerical analysis
In this part, some numerical examples are used to verify the correctness and validity of the theoretical analysis.

The spatial convergence order
For Tables 1 and 2, the parameters are chosen as follows, τ = 0.01, T = 0.1, ε=0.14 and mesh steps h = 1 16 , 1 32 , 1 64 , 1 128 . The spatial convergence orders of relative error ê φ H 1 are close to 1, which is consistent with the convergence order obtained from theoretical analysis. Moreover, different θ and γ have little effect on the corresponding convergence order.

The temporal convergence order
For Tables 3 and 4, the parameters are chosen as follows, ε = 0.01, T = 0.1, h = τ =0.0625, 0.03125, 0.015625. The temporal convergence orders of relative error ê φ H 1 are close to 1, which is consistent with the convergence order obtained from theoretical analysis.

Energy dissipation
Let us test the energy dissipation of our proposed scheme. The energy functional Eq (1.13) of the modified Cahn-Hilliard-Hele-Shaw system Eqs (1.4)-(1.10) can be discreteized as Correspondingly, the modified energy of the fully discrete scheme Eqs (2.9)-(2.13) is defined as For the test, the parameters are chosen as follows: T = 5, τ = 0.001, h = 1 64 , ε = 0.4, γ = 0.5. In Figure 1, we can see that the energy functional is non-increasing for θ=0, 0.1, 1.

Spinodal Decomposition
In this part, we present the phase separation dynamics that is called spinodal decomposition in the modified Cahn-Hilliard-Hele-Shaw system. In the simulation, the computational domain is chosen as [0, 1] × [0, 1], the parameters are chosen as follows: ε = 0.05, γ = 0.45, τ = 0.0001. Then, let us take the initial condition where rand() ∈ [0, 1]. The process of coarsening is shown in the following figures. From figures 2-19, we can see that the contours of φ are gradually coarsened over time. However, the profiles obtained by different θ are similar at the same time T . From left to right, the coarsening processes of θ = 50, 200 are not obvious compared with the coarsening processes of θ = 0. We know the bigger θ can suppress the coarsening process.

Summary
In this paper, a decoupled scheme of the modified Cahn-Hilliard-Hele-Shaw system is studied. In our scheme, the velocity and pressure are decoupled, and a Possion equation is solved to update the pressure at each time step. Unconditional stability of the scheme in energy is proved. The convergence analysis are addressed in the frame of finite element method. Furthermore, the theoretical part is verified by several numerical examples. The results show that the numerical examples are consistent with the results of the theoretical part.