 # Numerical approximations and error analysis of the Cahn-Hilliard equation with reaction rate dependent dynamic boundary conditions

We consider numerical approximations and error analysis for the Cahn-Hilliard equation with reaction rate dependent dynamic boundary conditions (P. Knopf et. al., arXiv, 2020). Based on the stabilized linearly implicit approach, a first-order in time, linear and energy stable scheme for solving this model is proposed. The corresponding semi-discretized-in-time error estimates for the scheme are also derived. Numerical experiments, including the comparison with the former work, the convergence results for the relaxation parameter K→0 and K→∞ and the accuracy tests with respect to the time step size, are performed to validate the accuracy of the proposed scheme and the error analysis.

## 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, first introduced in , was originally utilized to describe the phase separation and de-mixing processes of binary mixtures. The standard Cahn-Hilliard equation can be written as follows:

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

where the parameter , () denotes a bounded domain whose boundary

with the unit outer vector field

. The function denotes the difference of two local relative concentrations, in order to describe the binary alloys. The regions with in the domain correspond to the pure phases of the materials, which are separated by a interfacial region whose thickness is proportional to .

In the Cahn-Hilliard equation, denotes 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 denotes the potential in . The classical choice of is the smooth double-well potential

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

which has a double-well structure with two minima at -1 and 1 and a local unstable maximum at 0.

Since the time-evolution of is confined in a bounded domain, suitable boundary conditions are needed. The classical choice is the homogeneous Neumann conditions:

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

where represents the outward normal derivative on . Obviously, the mass conservation law holds in the bulk (i.e., in ) with the no-flux boundary condition (1.4):

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

In addition, the time evolution of the bulk free energy (Eq. (1.2)) is decreasing with the boundary conditions (1.4) and (1.5), namely,

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

When some particular applications (for instance, the hydrodynamic applications such as contact line problems) are taken into consideration, it’s necessary to describe the short-range interactions between the mixture and the solid wall. However, the standard homogeneous Neumann conditions neglects the effects of the boundary to the bulk dynamics. Thus, several dynamic boundary conditions have been proposed and analysed in recent years, see for instance, [22, 27, 10, 12, 5, 6, 21, 17, 19, 18]. These dynamic boundary conditions are based on the system with added surface free energy [7, 8, 16]. The total free energy can be written as

 (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, denotes the thickness of the interfacial region on and the parameter is related to the surface diffusion. When , it corresponds to the moving contact line problem .

In the present work, we summarize three Cahn-Hilliard models with dynamic boundary conditions in detail. All the dynamic boundary conditions of the three models have a Cahn-Hilliard type structure. And they can be interpreted as an -gradient flow of the total free energy.

The first Cahn-Hilliard model with dynamic boundary conditions was proposed by G.R. Goldstein, A. Miranville, and G. Schimperna :

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

In the present work, we denote the model as the GMS model for convenience. Here, denotes the Laplace-Beltrami operator on . Note that this model describes the chemical reactions occurring on the boundary and the chemical potentials in the bulk and on the boundary are the same. Moreover, the dynamic boundary conditions ensure the conservation of the total mass (namely, the sum of the bulk and boundary mass):

 (1.11) ∫Ωϕ(t)dx+∫Γψ(t)dS=∫Ωϕ(0)dx+∫Γψ(0)dS,for all t∈[0,T],

and the energy dissipation law:

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

The second Cahn-Hilliard model with dynamic boundary conditions was proposed by C. Liu and H. Wu :

 (1.13) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ϕ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).

We denote it as the Liu-Wu model for short. Here, denotes the chemical potential on the boundary. The model assumes that there is no mass exchange between the bulk and the boundary, namely, . Different from the GMS model (), the chemical potential and are not directly coupled. Similarly, we can obtain the following mass conservation law:

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

indicating that the Liu-Wu model satisfies the mass conservation law in the bulk and on the boundary respectively. Moreover, the energy dissipation law (1.12) also holds for the Liu-Wu model. The readers can find the well-posedness results for the Liu-Wu model and the GMS model in  and  respectively.

Recently, Knopf 

proposed a new model, which can be interpreted as an interpolation between the Liu-Wu model and the GMS model. It reads as follows,

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

In the present work, we use the authors’ initials and refer it to be the KLLM model for convenience. Here, and represent the chemical potentials in and on , respectively. Notice that and are coupled by the Robin type boundary condition , where the positive parameter is the relaxation parameter. The constant can be interpreted as the reaction rate, since the reactions between the materials are described by . The well-posedness of the system (1.15) and convergence to the Liu-Wu model (as ) and the GMS model (as ) in both the weak and the strong sense have been investigated by Knopf .

The numerical approximations of the Cahn-Hilliard equation and its variants have already been well investigated. There exists extensive efficient techniques for the time discretization, such as the stabilized linearly implicit approach , the convex splitting approach [23, 13], the invariant energy quadratization (IEQ) method [28, 29] and the scalar auxiliary variable (SAV) method. Recently, there have been numerical approximations for the Cahn-Hilliard equation with dynamic boundary conditions, see for instance, [3, 4, 15, 9, 26]. Specifically, for the Liu-Wu model, the finite element scheme has been proposed in [26, 11], where the straightforward discretization based on piecewise linear finite element functions was utilized to simulate the model, and the corresponding nonlinear system was solved by Newton’s method. A recent contribution on the numerical analysis can be found in . For the KLLM model, we refer the readers to  for the finite element numerical approximations and numerical analysis. However, the backward implicit Euler method was used for time discretization in the above finite element schemes, where one needs to solve nonlinear systems at each time step. Recently, based on the stabilized linearly implicit approach, a linear and energy stable numerical scheme has been proposed for the Liu-Wu model  and the corresponding semi-discrete-in-time error estimates are carried out.

Inspired by the numerical scheme in , a first-order in time, linear and energy stable scheme for solving the KLLM model is proposed in the present work. The scheme allows us to simulate the KLLM model as well as the two limit models – the Liu-Wu model and the GMS model. Note that the scheme is highly efficient since one only needs to solve a linear equation at each time step. Numerical simulations are performed in the two-dimensional space to validate the accuracy and stability of the scheme by comparing with the former work. We also investigate the error estimates in semi-discrete-in-time for the scheme. To the best of the authors’ knowledge, the proposed scheme is the first linear numerical scheme to solve the KLLM model and it is the first work to give the corresponding semi-discrete-in-time error estimates.

The rest of the paper is organized as follows. We first present some notions and notation appearing in this article in Section 2. In Section 3, the stabilized scheme for the KLLM model and the energy stability are derived. The error estimates are constructed in Section 4. In Section 5, we present the numerical examples and illustrate the convergence results for and . The accuracy tests are also displayed in this section. Finally, the conclusion is 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 parameters and solutions, and use to say that there is a generic constant such that .

## 3 The Cahn-Hilliard equation with reaction rate dependent dynamic boundary conditions and its numerical scheme

In this section, we first summarize the mass conservation and the energy dissipation law of the KLLM model. Then we propose the stabilized linear numerical scheme and prove the discrete energy dissipation law.

Since is the phase-field order parameter in the bulk, denote its trace as the order parameter on the boundary. In the bulk , assume that is a locally conserved quantity that satisfies the continuity equation

 (3.1) ϕt+∇⋅(ϕu)=0,(x,t)∈Ω×(0,T),

where is the microscopic effect velocity.

We assume that there exists mass exchange between the bulk and the boundary , which is denoted by the flux . Assume that the mass flux is directly driven by differences between the chemical potentials in the sense that

 (3.2) K(J⋅n)=K(ϕu⋅n)=μ−μΓ,(x,t)∈Γ×(0,T),

where is a positive parameter describing the extent of mass exchange. Eq. (3.2) is the boundary condition of .

Assume that the boundary dynamics is characterized by a local mass conservation law analogous to (3.1), such that

 (3.3) ψt+∇Γ⋅(ψv)−J⋅n=0,(x,t)∈Γ×(0,T),

where denotes the microscopic effective tangential velocity field on the boundary . Assume that is a closed manifold, thus, there is no need to impose any boundary condition on .

The mass is conserved in the sense that

 (3.4) ∫Ωϕ(t)dx+∫Γψ(t)dS=∫Ωϕ(0)dx+∫Γψ(0)dS,∀t∈[0,T].

To this end, integrating (3.1) over , we have

 (3.5) ddt∫Ωϕ(⋅,t)dx+∫Γϕu⋅ndS=0,∀t∈(0,T),

and integrating (3.3) over , we have

 (3.6) ddt∫Γψ(⋅,t)dS−∫ΓJ⋅ndS=0,∀t∈(0,T).

Combining (3.5) with (3.6) and the flux , we obtain the total mass conservation law, see (3.4).

Then we show the energy law of the KLLM model, where the total free energy (sum of the bulk and surface free energies) is decreasing in time. Precisely, multiplying the first equation of (1.15) by and integrating over , we get

 (ϕt,μ)Ω=(Δμ,μ)Ω=(∂nμ,μ)Γ−∥∇μ∥2L2(Ω).

Since

 (ϕt,μ)Ω=(ϕt,−εΔϕ+1εF′(ϕ))Ω,
 (ϕt,−εΔϕ)Ω=−(ε∂nϕ,ϕt)Γ+ε2ddt(∫Ω|∇ϕ|2dx),
 (ϕt,1εF′(ϕ))Ω=ddt(∫Ω1εF(ϕ)dx),

we arrive that

 (3.7) ddt(∫Ω1εF(ϕ)dx+ε2∫Ω|∇ϕ|2dx)−(ε∂nϕ,ϕt)Γ=(∂nμ,μ)Γ−∥∇μ∥2L2(Ω).

Multiplying the boundary equation in (1.15) by and integrating over , we get

 (ψt,μΓ)Γ=(ΔΓμΓ,μΓ)Γ−(∂nμ,μΓ)Γ=−∥∇ΓμΓ∥2L2(Γ)−(∂nμ,μΓ)Γ.

Since

 (ψt,μΓ)Γ=(ψt,−δκΔΓψ+1δG′(ψ)+ε∂nϕ)Γ,
 (ψt,−δκΔΓψ)Γ=δκ2ddt(∫Γ|∇Γψ|2dS),
 (ψt,1δG′(ψ))Γ=ddt(∫Γ1δG(ψ)dS),

we arrive that

 (3.8) ddt(∫Γ1δG(ψ)dS+δκ2∫Γ|∇Γψ|2dS)+(ε∂nϕ,ψt)Γ=−(∂nμ,μΓ)Γ−∥∇ΓμΓ∥2L2(Γ).

Adding (3.7) and (3.8) together, we get

 (3.9) ddt(∫Ω1εF(ϕ)+ε2|∇ϕ|2dx+∫Γ1δG(ψ)+δκ2|∇Γψ|2dS) =−∥∇μ∥2L2(Ω)−∥∇ΓμΓ∥2L2(Γ)+(∂nμ,μ−μΓ)Γ =−∥∇μ∥2L2(Ω)−∥∇ΓμΓ∥2L2(Γ)−K∥∂nμ∥2L2(Γ).

Since , we arrive at

 ddt(∫Ω1εF(ϕ)+ε2|∇ϕ|2dx+∫Γ1δG(ψ)+δκ2|∇Γψ|2dS)≤0,

namely,

 ddt[Ebulk(ϕ)+Esurf(ψ)]≤0.

Now we present the numerical scheme for the KLLM model (namely, Eq. (1.15)). The scheme can be written as follows,

 (3.10) ϕn+1−ϕnτ=Δμn+1,in Ω, (3.11) μn+1=−εΔϕn+1+1εF′(ϕn)+s1(ϕn+1−ϕn),in Ω, (3.12) K∂nμn+1=μn+1Γ−μn+1,on Γ, (3.13) ϕn+1|Γ=ψn+1,on Γ, (3.14) ψn+1−ψnτ=ΔΓμn+1Γ−∂nμn+1,on Γ, (3.15) μ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.

###### Remark 3.1.

The parameters . And the stabilization terms and are added in the bulk and on the boundary to enhance the stability, respectively.

###### Remark 3.2.

For the Liu-Wu model, we need to modify Eq. (3.12) to be

 ∂nμn+1=0,on Γ,

and the last term in (3.14) vanishes. For the GMS model, Eq. (3.12) is modified to be

 μn+1|Γ=μn+1Γ,on Γ.

We have the energy stability as follows.

###### Theorem 3.3.

Assume that , , the scheme (3.10)-(3.15) is energy stable in the sense that

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

where

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

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

 (3.18) (ϕn+1−ϕnτ,μn+1)Ω=(Δμn+1,μn+1)Ω=(∂nμn+1,μn+1)Γ−∥∇μn+1∥2L2(Ω).

For the boundary integral term, by using (3.12), we have

 (∂nμn+1,μn+1)Γ=1K(μn+1Ω−μn+1,μn+1)Γ.

By using (3.11), we have

 (3.19) (ϕn+1−ϕnτ,μn+1)Ω=(ϕn+1−ϕnτ,−εΔϕn+1+1εF′(ϕn)+s1(ϕn+1−ϕn))Ω,

and

 (3.20) (ϕn+1−ϕnτ,−εΔϕn+1)Ω=−ε(∂nϕn+1,ϕn+1−ϕnτ)Γ+ε(∇ϕn+1,∇ϕn+1−∇ϕnτ)Ω.

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

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

By using (3.15), we have

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

and

 (3.23) (ψn+1−ψnτ,−δκΔΓψn+1)Γ=(∇Γψn+1−∇Γψnτ,δκ∇Γψn+1)Γ.

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

 (3.24) 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μn+1,μn+1)Γ−∥∇μn+1∥2L2(Ω)−∥∇Γμn+1Γ∥2L2(Γ)−(∂nμn+1,μn+1Γ)Γ =−∥∇μn+1∥2L2(Ω)−∥∇Γμn+1Γ∥2L2(Γ)−1K∥μn+1−μn+1Γ∥2L2(Γ),

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∥2L2(Ω) +(δκ∇Γψn+1,∇Γψn+1−∇Γψnτ)Γ+1δ(G′(ψn),ψn+1−ψnτ)Γ+s2τ∥ψn+1−ψn∥2L2(Γ) =ε(∇ϕn+1,∇ϕn+1−∇ϕnτ)Ω+1ε(F(ϕn+1)−F(ϕn)τ,1)Ω−12ε(F′′(η),(ϕn+1−ϕn)2τ)Ω +s1τ∥ϕn+1−ϕn∥2L2(Ω)+δκ(∇Γψn+1,∇Γψn+1−∇Γψnτ)Γ+1δ(G(ψn+1)−G(ψn)τ,1)Γ −12δ(G′′(ζ),(ψn+1−ψn)2τ)Γ+s2τ∥ψn+1−ψn∥2L2(Γ) =ε2τ(∥∇ϕn+1∥2L2(Ω)−∥∇ϕn∥2L2(Ω)+∥∇ϕn+1−∇ϕn∥2L2(Ω))+1ετ(F(ϕn+1)−F(ϕn),1)Ω +1τ(s1−12εF′′(η))∥ϕn+1−ϕn∥2L2(Ω)+δκ2τ(∥∇Γψn+1∥2L2(Γ)−∥∇Γψn∥2L2(Γ)+∥∇Γψn+1−∇Γψn∥2L2(Γ)) +1δτ(G(ψn+1)−G(ψn),1)Γ+1τ(s2−12δG′′(ζ))∥ψn+1−ψn∥2L2(Γ) =1τ[E(ϕn+1,ψn+1)−E(ϕn,ψn)]+ε2τ∥∇ϕn+1−∇ϕn∥2L2(Ω)+δκ2τ∥∇Γψn+1−∇Γψn∥2L2(Γ) +1τ(s1−12εF′′(η))∥ϕn+1−ϕn∥2L2(Ω)+1τ(s2−12δG′′(ζ))∥ψn+1−ψn∥2L2(Γ).

Thus, we have

 1τ[E(ϕn+1,ψn+1)−E(ϕn,ψn)]+ε2τ∥∇ϕn+1−∇ϕn∥2L2(Ω)+δκ2τ∥∇Γψn+1−∇Γψn∥2L2(Γ) +1τ(s1−12εF′′(η))∥ϕn+1−ϕn∥2L2(Ω)+1τ(s2−12δG′′(ζ))∥ψn+1−ψn∥2L2(Γ) =−∥∇μn+1∥2L2(Ω)−∥∇Γμn+1Γ∥2L2(Γ)−1K∥μn+1−μn+1Γ∥2L2(Γ)≤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.10)-(3.15) 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 numerical scheme (3.10)-(3.15). During the estimate, the mathematics induction is utilized and the trace theorem is applied to estimate the boundary terms.

Assume that the Lipschitz properties hold for and , and and are bounded:

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

which are necessary for error estimates.

###### Remark 4.1.

One example of the functionals and , satisfying the assumptions mentioned above, is the modified double-well potential:

 (4.3) F(ϕ)=G(ϕ)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩(ϕ−1)2ϕ>1,14(ϕ2−1)2−1≤ϕ<1,(ϕ+1)2ϕ≤−1.

The Lipschitz property holds for and and

 (4.4) maxϕ∈R|F′′(ϕ)|=maxψ∈R|G′′(ψ)|≤2.

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

 (4.5) ϕ(tn+1)−ϕ(tn)τ=Δμ(tn+1)+Rn+1ϕ,in