A Decoupled Energy Stable Numerical Scheme for the Modified Cahn–Hilliard–Hele–Shaw System with Logarithmic Potential

A decoupled unconditionally stable numerical scheme for the modied Cahn–Hilliard–Hele–Shaw system with logarithmic potential is proposed in this paper. Based on the convex-splitting of the associated energy functional, the temporal discretization of the scheme is given. e fractional step method is used to decouple the nonlinear modied Cahn–Hilliard equation from the pressure gradient. en, at each time step, one only needs to solve Poisson’s equation which is obtained by using an incremental pressure-stabilization technique. In terms of logarithmic potential, using the regularization procedure can make the domain extended from (− 1, 1) to (− ∞,∞) for F(φ). Finally, the proof of the energy stability and the calculation of the error estimate are presented. Numerical examples are recorded to illustrate the accuracy and performance for the proposed scheme.


Introduction
e Cahn-Hilliard-Hele-Shaw system (CH-HS) [1], which is used to describe the motion of a viscous uid between two closely spaced parallel plates, is an important mathematical model. e applications of the CH-HS model and its variant are abundant. Aritomo Shinozaki and Yoshitsugu Oono performed the simulation of spinodal decomposition of a binary incompressible uid in a Hele-Shaw cell [2] by a similar set of equations. Recently, the CH-HS system has been applied in the study of Sa man-Taylor instability [3] when a more viscous uid is displaced by a less viscous one resulting in complex pattern formation. is system also serves as a phase eld and di usion interface models of tumor growth [4][5][6]. When given thought to permeability/hydraulic conductivity, the CH-HS system, also known as Cahn-Hilliard-Darcy (CH-D) equation, is used to model multiphase ow in porous media. e CH-D system has been studied by many authors and one can refer to Refs. [7,8].
Let Ω ∈ R d (d 2, 3) be an open polygonal or polyhedral domain with a Lipschitz continuous boundary zΩ. e modi ed Cahn-Hilliard-Hele-Shaw system can be studied in the article as follows: with the convertive term, the equations (1d) and (1e) are considered as the Darcy equation with the elastic forcing term, and the equations (1f ) and (1g) are the initial and Neumman boundary conditions, respectively. Numerical methods are a crucial tool to study the dynamics described by the CH-HS system, which have been extensively investigated. For instance, an unconditionally stable and solvable finite difference scheme for the CH-HS equations [11] was presented by Wise et al., which was based on a convex splitting of the discrete Cahn-Hilliard(CH) energy and was semiimplicit. Guo et al. put forward an efficient and stable energy fully discrete local discontinuous Galerkin (LDG) method for the CH-HS system in Ref. [12]. Also, in literature [13], the error estimates of the mixed finite element method of the CH-HS system was analyzed by Chen Wenbin et al. We can refer to references [14][15][16][17].
e modified Cahn-Hilliard equation was introduced by Shahriari [18] to suppress the phase coarsening. A discontinuous Galerkin method coupled with convex splitting for the modified Cahn-Hilliard equation was proposed by Aristotelous et al. [19]. Other related researches for the modified Cahn-Hilliard equation can be seen in Refs. [20][21][22]. When the modified Cahn-Hilliard equation is coupled with the Darcy equation, there is the modified Cahn-Hilliard-Hele-Shaw equation. For the modified Cahn-Hilliard-Hele-Shaw system, Jia presented a coupled finite element method for solving the modified CH-HS system [23], in which the time discretization was based on the convex splitting of the energy functional in the modified CH equation, and the high-order nonlinear term and the linear term of the chemical potential were treated explicitly and implicitly, respectively.
About the potential function, a modified CH-HS system with double well potential was solved in Ref. [23]. We can find that no theoretical and numerical analysis has been presented for the modified CH-HS system with logarithmic potential although the CH-HS equation with logarithmic potential had appeared, such as the nonlocal CH-HS system with logarithmic potential [24]. A decoupling numerical method for solving the CH-HS system with logarithmic potential was proposed by Guo et al. [25], and a fully discrete scheme was proposed for solving the CH-HS system with logarithmic potential and concentration dependent mobility by Guo et al. [25]. erefore, we used the free energy density function of the logarithmic potential in this paper, which is defined as follows [26]: ∀a ∈ (0, 1), e logarithmic free energy functional is defined as and the CH-HS system is mass conservative, i.e., (ϕ(·, t), 1) � (ϕ 0 , 1) and energy dissipative: where (·, ·) is the standard L 2 inner product over Ω.
Here, we use the Elliott and Luckhaus regularization for problem with the logarithmic free energy function F(ϕ) replaced by the twice continuously differentiable function � F(ϕ) [27,28], where ∀κ ∈ (0, 1): (6) and the monotone function is given by Mathematical Problems in Engineering For κ ≤ 1/2 and ∀ϕ, � f(ϕ) ≤ L holds true, where L is a positive constant.
Hereinafter, f and F will be substituted by � F and � f in our analysis. However, for convenience, the �will be omitted.
We note that the equations in (1a)-1g are coupled, nonlinear, and numerically stiff with large spatial derivative over a small transition layer.
us, solving the modified CH-HS system numerically is challenging. In reference [23], Jia presented a coupled finite element method for solving the modified CH-HS system. In fact, compared with the method of coupling, the decoupling is more simple. Daozhi Han introduced a decoupled numerical scheme for the CH-HS system [29], an operator-splitting/fractional-step method is used to decouple the nonlinear Cahn-Hilliard equation from the update of pressure, and a pressure-stabilization technique is used in the update of pressure so that only Poisson's equation with constant coefficient needs to be solved at each time step. Using the method of decoupling makes the calculation easier and faster, which can reduce the cost of calculation.
In this paper, we propose a decoupled unconditionally stable numerical scheme for the modified Cahn-Hilliard-Hele-Shaw system with logarithmic potential among time discretization basing on the convex splitting of the energy functional dated back to eyre [30]. e nonlinear modified CH equation is fully decoupled from the pressure in the Darcy equation by using an operator splitting strategy, and the incremental pressure-stabilization technique [31] is used in pressure updates, so that we only need to solve Possion's equation at each time step. For the logarithmic potential, the use of the regularization procedure makes the domain for the regularized functional F(ϕ) extended from (− 1, 1) to (− ∞, ∞) so that the problem of the overflow on the boundary can be avoided. Finally, the stability and convergence of the scheme are presented.
e Cahn-Hilliard-Hele-Shaw system (CHHS) with logarithmic potential is a strong nonlinear system. Based on this reason, it is extremely difficult to find the exact solution of the system. To overcome this difficulty, we proposed a numerical method. e first-order scheme based on the convex-splitting method is introduced for the CHHS system. e fully discrete scheme exploits the first-order backward Euler method for time discretization and the mixed finite element method for spatial discretization. It is uniquely solvable and stable in energy. e rest of the article is organized as follows. e semidiscrete and fully discrete finite element schemes are reviewed in Section 2. en, we show that the fully discrete scheme is unconditionally stable in energy in Section 3. In Section 4, the detailed convergence analysis is given. Finally, some numerical examples show that the proposed scheme is convergent and efficient in Section 5. Some conclusions are drawn in Section 6.

The Discrete Numerical Scheme
e weak formulation of the modified Cahn-Hilliard-Hele-Shaw equation can be written as and f 2 (ϕ): � ϕ.
Let N be a positive integer and 0 � t 0 < t 1 < · · · < t N � T be a uniform partition of [0, T]. Let us denote by τ: � t n+1 − t n , n � 0, 1 · · · N, where τ is the time step size. For the modified Cahn-Hilliard-Hele-Shaw equation, we proposed the following numerical scheme which is discrete in time and continuous in space, i.e., seek ϕ n+1 , μ n+1 , ξ n+1 , p n+1 , such that Mathematical Problems in Engineering where the velocity is given by In fact, the velocity u n+1 in equation (10) is an intermediate velocity, and the correction of velocity is the same as the incremental pressure projection method for Navier-Stokes equation (31). So the velocity can be obtained by using the intermediate velocity and satisfies Combining (11) with (12), the Darcy equation will be obtained. By using the divergence operator to (12), the true velocity u n+1 can be eliminated. So, we have

e Fully Discrete Scheme.
Let T h � K { } be a conforming, shape-regular, globally quasi-uniform family of (14) where P k (x, y) is the space of polynomials of degree at most k ∈ Z + . In addition, we define . e fully discrete finite element formulation for the modified Cahn-Hilliard-Hele-Shaw To analyze our scheme, we make the following projection.
Lemma 1 (see Ref. [13]). P is linear, and given w ∈ L 2 (Ω), it follows that In particular, since P(w) ∈ W, Consequently, Definition 2 (see Ref. [13]). Defining where is the unique solution to Clearly P h ∈ W h . Furthermore, we have what as follows.
Lemma 2 (see Ref. [13]). P h is linear, and given any w ∈ L 2 (Ω), it follows that In particular, since P h (w) ∈ W h , Consequently, ere is an estimate for the difference between the projections P and P h . Lemma 3 (see Ref [13]). Suppose that w ∈ H q (Ω) with the compatible boundary conditions w · n � 0 on zΩ and en Lemma 4 (see Ref. [32]). Define the following variational problem, where where (·, ·) − 1,h defines an inner product on _ S h . Hence, for ∀ζ ∈ _ S h and ∀g ∈ S h , Moreover, the following Poincar e ′ -type estimate holds true: for some constants c > 0, independent of h, where

The Energy Stability
h be the unique solution of scheme (15), and define en, the following energy law holds true for any Proof. Utilizing the definition, we get one sees that the scheme (15) can be reformulated as Taking the test function (38) and (39), respectively, and using the identity 2(a, Summing up equations (42)-(44) and using Lemma 4, one has In view of the convexity and concavity about

and through Taylor expansion, we get
where Substituting (47) into (44) and utilizing the identity Taking the inner product of equation (36) en taking the test function q h � τ/cp n h in equation (40) and using the identity 2(a − b, b) Summing up equations (49)-(51), one can conclude that and using the Cauchy-Schwartz inequality, it can be written that Adding (53) into (51), one obtains us, the proof is completed.
h be the unique solution of scheme (15) for constant C > 0 that is independent of h and τ.
Proof. Using the definition of free energy (35), we easily see that (56) holds true. Summing equation (36) from i � 1 to i � n, we get which leads to (56) and (57) immediately.

The Error Estimate
For ϕ, μ, ξ, u , we make the following regularization assumption: where q ≥ 1 is the spatial approximation order given in the space S h . Define the backward difference operator δ τ ϕ n+1 : � ϕ n+1 − ϕ n /τ. For simplicity, the following notation is introduced: Definition 3 (see Ref. [32]). e Ritz projection satisfies the following estimate Lemma 5. Assuming that (ϕ, μ) is weak solution to equation (8), satisfying the regularity assumptions (62), we have the following estimate where C > 0 is a constant independent of h and τ.  (8) and (15), respectively. en, for any h, τ > 0, there exists C > 0 independent of h and τ, such that Proof. Subtracting equation (15) from (8) at t � n + 1, we have , and q h � ϵe n+1 p in equations (65)-(68), respectively, and adding them together and using Lemma 4, we get Mathematical Problems in Engineering where we denote Now, we estimate M i separately. Using Cauchy-Schwarz inequality, Poincar e � inequality, Young's inequality, and Lemma 5, we get (71) Using (33) in Lemma 4, Definition 3, and Young's inequality, we get the following estimate: Homoplastically, according to the Poincar e � inequality and equation (32) in Lemma 4, we obtain and by Definition 3, we can estimate M 3 as About the term M 4 , referring to the concept in Ref. [13], the following estimate is obtained where D: � ‖ϕ n h ‖ 4 L ∞ + 1. For M 5 , according to Lemma 4, Definition 3, Young's inequality, and Taylor extension ‖∇τδ τ ϕ(t)‖ ≤ Cτ 2 , the following inequality is established: where b(ϕ, u, ω): � (ϕu, ∇ω). Following the same idea as in Ref. [13], we know D: � ‖ϕ n h ‖ 4 L ∞ + 1 ≤ C. en, substituting (76) into equation (75) and multiplying the result by 2τ, we get Summing the above estimate (77) from i � 1 to i � n and applying the discrete Gronwall inequality, one can conclude that Mathematical Problems in Engineering us, the theorem is proved.

Numerical Results
In the part of numerical experiment, some numerical examples are given to verify the accuracy of the above theoretical analysis. In the following tests, for the phase field ϕ, chemical potential μ, velocity u, and the pressure p, we use the P1 finite element space. All examples are carried out by using the Freefem++ package [33].  Tables 3 and 4 with θ � 0, 0.0001 and c � 0.01, 0.001, respectively, which confirm the theoretical convergence of O(h).

Spinodal
Decomposition. Now, we study the phase separation dynamics, which is called spinodal decomposition in the modified Cahn-Hilliard-Hele-Shaw system with logarithmic potential. We select test parameters ϵ � 0.02, κ � 0.01, c � 0.01, and h � 1/32 for simulation in the do- e initial condition is taken as the randomly perturbed concentration fields as follows: where rand ∈ [0, 1]. We show the following coarsening process.
We demonstrate evolutions of coarsening dynamics for θ � 0, θ � 10, and θ � 200 with τ � 0.0001 in Figures 1-18. It can be seen that the particles obtained by different θ are similar when T � 0.0001. When θ � 0, the coarsening process gradually become clearer with the increase of T. We can obtain similar results for θ � 10 and θ � 200, so we do not report them here for the sake of brevity.

Energy
Dissipation. Now, we test the energy dissipation for our proposed scheme. e energy function (4) of the modified CH-HS system (1a)-1g can be discreteized as and the modified energy of the fully discrete scheme (15) is defined as We choose the final time T � 0.1, τ � 0.01, and ϵ � 0.1, a � 0.5, c � 0.01, κ � 0.01. e evolution of energy is shown in Figure 19, and it is nonincreasing apparently.
In Figure 19(a), we conclude that the evolution of energy is dissipative with θ � 0. Comparison of the energy with different θ is shown in Figure 19(b). We can easily obtain that the energy decreases faster with larger θ when time step τ � 0.02, which is consistent with the fact that the rate of coarsening dynamics can be improved largely.

Conclusion
is paper presented an efficient and unconditionally stable method for studying the modified Cahn-Hilliard-Hele-Shaw system with logarithmic potential, which incorporated with the first-order backward differentiation formula in time and the mixed finite element method in space. In order to update the pressure, it is significant to solve the linear system at each time step by using an incremental pressure-stabilization technique. By exploiting the convex splitting method, we could treat the nonlinear term in free energy for the logarithmic potential, and the use of the regularization procedure makes the domain for the functional F(ϕ) extended from (− 1, 1) to (− ∞, ∞) so that the problem of the overflow on the boundary can be avoided. en, we not only proved that the first-order discrete scheme is unconditionally stable in energy but also carried out a rigorous error analysis. Various numerical experiments were given to demonstrate the first-order temporal accuracy of numerical algorithm and stability of discrete energy consistent with theoretical analysis. Moreover, we illustrated the capability to capture the realistic phenomena, such as the spinodal decomposition, which revealed the status of internal phase separation.

Data Availability
All data generated or analyzed during this study are included in this published article.  Figure 19: e evolution of discrete energy. 20 Mathematical Problems in Engineering