Keywords

1 Introduction

Two-phase flow is omnipresent in many natural and industrial processes, especially for the petroleum industry, the two-phase flow is throughout the whole upstream production process including oil and gas recovery, transportation and refinery, e.g. [1].

As a critical component in the two-phase fluid system, the interface is usually considered as a free surface, its dynamics is determined by the usual Young-Laplace junction condition in classical interface tracking or reconstruction approaches, such as Level-set [2], volume-of-fluid [3] and even some advanced composite method like VOSET [4].

But when it traced back to 19th century, Van der Waals [5] provided a new alternative point of view that the interface has a diffuse feature, namely non-zero thickness. It can be implicitly characterized by scalar field, namely phase filed, taking constant values in the bulk phase areas and varying continuously but radically across a diffuse front. Within this thin transition domain, the fluids are mixed and store certain quantities of “mixing energy” in this region. Thus, unlike other methods proposed graphically, phase dynamics is derived from interface physics by energy variational approach regardless of the numerical solution, which gives rise to coupled nonlinearly well-posed system at partial differential equation continuous form that satisfies thermodynamically consistent energy dissipation laws. Then there is possibility for us to design numerical scheme preserving energy law in discrete form [6].

The diffuse interface approach excels in some respects of handling two-phase flow among other available methods. Firstly, it is based on the principle of energy minimization. Hence it can deal with moving contact lines problems and morphological changes of interface in a natural way effortlessly, such as droplet coalescence or break-up phenomena. Secondly, we can benefit from the simplicity of formulation, ease of numerical implementation without explicitly tracking or reconstructing interface, also the capability to explore essential physics at the interfacial domain. The accessibility of modeling various material properties or complex interface behaviors directly by introducing appropriate energy functions. Therefore, enforcing certain rheological fluid or modeling polymeric solutions or viscoelastic behaviors would be alluring feature naturally. For these benefits, the diffuse interface model attracted substantial academic attention in recent years, a great number of the advanced and cutting-edge researches are conducted corresponding to partial immiscible multi-components flow based on phase field theory [7] and thermodynamically consistent diffuse interface model for two-phase flow with thermo-capillary [8] et al.

The classical diffuse interface model for cases of two-phase incompressible viscous Newtonian fluids is known as the model H [9]. It has been successfully applied to simulate flows involving incompressible fluids with same densities for both phase components. This model is restricted to the matched density case using Boussinesq approximation. Unlike the matched density case, when it comes to the case with big density ratio, the incompressibility cannot guarantee mass conservation any longer in this model. Therefore, the corresponding diffuse interface model with the divergence free condition no longer preserve an energy law. Thus, a lot of further works have been done by (1998) Lowengrub [10], (2002) Boye [11], (2007) Ding [12], (2010) Shen [13] and most recently Benchmark computations were carried out by (2012) Aland [14].

Generally there are two kinds of approaches to deal with variable densities problem, one is that material derivative of momentum equation written in one kind form that takes density variations into consideration without resorting to mass conservation to guarantee stability of energy proposed by Guermond [15]. Another approach is proposed by Abels [16]. The approach introduces an intermediate velocity to decouple the Cahn-Hilliard equation and Navier-Stokes equation system in Minjeaud’s paper [17], and recently this approach is applied in Shen [18] to simulate the model in [16]. However, the schemes proposed in [17, 18] employ the intermediate velocity in the Cahn-Hilliard equation only, imposing the mass balance equation to ensure the discrete energy-dissipation law. Very recently, in Kou [19] the schemes that the intermediate velocity is applied in both mass balance equations and the momentum balance equation are developed to simulate the multi-component diffuse-interface model proposed in [20] to guarantee the consistency between the mass balance and energy dissipation. In this paper, we extend this treatment to the model in [16]. However, this extension is not trivial due to a crucial problem that Cahn-Hilliard equation is not equivalent to mass balance equation. In order to deal with this problem, a novel scheme applying the intermediate velocity in Navier-Stokes equation will be proposed in this paper.

The rest part of this paper is organized as follows. In Sect. 2 we introduce a diffuse interface model for two-phase flow with variable densities and viscosities in detail; In Sect. 3 we propose a brand new numerical scheme for solving the coupled Navier–Stokes and Cahn–Hilliard equation system based on this model; In Sect. 4 some numerical results are demonstrated in this part to validate this scheme comparing with benchmark and to exam the accuracy, discrete energy decaying tendency. Other cases and numerical performances will be investigated to show the robustness of the novel scheme.

2 Mathematical Formulation and Physical Model

The phase field diffuse interface model with variable densities and viscosities can be described through the following Cahn-Hilliard equation coupled with Navier-Stokes equation. An introduced phase field variable \( \phi \), namely order parameter, defined over the domain, identifies the regions occupied by the two fluids.

$$ \phi \left( {{\mathbf{x}},t} \right) = \left\{ {\begin{array}{*{20}c} 1 & {fluid \;1} \\ { - 1} & {fluid \;2} \\ \end{array} } \right. $$
(1)

With a thin smooth transition front of thickness \( \varepsilon \) bridging two fluids, the microscopic interactions between two kinds of fluid molecules rules equilibrium profiles and configurations of interface mixing layer neighboring level-set \( \Gamma_{t} = \left\{ {\phi :\left( {{\mathbf{x}},t} \right) = 0} \right\} \). For the situation of isotropic interactions, the following Ginzburg–Landau type of Helmholtz free energy functional is given by the classical self-consistent mean field theory in statistical physics [21]:

$$ W\left( {\phi ,\nabla \phi } \right) = \lambda \mathop \smallint \nolimits_{\Omega } \left( {\frac{1}{2}\left\| {\nabla \phi } \right\|^{2} + F\left( \phi \right)} \right)dx $$
(2)

The foremost term in right hand side represents the effect of mixing of interactions between the materials, and the latter one implies the trend of separation. Set the Ginzburg–Landau potential in the usual double-well form \( {\text{F}}\left( \phi \right) = \frac{{\phi^{2} - 1}}{{4\varepsilon^{2} }} \). \( \lambda \) means mixing energy density, \( \varepsilon \) is the capillary width of interface between two phases. If we focus on one-dimensional interface and assume that total diffusive mixing energy in this domain equals to traditional surface tension coefficient:

$$ \sigma = \uplambda \mathop \smallint \nolimits_{ - \infty }^{ + \infty } \left\{ {\frac{1}{2}\left( {\frac{d\phi }{dx}} \right)^{2} + F\left( \phi \right)} \right\}dx $$
(3)

The precondition that diffuse interface is at equilibrium is valid, then we can get the relationship among surface tension coefficient \( \sigma \), capillary width \( \varepsilon \) and mixing energy density \( \lambda \):

$$ \sigma = \frac{{\left( {2\sqrt 2 } \right) }}{3}\frac{\lambda }{\varepsilon } $$
(4)

The evolution of phase field (\( \phi ) \) is governed by the following Cahn-Hilliard equations:

$$ \left\{ {\begin{array}{*{20}c} {\phi_{t} + \nabla \cdot \left( {{\mathbf{u}}\phi } \right) = M\Delta \mu } \\ {\mu = \frac{\delta W}{\delta \phi } = f\left( \phi \right) - \Delta \phi } \\ \end{array} } \right. $$
(5)

Where \( \mu \) represents the chemical potential, namely \( \frac{\delta W}{\delta \phi } \) that indicates the variation of the energy \( W \) with respect to \( \phi \); the parameter \( M \) is a mobility constant related to the diffusivity of bulk phases and \( f\left( \phi \right) = F^{\prime}\left( \phi \right) \)

The momentum equation for the two-phase system is presented as the usual form

$$ \left\{ {\begin{array}{*{20}c} {\rho \left( {{\mathbf{u}}_{t} + \left( {{\mathbf{u}} \cdot \nabla } \right){\mathbf{u}}} \right) = \nabla \cdot \tau } \\ {\tau = \eta D\left( {\mathbf{u}} \right) - pI +\varvec{\tau}_{e} } \\ {\varvec{\tau}_{e} = - \lambda \left( {\nabla \phi \otimes \nabla \phi } \right)} \\ \end{array} } \right. $$
(6)

with the identical equation

$$ \begin{aligned} & \nabla \cdot \left( {\nabla {\mathbf{u}} \otimes \nabla {\mathbf{u}}} \right) = \left( {\Delta \phi - f\left( \phi \right)} \right) + \frac{1}{2}\nabla \left( {\left\| {\nabla \phi } \right\|^{2} + F\left( \phi \right)} \right) = - \mu \nabla \phi \\ & + \,\frac{1}{2}\nabla \left( {\left\| {\nabla \phi } \right\|^{2} + F\left( \phi \right)} \right) = \phi \nabla \mu + \frac{1}{2}\nabla \left( {\left\| {\nabla \phi } \right\|^{2} + F\left( \phi \right) - \phi \mu } \right) \\ \end{aligned} $$
(7)

The second term in Eq. (7) can be merged with the pressure gradient term \( p \), then the pressure of the momentum equation should be a modified pressure \( p^{*} \), It denotes that: \( p^{*} = p + \frac{1}{2}\nabla \left( {\left| {\nabla \phi } \right|^{2} + F\left( \phi \right) - \phi \mu } \right) \). The \( p \) is represented the modified one in the following contents for unity.

2.1 Case of Matched Density

The governing equations system:

$$ \left\{ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\phi_{t} + \nabla \cdot \left( {{\mathbf{u}}\phi } \right) = M\Delta \mu } \\ {\mu = f\left( \phi \right) - \Delta \phi } \\ \end{array} } \\ {{\mathbf{u}}_{t} + \left( {{\mathbf{u}} \cdot \nabla } \right){\mathbf{u}} = \nabla \cdot \eta D\left( {\mathbf{u}} \right) - \nabla p - \phi \nabla \mu } \\ {\nabla \cdot {\mathbf{u}}} \\ \end{array} } \right. $$
(8)

A set of appropriate boundary condition and initial condition is applied to the above system: no-slip boundary condition for momentum equation and the period boundary condition for the Cahn-Hilliard equation.

$$ {\mathbf{u}}|_{\partial \Omega } = 0 \;\; \left. {\frac{\partial \phi }{\partial n}} \right|_{\partial \Omega } = 0 \;\; \left. {\frac{\partial \mu }{\partial n}} \right|_{\partial \Omega } = 0 $$
(9)

Also the initial condition:

$$ {\mathbf{u}}|_{t = 0} = {\mathbf{u}}_{0} \phi |_{t = 0} = \phi_{0} $$
(10)

If the density contrast of two phases is relatively little, a common approach is to employ the Boussinesq approximation [22], replacing momentum equation in equation system (8) by

$$ \rho_{0} \left( {{\mathbf{u}}_{t} + \left( {{\mathbf{u}} \cdot \nabla } \right){\mathbf{u}}} \right) = \nabla \cdot \eta D\left( {\mathbf{u}} \right) - \nabla p - \phi \nabla \mu + \frac{\phi + 1}{2}\delta \rho {\mathbf{g}} $$
(11)

Set the background density \( \rho_{0} = \left( {\rho_{1} + \rho_{2} } \right) \) and term \( \frac{\phi + 1}{2}\delta \rho {\mathbf{g}} \) is an additional body force term in charge of the equivalent gravitational effect caused by density difference. Since the density \( \rho_{0} \) distributed everywhere in this field do not change respect to time. If the divergence of the velocity field \( \nabla \cdot {\mathbf{u}} = 0 \) holds. Then basic mass conservation \( \rho_{t} + \nabla \cdot \left( {\rho {\mathbf{u}}} \right) = 0 \) is a natural consequence of incompressibility.

By inner product operation of 1st 2nd and 3rd equation of system (8) with \( - \mu \), \( \phi_{t} \), \( {\mathbf{u}} \) respectively and summation of these three results, It is easily to conclude that system (8) admits the following energy dissipation law:

$$ \frac{d}{dt}\mathop \smallint \nolimits_{\Omega } \left( {\frac{1}{2}\left\| {\mathbf{u}} \right\|^{2} + \frac{\lambda }{2}\left\| {\nabla \phi } \right\|^{2} + \lambda F\left( \phi \right)} \right)dx = - \mathop \smallint \nolimits_{\Omega } \left( {\frac{\eta }{2}\left\| {D\left( u \right)} \right\|^{2} + M\left\| {\nabla \mu } \right\|^{2} } \right)dx $$
(12)

2.2 Case of Variable Density

Now we consider the case where the density ratio is so large that the Boussinesq approximation is no longer in effect. Here introduced a phase field diffuse interface model for incompressible two-phase flow with different densities and viscosity proposed by Abels [16].

$$ \left\{ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\phi_{t} + \nabla \cdot \left( {{\mathbf{u}}\phi } \right) = M\Delta \mu } \\ {\mu = f\left( \phi \right) - \Delta \phi } \\ \end{array} } \\ {\rho {\mathbf{u}}_{t} + \left( {\rho {\mathbf{u}} + \varvec{J}} \right) \cdot \nabla {\mathbf{u}} = \nabla \cdot \eta D\left( {\mathbf{u}} \right) - \nabla p - \phi \nabla \mu } \\ {\nabla \cdot {\mathbf{u}} = 0} \\ \end{array} } \right. $$
(13)

among them

$$ \varvec{J} = \frac{{\rho_{2} - \rho_{1} }}{2}M\nabla \mu $$
(14)

The density and viscosity is the function of phase parameter.

$$ \rho \left( \phi \right) = \frac{{\rho_{1} - \rho_{2} }}{2}\phi + \frac{{\rho_{1} + \rho_{2} }}{2} \eta \left( \phi \right) = \frac{{\eta_{1} - \eta_{2} }}{2}\phi + \frac{{\eta_{1} + \rho \eta_{2} }}{2} $$
(15)

The mass conservation property can be derived from Eqs. (14), (15) and (16).

$$ \rho_{t} + \nabla \cdot \left( {\rho {\mathbf{u}}} \right) + \nabla \cdot \varvec{J} = 0 $$
(16)

The NSCH governing system holds thermodynamically consistency and energy law. We can obtain the following energy dissipation law:

$$ \frac{d}{dt}\mathop \smallint \nolimits_{\Omega } \left( {\frac{\rho }{2}\left\| {\mathbf{u}} \right\|^{2} + \frac{\lambda }{2}\left\| {\nabla \phi } \right\|^{2} + \lambda F\left( \phi \right)} \right)dx = - \mathop \smallint \nolimits_{\Omega } \left( {\frac{\eta }{2}\left\| {D\left( u \right)} \right\|^{2} + \left\| {M\nabla \mu } \right\|^{2} } \right)dx $$
(17)

If we add the gravity in this domain, such as modeling topological evolution of a single bubble rising in a liquid column, the total energy must contain the potential energy. Then this energy dissipation law can be expressed as follow:

$$ \frac{{dE_{tot} }}{dt} = \frac{d}{dt}\mathop \smallint \nolimits_{\Omega } \left( {\frac{\rho }{2}\left\| {\mathbf{u}} \right\|^{2} + \frac{\lambda }{2}\left\| {\nabla \phi } \right\|^{2} + \lambda F\left( \phi \right) - \rho {\mathbf{g}}{\text{y}}} \right)dx \le 0 $$
(18)

3 Decoupled Numerical Scheme

3.1 Time Discretization

In matched density case, for simplicity of presentation, we will assume that \( \eta_{1} = \eta_{2} = \eta \). Given initial conditions \( \phi^{0} \), \( {\mathbf{u}}^{0} \), \( p^{0} \) we compute \( \phi^{{{\text{k}} + 1}} \), \( \mu^{{{\text{k}} + 1}} \), \( p^{{{\text{k}} + 1}} \), \( \widetilde{{\mathbf{u}}}^{{{\text{k}} + 1}} \), \( {\mathbf{u}}^{{{\text{k}} + 1}} \) for \( {\text{n}} \ge 0 \). Here is an additional term to the convective velocity introduced based on the idea from [17]. Then the intermediate velocity term \( \widehat{{\mathbf{u}}}^{\text{k}} = {\mathbf{u}}^{\text{k}} - \delta t\frac{{\phi^{k} \nabla \mu^{k + 1} }}{{\rho^{k} }} \) makes the Cahn-Hilliard equation and the Navier-Stokes equation decoupled fundamentally. The novel scheme can be described as below:

$$ \left\{ {\begin{array}{*{20}c} {\frac{{\phi^{k + 1} - \phi^{k} }}{\delta t} + \nabla \cdot \left( {\widehat{{\mathbf{u}}}^{\text{n}} \phi^{k + 1} } \right) = M\Delta \mu^{k + 1} } \\ {\mu^{k + 1} = \lambda \left( {f_{e} \left( {\phi^{k} } \right) + f_{c} \left( {\phi^{k + 1} } \right) - \Delta \phi^{k + 1} } \right).or. } \\ { \mu^{k + 1} = \lambda \left( {f\left( {\phi^{k} } \right) - \Delta \phi^{k + 1} } \right) + \frac{\lambda }{{ \varepsilon^{2} }}\left( {\phi^{k + 1} - \phi^{k} } \right) } \\ {{\mathbf{n}} \cdot \nabla \left. {\phi^{k + 1} } \right|_{\partial \Omega } = 0 \;\;\; {\mathbf{n}} \cdot \nabla \left. {\mu^{k + 1} } \right|_{\partial \Omega } = 0} \\ \end{array} } \right. $$
(19)

We add a stabilizing term \( \frac{\lambda }{{ \varepsilon^{2} }}\left( {\phi^{k + 1} - \phi^{k} } \right) \) or treat the term by convex splitting method \( f\left( \phi \right) = f_{e} \left( \phi \right) + f_{c} \left( \phi \right) \). Then the time step for computation will not be strictly limited in extreme range by the coefficient Capillary width \( \varepsilon \).

$$ \left\{ {\begin{array}{*{20}c} {\rho_{0} \frac{{\widetilde{{\mathbf{u}}}^{k + 1} - {\mathbf{u}}^{k} }}{\delta t} + \left( {{\mathbf{u}}^{k} \cdot \nabla } \right)\widetilde{{\mathbf{u}}}^{k + 1} = \nabla \cdot \eta_{0} D\left( {\widetilde{{\mathbf{u}}}^{k + 1} } \right) - \nabla p^{k} - \phi^{k} \nabla \mu^{k + 1} + F_{m} } \\ {\left. {\widetilde{{\mathbf{u}}}^{k + 1} } \right|_{\partial \Omega } = 0} \\ \end{array} } \right. $$
(20)

Then we can get the pressure by solving a constant coefficient Poisson equation and correct the velocity filed to satisfy divergence free condition.

$$ \left\{ {\begin{array}{*{20}c} {\rho_{0} \frac{{{\mathbf{u}}^{k + 1} - \widetilde{{\mathbf{u}}}^{k + 1} }}{\delta t} = - \nabla \left( {p^{k + 1} - p^{k} } \right)} \\ {\nabla \cdot {\mathbf{u}}^{k + 1} = 0} \\ {\left. {{\mathbf{u}}^{k + 1} } \right|_{\partial \Omega } = 0} \\ \end{array} } \right. $$
(21)

For variable density case, [15, 16] serve as incentive for the novel numerical scheme below. To deal with the variable densities and ensure numerical stability, we have to define a cut-off function \( \widehat{\phi } \) for the phase order parameter at first place.

$$ \phi = \left\{ {\begin{array}{*{20}c} \phi & {\left| \phi \right| \le 1} \\ {sign\left( \phi \right)} & {\left| \phi \right| > 1} \\ \end{array} } \right. $$
(22)

Given initial conditions \( \phi^{0} ,{\mathbf{u}}^{0} ,p^{0} ,\rho_{0} ,\eta_{0} \) we compute \( \phi^{{{\text{k}} + 1}} ,\mu^{{{\text{k}} + 1}} ,p^{{{\text{k}} + 1}} {\mathbf{u}}^{{{\text{k}} + 1}} ,\rho^{k + 1} , \) \( \eta^{k + 1} \) for \( {\text{n}} \ge 0 \).The discretization of the Cahn-Hilliard part is same with the matched density as Eq. (19).

We update the density and viscosity by cut-off function

$$ \left\{ {\begin{array}{*{20}c} {\rho \left( {\phi^{n + 1} } \right) = \frac{{\rho_{1} - \rho_{2} }}{2}\widehat{\phi }^{n + 1} + \frac{{\rho_{1} + \rho_{2} }}{2} } \\ {\eta \left( {\phi^{n + 1} } \right) = \frac{{\eta_{1} - \eta_{2} }}{2}\widehat{\phi }^{n + 1} + \frac{{\eta_{1} + \rho \eta_{2} }}{2}} \\ \end{array} } \right. $$
(23)

For the momentum equation part, we use

$$ \left\{ {\begin{array}{*{20}c} {\rho^{k} \frac{{{\mathbf{u}}^{k + 1} - {\mathbf{u}}^{k} }}{\delta t} + \left( {\rho^{\text{k}} \widehat{{\mathbf{u}}}^{k} + \varvec{J}^{k} } \right) \cdot \nabla {\mathbf{u}}^{k + 1} = \nabla \cdot \eta D\left( {{\mathbf{u}}^{k + 1} } \right) - \nabla p^{k} } \\ { - \phi^{k} \nabla \mu^{k + 1} + \frac{1}{2}\left( {\frac{{\rho_{1} + \rho_{2} }}{2}\nabla \cdot \widehat{{\mathbf{u}}}^{k} } \right){\mathbf{u}}^{k + 1} } \\ {\left. {{\mathbf{u}}^{k + 1} } \right|_{\partial \Omega } = 0} \\ \end{array} } \right. $$
(24)

together with

$$ \varvec{J}^{k} = \frac{{\rho_{2} - \rho_{1} }}{2}M\nabla \mu^{k} $$
(25)

For saving the computer time consuming and stability, we adopt schemes based on pressure stabilization.

$$ \left\{ {\begin{array}{*{20}c} {\Delta \left( {p^{k + 1} - p^{k} } \right) = \frac{\theta }{\delta t}\nabla \cdot {\mathbf{u}}^{n + 1} } \\ {\left. {\nabla \left( {p^{k + 1} - p^{k} } \right)} \right|_{\partial \Omega } = 0} \\ {\theta = \frac{1}{2}{ \hbox{min} }\left( {\rho_{2 } ,\rho_{1} } \right)} \\ \end{array} } \right. $$
(26)

Pressure stabilization at the initial stage might cause velocity without physical meaning because it cannot satisfy solenoidal condition strictly. If we use pressure correction method, we have to face the non-linear Poisson equation. Thus, the solution must cost much more time.

3.2 Spatial Discretization

For 2-D cases, the computational domain is \( \Omega = \left( {0,L_{x} } \right) \times \left( {0,L_{y} } \right) \), the staggered grid are used for spatial discretization. The cell centers are located on

$$ x_{i} = \left( {i - \frac{1}{2}} \right)h_{x} \;\; i = 1, \ldots ,n_{x} \;\;\;\;\;\;\; y_{i} = \left( {j - \frac{1}{2}} \right)h_{y} \;\; i = 1, \ldots ,n_{y} $$

Where \( h_{x} \) and \( h_{y} \) are grid spacing in \( x \) and \( y \) directions. \( n_{x } \), \( n_{y } \) are the number of grids along \( x \) and \( y \) coordinates respectively. In order to discretize the coupled Cahn-Hilliard and Navier-Stokes system, the following finite volume method is introduced (Fig. 1).

Fig. 1.
figure 1

The staggered grid based on finite volume method

Full size image
$$ \begin{array}{*{20}l} {U_{h} = \left. {\left( {x_{{i - \frac{1}{2}}} , y_{j} } \right)} \right|i = 1, 2, \ldots ,n_{x} + 1;} \hfill & {j = 1,{ 2}, \ldots ,n_{y} } \hfill \\ {V_{h} = \left. {\left( {x_{i} , y_{{j - \frac{1}{2}}} } \right)} \right|i = 1, 2, \ldots ,n_{x} ;} \hfill & {j = 1,{ 2}, \ldots ,n_{y} + 1} \hfill \\ {P_{h} = \left. {\left( {x_{{i - \frac{1}{2}}} , y_{j} } \right)} \right| i = 1,{ 2}, \ldots ,n_{x} ;} \hfill & {j = 1,{ 2}, \ldots ,n_{y} } \hfill \\ \end{array} $$

Where \( P_{h} \) is cell-centered space, \( U_{h} \) and \( V_{h} \) are edge-center space, \( \left( {\phi ,\mu ,p,\rho ,\eta } \right) \in P_{h } \), \( u \in U_{h} \), \( v \in V_{h} \). Some common differential and averaged operators are used to interpolation of these physical variables from one space to another space which are not discussed in detail here.

4 Numerical Results

4.1 Validation of the Novel Scheme

Case1: A Bubble Rising in Liquid in 2D Domain

There is a rectangular domain \( \Omega = \left( {0,1} \right) \times \left( {0,2} \right) \) filled with two-phase incompressible fluid and a lighter bubble (with density \( \rho_{1} \) and dynamic viscosity \( \eta_{1} \)) in a heavier medium (with density \( \rho_{2} \) and dynamic viscosity \( \eta_{2} \)) rises from a fixed position initially. As it described in benchmark test paper [23], the boundary condition imposed on the vertical walls and top wall are replaced by free slip condition. The physical parameters for case 1 follows Table 1.

Table 1. The physical parameters for numerical test case 1.
Full size table

The initial bubble is perfect round with radius \( {\text{r}} = 0.25 \) and its center is set at the point \( \left( {x_{0} ,y_{0} } \right) = \left( {0.5,1} \right) \). The initial profile of \( \phi \) is set as

$$ \phi \left( {x,y} \right) = - { \tanh }\,(\frac{1}{\sqrt 2 \varepsilon }\left( {\sqrt {\left( {x - x_{0} } \right)^{2} + \left( {y - y_{0} } \right)^{2} } - r} \right) $$
(27)

We must note these parameters have been through the non-dimensionalization. Mobility coefficient is an additional numerical parameter, which is not appeared in the sharp interface model. The value is chosen in a rational range for comparison of different spatial step and interface thickness. Furthermore, the interface thick \( \varepsilon \) is chosen proportional to \( h \). The energy density parameter \( \lambda \) can be calculated by the surface tension coefficient \( \sigma \) through the Eq. (4). The time step \( \delta t = 0.0001 \) (Fig. 2).

Fig. 2.
figure 2

Snapshots of the bubble evolution with velocity vector field (computed on h = 1/150 grid).

Full size image

The bubble shape at the time point t = 3.0 calculated by the novel scheme is compared with the solution from the benchmark paper by the level-set method [23] and the diffuse interface method [14] in Fig. 3(a). Different bubble shapes at the time t = 3.0 are compared from the coarsest grid (h = 1/50, \( \varepsilon = 0.02 \)) to the finest grid (h = 1/200, \( \varepsilon = 0.005 \)).

Fig. 3.
figure 3

Bubble shapes at t = 3 for the novel scheme comparing with the level-set and diffuse interface benchmark results provided in [14, 23] (a); bubble shapes at t = 3 solved by grid with different refinement level and different interface thickness (b).

Full size image

The shapes of bubble differ distinctly for different values of interface thickness \( \varepsilon \). But they seem to be convergent so that there is no significant differences for the finest grid and the case with \( \varepsilon = 0.008 \). We can also remark that the bubble shape from novel scheme is quite approximate to the benchmark level-set and diffuse interface results. But it is clearly not sufficient to only look at the bubble shapes, therefore we use some previously defined benchmark quantities to validate the new scheme rigorously.

I. Center of the Mass:

Various positions of points can be used to track the motion of bubbles. The most common way is to use the center of mass defined by

$$ y_{\text{c}} = \frac{{\smallint_{\phi > 0} ydx}}{{\smallint_{\phi > 0} 1dx}} $$
(28)

with y as the vertical coordinate of \( {\mathbf{x}}\left( {x,y} \right) \)

II. Rise Velocity

\( v \) is the vertical component of the bubble’s velocity \( {\mathbf{u}} \). Where \( \phi > 0 \) denotes the region that bubble occupies. The velocity is volume average velocity of bubble.

$$ v_{c} = \frac{{\smallint_{\phi > 0} vdx}}{{\smallint_{\phi > 0} 1dx}} $$
(29)

Combining the contours given in Fig. 3 and Figs. 4, 5 and 6 about the mass center and rise velocity of the lighter bubble before t = 3, a significant increase in rise velocity can be observed at the initial stage and then the velocity decrease slowly to a constant value when the time advances to t = 3 gradually. The shape of bubble will reach to a temporary steady state when the rise velocity keeps constant. The mass center of the bubble could be recognized as a linear function of time asymptotically after it is higher than y = 0.6.

Fig. 4.
figure 4

Center of the mass of the bubble for the novel scheme comparing with the level-set and diffuse interface benchmark result provided in [14, 23] (a); Solutions of center of mass from grids with different refinement level and different interface thickness (b).

Full size image
Fig. 5.
figure 5

Rise velocity of the bubble for the novel scheme comparing with the level-set and diffuse interface benchmark results (a); close-up of rise velocity at maximum values area (b).

Full size image
Fig. 6.
figure 6

Solutions of rise velocity of the bubble from grids with different refinement level and different interface thickness (a); close-up of rise velocity at maximum values area (b).

Full size image

Although the difference is visible for the coarsest grid and thickest interface on the velocity plot. It is obvious to see that results become closer and closer to the line corresponding to the computation on the finest grid with refinement. The results calculated by the new scheme shows good agreement with the level-set solution [23] and benchmark diffuse interface method [14].

The decaying trend of discrete energy in Fig. 7 confirms that the proposed scheme is energy stable. The whole system reaches the equilibrium state at t = 25.

Fig. 7.
figure 7

Energy dissipation of the whole system

Full size image

Figure 8 gives the evolution of free energy and kinetic energy respectively. At the early stage of bubble rising, kinetic energy and free energy rise dramatically, which come from part of the reduced gravity potential energy. Then the velocity keep constant to some extent at next stage. The bubble shape also change into a relative steady state. When the bubble almost touching the top lid. The kinetic energy gives a considerable decrease to the zero. Then the gas phase will evolve to a stratified state finally under the lead of the diffusion in the Cahn Hilliard equation.

Fig. 8.
figure 8

Variation of free energy and kinetic energy of the whole system

Full size image

Case2: Novel Scheme for Matched Density with Boussinesq Approximation

In the section, We simulate a physical problem with matched viscosities and a relative low density contrast \( \left( {\rho_{1} = 1 ;\rho_{2} = 10} \right) \) which ensures the Boussinnesq approximation is applicable. We set 2d bubble diameter \( d = 1.0 \), \( {\text{g}} = 1.0 \), \( \eta_{1} = \eta_{2} = 0.1 \), \( \lambda = 4 \times 10^{ - 3} \), Mobility = 0.02 and \( \varepsilon = 0.008 \). The incompressible fluid in a rectangular domain \( \Omega = \left( {0,2} \right) \times \left( {0,4} \right) \) with initially a lighter bubble center located on \( \left( {x_{0} ,y_{0} } \right) = \left( {1,1.5} \right) \). These dimensionless parameters are set according to the cases in [13]. The mesh size is \( 200 \times 400 \) and the boundary condition at the vertical wall is no-slip boundary condition (Fig. 9).

Fig. 9.
figure 9

The shape and displacement comparison of the bubble calculated by novel scheme with the results presented in [13] for matched density case.

Full size image

It is easy to find that the shape of the bubble and the vertical displacement at different time step is pretty similar with the reference contour provided in [13] by GUM/PSM method. The reference only have a contour line at certain \( \phi \). For a beautiful presentation of the integral contour calculated by the novel scheme here, we adjust interface thickness \( \varepsilon = 0.008 \). The novel scheme can be employed in the situation.

4.2 Robustness Test of the Novel Scheme

Case3: Examining the Performance for Big Density Contrast

We now consider two-phase incompressible fluid with the same initial condition and boundary condition with the case1. But density and viscosity contrast is much more violent in this case (Table 2 and Fig. 10).

Table 2. The physical parameters for numerical test case 3.
Full size table
Fig. 10.
figure 10

The shapes and displacement comparison of the bubble calculated by novel scheme with the benchmark level-set results and the contours at t = 3.0 provided by three common commercial software in [23]

Full size image

The break-up of bubble happen before the time t = 3.0 in the level-set benchmark. So it is not appropriate to compare the results from the new scheme with the contour. But from the results of some common commercial computing software [23]. It’s not that difficult to find the shape of bubble and the vertical displacement at t = 3.0 solved by our scheme is pretty close to them. Although it could be some slight diffusion on the interface caused by the Cahn-Hilliard system itself. The case shows the robustness of the novel scheme proposed in this paper. It can not only handle an extreme numerical situation with harsh density and viscosity ratio but get reliable results to some extent.

5 Concluding Remark

The numerical simulation and approximation of incompressible and immiscible two-phase flows with matched and variable densities and viscosities is the main topic in this paper. We proposed a brand new scheme for coupled diffuse interface model system with matched and variable densities and viscosities that satisfies the mass conservation and admits an energy law. Plenty of numerical experiments are carried out to illustrate the validation, accuracy compared with sharp interface method by the benchmark problem and to test the robustness of the new scheme for some extreme cases as well.