An energy stable $C^0$ finite element scheme for a quasi-incompressible phase-field model of moving contact line with variable density

TL;DR

This paper introduces an energy stable C0 finite element scheme for simulating two-phase flows with moving contact lines and variable density, ensuring mass conservation and energy stability.

Contribution

It develops a thermodynamically consistent phase-field model with a novel finite element scheme that is mass conserving and energy stable for quasi-incompressible flows.

Findings

Scheme is mass conservative and energy stable

Achieves 2nd order convergence with P1 elements

Achieves 3rd order convergence with P2 elements

Abstract

In this paper, we focus on modeling and simulation of two-phase flow with moving contact lines and variable density. A thermodynamically consistent phase-field model with General Navier Boundary Condition is developed based on the concept of quasi-incompressibility and the energy variational method. Then a mass conserving and energy stable C0 finite element scheme is developed to solve the PDE system. Various numerical simulation results show that the proposed schemes are mass conservative, energy stable and the 2nd order for P1 element and 3rd order for P2 element convergence rate in the sense of L2 norm.

Tables

Table 1. Table 1: L 2 superscript 𝐿 2 L^{2} norm of the error and convergence rate for velocity 𝐮 = ( u x , u y ) 𝐮 subscript 𝑢 𝑥 subscript 𝑢 𝑦 \mathbf{u}=(u_{x},u_{y}) , phase function c 𝑐 c , at time t = 0.2 𝑡 0.2 t=0.2 with density ratio ρ 1 : ρ 2 = 0.8 : 1 : subscript 𝜌 1 subscript 𝜌 2 0.8 : 1 \rho_{1}:\rho_{2}=0.8:1 viscosity ratio η 1 : η 2 = 1 : 1 : subscript 𝜂 1 subscript 𝜂 2 1 : 1 \eta_{1}:\eta_{2}=1:1 .

Space step $h$	P1 Element
	Err($u_x$)	Rate	Err($u_y$)	Rate	Err(c)	Rate
1/160	4.9e-2		2.6e-2		8.5e-2
1/226	1.9e-2	2.57	1.8e-2	1.46	5.3e-2	1.59
1/320	7.2e-3	2.64	6.4e-3	2.78	2.2e-2	2.38
Space step $h$	P2 Element
	Err($u_x$)	Rate	Err($u_y$)	Rate	Err(c)	Rate
1/80	2.9e-2		2.5e-2		7.4e-2
1/113	9.9e-3	2.93	7.5e-3	3.32	2.8e-2	2.63
1/160	3.0e-3	3.34	2.5e-3	2.97	1.0e-2	2.80

Table 2. Table 2: L 2 superscript 𝐿 2 L^{2} norm of the error and convergence rate for velocity 𝐮 = ( u x , u y ) 𝐮 subscript 𝑢 𝑥 subscript 𝑢 𝑦 \mathbf{u}=(u_{x},u_{y}) , phase function c 𝑐 c , at time T = 0.2 𝑇 0.2 T=0.2 with density ratio ρ 1 : ρ 2 = 0.1 : 10 : subscript 𝜌 1 subscript 𝜌 2 0.1 : 10 \rho_{1}:\rho_{2}=0.1:10 viscosity ratio η 1 : η 2 = 0.1 : 10 : subscript 𝜂 1 subscript 𝜂 2 0.1 : 10 \eta_{1}:\eta_{2}=0.1:10 .

Space step $h$	P1 Element
	Err($u_x$)	Rate	Err($u_y$)	Rate	Err(c)	Rate
1/160	8.8e-3		4.9e-3		2.2e-2
1/226	6.9e-3	1.27	3.8e-3	1.27	2.0e-2	1.12
1/320	2.2e-3	3.17	1.5e-3	2.59	7.8e-3	2.52
Space step $h$	P2 Element
	Err($u_x$)	Rate	Err($u_y$)	Rate	Err(c)	Rate
1/80	5.8e-3		3.6e-3		2.5e-2
1/113	4.0e-3	1.46	1.9e-3	1.86	1.1e-2	2.20
1/160	1.1e-3	3.83	6.3e-4	3.08	3.5e-3	3.27

Equations

$$\rho =\frac{M}{V}=\frac{M_1}{V}+\frac{M_2}{V}={\bar{\rho }}_1+{\bar{\rho }}_2.$$

$$\frac{1}{\rho }=\frac{V}{M}=\frac{V_1}{M}+\frac{V_2}{M}=\frac{c_1}{{\rho }_1}+\frac{c_2}{{\rho }_2}=\frac{c}{{\rho }_1}+\frac{1-c}{{\rho }_2},$$

$$\frac{\partial {\bar{\rho }}_i}{\partial t}+\nabla \cdot ({\mathbf{u}}_i{\bar{\rho }}_i)=0,i=1,2.$$

$$\rho \mathbf{u}={\bar{\rho }}_1{\mathbf{u}}_1+{\bar{\rho }}_2{\mathbf{u}}_2,$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\mathbf{u}\rho )=0.$$

$$\frac{d}{dt}{\int }_{V(t)}\rho cdx=-{\int }_{\partial V(t)}{\mathbf{j}}_c\cdot \mathbf{n}dS,$$

$$\frac{d}{dt}{\int }_{V(t)}\rho \mathbf{u}dx={\int }_{\partial V(t)}({\mathbf{\sigma }}_{\eta }+{\mathbf{\sigma }}_c)\cdot \mathbf{n}dS.$$

$$\displaystyle\left\{\begin{array}[]{l}\rho\frac{Dc}{Dt}=-\nabla\cdot\bm{j}_{c},\\ \frac{D\rho}{dt}+\rho\nabla\cdot\mathbf{u}=\frac{\partial\rho}{\partial t}+\nabla\cdot(\mathbf{u}\rho)=0,\\ \rho\frac{D\mathbf{u}}{Dt}=\nabla\cdot\bm{\sigma}_{\eta}+\nabla\cdot\bm{\sigma}_{c},\end{array}\right.$$

$$\nabla \cdot \mathbf{u}=-\frac{1}{\rho }\frac{D\rho }{Dt}=\frac{1}{{\rho }^2}\frac{d\rho }{dc}\left(\nabla \cdot {\mathbf{j}}_c\right).$$

$$\alpha =-\frac{1}{{\rho }^2}\frac{d\rho }{dc}=\frac{{\rho }_2-{\rho }_1}{{\rho }_1{\rho }_2},$$

$$\nabla \cdot \mathbf{u}=-\alpha \nabla \cdot {\mathbf{j}}_c.$$

$$\displaystyle\left\{\begin{array}[]{l}\mathbf{u}\cdot\mathbf{n}=0,~{}~{}\mathbf{u}^{s}\cdot\bm{\tau}_{i}=u^{s}_{\tau_{i}}=f_{\tau_{i}}\\ \frac{D_{\Gamma}c}{Dt}=J_{\Gamma},\\ \bm{j}_{c}\cdot\mathbf{n}=0,\\ \end{array}\right.$$

$$E^{tot}$$

$$G(c)=\frac{1}{4}{c}^2(1-c{)}^2,f_w(c)=-\frac{\sigma }{2}\cos ({\theta }_s)\sin (\frac{(2c-1)\pi }{2}),$$

$$\Delta ={\int }_{\Omega }2\eta (c)|{\mathbf{D}}_{\eta }{|}^{2}d\mathbf{x}+{\int }_{\Omega }\lambda (c)|\nabla \cdot \mathbf{u}{|}^{2}d\mathbf{x}+{\int }_{\Omega }\frac{1}{\mathcal{M}}|{\mathbf{j}}_c{|}^{2}d\mathbf{x}+{\int }_{\partial {\Omega }_w}\left(\frac{1}{{\mathcal{M}}_{\Gamma }}{J}_{\Gamma }^2+{\beta }_{\Gamma }|{\mathbf{u}}^s{|}^2\right)dS,$$

$$\frac{1}{\eta (c)}=\frac{c}{{\eta }_1}+\frac{(1-c)}{{\eta }_2},\frac{1}{{\beta }_{\Gamma }(c)}=\frac{c}{{\beta }_{{\Gamma }_1}}+\frac{(1-c)}{{\beta }_{{\Gamma }_2}},$$

$$\frac{d}{dt}{\int }_{\Omega }\rho (\mathbf{x},t)f(\mathbf{x},t)d\mathbf{x}={\int }_{\Omega }\rho \frac{Df}{Dt}d\mathbf{x}.$$

$$\frac{d{E}^{tot}}{dt}=\frac{d}{dt}{E}_{kin}+\frac{d}{dt}{E}_{mix}+\frac{d}{dt}{E}_w=I_1+I_2+I_3.$$

$$I_1$$

$$I_2$$

$$I_3=\frac{d}{dt}{\int }_{\partial {\Omega }_w}{f}_{w}dS={\int }_{\partial {\Omega }_w}\frac{d{f}_w}{dc}\frac{\partial c}{\partial t}dS.$$

$$\frac{d{E}^{tot}}{dt}$$

$$\tilde{\mu }$$

$$L(c)$$

$$\displaystyle\left\{\begin{array}[]{l}\bm{j}_{c}=-\mathcal{M}\nabla\tilde{\mu},\\ \bm{\sigma}_{\eta}=2\eta\bm{D}_{\eta}+\lambda\nabla\cdot\mathbf{u}I-pI=\eta(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})+\lambda\nabla\cdot\mathbf{u}I-pI,\\ \bm{\sigma}_{c}=-\lambda_{c}\gamma^{2}\rho\left(\nabla c\otimes\nabla c\right),\\ J_{\Gamma}=-\mathcal{M}_{\Gamma}L(c),\\ u^{s}_{\tau_{i}}=\beta_{\Gamma}^{-1}\bm{\tau}_{i}\cdot(-(\bm{\sigma}_{\eta}+\bm{\sigma}_{c})\cdot\mathbf{n}+\frac{df_{w}}{dc}\nabla_{\Gamma}c).\end{array}\right.$$

$$u_{{\tau }_i}^s={\beta }_{\Gamma }^{-1}{\mathbf{\tau }}_i\cdot (-{\mathbf{\sigma }}_{\eta }\cdot \mathbf{n}+L(c){\nabla }_{\Gamma }c).$$

$$\rho \frac{Dc}{Dt}=\nabla \cdot (\mathcal{M}\nabla \tilde{\mu }),$$

$$\tilde{\mu }={\lambda }_c\frac{dG}{dc}-\frac{{\lambda }_c{\gamma }^2}{\rho }\nabla \cdot (\rho \nabla c)+\alpha p,$$

$$\frac{D\rho }{Dt}+\rho \nabla \cdot \mathbf{u}=0,$$

$$\rho \frac{D\mathbf{u}}{Dt}=\nabla \cdot (2\eta {\mathbf{D}}_{\eta })+\nabla (\lambda \nabla \cdot \mathbf{u})-\nabla p-\nabla \cdot ({\lambda }_c{\gamma }^2\rho \nabla c\otimes \nabla c),$$

$$\displaystyle\left\{\begin{array}[]{l}\frac{D_{\Gamma}c}{Dt}=-\mathcal{M}_{\Gamma}L(c),\\ \nabla\tilde{\mu}\cdot\mathbf{n}=0,\\ \mathbf{u}\cdot\mathbf{n}=0,\\ u_{\tau_{i}}^{s}=\beta_{\Gamma}^{-1}(-(\mathbf{n}\cdot\bm{\sigma}_{\eta}\cdot\bm{\tau}_{i})+L(c)\partial_{\tau_{i}}c),i=1,2.\end{array}\right.$$

$$\frac{d{E}^{tot}}{dt}$$

$${\int }_{\Omega }\rho \frac{Dc}{Dt}\tilde{\mu }d\mathbf{x}=-{\int }_{\Omega }\mathcal{M}|\nabla \tilde{\mu }{|}^{2}d\mathbf{x},$$

Full text

1 Introduction

The modeling and simulation of moving contact lines (MCLs), where the interface of two or more immersible fluids intersects with a solid wall [1, 2], have attracted much attention in recent years. Applications of MCLs in industries and medical fields (for example, printing [3], spray cooling of surfaces [4], blood clot [5], microfluidics [6], surfactant [7, 8]) have motivated scientific interests and mathematical challenges on associated issues such as the stress singularity and contact angle hysteresis. In order to model the dynamics around the contact lines, various types of models and approximations have been developed, such as direct molecular dynamics simulations [9, 10, 11], phase-filed models [12, 13, 14, 15, 16, 17], microscopic–macroscopic hybrid model [18, 19], front tracking model [20, 21, 22, 8] and Lattice Boltzmann model [23]. For reviews of the current status of the MCLs problem, we refer to the articles [24] and [25].

Among those models, phase-field method (or diffusive interface method) [26, 27, 28] is one of the most popular and powerful methodologies. It has two main advantages. Firstly, it is easy to track the interface and numerically implement even if there are topological changes [29]. Secondly, it can be derived by energy-based variational approach [30, 31, 32]. As a result, the obtained system is compatible with the law of energy dissipation, which makes it possible to design efficient and energetically stable numerical schemes.

One of the main challenges in phase-field method is to model the immersible two-phase flow with different densities. When the density ratio between the two phases is small, it could be handled by the Boussinesq approximation [33]. However, it could not be extended to the case with a large density ratio due to its underlying assumption [32]. One key problem arises from the inconsistency between the mass conservation and the incompressibility especially near the diffusive interface region. It was first pointed out by Lowengrub [34] and later by Shen et al. [32, 35]. Two main approaches are proposed to overcome this difficulty: one is based on volume averaged velocity; the other is based on the mass averaged velocity. For the volume averaged velocity model, the incompressiblity is assumed everywhere including the interfacial region [36, 37, 38, 39, 40]. An thermodynamically consistent and frame invariant model was developed by Abels et al. [41], where the mass conservation equation is modified with a mass correction term. On the other hand, for the mass averaged velocity method, the mass conservation is assured instead of incompressibility. This naturally yields the quasi-incompressible Navier-Stoker-Cahn-Hilliard (q-NSCH) model [34, 42], which in fact leads to a slightly compressible mixture only inside the interfacial region.

In the present paper, we first rederive and generalize the thermodynamcially consistent q-NSCH model in [34] from a variational point of view by combining with the Energy Variational Approach (EnVarA) [31, 43, 44] and Onsager’s Variation Principle [8, 21, 22, 41]. It starts from two functionals for the total energy and dissipation, together with the kinematic equations based on physical laws of conservation. The specific forms of the fluxes and stresses in the kinematic equations could be obtained by taking the time derivative of the total energetic functional and comparing with the predefined dissipation functional. More details could be found in [45]. In addition to bulk energy and dissipation, the energy and dissipation on the boundary are introduced to model the dynamics of contact lines. Our energy variational approach consistently yields both the correct bulk equations (the q-NSCH system) and a modified General Navier-Stokes Boundary Condition (GNBC) for the case of mass averaged velocity. The density effect on the contact line is explicitly modeled compared with the traditional GNBC [14, 15, 46, 47, 48, 49] in the case of volume averaged velocity, where the effect is modeled implicitly by the bulk and boundary interactions.

The second goal of our paper is to design an efficient energy stable scheme for the obtained q-NSCH system with large density ratio. There are not many such schemes developed for the MCLs. For the incompressible NSCH system, the development of such schemes may be found in [16, 46, 47, 50, 51, 52, 53, 54, 49] and only a few of them [47, 49, 52] are for variable density MCL models using the volume-averaged velocity (satisfying the incompressible condition in the whole domain). We shall develop an energy stable scheme for our thermodynamically consistent variable density q-NSCH system using the mass-averaged velocity. Based on the author’s previous works [13, 42, 55], we design a mass conservative $C^{0}$ finite element method for the q-NSCH system with a consistent discrete energy law. Thanks to a $\Delta p$ term in the quasi-incompressible condition, which is similar to the pressure stabilization of pseudo-compressibility methods [56, 57, 58], q-NSCH system does not need to satisfy the Babuska-Brezzi inf-sup condition [58, 59, 60, 61]. This may be considered as another benefit of our quasi-incompressible NSCH system.

The rest of paper is organized as follows. In Section 2, we present the thermodynamically consistent derivation of the q-NSCH system and its non-dimensionlization. The $C^{0}$ finite element algorithm for the q-NSCH system and the energy stable analysis are shown in Section 3. Section 4 presents the numerical results, including the convergence case study, and the examples of moving droplets and rising bubbles.

2 Mathematical Model

2.1 Mass-averaged velocity and laws of conservation

We consider a complex mixture consisting of two phase fluids with different densities. The interface of two fluids intersects with the wall $\partial\Omega_{w}$ at the contact line $\Gamma_{w}$ (see Fig. 1 (a)). Around the interface, we choose a control volume $V(t)$, where there are two phases labeled by $i=1,2$ with volume $V_{i}$ and mass $M_{i}$ (see Fig. 1 (b)). If the local average density of each phase is denoted by $\bar{\rho}_{i}=M_{i}/V$ and pure phase density is denoted by $\rho_{i}=M_{i}/V_{i}$, then density of mixture is

[TABLE]

Let $c_{i}=\frac{M_{i}}{M}$ be the mass fraction of each phase. Then we have [62]

[TABLE]

where $c=c_{1}$ is adopted in the last equality.

Remark 2.1**.**

Note that according to the definition (2), the mixture density $\rho$ is almost constant everywhere except in the interfacial region.

If we assume those two fluids move with velocities $\mathbf{u}_{i}$ ($i=1,2$), then the mass conservation of each phase inside the control volume is

[TABLE]

Introducing the mass averaged velocity as

[TABLE]

and combining with Eq.(1) yields the conservation of mass for the mixture

[TABLE]

Next, with an arbitrary volume $V(t)\in\Omega$, laws of conservation state

[TABLE]

Here, the first equation is the conservation of phase-field function (phase 1) and $\bm{j}_{c}$ is the flux of phase-field function. The second equation is the conservation of momentum where $\bm{\sigma}_{\eta}$ is the viscous stress and $\bm{\sigma}_{c}$ is the extra stress induced by two-phase interface due to nonzero $\nabla c$.

Thanks to the Reynolds transport theory [42, 45], Eqs.(5)-(7) yield the following kinematic equations in the domain $\Omega$

[TABLE]

where $\frac{D}{Dt}=\frac{\partial}{\partial t}+\mathbf{u}\cdot\nabla$ is the material derivative.

By definition of $\rho=\rho(c(\mathbf{x},t))$ in Eq. (2), above equations yield the quasi-incompressibility condition [34, 42]

[TABLE]

In the present case, we denote

[TABLE]

then the quasi-incompressibility condition is written as

[TABLE]

Remark 2.2**.**

Equations (13)-(14) show that the quasi-incompressibilty condition depends on the density difference of two fluids. When the two phases have the same density, i.e. $\rho_{1}=\rho_{2}$, it will consistently degenerate to the incompressibility condition. It makes a difference when two fluids have large density ratio and near the interfacial region [34, 42, 62].

On the boundary of domain $\partial\Omega$, the following boundary conditions are used

[TABLE]

where $\mathbf{u}^{s}=\mathbf{u}_{\tau}-\mathbf{u}_{w}$ with $\mathbf{u}_{\tau}=\mathbf{u}-(\mathbf{u}\cdot\mathbf{n})\mathbf{n}$ is the fluid slip velocity with respect to the wall, $\frac{D_{\Gamma}c}{Dt}=\frac{\partial c}{\partial t}+\mathbf{u}\cdot\nabla_{\Gamma}c$ is the surface material derivative, the Allen-Cahn type boundary condition is used for $c$ and $\nabla_{\Gamma}=\nabla-\mathbf{n}(\mathbf{n}\cdot\nabla)$ is surface gradient on the boundary $\partial\Omega$. The quantities $f_{\tau_{i}}$ and $J_{\Gamma}$ are to be determined. During the derivation, we assume the solid wall is fixed, i.e. $\mathbf{u}^{s}=\mathbf{u}_{\tau}$.

2.2 Model derivation

Now we start to derive the exact forms of $\bm{j}_{c}$, $\bm{\sigma}_{\eta}$, $\bm{\sigma}_{c}$ in Eq. (11), $f_{\tau_{i}}$ and $J_{\Gamma}$ in Eq. (18) by using energy variational method.

The total energy consists of the kinetic energy, the phase mixing energy and the energy on solid wall boundary $\partial\Omega_{w}$

[TABLE]

together with

[TABLE]

where $\lambda_{c}$ is the mixing energy density, $\gamma$ is the capillary width of the interface, $\theta_{s}$ is static contact angle and $\sigma$ is surface tension. The mixing energy $E_{mix}$ represents the competition between a homogeneous bulk mixing energy density term $G(c)$ (‘hydrophobic’ part) that enforces total separation of the two phases into pure components, and a gradient distortional term $\frac{|{\nabla}c|^{2}}{2}$ (‘hydrophilic’ part) that represents the nonlocal interactions between two phases and penalizes spatial heterogeneity.

The dissipation functional is composed of the dissipation due to fluid friction and irreversible mixing of two phases in bulk and the dissipation on the boundary

[TABLE]

where $\lambda(c)$ and $\eta(c)$ are the two Lamé coefficients, $\bm{D}_{\eta}=(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})/2$ is the strain rate, $\mathcal{M}$ is mobility coefficient in bulk, $\mathcal{M}_{\Gamma}$ is mobility coefficient on the wall, $\beta_{\Gamma}(c)$ is wall friction coefficient. In the present paper, $\eta(c)$ and $\beta_{\Gamma}$ are approximated by

[TABLE]

where $\eta_{i}$ and $\beta_{\Gamma,i}$ with $i=1,2$ are coefficients of each phase.

During the derivation, the following lemma is frequently used.

Lemma 2.1**.**

For a continuous function $f(\mathbf{x},t)$, if the density $\rho$ satisfies the conservation law (5) in the domain $\Omega$ and $\mathbf{u}\cdot\mathbf{n}=0$ on the boundary $\partial\Omega$, then we have

[TABLE]

By taking the time derivative of the total energetic functional, we have

[TABLE]

For the first term in (22), using the last two equations in Eq.(11) yields

[TABLE]

where we have introduced a Lagrangian multiplier $p$ with respect to the quasi-compressibility condition (12) and have used the boundary conditions $\mathbf{u}\cdot\mathbf{n}=0$ and $\bm{j}_{c}\cdot\mathbf{n}=0$. For the second term in (22), using the first equation in Eq.(11) and last two boundary conditions in Eq. (18) yields

[TABLE]

where $\mu=\lambda_{c}\left(\frac{dG}{dc}-\frac{1}{\rho}\gamma^{2}\nabla\cdot(\rho\nabla c)\right)$. The detailed derivations of Eqs.(23) -(24) are given in Appendix A. The last term $I_{3}$ in (22) yields

[TABLE]

Combining Eqs.(23) -(25), we obtain the derivative of the energy functional

[TABLE]

where we have defined

[TABLE]

Using energy dissipation law ${dE^{tot}}/{dt}=-\Delta$ [43, 63] and comparing (26) with the predefined dissipation functional in Eq.(21) yield

[TABLE]

By the definition of $\bm{\sigma}_{\eta}$ and $\bm{\sigma}_{c}$, the slip boundary condition (last equation in $(\ref{eq23})$) could be further written in the GNBC format

[TABLE]

To summarize, we have the following model for the two-phase flow with variable density for three unknowns $c,\mathbf{u},p$, in domain $\Omega$,

[TABLE]

with boundary conditions on $\partial\Omega$

[TABLE]

where $L(c)$ and $\bm{\sigma}_{\eta}$ are defined in (28) and (34).

Remark 2.3**.**

Note that in the above boundary conditions (41), the density effect on the contact line dynamics is explicitly modeled both in the boundary dynamics of phase-field and in velocity slip boundary condition through $L(c)$ term.

It is worth noting that the above system satisfies the following energy dissipation law.

Theorem 2.2**.**

If $c,\mathbf{u},p$ are smooth solutions of above system (36)-(41), then the following energy law is satisfied:

[TABLE]

Proof: The main idea of the proof is obtained by multiplying the phase-field equation (36a) by $\tilde{\mu}$, multiplying the chemical potential equation (36b) by $\frac{dc}{dt}$, multiplying the mass conservation equation (36c) by $p$, multiplying the Navier-Stokes equation (36d) by $\mathbf{u}$, and summing them up.

Taking the inner product of the phase-field equation (36a) with $\tilde{\mu}$ results in the following equation

[TABLE]

where we used the boundary condition $\partial_{n}\tilde{\mu}=0$ in (41).

Multiplying the chemical potential (36b) by $\rho\frac{Dc}{Dt}$ yields

[TABLE]

Summing up the above two equations, we have

[TABLE]

Multiplying the Navier-stokes equation (36d) by $\mathbf{u}$ followed by integration by parts, the rate of change of kinetic energy is calculated as

[TABLE]

where we have used the definition $\frac{D_{\Gamma}c}{Dt}=\frac{\partial c}{\partial t}+\mathbf{u}_{\tau}\cdot\nabla_{\Gamma}c$, $\mathbf{u}_{\tau}=\mathbf{u}^{s}$ and boundary conditions in (41).

From the derivation of $I_{2}$ in Appendix A, we get

[TABLE]

Combining the equations (45)-(47) leads to the final energy dissipation law. ∎

Remark 2.4**.**

In the above derivation, we have neglected the external body force, for example the gravity. If the effect of gravity needs to be taken into consideration, an extra gravitational potential should be added to the total energy

[TABLE]

where $g$ is the gravitational constant and $z$ is the vertical position. Using the fact that [62]

[TABLE]

the conservation of momentum equation (36d) is changed to

[TABLE]

2.3 Non-dimensionalization and Reformulation

In the following parts of the article, we assume that $\lambda=-\frac{2\eta}{3}$ for simplicity. Now we introduce the dimensionless variables

[TABLE]

Here $L^{*}$, $U^{*}$, $\rho^{*}$, $\eta^{*}$ and $\mathcal{M}^{*}$ are the characteristic scales of length, velocity, density, viscosity, and mobility coefficient, which are defined as

[TABLE]

For convenience, the hat symbol will be removed in the dimensionless quantities, and the dimensionless system of (36,41) is given by

[TABLE]

with boundary conditions

[TABLE]

where

[TABLE]

and with dimensionless parameters

[TABLE]

If we define

[TABLE]

then the system (53) could be rewritten as

[TABLE]

If we define the Sobolev spaces as follows

[TABLE]

and then above system satisfies the following energy dissipation law.

Theorem 2.3**.**

If $(c,\mu,\mathbf{u},\bar{p})\in\bm{W}_{b},$ are smooth solutions of above system (58) with boundary conditions (54), then the following energy law is satisfied:

[TABLE]

The proof is similar as Theorem 2.2. Here we omit the details.

Remark 2.5**.**

When the walls move, i.e. $\mathbf{u}_{w}\neq 0$, the above energy dissipation law has an extra term induced by external energy input

[TABLE]

3 Numerical Scheme and Analysis

3.1 Time-discrete primitive method

In this section, we present the numerical method of system (58) with boundary conditions (54) in the primitive variable formulation. Let $\Delta t>0$ denote the time step, and assume $(c^{n},\bar{\mu}^{n},\mathbf{u}^{n},\bar{p}^{n})$ are the solutions at the time $t=n\Delta t$. We then find the solutions at time $t=(n+1)\Delta t$ are $\left(c^{n+1},\bar{\mu}^{n+1},\mathbf{u}^{n+1},\bar{p}^{n+1}\right)$ that satisfy

[TABLE]

with boundary conditions

[TABLE]

where $\tilde{\mu}$ is defined in (57b) and we have used the notations

[TABLE]

For above discretization, it satisfies the following properties.

Lemma 3.4**.**

([42]) If $c^{n+1}$ is the solution of above system (64)-(65), then we have

[TABLE]

And the system (64)-(65) yields mass conservation for each component of binary fluid

[TABLE]

Theorem 3.5**.**

If $\left(c^{n+1},\mu^{n+1},\mathbf{u}^{n+1},\bar{p}^{n+1}\right)\in\bm{W}_{b}$ are solutions of above system (64) with boundary conditions (65), then the following energy law is satisfied:

[TABLE]

where

[TABLE]

is the discretized total energy.

Proof: Taking inner product of the first equation (64a) with $\frac{\Delta t}{\beta}\bar{\mu}^{n+1}$ results in the following equation

[TABLE]

where we used the boundary condition (65b) and the definition of $\bar{\mu}$ in (57b).

Multiplying the second equation (64b) with $\frac{c^{n+1}-c^{n}}{\beta}$ yields [42]

[TABLE]

where we use the results in above Lemma and the boundary condition (65a).

Multiplying the Navier-Stokes equation (64c) with $\Delta t\mathbf{u}^{n+1}$, we have

[TABLE]

where we used the slip boundary condition (65e) and the tensor calculation in Appendix C.

For the last term in above equation, combining the definition of $L^{n+1/2}$ in (65c) and equation (65a) yields

[TABLE]

Then equation (74) could be rewritten as

[TABLE]

Multiplying the last equation (64d) with $\frac{\Delta t}{\beta}\bar{p}^{n+1}$ yields

[TABLE]

Summing up equations (72), (73), (76) and (77) results

[TABLE]

where we used the definition of $\bar{\mu}$ (57b) and slip boundary condition (65e). ∎

3.2 Fully-discrete $C^{0}$ finite element scheme

The fully-discrete $C^{0}$ finite element scheme for this time-discrete primitive scheme (64)-(65) is presented in the section. For simplicity, we only consider a two-dimensional case here. It is straightforward to extend the results to three-dimensional case. The domain $\Omega$ is a bounded domain with Lipschitz-continuous boundary $\partial\Omega$. Specifically, we denote $\partial\Omega_{w}$ as the solid wall where the slip boundary condition is used. Let $\bm{W}^{h}_{b}=H^{h}\times H^{h}\times\bm{U}^{h}_{b}\times P^{h}$ be the finite dimensional space of $\bm{W}_{b}$ based on a given finite element discretization of $\Omega$. If we assume that $\rho^{n}\in L^{\infty}(\Omega)$ and positive [55], then the weak form of semi-discrete system (64) with boundary conditions (65) is the following: finding $(c_{h}^{n+1},\bar{\mu}_{h}^{n+1},\mathbf{u}_{h}^{n+1},\bar{p}_{h}^{n+1})\in\bm{W}_{b}^{h}$, such that

[TABLE]

for any $(\psi_{h},\chi_{h},\mathbf{v}_{h},q_{h})\in\bm{W}_{0}^{h}$.

Lemma 3.6**.**

The fully discretized system (79) satisfies mass conservation for each component of binary fluid

[TABLE]

Proof: Setting $\psi_{h}=q_{h}=\rho_{h}^{n+1}$ in Eqs.(79a) and (79d), we have

[TABLE]

Adding the above two equations and using Eqs. (13) and (68b) yields

[TABLE]

Using the boundary condition $\mathbf{u}^{n+1}_{h}\cdot\mathbf{n}=0$, we have the conservation of total mass

[TABLE]

Choosing $\psi_{h}=q_{h}=\hat{\rho}_{h}^{n+1}c^{n+1}_{h}$ in Eqs.(79a) and (79d), similarly we have

[TABLE]

Choosing $\psi_{h}=1$ in Eq.(79a) yields

[TABLE]

Adding the above two equations and using the velocity boundary condition $\mathbf{u}^{n+1}_{h}\cdot\mathbf{n}=0$, we have

[TABLE]

Theorem 3.7**.**

If $(c_{h}^{n+1},\bar{\mu}_{h}^{n+1},\mathbf{u}_{h}^{n+1},\bar{p}_{h}^{n+1})$ are solutions of the above system (79), then the following energy law is satisfied:

[TABLE]

where

[TABLE]

is the discretized total energy.

It can be proved by choosing

[TABLE]

in Eq. (79) and following the proof of Theorem 3.5.

4 Simulation Results

In this section, we present some numerical simulations using the aforementioned algorithm. Three cases are considered: Couette flow, moving droplets in shear flow and rising bubbles to illustrate the convergence rate, the effect of contact angle and the quasi-incompressibility of two-phase flow with large density ratio, respectively. All of the numerical simulations in this part are based on the proposed finite element scheme and implemented with the FreeFem++ [64].

4.1 Convergence Study: Couette Flow

We start with convergence test using Couette flow with different density and viscosity [47, 49] as in Fig. 2. The domain size is $[0,0.6]\times[0,0.1]$ The top and bottom walls move oppositely with $\mathbf{u}_{w}=(1,0)^{T}$. We do the convergence study for two-phase fluids with both low and high density ratios.

For the case of low density ratio, the parameters are listed as follows:

[TABLE]

[TABLE]

[TABLE]

We first present the convergence study for P1 element with $h=1/160,1/226,1/320,1/640$ and P2 element with $h=1/80,1/113,1/160,1/320$. The results with $h=1/640$ and $1/320$ are used as the reference solutions for P1 and P2 elements, respectively.

For the case of high density ratio, the parameters are listed as follows:

[TABLE]

[TABLE]

[TABLE]

The convergence rate for both P1 element and P2 element are shown in Table 2. It illustrates the 2nd-order for P1 element and 3rd-order for P2 element convergence rate in the sense of $L^{2}$ norm.

The profile of interface and velocity fields around steady state are shown in Fig. 3.

The fluid velocities on the wall are shown in Fig.4. It shows that for both low and high density ratio case, P1 element could yield consistent contact velocity with P2 element.

In Fig.5, we check the $L^{2}$ norm of $\nabla\cdot\mathbf{u}$ with different $\epsilon$. The results confirm that as $\epsilon$ decreases, the solution converges to the sharp interface incompressible fluids.

In Fig. 6, we check the total mass convergence of each phase in Lemma 3.6 for both low and high density ratios. It confirms that P1 and P2 elements could preserve the mass very well in both cases.

Then we set the wall velocity $\mathbf{u}_{w}=(0,0)^{T}$ to check the evolution of the total free energy when there is no input energy from outside. It is shown in Fig.7 that the free energy decreases over time for both methods and two density ratios, indicating that our schemes are energy stable.

4.2 Contact Angle Effect: moving droplet

In this example, we show the dynamics of an oil droplet in water with shear flow. The density ratio is $\rho_{1}:\rho_{2}=0.8:1$ and viscosity ratio is $\eta_{1}:\eta_{2}=1:1$. The other parameters are as follows

[TABLE]

[TABLE]

The domain size is $[0,4]\times[0,0.5]$ with adaptive mesh and $\Delta t=4\times 10^{-4}$. The initial profile is set to be a half circle

[TABLE]

In Figs. 8 and 10, the profiles of droplets under shear flow at different time are presented. For the acute contact angle cae (Fig. 8 ), the droplet is elongated by the shear flow force and hydrophilic force on the wall. The distance between two contact points increases over time (see Fig. 11 black curve) as a spreading droplet. While for the obtuse case (see Fig. 10), the hydrophobic force induced the shrink of contact lines on the wall. The distance between two contact points keeps decreasing (see Fig.11 blue curve). With the help of shear force, the droplet eventually detaches from the wall around $t=0.1$ and get stabilized at the center of the flow. When the contact angle is $90^{\circ}$ (Fig. 9), the competition between wall attraction force and bulk shear force first elongates the droplet and finally breaks the bubble around time $t=0.15$.

4.3 Large Density Ratio: Rising Bubble

As a last example, we carry out numerical simulation of an air bubble raising in water. The density ratio is set to be $\rho_{1}:\rho_{2}=0.001:1$ and viscosity ratio is $\eta_{1}:\eta_{2}=0.01:1$. The domain size is $(x,y)\in[0,0.15]\times[0,0.15]$ with mesh size $h=1/540$ and timestep $\Delta t=2\times 10^{-4}$. Parameters are listed as follows

[TABLE]

[TABLE]

The initial profile is set to be a half circle with radius 0.05 and center at $(0.075,0)$:

[TABLE]

The snapshots of interfaces with velocity fields and $\nabla\cdot\mathbf{u}$ profiles for bubbles with acute contact angle $\theta_{s}=60^{\circ}$ and obtuse contact angle $\theta_{s}=120^{\circ}$ are presented in Figs. 12 and 13, respectively. When the angle is acute, the attractive (hydrophilic) force from the wall competes with the buoyancy force and break the bubble. While for the obtuse case, the wall repulsive (hydrophobic) fore enhances the bubble rising under buoyancy force. The $\nabla\cdot\mathbf{u}$ profiles confirm that the quasi-impressible property of two-phase fluid with different density only happens around the interface due to the slightly mixing [42].

In Fig. 14, we show the dynamics of rising velocity $V_{c}=\frac{\int_{\Omega}u_{y}cd\mathbf{x}}{\int_{\Omega}cd\mathbf{x}}$ of bubble with different static contact angles. The vertical dash lines are the time when bubbles break ($\theta_{s}=60^{\circ},90^{\circ}$) or fully detach ($\theta_{s}=120^{\circ}$) from wall. It shows that the hydrophobic bubble (Blue line $\theta_{s}=120^{\circ}$ ) has a larger acceleration to form a sealing bubble. At $t=0.0384$, the bubble fully detaches from the wall. For the hydrophilic bubbles, in the beginning, the velocity increases slowly due to the competition between the hydrophilic force, the surface tension and the buoyancy force. The bubble is stretched into a tear shape which induces a larger velocity around the narrow neck region (see Fig. 12). The maximum velocity is achieved around the break time because the instantaneous response of the surface tension to the large surface deformation.

5 Conclusion

In this paper, we first derived the q-NSCH system for MCLs with variable density by using energy variational method consistently. GNBC for mass-averaged velocity is obtained during the variation due to the boundary dissipation.

Then we designed an energy stable $C^{0}$ finite element scheme to solve the obtained system. We also proved that the fully discrete scheme is mass conservative for each phase. Thanks to the quasi-incompressible condition with $\Delta p$ term, the finite element space for Navier-Stokes equations do not need to satisfy the Babuska-Brezzi inf-sup condition, as in the case of the pressure stabilization method for the standard Navier-Stokes equations.

Three examples are investigated numerically. The Couette flow test illustrates the 2nd-order for P1 element and 3rd-order for P2 element convergence rate in the sense of $L^{2}$ norm and the energy decay of the scheme. The contact angle effect on the droplet is illustrated by moving droplet in shear flow. Finally, a rising bubble is simulated to confirm the ability of our scheme to handle large density ratio and the quasi-incompresiblity only happens around the interface.

Acknowledgments

This work was partially supported by the NSFC (grant numbers 11771040, 11861131004), NSERC (CA) and the Fields Institute. Lingyue Shen was partially supported by the Chinese Scholarship Council for studying at the University of Dundee.

Appendix A Energy Variation Details

For the first term $I_{1}$ in (22), using the last two equations in Eq.(11) yields

[TABLE]

where we have introduced a Lagrangian multiplier $p$ with respect to the constraint (14) and have used the boundary condition $\mathbf{u}\cdot\mathbf{n}=0$ and $\bm{j}_{c}\cdot\mathbf{n}=0$.

For the second term $I_{2}$ in (22), using the first equation in Eq.(11) and last two boundary conditions in Eq. (18) yields

[TABLE]

where $\mu=\lambda_{c}\left(\frac{dG}{dc}-\frac{1}{\rho}\gamma^{2}\nabla\cdot(\rho\nabla c)\right)$.

Appendix B Proof of Lemma 2.1

Proof:

[TABLE]

Appendix C Tensor Calculation

[TABLE]