FULLY-DISCRETE FINITE ELEMENT NUMERICAL SCHEME WITH DECOUPLING STRUCTURE AND ENERGY STABILITY FOR THE CAHN–HILLIARD PHASE-FIELD MODEL OF TWO-PHASE INCOMPRESSIBLE FLOW SYSTEM WITH VARIABLE DENSITY AND VISCOSITY

We construct a fully-discrete finite element numerical scheme for the Cahn–Hilliard phasefield model of the two-phase incompressible flow system with variable density and viscosity. The scheme is linear, decoupled, and unconditionally energy stable. Its key idea is to combine the penalty method of the Navier–Stokes equations with the Strang operator splitting method, and introduce several nonlocal variables and their ordinary differential equations to process coupled nonlinear terms. The scheme is highly efficient and it only needs to solve a series of completely independent linear elliptic equations at each time step, in which the Cahn–Hilliard equation and the pressure Poisson equation only have constant coefficients. We rigorously prove the unconditional energy stability and solvability of the scheme and carry out numerous accuracy/stability examples and various benchmark numerical simulations in 2D and 3D, including the Rayleigh–Taylor instability and rising/coalescence dynamics of bubbles to demonstrate the effectiveness of the scheme, numerically. Mathematics Subject Classification. 65N12, 65M12, 65M70. Received February 17, 2021. Accepted September 2, 2021.

However, the fact is that for the variable density model, as far as the author knows, the only numerical technique that can effectively handle the coupling of flow field and phase-field variable and meet the requirements of owning a fully decoupled structure and unconditional energy stability is the stabilization technique developed in [4,25,33]. However, the stabilization technique also has some obvious disadvantages. For example, one has to solve the Cahn-Hilliard equations with variable coefficients (related to 1/ ) at each time step due to the applied stabilization term. This inevitably leads to a higher computational cost than solving a constant-coefficient system. Moreover, in the momentum equation, the scheme developed in [4,25,33] linearizes the advection terms implicitly, which leads to further calculation costs in practice (see Rem. 3.7 for details and Fig. 2c for computational efficiency).
Therefore, on the basis of the scheme developed in [4,25,33], for the flow-coupled two-phase Cahn-Hilliard model with variable density and viscosity, we aim to construct a fully discrete finite element scheme to overcome the shortcomings of stabilization technique. More precisely, we expect that the resulting scheme can solve as many constant coefficient systems as possible, so as to reduce the computational cost as much as possible while maintaining linearity, and owning a decoupling structure and unconditional energy stability. To this end, we note that many nonlinear coupling items satisfy the so-called "zero-energy-contribution" feature (see Rem. 2.1). By combining this particular property with the operator splitting method, we introduce several nonlocal variables and design some special ordinary differential equations (ODE). These variables and ODEs play a vital role in obtaining unconditional energy stability and the decoupling structure of the numerical scheme. We further combine this new decoupling technique with several other efficient approaches, including the penalty method for the fluid momentum equation, the quadratization method for the nonlinear energy potential to arrive at a numerical scheme with the desired properties. The actual implementation of the scheme is also very simple. The calculation of each time step only includes solving several elliptic equations, among which the Cahn-Hilliard equation has constant coefficients, the fluid momentum equation has only two variable coefficient terms, and the pressure Poisson equation also has constant coefficients. The solvability and unconditional energy stability of the proposed scheme are given, and various 2D and 3D numerical examples are further simulated, including Rayleigh-Taylor instability and liquid bubble rising dynamics, to numerically demonstrate the stability and accuracy of the scheme.
We organize the rest of the article as follows. In Section 2, we introduce the variable density/viscosity twophase flow-coupled Cahn-Hilliard model. In Section 3, a numerical scheme is constructed and its unconditional energy stability is strictly proved. The detailed implementation process as well as the solvability of each step are also given. In Section 4, we carry out numerical examples to verify the stability and accuracy. In Section 5, some concluding remarks are given.

The model and its energy structure
The CH-NS model with variable density and viscosity for the two-phase fluid flow system was developed in [1]. Before introducing it, some notations are given here. The 2 inner product of any two functions ( ) and ( ) is denoted by ( , ) = ∫︀ Ω ( ) ( )d , and the 2 norm of ( ) is denoted by ‖ ‖ = ( , ) 1 2 . Suppose that the computational domain Ω ∈ R dim dim = 2, 3 is a smooth, open bounded and connected domain. To represent a two-phase immiscible fluid mixture, a phase-field variable ( , ) is introduced such that: with a thin, smooth transition region with the width ( ) ( ≪ 1) connecting these two distinct values. Note that ∫︀ Ω 1+ 2 d and ∫︀ Ω 1− 2 d can also be viewed as the volume fraction of the fluid 1 and 2, respectively. The density and the viscosity of fluid 1 and 2 are set as 1 , 1 and 2 , 2 , respectively. We assume 1 ≤ 2 without loss of generality. Then, the governing PDE system of the variable density/viscosity two-phase system constructed in [1] reads as: The the initial conditions read as 10) and the boundary conditions read as Here, u is the fluid velocity, is the pressure, ( ) is the density, and ( ) is the viscosity, is the chemical potential, is the surface tension parameter, ( ) = 3 − is the derivative of the nonlinear double-well potential ( ) = 1 4 ( 2 − 1) 2 , and n is the unit outward normal on the boundary Ω. The nonlinear term ∇ is the surface tension term that is derived from the stress tensor, (u · ∇)u and ∇ · (u ) are two fluid advection terms. Note that periodic boundary conditions are another commonly used boundary condition.
The governing system (2.2)-(2.5) holds the law of energy dissipation, which can be derived by performing standard energy estimate and shown by the following theorem.
is the total free energy of the system.
Proof. First, we multiply the 2 inner product of (2.2) with , of (2.3) with − , and perform integration by parts to get Second, we mulitply the 2 inner product of (2.4) with u, use integration by parts, and apply the divergence-free condition (2.5) to obtain Third, from (2.2), (2.7) and (2.8), we obtain By multiplying 1 2 u to (2.17), we derive Hence, we take the 2 inner product of (2.18) with u to derive Finally, by combining (2.14)-(2.16), and (2.19), and using Remark 2.1, we get (2.12).
Remark 2.1. During the process of deriving (2.12), many nonlinear terms are canceled due to integration by parts and the boundary condition u| Ω = 0, including I and II : (∇ · (u ), ) + ( ∇ , u) = 0, (2.20) III and IV : V and VI : Equations (2.20)-(2.22) can be regarded as the "zero-energy-contribution" property [36][37][38][39][40][41], namely, the contribution of all these nonlinear terms to the total free energy of the system is zero. We will use this property to design decoupling type numerical schemes in the next section.

Numerical scheme
We expect to construct an effective time marching scheme that can own the desired properties of decoupling, linearity, and unconditional energy stability. To overcome the challenges on how to discretize the nonlinear terms to obtain the energy stability, we design several nonlocal variables and their associated ODEs and reformulate the PDE system to an equivalent one.

Reformation to an equivalent PDE system
First, we introduce a nonlocal scalar variable ( ), which is defined as where and are two positive constants. The new variable is defined as the square root of the nonlinear double-well potential after extracting a quadratic term 2 2 . This quadratic extraction term can help to maintain the 1 stability of (see [43]). The constant is used to ensure the radicand positive since ( ) − 2 2 is always bounded from below. The use of a new variable to "quadratize" the nonlinear energy potential is the so-called IEQ or SAV method, see [8], which is an efficient method to obtain the linear and stable scheme for gradient flow type models.
Second, we introduce two nonlocal variables ( ) and ( ) and design two special ODEs for them that read as where u| Ω = 0. Note that the design of these ODEs is due the nonlinear terms mentioned in Remark 2.1 satisfying the so-called zero-energy-contribution property. By using the initial conditions and Remark 2.1, it is easy to see that these two ODEs are equivalent to = 0 and = 0 . Thus from their initial conditions, it is easy to see that ( ) = 1 and ( ) = 1 are their exact solutions.
Then, using the new variables , , , the CH-NS system (2.2)-(2.5) can be rewritten to the following equivalent form: where , , are given in (2.7)-(2.9), and ( ) = Since the variables , , are ODEs for time, the boundary conditions of the new system (3.4)-(3.10) are still (2.11), and the initial conditions read as (3.11) Note that there exist some modifications appear in the system (3.4)-(3.10). Here we highlight that even after these modifications, the new system is still equivalent to the original system (2.2)-(2.5). For the momentum equation (3.6), an extra term, 1 2 , are added. This extra term is actually zero from (2.18) and = 1. We also multiply some terms, e.g., ∇ and (u · ∇)u by or . Since = = 1, these modifications will not change the system.
The new system (3.4)-(3.10) also follows the law of energy dissipation, which is derived in the following theorem.
Proof. By multiplying the 2 inner product of (3.4) with , and of (3.5) with − , we get By taking the 2 inner product of (3.6) with u and using integration by parts and the divergence-free condition (3.7), we obtain d d (3.16) By multiplying (3.8) with 2 , (3.9) with , and (3.10) with , we obtain d d By combining the above six equations (3.14)- (3.19) and noting that the two terms marked with the same Roman numerals cancel each other out, we obtain the law of energy dissipation (3.12).
Remark 3.1. Note that we add the constant − − 1 in the above modified energy (3.13) to make it consistent with the original energy (2.13) for the continuous case. This is because in the continuous case, we have (3.1) and ( ) = ( ) = 1. Therefore, after adding the constant − − 1 in (3.13), this modified energy just becomes (2.13).

Numerical scheme
In this subsection, we develop a fully-discrete numerical scheme to solve the system (3.4)-(3.10) which is an equivalent system of (2.2)-(2.5). We denote > 0 as a time step size and = for 0 ≤ ≤ with = .
Let be the numerical approximation to the function (·, )| = . Some finite dimensional discrete subspaces are introduced as follows. Assuming that the polygonal/polyhedral domain Ω is discretized by a conforming and shape regular triangulation/tetrahedron mesh ℎ that is composed by open disjoint elements such thatΩ = ⋃︀ ∈ ℎ¯. We use to denote the space of polynomials of total degree at most and define the following finite element spaces: Besides, we assume the pair of spaces ( ℎ , ℎ ) satisfy the inf-sup condition [14]: where the constant only depends on Ω. A well known inf-sup stable pair ( ℎ , ℎ ) is the Taylor-Hood element [14]. The semi-discrete formulations of the system (3.4)-(3.7) in the weak form reads as: where , , is given in (2.7)-(2.9). By combining the penalty method for the Navier-Stokes equation, the Strang operator splitting method, we construct a fully-discrete scheme as follows: given u ℎ , ℎ , ℎ , , , , Step 2. Find u +1 ℎ ∈ ℎ , and +1 ∈ R such that Step Some explanations of the scheme (3.29)-(3.36) are given through the following remarks.
Remark 3.2. We adopt the penalty method given in [19,25,33] to discretize the Navier-Stokes equations with variable density. In this way, the computation of pressure is not only completely decoupled from the momentum equation, but also avoids solving the elliptic equation with 1/ as the variable coefficient, which can greatly reduce the calculation cost and increase the computational efficiency. Note that we choose −1 ℎ = 0 ℎ in practice. Remark 3.3. Note that we use the implicit and explicit combination method to handle nonlinear terms. For instance, the term (u , ∇ ) in (3.22) is discretized as +1 (u ℎ ℎ , ∇ ℎ ) in (3.29). Thus the scheme is linear.
Remark 3.4. We use the cut-off function to update the density +1 in (3.37), which gives which is quite important to derive the energy stability, see the proof in Theorem 3.2. In fact, the cut-off technique of is a commonly-used technique for the variable density phase-field model, either for the modeling, or the numerical algorithm, see also in [13,27,32,33,44].
Remark 3.5. In step 1, to obtain the intermediate velocityũ +1 ℎ , we split the surface tension term +1 ℎ ∇ ℎ from the momentum equation, which is the so-called first-order operator Strang-splitting method, see also in [4,25,33]. However, there are some slight differences in the discretization of this specific term between our developed scheme and the scheme developed in [4,25,33]. We refer to the detailed comparisons given in Remark 3.7.
The following theorem shows the unconditional energy stability held by the scheme (3.29)-(3.36).
Proof. By taking ℎ = 2 u +1 ℎ in (3.34) and using the following identity (3.42) By taking ℎ = 2ũ +1 in (3.32), we obtain By combining (3.42) and (3.43), we obtain By combining (3.45) and (3.46), we derive (3.47) where we use the following identity: We subtract (3.36) at ( + 1)-step and -step to obtain where we use the Cauchy-Schwarz inequality. Hence (3.50) implies Using the inequality 1 2 ( + ) 2 ≤ 2 + 2 , we obtain Therefore, by combining (3.47), (3.52) and (3.53), we obtain (3.55) (3.57) We multiply (3.31) with 4 +1 to get We multiply (3.33) with 2 +1 to get We multiply (3.35) with 2 +1 to get (3.60) By combining (3.55)-(3.60), we derive Remark 3.6. The energy stability is formally derived for the numerical scheme developed in this paper. Since energy stability guarantees the stability of +1 in 1 , it is expected that the relevant error estimates follow the same reasoning method. Following the error analysis carried out for the different scheme of the Navier-Stokes coupled Cahn-Hilliard phase-field model for the constant density case, cf. [6], the error analysis associated with the developed scheme in this article will be performed in our future work.

Implementation of step 2:
We now give the detailed implementation process of the scheme (3.34) and (3.35) in step 2.
First, by using the nonlocal variable +1 , the unknown velocity field u +1 ℎ can be written as a linear split form that reads as We replace u +1 ℎ using (3.82) and then decompose the obtained equation into the following two sub-equations according to +1 , namely, (3.83) Note that the variable coefficients +1 ℎ + ℎ and +1 ℎ are all positive from (3.38) due to the application of the cut-off technique, which implies these variable-coefficient equations are well-posed. After solving these variable-coefficient elliptic equations, we obtain u +1 1ℎ , u +1 2ℎ . Second, we solve the auxiliary variable +1 from (3.35). By applying the split form of the variable u +1 ℎ given in (3.82) to (3.35), we obtain the following formulation: The solvability of (3.84) can be verified by showing 1 − 2 ̸ = 0, which is obtained by the energy estimate. By taking ℎ = u +1 2ℎ in the second equation of (3.83) and using (3.38), we derive Hence, the total cost at each time step only includes three constant-coefficient elliptic systems ((3.67)-(3.70)); two variable coefficient elliptic equations (3.83); and one pressure Poisson equation (3.36). It can be seen that each step is the realization of decoupling type, which means that the calculation is Remark 3.7. We recall a fully-decoupled type scheme has been developed in [4,25,33]. Here, we give the detailed scheme so that we can compare the scheme developed in [4,25,33] and the scheme (3.29)-(3.36) developed in this article. For simplicity, we only discretize the related terms while other terms remain continuously so that the readers can see the detailed discretization more clearly.
First, compared with (3.32), the scheme given in [4,25,33] discretizes it as Therefore, we getũ The termũ +1 is used as the advective velocity in the Cahn-Hilliard equation which is discretized as By replacingũ +1 in the above using (3.87), we get Therefore, the scheme in [4,25,33] leads to solve the Cahn-Hilliard equation with variable coefficients with a factor 1 , which is very costly in computational costs. Second, for the discretization of the nonlinear terms in the momentum equation, compared with (3.34), the scheme in [4,25,33] discretize them as Although this discrete method can obtain unconditional energy stability, in practice, a large number of variable coefficients bring up unpredictable effects on the condition number of the linear system, leading to higher computational costs. This can be seen in Figure 2c where we compare the average number of iterations per time step needed by our developed scheme (3.29)-(3.36) and the scheme developed in [4,25,33]. Some further works on the convergence analysis of the proposed scheme will be performed by following the lines of [9,10,26]

Numerical examples
In this section, we carry out the developed algorithm (3.29)-(3.36) using the splitting technique through scalar variables (denoted as DSS for short) to investigate the accuracy and energy stability numerically, and also implement several benchmark problems including the Rayleigh-Taylor instability problems and the dynamics of rising droplets.
We set the computational domain as a rectangular region, and use the Taylor-Hood element (see [14]) for ℎ and ℎ that satisfies inf-sup condition and set the finite element spaces (3.20) with 1 = 1, 2 = 2. We set the initial conditions as  To verify the temporal convergence order, we fix mesh size ℎ = 1 256 so that the grid size is sufficiently small and the spatial discretization errors are negligible compared with the time discretization error. Since the exact solutions are not known, we choose the numerical solutions obtained with a very tiny time step size = 1 − 9 computed by the scheme DSS as the exact solution approximately for computing errors. In Figure 1a, the 2errors between the numerical solution and the exact solution at = 0.5 are plotted (the exact solution of and is 1). We observe that all variables follow the first-order time accuracy.

Accuracy and stability test
To verify the spatial convergence order, in Figure 1b, we plot the errors of DSS for various spatial grid size ℎ. We choose sufficiently small ( = 1 − 9) so that the errors are only dominated by the spatial discretization error. We use the numerical solutions obtained with a very tiny mesh size = 1 − 9, 1 ℎ = 1 512 computed by the scheme DSS as the exact solution approximately for computing errors. The errors in various norms between the numerical solution and the exact solution at = 0.2 are plotted. We can see that the second-order convergence rates are followed by the 1 -error of the velocity, 2 -error of the pressure , and 2 -error of the phase-field variable , while third-order convergence rates are observed for the 2 -error of the velocity. These results are in full agreements with the theoretical expectation of accuracy for 2/ 1 element for (u, ) and 1 element for .
We further perform the energy stability tests by plotting the total free energy (2.13) evolved over time in Figure 2a by using various time steps. The decays of all energy curves illustrates the energy stability. In Figure 2b, using = 0.01, we plot the time evolution of the total free energy (2.13) (in the original form) and (3.40) (in the modified discrete form). The overlapping of these two energy curves show that these two energies are consistent numerically.
In Figure 2c, we compare the average number of iterations per time step needed by DSS and that needed by the scheme developed in [4,25,33] (see Rem. 3.7, denoted by Stab for short). Note that the time of each iteration is basically the same as the time to solve a Poisson equation with constant coefficients. The comparison shows that the efficiency of DSS is much higher than that of Stab, because DSS needs to solve far fewer variablecoefficient equations. In Figure 2d, we compare the 2 error of when using various stabilization parameters . When = 0, the scheme is not stable for larger time steps, so the error points are missed. When = 4, the accuracy curve is not only very close to the results of = 0 for small time steps, but also stable for all other tested big time steps, and the entire curve shows the first-order convergence rate. When = 20 and = 100, we can see that even though the scheme is always stable, the accuracy is much worse than the case of = 4. Therefore, we adopt = 4 in all our computations.

Rayleigh-Taylor instability
Rayleigh-Taylor interface instability is a benchmark phenomenon of two-phase flow with variable density, see [11,18,24,34]. The classic experimental setup is that the system consists of two layers of fluids with different densities. Under the action of gravity, slight interference caused by potential energy will cause the heavier fluid to sink, and the lighter fluid will rise, resulting in irregularities (ripples) caused by downward movement will quickly expand into a series of "RT fingers" (shown in Figs. 3 and 5).    We set the 2D computational domain as ( , ) ∈ Ω = [0, ]×[0, 2 ] and add the gravity field to the momentum equation, which reads as where g = (0, 0 ) and 0 is the gravity constant. We rescale the variables/parameters by the following way: We set the initial conditions of a perturbed two-layer fluid as follows, , u 0 = (0, 0), 0 = 0. Therefore, the Reynolds number is 1000.
In the first simulation, we set the density of the heavier fluid as 2 = 3 to investigate the dynamical changes of the fluid interface when the density ratio of the two fluids is not much different. The so-called Atwood number is given by Tryggvason's definition in [34], which is = max− min max+ min . For this case, = 0.5. Snapshots of the profiles of at various times are shown in Figure 3. While the upper heavier fluid flows to the lower layer under the action of the gravity field, a roll-like structure is formed. In the initial stage of the production of the rolling structure, as the main vortex scrolls, some small branches (R-T fingers) are formed, see = 3.5, and more complex fingering patterns are produced gradually. We further plot the profiles of the velocity field at various times in Figure 4. In the second simulation, we increase the density of the upper heavier fluid as 2 = 100, namely, = 0.9802. From the snapshots of at different times shown in Figure 5 and the velocity field profiles shown in Figure 6, we observe that more complex vortices and R-T fingers appear during the formation of the rolling structure. These obtained simulations are qualitatively consistent with numerical simulations and experimental results using different models and different numerical methods, see [11,18,24,34].

Dynamics of rising droplets
In this numerical example, we examine the rising phenomenon of 3D light droplets immersed in a fluid with a larger density due to the action of the gravitational field. The momentum equation (4.3) equipped with the gravity force (g = (0, 0, 0 )) is still adopted. The non-dimensionalized method is also consistent with the previous example. We first study the dynamic changes of the rising deformation of a single droplet. The initial droplet is set at the bottom of the domain with the following formulation: Snapshots of the isosurface { = 0} at different times are plotted in Figure 7. Due to the influence of gravity, the shape of the droplet will change dramatically as it rises gradually, from the initial spherical shape to a flat cap shape. Note that the results of the deformation we obtained are consistent with the physical experimental results given in [3] and the simulation results of different models given in [2]. Next, we set the initial conditions to two non-contact droplets located at different latitudes to study the rising dynamics of multiple droplets under the influence of gravity and examine their mutual influence. The model parameters are still given in (4.7).
In the first simulation, the initial two droplets of the same size are placed at different heights, one of which is suspended directly above the other. This is the so-called coaxial coalescence situation, see [5]. The formulation  Figure 8a, where we can see that the deformations of the upper and lower droplets under the influence of gravity are completely different. At the initial stage, both droplets will change from a spherical shape to a cap shape, but the height of the cap formed by the upper droplet is very short, and the cross section of the cap becomes very large. However, the height of the cap formed by the droplet below is very large, and the cross section becomes relatively small. And the moving speeds of the upper and lower droplets are also different. The lower droplet moves faster, so it will catch up with the higher droplet, and then they merge together (approximately at = 1.4). For comparison, we show the experimental results of [5] in Figure 8b which can be seen to be very consistent with our numerical simulation qualitatively.
In the second simulation, we also use two droplets of the same size. and place them at different heights, but one droplet is in a staggered position on the other droplet. This is the so-called oblique coalescence case, cf. [5]. The initial conditions are still given in (4.9) with 1 = 0.3, 2 = 0.5. The isosurface of { = 0} at different times are plotted in Figure 9a. We can see that when the two droplets start to rise upward, the obvious difference Figure 9. (a) The dynamical motion of two 3D rising droplets (oblique coalescence case) in the gravity field where 1 = 1 and 2 = 1000. Profiles of the isosurfaces of the interface { = 0} are taken at = 0, 1.7, 1.9, 2.1, 2.3, and 2.9; (b) the experimental results by Brereton and Korotney in [5].
from the coaxial simulation is that their staggered position makes the two droplets become a hat shape with non-symmetrical structure. As they rise gradually, the lower droplet moves faster, while the upper droplet moves slowly, so the lower droplet will eventually be captured by the upper droplet ("fish-capture" phenomenon), and then the two droplets will be completely fuse together to form a large, obviously inclined cap-shaped drop. Similarly, we also compare with the experimental results given [5] shown in the Figure 9b. It can be seen that for the oblique coalescence case, our numerical simulation is also very consistent with the experimental results qualitatively.

Concluding remarks
The purpose of this paper is to construct a new fully discrete finite element numerical scheme for solving the Cahn-Hilliard-Navier-Stokes phase-field model with variable density and viscosity. With the help of multiple nonlocal variables and some specially designed ODEs, this scheme has the characteristics of linearity, full decoupling, and easy implementation. The calculation cost of each step only includes solving a few linear elliptic equations, and most of them are constant coefficients. We strictly show the solvability and energy stability and give the implementation process in detail. Various 2D and 3D numerical examples are performed to demonstrate the accuracy and efficiency. Some benchmark simulations, such as Rayleigh-Taylor instability and drop rising dynamics, are also carried out, and the good agreement with physical experiments proves the effectiveness of the scheme.