A fractional step lattice Boltzmann model for two phase flows with large density differences

TL;DR

This paper introduces a fractional step lattice Boltzmann model for accurately simulating two-phase flows with large density differences, improving interface tracking and numerical stability.

Contribution

The proposed method combines fractional step schemes with correction terms and advanced collision models to enhance accuracy and stability in high density ratio two-phase flow simulations.

Findings

Successfully simulates flows with density ratios up to 1000.

Accurately predicts interface deformation and flow dynamics.

Validated against analytical solutions and literature data.

Abstract

In this paper, a fractional step lattice Boltzmann method is proposed to model two-phase flows with large density differences by solving Cahn-Hilliard phase-field equation and the incompressible Navier-Stokes equations.In order to maintain a hyperbolic tangent property of the interface profile and conserve the volume, an interfacial profile correction term and a flux correction term are added into the original Cahn-Hilliard equation respectively. By using a fractional step scheme, the modified Cahn-Hilliard equation is split into two sub-equations. One is solved in the framework of lattice Boltzmann equation method. The other is solved by the finite difference method. Compared with the previous lattice Boltzmann methods, the proposed method is able to maintain the order parameter within a physically meaningful range, which is conductive to track the interface accurately. In addition, the multi-relaxation-time collision model and a high-order compact selective filter operation are employed to enhance the numerical stability. The proposed method can simulate two-phase fluid flows with the density ratio up to $1000$. In order to validate the accuracy and capability of the method, several benchmark problems, including single vortex deform of a circle, translation of a drop, Laplace-Young law, capillary wave and rising bubble with large density ratios, are presented. The results are in good agreement with the analytical solutions and the data in the literature for the investigated benchmarks.

Tables

Table 1. Table 1: The Relative maximum errors with various λ 𝜆 \lambda .

$\lambda$	$U=0.01$				$U=0.05$
	$\sigma =0.01$	$\sigma =0.001$	$\sigma =0.0001$		$\sigma =0.01$	$\sigma =0.001$	$\sigma =0.0001$
$0.0$	0.06314	0.13071	0.20146		0.10459	0.19082	0.21659
$0.001$	0.01638	0.01384	0.01708		0.06907	0.09880	0.10610
$0.005$	0.00002	0	0.00254		0.01356	0.01665	0.01805
$0.01$	0	0.00004	0.00168		0.00005	0.00058	0.00075
$0.05$	0	–	–		0	0.00088	0.00385

Table 2. Table 2: The Relative minimum errors with various λ 𝜆 \lambda .

$\lambda$	$U=0.01$				$U=0.05$
	$\sigma =0.01$	$\sigma =0.001$	$\sigma =0.0001$		$\sigma =0.01$	$\sigma =0.001$	$\sigma =0.0001$
$0.0$	0.04495	0.10503	0.15960		0.08267	0.15024	0.17401
$0.001$	0.01276	0.01210	0.01552		0.05657	0.08323	0.09033
$0.005$	0.00002	0	0.00274		0.01183	0.01506	0.01631
$0.01$	0	0.00006	0.00166		0.00005	0.00053	0.00070
$0.05$	0	–	–		0.00002	0.00072	0.03209

Table 3. Table 3: Relative errors in the volume of a droplet with different radii at t = 2 × 10 5 ​ δ ​ t 𝑡 2 superscript 10 5 𝛿 𝑡 t=2\times 10^{5}\delta t .

Raidus	10	20	30	40
Original CH	3.9344	0.6426	0.2848	0.1596
Present	0	0	0	0

Table 4. Table 4: Physical parameters and dimensionless numbers

Test	Eo	Re	$U_g$	${\rho }_1$	${\rho }_2$	${\mu }_1$	${\mu }_2$
1	10	35	0.0024	1000	100	10	1
2	125	35	0.0019	1000	1	10	0.1

	SF-5	SF-7	SF-9	SF-11
$d_0$	$6/16$	$5/16$	0.243527493120	0.215044884112
$d_1$	$-4/16$	$-15/64$	-0.204788880640	-0.187772883589
$d_2$	$1/16$	$3/32$	0.120007591680	0.123755948787
$d_3$		$-1/64$	-0.045211119360	-0.059227575576
$d_4$			0.008228661760	0.018721609157
$d_5$				-0.002999540835

Equations

$${\partial }_t\varphi +\nabla \cdot (\varphi \mathbf{u})=M{\nabla }^2{\mu }_{\varphi },$$

$$\Psi (\varphi )=\beta ({\varphi }^2-1{)}^2+\frac{\kappa }{2}{\left|\nabla \varphi \right|}^2,$$

$${\mu }_{\varphi }\equiv \frac{\delta \Psi }{\delta \varphi }=4\beta ({\varphi }^3-\varphi )-\kappa {\nabla }^2\varphi .$$

$${\varphi }^{eq}(\mathbf{r})=\tanh \left(\frac{2\mathbf{r}}{W}\right),$$

$${\partial }_t\varphi +\nabla \cdot (\varphi \mathbf{u})=M{\nabla }^2{\mu }_{\varphi }+\nabla \cdot {\mathbf{J}}_{\varphi }+\nabla \cdot \mathbf{q},$$

$${\partial }_t\varphi +\nabla \cdot (\varphi \mathbf{u})$$

$${\partial }_t\varphi$$

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

$${\partial }_t(\rho \mathbf{u})+\nabla \cdot (\rho \mathbf{u}\mathbf{u})$$

$$\rho ={\rho }_1\frac{1+\varphi }{2}+{\rho }_2\frac{1-\varphi }{2},$$

$${\mathbf{F}}_s=\kappa \nabla \cdot (|\nabla \varphi {|}^2-\nabla \varphi \otimes \nabla \varphi ).$$

$${\partial }_t\mathbf{u}+\nabla \cdot (\mathbf{u}\mathbf{u})=-\nabla p+\nabla \cdot 2\nu D(\mathbf{u})+\frac{{\mathbf{F}}_p+{\mathbf{F}}_s+{\mathbf{F}}_{\nu }+{\mathbf{F}}_m+{\mathbf{F}}_g}{\rho },$$

$$f_i(\mathbf{x}+{\mathbf{c}}_{\mathbf{i}}\delta t,t+\delta t)-f_i(\mathbf{x},t)=-{\Lambda }_{ij}^f(f_j(\mathbf{x},t)-f_j^{eq}(\mathbf{x},t))+\left(\mathbf{I}-\frac{{\Lambda }_{ij}^f}{2}\right)F_j\delta t,$$

$$f_i^{eq}=\{\begin{array}{rlrl}& \varphi +({\omega }_0-1)\eta \mu -{\omega }_0\varphi \frac{{\mathbf{u}}^2}{2c_s^2},& & \text{i=0}\\ & {\omega }_i\eta \mu +{\omega }_i\varphi \left(\frac{{\mathbf{c}}_i\cdot \mathbf{u}}{c_s^2}+\frac{({\mathbf{c}}_i\cdot \mathbf{u}{)}^2}{2c_s^4}-\frac{{\mathbf{u}}^2}{2c_s^2}\right),& & \text{i=0}\end{array}$$

$${\mathbf{c}}_i=\{\begin{array}{rlrl}& (0,0)e,& & \text{i=0}\\ & (\cos [(i-1)\pi /2],\sin [(i-1)\pi /2])e,& & \text{i=1,…,4}\\ & \sqrt{2}(\cos [(i-5)\pi /2+\pi /4],\sin [(i-5)\pi /2+\pi /4])e,& & \text{i=5,…,8},\end{array}$$

$${\Lambda }_{ij}^f={\text{M}}^{-1}{\text{S}}^f\text{M},$$

$$\bm{\text{M}}=\left[\begin{array}[]{ccccccccc}1&1&1&1&1&1&1&1&1\\ -4&-1&-1&-1&-1&2&2&2&2\\ 4&-2&-2&-2&-2&1&1&1&1\\ 0&1&0&-1&0&1&-1&-1&1\\ 0&-2&0&2&0&1&-1&-1&1\\ 0&0&1&0&-1&1&1&-1&-1\\ 0&0&-2&0&2&1&1&-1&-1\\ 0&1&-1&1&-1&0&0&0&0\\ 0&0&0&0&0&1&-1&1&-1\end{array}\right],$$

$${\text{S}}^f=\text{diag}(s_0^f,s_1^f,s_2^f,s_3^f,s_4^f,s_5^f,s_6^f,s_7^f,s_8^f),$$

$$F_i={\omega }_i\frac{\varphi {\mathbf{c}}_i\cdot \left({\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}\right)}{c_s^2}.$$

$$m=\eta c_s^2({\tau }_f-0.5)\delta t.$$

$$\bar{\varphi }=\sum _i{f}_i.$$

$${\varphi }^{(1)}$$

$${\varphi }^{(2)}$$

$$\varphi$$

$$\frac{\partial {\varphi }^{\ast }}{\partial x}$$

$$\frac{\partial {\varphi }^{\ast }}{\partial y}$$

$$g_i(\mathbf{x}+{\mathbf{c}}_{\mathbf{i}}\delta t,t+\delta t)-g_i(\mathbf{x},t)=-{\Lambda }_{ij}^g(g_j(\mathbf{x},t)-g_j^{eq}(\mathbf{x},t))+\left(\mathbf{I}-\frac{{\Lambda }_{ij}^g}{2}\right)G_j{\delta }_t,$$

$$g_i^{eq}={\omega }_i\left(\frac{p}{c_s^2}+\frac{{\mathbf{c}}_i\cdot \mathbf{u}}{c_s^2}+\frac{({\mathbf{c}}_i\cdot \mathbf{u}{)}^2}{2c_s^4}-\frac{{\mathbf{u}}^2}{2c_s^2}\right).$$

$$G_i={\omega }_i\frac{1}{\rho }\left(\frac{{\mathbf{c}}_i\cdot \mathbf{F}}{c_s^2}+\frac{\mathbf{u}(\mathbf{F}-\nabla p):({\mathbf{c}}_i{\mathbf{c}}_i-c_s^2\mathbf{I})}{c_s^4}\right),$$

$$p=c_s^2\sum _i{g}_i,$$

$$\mathbf{u}=\sum _i{g}_i{\mathbf{c}}_i+\frac{\mathbf{F}}{2\rho }\delta t.$$

$$\nabla \mathbf{u}+\nabla {\mathbf{u}}^T=\frac{(s_7^g-s_1^g)(\mathbf{F}\cdot \mathbf{u})\mathbf{I}-s_7^g(\mathbf{F}\mathbf{u}+\mathbf{u}\mathbf{F})}{2\rho c_s^2}-\frac{{\sum }_i{\mathbf{c}}_i{\mathbf{c}}_i({\mathbf{M}}^{-1}{\mathbf{S}}^g\mathbf{M}{)}_{ij}(f_j-f_j^{eq})}{c_s^2\delta t}.$$

$$\mu ={\mu }_1\frac{1+\varphi }{2}+{\mu }_2\frac{1-\varphi }{2},$$

Full text

††preprint: APS/123-QED

I Introduction

Over the years, the phase-field model has emerged as a reliable and versatile method for multiphase flows anderson1998diffuse ; ding2007diffuse ; lee2010lattice ; yue2004diffuse ; jacqmin1999calculation ; li2016multi . In this framework, the bulk fluid is represented as an order parameter with constant values. The interface between two phase fluids is described as a finite transition region where phase variable can continuously vary between their values in the bulk fluids. This leads to the two main advantages of diffuse-interface model. One advantage is natural capability of capturing complex topological changes. Another feature is its rigorous thermodynamic basis, such as, most models satisfy a nonlinear stability relationship by using a non-convex free-energy function. The evolution of the order variable is governed by the phase-field equations, including Cahn-Hilliard (CH) and Allen Cahn (AC) equations. To model two-phase flows, the phase field equation should couple with a fluid flow governed by the modified Navier-Stokes (NS) equations with a surface force, leading to the so-called Navier-Stokes/Cahn-Hilliard (NSCH) or Navier-Stokes/Allen-Cahn (NSAC) system. Based on the phase-field theory for multiphase flows, a great number of numerical methods have been developed jacqmin1999calculation ; yue2004diffuse ; fakhari2010phase ; zu2013phase ; shen2015decoupled . Among these methods, the phase-field-based lattice Boltzmann equation (LBE) method has received extensive attention because of its some features chen1998lattice ; lee2005stable ; fakhari2010phase ; zu2013phase ; liang2014phase .

The LBE method is based on simplified mesoscopic kinetic equation with discrete-velocity distribution function instead of discretizations of the macroscopic governing equations. However, the macroscopic governing equations can be recovered from the lattice Boltmzann equation through the Champan-Enskog expansion only if the equilibrium function and the source term are designed properly. At present, most of the LBE methods for two phase flows are developed based on the NSCH system. To correctly recover the desired NSCH equations, great efforts have been made zheng2005lattice ; lee2005stable ; zu2013phase ; liang2014phase ; inamuro2016improved . For example, Zheng et al.zheng2005lattice added a spatial difference term of the distribution function as the source term to recovery the CH equation correctly. Liang et al.liang2014phase proposed a time-derivative term related to the order parameter and the velocity in the source term to ensure the correct CH equation. Li et al. li2012additional presented an additional interfacial force to recover the target momentum equation correctly by considering the continuity equation being nonzero near the interface. However, it is still not satisfactory for the capability of capturing the interfaces and modeling two-phase flows with large density contrasts. In particular, it is a challenging task for the NSCH system because the CH equation is a fourth-order partial differential equation and the system of NSCH equations is very stiff. To model large-density-ratio flows, several approaches have been developed. For example, Wang et al. wang2015lattice proposed a LBE flux solver for two-phase flows with large density ratios, in which the CH equation is solved by a stable high-order WENO difference scheme. Lee et al. lee2005stable proposed a stable mixing difference scheme for the force terms in their model to achieve large density ratios. Inamuro et al inamuro2004lattice proposed a free-energy method for large density ratio flow, in which an additional pressure poisson equation is solved iteratively to enforce the divergence-free condition. Kim et al. kim2015lattice employed a filtering operation for the pressure field to enhance the numerical stability in addition to improving the divergence-free velocity by reducing the incompressibility error. Although these methods can simulate the large-density-ratio flows, another important problem has not received much attention. The value of the order parameter defined as the difference between the volume fractions of two phases should theoretically be $[-1,1]$. However, the numerical solutions of the CH equation have maxima and minima values that are larger than $\phi=1$ and smaller than $\phi=-1$, respectively, due to flow and diffusion. These extreme values can give rise to nonphysical fluid characteristics. In particular, for multiphase flows with large density ratios, this phenomenon is more pronounced leading to numerical instability. In order to address this problem, a cutoff technique for the order parameter is widely used once the order parameter exceeds the theoretical maximum and minimum values inamuro2004lattice ; ren2016improved ; shen2015decoupled . However, this method violates the mass conservation. Recently, Li et al. li2016phase proposed a modified CH equation with an interfacial profile correction term based on the equilibrium profile obtained in the thermodynamically derived phase-field model. The added term is used to enforce the phase-field profile to be a hyperbolic tangent profile. As a result, the values of order parameter in the bulk of each component are almost $\pm 1$. Based on the work of Li et al., Zhang and Ye et al zhang2017flux further proposed a flux-corrected term by considering removing the bulk diffusion motion which may contribute to the loss of volume. And the modified CH equation performed a good performance in terms of eliminating the shrinkage and coarsening effects. Although the hydrodynamic equation is not considered in their work, this modified CH equation is also potential for studying two-phase fluid flows with large density ratios.

In addition, the Allen Cahn phase field equation is also proposed for capturing the interface of two phases chiu2011conservative ; sun2007sharp ; geier2015conservative . Compared with the CH equation with a fourth-order derivative, the AC equation is easier to be solved because only second-order partial differential equation are involved. By recovering the AC and NS equations though the Champman-Enskog analysis, several LBE methods of multiphase flows with large density ratios have been proposed liang2018phase ; fakhari2017diffuse ; joshi2018positivity .

In this work, our objective is to develop a LBE method to model two-phase flows with large density ratios based on the NSCH system. To do this, the modified CH equation mentioned above with a convection term is employed instead of the original convection CH equation. The resulting equation is solved by a fractional step scheme efficiently. To improve the numerical stability, the multi-relaxation-time (MRT) collision model is employed and the source term distribution function is redesigned to recover the governing equations correctly. In addition, a filtering operation is applied to the pressure field and velocity field as well.

The rest of the present paper is organized as follows. In section 2, we describe the modified phase-field equation and Navier-Stokes equations of two incompressible fluids with different densities. In section 3, the MRT-LBE model is constructed to recover the macroscopic governing equations correctly. The proposed model will be verified in section 4. Finally, a brief summary is given in Section 5.

II Governing equations

II.1 The modified Cahn-Hilliard equation

In the diffuse interface model, the evolution of an order parameter $\phi$ is generally governed by a convective CH equation,

[TABLE]

where $\bm{u}$ is the velocity, $M$ is the mobility and $\mu_{\phi}$ is the chemical potential defined below. The CH equation is derived by minimizing the free energy of the system, which can be described by the Helmholtz free energy functional abels2012thermodynamically ,

[TABLE]

where $\beta$ and $\kappa$ are constant that depend on the surface tension $\sigma$ and the interface thickness W, i.e., $\kappa=\frac{3}{8}\sigma W$ and $\beta=\frac{3\sigma}{4W}$. A variational procedure applied to the total free energy will yield the chemical potential jacqmin2000contact ; yan2007lattice ,

[TABLE]

Consider a one-dimensional interface at the equilibrium condition, one can obtain the equilibrium profile for the order parameter,

[TABLE]

where $\bm{r}$ is the coordinate normal to the interface. Based on Eq. (4), the interface profile between the two phases should be a hyperbolic tangent profile at equilibrium. Motivated by this idea, Li et al. li2016phase proposed an interface correction term to maintain the hyperbolic tangent profile of the interface li2016phase . More recently, a flux correction term is developed to maintain conserve the volume zhang2017flux . By considering these two correction terms, the original convection CH equation can be modified to the following form,

[TABLE]

where $\bm{J}_{\phi}=\left(|\nabla\phi|-\frac{2(1-\phi^{2})}{W}\right)\bm{n}$ is used to enforce the interface profile to be $\tanh\left(\frac{2r}{W}\right)$, $\bm{q}=(-M\nabla\mu_{\phi}\cdot\bm{n})\bm{n}$ is used to remove the bulk diffusion, and $\bm{n}=\frac{\nabla\phi}{|\nabla\phi|}$ is an outward pointing normal vector. Instead of directly solving Eq (5) at each time step, a fractional step scheme is employed, which consists of the following two sequential steps:

[TABLE]

where $\bm{J}_{d}=\bm{J}_{\phi}+\bm{q}$. Eq (6) is solved by the LBE method and Eq (7) can be efficiently solved by the finite difference method, which will be presented below.

II.2 The incompressible Navier-Stokes equations

We here consider two-phase flows of immiscible incompressible fluids with different densities. The governing equations in the conservation form are as follows abels2017diffuse ; abels2012thermodynamically ,

[TABLE]

with

[TABLE]

where $\mu$ is the dynamic viscosity, $D(\bm{u})=\frac{1}{2}(\nabla\bm{u}+\nabla\bm{u}^{T})$ denotes the rate of deformation tensor, $\bm{I}$ is Kronecker delta, $\rho_{1}$ and $\rho_{2}$ are the component densities, $p_{h}$ is the pressure, $\bm{J}=\frac{\partial\rho}{\partial\phi}M\nabla\mu_{\phi}$ is the flux due to the diffusion, $\bm{F}_{s}=-\kappa\nabla\cdot(\nabla\phi\otimes\nabla\phi)$ describes the capillary force exerted to the fluids by the interface, $\bm{F}_{g}=\rho\bm{g}$ is the gravity force with $\bm{g}$ representing the gravitational acceleration. The above system is thermodynamically consistent and satisfies an energy dissipation law. In addition to the surface force above, several forms of the surface force have been used in the literature kim2005continuous . In this study, we employ the following formulation starovoitov1994model ,

[TABLE]

Then, the pressure is redefined as $p^{\prime}=p_{h}+\kappa|\nabla\phi|^{2}$.

The momentum equaiton can also be expressed as the following equivalent nonconservative form

[TABLE]

where $\nu=\mu/\rho$ is the kinematic viscosity, $\bm{F}_{p}=-p\nabla\rho$, $\bm{F}_{\nu}=\nu(\nabla\bm{u}+\nabla\bm{u}^{T})\cdot\nabla\rho$ , $\bm{F}_{m}=\bm{J}\cdot\nabla\bm{u}$ and $p=p^{\prime}/\rho$.

III Phase-field-Based Lattice Boltzmann method for incompressible two-phase flows

III.1 The MRT-LBE model for the Cahn-Hilliard equation

The standard form of discrete Boltzmann equation with MRT collision model for the CH equation can be written as

[TABLE]

where $f_{i}(\bm{x},t)$ is the distribution functions for the order parameter field with a discrete velocity $\bm{c}_{i}$ at point $\bm{x}$ and time $t$, $f_{i}^{eq}(\bm{x},t)$ is the corresponding equilibrium state, $\delta t$ is the time step, $\Lambda_{ij}^{f}$ is the collision operate defined below, $\bm{I}$ is the unit matrix and ${F}_{j}$ is the source term. The equilibrium distribution functions $f_{i}^{eq}(\bm{x},t)$ is defined as zhang2018high

[TABLE]

where $\omega_{i}$ is the weighting coefficient, $c_{s}$ is the sound speed and $\eta$ is an adjustable parameter for a given mobility. Because of only two-dimensional simulations being considered, the two-dimensional nine-velocity (D2Q9) structures is used to solve flow fields. The discrete particle velocity $\bm{c}_{i}$ in D2Q9 can be defined as

[TABLE]

where $e=\delta x/\delta t=\delta y/\delta t$ with $\delta x$ and $\delta y$ representing the lattice spacing in the x-direction and y-direction respectively. In this study, both $\delta x$ ($\delta y$) and $\delta t$ are fixed at unity. The corresponding weighting coefficient in each direction is given by $\omega_{0}=4/9$, $\omega_{1,\ldots,4}=1/9$, $\omega_{1,\ldots,8}=1/36$ and $c_{s}=1/\sqrt{3}$.

In the MRT-LBE method, the collision with the forcing term effect is executed in momentum space by multiplying a transformation matrix M while the propagation of the MRT model remains in velocity space. Thus, the collision operate can be written as,

[TABLE]

with

[TABLE]

Here $\bm{\text{S}}^{f}$ is a diagonal relaxation matrix,

[TABLE]

where each element represents the inverse of a relaxation time for the distribution function $f_{i}$ as it is relaxed to the equilibrium distribution function in moment space. For the stability, $0\leq s_{i}^{f}\leq 2$ must be satisfied. The parameters $s_{3}^{f}$ and $s_{5}^{f}$ are related to the relaxation time $\tau_{f}$ in the single-relaxation-time (SRT) model, i.e., $s_{3}^{f}=s_{5}^{f}=1/\tau_{f}$. If all elements equal to each other, the MRT model will reduce to the SRT model. The source term $F_{i}$ in Eq.(12) is given by zhang2018high

[TABLE]

Though the Chapman-Enskog analysis, the mobility is defined as

[TABLE]

The order parameter can be calculated by

[TABLE]

In general, the values of $\bar{\phi}$ not strictly attain in $[-1,1]$. To maintain the physical meaningful values of the order parameter, Eq. (7) is solved by the finite difference method thereafter. To be specific, a third-order total variation diminishing Runge-Kutta scheme is used for the temporal discretization,

[TABLE]

where $L(\phi)=\nabla\cdot\bm{J}_{d}$, $\phi^{(1)}$ and $\phi^{(2)}$ are the intermediate values of the order parameter. Here, the four-order center difference scheme is applied to discretize the right hand terms,

[TABLE]

where $\phi^{*}$ represents $\bar{\phi},\phi^{(1)}$ and $\phi^{(2)}$.

III.2 The MRT-LBE model for the incompressible Navier-Stokes equations

The discrete Boltzmann equation with MRT collision model for hydrodynamic equations can be written as he1999lattice

[TABLE]

where $g_{i}(\bm{x},t)$ is a pressure distribution function with a discrete velocity $\bm{c}_{i}$ at point $\bm{x}$ and time $t$, $\Lambda_{ij}^{g}$ is the collision operator for the pressure distribution function, $g_{i}^{eq}(\bm{x},t)$ is the corresponding equilibrium distribution function and $G_{j}$ is the source term. Similar to Eq. (15), $\Lambda_{ij}^{g}$ can be expressed as $\bm{\text{M}}^{-1}\bm{\text{S}}^{g}\bm{\text{M}}$, where $\bm{\text{S}}^{g}=\text{diag}(s_{0}^{g},s_{1}^{g},s_{2}^{g},s_{3}^{g},s_{4}^{g},s_{5}^{g},s_{6}^{g},s_{7}^{g},s_{8}^{g})$. Here $s_{7}^{g}$ and $s_{8}^{g}$ is related to the relaxation time $\tau_{g}$ in the SRT version, i.e., $s_{7}^{g}=s_{8}^{g}=1/\tau_{g}$. The equilibrium distribution function can be defined as

[TABLE]

Note that the first moment of the equilibrium distribution functions above is based on the fluid velocity $\bm{u}$, not the momentum $\rho\bm{u}$. Compared with the momentum-based formulation, the velocity-based formulation can remove differentiation-by-part error for the divergence-free condition kim2015lattice ; zu2013phase . To recover the momentum equation exactly though the Champan-Enskog expansion, the source term $G_{i}$ in Eq. (25) is defined as guo2002discrete ,

[TABLE]

where $\bm{F}=\bm{F}_{p}+\bm{F}_{s}+\bm{F}_{\nu}+\bm{F}_{m}+\bm{F}_{g}$. Taking the zero- and first order moments of the distribution function $g_{i}$, the macroscopic quantities $\bm{u}$ and $p$ are calculated by

[TABLE]

[TABLE]

Note that the deviatoric stress tensor in $\bm{F}_{\nu}$ can be obtained in terms of the distribution function in addition to using the finite difference calculation. For the MRT model, $\nabla\bm{u}+\nabla\bm{u}^{T}$ can be calculated by

[TABLE]

However, it is clear that obtaining the velocity in Eq. (29) requires implicit calculation. Instead of implicit calculation, this equation is solved iteratively until $\sum|\bm{u}^{n+1}-\bm{u}^{n}|/\sum|\bm{u}^{n}|<\zeta$ is satisfied in the whole domain, where $\zeta$ is a given error threshold. And we found that only few steps are required to maintain numerical stability.

The viscosity $\mu$ is given by

[TABLE]

where $\mu_{1}$ and $\mu_{2}$ are the viscosities of the two fluids, respectively. Once the viscosity of the fluid is obtained, the relaxation time is computed by $\tau_{g}=\mu/(\rho c_{s}^{2}\delta t)+0.5$. Through the Chapman-Enskog analysis, the governing equation in Eq. (11) can be recovered correctly. (see Appendix. A for details.)

The force terms in Eq. (29) involve the first-order derivatives and the Laplacian, which need to be carefully discretized. Like most used, the following isotropic discretization is adopted,

[TABLE]

where $\Phi$ is an arbitrary function. In addition, the finite-difference schemes may generate grid-to-grid oscillations leading to numerical instability. A spatial selective filter is therefore used to eliminate grid-to-grid oscillations tam1993dispersion ; ricot2009lattice . After each time step, the variable $\bm{\chi}=(p,\bm{u})$ is filtered by the following high-order selective filter chao2011filter ; berland2007high ,

[TABLE]

where $P$ and $Q$ are the number of points of the stencil, $D$ is the space dimension, $a$ is a constant between 0 and 1, which describes the intensity of the filter and $d_{m}=d_{-m}$ is a filter coefficient. The filter coefficients are given in Appendix B.

IV Numerical Results and discussion

In this section, we will validate the proposed MRT-LBE method through several numerical benchmarks, including single vortex deformation of a circle, translation of a drop, Laplace-Young law, capillary wave and a rising bubble with large density ratio. The results of each test will compared with the analytical solutions and the data in existing literature.

IV.1 Single vortex deformation of a circular droplet

In order to validate the proposed method the performance of tracking the interface, the problem of Rider and kothe (1995) is simulated in a square domain $L\times L$ balcazar2014finite . A circular droplet with radius $R$ initially is placed at $(0.5L,0.75L)$. The flow field is given as balcazar2014finite ; liang2014phase ,

[TABLE]

where $T=nL/U_{0}$ with $n$ being a positive integer. Based on Eq. (34), the circular droplet will undergo the largest deformation at $T$ and come back to its initial position at $2T$. In the simulations, the parameters are set as $L=400,\sigma=0.001,W=4,\tau_{g}=0.5+\sqrt{3}/6,\eta=2,U_{0}=0.02,\lambda=0.05$. The evolution of the interface advection and deformation obtained by the present method is shown in Fig. 1d. It can be found that the interface is stretched into a very long filament that spirals around the center of the domain, which compares very well with that obtained using other high-order methods balcazar2014finite ; sheu2009development . In particular, the values of the order parameter stay in the physically reasonable interval [-1,1] in the whole process. Fig. 2 shows the result of both the present method and the original CH method at $T$. It can be seen that the interface shapes predicted by both methods are significantly different at the tail. However, the results of the present method are consistent with those in desjardins2008accurate . This is because the present method attempts to maintain the volume conservation, which leads to the formation of drops of the size of a few grids. Fig. 3 shows the initial and the final profile of the interface advected to periods $2T$ with $n=2,4$. From Fig. 3, we can see that the shapes of the droplet predicted by both methods can return to its original position at $2T$. However, the interface profile predicted by the present method returns to its original location of the interface with greater accuracy than the one of the original CH method. The wiggles appearing are mainly due to the stringent velocity field and rapid shift in its direction at $T/2$ and $3/2T$. Finally, we compare the evolution of the normalized area of the circle over time for both methods, as shown in Fig. 4. At $2T$, most of the area of the droplet in the present method can be recovered and the maximum volume loss of the circular droplet is $0.346\%$. In comparison, the area of the droplet in the original CH method has significant volume loss and the maximum volume loss of the circular droplet is $9.26\%$.

IV.2 Effect of parameter $\lambda$

To investigate the effect of parameter $\lambda$ in Eq. (5) on the interface between two phases, we consider a translation of a circular drop of radius $R$ with a constant velocity on the domain $L\times 3L$. The initial shape is set as

[TABLE]

The exact solution at time t can be obtained by $\phi(x,y)=\tanh(\frac{2(R-\sqrt{(x-50-ut)^{2}+(y-50)^{2}})}{W})$. In the simulations, the velocity is set as $\bm{u}=(u,v)=(0.01,0)$. The other parameters are given by $L=100$, $\tau_{g}=0.5+\sqrt{3}/6$, $W=4$, $R=26$. We run the simulation up to $t=4000\delta t$. Comparisons between the numerical results of the modified CH equation with different values of $\lambda$ are shown in Fig. 5.

It can be found that the values of $\phi$ obtained by the original CH method ($\lambda=0$) are slightly smaller than $-1$ on the left and larger than $1$ on the right. By contrast, the interfaces shape obtained by the present model with $\lambda=0.05$ are in excellent agreement with the exact solutions. In addition, the location of the interface predicted by the original CH method deviates from that of the exact solution while the present method can track the interface location correctly.

To furthermore quantitatively investigate the effect of $\lambda$ on the results, we define the relative maximum and minimum errors of the order parameter as follows

[TABLE]

Then, we repeat the above test with $\lambda=0,0.001,0.005,0.01,0.05$. The relative maximum and minimum errors of the order parameter are presented in Tables. 1 and  2. It can be found that the relative errors decrease significantly as $\lambda$ increases, which means that the numerical values of the order parameter are gradually adjusted to $-1$ or $1$. However, for larger $\lambda$, relative maximum and minimum errors of the order parameter cannot continue to decrease. A tiny deviation still exists. In addition, a larger $\lambda$ may give rise to numerical instability. Therefore, for the rest of the paper, we will use $\lambda=0.005$ for presentation unless otherwise stated.

IV.3 Laplace-Young law

All of the tests are carried out by a given velocity field and the hydrodynamic equation is not considered. To validate the accuracy of the proposed model coupled with the hydrodynamic equation, we carry out the Laplace-Young law test. Initially, a static circular droplet with radius $R$ is located in the center $(x_{c},y_{c})$ of a domain size $L\times L$. Periodic boundary conditions are applied on all the boundaries, and the order parameter is given by

[TABLE]

Theoretically, in the equilibrium state, the pressure difference ($\Delta p$) between the inside and outside droplet is equivalent to the ratio of the surface tension and the radius of the droplet, i.e., $\Delta P=\sigma/R$. The density and dynamic viscosity of the droplet are fixed as $\rho_{1}=1000,\rho_{2}=1,\mu_{1}=10,\mu_{2}=0.1$. In Eq. (33), the coefficients of ten-order filter (SF-11) are used and the intensity of the filter is fixed at $a=0.5$ unless otherwise stated. And the other parameters are taken as: $L=100,x_{c}=y_{c}=50,\lambda=0.005,\tau_{g}=1/2+\sqrt{3}/6,\eta=2$ and $W=4$. Fig. 6a shows the relationship between the pressure differences and the reciprocal of different radii of the droplet with three surface tensions ($\sigma=0.01,0.005,0.001$) measured after $2\times 10^{5}\delta t$. The pressure differences in Fig. 6a show good agreement with the Laplace-Law, which verifies the accuracy of the proposed model in modeling the interface tension with large density ratios.

In addition, the volume of the droplet with different radius (R=10, 20, 30, 40) are measured to validate the volume-conserved property for the present method. To quantitative analysis, the percentage error in the droplet volume is calculated by $\text{Err}=((\int_{\phi(\bm{x},0)>0}1d\bm{x}-\int_{\phi(\bm{x},t)>0}1d\bm{x})/\int_{\phi(\bm{x},0)>0}1d\bm{x})\times 100$, where $\phi(\bm{x},0)$ and $\phi(\bm{x},t)$ are the order parameter at $t=0$ and $t=2\times 10^{5}\delta t$, respectively. As the numerical instability occurs for the original CH method with large density contrast, the parameters in this case are set as, $L=200,\rho_{1}=10,\rho_{2}=1,\sigma=0.001,\tau_{f}=0.6$. The results are summarized in Table 3. It can be observed that the volume loss of the droplet obtained by the original CH method increases with the decreasing radius of the droplet while the volume of the droplet obtained by the present method conserve volume exactly.

IV.4 Capillary wave

To further validate the present method, a capillary wave problem prosperetti1981motion ; shao2015hybrid is carried out. As shown in Fig. 7, a heavier fluid (labeled as fluid 2) is placed under a lighter fluid (labeled as fluid 1) with a small perturbation $y=H/2+H_{0}\cos(\kappa x)$ on the interface, where $H$ is the height of the computation domain, $H_{0}$ is the initial perturbation amplitude and $\kappa$ is the wave number. Under the influence of surface tension and viscosity, the capillary wave with a decaying amplitude oscillates until the fluid is at rest. According to the work of Prosperetti prosperetti1981motion , the evolution of interface wave amplitude $h(t)$ is given by

[TABLE]

where $\beta$ is defined as $\rho_{1}\rho_{2}/(\rho_{1}+\rho_{2})^{2}$, $\nu$ is the kinematic viscosity, $\omega_{0}^{2}=(\sigma\kappa^{3})/(\rho_{2}+\rho_{1})$ and $\mbox{erfc}(z_{i})$ is the complementary error function of complex $z_{i}$, which is the solution of the following equation

[TABLE]

and $Z_{i}$ is defined as

[TABLE]

For this case, the computation domain is set as $L\times H=1\times 3$, the kinetic viscosity is $0.01$ and the surface tension is $0.02$. The initial amplitude of the perturbation wave is $H_{0}=0.01H$. The wave number is set as $\kappa=2\pi/L$. No-slip boundary conditions are imposed on the top and bottom wall while periodic boundary conditions on the side directions. In the simulations, a grid size of $300\times 900$ is applied to discretize the computational domain. The other parameters are set as $\tau_{g}=0.5+\sqrt{3}/6$, $W=4$, $\lambda=0.005$. In addition, the non-dimensional time is normalized by $T=t\omega_{0}$. Fig. 8 shows the comparison of the numerical results and the theoretical predictions for $\rho_{2}/\rho_{1}=1,100,1000$. It can be seen that the numerical results are in good agreements with the theoretical solution. However, the attenuation of the amplitude in the theoretical is faster than the one predicted from the simulation. This delay has been pointed out in denner2016frequency that the nonlinear effects are negligible due to the linearized Navier-Stokes equations assumption. Overall, this validates the accuracy of the proposed model in capturing the interface of immiscible two-phase fluids.

IV.5 Rising bubble

In this section, we examine the present model by simulating a single bubble rising in a liquid column and compare the results with the data in the literature. For this problem, a numerical benchmark configuration has been proposed for two-dimensional bubble dynamics, including a small density and viscosity ratio case (labeled as test 1) and a large density and viscosity ratio case (labeled as test 2) hysing2009quantitative ; aland2012benchmark . Different numerical approaches have been quantitatively compared, such as Eulerian level set finite element methods, arbitrary Lagrangian-Eulerian moving grid approaches, and diffuse interface methods. The schematic diagram of this problem is depicted in Fig. 9. A bubble of radius $R$ is placed at (0.5m,0.5m) in a rectangular domain $1\text{m}\times 2\text{m}$. The hydrodynamics are governed by two non-dimensional numbers: the Reynolds number (Re) and the Eotvos number (Eo),

[TABLE]

where $g$ is the gravitational acceleration, $U_{g}=\sqrt{g2R}$ is the gravitational velocity. No-slip conditions are applied on the upper and lower walls and periodic conditions are used on the left and right boundaries. In order to quantify the dynamics of the bubble during its course of rise, the following benchmark quantities are measured, including rise velocity and center of mass, which are defined as,

[TABLE]

[TABLE]

where $v$ is the velocity component in the vertical direction.

The fluid and physical parameters for both tests are listed in table  4. In the simulations, a grid size of $240\times 480$ is applied to discretize the flow domain. The gravitational acceleration and surface tension are determined by the non-dimensional numbers.

Fig. 10 compares the bubble shapes predicted by the present method at the final time (t=3s) with the benchmark solutions of Aland and Voigt aland2012benchmark in terms of test 1. It can be observed that the overall interface shapes predicted by the present method agree well with those of published data. Furthermore, Fig. 11 compares the center of mass and the rising velocity for the rising bubble with the benchmark solutions of Aland and Voigt aland2012benchmark . It can be found that the results obtained by the present method are almost in consistent with the results from the literature.

We also carried out the simulation with a density ratio 1000 and viscosity ratio of 100 in terms of test 2. Both the bubble shapes predicted by the present method and those of published data at $t=3s$ are shown in Fig. 12. From Fig. 12, the shape of bubble experiences significant topology change. And there are visible differences between all of the three models. The thin filaments break up for the results of Hysing et al hysing2009quantitative while no break off occurs for both the results of S. Aland and A. Voigit et al.aland2012benchmark and the present method. However, the overall interface shapes and positions are quite similar. This difference may be caused by the grid resolution because a appropriate thickness of the interface is required for the diffuse interface approximations. As the interface thickness tends to zero, the differences between the results of the three methods should be rather small. In addition, the center of mass and rise velocity over time are shown in Fig. 13. Again, good agreements with the published results are demonstrated. Based on the above results, the accuracy and performance of the present method are confirmed.

V CONCLUSIONS

In this paper, a MRT-LBE method for multiphase flow with large density ratios is developed based on the phase-field theory. To accurately capture the interface, the modified CH equation with the interface correction term and the flux correction term is solved by a fractional scheme. In addition, a high-order compact selective filter is applied to macroscopic quantities (pressure and velocity) to improve the numerical stability. Thereafter, several numerical benchmarks for two-phase flows are performed to evaluate the accuracy and stability of the proposed method. First, the results of both the single vortex deform of a circular droplet and the translation of a drop have been compared fairly well with the theoretical interface shape and location. The accuracy of the proposed method for capturing the interface is demonstrated. By coupling the hydrodynamic equations, the present method is further assessed by the Laplace-Young law and the capillary wave problems. The results are in good agreement with the theoretical solutions. Finally, a rising bubble with large density ratios are simulated to examine the performance of the present method, which demonstrates the capability of the present method for problems with complex interface deformation and large density ratios.

ACKNOWLEDGEMENTS

Appendix A CHAPMAN-ENSKOG ANALYSIS OF THE MRT-LBE MODE FOR THE NAVER-STOKES EQUATIONS

To recover the incompressible immiscible two-phase governing equations from Eq. (25), the following multi-scale expansions are introduced,

[TABLE]

where $\epsilon$ is a small expansion parameter. By applying the Taylor expansion to Eq. (25), one can obtain

[TABLE]

where $D_{i}=\partial_{t}+c_{i}\cdot\nabla$. Substituting Eq. (44) into (45) and taking the zeroth, first-, and second-order in $\epsilon$ leads to

[TABLE]

where $D_{1i}=\partial_{t_{1}}+\bm{c}_{i}\nabla$. Multiplying the matrix M on both sides of Eqs. (46-48), one can obtain the equations in the moment space,

[TABLE]

where $\bm{m}_{g}=\bm{\text{M}}g$, $\bm{\tilde{G}}=\bm{\text{M}}G$, $\bm{\tilde{D}}_{1}=\bm{\text{M}}\bm{D}_{1}\bm{\text{M}}^{-1}$. According to Eqs. (26) and (27), the elements of $\bm{m}_{g}$ and $\bm{\tilde{G}}$ can be expressed as,

[TABLE]

[TABLE]

[TABLE]

Substituting Eqs. (52-54) into Eq. (50), the elements of Eq. (50) with $i=0,3,5$ can be written as,

[TABLE]

and the elements of Eq. (50) with $i=1,7,8$ can be written as,

[TABLE]

Similarly, substituting Eqs. (52-54) into Eq. (51), the elements of Eq. (51) with $i=0,3,5$ can be written as,

[TABLE]

[TABLE]

[TABLE]

where $\nu=c_{s}^{2}(\frac{1}{s_{7,8}^{g}}-\frac{1}{2})\delta t$ represents the kinetic viscosity and $\xi=c_{s}^{2}(\frac{1}{s_{1}^{g}}-\frac{1}{2})\delta t$ represents the bulk viscosity. With the help of Eq. (56), Eqs. (58) and (59) can be reduced to

[TABLE]

By assembling Eq. (55) at the $\epsilon$ scale and Eq. (60) at the $\epsilon^{2}$ scale, the recovered macroscopic equations are

[TABLE]

[TABLE]

In the nearly incompressible limit, the time derivative of the pressure is small and the divergence-free condition is approximately satisfied. Thus Eq (25) can exactly recover the hydrodynamic equations.

Appendix B Selective filters

Coefficients of the explicit centered selective filters (SF):