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:
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:
where denotes the potential in . The classical choice of is the smooth double-well potential
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:
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):
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
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 :
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):
and the energy dissipation law:
The second Cahn-Hilliard model with dynamic boundary conditions was proposed by C. Liu and H. Wu :
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:
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,
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.
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 :
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
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
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
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
To this end, integrating (3.1) over , we have
and integrating (3.3) over , we have
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
we arrive that
Multiplying the boundary equation in (1.15) by and integrating over , we get
we arrive that
Since , we arrive at
Now we present the numerical scheme for the KLLM model (namely, Eq. (1.15)). The scheme can be written as follows,
Here, is an arbitrary and fixed time, is the number of time steps and is the step size.
The parameters . And the stabilization terms and are added in the bulk and on the boundary to enhance the stability, respectively.
We have the energy stability as follows.
By taking inner product of (3.10) with in , we have
For the boundary integral term, by using (3.12), we have
By using (3.11), we have
By using (3.15), we have
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:
which are necessary for error estimates.
One example of the functionals and , satisfying the assumptions mentioned above, is the modified double-well potential:
The Lipschitz property holds for and and
The PDE system (1.15) can be rewritten as the following truncated form,