The error analysis for the Cahn-Hilliard phase field model of two-phase incompressible flows with variable density

: In this paper, we consider the numerical approximations of the Cahn-Hilliard phase field model for two-phase incompressible flows with variable density. First, a temporal semi-discrete numerical scheme is proposed by combining the fractional step method (for the momentum equation) and the convex-splitting method (for the free energy). Second, we prove that the scheme is unconditionally stable in energy. Then, the L 2 convergence rates for all variables are demonstrated through a series of rigorous error estimations. Finally, by applying the finite element method for spatial discretization, some numerical simulations were performed to demonstrate the convergence rates and energy dissipations.

This model is also called the Cahn-Hilliard-Navier-Stokes model when the density is constant, which has many practical applications in physical and engineering, such as wetting, coating, and painting.By considering the influence of variable density, this model has broad applicability, which include highly stratified flows, interfaces between fluids of different densities and some problems of inertial confinement.
The phase field method which has a wide range of applications is one of the main methods to deal with the fluid interface in two-phase flow modeling; see [2][3][4][5] and the references therein.It was initially developed to simulate solid-liquid phase transitions, where the interface is treated as a thin, smooth transition layer [6][7][8] to remove the singularities between two phases.The basic framework is to use phase field variables to represent the volume fraction of fluid components and then adopt a variational form to derive the model.Recently, it has increasingly attracted researchers' interest, mainly because the phase field method is superior to other available methods in some aspects of twophase flow [9].Various material properties or complex interface behaviors can be simulated directly by introducing suitable energy functions.Numerical solutions play a crucial role in their study and applications because the analytic solutions are usually not available.
A few works have been devoted to the design, analysis and implementation of numerical schemes for the phase field model, although this model is very classical and canonical.We briefly review the available methods used in the phase field model.It is worth noting that there are projection/gauge/penalty methods [10][11][12], scalar auxiliary variables [13,14], linear stability [15,16], convex splitting [17][18][19], invariable energy quadratization(IEQ) [20,21], nonlinear quadratic [22], exponential time differencing [23], etc.In practical applications, we often couple the flow-field equation with the phase field equation.Typically two-phase incompressible flow models are coupled to phase-field models.The Cahn-Hilliard model is taken into account in this paper, because it is effective in the following two aspects: (i) Cahn-Hilliard models can accurately conserve the volume and dynamics; (ii) the equation is one of the most important models in mathematical physics.Because of these reasons, here we use the Cahn-Hilliard equation developed in [24] to couple the incompressible flows with variable density.
There also have been a lot of works on numerical approximates for the Cahn-Hilliard phase field model for incompressible two-phase flows with variable density.Hohenberg and Halperin proposed the model in [25] to simulate two incompressible viscous fluids with constant density.In [26], Gurtin et al. obtained the equal model by using the framework of rational continuum mechanics.A fully adaptive energy stabilization scheme is proposed in [27].An efficient Picard iteration procedure was designed in [28] to further decouple the model.In the last years, many authors have been concerned with designing incompressible two-phase flow models with variable densities.Several efficient and energy-stable time discretization schemes for the coupled nonlinear Cahn-Hilliard phase field system with variable density are constructed; see [29].Yang and Dong presented an energy-stable scheme in [30] for the numerical approximation of the two-phase governing equations with variable density and viscosity for the two fluids by introducing a scalar-valued variable related to the total of the kinetic energy and the potential free energy.A second-order accurate, coupled, energy-stable schemeis proposed in [31], where the Crank-Nicolson method and the IEQ method were used.In [32], the conservation scheme of the first-order energy law was established, in which the Cahn-Hilliard solver was used to decouple from the two-phase incompressible flows solver through the use of the fractional step method.Ye et al. [33] have designed a fully-decoupled type scheme to solve the Cahn-Hilliard phase field model for a two-phase incompressible fluid flow system with constant density.They only give a detailed practical implementation method and also prove the solvability.None of the various existing schemes have been subjected to error analysis, where the main difficulty lies in the delicate treatment of a several of nonlinear terms.Rigorous error estimates of models with variable densities, using an optimal order error bound, may seem to be a difficult prospect, but is a very interesting direction for future research.
In this paper, we finally arrive at an unconditionally stable in energy, first-order time-accurate scheme for the incompressible Cahn-Hilliard two phase flows with variable density by coming up with a fractional step method.This method has an advantage over the projection method in that the original boundary conditions of the problem can be implemented in all substeps of the scheme.The popular approach to discretizing the Cahn-Hilliard phase field model (1.1b) and (1.1c) in time is based on the convex-splitting of the free energy functional, i.e., an idea that can be traced back to [32].In the convex-splitting framework, one treats the contribution from the convex part implicitly and the contribution from the concave part explicitly.This treatment promotes the energy stability of the scheme and this property is unconditional in terms of time steps.We also give a rigorous proof of the convergence results and error estimates in the theoretical analysis.The main contribution of this paper is a rigorous error analysis, particularly under the condition that energy stability is available.To the best of the authors' knowledge, the proof developed in this article is the first to have the description of error estimation.The accuracy and stability are also demonstrated through the simulation of various numerical examples, where the challenge is in creating an efficient and easy to implement numerical scheme that preserves the energy dissipation law.
The rest of this article is organized as follows.In Section 2, we construct an efficient time discrete scheme for variable density and derive unconditional energy stability.In Section 3, the error analysis of the semi-discrete scheme in time is provided.Some numerical experimentations are given in Section 4. Finally, conclusions are drawn in Section 5.

The discrete scheme
For the sake of simplicity, some notations are needed for the following content.We assume that the domain Ω ∈ R 2 is open, sufficiently smooth and bounded.For any two functions ϕ(x) and ψ(x), their L 2 inner product on Ω is denoted by (ϕ, ψ) = Ω ϕ(x)ψ(x)dx, and the L 2 norm of ϕ(x) is denoted by ∥ϕ∥ = (ϕ, ϕ) 1 2 .Let τ > 0 be the time step size and set t n = nτ for 0 ≤ n ≤ N with T = Nτ.Moreover, we introduce the following spaces, where n is the outward normal vector of ∂Ω.Next, we reformulate the system (1.1) as follows: ) ) ) Remark 1.In (2.1d), the term 1 4 ρ n+1 (∇ • u n ) ũn+1 is added to obtain the unconditional stability, as it is 0 if ∇ • u n = 0. Theorem 1. (S tability o f ρ) For any τ > 0 and any sequence {u n } n=0,...,N satisfying the boundary condition u n • n = 0 on ∂Ω, the solution {ρ n } n=1,...,N to (2.1a) satisfies (2.2) Proof.Testing (2.1a) with 2τρ n+1 gives Owing to the boundary conditions on u n , we note that (2.4) Thus, we get Summing all indices n ranging 0 to N − 1, the proof is completed.□ Next, the stability of (2.1b)-(2.1e)will be proved in the theorem below.Moreover, since the kinetic energy of the fluid is 1

√
ρ n u n 2 , it is more suitable to establish a bound based on √ ρ n u n 2 than on the velocity itself.For simplicity, let us say that σ n = √ ρ n for all 1 ≤ n ≤ N and σ 0 = √ ρ 0 .
Theorem 2. (S tability o f energy) For any τ > 0, (2.1b)-(2.1e)satisfy the the conditions of following energy estimates: Proof.Multiplying (2.1b) by 2τλw n+1 and integrating over Ω, we get where we use the following identity: By taking the inner product of (2.1d) with 2τ ũn+1 , we have (2.9) (2.10) Taking into account the boundary condition on u n and using integration by parts, we have (2.11) Summing the above inequality, we arrive at (2.12) Adding up the above inequality from n = 0 to N − 1, we obtain Theorem 2.

□
The bound on the pressure p is proved in the following theorem.
Theorem 3. (S tability o f p) For any τ > 0, the solution p n+1 to (2.1e) satisfies the following inequality: Proof.Under the inf-sup condition, there exists a positive constant β such that (2.15) Given that ∥ρ n ∥ ≤ ∥ρ 0 ∥ for all 1 ≤ n ≤ N, and by using Hölder's inequality, we have Then, by the Sobolev embedding inequality ∥v∥ L 6 ≤ C∥∇v∥ L 2 for any v ∈ V, we get (2.17) Substituting the above inequalities into Eq (2.14), we obtain (2.18) And by using Theorem 2.2, we get the desired result.□

Temporal error estimates
In this section, we will give the time error estimates and show that the scheme has a first-order convergence rate.Although we verified that the scheme (2.1) is unconditionally stable in the previous chapter, we need to make the following assumptions [34,35] when conducting temporal error analysis: where χ is a number in (0, ρ min 0 ].We assume that the exact solution (ρ, u, ϕ, w, p) is sufficiently smooth.To be more precise, We denote In (1.1), taking t = t n+1 and subtracting from (2.1), we get the following error equations: where If the exact solution is sufficiently smooth, it is easy to establish the following estimate of the truncation error.
Lemma 1.Under the regularity assumptions given by (3.3), the truncation errors satisfy: Proof.By using the integral residual of the Taylor formula, we have By Hölder's inequality, we can derive Similarly, we can prove the inequality of R n+1 ϕ .For R n+1 u ,we can rewrite Using R n+1 ρ estimation and Hölder's inequality can yield the result for R n+1 u .□ We introduce the following Gronwall's inequality, which will frequently be used in error estimates.
Lemma 2. Let a k , b k , c k and γ k , for integers k ≥ 0, be the nonnegative numbers such that Suppose that τγ k < 1, for all k, and set σ k = (1 − τγ k ) −1 .Then, We verify Lemma 5 by the following lemma.
Then, for n< T τ − 1, we have (3.9) Proof.For G n+1 c , we use (3.8) to conclude that where we have used the a priori bound ∥ϕ n ∥ H 1 ≤ C implied by the stability result given by Theorem 2.
Using the H 2 regularity results for elliptic equations, we conclude that from (3.4c), we know that thus,

□
The error estimate for the discrete density ρ n+1 is derived in the following lemma. (3.12) Using ∇ • u (t n+1 ) = 0 in Ω and u (t n+1 ) = 0 on ∂Ω, we have By using Young's inequality, the Cauchy-Schwarz inequality and Lemma 1, we have where ε > 0 is a sufficiently small constant.Then, by the Sobolev inequality and Young inequality, we get (3.14) Similarly, the last term can be estimated as follows: If ε is sufficiently small such that ετ ≤ 1 6 , substituting the estimates of K i (1 ≤ i ≤ 4) into (3.12),we have e n+1 Thanks to ∇ • u n = 0 in Ω, we have Therefore, we get where the Sobolev embedding inequality ∥v∥ L 6 ≤ C∥∇v∥ L 2 for any v ∈ V is used.
For A 4 , we have and (3.25) From Lemma 1, we have From Lemma 3, we obtain where we use this identity We denote Taking the gradient of Gn+1 , we get In view of (3.10) and the bound ∥ϕ n ∥ H 1 < C, we have   (3.37) Substituting the above estimates into (3.18),we have ( Adding up from 0 to N − 1, and applying Gronwall's inequality, we infer that

□
Lemma 5 shows that the fractional scheme converges at a rate of O(τ ).However, since we use the first-order backward Euler method of time discretization, it is not optimal from the perspective of theoretical analysis.Next, we will improve the convergence speed to the first order.Lemma 6. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)-(3.2) are valid.For sufficiently small τ, there are the following error estimates: (3.39) Proof.Taking the sum of (3.4d) and (3.4e), we get Using integration by parts, we have we rewrite the L 3 as From (3.29), combining the two inequalities above, we deduce that for sufficiently small τ, taking the sum of (3.55) from 0 to N-1 and using the discrete Gronwall inequality, we can obtain Lemma 6. □ Theorem 4. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)-(3.2) are valid.For sufficiently small τ, there are the following error estimates: for all 1 ≤ n ≤ N.
Tables 1 and 2 verify that (3.4) is the first-order convergence rate O(τ) of (ϕ, σu, u, ρ), which is consistent with the conclusion obtained from theoretical analysis.It is only proved in Theorem 3.2 that the pressure is the half-order convergence rate O(τ 1 2 ) because of technical reasons.However, the numerical results on p in Table 1 still reach the first-order optimal convergence rate O(τ).Tables 3 and  4 show the convergence rate with another set of parameters.  1 shows the evolution of the total energy at τ = 0.02.The downward trend of the energy curve confirms that our scheme is unconditionally energy stable.We also see a downward trend in energy when using different parameters.The energy curves for different time steps are shown in Figure 2 as a result of keeping the other parameters unchanged.It can been found that the curves are very similar, which means that the scheme is robust against different time steps.

Conclusions
To solve the Cahn-Hilliard phase field model for two-phase incompressible flows with variable density, we have designed a novel time marching scheme, which can significantly improve the calculation efficiency.The method is efficient because we decoupled the pressure from the velocity and phase field.We have also proved the unconditional energy stability, presented the error analysis and provided various numerical examples to demonstrate the stability and accuracy of the scheme.In addition, the decoupling method developed in this paper is universally applicable, and this method is always applicable for the generation of an effective fully decoupling scheme.

Use of AI tools declaration
The authors declare they have not used Artificial Intelligence tools (AI) in the creation of this article.