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

1.   Introduction

In this paper, we focus on the following hydrodynamically consistent Cahn-Hilliard phase field model for incompressible two-phase flows with variable density

where $\phi$, $w$, $\boldsymbol{u}$ and $p$ are the phase field function, the chemical potential, the velocity of flow and pressure respectively. Moreover, $\rho = \frac{1}{2}(\rho_1+\rho_2)+\frac{\phi}{2}(\rho_1-\rho_2)$ is the density, where $\rho_1$ and $\rho_2$ are the densities of the two fluids. The parameters $\gamma$, $\varepsilon$, $\eta$ and $\lambda$ are the mobility parameter related to the relaxation time scale, the interface thickness ^[1], the viscosity of the field, and the mixing energy density respectively. The system (1.1) is supplemented with appropriate boundary and initial conditions which are given as

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 two-phase 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: (ⅰ) Cahn-Hilliard models can accurately conserve the volume and dynamics; (ⅱ) 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.

2.   The discrete scheme

For the sake of simplicity, some notations are needed for the following content. We assume that the domain $\Omega\in R^2$ is open, sufficiently smooth and bounded. For any two functions $\phi(x)$ and $\psi(x)$, their $L^2$ inner product on $\Omega$ is denoted by $(\phi, \psi) = \int_\Omega \phi(x)\psi(x)dx$, and the $L^2$ norm of $\phi(x)$ is denoted by $\Vert\phi\Vert = (\phi, \phi)^{\frac{1}{2}}$. Let $\tau > 0$ be the time step size and set $t^n = n\tau$ for $0\leq n \leq N$ with $T = N\tau$. Moreover, we introduce the following spaces,

where ${\bf{n}}$ is the outward normal vector of $\partial\Omega$. Next, we reformulate the system (1.1) as follows:

Remark 1. In (2.1d), the term $\frac{1}{4} \rho^{n+1}\left(\nabla \cdot \boldsymbol{u}^{n}\right) \tilde{\boldsymbol{u}}^{n+1}$ is added to obtain the unconditional stability, as it is 0 if $\nabla\cdot \boldsymbol{u}^n = 0$.

Theorem 1. $(Stability\; of\; \rho)$ For any $\tau > 0$ and any sequence $\left\{\boldsymbol{u}^{n}\right\}_{n = 0, \ldots, N}$ satisfying the boundary condition $\boldsymbol{u}^{n} \cdot \boldsymbol{n} = 0$ on $\partial \Omega$, the solution $\left\{\rho^{n}\right\}_{n = 1, \ldots, N}$ to (2.1a) satisfies

Proof. Testing (2.1a) with $2 \tau \rho^{n+1}$ gives

Owing to the boundary conditions on $\boldsymbol{u}^{n}$, we note that

Thus, we get

Summing all indices $\mathrm{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 $\frac{1}{2}\left\| \sqrt{\rho^n}\boldsymbol{u}^n\right\|^2$, it is more suitable to establish a bound based on $\left\| \sqrt{\rho^n}\boldsymbol{u}^n\right\|^2$ than on the velocity itself. For simplicity, let us say that $\sigma^n = \sqrt{\rho^n}$ for all $1\leq n\leq N$ and $\sigma_0 = \sqrt{\rho_0}$.

Theorem 2. $(Stability\; of\; energy)$ For any $\tau > 0$, (2.1b)–(2.1e) satisfy the the conditions of following energy estimates:

Proof. Multiplying (2.1b) by $2 \tau \lambda w^{n+1}$ and integrating over $\Omega$, we get

Multiplying (2.1c) by $2 \lambda\left(\phi^{n+1}-\phi^{n}\right)$ yields

where we use the following identity:

Testing (2.1a) with $\tau\vert\tilde{\boldsymbol{u}}^{n+1}\vert^{2}$ leads to

By taking the inner product of (2.1d) with $2 \tau \tilde{\boldsymbol{u}}^{n+1}$, we have

Testing (2.1e) with $2 \tau \boldsymbol{u}^{n+1}$ yields

Taking into account the boundary condition on $\boldsymbol{u}^{n}$ and using integration by parts, we have

Summing the above inequality, we arrive at

Adding up the above inequality from $n = 0$ to $N-1$, we obtain Theorem 2.3. □

The bound on the pressure $p$ is proved in the following theorem.

Theorem 3. $(Stability\; of\; p)$ For any $\tau > 0$, the solution $p^{n+1}$ to (2.1e) satisfies the following inequality:

Proof. Under the inf-sup condition, there exists a positive constant $\beta$ such that

Testing (2.14) with all $v \in V$ leads to

Given that $\left\|\rho^{n}\right\| \leq\left\|\rho_{0}\right\|$ for all $1 \leq n \leq N$, and by using Hölder's inequality, we have

Then, by the Sobolev embedding inequality $\|v\|_{L^{6}} \leq C\|\nabla v\|_{L^{2}}$ for any $v \in V$, we get

Substituting the above inequalities into Eq (2.14), we obtain

And by using Theorem 2.2, we get the desired result. □

3.   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 $\chi$ is a number in $(0, \rho_0^{min}]$.

We assume that the exact solution $(\rho, \boldsymbol{u}, \phi, 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:

for all $0 \leq n \leq N-1$.

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_\phi^{n+1}$. For $R_u^{n+1}$, we can rewrite

Using $R_\rho^{n+1}$ estimation and Hölder's inequality can yield the result for $R_u^{n+1}$. □

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 $\gamma_{k}$, for integers $k\geq0$, be the nonnegative numbers such that

Suppose that $\tau\gamma_{k} < 1$, for all k, and set $\sigma_{k} = (1-\tau\gamma_{k})^{-1}$. Then,

We verify Lemma 5 by the following lemma.

Lemma 3. Define

Then, for n$< \frac{T}{\tau}-1$, we have

Proof. For $\left\|G_{c}^{n+1}\right\|$, we use (3.8) to conclude that

where we have used the $a$ priori bound $\left\|\phi^{n}\right\|_{H^{1}} \leq 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 $\rho^{n+1}$ is derived in the following lemma.

Lemma 4. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)–(3.2) hold. Then, we have

for all $0 \leq n \leq N-1$.

Proof. Multiplying (3.4a) by $2 \tau e_{\rho}^{n+1}$ and integrating over $\Omega$, we have

Using $\nabla \cdot \boldsymbol{u}\left(t_{n+1}\right) = 0$ in $\Omega$ and $\boldsymbol{u}\left(t_{n+1}\right) = 0$ on $\partial \Omega$, we have

By using Young's inequality, the Cauchy-Schwarz inequality and Lemma 1, we have

where $\varepsilon > 0$ is a sufficiently small constant.

Then, by the Sobolev inequality and Young inequality, we get

Similarly, the last term can be estimated as follows:

If $\varepsilon$ is sufficiently small such that $\varepsilon \tau \leq \frac{1}{6}$, substituting the estimates of $K_{i}\ (1 \leq i \leq 4)$ into (3.12), we have

Using the discrete Gronwall inequality, we obtain the desired result. □

Lemma 5. 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 $\tau$, there are the following error estimates:

Proof. Let us multiply (3.4a) by $\tau\vert\tilde{e}_{u}^{n+1}\vert^{2}$, (3.4b) by $2 \tau e_{\phi}^{n+1}$ and $2 \tau \lambda e_{w}^{n+1}$, (3.4c) by $2\tau\gamma e_{w}^{n+1}$ and $-2 \lambda(e_{\phi}^{n+1}-e_{\phi}^{n})$, (3.4d) by $2 \tau \tilde{e}_{u}^{n+1}$ and (3.4e) by $2 \tau e_{u}^{n+1}$. Summing up all of the above equations, we have

Thanks to $\nabla \cdot {\bf u}^{n} = 0$ in $\Omega$, we have

Therefore, we get

where the Sobolev embedding inequality $\|v\|_{L^{6}} \leq C\|\nabla v\|_{L^{2}}$ for any $v \in V$ is used.

For $A_4$, we have

Using the $Poincar\acute{e}$ inequality and Cauchy-Schwarz inequality, we can deal with $A_6$, $A_7$ and $A_9$:

and

From Lemma 1, we have

From Lemma 3, we obtain

and

where we use this identity

We denote

Similar to the method estimated from Lemma 3, we can obtain

Taking the gradient of $\tilde{G}^{n+1}$, we get

Since $H^{2}(\Omega) \subset L^{\infty}(\Omega)$ and by utilizing the bound of $\left\|\nabla e_{\phi}^{n+1}\right\|_{L^{2}}$ implied by Theorem 2, we conclude that

In view of (3.10) and the bound $\left\|\phi^{n}\right\|_{H^{1}} < C$, we have

Dealing with the penultimate term of (3.28) and in view of (3.31) and (3.33), we get

Then, we estimate the last term of (3.28),

Combining (3.35) and (3.36), we can obtain

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(\tau^{\frac{1}{2}})$. 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 $\tau$, there are the following error estimates:

Proof. Taking the sum of (3.4d) and (3.4e), we get

Let us multiply (3.4a) by $\tau\tilde{e}_{u}^{n+1}$, (3.4b) by $2 \tau e_{\phi}^{n+1}$ and $2 \tau \lambda e_{w}^{n+1}$, (3.4c) by $2\tau\gamma e_{w}^{n+1}$ and $-2 \lambda(e_{\phi}^{n+1}-e_{\phi}^{n})$ and (3.40) by $2 \tau e_{u}^{n+1}$. Summing up all of the above equations, we obtain

Using integration by parts, we have

we rewrite the $L_{3}$ as

then,

Lemma 5 shows that

According to Young's inequality and the Cauchy-Schwarz inequality, we can deduce that

Therefore, we have

Using the embedding inequality and Cauchy-Schwarz inequality, we can deal with $M_4$–$M_9$

and

Next, we estimate $M_{10} + M_{11}$ and $M_{12} + M_{13}$ as follows:

and

From (3.29), combining the two inequalities above, we deduce that

Plugging the above inequality into (3.41) gives

for sufficiently small $\tau$, 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 $\tau$, there are the following error estimates:

for all $1 \leq n \leq N$.

Proof. From Lemmas 5 and 6, we obtain

Thus, we only prove that

In fact, we have

□

Theorem 3.6 states that both $\sigma^{n}\boldsymbol{u}^{n}, \; \rho^{n}, \; \boldsymbol{u}^{n}$, $\phi^{n}$ and $w^{n}$ are order 1 approximations to $\sigma\boldsymbol{u}, \; \rho, \; \boldsymbol{u}$, $\phi$ and $w$ in $l^{\infty}\left(L^{2}(\Omega)\right)$, $l^{\infty}\left(L^{2}(\Omega)\right)$, $l^{2}\left(H_{0}^{1}(\Omega)\right)$, $l^{\infty}\left(H_0^{1}(\Omega)\right)$ and $l^{\infty}\left(H_0^{1}(\Omega)\right)$, respectively. Finally, we can obtain order $\frac{1}{2}$ error estimates for $p$ approximation in $l^{\infty}\left(L^{2}(\Omega)\right)$.

Theorem 5. Under the assumptions in Theorem 4, the following holds true:

Proof. Let us rewrite (3.40) as

To prove the theorem, we introduce the discrete inf-sup condition, i.e.,

Then, we can constrain the products of the right-hand side of (3.61) with an arbitrary $v \in V$ as follows:

and

For the second term on the right-hand side, using $\|v\|_{L^{3}} \leq\|v\|^{\frac{1}{2}}\|\nabla v\|^{\frac{1}{2}}$, we have

By the split method, we arrive at

Utilizing H$\ddot{o}$lder's inequality and Young's inequality, we derive

By adding up the inequalities and incorporating (3.62), we get the desired result. □

4.   Numerical results

In this section, we will present some numerical experiments to prove the validity and accuracy of our method. In the following simulation, for phase field $\phi$, chemical potential $w$ and pressure $p$, we take the $P1$ finite element space (continuous piecewise linear), and for fluid velocity $\boldsymbol{u}$, we take the $P2$ finite element space. All experiments were conducted in Freefem++. We fixed $\eta$ = 0.8, $\lambda$ = 0.7, $\gamma$ = 0.0006, $\varepsilon$ = 0.1, $T$ = 0.1, $\rho_1$ = 1 and $\rho_2$ = 3. The computational domain and the initial conditions were taken as

Tables 1 and 2 verify that (3.4) is the first-order convergence rate $O(\tau)$ of $(\phi, \sigma\boldsymbol{u}, \boldsymbol{u}, \rho)$, 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(\tau^{\frac{1}{2}})$ because of technical reasons. However, the numerical results on $p$ in Table 1 still reach the first-order optimal convergence rate $O(\tau)$. Tables 3 and 4 show the convergence rate with another set of parameters.

Figure 1 shows the evolution of the total energy at $\tau = 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.

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

Acknowledgments

This work was supported by the Research Project Supported by Shanxi Scholarship Council of China (No. 2021-029), Shanxi Provincial International Cooperation Base and Platform Project (202104041101019), Shanxi Province Natural Science Research (202203021211129), Shanxi Province Natural Science Research (No. 202203021212249) and Special/Youth Foundation of Taiyuan University of Technology (No. 2022QN101).

Conflict of interest

All authors declare that they have no competing interests in this paper.