 # Numerical Approximations and Error Analysis of the Cahn-Hilliard Equation with Dynamic Boundary Conditions

We consider the numerical approximations of the Cahn-Hilliard equation with dynamic boundary conditions (C. Liu et. al., Arch. Rational Mech. Anal., 2019). We propose a first-order in time, linear and energy stable numerical scheme, which is based on the stabilized linearly implicit approach. The energy stability of the scheme is proved and the semi-discrete-in-time error estimates are carried out. Numerical experiments, including the comparison with the former work, the accuracy tests with respect to the time step size and the shape deformation of a droplet, are performed to validate the accuracy and the stability of the proposed scheme.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The Cahn-Hilliard equation is one of the most fundamental models which describe the phase separation processes of binary mixtures. The classical Cahn-Hilliard equation, first introduced in , can be written as follows:

 (1.1) ⎧⎨⎩ϕt=Δμ,in%  Ω×(0,T),μ=−εΔϕ+1εF′(ϕ),in Ω×(0,T),

where and () is a bounded domain with the smooth boundary . The phase-field order parameter represents the difference of two local relative concentrations to describe the binary mixtures. In the domain , correspond to the pure phases of the materials, which are separated by a interfacial region whose thickness is proportional to the parameter . represents the chemical potential in , which can be expressed as the Fréchet derivative of the bulk free energy:

 (1.2) Ebulk(ϕ)=∫Ωε2|∇ϕ|2+1εF(ϕ)dx,

where stands for the bulk potential, which usually has a double-well structure with two minima at -1 and 1 and a local unstable maximum at 0. A classical choice is the regular double-well potential

 (1.3) F(x)=14(x2−1)2,x∈R.

When the time-evolution of is confined in a bounded domain, suitable boundary conditions should be considered for the system (1.1). The homogeneous Neumann conditions are the classical boundary conditions:

 (1.4) ∂nμ=0,on Γ×(0,T),
 (1.5) ∂nϕ=0,on Γ×(0,T),

where

denotes the unit outer normal vector and

denotes the outward normal derivative on . The Cahn-Hilliard equation with the boundary conditions (1.4) and (1.5) can be viewed as an -gradient flow of the bulk free energy.

The no-flux boundary condition (1.4) guarantees the conservation of mass in the bulk (i.e., in ):

 (1.6) ∫Ωϕ(t)dx=∫Ωϕ(0)dx,t∈[0,T].

Moreover, the boundary conditions (1.4) and (1.5) imply that the bulk free energy (Eq. (1.2)) is decreasing with respect to time, namely,

 (1.7) ddtEbulk(ϕ(t))+∫Ω|∇μ|2dx=0,t∈(0,T).

However, the Cahn-Hilliard equation with homogeneous Neumann conditions neglects the effects of the boundary to the bulk dynamics. Thus, it is not suitable for some applications (for instance, hydrodynamic applications such as contact line problems). In order to describe the effective interactions between the solid wall and the binary mixture, physicists added the suitable surface free energy functional into the system [4, 5, 14]:

 (1.8) Etotal(ϕ)=Ebulk(ϕ)+Esurf(ϕ),
 (1.9) Esurf(ϕ)=∫Γδκ2|∇Γϕ|2+1δG(ϕ)dS,

where represents the tangential or surface gradient operator on , is the surface potential, the parameter is related to the surface diffusion and denotes the thickness of the interfacial region on . When , it corresponds to the moving contact line problem . Recently, for the total free energy (1.8), various dynamic boundary conditions for the Cahn-Hilliard equation have been proposed and investigated, see for instance, see [14, 9, 17, 15, 16], and references therein.

In the present work, the Cahn-Hilliard equation with the dynamic boundary conditions, which was derived by an energetic variational approach by Liu and Wu (Liu-Wu model, for short) , is considered. It reads as follows:

 (1.10) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ϕt=Δμ,in%  Ω×(0,T),μ=−εΔϕ+1εF′(ϕ),in Ω×(0,T),∂nμ=0,on Γ×(0,T),ϕ|Γ=ψ,on Γ×(0,T),ψt=ΔΓμΓ,on Γ×(0,T),μΓ=−δκΔΓψ+1δG′(ψ)+ε∂nϕon Γ×(0,T),

where denotes the Laplace-Beltrami operator on . The dynamic boundary conditions (with , ) turns out to be a surface Cahn-Hilliard type equation for the trace of on , coupled with the bulk evolution in terms of . The existence and uniqueness of weak and strong solutions of the Liu-Wu model have been established in . A different approach to construct the weak solutions of the Liu-Wu model is proposed in .

The numerical approximations of the Cahn-Hilliard equation and its variants have been intensively investigated. The stabilized linearly implicit approach , the approaches based on the convex-concave splitting [19, 10], the invariant energy quadratization [23, 24] and the scalar auxiliary variable (SAV)  method are efficient techniques for the time discretization. Recently, there have been numerous contributions on the numerical approximation of the Cahn-Hilliard equation with dynamic boundary conditions [2, 3, 13, 7, 22]. For the numerical approximations of Liu-Wu model, the first finite element scheme was proposed in  and the corresponding numerical results were presented in , where the straightforward discretization based on piecewise linear finite element functions was utilized to simulate Liu-Wu model, and the corresponding nonlinear system was solved by Newton’s method. A recent contribution on the numerical analysis for the Liu-Wu model can be found in , where a different discrete scheme was proposed and the connection between and the chemical potentials was investigated. However, the backward implicit Euler method was used for time discretization in the above discrete schemes, which lead to nonlinear systems at each time step.

In the present work, a first-order in time, linear and energy stable scheme for solving the Liu-Wu model is proposed based on the stabilized linearly implicit approach. At each time step, one only needs to solve one linear equation and thus, the scheme is highly efficient. The energy stability of the scheme is proved and various numerical simulations in two-dimensional spaces are performed to validate the accuracy and stability of the scheme by comparing with the former work. The error estimates in semi-discrete-in-time for the scheme are also carried out. To the best of the authors’ knowledge, the proposed scheme in this paper is the first linear and energy stable scheme for solving the Liu-Wu model and it is the first work to give the semi-discrete-in-time error estimates for the model.

The rest of the paper is organized as follows. In Section 2, we first recall some notions and notation appearing in this article. In Section 3, a simple derivation of Liu-Wu model and the stabilized scheme with the energy stability are derived. In Section 4, we construct the error estimates. The accuracy tests and numerical examples are presented in Section 5. Finally, some concluding remarks are presented in Section 6.

## 2 Preliminaries

Before giving the stabilized scheme and the corresponding error analysis, we make some definitions in this section. The norm and inner product of and are denoted by , and , respectively. The usual norm in and are denoted by and respectively.

We consider a finite time interval and a domain (), which is a bounded domain with sufficient smooth boundary and is the unit outer normal vector on .

Let be the time step size. For a sequence of functions in some Hilbert space , we denote the sequence by and define the following discrete norm for :

 (2.1) ∥fτ∥l∞(E)=max0≤n≤N(∥fn∥E).

We denote by a generic constant that is independent of but possibly depends on the data and the solution, and use to say that there is a generic constant such that .

## 3 Derivation of the Cahn-Hilliard equation with dynamic boundary conditions and its numerical scheme

In this section, we first propose a simple derivation of the Liu-Wu model, indicating that it satisfies the energy dissipation law and mass conservation.

Since is the phase-field order parameter in the bulk, denote its trace as the order parameter on the boundary. When the mass conservation holds true in the bulk and on the boundary respectively, and satisfy the following continuity equations :

 (3.1) ϕt+∇⋅(ϕu)=0,in Ω×(0,T), ψt+∇Γ⋅(ϕv)=0,on Γ×(0,T),

where is the microscopic effect velocity and is the microscopic effective tangential velocity field on the boundary. Assume that there is no mass exchange between the bulk and the boundary, satisfies the following boundary condition:

 (3.2) u⋅n=0,on Γ×(0,T).

Since the boundary is closed, there is no need to impose any boundary condition on .

For the total free energy (1.8), we consider the following energy dissipation law:

 (3.3) ddtEtotal(t)=−Dtotal(t),

which is based on the first and second laws of thermodynamics . And the rate of energy dissipation also consists of two contributions from the bulk and the boundary, namely,

 (3.4) Dtotal(t)=Dbulk(t)+Dsurf(t).

Assume that

 (3.5) Dbulk(t)=∫Ωϕ2|u|2dx,Dsurf(t)=∫Γψ2|v|2dS,

and substitute (1.2) and (1.9) into (3.3), we obtain

 (3.6) ddt[∫Ωε2|∇ϕ|2+1εF(ϕ)dx+∫Γδκ2|∇Γϕ|2+1δG(ϕ)dS] =−∫Ωϕ2|u|2dx−∫Γψ2|v|2dS.

The left part of (3.6) can be written as

 (3.7) ddt[∫Ωε2|∇ϕ|2+1εF(ϕ)dx+∫Γδκ2|∇Γϕ|2+1δG(ϕ)dS] =∫Ω[−εΔϕ+1εF′(ϕ)]ϕtdx+∫Γ[ε∂nϕ−δκΔΓψ+1δG′(ψ)]ψtdS =∫Ωϕ∇[−εΔϕ+1εF′(ϕ)]⋅udx+∫Γψ∇Γ[ε∂nϕ−δκΔΓψ+1δG′(ψ)]⋅vdS =∫Ωϕ∇μ⋅udx+∫Γψ∇Γ[ε∂nϕ+μs]⋅vdS.

Here, and are the chemical potentials in and on , respectively, which can be expressed as the Fréchet derivative of the bulk free energy and the surface free energy , namely,

 (3.8) μ=δEbulkδϕ=−εΔϕ+1εF′(ϕ) μs=δEsurfδψ=−δκΔΓψ+1δG′(ψ).

From (3.6) and (3.7), we obtain

 (3.9) ϕu=−∇μ, ψv=−∇Γ(ε∂nϕ+μs).

Substitute (3.9) into (3.1), we obtain the Liu-Wu model (Eq. (1.10)) with .

Obviously, the Liu-Wu model satisfies the mass conservation law in the bulk and on the boundary:

 (3.10) ∫Ωϕ(t)dx=∫Ωϕ(0)dxand∫Γψ(t)dS=∫Γψ(0)dS,t∈[0,T],

Moreover, from the energy dissipation law, it is easy to see that the total free energy is decreasing in time:

 (3.11) ddtEtotal(ϕ,ψ)=−∥∇μ∥2Ω−∥∇ΓμΓ∥2Γ≤0.

Now we present the stabilized scheme for the Cahn-Hilliard equation with dynamic boundary conditions (namely, (1.10)). The scheme can be written as follows,

 (3.12) ϕn+1−ϕnτ=Δμn+1,in Ω, (3.13) μn+1=−εΔϕn+1+1εF′(ϕn)+s1(ϕn+1−ϕn),in Ω, (3.14) ∂nμn+1=0,on Γ, (3.15) ϕn+1|Γ=ψn+1,on Γ, (3.16) ψn+1−ψnτ=ΔΓμn+1Γ,on Γ, (3.17) μn+1Γ=−δκΔΓψn+1+1δG′(ψn)+ε∂nϕn+1+s2(ψn+1−ψn),on Γ.

Here, is an arbitrary and fixed time, is the number of time steps and is the step size. We have the energy stability as follows.

###### Theorem 3.1.

Assume that , , the scheme (3.12)-(3.17) is energy stable in the sense that

 (3.18) E(ϕn+1,ψn+1)−E(ϕn,ψn)τ ≤−∥∇μn+1∥2Ω−∥∇Γμn+1Γ∥2Γ,

where

 (3.19) E(ϕn,ψn)=∫Ω1εF(ϕn)+ε2|∇ϕn|2dx+∫Γ1δG(ψn)+δκ2|∇Γψn|2dS
###### Proof.

By taking inner product of (3.12) with in , we have

 (3.20) (ϕn+1−ϕnτ,μn+1)Ω=(Δμn+1,μn+1)Ω=−∥∇μn+1∥2Ω.

By using (3.13), we have

 (3.21) (ϕn+1−ϕnτ,μn+1)Ω =(ϕn+1−ϕnτ,−εΔϕn+1+1εF′(ϕn)+s1(ϕn+1−ϕn))Ω,
 (3.22) (ϕn+1−ϕnτ,−εΔϕn+1)Ω =−ε(∂nϕn+1,ϕn+1−ϕnτ)Γ+ε(∇ϕn+1,∇ϕn+1−∇ϕnτ)Ω.

For the boundary integral term in (3.22), by taking the inner product of (3.16) with on , we obtain

 (3.23) (ψn+1−ψnτ,μn+1Γ)Γ=(ΔΓμn+1Γ,μn+1Γ)Γ=−∥∇Γμn+1Γ∥2Γ.

By using (3.17), we have

 (3.24) (ψn+1−ψnτ,μn+1Γ)Γ =(ψn+1−ψnτ,−δκΔΓψn+1+1δG′(ψn) +ε∂nϕn+1+s2(ψn+1−ψn))Γ,
 (3.25) (ψn+1−ψnτ,−δκΔΓψn+1)Γ=(∇Γψn+1−∇Γψnτ,δκ∇Γψn+1)Γ.

To handle the nonlinear term associated with and in (3.21) and (3.24), we need the following identities

 (3.26) F′(ϕn)(ϕn+1−ϕn) =F(ϕn+1)−F(ϕn)−F′′(η)2(ϕn+1−ϕn)2, G′(ϕn)(ϕn+1−ϕn) =G(ϕn+1)−G(ϕn)−G′′(ζ)2(ϕn+1−ϕn)2.

Combining the equations mentioned above, we get

 (ϕn+1−ϕnτ,μn+1)Ω+(ψn+1−ψnτ,μn+1Γ)Γ=−∥∇μn+1∥2Ω−∥∇Γμn+1Γ∥2Γ,

and

 (ϕn+1−ϕnτ,μn+1)Ω+(ψn+1−ψnτ,μn+1Γ)Γ =ε(∇ϕn+1,∇ϕn+1−∇ϕnτ)Ω+1ε(F′(ϕn),ϕn+1−ϕnτ)Ω+s1τ∥ϕn+1−ϕn∥2Ω +(δκ∇Γψn+1,∇Γψn+1−∇Γψnτ)Γ+1δ(G′(ψn),ψn+1−ψnτ)Γ+s2τ∥ψn+1−ψn∥2Γ =ε(∇ϕn+1,∇ϕn+1−∇ϕnτ)Ω+1ε(F(ϕn+1)−F(ϕn)τ,1)Ω−12ε(F′′(η),(ϕn+1−ϕn)2τ)Ω +s1τ∥ϕn+1−ϕn∥2Ω+δκ(∇Γψn+1,∇Γψn+1−∇Γψnτ)Γ+1δ(G(ψn+1)−G(ψn)τ,1)Γ −12δ(G′′(ζ),(ψn+1−ψn)2τ)Γ+s2τ∥ψn+1−ψn∥2Γ =ε2τ(∥∇ϕn+1∥2Ω−∥∇ϕn∥2Ω+∥∇ϕn+1−∇ϕn∥2Ω)+1ετ(F(ϕn+1)−F(ϕn),1)Ω +1τ(s1−12εF′′(η))∥ϕn+1−ϕn∥2Ω+δκ2τ(∥∇Γψn+1∥2Γ −∥∇Γψn∥2Γ+∥∇Γψn+1−∇Γψn∥2Γ)+1δτ(G(ψn+1)−G(ψn),1)Γ +1τ(s2−12δG′′(ζ))∥ψn+1−ψn∥2Γ =1τ[E(ϕn+1,ψn+1)−E(ϕn,ψn)]+ε2τ∥∇ϕn+1−∇ϕn∥2Ω +δκ2τ∥∇Γψn+1−∇Γψn∥2Γ+1τ(s1−12εF′′(η))∥ϕn+1−ϕn∥2Ω +1τ(s2−12δG′′(ζ))∥ψn+1−ψn∥2Γ.

Thus, we have

 1τ[E(ϕn+1,ψn+1)−E(ϕn,ψn)]+ε2τ∥∇ϕn+1−∇ϕn∥2Ω +δκ2τ∥∇Γψn+1−∇Γψn∥2Γ+1τ(s1−12εF′′(η))∥ϕn+1−ϕn∥2Ω +1τ(s2−12δG′′(ζ))∥ψn+1−ψn∥2Γ=−∥∇μn+1∥2Ω−∥∇Γμn+1Γ∥2Γ≤0.

Therefore, under the conditions that

 s1≥12εmaxξ∈RF′′(ξ)

and

 s2≥12δmaxη∈RG′′(η),

we have

 1τ[E(ϕn+1,ψn+1)−E(ϕn,ψn)]≤0,

namely, the scheme (3.12)-(3.17) is energy stable. ∎

## 4 Error estimates for the stabilized semi-discrete scheme

In this section, we establish the error estimates for the phase functions and for the stabilized scheme (3.12)-(3.17).

Assume that the Lipschitz properties hold for the derivatives of and ,

 (4.1) maxϕ∈R|F′′(ϕ)|≤L1,
 (4.2) maxψ∈R|G′′(ψ)|≤L2,

which are necessary for error estimates.

The PDE system (1.10) can be rewritten as the following truncated form,

 (4.3) ϕ(tn+1)−ϕ(tn)τ=Δμ(tn+1)+Rn+1ϕ,in Ω, (4.4) μ(tn+1)=−εΔϕ(tn+1)+1εF′(ϕ(tn))+s1(ϕ(tn+1)−ϕ(tn))+Rn+1μ,in Ω, (4.5) ∂nμ(tn+1)=0,on Γ, (4.6) ϕ(tn+1)|Γ=ψ(tn+1),on Γ, (4.7) ψ(tn+1)−ψ(tn)τ=ΔΓμΓ(tn+1)+Rn+1ψ,on Γ, μΓ(tn+1)=−δκΔΓψ(tn+1)+1δG′(ψ(tn))+ε∂nϕ(tn+1) (4.8) +s2(ψ(tn+1)−ψ(tn))+Rn+1Γ,on Γ,

where

 (4.9) Rn+1ϕ=ϕ(tn+1)−ϕ(tn)τ−ϕt(tn+1),
 (4.10) Rn+1ψ=ψ(tn+1)−ψ(tn)τ−ψt(tn+1),
 (4.11) Rn+1μ=1εF′(ϕ(tn+1))−1εF′(ϕ(tn))−s1(ϕ(tn+1)−ϕ(tn)),
 (4.12) Rn+1Γ=1δG′(ψ(tn+1))−1δG′(ψ(tn))−s2(ψ(tn+1)−ψ(tn)).

We assume that the exact solutions of the system (1.10) are sufficiently smooth. From the Taylor expansion, it’s easy to prove that

###### Lemma 4.1.

The truncation errors satisfy

 (4.13) ∥Rϕ,τ∥l∞(H1(Ω))+∥Rμ,τ∥l∞(H1(Ω))≲τ, ∥Rϕ,τ∥l∞(L2(Ω))+∥Rμ,τ∥l∞(L2(Ω))≲τ, ∥Rψ,τ∥l∞(H1(Γ))+∥RΓ,τ∥l∞(H1(Γ))≲τ, ∥Rψ,τ∥l∞(L2(Γ))+∥RΓ,τ∥l∞(L2(Γ))≲τ.

Thus we can establish the estimates for the stabilized scheme as follows.

###### Theorem 4.2.

Provided that the exact solutions are sufficiently smooth, the solution of the scheme (3.12)-(3.17) satisfy the following error estimates

 (4.14) ∥eϕ,τ∥l∞(H1(Ω))+∥eψ,τ∥l∞(H1(Γ))≲τ, ∥eϕ,τ∥l∞(L2(Ω))+∥eψ,τ∥l∞(L2(Γ))≲τ.

Here, the error functions are defined as

 (4.15) enϕ=ϕ(tn)−ϕn,enμ=μ(tn)−μn, enψ=ψ(tn)−ψn,enΓ=μΓ(tn)−μnΓ,

, and the corresponding sequence of error functions are denoted as , , and .

###### Proof.

By subtracting (4.3)-(4) from the corresponding scheme (3.12)-(3.17), we derive the error equations as follows,

 (4.16) 1τ(en+1ϕ−enϕ)=Δen+1μ+Rn+1ϕ,in Ω, (4.17) en+1μ=−εΔen+1ϕ+1ε(F′(ϕ(tn))−F′(ϕn))+s1(en+1ϕ−enϕ)+Rn+1μ,in Ω, (4.18) ∂nen+1μ=0,on Γ, (4.19) en+1ϕ|Γ=en+1ψ,on Γ, (4.20) 1τ(en+1ψ−enψ)=ΔΓen+1Γ+Rn+1ψ,on Γ, en+1Γ=−δκΔΓen+1ψ+1δ(G′(ψ(tn))−G′(ψn))+ε∂nen+1ϕ (4.21) +s2(en+1ψ−enψ)+Rn+1Γ,on Γ.

By taking the inner product of (4.16) with in , we obtain

 (en+1ϕ−enϕ,en+1μ)Ω+τ∥∇en+1μ∥2Ω=τ(Rn+1ϕ,en+1μ)Ω.

By taking the inner product of (4.16) with in , we obtain

 ε2(∥en+1ϕ∥2Ω−∥enϕ∥2Ω+∥en+1ϕ−enϕ∥2Ω) =−ετ(∇en+1μ,∇en+1ϕ)Ω+ετ(Rn+1ϕ,en+1ϕ)Ω,

where the boundary terms vanish due to . By taking the inner product of (4.17) with in , we obtain

 −(en+1μ,en+1ϕ−enϕ)Ω+ε2(∥∇en+1ϕ∥2Ω−∥∇e