Mass Conservative and Energy Stable Finite Difference Methods for the Quasi-incompressible Navier-Stokes-Cahn-Hilliard system: Primitive Variable and Projection-Type Schemes

TL;DR

This paper introduces two finite difference schemes for the quasi-incompressible Navier-Stokes-Cahn-Hilliard system that are mass conservative and energy stable, with efficient solvers validated through various fluid flow simulations.

Contribution

The paper presents two novel finite difference methods that ensure mass conservation and energy stability for the q-NSCH system, along with an efficient multigrid solver.

Findings

Mass conserved up to 10^-8 in simulations.

Energy decreases monotonically over time.

Methods accurately reproduce analytical and previous numerical results.

Abstract

In this paper we describe two fully mass conservative, energy stable, finite difference methods on a staggered grid for the quasi-incompressible Navier-Stokes-Cahn-Hilliard (q-NSCH) system governing a binary incompressible fluid flow with variable density and viscosity. Both methods, namely the primitive method (finite difference method in the primitive variable formulation) and the projection method (finite difference method in a projection-type formulation), are so designed that the mass of the binary fluid is preserved, and the energy of the system equations is always non-increasing in time at the fully discrete level. We also present an efficient, practical nonlinear multigrid method - comprised of a standard FAS method for the Cahn-Hilliard equation, and a method based on the Vanka-type smoothing strategy for the Navier-Stokes equation - for solving these equations. We test the scheme in the context of Capillary Waves, rising droplets and Rayleigh-Taylor instability. Quantitative comparisons are made with existing analytical solutions or previous numerical results that validate the accuracy of our numerical schemes. Moreover, in all cases, mass of the single component and the binary fluid was conserved up to 10 to -8 and energy decreases in time.

Tables

Table 1. Table 1 : Convergence test for the primitive method and projection method.

	Primitive		Projection
$m$	error	rate	error	rate
16 & 32	6.063$\times {10}^{-3}$	—	6.409$\times {10}^{-3}$	—
32 & 64	2.310$\times {10}^{-3}$	1.40	2.603$\times {10}^{-3}$	1.30
64 & 128	6.812$\times {10}^{-4}$	1.76	7.814$\times {10}^{-4}$	1.73
128 & 256	1.622$\times {10}^{-4}$	2.06	1.917$\times {10}^{-4}$	2.03

Table 2. Table 2 : Average computational cost for the Primitive method and Projection method.

	Ave. CPU time per time step (sec)
	Primitive	Projection
16	2.913$\times {10}^{-2}$	1.075$\times {10}^{-2}$
32	1.659$\times {10}^{-1}$	3.816$\times {10}^{-2}$
64	2.659$\times {10}^{-1}$	1.149$\times {10}^{-1}$
128	1.027$\times {10}^0$	4.170$\times {10}^{-1}$
256	2.622$\times {10}^0$	1.246$\times {10}^0$

Equations

$$\rho {\mathbf{u}}_t+\rho \mathbf{u}\cdot \nabla \mathbf{u}=$$

$$\displaystyle+\boldsymbol{\nabla}\big{(}\lambda(c)(\boldsymbol{\nabla}\cdot\boldsymbol{u})\big{)}-\rho g\boldsymbol{j},$$

$$\nabla \cdot \mathbf{u}=$$

$$\rho c_t+\rho \mathbf{u}\cdot \nabla c=$$

$${\mu }_c=$$

$$\mathbf{u}{|}_{\partial \Omega }={\mathbf{u}}_{\partial \Omega },$$

$$\mathbf{n}\cdot \nabla c{|}_{\partial \Omega }=\mathbf{n}\cdot \nabla {\mu }_c{|}_{\partial \Omega }=0,$$

$$\alpha \rho c_t+\alpha \rho \mathbf{u}\cdot \nabla c=\nabla \cdot \mathbf{u}$$

$${\rho }_t+\nabla \cdot (\rho \mathbf{u})=0,$$

$$\rho {\mathbf{u}}_t+\rho \mathbf{u}\cdot \nabla \mathbf{u}=$$

$$\displaystyle+\cfrac{1}{3Re}\boldsymbol{\nabla}\big{(}\mu(c)(\boldsymbol{\nabla}\cdot\boldsymbol{u})\big{)}-\cfrac{\rho}{Fr}\boldsymbol{j},$$

$$\nabla \cdot \mathbf{u}=$$

$$\rho c_t+\rho \mathbf{u}\cdot \nabla c=$$

$${\mu }_c=$$

$$\rho {\mathbf{u}}_t+\rho (\mathbf{u}\cdot \nabla )\mathbf{u}+\frac{1}{2}{\rho }_t\mathbf{u}+\frac{1}{2}\nabla \cdot (\rho \mathbf{u})\mathbf{u}=-\frac{1}{M}\nabla \bar{p}+\frac{1}{M}\rho {\bar{\mu }}_c\nabla c$$

$$\displaystyle+\cfrac{1}{Re}\boldsymbol{\nabla}\cdot(\mu(c)\boldsymbol{\nabla}\boldsymbol{u})+\cfrac{1}{3Re}\boldsymbol{\nabla}\big{(}\mu(c)(\boldsymbol{\nabla}\cdot\boldsymbol{u})\big{)}-\cfrac{\rho}{Fr}\boldsymbol{j},$$

$$\displaystyle\boldsymbol{\nabla}\cdot\boldsymbol{u}=\cfrac{\alpha}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c)\boldsymbol{\nabla}\bar{\mu}_{c}\big{)}+\cfrac{\alpha^{2}}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c)\boldsymbol{\nabla}\bar{p}\big{)},$$

$$\displaystyle\rho c_{t}+\rho(\boldsymbol{u}\cdot\boldsymbol{\nabla})c=\cfrac{1}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c)\boldsymbol{\nabla}\bar{\mu}_{c}\big{)}+\cfrac{\alpha}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c)\boldsymbol{\nabla}\bar{p}\big{)},$$

$$\rho {\bar{\mu }}_c=\frac{M\eta }{ϵWe}\rho f(c)+\frac{M\eta }{ϵWe}\frac{\partial \rho }{\partial c}F(c)+\frac{M\eta ϵ}{We}\frac{\partial \rho }{\partial c}\frac{1}{2}(\nabla c\cdot \nabla c)-\frac{M\eta ϵ}{We}\nabla \cdot (\rho \nabla c).$$

$$\bar{p}$$

$${\int }_{\Omega }{\rho }_t\mathrm{d}\mathbf{x}$$

$${\rho }_t+\nabla \cdot (\rho \mathbf{u})=0.$$

$$\displaystyle(\rho c)_{t}+\boldsymbol{\nabla}\cdot(\rho\boldsymbol{u}c)=\cfrac{1}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c)\boldsymbol{\nabla}\bar{\mu}_{c}\big{)}+\cfrac{\alpha}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c)\boldsymbol{\nabla}\bar{p}\big{)}=\cfrac{1}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c)\boldsymbol{\nabla}\mu_{c}\big{)}.$$

$$\frac{dE}{dt}=$$

$$=$$

$$\displaystyle\cfrac{d}{dt}\bigg{(}\cfrac{1}{2}\rVert\sqrt{\rho}\boldsymbol{u}\rVert^{2}_{L^{2}}+\int_{\Omega}\big{(}\cfrac{1}{Fr}\rho y\big{)}{\rm d}\boldsymbol{x}\bigg{)}$$

$$-\frac{1}{Re}\Vert \sqrt{\mu (c)}\nabla \mathbf{u}{\Vert }_{L^2}^2-\frac{1}{3Re}\Vert \sqrt{\mu (c)}\nabla \cdot \mathbf{u}{\Vert }_{L^2}^2,$$

$${\int }_{\Omega }\rho \mathbf{u}\cdot \mathbf{j}$$

$$\displaystyle\int_{\Omega}\bigg{(}\rho(\boldsymbol{u}\cdot\boldsymbol{\nabla})\cfrac{1}{2}(\boldsymbol{u}\cdot\boldsymbol{u})+\cfrac{1}{2}(\boldsymbol{u}\cdot\boldsymbol{u})\boldsymbol{\nabla}\cdot(\rho\boldsymbol{u})\bigg{)}{\rm d}\boldsymbol{x}=0.$$

$$\displaystyle\int_{\Omega}0~{}{\rm d}\boldsymbol{x}=\int_{\Omega}\bigg{(}-\cfrac{1}{M}(\boldsymbol{\nabla}\cdot\boldsymbol{u})\bar{p}-\cfrac{\alpha}{PeM}m(c)\boldsymbol{\nabla}\bar{\mu}_{c}\cdot\boldsymbol{\nabla}\bar{p}-\cfrac{\alpha^{2}}{PeM}m(c)\boldsymbol{\nabla}\bar{p}\cdot\boldsymbol{\nabla}\bar{p}\bigg{)}{\rm d}\boldsymbol{x}.$$

$$0$$

$$\displaystyle-\cfrac{\alpha}{PeM}m(c)\boldsymbol{\nabla}\bar{\mu}_{c}\cdot\boldsymbol{\nabla}\bar{p}\bigg{)}{\rm d}\boldsymbol{x}.$$

$$\displaystyle\cfrac{d}{dt}\bigg{(}\cfrac{\eta\epsilon}{2We}\rVert\sqrt{\rho}~{}\boldsymbol{\nabla}c\rVert^{2}_{L^{2}}+\int_{\Omega}\big{(}\cfrac{\eta}{\epsilon We}\rho F(c)\big{)}{\rm d}\boldsymbol{x}\bigg{)}=\int_{\Omega}\cfrac{1}{M}\rho c_{t}\bar{\mu}_{c}{\rm d}\boldsymbol{x}.$$

$$\displaystyle\rho^{n}\cfrac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^{n}}{\delta t}+\rho^{n}(\boldsymbol{u}^{n}\cdot\boldsymbol{\nabla})\boldsymbol{u}^{n+1}+\cfrac{1}{2}\boldsymbol{u}^{n+1}\big{(}\cfrac{\rho^{n+1}-\rho^{n}}{\delta t}+\boldsymbol{\nabla}\cdot(\rho^{n}\boldsymbol{u}^{n})\big{)}=-\cfrac{1}{M}\boldsymbol{\nabla}\bar{p}^{n+1}$$

$$\displaystyle+\cfrac{1}{M}\rho^{n}{\bar{\mu}_{c}}^{n+1}\boldsymbol{\nabla}c^{n}+\cfrac{1}{Re}\boldsymbol{\nabla}\cdot(\mu(c^{n})\boldsymbol{\nabla}\boldsymbol{u}^{n+1})+\cfrac{1}{3Re}\boldsymbol{\nabla}\big{(}\mu(c^{n})(\boldsymbol{\nabla}\cdot\boldsymbol{u}^{n+1})\big{)}-\cfrac{\rho^{n}}{Fr}\boldsymbol{j},$$

$$\displaystyle\boldsymbol{\nabla}\cdot\boldsymbol{u}^{n+1}=\cfrac{\alpha}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c^{n})\boldsymbol{\nabla}{\bar{\mu}_{c}}^{n+1}\big{)}+\cfrac{\alpha^{2}}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c^{n})\boldsymbol{\nabla}\bar{p}^{n+1}\big{)},$$

$$\displaystyle\rho^{n+1}\cfrac{c^{n+1}-c^{n}}{\delta t}+\rho^{n}\boldsymbol{u}^{n+1}\cdot\boldsymbol{\nabla}c^{n}=\cfrac{1}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c^{n})\boldsymbol{\nabla}{\bar{\mu}_{c}}^{n+1}\big{)}$$

$$\displaystyle+\cfrac{\alpha}{Pe}\boldsymbol{\nabla}\cdot\big{(}m(c^{n})\boldsymbol{\nabla}\bar{p}^{n+1}\big{)},$$

$${\rho }^{n+1}{{\bar{\mu }}_c}^{n+1}=\frac{M\eta }{ϵWe}{\rho }^{n+\frac{1}{2}}g(c^{n+1},c^n)+\frac{M\eta }{ϵWe}{F}^{n+\frac{1}{2}}(c)r(c^{n+1},c^n)$$

$$+\frac{ϵ\eta M}{2We}(\nabla c\cdot \nabla c{)}^{n+\frac{1}{2}}r(c^{n+1},c^n)-\frac{ϵ\eta M}{We}\nabla \cdot ({\rho }^{n+\frac{1}{2}}\nabla c^{n+\frac{1}{2}}),$$

$$\displaystyle g(c^{n+1},c^{n})=\cfrac{1}{4}\big{(}c^{n+1}(c^{n+1}-1)+c^{n}(c^{n}-1)\big{)}(c^{n+1}+c^{n}-1)$$

Full text

††journal: Computer Methods in Applied Mechanics and Engineering

1 Introduction

Phase-field, or diffuse-interface models [4, 29], have now emerged as a powerful method to simulate many types of multiphase flows, including drop coalescence, break-up, rising and deformations in shear flows [19, 26, 30, 31, 34, 52], contact line dynamics [14, 15, 27], thermocapillary effects [18, 20], and tumor growth [24, 38]. Phase-field model are based on models of fluid free energy which goes back to the work of Cahn and collaborators [6, 7]. The basic idea is to introduce a phase variable (order parameter) to characterize the different phases that varies continuously over thin interfacial layers and is mostly uniform in the bulk phases. Sharp interfaces are then replaced by the thin but nonzero thickness transition regions where the interfacial forces are smoothly but locally distributed in the bulk fluid. One set governing equations for the whole computational domain can be derived variationally from the free energy, where the order parameter fields satisfy an advection-diffusion equation (usually the advective Cahn-Hilliard equations) and is coupled to the Navier-Stokes equations through extra reactive stresses that mimic surface tension.

The classical phase-field model, the Model H [25], was initially developed for simulating a binary incompressible fluid where components are density matched, and was later generalized for simulating binary incompressible fluids with variable density components [1, 2, 5, 11, 13, 28, 36, 41, 43], in which some models, however, do not satisfy the Galilean invariance or are not thermodynamic consistency. As the phase-field model can be derived through a variational procedure, thermodynamic consistency of the model equations can serve as a justification for the model. In addition, this approach ensures the model compatible with the laws of thermodynamics, and to have a strict relaxational behavior of the free energy, hence the models are more than a phenomenological description of an interfacial problem. Lowengrub and Truskinovsky [36] and Abels et al. [1] extended the Model H to a thermodynamically consistent model for variable density using two different modelling assumptions on the phase variable (mass concentration [36] or volume fraction [1]) and the velocity field (mass averaged [36] or volume averaged velocity [1]). Although the two models are developed to represent the same type of flow dynamics, the resulting equations have significant differences due to the underlying modeling choices. In particular, the quasi-incompressible NSCH model (q-NSCH) developed by Lowengrub and Truskinovsky [36] adopts a mass-averaged velocity, and the fluids are mixing at the interfacial region which generates the changes in density. Such a system was called quasi-incompressible, which leads to a (generally) non-solenoidal velocity field ($\boldsymbol{\nabla}\cdot\boldsymbol{u}\neq 0$ but was given through the quasi-incompressibility condition) and an extra pressure term appears in the Cahn-Hilliard equation comparing to Model H. In the model of Abels et. al. [1], a solenoidal (divergence-free) velocity field is obtained due to the volume-averaged mixture velocity modeling assumption. However the mass conservation equation of their model is modified by adding a mass correction term. Most recently, another quasi-incompressible phase-field model [51] was developed to study the binary fluid with variable density, where the volume fraction is employed as the phase variable leading to a different free energy.

Solving the q-NSCH model is quite a challenging problem. The CH equation is a fourth order nonlinear parabolic PDE, which contains an extra pressure term; the solution of the phase variable varies sharply through the thin diffuse-interface region where the velocity field is non-solenoidal; the variable density is a non-linear function of the phase variable; the NS and CH equations are strongly coupled, which further increases the mathematical complexity of the model and that makes it difficult to design provably stable numerical schemes. Recently, it has been reported that thermodynamic consistency can serve as not only a critical justification for the phase-field modeling, but also an important criterion for the design of numerical methods. When the thermodynamic consistency is preserved at the discrete level, it guarantees the energy stability of the numerical method and also the accuracy of the solution, especially for the case where a rapid change or a singularity occurs in the solution, such as occurs in non-Newtonian hydrodynamic systems [32, 33]. Therefore it is highly desirable to design such an energy stable method for the q-NSCH model, which dissipates the energy (preserve thermodynamic consistency) at the discrete level. Many time-discrete or fully discrete level energy stable methods [8, 21, 23, 28, 33, 35] have been presented for the other types of NSCH models for binary incompressible fluid with the solenoidal velocity field. However for the q-NSCH model presented by Lowengrub and Truskinovsky [36] or the other quasi-incompressible type models with the non-solenoidal velocity, relatively few time-discrete energy stable methods are available [16, 43, 19]. Very recently, a $C^{0}$ finite element method for the q-NSCH system with a consistent discrete energy law was presented by Guo et. al. [19]. where interface topological transitions are captured and the quasi-incompressibility is handled smoothly. At the fully discrete level, however, there are no available energy stable numerical methods for the q-NSCH model.

Another important criterion for the method design is to guarantee mass conservation of the binary incompressible fluid at the fully discrete level. Due to appearance of the diffusion term and numerical dissipation introduced in discretization of the convective term in the Cahn-Hilliard equation, the total mass of the binary fluid is usually not preserved exactly. This phenomenon has been reported in several works [47, 53], where the phase-field models were used to study the binary incompressible fluid. The mass loss can get even worse for binary fluid flows with large density ratios, as large numerical dissipation is needed to obtain a stable solution. To handle the issue of mass loss, usually the fine grids and small thickness of diffuse-interfaces are used in the phase-field model to improve mass conservation [3, 12]. Another way to compensate the mass loss is to add the extra mass correction terms into the CH equation [47]. However, this may need additional efforts for correcting the order parameter at each time step, and the energy stability of the numerical methods can be hardly maintained due to the artificial mass correction term.

In the present paper, we develop two fully mass conservative, energy stable, staggered grid finite difference methods for the q-NSCH model. The temperal and space discretization for both methods are so designed that the mass for the binary fluid is preserved naturally at the fully discrete level, and the extra artificial mass correction is not required. Moreover, both methods are energy stable at the fully discrete level. Our first method, primitive method, that uses a primitive variable formulation is based on the Vanka-type smoother [46], where the momentum and continuity equations are coupled implicitly and the velocities and pressures are updated simultaneously in a linear sense. Our second method that uses a projection-type formulation was originated from the work by A. Chorin [9] for solving the Navier-Stokes equations. The key advantage of the projection method is its efficiency such that the computations of the velocity and the pressure fields are decoupled. Our projection method differs from the traditional projection method in that the latter method usually uses the pressure to project the intermediate velocity onto a space of divergence-free velocity field [28], whereas we enforce the quasi-incompressible condition instead. To the best of the author’s knowledge, for the q-NSCH model, our two finite difference methods which preserve the mass and meanwhile are energy stable at the fully discrete level are new. To solve the schemes efficiently, we design a practical nonlinear multigrid solver - comprised of a standard FAS method for the CH equation, and a method based on the Vanka-type smoothing strategy for the NS equation - for solving these equations.

The rest of the paper is organized as follows: in §2, we introduce the q-NSCH model and its non-dimensionaliz-ation. For convenience of the numerical design later, we reformulate the system and show the mass conservation and the energy law (thermodynamic consistency) of the reformulated model in §3. In §4 we present our two numerical methods, namely the primitive method and projection method, at the time-discrete level, and we demonstrate the mass conservation property and energy stability of both methods. In §5, we introduce some basic definitions and notations for the finite difference discretization on a staggered grid. In §6, we descretize in space and we show that for both methods the mass conservation and energy stability can be achieved at the fully discrete level. In §7, we introduce a non-linear multigrid solver for our fully discrete numerical schemes. In §8, we test our methods and compare with the existing results. In §9 we present a convergence test for both methods. §10 is the conclusion. Moreover, the extra notations for the finite difference discretization and some useful propositions are listed in the Appendix A. The multigrid solver is briefly introduced in B.

2 Quasi-Incompressible NSCH System

2.1 Dimensional System Equations

As derived in [36], the q-NSCH system governing a binary incompressible fluid with variable density and viscosity is

[TABLE]

Here $\boldsymbol{u}$ is the velocity, $p$ is the pressure, $g$ is the gravity, $c$ is the phase variable (mass concentration), $\mu_{c}$ is the chemical potential, $f(c)=F^{\prime}(c)$, and $F(c)=c^{2}(c-1)^{2}/4$ is the double-well potential, $m(c)=\sqrt{c^{2}(1-c)^{2}}$ is the variable mobility, $\rho=\rho(c)=\rho_{1}\rho_{2}/((\rho_{2}-\rho_{1})c+\rho_{1})$ is the variable density for the binary fluid, $\rho_{1}$ and $\rho_{2}$ are the constant densities for the two incompressible fluids, $\alpha=(\rho_{2}-\rho_{1})/\rho_{1}\rho_{2}$ is a constant such that $\rho^{\prime}=-\alpha\rho^{2}$, $\mu(c)=\mu_{1}\mu_{2}/((\mu_{2}-\mu_{1})c+\mu_{1})$ is the variable viscosity, $\mu_{1}$ and $\mu_{2}$ are the constant viscosities of the two fluids, $\lambda(c)=-\frac{2}{3}\mu(c)$, $\epsilon$ is a small parameter that is related to the thickness of the diffuse-interface, $\sigma$ is the surface tension from the sharp interface model, $\eta$ is a ratio parameter that relates the phase-field model and sharp interface model [18].

The no-slip boundary condition is imposed for the velocity field

[TABLE]

and the Neumann boundary conditions are imposed for the phase-field variables,

[TABLE]

where $\boldsymbol{n}$ is the normal vector pointing out of the physical domain $\Omega$.

Note that, by multiplying (2) with $1/\alpha$ and substituting into (3), we obtain,

[TABLE]

Multiplying the above equation by $-\alpha\rho$ and using the definition of $\alpha$, we obtain

[TABLE]

which shows that the q-NSCH system (1)-(4) satisfies the mass conservation. Note that this equation is also required to reformulate the system equation in §3.

2.2 Non-dimentionalization

Let $L_{*}$ and $U_{*}$ denote the characteristic scales of length and velocity, we then introduce the dimensionless independent variables: $\hat{x}=\boldsymbol{x}/L_{*}$, $\hat{\boldsymbol{u}}=\boldsymbol{u}/U_{*}$, $\hat{t}=tU_{*}/L_{*}$, and the following natural scaling of the dependent variables: $\hat{p}=\rho_{*}{\mu_{c}}_{*}p_{l}$, $\hat{\rho}=\rho_{l}/\rho_{*}$, $\hat{\mu}=\mu_{l}/\mu_{*}$, $\hat{\mu}_{c}=\mu_{c}/{\mu_{c}}_{*}$, $\hat{\eta}=\eta_{*}/\rho_{l}$, $\hat{\epsilon}=\epsilon/L_{*}$ and $\bar{m}(c)=m_{l}/m_{*}(c)$, where the subscripts denote characteristic quantities. Omitting the hat notation, the non-dimensional q-NSCH model is

[TABLE]

where $M=U_{*}^{2}/{\mu_{c}}_{*}$ is an analogue of the squared Mach number measuring the relative strength of the kinetic energy to chemical energy [36], $We=\rho_{*}U_{*}^{2}L_{*}/\sigma$ is the Weber number, $Re=\rho_{*}L_{*}U_{*}/\mu_{l}$ is the classical Reynolds number, $Fr=U_{*}^{2}/gL_{*}$ is the Froude number, $Pe=\rho U_{*}L_{*}/m_{*}{\mu_{c}}_{*}$ is the diffusional Peclet number.

Note that the sharp-interface limit analysis is carried out in [36, 18] to show the convergence of the q-NSCH model. In particular, as the thickness of the diffuse interface approached to 0 ($\epsilon\rightarrow 0$), the q-NSCH reduces to the classical sharp-interface model for binary incompressible fluids. We will show this convergence property through the numerical simulations in §8.

3 Reformulation of the system equations

For convenience of the numerical scheme design later, we follow the same strategy used in [19] to reformulate the q-NSCH system (9)-(12) to obtain

[TABLE]

Here we have defined a new pressure $\bar{p}$ and a new chemical potential $\bar{\mu}_{c}$, such that

[TABLE]

Note that a zero term $\frac{1}{2}\rho_{t}\boldsymbol{u}+\frac{1}{2}\boldsymbol{\nabla}\cdot(\rho\boldsymbol{u})\boldsymbol{u}$ has been added to the momentum equation, which can be seen as the multiplication of the continuity equation (8) and the velocity $\boldsymbol{u}$. The reason is that the test functions in deriving the energy law can be made much simpler (See also Remark 3.3 below) and thus make it easier to achieve the energy stable finite difference method.

We now show that, in the time-continuous and space-continuous level, the reformulated non-dimensional system (13)-(16) satisfies the Mass conservation for the binary fluid and also the fluid components. In order to show these properties, we will only consider the homogeneous boundary condition for $\boldsymbol{u}_{\partial\Omega}=0$, such that, the boundary terms that are originated from the integration by parts can be dropped by using the homogeneous boundary conditions.

Theorem 3.1**.**

The non-dimensional q-NSCH system (13)-(16) preserve the mass of the binary fluid $\rho$ and the fluid components $\rho c$, ${\rm i.e}$,

[TABLE]

Proof.

Using the same scheme for (7), we obtain the continuity equation of the system

[TABLE]

Taking integration of Eq.(19) over $\Omega$ with the help of the homogeneous boundary condition of $\boldsymbol{u}$ in Eq.(5), we obtain the mass conservation for the binary fluid $\int_{\Omega}\rho_{t}{\rm d}\boldsymbol{x}=0$. Multiplying Eq.(19) by $c$, and adding to Eq.(15), we obtain

[TABLE]

Taking integration of Eq.(20) over $\Omega$ with the help of the homogeneous boundary condition (6), we obtain the mass conservation for the single fluid $\int_{\Omega}(\rho c)_{t}{\rm d}\boldsymbol{x}=0$. ∎

Theorem 3.2**.**

The non-dimensional q-NSCH system (13)-(16) is energy stable, namely the system equations satisfy the following energy dissipation law:

[TABLE]

where $E$ is the total energy of the binary fluid, and $||\cdot||_{L^{2}}$ denotes the norm of $L^{2}(\Omega)$ in Sobolev spaces.

Proof 1**.**

Multiplying Eq.(13) by $\boldsymbol{u}$, using integration-by-parts and dropping boundary terms, we obtain

[TABLE]

where we have used the continuity equation (19), the homogeneous boundary condition (5), the following identities

[TABLE]

Multiplying Eq.(14) by $\bar{p}/M$ and using integration-by-parts, we obtain

[TABLE]

Multiplying Eq.(15) by $\bar{\mu}_{c}/M$ and using integration-by-parts, we obtain

[TABLE]

Multiplying Eq.(16) by $c_{t}/M$ and using integration-by-parts, we obtain

[TABLE]

Summing up the four equations, (22) and (25)-(27), we obtain the energy dissipation law (21) of the continuous system equations.

Remark 3.3**.**

Note that the original non-dimensional q-NSCH system (9)-(12) also satisfies the energy law (21), which requires, however, much more complicated test functions to derive the energy law. Our reformulation makes the derivation easier with relatively simpler testing functions being used.

4 Time Discrete, Mass Conservative and Energy Stable Schemes

We now present two time-discrete methods, namely the primitive method (using the primitive variable formulation), and the projection method (using a projection-type formulation). Both methods are mass conservative and energy stable. In the projection method, we show that there is a pressure-Poisson equation naturally occurring in the reformulated system equation that can be used for solving the pressure. This differs from the traditional projection method in that the latter requires constructing an extra pressure-Poisson equation for solving the pressure. Here we present semi-discrete schemes that motivate the fully discrete schemes that we exhibit in later sections.

4.1 Time-Discrete Primitive Method

We first present the following numerical method in the primitive variable formulation. Let $\delta t>0$ denote the time step, and assume $\boldsymbol{u}^{n}$, $\bar{p}^{n}$, $c^{n}$, $\bar{\mu}_{c}^{n}$ are the solution at the time $t=n\delta t$, We then find the solutions at time $t=(n+1)\delta t$ are $\boldsymbol{u}^{n+1}$, $\bar{p}^{n+1}$, $c^{n+1}$, $\bar{\mu}_{c}^{n+1}$ that satisfy

[TABLE]

where $\rho^{n+\frac{1}{2}}=(\rho^{n+1}+\rho^{n})/2$, $c^{n+\frac{1}{2}}=(c^{n+1}+c^{n})/2$, $F^{n+\frac{1}{2}}(c)=(F(c^{n+1})+F(c^{n}))/2$, $(\boldsymbol{\nabla}c\cdot\boldsymbol{\nabla}c)^{n+\frac{1}{2}}=(\boldsymbol{\nabla}c^{n+1}\cdot\boldsymbol{\nabla}c^{n+1}+\boldsymbol{\nabla}c^{n}\cdot\boldsymbol{\nabla}c^{n})/2$ are the temporal average, and

[TABLE]

is an approximation to the nonlinear function $F^{\prime}(c)=f(c)=c(c-1)(c-1/2)$. Here we note the identity,

[TABLE]

[TABLE]

is an approximation of the nonlinear function $\partial\rho/\partial c=-\alpha\rho^{2}$, which satisfies the following identity:

[TABLE]

Note that the function $g(c^{n+1},c^{n})$ and $r(c^{n+1},c^{n})$ are critical for the achievement of the mass conservation and energy stability of the numerical schemes.

Using the following boundary conditions

[TABLE]

we now show that the time-discrete primitive method (28)-(31) satisfies the following properties:

Theorem 4.4**.**

The time-discrete primitive method (28)-(31) is mass conservative for the binary fluid and the fluid components, ${\rm i.e}$,

[TABLE]

Here we use $||\cdot||_{L^{2}}$ to denote the norm of $L^{2}(\Omega)$ in Sobolev spaces, and $(\cdot,\cdot)_{L^{2}}$ denotes the inner product in $L^{2}(\Omega)$.

Proof 2**.**

Multiplying (29) by $1/\alpha$ and substituting into (30), we obtain

[TABLE]

Multiplying the above by $-\alpha\rho^{n}$, we obtain the continuity equation at the time discrete level

[TABLE]

where we have used the identity (34) and (35). Integrating of Eq.(39) over $\Omega$, thanks to the boundary condition (36), we obtain the mass conservation for the binary fluid $(\rho^{n+1},1)_{L^{2}}=(\rho^{n},1)_{L^{2}}$. Multiplying (39) by $c^{n}$, and adding to Eq.(30), we obtain

[TABLE]

Integrating of Eq.(40) over $\Omega$, thanks to the definition (17) and boundary condition (36), we obtain the mass conservation for the single fluid $(\rho^{n+1}c^{n+1},1)_{L^{2}}=(\rho^{n}c^{n},1)_{L^{2}}$.

Theorem 4.5**.**

The time-discrete primitive method (28)-(31) is energy stable, ${\rm i.e}$,

[TABLE]

where $E^{n+1}$ is the total energy at the time discrete level.

Proof.

Multiplying Eq.(28) by $\delta t\boldsymbol{u}^{n+1}$, using integration-by-parts and dropping the boundary terms, we obtain

[TABLE]

where we have used the homogeneous boundary condition (36), the identities (23), (24), (39) and

[TABLE]

Note that for all the following derivations, the boundary terms originated from the integration-by-parts can be dropped by using the homogeneous boundary condition (36).

Multiplying Eq.(29) by $\delta t\bar{p}^{n+1}/M$ and using integration-by-parts, we obtain

[TABLE]

Multiplying Eq.(30) by $\delta t{\bar{\mu}_{c}}^{n+1}/M$ and using integration-by-parts, we obtain

[TABLE]

Multiplying (31) by $(c^{n+1}-c^{n})/M$ and using integration-by-parts, we obtain

[TABLE]

Summing up the four relations, (42) and (44)-(46), we obtain the energy stability (41) for time-discrete primitive method. ∎

4.2 Time-Discrete Projection Method

To design an efficient projection-type methods, we follow the projection formulation to decouple the computation of the velocity $\boldsymbol{u}$ and pressure $p$. In particular, an intermediate velocity that does not satisfy the quasi-incompressibility constraint (14) is computed at each time step, the pressure is then used to correct the intermediate velocity to get the next updated velocity that satisfies the quasi-incompressible constraint. Our projection method differs from traditional projection methods in that traditional methods usually use the pressure to project the intermediate velocity onto a space of divergence-free velocity fields (See [28] as example). Moreover, in our projection method, a pressure-Poisson equation (49) naturally occurs in the reformulated system equations, and an extra pressure-Poisson equation is not required, which, however, is usually compulsory in traditional projection methods (This can be done applying the divergence operator to Eq.(48), see [42] as a review). Our projection method for the q-NSCH system (13)-(16) is the following: given $\tilde{\boldsymbol{u}}^{n}$, $\boldsymbol{u}^{n}$, $\bar{p}^{n}$, $c^{n}$, $\bar{\mu}_{c}^{n}$, find the solution $\tilde{\boldsymbol{u}}^{n+1}$, $\boldsymbol{u}^{n+1}$, $\bar{p}^{n+1}$, $c^{n+1}$, $\bar{\mu}_{c}^{n+1}$ satisfying

[TABLE]

with the following boundary conditions

[TABLE]

Here the intermediate velocity $\tilde{\boldsymbol{u}}$ is solved first in (47), then the pressure $\bar{p}$ is solved in (49) and is used to correct $\tilde{\boldsymbol{u}}$ to obtain the velocity $\boldsymbol{u}$ that satisfies the quasi-incompressible constraint through the projection equation (48).

Theorem 4.6**.**

The time-discrete projection scheme (47)-(51) is mass conservative for the binary fluid and single fluid, ${\rm i.e}$,

[TABLE]

Proof.

Using the same strategy that used for Theorem 4.4, we obtain the mass conservation equation from our time-discrete projection method:

[TABLE]

where we have used the identity (35). Integrating of Eq.(54) over $\Omega$, thanks to the homogeneous boundary conditions (52), we obtain the mass conservation for the binary fluid $(\rho^{n+1},1)_{L^{2}}=(\rho^{n},1)_{L^{2}}$. Again, multiplying Eq.(54) by $c^{n+1}$, and adding to Eq.(50) we obtain

[TABLE]

Integrating of Eq.(55) over $\Omega$, thanks to the boundary condition (52), we obtain the mass conservation for fluid components $(\rho^{n+1}c^{n+1},1)_{L^{2}}=(\rho^{n}c^{n},1)_{L^{2}}$. ∎

Theorem 4.7**.**

The time-discrete projection scheme (47)-(51) is energy stable, ${\rm i.e}$,

[TABLE]

where $E^{n+1}$ is the total energy at time discrete level.

Remark 4.8**.**

Here we omit the details of the proof, as the derivations here are similar with the proof for the primitive method in many aspects. The primary differences are that in the projection method, we have one more projection equation (48), and the mass conservation equation (54) here is slightly different to that of the primitive methods, moreover the test functions are different as well. In particular, to show the energy stability of the projection method, we multiply Eq.(47) by $\delta t\tilde{\boldsymbol{u}}^{n+1}$, Eq.(48) by $\delta t\boldsymbol{u}^{n+1}$, Eq.(49) by $\delta t\bar{p}^{n+1}/M$, Eq.(50) by ${\delta t\bar{\mu}_{c}}^{n+1}/M$ and Eq.(51) by $(c^{n+1}-c^{n})/M$. After using the integration-by-parts with the homogeneous boundary conditions, we sum up the resulted relations to obtain the energy stability for the projection method (56).

Remark 4.9**.**

Both numerical methods are highly coupled and non-linear, however we never observe a problem with existence and uniqueness of the solution in all our extensive numerical experiments. Here we refer to [22] for some results and analysis methods about this issue.

5 Finite Difference Discretization on Staggered Grid

Before we present our fully discrete finite difference schemes, we first show some basic definitions and notations for the finite difference discretization on a staggered grid. Here we use the notation and results for cell-centered functions from [49, 10, 50]. Let $\Omega=(0,L_{x})\times(0,L_{y})$, with $L_{x}=m_{1}\cdot h$ and $L_{y}=m_{1}\cdot h$, where $m_{1}$ and $m_{2}$ are positive integers and $h>0$ is the spatial step size. For simplicity we assume that $L_{x}=L_{y}$. Consider the following four sets

[TABLE]

where $x_{i+\frac{1}{2}}=i\cdot h$ and $x_{i}=(i-\frac{1}{2})\cdot h$. Here $E_{m_{1}}$ and $E_{\overline{m}_{1}}$ are called the uniform partition of $[0,L_{x}]$ of size $m_{1}$, and its elements are called edge-centered points. The two points belonging to $E_{\overline{m}_{1}}$$E_{m_{1}}$are called ghost points. The elements of$C_{m_{1}}$and$C_{\overline{m}{1}}$are called cell-centered points. Again, the two points belonging to$C{\overline{m}{1}}\C{m_{1}}$are called ghost points. Analogously, the sets$E_{m_{2}}$and$E_{\overline{m}{2}}$contain the edge-centered points, and$C{m_{2}}$and$C_{\overline{m}{2}}$contain the cell-centered points of the interval$[0,L{y}]$.

We then define the following function spaces

[TABLE]

for cell-centered functions, vertex-centered functions, east-west edge-centered functions and north-south edge-centered functions respectively. Due to the different locations of the functions, we define several average and difference operators as follows:

[TABLE]

We also define an average operator $\boldsymbol{A}=(\begin{smallmatrix}A_{x}&0\\ 0&A_{y}\end{smallmatrix})$ and the following divergence operator:

[TABLE]

We refer the reader to A.1 and A.2 for a description of our notations for the above spaces and operators. Moreover $(\cdot,\cdot)_{2}$, and $[\cdot,\cdot]_{ew}$, $[\cdot,\cdot]_{ns}$, and $\langle\cdot,\cdot\rangle_{vc}$ denote the fully discrete inner product of the cell-centered, edge-centered and vertex-centered variables respectively which are defined in A.3.

Note that in this paper, the cell-centered functions are the phase variable $c$, chemical potential $\mu_{c}$, $\bar{\mu}_{c}$, and pressure $p$ and $\bar{p}$, the east-west edge-centered function is the x-component of the velocity, $u$ and $\tilde{u}$(for projection method), and the north-south edge-centered function is the y-component of the velocity, $v$ and $\tilde{v}$ (for the projection method).

6 Fully Discrete Mass Conservative, Energy Stable Schemes

In this section we describe and analyze our two staggered grid finite difference schemes for q-NSCH model. We show that the property of mass conservation and energy stability can be achieved at the fully discrete level for both schemes.

6.1 Fully Discrete Primitive Method

The fully-discrete scheme for the primitive method (28)-(31) is the following: Let $\delta t>0$ represent the time step, and the grid functions $c^{n},{\bar{\mu}}_{c}^{n},\bar{p}^{n}\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$, $u^{n}\in\mathcal{E}^{ew}_{m_{1}\times m_{2}}$ and $v^{n}\in\mathcal{E}^{ns}_{m_{1}\times m_{2}}$, and $\boldsymbol{u}^{n}=(u^{n},v^{n})$ be the solution at time $t=n\delta t$, find $c^{n+1},{\bar{\mu}}^{n+1}_{c},\bar{p}^{n+1}\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$, $u^{n+1}\in\mathcal{E}^{ew}_{m_{1}\times m_{2}}$, $v^{n+1}\in\mathcal{E}^{ns}_{m_{1}\times m_{2}}$, and $\boldsymbol{u}^{n+1}=(u^{n+1},v^{n+1})$ at $t=(n+1)\delta t$ such that:

[TABLE]

In Eq.(62), we let

[TABLE]

In Eq.(64) we let

[TABLE]

Note that the special discretization for the advection terms $a_{x}(\rho^{n}_{ew}u^{n+1}D_{x}c^{n})$ and $a_{y}(\rho^{n}_{ns}v^{n+1}D_{y}c^{n})$ in (69), and the surface tension terms $A_{x}(\rho^{n}_{ew}{\bar{\mu}_{c}}^{n+1})D_{x}c^{n}$ and $A_{y}(\rho^{n}_{ns}{\bar{\mu}_{c}}^{n+1})D_{y}c^{n}$ in (68), are introduced in A.4 and A.5 respectively. These discretizations are critical for deriving the fully mass conservation and energy stability of our primitive method. Moreover, $\lVert\boldsymbol{\nabla}_{D}c^{n+\frac{1}{2}}\rVert_{2}=\lVert\boldsymbol{\nabla}_{D}c^{n+1}\rVert_{2}+\lVert\boldsymbol{\nabla}_{D}c^{n}\rVert_{2}$ is the temperal average of the norm $\lVert\boldsymbol{\nabla}_{D}c\rVert_{2}$, which is described in A.7. The following expression for $m(c^{n+1})$ is used for the computations:

[TABLE]

We assume the cell-centered functions satisfy the following homogeneous Neumann boundary conditions

[TABLE]

and the velocity $\boldsymbol{u}^{n+1}=(u^{n+1},v^{n+1})$ satisfies the no-slip boundary condition

[TABLE]

A detailed description of the discrete boundary conditions is provided in A.6.

Theorem 6.10**.**

The fully discrete primitive scheme (62)-(65) is mass conservative for the binary fluid, ${\rm i.e}$,

[TABLE]

Proof.

Multiply Eq.(63) by $1/\alpha$ and substituting into Eq.(64), we obtain

[TABLE]

Multiplying the above equation by $-\delta t~{}\alpha\rho^{n}$, we obtain the continuity equation at the fully discrete level

[TABLE]

where we have used the relation (35) and the following identity

[TABLE]

Note that the definition of $\rho_{ew}$, $\rho_{ns}$ and a detailed derivation of (76) are given in A.4. Applying summation-by-parts to Eq.(75) and utilizing the homogeneous boundary conditions (72), we obtain the fully discrete mass conservation for the binary fluid:

[TABLE]

∎

Theorem 6.11**.**

The fully discrete primitive scheme (62)-(65) is energy stable at the fully discrete level, ${\rm i.e}$,

[TABLE]

where $E^{n+1}_{h}$ is the total energy of the system at the fully discrete level. Here all the norms are defined by Eqs.(161)-(163). Here $\rVert\cdot\rVert_{2}$ is the fully discrete norm that is defined in A.7.

Proof.

Multiplying Eq.(62) by $\delta t\boldsymbol{u}^{n+1}=\delta t(u^{n+1},v^{n+1})$ in x and y direction respectively, using summation-by-parts equations and dropping the boundary terms, we obtain

[TABLE]

where a special discretization for surface tension terms $A_{x}(\rho^{n}_{ew}{\bar{\mu}_{c}}^{n+1})D_{x}c^{n}$ and $A_{y}(\rho^{n}_{ns}{\bar{\mu}_{c}}^{n+1})D_{y}c^{n}$, and the corresponding derivations are introduced in A.5. The various summation-by-parts equations we have used here are Eqs.(119)-(133). Note that in all the derivations throughout this theorem, the boundary terms that originated from summation-by-parts can be eliminated by utilizing the homogeneous boundary conditions (71) and (72). We have also used the fully discrete mass conservation (77) and the following identity

[TABLE]

Multiplying Eq.(63) by $\delta t{\bar{p}}^{n+1}$ and using the summation-by-parts Eqs.(132) and (133), we obtain

[TABLE]

Multiplying Eq.(64) by $\delta t{\bar{\mu}_{c}}^{n+1}$ and using the summation-by-parts Eqs.(132) and (133), we obtain

[TABLE]

Again the special discretization for advection terms $a_{x}(\rho^{n}_{ew}D_{x}c^{n}u^{n+1}))$ and $a_{y}(\rho^{n}_{ns}D_{y}c^{n}v^{n+1})$, and the corresponding derivations are introduced in A.4. Multiplying Eq.(65) by $({c}^{n+1}-{c}^{n})/M$ and using the norm definition (157), we obtain

[TABLE]

Summing up the four relations (79), (81)-(83), we obtain the energy stability (78) of the primitive method at the fully discrete level. Here all the norms are defined by Eqs.(161)-(163). ∎

6.2 Fully Discrete Projection Method

The fully-discrete scheme for projection method (47)-(51) is the following: given $c^{n}$, $\bar{\mu}^{n}_{c}$, $\bar{p}^{n}\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$, $u^{n}\in\mathcal{E}^{ew}_{m_{1}\times m_{2}}$, $v^{n}\in\mathcal{E}^{ns}_{m_{1}\times m_{2}}$, and $\boldsymbol{u}^{n}=(u^{n},v^{n})$ at time $t=n\delta t$, find grid functions $c^{n+1},~{}\bar{\mu}^{n+1}_{c},~{}\bar{p}^{n+1}\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$, $u^{n+1}$, $\tilde{u}^{n+1}\in\mathcal{E}^{ew}_{m_{1}\times m_{2}}$, $v^{n+1}$, $\tilde{v}^{n+1}\in\mathcal{E}^{ns}_{m_{1}\times m_{2}}$, and $\boldsymbol{u}^{n+1}=(u^{n+1},v^{n+1})$ at time $t=(n+1)\delta t$:

[TABLE]

Note that the terms $\rho^{n}\boldsymbol{u}^{n}\cdot\boldsymbol{\nabla}\tilde{\boldsymbol{u}}^{n+1}$ and $\boldsymbol{\nabla}\cdot(\rho^{n}\boldsymbol{u}^{n})\boldsymbol{u}^{n+1}/2$ in Eq.(84), the term $\rho^{n+1}{\bar{\mu}_{c}}^{n+1}\boldsymbol{\nabla}c^{n+1}$ in Eq.(85), and the term $\rho^{n+1}\boldsymbol{u}^{n+1}\cdot\boldsymbol{\nabla}c^{n+1}$ in Eq.(87) are defined analogously as in Eqs.(66)-(69), where the main difference is that in the current method, the above three terms have different upper subscript representing the solution at the different time step. Again, $\lVert\boldsymbol{\nabla}_{D}c\rVert_{2}^{n+\frac{1}{2}}=\lVert\boldsymbol{\nabla}_{D}c^{n+1}\rVert_{2}+\lVert\boldsymbol{\nabla}_{D}c^{n}\rVert_{2}$ is the temperal average. Once again, the reader is referred to the Appendix and [50, 49, 10] for a description of the finite difference notation used here. We assume the cell-centered functions satisfy the following Neumann boundary conditions

[TABLE]

and the intermediate velocity $\tilde{\boldsymbol{u}}^{n+1}=(\tilde{u}^{n+1},\tilde{v}^{n+1})$ satisfies the no-slip boundary condition

[TABLE]

and the velocity $\boldsymbol{u}^{n+1}=(u^{n+1},v^{n+1})$ satisfies the following boundary condition

[TABLE]

A detailed description for the boundary condition is provided in A.6.

Theorem 6.12**.**

The scheme (84)-(88) is mass conservative for the two-phase fluid, ${\rm i.e}$,

[TABLE]

Proof.

Multiplying Eq.(86) by $1/\alpha$ and substituting into (87), we obtain

[TABLE]

Again the definition of $\rho_{ew}$, $\rho_{ns}$ and a detailed derivation of (76) are given in A.4. Multiplying the above equation by $-\alpha\rho^{n}$ and using the same treatment in Eqs.(76) and (77), we obtain the mass conservation of the binary fluids from our projection method:

[TABLE]

Using summation by parts, we obtain

[TABLE]

where we have shown that our projection method preserves the mass of the binary fluid at the fully discrete level. ∎

Theorem 6.13**.**

The fully discrete projection scheme (84)-(88) is energy stable, ${\rm i.e}$,

[TABLE]

where $E^{n+1}_{h}$ is the total energy of the system at the fully discrete level.

Remark 6.14**.**

Here we omit the details of the proof, as the derivations here are similar with the proof for the primitive method in many aspects. The primary differences is that in the projection method has one more projection equation (85), and the test functions are different. In particular, to show the energy stability of the projection method, we multiply Eq.(84) by $\delta t\tilde{\boldsymbol{u}}^{n+1}$, Eq.(85) by $\delta t\boldsymbol{u}^{n+1}$, Eq.(86) by $\delta t\bar{p}^{n+1}/M$, Eq.(87) by ${\delta t\bar{\mu}_{c}}^{n+1}/M$ and Eq.(88) by $(c^{n+1}-c^{n})/M$. After using the summation-by-parts with the homogeneous boundary conditions, we sum up the resulted relations to obtain the energy stability for the projection method (96).

7 Multigrid Solver

In this paper, we present an efficient nonlinear FAS multigrid solver for our schemes. The solver is motivated by that described for the Cahn-Hilliard-Brinkman scheme in an existing paper [10], where a finite difference method in primitive variable formulations is used. The primary difference is that in the present paper, the multigrid solver for our primitive method is designed for a much more complicated and highly non-linear problem comprised of a full Navier-Stokes equation and Cahn-Hilliard equation with variable density.

7.1 Primitive Methods

In the multigrid solver for our the primitive method (62)-(65), the smoothing operators for the Cahn-Hilliard equation and Navier-Stoke equation are decoupled. Specifically, for each grid cell $(i,j)$, we perform the following steps:

Update $c^{k+1}_{i,j}$, $\bar{\mu}^{k+1}_{c~{}i,j}$ using a non-linear Gauss-Seidel method on the CH equations (64) and (65). 2. 2.

Update the five variables $u^{k+1}_{i\pm\frac{1}{2},j}$, $v^{k+1}_{i,j\pm\frac{1}{2}}$, and $\bar{p}^{k+1}_{i,j}$ using a Vanka-type smoothing strategy [39, 45, 46, 48] on the NS equations (62)-(63) with the updated values for $c^{k+1}_{i,j}$, $\bar{\mu}^{k+1}_{c~{}i,j}$.

Here $k$ stands for the iteration step at the current time step. Note that $u$ and $v$ are edge-centered variables and this contributes to the complication of the method. Here we omit the details for the relaxation, and we refer to the [49] as a description for the Vanka-type smoother for the fluid equation in primitive variable formulation. Note that the smoother operator for CH equation and the Vanka-type smoother for NS equation is performed in the RedBlack order.

7.2 Projection Methods

For the projection method (84)-(88), we first perform a relaxation on the Cahn-Hilliard equations (87) and (88), which is the same as the primitive method; and then relax the flow equation (84) to obtain the intermediate velocity by using a Vanka-type smoother. Note that the smoother used here differs from that of the primitive method in that in the present case the pressure is not updated together with the four intermediate velocity variables. We next relax the mass conservation equation (86) to obtain the pressure, and finally we update the velocity through the projection equation (85).

Specifically, in the proposed smoother, for each grid cell $(i,j)$, we perform the following steps:

Update $c^{k+1}_{i,j}$, $\bar{\mu}^{k+1}_{c~{}i,j}$ by using a nonlinear Gauss-Seidel method on CH equations (87) and (88). 2. 2.

Update the four intermediate velocity variables $\tilde{u}^{k+1}_{i\pm\frac{1}{2},j}$, $\tilde{v}^{k+1}_{i,j\pm\frac{1}{2}}$ using a Vanka-type smoothing strategy on the fluid equation (84), with the updated $c^{k+1}_{i,j}$, $\bar{\mu}^{k+1}_{c~{}i,j}$. 3. 3.

Update the pressure $\bar{p}^{k+1}_{i,j}$ by using a nonlinear Gauss-Seidel method on the the mass conservation equation (86) with the updated $c^{k+1}_{i,j}$, $\bar{\mu}^{k+1}_{c~{}i,j}$ and $\tilde{u}^{k+1}_{i\pm\frac{1}{2},j}$, $\tilde{v}^{k+1}_{i,j\pm\frac{1}{2}}$. 4. 4.

Update the four velocity variables $u^{k+1}_{i\pm\frac{1}{2},j}$ and $v^{k+1}_{i,j\pm\frac{1}{2}}$ through the projection equation (85) with the updated $c^{k+1}_{i,j}$, $\bar{\mu}^{k+1}_{c~{}i,j}$, $\tilde{u}^{k+1}_{i\pm\frac{1}{2},j}$, $\tilde{v}^{k+1}_{i,j\pm\frac{1}{2}}$ and $\bar{p}^{k+1}_{i,j}$.

Here we omit the full details of the implementation and smoothing strategy and refer the reader to [10] for the remaining details of the solver. Moreover, for both methods, we use a standard FAS V-cycle approach that can be found in [10].

8 Numerical example

In this section we investigate the performance of our numerical schemes by solving several test problems. For the advantage of using the finite difference method, we will focus only on rectangular domains. Due to the page limit, only two figures are shown in the B to illustrate the mass conservation of our methods for Case 1 in Example 1.

8.1 Capillary Wave

The first test is the damping of a sinusoidal, capillary wave, which takes into account the surface tension, gravity, and two phase flows with variable density. In [40], an analytical solution was found in the case of the small-amplitude waves on an interface between incompressible viscous fluids in an infinite domain. To simulate this problem, we choose a computational domain, $\Omega=\{(x,y):0\leqslant x\leqslant 1,0\leqslant y\leqslant 1\}$. We assume that the equilibrium position of the interface coincides with $x$ axis, and the capillary wave-length equals to the length of the domain in $x$-direction. We further assume zero initial velocity, and that the initial profile of the interface given by

[TABLE]

with the perturbation $\tilde{y}(x)=0.5-H_{0}~{}{\rm{cos}}kx$, $k={2\pi}/{\lambda_{w}}=2\pi$ and the initial amplitude of the perturbation wave $H_{0}=0.01$. We set the gravity $g=1$ and surface tension $\sigma=1$. Moreover, the ratio parameter $\eta$ that relates the sharp interface model and phase field model is determined through the following equation [18]:

[TABLE]

To test our schemes, two cases with different density and viscosity ratios are considered. In Case 1 and 2, we choose the following values for kinematic viscosities and densities for the two fluids respectively:

[TABLE]

The other non-dimensional parameters are set as

[TABLE]

which is corresponding to the asymptotic analysis of the q-NSCH model [36]. Periodic conditions are imposed on the left and right boundary for the velocity $\boldsymbol{u}$, phase-field function $c$ and chemical potential $\mu_{c}$. At the upper and lower boundaries, we impose the no-slip boundary condition for the velocity, and no-flux boundary conditions for the phase-field functions $c$, $\mu_{c}$. The time step is set as $\delta t=10^{-3}$. For each case, we use two values of $\epsilon=0.005,~{}{\rm and}~{}0.0025$ with the corresponding grid size $[256\times 256]$ and $[512\times 512]$ respectively. Due to the sharp interface analysis [18], the numerical results of this phase-field model approaches to that of the sharp interface model as the value of $\epsilon$ decreases. Both schemes are computed and the numerical results are compared with the analytical solution. Figure 1 shows the capillary wave amplitude for Case 1 with density ratio $1:10$, where for both schemes, the numerical results all agree well with the analytical solution. Moreover, as $\epsilon$ decreases, the numerical results converge to the analytical solutions. In B, Figures 11 and 12 show the time evolution of the mass of single component $\rho c$ and binary fluids $\rho$, where it can be observed that both methods preserve the mass well. In particular, the mass can be conserved up to $10^{-10}$ with the primitive method, which performs slightly better than the projection method that preserves the mass up to $10^{-9}$. The similar results can be observed in Figure 2 for the Case 2 with density ratio $1:1000$. Moreover, it has been confirmed that the mass of the single component and the binary fluids are preserved up to $10^{-9}$ by primitive method and up to $10^{-8}$ by the projection method for Case 2. In Figure 3, we show the energy dissipation for both cases with both methods. As predicted by the Theorem 6.11 and 6.13, the energy decreases for both methods, exhibiting a similar dissipation way.

8.2 Rising Droplets

As a second test, we simulate the dynamics of rising droplets. The test setup is taken from in [37]. In particular, the computational domain $\Omega=[0,1]\times[0,2]$ is filled with the heavier fluid ($c=0$) and a initially circular shaped lighter fluid ($c=1$) is placed inside. The initial drop has a radius of 0.25 and is centered at $[0.5,0.5]$. This leads to the initial profile of the interface given by

[TABLE]

with $r=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}$ and $R_{0}=0.25$. The parameters of the outer fluid are $\rho_{2}=1000$, $\mu_{2}=10$, and $\rho_{1}=100$, $\mu_{1}=1$ for the drop fluid. The gravity is $\boldsymbol{g}=(0,0.98)$, and the surface tension is $\sigma=24.5$, which lead to the values of following non-dimensional parameters

[TABLE]

As in the previous example, periodic conditions are imposed on the left and right boundary for the velocity $\boldsymbol{u}$, pressure $\bar{p}$, phase-field function $c$ and chemical potential $\mu_{c}$, and the no-slip boundary condition for the velocity, and no-flux boundary conditions for the phase-field functions $c$, $\mu_{c}$ are imposed at the upper and lower boundaries. The time step is $\delta t=10^{-3}$ for the two methods. Moreover, to test the convergence of the diffuse-interface, we use two values of $\epsilon=0.005~{}{\rm and}~{}0.0025$, which corresponds to the grid size $[256\times 256]$ and $[512\times 512]$ respectively.

Because the droplet is lighter than the surrounding fluid, the droplet rises. For a rigorous estimate of the accuracy of the simulation, we calculated the rising velocity that is determined by:

[TABLE]

where $v$ is the second (vertical) component of the velocity $\boldsymbol{u}$. Moreover, to show the quasi-incompressibility of the q-NSCH model, we calculate the divergence of velocity $\boldsymbol{\nabla}_{d}\cdot\boldsymbol{u}$ at the fully discrete level.

Snapshots of the deformed droplet interfaces and the $\boldsymbol{\nabla}\cdot\boldsymbol{u}$ (quasi-incompressibility) are presented in Figure 4, where we observe that the drop deforms slowly, resulting in a mushroom shape. Recall that the divergence-free condition does not hold for quasi-incompressible fluids with different densities because the fluids may mix slightly across the interface. The two incompressible fluids can be compressible across the interface where the two components are mixed. It can be observed that the fluid is incompressible ($\boldsymbol{\nabla}\cdot\boldsymbol{u}=0$) almost everywhere except along the moving interface. Near the interface, waves of expansion ($\boldsymbol{\nabla}\cdot\boldsymbol{u}>0$ ) and compression ($\boldsymbol{\nabla}\cdot\boldsymbol{u}<0$ ) are observed. Figure 5 shows the droplet shapes at the final time ($t=3$), where we observe that the droplet shapes differ clearly for different values of $\epsilon$ but seem to converge so that there is no big difference for the finest values $\epsilon=0.0025$ and the result obtained from the [37] by using a sharp interface model. Figure 6 plots the rising velocity of our numerical relusts and the result obtained from [37], where the agreement improves as $\epsilon\rightarrow 0$. In Figure 7, we show the energy dissipation of the binary fluid system obtained from both methods by using different values of $\epsilon$. Note that, to show the energy dissipation converges as $\epsilon$ is decreased, we also compute the example by using a even smaller $\epsilon=0.00125$. As expected, the energy decreases and yields very similar way for both methods. it can be observed that the It has been confirmed that the mass of the single component and the binary fluids are preserved up to $10^{-10}$ by primitive method and up to $10^{-9}$ by the projection method.

8.3 Rayleigh-Taylor Instability

Our last test is the Rayleigh-Taylor instability which would occur for any perturbation along the interface between a heavy fluid ($c=0$) on top of a light fluid ($c=1$), and is characterised by the density difference between the two fluids. The instability is characterized by the non-dimensional parameter Atwood ratio, that $At=(\rho_{A}-\rho_{B})/(\rho_{A}+\rho_{B})$. The initial growth and long-time evolution of Rayleigh-Taylor instability was investigated by Tryggvason [44] for inviscid incompressible flows with zero surface tension, at $At=0.50$. Guermond et al [17] studied this instability at the same value of At but accounted for viscous effects. Ding [11] studied this instability problem by using the a different phase-field model where the velocity satisfies the divergence free constraint. We validate our code here by investigating the same problem as Guermond et al [17], i.e., at $At=0.50$ and $Re(=\rho_{l}L^{3/2}g^{1/2}/\mu_{l})=3000$, with the initial interface being located in a rectangular domain $[0,L]\times[0,4L]$ at $\tilde{y}(x)=2L+0.1L\rm{cos}(2\pi x/L)$, which represents a planar interface superimposed by a perturbation of wave number $k=1$ and amplitude $0.1L$. Here we set $L=1$, and we take

[TABLE]

Here we set the gravitational acceleration $g=1$ and surface tension $\sigma=0$. In the present case of zero surface tension, the Cahn-Hilliard equation simply amounts to interface tracking only. The non-dimensional parameters are set as

[TABLE]

Periodic conditions are imposed on the left and right boundary for $\boldsymbol{u}$, $c$ and $\mu_{c}$. At the upper and lower boundaries, we impose the no-slip boundary condition for $\boldsymbol{u}$, and no-flux boundary conditions for $c$ and $\mu_{c}$. We use two values of $\epsilon=0.005,~{}{\rm and}~{}0.0025$ which corresponds to the grid size $[128\times 512]$ and $[256\times 1024]$. We set time-step $\delta t=10^{-3}$. Results are presented in Figure 8 in terms of the y-coordinate of the top of the rising fluid and the bottom of the falling fluid, together with the corresponding previous results of Tryggvason [44], Guermond et al [17] (sharp interface models) and Ding [11] (another phase-field model). For both methods, good agreement is observed with these results. As the value of $\epsilon$ decreases, our numerical resutls converge to the other numerical results. The evolution of the interface of our numerical results (projection with $\epsilon=0.0025$) and the results in [11] are shown in Figure 9, in which the rolling-up of the falling fluid can be clearly seen. At the early time, two counter-rotating vortices are formed along the sides of the falling filament and grow with time. To a certain extent, the two vortices are shed and a pair of secondary vortices occurs at the tails of the roll-ups. Our results agree with those obtained in [11]. However, comparing to their results, the small structures around the vortices are preserved and can be observed more clearly due to the mass conservative property of our numerical methods. The time evolution of the energy is shown in Figure 10, which decreases as expected. It has been confirmed that the mass of the single component and the binary fluids are preserved up to $10^{-12}$ by primitive method and up to $10^{-10}$ by the projection method.

9 Convergence test

To show our two methods are both first order accurate in time and second order accurate in space, we carry out a convergence test by considering the Cauchy sequence. We compute our q-NSCH system with the following function as the initial condition for the phase variable $c$,

[TABLE]

and zero for all the other variables, including $\mu$, $p$ and $\boldsymbol{u}$. The computational domain is $[0,1]\times[0,1]$, and the homogeneous Neumann boundary conditions for $c$ $\mu$ and $p$, and $\boldsymbol{u}$ are applied on the boundary. We refine the mesh and time step according to the schedule

[TABLE]

Here $m$ is the grid points in both $x$ direction and $y$ direction, $\Delta t$ is the time step. To compare solutions on different grid resultions, we push the solution at the coarse grid up to the next fine grid to calculate the difference in $L_{2}$ norm, and then obtain the convergence rate which are shown in Table 1. The second order convergence rate are achieved for both methods, which indicates that both of our methods are second order accurate in space and first order accurate in time. Moreover, we also present the average computational cost for each time step of both methods with different grid resolutions in Table 2, which indicates that the Projection method is much more efficient than the Primitive method. All the tests are carried out on a desktop with 4.0GHz AMD(R) FX(TM)-8350 processor.

10 Conclusion

In this paper we presented and analyzed two fully mass conservative and energy stabel finite difference schemes for the q-NSCH sytem governing the binary incompressible fluid flows with variable density and viscosity. At the continuous level, we reformulated the system equations and show that the system conserves both the component and binary fluid mass and the energy is non-increasing due to the energy dissipation law that underlies the model. Based on the reformulated system, we introduced two time-discrete and fully discrete numerical methods using primitive variable and projection-type formulation. Both schemes are fully mass conserving and the extra mass correction is not necessary. Moreover, the fully discrete energy stability are achieved for both methods, where the enregy functionals are always non-increasing. In particular, our projection method differs from the traditional projection methods in two ways:

(a) due to the quasi-incompressibility of the q-NSCH model, the pressure here is used to correct the intermediate velocity to get the updated velocity that satisfies the quasi-incompressible constraint, whereas in most of the existing projection methods the intermediate velocity is projected onto a space of divergence-free velocity field due to the corresponding model constraints (See [28] as example);

(b) in our projection method, a pressure-Poisson equation (continuity equation (49)) naturally occurs in the reformulated system equations, whereas, in most of the projection methods, the extra construction of the pressure-Poisson equation is required (this is usually done by applying the divergence operator to the decoupled equation (48)); We also present an efficient nonlinear FAS multigrid method for each method, which is motivated by a Vanka-type smoothing strategy for the Cahn-Hilliard-Brinkman equation [10].

Three numerical examples are investigated numerically, including the Capillary Wave, Rayleigh-Taylor instability and rising droplets. Quantitative comparisons are made with the existing analytical solution or the previous numerical results to validate the accuracy of our numerical schemes. Moreover it has been confirmed that the mass of the single fluid and the binary fluids are preserved up to $10^{-8}$, and the energy is always non-increasing for all the examples. The convergence property of the q-NSCH is investigated as well. In particular, our numerical results converge to the analytical or numerical solutions of the other sharp interface models as the thickness of the diffuse-interface decreases. Moreover, in the rising droplet example, we show that the quasi-incompressibility ($\boldsymbol{\nabla}\cdot\boldsymbol{u}\neq 0$) is captured smoothly together with the move interface which indicates that the quasi-incompressibility does not give any problems to our schemes.

Acknowledgement

ZG gratefully acknowledge partial support from the $150^{\rm th}$ Anniversary Postdoctoral Mobility Grant of London Mathematical Society (PMG14-15 09). JL and ZG gratefully acknowledge partial support from National Science Foundation Grants NSF-DMS-1719960, NSF-DMS-1522775 and the National Institute of Health grant P50GM76516 for the Center of Excellence in Systems Biology at the University of California, Irvine. PL is partially supported by the National Natural Science Foundation of China (No. 91430106) and the Fundamental Research Funds for the Central Universities (No. 06500073). SMW gratefully acknowledges support from a grant from the National Science Foundation (NSF-DMS 1418692) and partial support from the Mathematics Department at the University of California, Irvine. Part of this work was completed during a sabbatical visit (03/2015-06/2015) spent at that institution.

Appendix A Finite Difference Discretization on a Staggered Grid

A.1 Basic Definitions and Properties

Here we use the notations and results for cell-centered functions from [49, 10, 50]. The reader is directed there for complete details. We begin with definitions of grid functions and difference operators needed for our discretization of a two-dimensional staggered grid. Throughout this appendix, we use the following symbols to denote the cell-centered, edge-centered and vertex-centered functions, such that

[TABLE]

Here we use the function spaces with over-lined subscript $\overline{m}_{1}$ and $\overline{m}_{2}$ to denote the space that include the ghost points in the $x$ and $y$ direction respectively. In component form, we define

[TABLE]

where $x_{i}=(i-\frac{1}{2})\cdot h$, $y_{j}=(j-\frac{1}{2})\cdot h$, and $i$ and $j$ can take on integer values.

A.2 Average and Difference Operators

We define the edge-to-center average and difference operators $a_{x},~{}d_{x}:\mathcal{E}^{ew}_{m_{1}\times m_{2}}\rightarrow\mathcal{C}_{m_{1}\times m_{2}}$ and $a_{y},~{}d_{y}:\mathcal{E}^{ns}_{m_{1}\times m_{2}}\rightarrow\mathcal{C}_{m_{1}\times m_{2}}$ component-wise via

[TABLE]

for $i=1,\cdots,m_{1}$ and $j=1,\cdots,m_{2}$.

The center-to-edge average and difference operators, $A_{x}$, $D_{x}$ : $\mathcal{C}_{\overline{m}_{1}\times m_{2}}\rightarrow\mathcal{E}^{ew}_{m_{1}\times m_{2}}$ and $A_{y}$, $D_{y}$ : $\mathcal{C}_{m_{1}\times\overline{m}_{2}}\rightarrow\mathcal{E}^{ns}_{m_{1}\times m_{2}}$ are defined component-wise as

[TABLE]

for $i=0,\cdots,m_{1}$ and $j=0,\cdots,m_{2}$.

The center-to-vertex average operator $\mathcal{A}$ : $\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}\rightarrow\mathcal{V}_{m_{1}\times m_{2}}$ is defined as

[TABLE]

for $i=0,\cdots,m_{1}$ and $j=0,\cdots,m_{2}$.

The edge-to-vertex average and difference operators, $\mathcal{A}_{x}$, $\mathcal{D}_{x}$ : $\mathcal{E}^{ns}_{\overline{m}_{1}\times m_{2}}\rightarrow\mathcal{V}_{m_{1}\times m_{2}}$, and $\mathcal{A}_{y}$, $\mathcal{D}_{y}$ : $\mathcal{E}^{ew}_{m_{1}\times\overline{m}_{2}}\rightarrow\mathcal{V}_{m_{1}\times m_{2}}$ are defined as

[TABLE]

for $i=0,\cdots,m_{1}$ and $j=0,\cdots,m_{2}$. The vertex-to-edge average and difference operators, $\mathfrak{A}_{x}$, $\mathfrak{D}_{x}$ : $\mathcal{V}_{m_{1}\times m_{2}}\rightarrow\mathcal{E}^{ns}_{m_{1}\times m_{2}}$, and $\mathfrak{A}_{y}$, $\mathfrak{D}_{y}$ : $\mathcal{V}_{m_{1}\times m_{2}}\rightarrow\mathcal{E}^{ew}_{m_{1}\times m_{2}}$ are defined as

[TABLE]

for $i=0,\cdots,m_{1}$ and $j=1,\cdots,m_{2}$.

A.3 Weighted Inner-Products

Based on the above definitions, we define the following 2D weighted grid inner-products:

[TABLE]

where $\phi\in\mathcal{C}_{m_{1}\times m_{2}}$, $u,\gamma\in\mathcal{E}^{ew}_{m_{1}\times m_{2}}$, $v,\omega\in\mathcal{E}^{ns}_{m_{1}\times m_{2}}$, and $f,g\in\mathcal{V}_{m_{1}\times m_{2}}$. We also define the following combined 2D weighted grid inner-products:

[TABLE]

We also define the one-dimensional inner-products for the edge-centered functions, or cell-centered functions, or the multiple combination of the edge-centered and cell-centered functions. Here for simplicity, we only introduce the one-dimensional inner-product for the edge-centered functions, the others combinations can be defined analogously:

[TABLE]

Here the first is defined for $u,\gamma\in\mathcal{E}^{ew}_{m_{1}\times m_{2}}$, and the second for $v,\omega\in\mathcal{E}^{ns}_{m_{1}\times m_{2}}$. $*$ indicates the sum of the functions in the direction along which the one-dimensional inner-product acts. Note that, throughout this section, all the boundary terms that originated from the summation-by-parts can be eliminated by using the homogeneous Neumann conditions for the cell-centered variables, and the no-slip boundary conditions for the edge centered variables. The results are also valid for the periodic boundary conditions for the cell-centered or edge-centered variables case. For the above definitions, we obtain the following results:

Proposition A.15**.**

(Summation-by-parts) if $\phi\in\mathcal{C}_{\overline{m}_{1}\times m_{2}}$, $u\in\mathcal{E}^{ew}_{m_{1}\times m_{2}}$ and $v\in\mathcal{E}^{ns}_{m_{1}\times m_{2}}$ then

[TABLE]

Proposition A.16**.**

Let $\phi\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$, $u,\gamma\in\mathcal{E}^{ew}_{\overline{m}_{1}\times m_{2}}$ and $v,\omega\in\mathcal{E}^{ns}_{m_{1}\times\overline{m}_{2}}$. Then

[TABLE]

Proposition A.17**.**

Let $\phi\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$, $u\in\mathcal{E}^{ew}_{\overline{m}_{1}\times\overline{m}_{2}}$ and $v\in\mathcal{E}^{ns}_{\overline{m}_{1}\times\overline{m}_{2}}$. Then

[TABLE]

For the sake of simplicity, we omit all the boundary terms of that originated in the following Propositions.

Proposition A.18**.**

Let $\phi\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$, $u,\gamma\in\mathcal{E}^{ew}_{\overline{m}_{1}\times\overline{m}_{2}}$ and $v,\omega\in\mathcal{E}^{ns}_{\overline{m}_{1}\times\overline{m}_{2}}$. Then

[TABLE]

Proposition A.19**.**

Let $\phi\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$, $u,\gamma\in\mathcal{E}^{ew}_{\overline{m}_{1}\times\overline{m}_{2}}$ and $v,\omega\in\mathcal{E}^{ns}_{\overline{m}_{1}\times\overline{m}_{2}}$. Then

[TABLE]

Proposition A.20**.**

Let $\phi\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$, $u\in\mathcal{E}^{ew}_{\overline{m}_{1}\times\overline{m}_{2}}$ and $v\in\mathcal{E}^{ns}_{\overline{m}_{1}\times\overline{m}_{2}}$. Then

[TABLE]

Proposition A.21**.**

Let $\phi,\psi,\zeta\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$. Then

[TABLE]

A.4 Special Treatment for the Advection Term in Cahn-Hilliard equation

To let our numerical schemes satisfy the mass conservation together with the energy dissipation law at the fully discrete level, we have employed a special treatment for the advection term in the Cahn-Hilliard equation (64) in primitive method and (87) in projection method, such that:

[TABLE]

for $i=1,\dots,m_{1}$ and $j=1,\dots,m_{2}$. Note that, for the sake of simplicity, we omit the upper subscript $n$ or $n+1$ that represent the time step for all the variables. These two terms are so designed that the primitive method (62)-(65) and the projection-type method (84)-(88) can satisfy the mass conservation property and the energy stability at the fully discrete level. In particular, to prove the identities (76) and (94), we show that

[TABLE]

Here we have used the definition of the variable density $\rho$, and the following identity which is analog to (35) in the temporal discretization:

[TABLE]

A.5 Special Treatment for the Surface Tension

To let our numerical schemes satisfy the mass conservation together with the energy dissipation law at the fully discrete level, we have employed a special treatment for the surface tension term in the momentum equations (62) of our primitive method, and in the projection equations (85) of our projection-type method, such that

[TABLE]

for $i=1,\dots,m_{1}$ and $j=1,\dots,m_{2}$. Note that, for the sake of simplicity, we omit the upper subscript $n$ or $n+1$ that represent the time step for all the variables. These two terms are so designed that our two numerical schemes can satisfy the energy dissipation law at the fully discrete level. In particular, when we multiply the momentum equation (62) or (85) by $\boldsymbol{u}=(u,v)$, the surface tension term can be written as

[TABLE]

Note that, throughout this section, all the boundary terms that originated from the summation-by-parts can be eliminated by using the homogeneous Neumann conditions for the cell-centered variables, and the no-slip boundary conditions for the edge centered variables. The results are also valid for the periodic boundary conditions for the cell-centered or edge-centered variables case.

A.6 Boundary Conditions

A.6.1 Cell-Centered Functions

In this paper we use the Neumann or periodic boundary condition for the cell-centered functions. Specifically, we shall say the cell-centered function $\phi\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$ satisfies homogeneous Neumann boundary conditions if and only if

[TABLE]

We use the notation $\boldsymbol{n}\cdot\boldsymbol{\nabla}_{D}\phi=0$ to indicate that $\phi$ satisfies (144) and (145). The cell-centered function $\phi\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$ satisfies periodic boundary conditions if and only if

[TABLE]

A.6.2 Edge-Centered Functions

We use the no-slip or periodic boundary conditions for the edge-centered functions. We shall say the velocity $\boldsymbol{u}=(u,v)$ (for $u\in\mathcal{E}^{ew}_{m_{1}\times m_{2}}$, $v\in\mathcal{E}^{ns}_{m_{1}\times m_{2}}$) satisfies the no-slip boundary conditions $\boldsymbol{u}|_{\Omega}=0$ if and only if

[TABLE]

We shall say the the velocity $\boldsymbol{u}=(u,v)$ (for $u\in\mathcal{E}^{ew}_{m_{1}\times m_{2}}$, $v\in\mathcal{E}^{ns}_{m_{1}\times m_{2}}$) satisfies the boundary condition $\boldsymbol{n}\cdot\boldsymbol{u}|_{\Omega}=0$ if and only if

[TABLE]

We shall say the edge-centered function $u\in\mathcal{E}^{ew}_{\overline{m}_{1}\times m_{2}}$ ($v\in\mathcal{E}^{ns}_{m_{1}\times\overline{m}_{2}}$) satisfies the periodic boundary conditions on the east-west (north-south) boundaries if and only if

[TABLE]

A.7 Norms

We define the following norms for cell-centered functions. If $\phi\in\mathcal{C}_{m_{1}\times m_{2}}$, then $\lVert\phi\rVert_{2}:=\sqrt{h^{2}(\phi,\phi)_{2}}$. For $\phi\in\mathcal{C}_{\overline{m}_{1}\times\overline{m}_{2}}$, we define the following norms

[TABLE]

We also define the following norms for the vector of the edge-centered functions $\boldsymbol{u}=(u,v)$ together with the cell-centered function $\phi$ , where $u\in\mathcal{E}^{ew}_{m_{1}\times m_{2}}$, $v$ $\in\mathcal{E}^{ns}_{m_{1}\times m_{2}}$ and $\phi$ $\in\mathcal{C}_{m_{1}\times m_{2}}$, such that

[TABLE]

Appendix B Figures for §8.1