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

Authors: Xianxia Wang
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

A decoupled unconditionally stable numerical scheme for the modified 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. The fractional step method is used to decouple the nonlinear modified Cahn–Hilliard equation from the pressure gradient. Then, 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.


---

## Body

## 1. Introduction

The Cahn–Hilliard–Hele–Shaw system (CH-HS) [1], which is used to describe the motion of a viscous fluid between two closely spaced parallel plates, is an important mathematical model. The 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 fluid in a Hele–Shaw cell [2] by a similar set of equations. Recently, the CH-HS system has been applied in the study of Saffman–Taylor instability [3] when a more viscous fluid is displaced by a less viscous one resulting in complex pattern formation. This system also serves as a phase field and diffusion interface models of tumor growth [4–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 flow in porous media. The CH-D system has been studied by many authors and one can refer to Refs. [7, 8].LetΩ∈Rd (d=2,3) be an open polygonal or polyhedral domain with a Lipschitz continuous boundary ∂Ω. The modified Cahn–Hilliard–Hele–Shaw system can be studied in the article as follows:(1a)∂tϕ+∇⋅ϕu=ϵΔμ,inΩT:=Ω×0,T,(1b)μ=1ϵfϕ−ϵΔϕ+ξ,inΩT,(1c)−Δξ=θϕ−ϕ¯0,inΩT,(1d)u=−∇p+γϕ∇μ,inΩT,(1e)∇⋅u=0,inΩT,(1f)ϕ|t=0=ϕ0,∀x∈Ω,(1g)u⋅n=0,∂nϕ=∂nμ=0,on∂Ω×0,T,where ϕ is the concentration field, u is the velocity, μ is the chemical potential, ξ is an auxiliary variable, p is the pressure, ϕ¯0=1/Ω∫Ωϕ0xdx, n is the unit outer normal of the boundary ∂Ω, and γ>0, ϵ>0, θ≥0 are physical constants. In particular, the system (1a)–1g is a classical CH-HS system when θ=0 (and further ξ=0, u=0) [9, 10]. Noting that the equations (1a) and (1b) in the system (1a)–1g can be seen as the modified Cahn–Hilliard(CH) equation 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–17].The 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–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]. Therefore, we used the free energy density function of the logarithmic potential in this paper, which is defined as follows [26]:∀a∈0,1,(2)Fϕ=a21+ϕln1+ϕ+1−ϕln1−ϕ+121−ϕ2,ϕ∈−1,1,(3)fϕ=F′ϕ=a2ln1+ϕ+ln1−ϕ−ϕ.The logarithmic free energy functional is defined as(4)Eϕ=∫Ωϵ22∇ϕ2+Fϕdx, and the CH-HS system is mass conservative, i.e., ϕ⋅,t,1=ϕ0,1 and energy dissipative:(5)dEdt=−ϵ∫Ω∇μ2dx−1γ∫Ωu2dx≤0,where ⋅,⋅ is the standard L2 inner product over Ω.Here, we use the Elliott and Luckhaus regularization for problem with the logarithmic free energy functionFϕ replaced by the twice continuously differentiable function Fˇϕ[27, 28], where ∀κ∈0,1:(6)Fˇϕ≔a21+ϕln1+ϕ+a4κ1−ϕ2+a21−ϕlnκ−aκ4+121−ϕ2,ϕ≥1−κ,a21+ϕln1+ϕ+1−ϕln1−ϕ+121−ϕ2,|ϕ|<1−κ,a21−ϕln1−ϕ+a4κ1+ϕ2+a21+ϕlnκ−aκ4+121−ϕ2,ϕ≤−1+κ,,and the monotone function is given by(7)fˇϕ:=Fˇϕ=a2ln1+ϕ+a2−a2κ1−ϕ−a2lnκ−ϕ,ϕ≥1−κ,a2ln1+ϕ−ln1−ϕ−ϕ,|ϕ|<1−κ,−a2ln1−ϕ−a2+a2κ1+ϕ+a2lnκ−ϕ,ϕ≤−1+κ..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. Thus, 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]. The 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.The 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. The first-order scheme based on the convex-splitting method is introduced for the CHHS system. The 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.The rest of the article is organized as follows. The semidiscrete and fully discrete finite element schemes are reviewed in Section2. Then, 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.


## 2. The Discrete Numerical Scheme

### 2.1. The Simi-Discrete Scheme

L2Ω is a square integrable function space, and H1Ω is the Sobolev space. Let us denote u,v=∫Ωuxvxdx, ϕ=ϕL2Ω=ϕ,ϕ1/2, ϕH1Ω=∫ΩDu2dx1/2, and ϕH1Ω=∫Ωϕ2dx+∫ΩDϕ2dx1/2. The weak formulation of the modified Cahn–Hilliard–Hele–Shaw equation can be written as(8)∂tϕ,ω−ϕu,∇ω+ϵ∇μ,∇ω=0,∀ω∈H1Ω,μ,υ−1ϵfϕ,υ−ϵ∇ϕ,∇υ−ξ,υ=0,∀υ∈H1Ω,∇ξ,∇φ−θϕ−ϕ¯0,φ=0,∀φ∈H1Ω,∇p+γϕ∇μ,∇q=0,∀q∈H1Ω,,with the initial condition ϕt=0=ϕ0, where fϕ=f1ϕ−f2ϕ(9)f1ϕ:=a2ln1+ϕ+a2−a2κ1−ϕ−a2lnκ,ϕ≥1−κ,a2ln1+ϕ−ln1−ϕ,ϕ<1−κ,−a2ln1−ϕ−a2+a2κ1+ϕ+a2lnκ,ϕ≤−1+κ,,and f2ϕ:=ϕ.LetN be a positive integer and 0=t0<t1<⋅⋅⋅<tN=T be a uniform partition of 0,T. Let us denote by τ: = tn+1−tn, 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,pn+1, such that(10)ϕn+1−ϕnτ,ω−ϕnun+1,∇ω+ϵ∇μn+1,∇ω=0,∀ω∈H1Ω,μn+1,υ−1ϵf1ϕn+1−f2ϕn,υ−ϵ∇ϕn+1,∇υ−ξn+1,υ=0,∀υ∈H1Ω,∇ξn+1,∇φ−θϕn+1−ϕ¯0,φ=0,∀φ∈H1Ω,∇pn+1−pn,∇q=un+1,∇q,∀q∈H1Ω,,where the velocity is given by(11)un+1=−∇pn+γϕn∇μn+1.In fact, the velocityun+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(12)u˜n+1−un+1=−∇pn+1−∇pn,∇⋅u˜n+1=0.Combining (11) with (12), the Darcy equation will be obtained. By using the divergence operator to (12), the true velocity un+1 can be eliminated. So, we have(13)∇⋅un+1=∇⋅∇pn+1−∇pn.


### 2.2. The Fully Discrete Scheme

LetTh=K be a conforming, shape-regular, globally quasi-uniform family of triangulations of Ω. h=max0≤i≤Nhi, where hi is the grid size. Sh is a piecewise continuous finite element space, which is defined as(14)Sh=vh∈C0ΩvhK∈Pkx,y,K∈Th⊂H1Ω,where Pkx,y is the space of polynomials of degree at most k∈Z+. In addition, we define L02: = ϕ∈L2Ω|ϕ,1=0 and S˙h: = Sh∩L02Ω. The fully discrete finite element formulation for the modified Cahn–Hilliard–Hele–Shaw equation is introduced; for any ωh,υh,φh,qh∈Sh×Sh×Sh×S˙h, find ϕhn+1,μhn+1∈Sh×Sh and ξhn+1,phn+1∈S˙h×S˙h such that(15)ϕhn+1−ϕhnτ,ωh+ϕhn∇phn+γϕhn∇μhn+1,∇ωh+ϵ∇μhn+1,∇ωh=0,μhn+1,υh−1ϵf1ϕhn+1−f2ϕhn,υh−ϵ∇ϕhn+1,∇υh−ξhn+1,υh=0,∇ξhn+1,∇φh−θϕhn+1−ϕ¯0,φh=0,∇phn+1+γϕhn∇μhn+1,∇qh=0,where ϕh0:=Rhϕ0.The operatorRh: ϕ∈H1Ω⟶Sh is the Ritz projection and satisfies(16)∇Rhϕ−ϕ,∇χ=0,∀χ∈Sh,Rhϕ−ϕ,1=0.To analyze our scheme, we make the following projection.Definition 1 (see Ref. [13]).
Define(17)W:=u∈L2Ω|u,∇q=0,∀q∈H1Ω,and the projection P:=L2Ω⟶W is defined via(18)Pw=∇p+w,where p∈H˙1Ω:=ϕ∈H1Ω|ϕ,1=0 is the unique solution to(19)∇p+w,∇q=0,∀q∈H1Ω.
Clearly,Pw∈W for any w∈L2Ω. Furthermore, we have what as follows.Lemma 1 (see Ref. [13]).
P is linear, and given w∈L2Ω, it follows that(20)Pw−w,v=0,∀v∈W.
In particular, sincePw∈W,(21)Pw−w,Pw=0.
Consequently,(22)Pw≤w,∀w∈L2Ω.Definition 2 (see Ref. [13]).
Defining(23)Wh:=uh∈L2Ω|uh,∇qh=0,∀qh∈Sh,we observe that W∈Wh. The projection Ph:=L2Ω⟶Wh is defined via(24)Phw=∇ph+w,where ph∈Sh. is the unique solution to(25)∇ph+w,∇qh=0,∀qh∈Sh.
ClearlyPh∈Wh. Furthermore, we have what as follows.Lemma 2 (see Ref. [13]).
Ph is linear, and given any w∈L2Ω, it follows that(26)Phw−w,v=0,∀v∈Wh.
In particular, sincePhw∈Wh,(27)Phw−w,Phw=0.Consequently,(28)Phw≤w,∀w∈L2Ω.
There is an estimate for the difference between the projectionsP and Ph.Lemma 3 (see Ref [13]).
Suppose thatw∈HqΩ with the compatible boundary conditions w⋅n=0 on ∂Ω and p∈Hq+1Ω, where(29)∇p=Pw−w.
Then(30)Phw−Pw≤ChqpHq+1.Lemma 4 (see Ref. [32]).
Define the following variational problem, givenζ∈S˙h, find Thζ∈S˙h such that(31)∇Thζ,∇χ=ζ,χ,where Th:S˙h⟶S˙h is an invertible linear operator. Suppose ζ,ψ∈S˙h and set(32)ζ,ψ−1,h:=∇Thζ,∇Thψ=ζ,Thψ=Thζ,ψ,where ⋅,⋅−1,h defines an inner product on S˙h. Hence, for ∀ζ∈S˙h and ∀g∈Sh,(33)ζ,g≤ζ−1,h∇g.
Moreover, the following Poincare′-type estimate holds true: for some constants c>0, independent of h, where(34)ζ−1,h≤cζ,∀ζ∈S˙h.


## 2.1. The Simi-Discrete Scheme

L2Ω is a square integrable function space, and H1Ω is the Sobolev space. Let us denote u,v=∫Ωuxvxdx, ϕ=ϕL2Ω=ϕ,ϕ1/2, ϕH1Ω=∫ΩDu2dx1/2, and ϕH1Ω=∫Ωϕ2dx+∫ΩDϕ2dx1/2. The weak formulation of the modified Cahn–Hilliard–Hele–Shaw equation can be written as(8)∂tϕ,ω−ϕu,∇ω+ϵ∇μ,∇ω=0,∀ω∈H1Ω,μ,υ−1ϵfϕ,υ−ϵ∇ϕ,∇υ−ξ,υ=0,∀υ∈H1Ω,∇ξ,∇φ−θϕ−ϕ¯0,φ=0,∀φ∈H1Ω,∇p+γϕ∇μ,∇q=0,∀q∈H1Ω,,with the initial condition ϕt=0=ϕ0, where fϕ=f1ϕ−f2ϕ(9)f1ϕ:=a2ln1+ϕ+a2−a2κ1−ϕ−a2lnκ,ϕ≥1−κ,a2ln1+ϕ−ln1−ϕ,ϕ<1−κ,−a2ln1−ϕ−a2+a2κ1+ϕ+a2lnκ,ϕ≤−1+κ,,and f2ϕ:=ϕ.LetN be a positive integer and 0=t0<t1<⋅⋅⋅<tN=T be a uniform partition of 0,T. Let us denote by τ: = tn+1−tn, 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,pn+1, such that(10)ϕn+1−ϕnτ,ω−ϕnun+1,∇ω+ϵ∇μn+1,∇ω=0,∀ω∈H1Ω,μn+1,υ−1ϵf1ϕn+1−f2ϕn,υ−ϵ∇ϕn+1,∇υ−ξn+1,υ=0,∀υ∈H1Ω,∇ξn+1,∇φ−θϕn+1−ϕ¯0,φ=0,∀φ∈H1Ω,∇pn+1−pn,∇q=un+1,∇q,∀q∈H1Ω,,where the velocity is given by(11)un+1=−∇pn+γϕn∇μn+1.In fact, the velocityun+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(12)u˜n+1−un+1=−∇pn+1−∇pn,∇⋅u˜n+1=0.Combining (11) with (12), the Darcy equation will be obtained. By using the divergence operator to (12), the true velocity un+1 can be eliminated. So, we have(13)∇⋅un+1=∇⋅∇pn+1−∇pn.


## 2.2. The Fully Discrete Scheme

LetTh=K be a conforming, shape-regular, globally quasi-uniform family of triangulations of Ω. h=max0≤i≤Nhi, where hi is the grid size. Sh is a piecewise continuous finite element space, which is defined as(14)Sh=vh∈C0ΩvhK∈Pkx,y,K∈Th⊂H1Ω,where Pkx,y is the space of polynomials of degree at most k∈Z+. In addition, we define L02: = ϕ∈L2Ω|ϕ,1=0 and S˙h: = Sh∩L02Ω. The fully discrete finite element formulation for the modified Cahn–Hilliard–Hele–Shaw equation is introduced; for any ωh,υh,φh,qh∈Sh×Sh×Sh×S˙h, find ϕhn+1,μhn+1∈Sh×Sh and ξhn+1,phn+1∈S˙h×S˙h such that(15)ϕhn+1−ϕhnτ,ωh+ϕhn∇phn+γϕhn∇μhn+1,∇ωh+ϵ∇μhn+1,∇ωh=0,μhn+1,υh−1ϵf1ϕhn+1−f2ϕhn,υh−ϵ∇ϕhn+1,∇υh−ξhn+1,υh=0,∇ξhn+1,∇φh−θϕhn+1−ϕ¯0,φh=0,∇phn+1+γϕhn∇μhn+1,∇qh=0,where ϕh0:=Rhϕ0.The operatorRh: ϕ∈H1Ω⟶Sh is the Ritz projection and satisfies(16)∇Rhϕ−ϕ,∇χ=0,∀χ∈Sh,Rhϕ−ϕ,1=0.To analyze our scheme, we make the following projection.Definition 1 (see Ref. [13]).
Define(17)W:=u∈L2Ω|u,∇q=0,∀q∈H1Ω,and the projection P:=L2Ω⟶W is defined via(18)Pw=∇p+w,where p∈H˙1Ω:=ϕ∈H1Ω|ϕ,1=0 is the unique solution to(19)∇p+w,∇q=0,∀q∈H1Ω.
Clearly,Pw∈W for any w∈L2Ω. Furthermore, we have what as follows.Lemma 1 (see Ref. [13]).
P is linear, and given w∈L2Ω, it follows that(20)Pw−w,v=0,∀v∈W.
In particular, sincePw∈W,(21)Pw−w,Pw=0.
Consequently,(22)Pw≤w,∀w∈L2Ω.Definition 2 (see Ref. [13]).
Defining(23)Wh:=uh∈L2Ω|uh,∇qh=0,∀qh∈Sh,we observe that W∈Wh. The projection Ph:=L2Ω⟶Wh is defined via(24)Phw=∇ph+w,where ph∈Sh. is the unique solution to(25)∇ph+w,∇qh=0,∀qh∈Sh.
ClearlyPh∈Wh. Furthermore, we have what as follows.Lemma 2 (see Ref. [13]).
Ph is linear, and given any w∈L2Ω, it follows that(26)Phw−w,v=0,∀v∈Wh.
In particular, sincePhw∈Wh,(27)Phw−w,Phw=0.Consequently,(28)Phw≤w,∀w∈L2Ω.
There is an estimate for the difference between the projectionsP and Ph.Lemma 3 (see Ref [13]).
Suppose thatw∈HqΩ with the compatible boundary conditions w⋅n=0 on ∂Ω and p∈Hq+1Ω, where(29)∇p=Pw−w.
Then(30)Phw−Pw≤ChqpHq+1.Lemma 4 (see Ref. [32]).
Define the following variational problem, givenζ∈S˙h, find Thζ∈S˙h such that(31)∇Thζ,∇χ=ζ,χ,where Th:S˙h⟶S˙h is an invertible linear operator. Suppose ζ,ψ∈S˙h and set(32)ζ,ψ−1,h:=∇Thζ,∇Thψ=ζ,Thψ=Thζ,ψ,where ⋅,⋅−1,h defines an inner product on S˙h. Hence, for ∀ζ∈S˙h and ∀g∈Sh,(33)ζ,g≤ζ−1,h∇g.
Moreover, the following Poincare′-type estimate holds true: for some constants c>0, independent of h, where(34)ζ−1,h≤cζ,∀ζ∈S˙h.


## 3. The Energy Stability

Theorem 1.
Letϕhn+1,μhn+1,ξhn+1,phn+1 be the unique solution of scheme (15), and define(35)Ξϕhn+1:=Eϕhn+1+θ2ϕhn+1−ϕ¯0−1,h2+τ2γ∇phn+12.
Then, the following energy law holds true for anyτ,h,ϵ>0, where(36)Ξϕhn+1+ε2∇ϕhn+1−ϕhn2+θ2ϕhn+1−ϕhn−1,h2,+τε∇μhn+12+τ2γuhn+12≤Ξϕhn.Proof.
Utilizing the definition, we get(37)uhn+1=−∇phn+γϕhn∇μhn+1,one sees that the scheme (15) can be reformulated as(38)ϕhn+1−ϕhnτ,ωh−ϕhnuhn+1,∇ωh+ϵ∇μhn+1,∇ωh=0,(39)μhn+1,υh−1ϵf1ϕhn+1−f2ϕhn,υh−ϵ∇ϕhn+1,∇υh−ξhn+1,υh=0,(40)∇ξhn+1,∇φh−θϕhn+1−ϕ¯0,φh=0,(41)∇phn+1−phn,∇qh=uhn+1,∇qh.
Taking the test functionωh=τμhn+1, υh=−ϕhn+1−ϕhn, φh=−Thϕhn+1−ϕhn in equations (38) and (39), respectively, and using the identity 2a,a−b=a2−b2+a−b2, we obtain(42)ϕhn+1−ϕhn,μhn+1+τϵ∇μhn+12−τϕhnuhn+1,∇μhn+1=0,(43)−μhn+1,ϕhn+1−ϕhn+1ϵf1ϕhn+1−f2ϕhn,ϕhn+1−ϕhn,+ϵ2∇ϕhn+12−∇ϕhn2+∇ϕhn+1−ϕhn2+ξhn+1,ϕhn+1−ϕhn=0,(44)−∇ξhn+1,∇Thϕhn+1−ϕhn+θϕhn+1−ϕ¯0,Thϕhn+1−ϕhn=0.
Summing up equations (42)–(44) and using Lemma 4, one has(45)τϵ∇μhn+12−τϕhnuhn+1,∇μhn+1+1ϵf1ϕhn+1−f2ϕhn,ϕhn+1−ϕhn,+ϵ2∇ϕhn+12−∇ϕhn2+∇ϕhn+1−ϕhn2+θϕhn+1−ϕ¯0,ϕhn+1−ϕhn−1,h=0.
In view of the convexity and concavity aboutfϕ=f1ϕ−f2ϕ, for all ϕ, there are f1′ϕ>0, f2′ϕ=1>0, and through Taylor expansion, we get(46)F1ϕhn−F1ϕhn+1=f1ϕhn+1ϕhn−ϕhn+1+12f1′η1ϕhn−ϕhn+12,≥f1ϕhn+1ϕhn−ϕhn+1,(47)F2ϕhn+1−F2ϕhn=f2ϕhnϕhn+1−ϕhn+12f2′η2ϕhn+1−ϕhn2,≥f2ϕhnϕhn+1−ϕhn,where η1∈ϕhn,ϕhn+1, η2∈ϕhn,ϕhn+1, we have(48)Fϕhn+1−Fϕhn,1≤f1ϕhn+1−f2ϕhn,ϕhn+1−ϕhn.
Substituting (47) into (44) and utilizing the identity 2a,a−b=a2−b2+a−b2 gives(49)τϵ∇μhn+12−τϕhnuhn+1,∇μhn+1+1ϵFϕhn+1−Fϕhn,1,+θ2ϕhn+1−ϕ¯0−1,h2−ϕhn−ϕ¯0−1,h2+ϕhn+1−ϕhn−1,h2,+ϵ2∇ϕhn+12−∇ϕhn2+∇ϕhn+1−ϕhn2≤0.
Taking the inner product of equation (36) with τ/γuhn+1, we get(50)τγuhn+12+τγ∇phn,uhn+1=−τϕhn∇μhn+1,uhn+1.
Then taking the test functionqh=τ/γphn in equation (40) and using the identity 2a−b,b=a2−b2−a−b2, we get(51)τ2γ∇phn+12−∇phn2−∇phn+1−∇phn2=τγ∇phn,uhn+1.
Summing up equations (49)–(51), one can conclude that(52)τϵ∇μhn+12+1ϵFϕhn+1−Fϕhn,1+ϵ2∇ϕhn+12−∇ϕhn2+∇ϕhn+1−ϕhn2,+θ2ϕhn+1−ϕ¯0−1,h2−ϕhn−ϕ¯0−1,h2+ϕhn+1−ϕhn−1,h2+τγuhn+12,+τ2γ∇phn+12−∇phn2−∇phn+1−∇phn2≤0.
Then, resettingqh=τ/2γphn+1−phn in (40), we get(53)τ2γ∇phn+1−phn2=τ2γ∇phn+1−phn,uhn+1,and using the Cauchy–Schwartz inequality, it can be written that(54)τ2γ∇phn+1−phn2≤τ2γuhn+12.
Adding (53) into (51), one obtains(55)τϵ∇μhn+12+1ϵFϕhn+1−Fϕhn,1+ϵ2∇ϕhn+12−∇ϕhn2+∇ϕhn+1−ϕhn2,+θ2ϕhn+1−ϕ¯0−1,h2−ϕhn−ϕ¯0−1,h2+ϕhn+1−ϕhn−1,h2+τ2γuhn+12,+τ2γ∇phn+12−∇phn2≤0.
Thus, the proof is completed.Corollary 1.
Letϕhn+1,μhn+1,ξhn+1,phn+1 be the unique solution of scheme (15). Assume that Ξϕh0≤C0, we have(56)max0<i≤n∇ϕhi+12+ϕhi+1−ϕ¯0−1,h2+∇phi+12≤C,(57)∑i=0n∇ϕhi+1−ϕhi2+ϕhi+1−ϕhi−1,h2≤C,(58)τ∑i=0n∇μhi+12+uhi+12≤C,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(59)Ξϕhn+1+ϵ2∑i=0n∇ϕhi+1−ϕhi2+θ2∑i=0nϕhi+1−ϕhi−1,h2,+τϵ∑i=0n∇μhi+12+τ2γ∑i=0nuhi+12≤Ξϕh0≤C,which leads to (56) and (57) immediately.


## 4. The Error Estimate

Forϕ,μ,ξ,u, we make the following regularization assumption:(60)ϕ∈H10,T;Hq+1Ω∩L∞0,T;H1Ω∩L∞0,T;Hq+1Ω,μ∈L∞0,T;H1Ω∩L20,T;Hq+1Ω,ξ∈L20,T;Hq+1Ω,u∈L∞0,T;HqΩ,ϕ∇μ∈L∞0,T;HqΩ,where q≥1 is the spatial approximation order given in the space Sh.Define the backward difference operatorδτϕn+1:=ϕn+1−ϕn/τ. For simplicity, the following notation is introduced:(61)e^ϕn+1:=Rhϕn+1−ϕhn+1,e˜ϕn+1:=ϕn+1−Rhϕn+1,e^μn+1:=Rhμn+1−μhn+1,e˜μn+1:=μn+1−Rhμn+1,e^ξn+1:=Rhξn+1−ξhn+1,e˜ξn+1:=ξn+1−Rhξn+1,e^pn+1:=Rhpn+1−phn+1,e˜pn+1:=pn+1−Rhpn+1,σϕn+1:=δτRhϕn+1−∂tϕn+1.Definition 3 (see Ref. [32]).
The Ritz projection satisfies the following estimate(62)ϕ−Rhϕ+hϕ−RhϕH1Ω≤Chq+1ϕHq+1Ω.Lemma 5.
Assuming thatϕ,μ is weak solution to equation (8), satisfying the regularity assumptions (62), we have the following estimate(63)σϕn+12≤Ch2q+2+Cτ2,∀t∈τ,T,where C>0 is a constant independent of h and τ.Theorem 2.
Assume thatϕ,μ,ξ,p and ϕhn+1,μhn+1,ξhn+1,phn+1 are the solutions of (8) and (15), respectively. Then, for any h,τ>0, there exists C>0 independent of h and τ, such that(64)∑i=0n2τϵ∇e^μi+12+ϵ∇e^ϕn+12+θ∇e^ϕn+1−1,h2,+∑i=0n2τϵ∇e^pi+12+∑i=0n4τγPhϕhi∇e^μi+12≤Cτ2+Ch2q.Proof.
Subtracting equation (15) from (8) at t=n+1, we have(65)δτe^ϕn+1,ωh+ϵ∇e^μn+1,∇ωh−σϕn+1,ωh,+ϕn+1∇pn+1+γϕn+1∇μn+1,∇ωh−ϕhn∇phn+γϕhn∇μhn+1,∇ωh=0,(66)e^μn+1,υh+e˜μn+1,υh−1ϵf1ϕn+1−f2ϕn+1,υh+1ϵf1ϕhn+1−f2ϕhn,υh,−ϵ∇e^ϕn+1,∇υh−e˜ξn+1,υh−e^ξn+1,υh=0,(67)∇e^ξn+1,∇φh+∇e˜ξn+1,∇φh−θe^ϕn+1,φh−θe˜ϕn+1,φh=0,(68)∇e^pn+1,∇qh+γϕn+1∇μn+1−ϕhn∇μhn+1,∇qh=0.
Settingωh=e^μn+1, υh=−δτe^ϕn+1, φh=−Thδτe^ϕn+1, and qh=ϵe^pn+1 in equations (65)–(68), respectively, and adding them together and using Lemma 4, we get(69)ϵ∇e^μn+12+ϵ2τ∇e^ϕn+12−∇e^ϕn2+∇e^ϕn+1−∇e^ϕn2,+ϵ∇e^pn+12+θ2τe^ϕn+1−1,h2−e^ϕn−1,h2+e^ϕn+1−e^ϕn−1,h2,=σϕn+1,e^μn+1+e˜μn+1,δτe^ϕn+1−θe˜ϕn+1,Thδτe^ϕn+1,−ϕn+1∇pn+1+γϕn+1∇μn+1,∇e^μn+1+ϕhn∇phn+γϕhn∇μhn+1,∇e^μn+1,−1ϵf1ϕn+1−f2ϕn+1,δτe^ϕn+1+1ϵf1ϕhn+1−f2ϕhn,δτe^ϕn+1,−ϵγϕn+1∇μn+1−ϕhn∇μhn+1,∇e^pn+1,=M1+M2+M3+M4+M5+M6,where we denote(70)M1=σϕn+1,e^μn+1,M2=e˜μn+1,δτe^ϕn+1,M3=−θe˜ϕn+1,Thδτe^ϕn+1,M4=−ϕn+1∇pn+1+γϕn+1∇μn+1,∇e^μn+1+ϕhn∇phn+γϕhn∇μhn+1,∇e^μn+1,M5=−1ϵf1ϕn+1−f2ϕn+1,δτe^ϕn+1+1ϵf1ϕhn+1−f2ϕhn,δτe^ϕn+1,M6=−ϵγϕn+1∇μn+1−ϕhn∇μhn+1,∇e^pn+1.
Now, we estimateMi separately. Using Cauchy–Schwarz inequality, Poincar e´ inequality, Young’s inequality, and Lemma 5, we get(71)M1≤σϕn+1,e^μn+1,=σϕn+1,e^μn+1−e^μn+1¯,≤σϕn+1e^μn+1−e^μn+1¯,≤σϕn+1∇e^μn+1,≤2σϕn+12+18∇e^μn+12,≤Ch2q+2+Cτ2+C∇e^μn+12.
Using (33) in Lemma 4, Definition 3, and Young’s inequality, we get the following estimate:(72)M2≤e˜μn+1,δτe^ϕn+1,≤∇e˜μn+1δτe^ϕn+1−1,h,≤Ch2q+Cδτe^ϕn+1−1,h2.
Homoplastically, according to the Poincare´ inequality and equation (32) in Lemma 4, we obtain(73)Thδτe^ϕn+12≤C∇Thδτe^ϕn+12=Cδτe^ϕn+1−1,h2,and by Definition 3, we can estimate M3 as (74)M3≤θe˜ϕn+1,Thδτe^ϕn+1,≤θe˜ϕn+1Thδτe^ϕn+1,≤θ2e˜ϕn+12+14δτe^ϕn+1−1,h2,≤Ch2q+2+Cδτe^ϕn+1−1,h2.
About the termM4, referring to the concept in Ref. [13], the following estimate is obtained(75)M4≤CDτ2+h2q+ϵ4∇e^μn+12+CD∇e^ϕn2−γPhϕhn∇e^μn+12,where D:=ϕhnL∞4+1.
ForM5, according to Lemma 4, Definition 3, Young’s inequality, and Taylor extension ∇τδτϕt≤Cτ2, the following inequality is established:(76)M5=−1ϵf1ϕn+1−f2ϕn+1,δτe^ϕn+1+1ϵf1ϕhn+1−f2ϕhn,δτe^ϕn+1,=−1ϵf1ϕn+1−f1ϕhn+1,δτe^ϕn+1+1ϵf2ϕn+1−f2ϕhn,δτe^ϕn+1,≤1ϵf1′ςϕn+1−ϕhn+1,δτe^ϕn+1+1ϵϕn+1−ϕhn,δτe^ϕn+1,≤Lϵ∇ϕn+1−ϕhn+1δτe^ϕn+1−1,h+1ϵ∇ϕn+1−ϕhnδτe^ϕn+1−1,h,≤C∇ϕn+1−ϕhn+12+14δτe^ϕn+1−1,h2+C∇ϕn+1−ϕhn2+14δτe^ϕn+1−1,h2,≤Ch2q+Cτ2+C∇e^ϕn+12+C∇e^ϕn2+Cδτe^ϕn+1−1,h2,where we note f1′ς=L.
Using Definitions1 and 2, Lemma 3, Taylor extension ∇τδτpn+1≤Cτ2, the Cauchy–Schwarz inequality, and Young’s inequality, one obtains(77)M6=−ϵγϕn+1∇μn+1−ϕhn∇μhn+1,∇e^pn+1,=−ϵγ∇pn+1+ϕn+1∇μn+1,∇e^pn+1+ϵγ∇pn+1−∇pn,∇e^pn+1,+ϵγ∇pn−∇phn,∇e^pn+1+ϵγ∇phn+ϕhn∇μhn+1,∇e^pn+1,=ϵγ∇pn+1−∇pn,∇e^pn+1+ϵγ∇pn−∇phn,∇e^pn+1,≤C∇pn+1−pn2+14∇e^pn+12+C∇pn−phn2+14∇e^pn+12,≤Cτ2+Ch2q+C∇e^pn+12.
Combining equations (69)–(75), we have(78)ϵ∇e^μn+12+ϵ2τ∇e^ϕn+12−∇e^ϕn2+∇e^ϕn+1−∇e^ϕn2+ϵ∇e^pn+12,+θ2τe^ϕn+1−1,h2−e^ϕn−1,h2+e^ϕn+1−e^ϕn−1,h2+γPhϕhn∇e^μn+12,≤Cτ2+Ch2q+CDτ2+h2q+C∇e^μn+12+C∇e^ϕn+12+C∇e^ϕn2+CD∇e^ϕn2+C∇e^pn+12+Cδτe^ϕn+1−1,h2.
The next step is to estimate the normδτe^ϕn+1−1,h2. Setting ωh=Thδτe^ϕn+1 in equation (63) and using the same concept as M4, we obtain(79)δτe^ϕn+1−1,h2=−ϵ∇e^μn+1,∇Thδτe^ϕn+1+σϕn+1,Thδτe^ϕn+1,+bϕn+1,un+1,Thδτe^ϕn+1−bϕhn,uhn+1,Thδτe^ϕn+1,≤ϵ∇e^μn+1∇Thδτe^ϕn+1+σϕn+1Thδτe^ϕn+1,+be˜ϕn+1,un+1,Thδτe^ϕn+1+bτδτRhϕn+1,un+1,Thδτe^ϕn+1,+be^ϕn,un+1,Thδτe^ϕn+1+bϕhn,un+1−uhn+1,Thδτe^ϕn+1,≤ϵ2∇e^μn+12+14δτe^ϕn+1−1,h2+σϕn+12+14δτe^ϕn+1−1,h2,+CDτ2+h2q−γPhϕhn∇e^μn+12+CD∇e^ϕn2+ϵ4δτe^ϕn+1−1,h2,≤Ch2q+2+Cτ2+CDτ2+h2q−γPhϕhn∇e^μn+12,+C∇e^μn+12+CD∇e^ϕn2+Cδτe^ϕn+1−1,h2,where bϕ,u,ω:=ϕu,∇ω.
Following the same idea as in Ref. [13], we know D:=‖ϕhn‖L∞4+1≤C. Then, substituting (76) into equation (75) and multiplying the result by 2τ, we get(80)2τϵ∇e^μn+12+ϵ∇e^ϕn+12−∇e^ϕn2+∇e^ϕn+1−∇e^ϕn2+2τϵ∇e^pn+12+θe^ϕn+1−1,h2−e^ϕn−1,h2+e^ϕn+1−e^ϕn−1,h2+4τγPhϕhn∇e^μn+12≤Cττ2+Cτh2q+C∇e^μn+12+C∇e^ϕn+12+C∇e^ϕn2+C∇e^pn+12.
Summing the above estimate (77) from i=1 to i=n and applying the discrete Gronwall inequality, one can conclude that(81)∑i=0n2τϵ∇e^μi+12+ϵ∇e^ϕn+12+θ∇e^ϕn+1−1,h2,+∑i=0n2τϵ∇e^pi+12+∑i=0n4τγPhϕhi∇e^μi+12≤Cτ2+Ch2q.
Thus, the theorem is proved.


## 5. 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].


### 5.1. Convergence and Energy Dissipation

We consider the modified Cahn–Hilliard–Hele–Shaw system given in domain0,1×0,1 with the initial condition(82)ϕ0=0.24cos2πxcos2πy+0.4cosπxcos3πy.


#### 5.1.1. The Order of Spatial Convergence

Here, we discuss the spatial convergence order ofϕ. The fixed test parameters ϵ=0.1, a=0.01, κ=0.01, τ=0.01, T=0.1, and the varying grid step size h=1/16,1/32,1/64,1/128 are used here. The errors and convergence rates in the H1 norm are presented in Tables 1 and 2 with θ=0,0.00005 and γ=0.1,0.001, respectively, which confirm the theoretical convergence of Oh.Table 1
The errors and convergence rates ofϕ with θ=0.
γ=0.1he^ϕH1rateγ=0.001e^ϕH1Rate1/161.17651.132231/320.5192681.179940.4966591.188841/640.259730.9994670.2531560.972231/1280.1288221.011640.126361.00249Table 2
The errors and convergence rates ofϕ with θ=0.00005.
γ=0.1he^ϕH1Rateγ=0.001e^ϕH1Rate1/161.176531.132291/320.5192591.180010.4966451.188961/640.2597320.9994270.2531890.9720011/1280.1288391.011460.1263121.00323


#### 5.1.2. The Order of Temporal Convergence

Similarly, the temporal convergence order ofϕ is discussed as follows. The test parameters are chosen as ϵ=0.01, a=0.01, κ=0.01, τ=0.01, T=0.1, and time step size τ=h=1/64,1/128,1/256. The errors and convergence rates in the H1 norm are presented in Tables 3 and 4 with θ=0,0.0001 and γ=0.01,0.001, respectively, which confirm the theoretical convergence of Oh.Table 3
The errors and convergence rates ofϕ with θ=0.
γ=0.01τe^ϕH1Rateγ=0.001e^ϕH1Rate1/642.549112.223371/1281.182251.108461.009541.139061/2560.5813861.023970.4878481.04919Table 4
The errors and convergence rates ofϕ with θ=0.0001.
γ=0.01τe^ϕH1Rateγ=0.001e^ϕH1Rate1/642.549132.223361/1281.182221.10851.009541.139051/2560.5813851.023930.4878461.04919


### 5.2. 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, γ=0.01, and h=1/32 for simulation in the domain 0,1×0,1. The initial condition is taken as the randomly perturbed concentration fields as follows:(83)ϕ0=2∗rand−1,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.Figure 1
T = 0.0001, θ = 0.Figure 2
T = 0.0001, θ = 10.Figure 3
T = 0.0001, θ = 200.Figure 4
T = 0.0005, θ = 0.Figure 5
T = 0.0005, θ = 10.Figure 6
T = 0.0005, θ = 200.Figure 7
T = 0.001, θ = 0.Figure 8
T = 0.001, θ = 10.Figure 9
T = 0.001, θ = 200.Figure 10
T = 0.005, θ = 0.Figure 11
T = 0.005, θ = 10.Figure 12
T = 0.005, θ = 200.Figure 13
T = 0.01, θ = 0.Figure 14
T = 0.01, θ = 10.Figure 15
T = 0.01, θ = 200.Figure 16
T = 0.02, θ = 0.Figure 17
T = 0.02, θ = 10.Figure 18
T = 0.02, θ = 200.


#### 5.2.1. Energy Dissipation

Now, we test the energy dissipation for our proposed scheme. The energy function (4) of the modified CH-HS system (1a)–1g can be discreteized as(84)Eϕhn+1=∫Ωϵ22∇ϕhn+12+Fϕhn+1dx,and the modified energy of the fully discrete scheme (15) is defined as(85)Ξϕhn+1:=Eϕhn+1+θ2ϕhn+1−ϕ¯0−1,h2+τ2γ∇phn+12.We choose the final timeT=0.1, τ=0.01, and ϵ=0.1, a=0.5, γ=0.01, κ=0.01. The evolution of energy is shown in Figure 19, and it is nonincreasing apparently.Figure 19
The evolution of discrete energy.
(a)(b)In Figure19(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.


## 5.1. Convergence and Energy Dissipation

We consider the modified Cahn–Hilliard–Hele–Shaw system given in domain0,1×0,1 with the initial condition(82)ϕ0=0.24cos2πxcos2πy+0.4cosπxcos3πy.


### 5.1.1. The Order of Spatial Convergence

Here, we discuss the spatial convergence order ofϕ. The fixed test parameters ϵ=0.1, a=0.01, κ=0.01, τ=0.01, T=0.1, and the varying grid step size h=1/16,1/32,1/64,1/128 are used here. The errors and convergence rates in the H1 norm are presented in Tables 1 and 2 with θ=0,0.00005 and γ=0.1,0.001, respectively, which confirm the theoretical convergence of Oh.Table 1
The errors and convergence rates ofϕ with θ=0.
γ=0.1he^ϕH1rateγ=0.001e^ϕH1Rate1/161.17651.132231/320.5192681.179940.4966591.188841/640.259730.9994670.2531560.972231/1280.1288221.011640.126361.00249Table 2
The errors and convergence rates ofϕ with θ=0.00005.
γ=0.1he^ϕH1Rateγ=0.001e^ϕH1Rate1/161.176531.132291/320.5192591.180010.4966451.188961/640.2597320.9994270.2531890.9720011/1280.1288391.011460.1263121.00323


### 5.1.2. The Order of Temporal Convergence

Similarly, the temporal convergence order ofϕ is discussed as follows. The test parameters are chosen as ϵ=0.01, a=0.01, κ=0.01, τ=0.01, T=0.1, and time step size τ=h=1/64,1/128,1/256. The errors and convergence rates in the H1 norm are presented in Tables 3 and 4 with θ=0,0.0001 and γ=0.01,0.001, respectively, which confirm the theoretical convergence of Oh.Table 3
The errors and convergence rates ofϕ with θ=0.
γ=0.01τe^ϕH1Rateγ=0.001e^ϕH1Rate1/642.549112.223371/1281.182251.108461.009541.139061/2560.5813861.023970.4878481.04919Table 4
The errors and convergence rates ofϕ with θ=0.0001.
γ=0.01τe^ϕH1Rateγ=0.001e^ϕH1Rate1/642.549132.223361/1281.182221.10851.009541.139051/2560.5813851.023930.4878461.04919


## 5.1.1. The Order of Spatial Convergence

Here, we discuss the spatial convergence order ofϕ. The fixed test parameters ϵ=0.1, a=0.01, κ=0.01, τ=0.01, T=0.1, and the varying grid step size h=1/16,1/32,1/64,1/128 are used here. The errors and convergence rates in the H1 norm are presented in Tables 1 and 2 with θ=0,0.00005 and γ=0.1,0.001, respectively, which confirm the theoretical convergence of Oh.Table 1
The errors and convergence rates ofϕ with θ=0.
γ=0.1he^ϕH1rateγ=0.001e^ϕH1Rate1/161.17651.132231/320.5192681.179940.4966591.188841/640.259730.9994670.2531560.972231/1280.1288221.011640.126361.00249Table 2
The errors and convergence rates ofϕ with θ=0.00005.
γ=0.1he^ϕH1Rateγ=0.001e^ϕH1Rate1/161.176531.132291/320.5192591.180010.4966451.188961/640.2597320.9994270.2531890.9720011/1280.1288391.011460.1263121.00323


## 5.1.2. The Order of Temporal Convergence

Similarly, the temporal convergence order ofϕ is discussed as follows. The test parameters are chosen as ϵ=0.01, a=0.01, κ=0.01, τ=0.01, T=0.1, and time step size τ=h=1/64,1/128,1/256. The errors and convergence rates in the H1 norm are presented in Tables 3 and 4 with θ=0,0.0001 and γ=0.01,0.001, respectively, which confirm the theoretical convergence of Oh.Table 3
The errors and convergence rates ofϕ with θ=0.
γ=0.01τe^ϕH1Rateγ=0.001e^ϕH1Rate1/642.549112.223371/1281.182251.108461.009541.139061/2560.5813861.023970.4878481.04919Table 4
The errors and convergence rates ofϕ with θ=0.0001.
γ=0.01τe^ϕH1Rateγ=0.001e^ϕH1Rate1/642.549132.223361/1281.182221.10851.009541.139051/2560.5813851.023930.4878461.04919


## 5.2. 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, γ=0.01, and h=1/32 for simulation in the domain 0,1×0,1. The initial condition is taken as the randomly perturbed concentration fields as follows:(83)ϕ0=2∗rand−1,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.Figure 1
T = 0.0001, θ = 0.Figure 2
T = 0.0001, θ = 10.Figure 3
T = 0.0001, θ = 200.Figure 4
T = 0.0005, θ = 0.Figure 5
T = 0.0005, θ = 10.Figure 6
T = 0.0005, θ = 200.Figure 7
T = 0.001, θ = 0.Figure 8
T = 0.001, θ = 10.Figure 9
T = 0.001, θ = 200.Figure 10
T = 0.005, θ = 0.Figure 11
T = 0.005, θ = 10.Figure 12
T = 0.005, θ = 200.Figure 13
T = 0.01, θ = 0.Figure 14
T = 0.01, θ = 10.Figure 15
T = 0.01, θ = 200.Figure 16
T = 0.02, θ = 0.Figure 17
T = 0.02, θ = 10.Figure 18
T = 0.02, θ = 200.


### 5.2.1. Energy Dissipation

Now, we test the energy dissipation for our proposed scheme. The energy function (4) of the modified CH-HS system (1a)–1g can be discreteized as(84)Eϕhn+1=∫Ωϵ22∇ϕhn+12+Fϕhn+1dx,and the modified energy of the fully discrete scheme (15) is defined as(85)Ξϕhn+1:=Eϕhn+1+θ2ϕhn+1−ϕ¯0−1,h2+τ2γ∇phn+12.We choose the final timeT=0.1, τ=0.01, and ϵ=0.1, a=0.5, γ=0.01, κ=0.01. The evolution of energy is shown in Figure 19, and it is nonincreasing apparently.Figure 19
The evolution of discrete energy.
(a)(b)In Figure19(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.


## 5.2.1. Energy Dissipation

Now, we test the energy dissipation for our proposed scheme. The energy function (4) of the modified CH-HS system (1a)–1g can be discreteized as(84)Eϕhn+1=∫Ωϵ22∇ϕhn+12+Fϕhn+1dx,and the modified energy of the fully discrete scheme (15) is defined as(85)Ξϕhn+1:=Eϕhn+1+θ2ϕhn+1−ϕ¯0−1,h2+τ2γ∇phn+12.We choose the final timeT=0.1, τ=0.01, and ϵ=0.1, a=0.5, γ=0.01, κ=0.01. The evolution of energy is shown in Figure 19, and it is nonincreasing apparently.Figure 19
The evolution of discrete energy.
(a)(b)In Figure19(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.


## 6. Conclusion

This 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 functionalFϕ extended from −1,1 to −∞,∞ so that the problem of the overflow on the boundary can be avoided. Then, 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.



---

*Source: 7104389-2022-09-26.xml*