# A ternary phase-field model for two-phase flows in complex geometries

Chengjie Zhan  Zhenhua Chai<EMAIL_ADDRESS>Baochang Shi  School of
Mathematics and Statistics, Huazhong University of Science and Technology,
Wuhan 430074, China  Institute of Interdisciplinary Research for Mathematics
and Applied Science, Huazhong University of Science and Technology, Wuhan
430074, China  Hubei Key Laboratory of Engineering Modeling and Scientific
Computing, Huazhong University of Science and Technology, Wuhan 430074, China

###### Abstract

In this work, a ternary phase-field model for two-phase flows in complex
geometries is proposed. In this model, one of the three components in the
classical ternary Cahn-Hilliard model is considered as the solid phase, and
only one Cahn-Hilliard equation with degenerate mobility needs to be solved
due to the condition of volume conservation, which is consistent with the
standard phase-field model with a single-scalar variable for two-phase flows.
To depict different wetting properties at the complex fluid-solid boundaries,
the spreading parameters in ternary phase-field model are determined based on
the Young’s law, in which the liquid-solid surface tension coefficient is
assumed to be a linear function of gas-liquid surface tension coefficient and
related to the contact angle and the minimum curvature of the solid surface.
In addition, to achieve a high viscosity in the solid phase and preserve the
velocity boundary conditions on the solid surface, the phase-field variable of
the solid phase is also used to derive the modified Navier-Stokes equations.
To test the present model, we further develop a consistent and conservative
Hermite-moment based lattice Boltzmann method where an adjustable scale factor
is introduced to improve the numerical stability, and conduct the numerical
simulations of several benchmark problems. The results illustrate that present
model has the good capability in the study of the two-phase flows in complex
geometries.

###### keywords:

Ternary phase-field model , two-phase flows , complex geometries , lattice
Boltzmann method

##  1 Introduction

The two-phase flows in complex geometries are usually observed in both nature
and engineering, and also become the fundamental problems in the fields of the
enhanced oil recovery [1, 2], the geological CO2 sequestration [3], the
advanced manufacture [4, 5], and so on. Due to the difficulty in exploring the
physical mechanism and the expensive cost of experimental approach, the
numerical simulation has become an effective method in the study of such
complex problems [6]. To describe two immiscible fluid phases, the Cahn-
Hilliard (CH) equation [7], as one of the popular used phase-field model [8,
9], has been widely applied due to its properties of the mass conservation and
thermodynamic consistency [10].

When we consider the fluid flows in complex geometries, how to implement the
boundary conditions is usually a tough issue. To solve this problem, some
approaches have been proposed, including the immersed boundary method [11,
12], the fictitious domain method [13, 14], the volume penalization method
[15, 16], the smoothed profile method [17, 18], and the diffuse domain method
[19, 20]. In these methods, a larger and regular domain containing the
original complex geometries and the fluid-solid boundaries is first
introduced, and simultaneously a sharp/diffuse indicator function at
Cartesian/Lagrangian grids is used to distinguish the interface between fluid
and solid phases. Then the boundary conditions at fluid-solid interface are
treated as the source terms and/or external force terms into the reformulated
governing equations or numerical schemes, and thus the difficulty in directly
treating the complex boundary conditions can be avoided.

For the multiphase problems, the wetting property also has an obvious
influence on the fluid dynamics. However, the application of the wetting
boundary condition usually needs to use the local normal vector of the solid
boundaries, which is difficult to be determined for the complex geometries. Up
to now, some of the above-mentioned methods [21, 22, 23, 24, 25, 26, 27] have
been successfully applied to investigate the multiphase flows with wetting
boundary conditions, where a smooth indicator function for the solid phase is
naturally combined with the diffuse-interface methods for two-phase flows.
However, in these methods, how to reconstruct the governing equations or
numerical schemes is also an important issue, and some additional treatments
are needed.

From another point of view, the two-phase flows in complex geometries can be
seen as the three-component flows in a larger and regular domain, and one of
which is considered as the solid phase. Actually, in the phase-field method,
Boyer and Lapuerta [28] designed a consistent free energy for the ternary CH
model, which can reduce to the standard form for two-phase cases when one
phase disappears. Although there have been many phase-field models for
multiphase flows [29, 30, 31] up to now, we only focus on this ternary-fluid
model in present work. Based on the model in Ref. [28], Yi [32] proposed a
phase-field method for the moving solid object in two-phase flows, where the
rigidity constraint of the solid phase and the no-slip boundary condition on
the solid-fluid interface is imposed by attributing a high viscosity to the
solid phase. This work aims to handle complex two-phase flows with a moving
solid object, but unfortunately, their results show that the solid phase has
some deformation. What is more, to consider a hydrophilic or hydrophobic solid
object, the surface tension coefficients of the gas/liquid/solid interfaces
are artificially given without any discussion. Later, Rohde and Wolff [33]
presented a ternary phase-field model for two-phase flow with precipitation
and dissolution. In their model, when they consider the wetting boundary
condition, the surface tension coefficients related to the solid phase are
assumed to be much larger than the surface tension coefficient between two
fluid phases, such that the mechanical balance condition at the triple points
can be preserved as far as possible. Additionally, to avoid the restriction of
the maximum principle for fourth-order CH equation, the free energy is
designed to enforce a relaxed constraint $\left(-\delta,1+\delta\right)$ for
the order parameter by using an unbounded potential function, instead of the
original range restriction $\left[0,1\right]$. Furthermore, the Navier-Stokes
(NS) equations for the fluid flows in their model are also more reasonable
through considering some conditions of consistency. However, in their work,
the values of the surface tension coefficients are not given, and the capacity
in predicting contact angle has also not been tested. More recently, Panter et
al. [34] used a ternary phase-field model [28] to investigate the capillary
rise within rough structures, in which one of three order parameters is
defined as the rough structures to smooth and curve the solids. However, in
their model, the wetting property of the rough structures seems not to be
considered, and the solid surface at the base of the pillars is still treated
via a boundary condition, rather than a diffuse interface.

In this work, to overcome the problems mentioned in above literature, we
present some improvements and propose a mathematical model for two-phase flows
in complex geometries. In the present model, a specific expression of the
liquid-solid surface tension coefficient is first designed, which is more
accurate in predicting the contact angles. Then a degenerate mobility is used
to overcome the mass leakage among different phases in the CH model [35].
Finally, the consistent and conservative NS equations is adopted to describe
the fluid flows, and the volume fraction of the solid phase is used to achieve
the high viscosity in solid phase and ensure the velocity boundary conditions
on solid interface. To numerically solve the present model, the lattice
Boltzmann method [36, 37, 38] is applied, which is an efficient numerical tool
in the simulation of complex fluid systems, especially for the multiphase
flows due to the features of the easy implementation of boundary conditions
and fully parallel algorithm [37].

The rest of this paper is organized as follows. In Section 2, the original
three-component CH model is first introduced, and then the model for two-phase
flows in complex geometries and the relation of surface tension coefficients
are designed, followed by the consistent and conservative NS equations. In
Section 3, the consistent and conservative LB method for the proposed phase-
field model is developed. In Section 4, the numerical validations are
performed by several benchmark tests, and the problems of displacement in a
complex channel and a porous medium are studied. Finally, some conclusions are
summarized in Section 5.

##  2 Ternary phase-field model for two-phase flows in complex geometries

###  2.1 The Cahn-Hilliard model for three-component flows

The mix free energy for three-component flows can be expressed as [28]

$E=\int_{\Omega}f_{b}\left(\phi_{0},\phi_{1},\phi_{2}\right)+\frac{3D}{8}\sum_{p=0}^{2}\gamma_{p}|\nabla\phi_{p}|^{2}\mathrm{d}\Omega,$
(1)

where $\phi_{p}$ is the phase-field variable of phase $p$ satisfying
$\phi_{0}+\phi_{1}+\phi_{2}=1$, $D$ is the interface thickness, and the
spreading coefficient $\gamma_{p}$ is related to the pairwise symmetric
surface tension coefficient $\sigma_{pq}$,

$\displaystyle\gamma_{0}$ $\displaystyle=\sigma_{01}+\sigma_{02}-\sigma_{12},$
(2) $\displaystyle\gamma_{1}$
$\displaystyle=\sigma_{01}+\sigma_{12}-\sigma_{02},$ $\displaystyle\gamma_{2}$
$\displaystyle=\sigma_{02}+\sigma_{12}-\sigma_{01}.$

$f_{b}\left(\phi_{0},\phi_{1},\phi_{2}\right)$ is the bulk free energy
density, and is given by [28]

$f_{b}\left(\phi_{0},\phi_{1},\phi_{2}\right)=\frac{12}{D}\left[\frac{\gamma_{0}}{2}\phi_{0}^{2}\left(1-\phi_{0}\right)^{2}+\frac{\gamma_{1}}{2}\phi_{1}^{2}\left(1-\phi_{1}\right)^{2}+\frac{\gamma_{2}}{2}\phi_{2}^{2}\left(1-\phi_{2}\right)^{2}\right].$
(3)

The chemical potential of $p$-th fluid can be calculated by the variation of
mix free energy,

$\mu_{p}=\frac{\delta E}{\delta\phi_{p}}=\frac{\partial
f_{b}}{\partial\phi_{p}}-\frac{3D}{4}\gamma_{p}\nabla^{2}\phi_{p}.$ (4)

Then one can obtain the following CH equation through minimizing the mix free
energy,

$\frac{\partial\phi_{p}}{\partial
t}+\nabla\cdot\left(\phi_{p}\mathbf{u}\right)=\nabla\cdot
M_{p}\nabla\bar{\mu}_{p},\quad p=0,1,2,$ (5)

where $\mathbf{u}$ is the fluid velocity, $M_{p}$ is the mobility,
$\bar{\mu}_{p}=\mu_{p}+L$ with $L$ being a Lagrangian multiplier used to
ensure the volume conservation of phase-field variables. To determine the
expression of $L$, one can assume
$M=M_{0}\gamma_{0}=M_{1}\gamma_{1}=M_{2}\gamma_{2}$ and sum all CH equations
with $p=0,1,2$,

$0=\frac{\partial 1}{\partial
t}+\nabla\cdot\left(1\mathbf{u}\right)=\nabla\cdot
M\sum_{p=0}^{2}\left(\frac{\mu_{p}}{\gamma_{p}}+\frac{L}{\gamma_{p}}\right),$
(6)

where the conservative condition of $\phi_{0}+\phi_{1}+\phi_{2}=1$ is used.
From the above equation, $L$ can be given by

$L=-\frac{\gamma^{T}}{3}\sum_{q=0}^{2}\frac{\mu_{q}}{\gamma_{q}},\quad\frac{3}{\gamma^{T}}=\frac{1}{\gamma_{0}}+\frac{1}{\gamma_{1}}+\frac{1}{\gamma_{2}}.$
(7)

With the above expression of $L$, the CH model can be rewritten as

$\frac{\partial\phi_{p}}{\partial
t}+\nabla\cdot\left(\phi_{p}\mathbf{u}\right)=\nabla\cdot\left(\frac{M}{\gamma_{p}}\nabla\bar{\mu}_{p}\right),$
(8a) $\bar{\mu}_{p}=\frac{\gamma^{T}}{3}\sum_{q\neq
p}\frac{1}{\gamma_{q}}\left(\mu_{p}-\mu_{q}\right),\quad p=0,1,2.$ (8b)

###  2.2 The Cahn-Hilliard model for two-phase flow in complex geometries

For the two-phase flows in complex geometries, the system can be viewed as a
three-phase case and be governed by the three-component CH model in Section
2.1 with some extra limitations. In this work, we set phase 0 to be the solid
phase, and thus $\phi_{0}$ becomes a smoothed indicator function to identify
the solid phase ($\phi_{0}=1$) and fluid phase ($\phi_{0}=0$).

To include the effect of the wetting property on the solid surface in a three-
phase system, the Young’s law [39] can be applied, i.e.,
$\sigma_{02}=\sigma_{01}+\sigma_{12}\cos\theta$ where $\theta$ is the contact
angle of fluid phase 1 to the solid phase 0 (see the schematic in Fig. 1).
Additionally, following the work of Rohde and Wolff [33], the mechanical
balance of surface tension coefficients and the angles at the triple points
can be given by

$\frac{\sin\alpha}{\sigma_{12}}=\frac{\sin\beta}{\sigma_{01}}=\frac{\sin\theta}{\sigma_{02}},$
(9)

where $\alpha+\beta+\theta=2\pi$. When the phase 0 is considered as solid
phase, we need to set $\alpha=\pi$, and $\sin\alpha=0$. Thus, the surface
tension coefficients need to satisfy $\sigma_{12}\ll\sigma_{01},\sigma_{02}$
[33]. Actually, under this condition, the terms with large value of
$\gamma_{0}$ in the mix free energy (1) can be seen as the penalties for solid
phase 0, and there should be an optimal penalizing parameter to give the most
reasonable result. However, as far as we know, there are no specific
expressions of the spreading coefficients or fluid/solid surface tension
coefficients at present. For this reason, in this work, we assume
$\sigma_{01}=k\sigma_{12}$ and give the following expression of $k$,

$k=\frac{9}{1+\kappa}\left[1+2\mathrm{sgn}\left(\theta-\frac{\pi}{2}\right)\sin^{2}\left(\frac{\pi}{4}\cos\theta\right)\right],$
(10)

which is a characteristic parameter related to the curvature and the contact
angle of the solid surface. Here $\kappa$ is the minimum dimensionless
curvature of a solid structure and $\mathrm{sgn}\left(x\right)$ is a sign
function. Specially, $\kappa=0$ denotes the case of an ideal solid plate. With
the above formula, the parameter $\gamma_{p}$ can be rewritten as

$\displaystyle\gamma_{0}$
$\displaystyle=\left(2k+\cos\theta-1\right)\sigma_{12},$ (11)
$\displaystyle\gamma_{1}$
$\displaystyle=\left(1-\cos\theta\right)\sigma_{12},$
$\displaystyle\gamma_{2}$
$\displaystyle=\left(1+\cos\theta\right)\sigma_{12}.$

Moreover, to eliminate the defects of the CH equation in terms of the
boundedness and phase mass conservation, a degenerate mobility
$M=M_{cst}+M_{deg}f\left(\bm{\phi}\right)$ is applied in this work, where
$M_{cst}$ and $M_{deg}$ are two positive constants and
$f\left(\bm{\phi}\right)$ is a function with
$\bm{\phi}=\left(\phi_{0},\phi_{1},\phi_{2}\right)$. Here it should be noted
that the form of mobility does not influence the formal derivation of the
present model.

Figure 1: Three-phase contact diagram.

######  Remark 1.

Compared to the original CH model for three-phase flows, the present model for
two-phase flows in complex geometries needs to set the values of $\sigma_{01}$
and $\sigma_{02}$ according to Eq. (10) and Young’s law.

######  Remark 2.

For the two phase flows in stationary geometries, only $\phi_{1}$ governed by
the CH model (8) needs to be solved and $\phi_{2}$ can be calculated by the
conservation of volume, which is consistent with the standard phase-field
equation with a single-scalar variable for two-phase flows.

######  Remark 3.

The phase-field variable $\phi_{0}$ for solid phase is usually not governed by
the CH model (8) with $p=0$, and it does not need to be evolved for stationary
geometries. In addition, the evolution of $\phi_{0}$ needs to be included when
the problems including solid moving, deformation or phase transition are
considered, for instance, the phase-field equations for dendritic growth [40,
41].

###  2.3 The consistent and conservative Navier-Stokes equations for flow
field

To depict the fluid flows, the following consistent and conservative NS
equations are adopted,

$\nabla\cdot\mathbf{u}=0,$ (12a)
$\frac{\partial\left(\rho\mathbf{u}\right)}{\partial
t}+\nabla\cdot\left[\left(\rho\mathbf{u}-\mathbf{S}\right)\mathbf{u}\right]=-\nabla
P+\nabla\cdot\frac{\mu}{1-\phi_{0}}\left[\nabla\mathbf{u}+\left(\nabla\mathbf{u}\right)^{\top}\right]+\mathbf{F}_{s}+\phi_{0}\rho\mathbf{f}+\mathbf{F}_{b},$
(12b)

where $\rho$ is the density, $P$ is the pressure.
$\mu/\left(1-\phi_{0}\right)$ is the effective dynamic viscosity with
$\left(1-\phi_{0}\right)$ being used to give a high viscosity in solid phase
and a normal viscosity in fluid region with $\phi_{0}=0$. $\mathbf{S}$ is the
extra mass flux across the interface of two fluid phases, and is given by
$\mathbf{S}=\sum_{p=0}^{2}\rho_{p}M\nabla\bar{\mu}_{p}/\gamma_{p}$ for the CH
equation (8). $\mathbf{F}_{s}=\sum_{p=0}^{2}\mu_{p}\nabla\phi_{p}$ is the
surface tension force, and $\mathbf{F}_{b}$ represents the body force, such as
gravity. The velocity boundary condition is imposed by the force
$\phi_{0}\rho\mathbf{f}$, and the form of fluid-structure interaction
$\mathbf{f}$ will be stated in the next section.

In addition, the density and dynamic viscosity are usually set as the linear
functions of the phase-field variables,

$\rho=\sum_{p=0}^{2}\rho_{p}\phi_{p},\quad\mu=\sum_{p=0}^{2}\mu_{p}\phi_{p},$
(13)

where $\rho_{p}$ and $\mu_{p}$ are the density and viscosity of the pure
$p$-th fluid.

##  3 Lattice Boltzmann method for the three-component phase-field model

In this section, we will develop the LB method to solve the mathematical model
of three-component flows under the general frameworks [42, 43]. However, it
should be noted that the degenerate mobility in CH equation is a function of
the phase-field variables and may become zero, which may cause the numerical
instability in the standard LB model. To overcome this problem, we rewrite the
CH equation into the following form with a constant diffusion coefficient,

$\frac{\partial\phi_{p}}{\partial
t}+\nabla\cdot\left(\phi_{p}\mathbf{u}\right)=\nabla\cdot\frac{M_{cst}+M_{deg}}{\gamma_{p}}\left[\nabla\bar{\mu}_{p}-M_{deg}\frac{1-f\left(\bm{\phi}\right)}{M_{cst}+M_{deg}}\nabla\bar{\mu}_{p}\right].$
(14)

For the traditional three-phase flows, the so-called _3f_ model composed of
the equations with $p=0$ and 1 needs to be solved, while in the _2f1s_ model
for two-phase flows in complex geometries, only the CH equation with $p=1$
needs to be considered.

In the LB method, the evolution equation for CH equation (14) can be given by

$f_{i}^{p}\left(\mathbf{x}+\mathbf{c}_{i}\Delta t,t+\Delta
t\right)=f_{i}^{p}\left(\mathbf{x},t\right)-\Lambda_{ij}^{p}\left[f_{j}^{p}\left(\mathbf{x},t\right)-f_{j}^{p,eq}\left(\mathbf{x},t\right)\right]+\Delta
t\left(\delta_{ij}-\frac{\Lambda_{ij}^{p}}{2}\right)R_{j}^{p},$ (15)

where $f_{i}^{p}\left(\mathbf{x},t\right)$ is the distribution function of
phase-field variable $\phi_{p}$ in the direction of $i$ ($i=0,1,\cdots,Q-1$
with $Q$ being the number of discrete velocity directions) at position
$\mathbf{x}$ and time $t$. $\mathbf{c}_{i}$ is the discrete velocity, $\Delta
t$ is the time step, and $(\Lambda_{ij}^{p})$ represents the collision matrix.
In addition, the distribution functions $f_{i}^{p,eq}$ and $R_{i}^{p}$ are
designed as

$f_{i}^{p,eq}=\begin{cases}\phi_{p}+\left(\omega_{i,p}-1\right)\eta_{p}\bar{\mu}_{p},&i=0,\\\
\omega_{i,p}\eta_{p}\bar{\mu}_{p}+\omega_{i}\mathbf{c}_{i}\cdot\phi_{p}\mathbf{u}/c_{s,p}^{2},&i\neq
0,\end{cases}\quad
R_{i}^{p}=\omega_{i,p}\left[\frac{\mathbf{c}_{i}\cdot\partial_{t}\left(\phi_{p}\mathbf{u}\right)}{c_{s,p}^{2}}+M_{deg}\frac{1-f\left(\bm{\phi}\right)}{M_{cst}+M_{deg}}\mathbf{c}_{i}\cdot\eta_{p}\nabla\bar{\mu}_{p}\right],$
(16)

where $\omega_{i,p}$ is the weight coefficient, $\eta_{p}$ is an adjustable
parameter, $c_{s,p}^{2}=d_{0}^{p}c^{2}$ is the lattice sound speed with
$c=\Delta x/\Delta t$, $d_{0}^{p}$ is a scale factor and $\Delta x$ is the
lattice spacing. The CH equation (14) can be correctly recovered from above LB
model through the direct Taylor (or Chapman-Enskog) expansion with
$\left(M_{cst}+M_{deg}\right)/\gamma_{p}=\left(1/s_{1}^{p}-0.5\right)\eta_{p}c_{s,p}^{2}\Delta
t$, where $s_{1}^{p}$ is the relaxation parameter related to an eigenvalue of
the collision matrix $(\Lambda_{ij}^{p})$.

We now consider the LB model for NS equations, and compared to our previous
works [44, 45], an adjustable scale factor is introduced to make the present
model more flexible. The evolution equation of the LB model reads

$g_{i}\left(\mathbf{x}+\mathbf{c}_{i}\Delta t,t+\Delta
t\right)=g_{i}\left(\mathbf{x},t\right)-\Lambda_{ij}\left[g_{j}\left(\mathbf{x},t\right)-g_{j}^{eq}\left(\mathbf{x},t\right)\right]+\Delta
t\left(\delta_{ij}-\frac{\Lambda_{ij}}{2}\right)F_{j}\left(\mathbf{x},t\right),$
(17)

where $g_{i}\left(\mathbf{x},t\right)$ is the distribution function for fluid
field, and the corresponding equilibrium distribution function
$g_{i}^{eq}\left(\mathbf{x},t\right)$ is given by

$g_{i}^{eq}=\lambda_{i}+\omega_{i}\frac{\mathbf{c}_{i}\cdot\rho\mathbf{u}}{c_{s}^{2}}+\omega_{i}\left[\frac{\left(c_{i\alpha}c_{i\alpha}-c_{s}^{2}\right)\left(\rho
u_{\alpha}u_{\alpha}-S_{\alpha}u_{\alpha}\right)}{c^{2}c_{s}^{2}-c_{s}^{4}}+\frac{c_{i\alpha}c_{i\bar{\alpha}}}{c_{s}^{4}}\left(\rho
u_{\alpha}u_{\bar{\alpha}}-\frac{S_{\alpha}u_{\bar{\alpha}}+u_{\alpha}S_{\bar{\alpha}}}{2}\right)\right],$
(18)

with $\lambda_{0}=\left(\omega_{0}-1\right)P/c_{s}^{2}+\rho_{0}$, and
$\lambda_{i}=\omega_{i}P/c_{s}^{2}\,(i\neq 0)$. $\omega_{i}$ and $c_{s}$ are
also the weight coefficient and lattice sound speed, which can be different
from those of LB model for CH equation. Here the two terms in the square
bracket correspond to the diagonal and non-diagonal parts of the momentum
flux, respectively. $\alpha=1,2,\cdots,d$, $\alpha<\bar{\alpha}\leq d$ with
$d$ being the dimensionality. The distribution function of force term is
expressed as

$F_{i}=\omega_{i}\left[\mathbf{u}\cdot\nabla\rho+\frac{\mathbf{c}_{i}\cdot\left(\mathbf{F}+\phi_{0}\rho\mathbf{f}\right)}{c_{s}^{2}}+\frac{\left(c_{i\alpha}c_{i\alpha}-c_{s}^{2}\right)M_{\alpha\alpha}^{2F}}{c^{2}c_{s}^{2}-c_{s}^{4}}+\frac{c_{i\alpha}c_{i\bar{\alpha}}M_{\alpha\bar{\alpha}}^{2F}}{c_{s}^{4}}\right],$
(19)

where $\mathbf{F}=\mathbf{F}_{s}+\mathbf{F}_{b}+\mathbf{F}_{c}$,
$\mathbf{F}_{c}=\nabla\cdot\left(\mathbf{S}\mathbf{u}-\mathbf{u}\mathbf{S}\right)/2$
is a corrected force term to exactly recover the consistent momentum equation
(12b), and the second-order tensor $\mathbf{M}^{2F}$ is designed as

$\mathbf{M}^{2F}=\partial_{t}\left(\rho\mathbf{u}\mathbf{u}-\frac{\mathbf{S}\mathbf{u}+\mathbf{u}\mathbf{S}}{2}\right)+c_{s}^{2}\left(\mathbf{u}\nabla\rho+\nabla\rho\mathbf{u}\right)+\left(c^{2}-3c_{s}^{2}\right)\left(\mathbf{u}\cdot\nabla\rho\right)\mathbf{I}.$
(20)

Through some asymptotic analysis methods, the NS equations (12) can be
recovered with the following relation [43],

$\frac{\mu}{\rho\left(1-\phi_{0}\right)}=\left(\frac{1}{s_{2a}}-\frac{1}{2}\right)\frac{c^{2}-c_{s}^{2}}{2}\Delta
t=\left(\frac{1}{s_{2b}}-\frac{1}{2}\right)c_{s}^{2}\Delta t,$ (21)

where $s_{2a}$ and $s_{2b}$ are two relaxation parameters and will be stated
below.

Finally, to update the evolution of the LB method, the macroscopic phase-field
variable, velocity and pressure are computed by

$\phi_{p}=\sum_{i}f_{i},$ (22)
$\rho\mathbf{u}^{*}=\sum_{i}\mathbf{c}_{i}g_{i}+\frac{\Delta
t}{2}\mathbf{F},\quad\mathbf{u}=\mathbf{u}^{*}+\frac{\Delta
t}{2}\phi_{0}\mathbf{f},$ (23)
$p=\frac{c_{s}^{2}}{1-\omega_{0}}\left[\sum_{i\neq
0}g_{i}+\left(\frac{1}{2}+H\right)\Delta
t\mathbf{u}\cdot\nabla\rho-\omega_{0}\frac{\left(\rho\mathbf{u}-\mathbf{S}\right)\cdot\mathbf{u}}{c^{2}-c_{s}^{2}}+\frac{K}{c^{2}}\frac{\Delta
t}{2}\partial_{t}\left(\left(\rho\mathbf{u}-\mathbf{S}\right)\cdot\mathbf{u}\right)\right],$
(24)

where $\mathbf{u}^{*}$ is the velocity without considering the fluid-solid
interaction, $\mathbf{f}$ can be discretized by
$\mathbf{f}=\left(\mathbf{u}_{w}-\mathbf{u}^{*}\right)/\Delta t$ with
$\mathbf{u}_{w}$ being the velocity of solid point. $K=1-d_{0}$ and
$H=K\left[\left(1-d_{0}\right)/s_{0}-\left(2-4d_{0}\right)/s_{2a}+\left(1-3d_{0}\right)/2\right]$
for two-dimensional (2D) problem, while $K=1-2d_{0}$ and
$H=K\left[\left(1-d_{0}\right)/s_{0}-\left(3-7d_{0}\right)/s_{2a}+\left(1-3d_{0}\right)\right]$
for three-dimensional (3D) case, $d_{0}$ is the adjustable scale factor and
$s_{0}$ is the relaxation parameter corresponding to the zero-order moment of
the distribution function. $\partial_{t}\left(\rho\mathbf{uu}\right)$ is
usually approximated by $\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$, other
temporal derivatives are discretized by the explicit Euler scheme, and the
gradient and Laplacian operators are computed by the second-order isotropic
central schemes, as shown in the previous work [44]. Here it is worth noting
that when $d_{0}=1/3$, the present model would reduce to the usual LB model
for the two-phase problems.

##  4 Numerical simulations and discussion

In the framework of LB method, there are some LB models based on different
collision matrices, such as the classical MRT-LB model [46, 47, 48], the
cascaded or central-moment LB model [49, 50], the Hermite moment LB model
[51], the central-Hermite moment LB model [52, 53] and so on. However, under
the unified framework of modeling of the MRT-LB model in Refs. [42, 43],
different forms of the collision matrices can be converted into each other
through a specific lower triangular matrix with the diagonal element of unity.
Therefore, when conducting the collision step of LB method in moment space to
improve the efficiency of the algorithm, one can use the natural moment space,
and multiply a specific lower triangular matrix on the relaxation matrix to
apply the corresponding collision matrix. For example, the collision matrix of
the Hermite moment LB model $\bm{\Lambda}=\mathbf{H}^{-1}\mathbf{KH}$
($\mathbf{H}$ is the Hermite transformation matrix and $\mathbf{K}$ is the
diagonal relaxation matrix) can be written as

$\bm{\Lambda}=\left(\mathbf{NM}\right)^{-1}\mathbf{KNM}=\mathbf{M}_{0}^{-1}\left(\mathbf{C}_{d}^{-1}\mathbf{N}^{-1}\mathbf{KN}\mathbf{C}_{d}\right)\mathbf{M}_{0}=\mathbf{M}_{0}^{-1}\mathbf{S}\mathbf{M}_{0},$
(25)

where $\mathbf{M}=\mathbf{C}_{d}\mathbf{M}_{0}$ is the natural transformation
matrix with $\mathbf{C}_{d}$ being a diagonal matrix formed by the powers of
the lattice speed $c$, $\mathbf{N}$ is the lower triangular matrix associated
with $\mathbf{H}$ and $\mathbf{M}$, and $\mathbf{S}$ is a new lower triangular
relaxation matrix.

In this work, we will apply the above Hermite moment LB model to perform some
2D and 3D simulations, and the specific forms of weight coefficients
$\omega_{i,p}$ and $\omega_{i}$, matrices $\mathbf{K}$, $\mathbf{C}_{d}$,
$\mathbf{M}_{0}$, and $\mathbf{N}$ in D2Q9 and D3Q15 lattice structures are
given in A. In the numerical simulations, the lattice units are used, and are
listed in Table 1. The half-way bounce-back scheme [54] is applied for the no-
flux scalar and no-slip velocity boundary conditions, and the non-equilibrium
extrapolation scheme [55] for the velocity or pressure boundary conditions at
inlet and/or outlet. Additionally, we apply the degenerate mobility by setting
$f\left(\bm{\phi}\right)=\prod_{p=0}^{2}\left(1-\phi_{p}\right)$ for _3f_
model and
$f\left(\bm{\phi}\right)=\phi_{1}\left(1-\phi_{1}\right)\left(1-\phi_{0}\right)$
for _2f1s_ model. In the fluid region with $\phi_{0}=0$, the present
$f\left(\bm{\phi}\right)$ in _2f1s_ model can reduce to the classical
degenerate mobility, which indicates that the boundedness of the solution of
CH equation for two-phase flows can be preserved [56].

Table 1: Units in lattice Boltzmann method Variable |  | Mass |  | Length |  | Time
---|---|---|---|---|---|---
Unit |  | $\mathrm{m}\mathrm{u}$ |  | $\mathrm{l}\mathrm{u}$ |  | $\mathrm{t}\mathrm{s}$

###  4.1 The results of two models for a simple case

As mentioned earlier, there are some differences between the _3f_ model and
_2f1s_ model, i.e., the values of solid/liquid surface tension coefficients,
the number of the CH equation, the expressions of the degenerate mobility and
the effective viscosity, and the first two of them are the key limitations
from _3f_ model to _2f1s_ model. We first show the numerical results of these
two models through performing a simple test in the 2D domain of
$[-L_{0}/2,L_{0}/2]\times[-L_{0}/2,L_{0}/2]$ with $L_{0}=300$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$. The periodic boundary condition is
applied at the left and right boundaries, while the no-flux boundary condition
is imposed on the upper and lower walls. A droplet or a neutral particle with
the radius $R_{0}=50$ $\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$ is placed at
the interface of a binary fluid layer, and the phase-field variables are
initialized by

$\displaystyle\phi_{0}\left(x,y\right)$
$\displaystyle=\frac{1}{2}+\frac{1}{2}\tanh\frac{R_{0}-\sqrt{x^{2}+y^{2}}}{D/2},$
(26) $\displaystyle\phi_{1}\left(x,y\right)$
$\displaystyle=\left[1-\phi_{0}\left(x,y\right)\right]\left(\frac{1}{2}-\frac{1}{2}\tanh\frac{2y}{D}\right).$

When the gravity is not considered, the droplet will form a lens at the
equilibrium state dominated by the surface tension, while the neutral particle
will maintain the initial state.

In our simulations, we set some physical parameters as $M_{cst}=0$,
$M_{deg}=0.1$, $\rho_{0}:\rho_{1}:\rho_{2}=5:10:1$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{l}\mathrm{u}^{3}\mathrm{]}$,
$\mu_{0}:\mu_{1}:\mu_{2}=0.05:0.1:0.01$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{(}\mathrm{l}\mathrm{u}\,\mathrm{t}\mathrm{s}\mathrm{)}\mathrm{]}$,
$D=5$ $\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$ and $\sigma_{12}=0.001$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{t}\mathrm{s}^{2}\mathrm{]}$.
$\sigma_{01}=\sigma_{02}=\sigma_{12}$ is adopted for the _3f_ model, while Eq.
(10) is used for _2f1s_ model. The final shapes of droplet and neutral
particle are plotted in Fig. 2. From this figure, one can observe that these
two models show a good performance for the multiphase problems.

(a) Droplet based on _3f_ model

(b) Neutral particle based on _2f1s_ model

Figure 2: The equilibrium shapes of the droplet and the neutral particle.
Figure 3: The equilibrium shapes of the droplet on an ideal wall with
different values of $\sigma_{01}$.

(a) $\theta=60^{\circ}$

(b) $\theta=120^{\circ}$

Figure 4: The equilibrium shapes of the droplet on an ideal wall with
different initial radii ($R_{0}=40,50,60$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$).

(a) $\theta=30^{\circ}$

(b) $\theta=45^{\circ}$

(c) $\theta=135^{\circ}$

(d) $\theta=150^{\circ}$

Figure 5: The equilibrium shapes of the droplet on an ideal wall at different contact angles. Table 2: The normalized droplet heights at different contact angles. Contact angle | $30^{\circ}$ | $45^{\circ}$ | $60^{\circ}$ | $75^{\circ}$ | $90^{\circ}$ | $105^{\circ}$ | $120^{\circ}$ | $135^{\circ}$ | $150^{\circ}$
---|---|---|---|---|---|---|---|---|---
Analytical | 0.5427 | 0.6822 | 0.8008 | 0.9013 | 1.0000 | 1.0985 | 1.1830 | 1.2628 | 1.3393
Numerical | 0.5579 | 0.6871 | 0.7996 | 0.9027 | 1.0000 | 1.0933 | 1.1825 | 1.2660 | 1.3389
Relative error | 2.72% | 0.72% | 0.15% | 0.16% | 0 | 0.48% | 0.04% | 0.25% | 0.03%

###  4.2 A droplet spreading on an ideal wall

To test the accuracy of the present _2f1s_ model, here we consider a droplet
spreading on an idea wall, and the solid plate with a finite width $H_{0}=10$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$ is absorbed into the computational
domain $[-L_{0}/2,L_{0}/2]\times[-H_{0},9H_{0}]$ with $L_{0}=300$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$. A semicircular droplet with the
radius $R_{0}$ is placed on the center of the wall, the periodic boundary
condition is applied in $x$-direction and no-flux boundary condition is
imposed on the top and bottom boundaries. Initially, the phase-field variables
are given by

$\displaystyle\phi_{0}\left(x,y\right)$
$\displaystyle=\frac{1}{2}-\frac{1}{2}\tanh\frac{2y}{D},$ (27)
$\displaystyle\phi_{1}\left(x,y\right)$
$\displaystyle=\left[1-\phi_{0}\left(x,y\right)\right]\left(\frac{1}{2}+\frac{1}{2}\tanh\frac{R_{0}-\sqrt{x^{2}+y^{2}}}{D/2}\right).$

In the simulations, some physical parameters are fixed as
$\rho_{0}:\rho_{1}:\rho_{2}=1$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{l}\mathrm{u}^{3}\mathrm{]}$,
$\mu_{0}=\mu_{1}=\mu_{2}=0.01$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{(}\mathrm{l}\mathrm{u}\,\mathrm{t}\mathrm{s}\mathrm{)}\mathrm{]}$
and $\sigma_{12}=0.001$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{t}\mathrm{s}^{2}\mathrm{]}$.
We first give a comparison of the results under different values of
$\sigma_{01}$ in Fig. 3. From this figure, one can observe that the spreading
of the droplet is not enough with a small $k=\sigma_{01}/\sigma_{12}$, while a
large value of $k$ leads to an inaccurate result, although the mechanical
balance condition (9) is preserved better. We believe that there should be an
optimal value of $k$ to give the best result, and the proposed equation (10)
seems to present a good performance, as shown in Fig. 3. We note that in this
case, the curvature of the solid plate is zero, and the parameter $k$ of Eq.
(10) is only related to the contact angle. To further test the accuracy of Eq.
(10), the droplet spreading on the solid plate with different initial radii
are considered. As seen from Fig. 4, the present numerical results agree well
with the analytical solutions. We also plot the results at a large range of
the contact angles in Fig. 5 where the initial radius of the droplet is 50
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$, and find that the finial shapes of
the droplet are exactly determined by the given contact angles. To give a
quantitative comparison, we also measure the height of the droplet from the
equilibrium shape in Table 2, and it is found that the numerical data are in
good agreement with the following analytical solution,

$H=R-R\cos\theta,\quad
R=R_{0}\sqrt{\frac{\pi/2}{\theta-\sin\theta\cos\theta}}.$ (28)

Here it should be noted that the droplet spreading length on the solid surface
is not measured because the diffuse-interface method for multiphase system
usually loses a certain accuracy at triple points.

(a) $\theta=45^{\circ}$

(b) $\theta=135^{\circ}$

Figure 6: The equilibrium shapes of the droplet on a cylinder surface with
different forms of mobility.

(a) $R_{0}=75$ $\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$, $R_{s}=50$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$

(b) $R_{0}=50$ $\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$, $R_{s}=50$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$

(c) $R_{0}=50$ $\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$, $R_{s}=75$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$

Figure 7: The equilibrium shapes of the droplet on a cylinder surface with
different radius ratios and $\theta=90^{\circ}$.

(a) $\theta=30^{\circ}$

(b) $\theta=60^{\circ}$

(c) $\theta=120^{\circ}$

(d) $\theta=150^{\circ}$

Figure 8: The equilibrium shapes of the droplet on a cylinder surface at
different contact angles. Figure 9: The normalized droplet heights at
different contact angles.

###  4.3 A droplet spreading on a cylinder surface

In this part, we consider a droplet spreading on a cylinder surface in the
square domain $[-L_{0}/2,L_{0}/2]\times[-L_{0}/4,3L_{0}/4]$ with $L_{0}=300$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$, and the other physical parameters
are the same as those in the previous section. In the computational domain, a
solid cylinder with the radius $R_{s}$ is located at the position
$\left(0,0\right)$, and a droplet with radius $R_{0}$ is placed on the
cylinder surface with the following initial phase-field variables,

$\displaystyle\phi_{0}\left(x,y\right)$
$\displaystyle=\frac{1}{2}+\frac{1}{2}\tanh\frac{R_{s}-\sqrt{x^{2}+y^{2}}}{D/2},$
(29) $\displaystyle\phi_{1}\left(x,y\right)$
$\displaystyle=\left[1-\phi_{0}\left(x,y\right)\right]\left[\frac{1}{2}+\frac{1}{2}\tanh\frac{R_{0}-\sqrt{x^{2}+\left(y-R_{s}\right)^{2}}}{D/2}\right].$

For this problem, the curvature of the cylinder is $1/R_{s}$, and thus we set
the dimensionless curvature to be $\kappa=R_{0}/R_{s}$. Usually, a constant
mobility in the CH equation for two-phase flows cannot preserve the
boundedness of the phase-field variable and the phase mass conservation.
However, a degenerate mobility can overcome these defects [56]. Here we
present a comparison between different forms of mobility through considering
the droplet spreading on a cylinder surface, and show the results in Fig. 6.
From this figure, one can find that the constant part of mobility causes the
obvious shrinkage of the droplet, and there is a mass leakage between the two
fluid phases, although the total mass is conserved. For this reason, the
degenerate mobility with $M_{cst}=0$ is used in the following simulations. On
the other hand, we also test Eq. (10) by considering the effect of
dimensionless curvature $\kappa$. As shown in Fig. 7, the numerical results
also match well with the analytical solutions. Additionally, we also present
the equilibrium shapes of the droplet at a large range of contact angles in
Fig. 8, and find that the droplets experience the greater deformations,
compared to those on the plate walls in Section 4.2. To conduct a quantitative
comparison, we further consider the maximum height of the equilibrium droplet
$H_{max}$, and it can be obtained by the initial shape under the condition of
equal area,

$H_{max}=r+k_{s},\quad\pi
r^{2}-r^{2}\cos^{-1}\left(\frac{k_{s}^{2}+r^{2}-R_{s}^{2}}{2k_{s}r}\right)-R_{s}^{2}\cos^{-1}\left(\frac{k_{s}^{2}-r^{2}+R_{s}^{2}}{2k_{s}R_{s}}\right)+k_{s}\sqrt{R_{s}^{2}-\left(\frac{k_{s}^{2}-r^{2}+R_{s}^{2}}{2k_{s}}\right)}=S_{0},$
(30)

where $S_{0}$ is the initial area, $k_{s}$ is related to the finial radius $r$
and the contact angle $\theta$,

$k_{s}=\sqrt{r^{2}+R_{s}^{2}-2rR_{s}\cos\theta}.$ (31)

Figure 9 plots the corresponding normalized maximum height of the droplet, and
it is found that the present results agree well with the analytical solutions
and previous data [57, 45].

Figure 10: The configuration of a droplet passing through $2\times 4$
cylindrical obstacles.

###  4.4 A droplet passing through $2\times 4$ cylindrical obstacles

We now study the problem of a droplet passing through $2\times 4$ cylindrical
obstacles with a large topological interface change. The cylindrical obstacles
with the radius $R_{s}=25$ $\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$ and
center distance $4R_{s}$ are placed at the central of the channel, the droplet
is initialized on the left side of the obstacles, and is driven by the left
inlet and right outlet velocity $U_{0}$ (see the schematic in Fig. 10). The
physical parameters are set to be $M_{deg}=0.01$,
$\rho_{0}:\rho_{1}:\rho_{2}=1$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{l}\mathrm{u}^{3}\mathrm{]}$,
$\mu_{0}=\mu_{1}=\mu_{2}=0.01$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{(}\mathrm{l}\mathrm{u}\,\mathrm{t}\mathrm{s}\mathrm{)}\mathrm{]}$,
$\sigma_{12}=0.001$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{t}\mathrm{s}^{2}\mathrm{]}$,
$\theta=170^{\circ}$, and Capillary number $Ca=\mu_{1}U_{0}/\sigma_{12}=0.01$.
We conduct some simulations with the channel length $L=2400$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$ and width $H=200$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$, and compare present results with
the available data [58] in Fig. 11. From this figure, one can see that present
results show the similar dynamic behaviors to those reported in the previous
work. Here it should be noted that the initial size of the droplet is not
clearly specified in Ref. [58], and there may be a difference in the droplet
size. From this test, one can conclude that the present model has a good
capability in the study of two-phase flows in complex geometries, and in the
following, it will be adopted to investigate two displacement problems.

Figure 11: The snapshots of the droplet passing through $2\times 4$
cylindrical obstacles [Results reported in Ref. [58]: (a), (c), (e), (g),
present results: (b), (d), (f), (h)].

###  4.5 Displacement in a Tesla valve

Tesla valve is a kind of check valve with no-moving-part, invented by Nicolas
Tesla in 1920 [59]. This special pipe structure usually presents different
flow patterns when the fluid flows in forward and reverse directions, and can
be considered as a fluidic diode. Here we consider the displacement of two-
phase flow in a 4 stages 3D T45C Tesla valve, and present the 2D projected
structure of 1 stage T45C Tesla valve in Fig. 12 where the pipe width is
$d_{s}=15$ $\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$ with 1
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$$=4\times 10^{-4}$
$\mathrm{[}\mathrm{m}\mathrm{]}$, and other structure parameters can be found
in Ref. [60]. The material properties of the driving fluid (phase 1) and
displaced fluid (phase 2) are set to be the same for simplicity, i.e.,
$\rho_{1}=\rho_{2}=1$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{l}\mathrm{u}^{3}\mathrm{]}$,
$\mu_{1}=\mu_{2}=0.1$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{(}\mathrm{l}\mathrm{u}\,\mathrm{t}\mathrm{s}\mathrm{)}\mathrm{]}$,
and other parameters are fixed as $D=4$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$, $\sigma_{12}=0.004$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{t}\mathrm{s}^{2}\mathrm{]}$
and $\theta=60^{\circ}$. In our simulations, the computational domain is
$\left(462+2\Delta\right)\times 20\times 136$
$\mathrm{[}\mathrm{l}\mathrm{u}^{3}\mathrm{]}$ with $\Delta=30$ being the
number of fluid layers added to the two sides of the domain in $x$-direction,
as shown in Fig. 13.

Initially, the Tesla valve is full of the displaced fluid except for the
driving fluid in the inlet section with the distance $d_{s}$ (see Fig. 13),
and the fluids are driven by a pressure drop $\Delta P$ from inlet to outlet.
We set $\Delta P=0.1$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{(}\mathrm{l}\mathrm{u}\,\mathrm{t}\mathrm{s}^{2}\mathrm{)}\mathrm{]}$
for both forward and reverse two-phase displacements, and display the
numerical results in Figs. 14 and 15. From these figures, one can observe that
the driving fluid gradually fills the pipe under the action of the pressure
drop in both forward and reverse cases. However, there are some differences in
the middle process. In the forward case, when the driving fluid reaches to the
bifurcation of the pipe, it tends to pass the straight pipe and the fluid in
straight pipe goes faster than the another part in bend pipe. In contrast, in
the reverse displacement, the fluids seem to have the same velocities in the
two bifurcated pipes in each stage valve, and the driving fluid has the same
tip positions in the two bifurcated pipes in $x$-direction. Additionally, one
can also find that the trapped displaced fluid at the early process of the
displacement disappears at the later process, which seems contrary to common
sense. To explain this phenomenon, we plot the vertical profile of the pipe at
$x=300$ $\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$ and $t=1.5\times 10^{5}$
$\mathrm{[}\mathrm{t}\mathrm{s}\mathrm{]}$ in Fig. 16. From this figure, one
can observe that the driving fluid is not filled with the square profile of
the pipe due to the larger inertia force caused by the pressure drop. As a
result, the displaced fluid cannot be completely trapped, and the trapped
parts can be displaced from four corners as time goes on.

Figure 12: The 2D projected structure of the one stage T45C Tesla valve.

(a) Forward flow

(b) Reverse flow

Figure 13: The 3D structure of the 4 stages T45C Tesla valve and the
initialization of driving phase.

(a) $t=5\times 10^{4}$ $\mathrm{[}\mathrm{t}\mathrm{s}\mathrm{]}$

(b) $t=1\times 10^{5}$ $\mathrm{[}\mathrm{t}\mathrm{s}\mathrm{]}$

(c) $t=1.5\times 10^{5}$ $\mathrm{[}\mathrm{t}\mathrm{s}\mathrm{]}$

(d) $t=2\times 10^{5}$ $\mathrm{[}\mathrm{t}\mathrm{s}\mathrm{]}$

Figure 14: The snapshots of the forward displacement in the 4 stages T45C
Tesla valve.

(a) $t=5\times 10^{4}$ $\mathrm{[}\mathrm{t}\mathrm{s}\mathrm{]}$

(b) $t=1\times 10^{5}$ $\mathrm{[}\mathrm{t}\mathrm{s}\mathrm{]}$

(c) $t=1.5\times 10^{5}$ $\mathrm{[}\mathrm{t}\mathrm{s}\mathrm{]}$

(d) $t=2\times 10^{5}$ $\mathrm{[}\mathrm{t}\mathrm{s}\mathrm{]}$

Figure 15: The snapshots of the reverse displacement in the 4 stages T45C
Tesla valve.

(a) Forward flow

(b) Reverse flow

Figure 16: The profiles of the displacements at $x=300$
$\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$ and $t=1.5\times 10^{5}$
$\mathrm{[}\mathrm{t}\mathrm{s}\mathrm{]}$. Here the blue solid lines are the
boundaries of pipes and red dotted lines represent the interfaces of two fluid
phases. Figure 17: The normalized distance of the tip of the driving fluid
from the inlet.

The influence of the viscosity ratio ($M_{\mu}=\mu_{1}/\mu_{2}$) is further
investigated, and $\mu_{1}$ is fixed at 0.1. We focus on two cases with
$M_{\mu}=0.2,1$ and 5, and show the results of both forward and reverse
displacements in Fig. 17. As shown in this figure, the normalized distances of
the tip of the driving fluid from inlet in forward and reverse cases are
basically the same as each other under the present pressure drop. In addition,
with the decrease of $M_{\mu}$, the displaced fluid is more difficult to be
displaced due to a larger viscous resistance, and thus the value of the
distance $d_{tip}^{*}$ takes more time to reach 1.

###  4.6 Displacement in a porous medium

We finally consider the displacement in a 3D porous medium, and the properties
of fluids are the same as those in Section 4.5. Figure 18(a) shows a sample of
the poorly sorted unconsolidated fluvial sandpack with a porosity
$\varepsilon_{o}=0.3585$ [61], and the original $256\times 256\times 256$ data
is firstly compressed into the size of $128\times 128\times 128$
$\mathrm{[}\mathrm{l}\mathrm{u}^{3}\mathrm{]}$ to adapt to the memory
limitation of the computing platform, as shown in Fig. 18(b). To apply the
present diffuse-interface method, the compressed data needs to be smoothed by
the finite-time evolution of a standard CH equation for two phases, and Fig.
18(c) plots the smoothed data after 200 iterations. In addition, the periodic
boundary condition with a pressure drop [62] is imposed on the inlet and
outlet boundaries in $x$-direction and no-flux boundary condition is applied
on other boundaries, which preserves the mass conservation of the system.

Figure 18: A sample of the poorly sorted unconsolidated fluvial sandpack [(a)
The original $256\times 256\times 256$ micro-CT image with porosity
$\varepsilon_{o}=0.3585$, (b) the compressed $128\times 128\times 128$ data
with porosity $\varepsilon_{c}=0.3024$, (c) the smoothed $128\times 128\times
128$ data with porosity $\varepsilon_{s}=0.3068$].

The fluids are driven by a pressure drop $\Delta p=5\times 10^{-3}$
$\mathrm{[}\mathrm{m}\mathrm{u}\mathrm{/}\mathrm{(}\mathrm{l}\mathrm{u}\,\mathrm{t}\mathrm{s}^{2}\mathrm{)}\mathrm{]}$,
and some simulations are carried out to predict the relative permeabilities of
the two-phase fluids in the porous medium. When the flows reach the steady
state, the relative permeability of a specific fluid can be calculated by

$K_{p}=\frac{\int_{\bar{\Omega}}u_{x,p}\mathrm{d}\bar{\Omega}}{\int_{\bar{\Omega}}u^{\prime}_{x,p}\mathrm{d}\bar{\Omega}},\quad
p=1,2,$ (32)

where $\bar{\Omega}$ is the pore spaces ($\phi_{0}\leq 0.5$) of the porous
medium, $u_{x,p}$ and $u^{\prime}_{x,p}$ are the $x$-component velocities of
fluid $p$ in two-phase flow and single-phase flow, respectively.

Figure 19: A example of the initial distribution of wetting fluid with
$R_{0}=48$ $\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$.

(a) $S_{1}=0.1668$

(b) $S_{1}=0.2160$

(c) $S_{1}=0.2690$

(d) $S_{1}=0.3908$

(e) $S_{1}=0.5367$

(f) $S_{1}=0.7128$

Figure 20: The finial distributions of wetting fluid in the porous medium at
different saturations. Figure 21: The relative permeabilities of wetting fluid
($K_{1}$) and non-wetting fluid ($K_{2}$) in the porous medium at
$\theta=60^{\circ}$. Figure 22: The relative permeabilities of wetting fluid
($K_{1}$) and non-wetting fluid ($K_{2}$) in the porous medium at
$S_{1}=0.4605$.

(a) $M_{\mu}=0.02$

(b) $M_{\mu}=1$

(c) $M_{\mu}=50$

Figure 23: The finial distributions of wetting fluid in the porous medium with
different viscosity ratios at $S_{1}=0.4605$. Figure 24: The relative
permeabilities of wetting fluid ($K_{1}$) and non-wetting fluid ($K_{2}$) in
the porous medium at $S_{1}=0.4605$.

(a) $\theta=30^{\circ}$

(b) $\theta=60^{\circ}$

(c) $\theta=90^{\circ}$

Figure 25: The finial distributions of wetting fluid in the porous medium with
different contact angles at $S_{1}=0.4605$.

We first investigate the effect of the saturation on the relative
permeabilities under the condition of $\theta=60^{\circ}$. The saturation of
wetting fluid $S_{1}$ is determined by the radius $R_{0}$ and the following
initialization of $\phi_{1}$,

$\phi_{1}\left(x,y,z\right)=\left[1-\phi_{0}\left(x,y,z\right)\right]\left[\frac{1}{2}+\frac{1}{2}\tanh\frac{R_{0}-\sqrt{\left(y-Ly/2\right)^{2}+\left(z-Lz/2\right)^{2}}}{D/2}\right],$
(33)

and Fig. 19 shows an example of the initial distribution of wetting fluid with
$R_{0}=48$ $\mathrm{[}\mathrm{l}\mathrm{u}\mathrm{]}$. After some simulations,
we plot the finial distributions of wetting fluid under different saturations
(or radii) in Fig. 20. From this figure, one can observe that when the
saturation of wetting fluid is increased, the wetting fluid flows into more
pores and the connectivity of wetting phase gradually increases. To
quantitatively analyze the present results, we also calculate the relative
permeabilities at different saturations, and plot them in Fig. 21. As shown in
this figure, the relative permeability of the wetting fluid increases with the
increase of $S_{1}$, while that of the non-wetting fluid decreases. Actually,
when the value of $S_{1}$ is small, as shown in Fig. 20(a), the wetting fluid
has fewer pathways in $x$-direction, which causes a low Darcy velocity. As a
result, the relative permeability of wetting fluid is small. In contrast, when
$S_{1}$ is high, the wetting fluid fills into the most of the pores, and forms
many flow paths from inlet to outlet such that the relative permeability of
wetting fluid becomes large. In addition, when the wetting fluid saturation is
$50\%$, the relative permeability of the wetting fluid is larger than that of
non-wetting fluid because of the effect of the wettability and the internal
distribution of wetting fluid in the porous medium.

We then consider the influence of the viscosity ratio
$M_{\mu}=\mu_{1}/\mu_{2}$ with $\mu_{1}=0.1$, and the saturation of the
wetting fluid is set as $S_{1}=0.4605$ ($R_{0}=48$). As seen from Fig. 22,
when the viscosity ratio is varied from $M_{\mu}=0.02$ to $M_{\mu}=50$, the
flow resistance exerted on the non-wetting fluid decreases, which results in
the increase of the relative permeability $K_{2}$. However, the variation of
the relative permeability $K_{1}$ shows some differences. When $M_{\mu}<1$,
the viscosity of non-wetting fluid decreases with the increase of $M_{\mu}$,
and thus the wetting fluid is easier to drive the non-wetting fluid and has a
better connectivity [see the final shapes in Figs. 23(a) and (b)]. However,
when $\mu_{2}<\mu_{1}$, the increase of $K_{1}$ is small because of the
dominant flow formed by the non-wetting fluid [see Figs. 23(b) and (c)].

The effect of the wettability is also discussed at $S_{1}=0.4605$. Figure 24
plots the relative permeability $K_{p}$ as a function of the contact angle.
From this figure, one can find that the contact angle has a significant effect
on the relative permeabilities. With the increase of the contact angle, the
relative permeability of wetting fluid increases a little, while that of non-
wetting fluid decreases a lot, this is because of the weakening of the
adhesive force between the wetting fluid and solid surface and the
strengthening of the adhesive force between the non-wetting fluid and solid
surface. We also show the distributions of wetting fluid at different contact
angles in Fig. 25, and find that when the solid surface become more
hydrophilic, the wetting fluid prefers to flow into small pores to minimize
the surface energy.

##  5 Conclusions

In this paper, we first proposed a ternary phase-field model for two-phase
flows in complex geometries. In this model, a specific expression of the
liquid-solid surface tension coefficient is designed, and it is assumed to be
related to the curvature and the contact angle of the solid surface. Then the
reformulated consistent NS equations is used to describe the fluid flows, in
which the phase-field variable for solid phase is introduced to achieve a high
viscosity in solid phase and preserve the velocity boundary conditions on
solid surface. To solve the present CH and NS equations for two-phase flows in
complex geometries, the Hermite-moment based LB method is further developed,
in which the collision step is conducted in the natural moment space with a
lower triangular relaxation matrix, and an adjustable scale factor is also
introduced to improve the numerical stability. To test the present MRT-LB
method, several benchmark problems are considered, and the numerical results
are in good agreement with the analytical solutions. Finally, the present
model is applied to study two 3D displacement problems in a complex channel
and a porous medium, and the simulation results illustrate the capacity of the
present model in the study of two-phase flows in complex geometries. In the
future, the present model will be extended to investigate the problems of two-
phase flows with moving solid objects or solid-liquid phase transition.

## Acknowledgments

This research has been supported by the National Natural Science Foundation of
China under Grants No. 12072127 and No 51836003, and the Interdiciplinary
Research Program of HUST (2023JCYJ002). The computation was completed on the
HPC Platform of Huazhong University of Science and Technology.

##  Appendix A The weight coefficients and transformation matrices of D2Q9 and
Q3Q15 lattice structures

In this appendix, we list the related weight coefficients and transformation
matrices of the D2Q9 and D3Q15 lattice structures. Here we set
$d_{0}^{p}=d_{0}$ and $\omega_{i,p}=\omega_{i}$ for simplicity, and
$d_{0}=0.5$ is used in our numerical simulations.

D2Q9: $\omega_{0}=1-2d_{0}+d_{0}^{2}$,
$\omega_{1-4}=\left(d_{0}-d_{0}^{2}\right)/2$, and $\omega_{5-8}=d_{0}^{2}/4$.
$\mathbf{C}_{d}=\mathbf{diag}\left(1,c,c,c^{2},c^{2},c^{2},c^{3},c^{3},c^{4}\right)$,
$\mathbf{K}=\mathbf{diag}\left(s_{0},s_{1},s_{1},s_{2a},s_{2a},s_{2b},s_{3},s_{3},s_{4}\right)$,

$\mathbf{M}_{0}=\begin{pmatrix}[r]1&1&1&1&1&1&1&1&1\\\ 0&1&0&-1&0&1&-1&-1&1\\\
0&0&1&0&-1&1&1&-1&-1\\\ 0&1&0&1&0&1&1&1&1\\\ 0&0&1&0&1&1&1&1&1\\\
0&0&0&0&0&1&-1&1&-1\\\ 0&0&0&0&0&1&-1&-1&1\\\ 0&0&0&0&0&1&1&-1&-1\\\
0&0&0&0&0&1&1&1&1\\\
\end{pmatrix},\quad\mathbf{N}=\begin{pmatrix}[r]1&0&0&0&0&0&0&0&0\\\
0&1&0&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0&0\\\ -c_{s}^{2}&0&0&1&0&0&0&0&0\\\
-c_{s}^{2}&0&0&0&1&0&0&0&0\\\ 0&0&0&0&0&1&0&0&0\\\
0&-c_{s}^{2}&0&0&0&0&1&0&0\\\ 0&0&-c_{s}^{2}&0&0&0&0&1&0\\\
c_{s}^{4}&0&0&-c_{s}^{2}&-c_{s}^{2}&0&0&0&1\\\ \end{pmatrix}.$ (34)

D3Q15: $\omega_{0}=1-3d_{0}+2d_{0}^{2}$,
$\omega_{1-6}=\left(d_{0}-d_{0}^{2}\right)/2$, and
$\omega_{7-14}=d_{0}^{2}/8$.
$\mathbf{C}_{d}=\mathbf{diag}\left(1,c,c,c,c^{2},c^{2},c^{2},c^{2},c^{2},c^{2},c^{3},c^{3},c^{3},c^{3},c^{4}\right)$,
$\mathbf{K}=\mathbf{diag}\left(s_{0},s_{1},s_{1},s_{1},s_{2a},s_{2a},s_{2a},s_{2b},s_{2b},s_{2b},s_{3},s_{3},s_{3},s_{3},s_{4}\right)$,

$\mathbf{M}_{0}=\setcounter{MaxMatrixCols}{15}\begin{pmatrix}[r]1&1&1&1&1&1&1&1&1&1&1&1&1&1&1\\\
0&1&-1&0&0&0&0&1&-1&1&-1&1&-1&1&-1\\\ 0&0&0&1&-1&0&0&1&1&-1&-1&1&1&-1&-1\\\
0&0&0&0&0&1&-1&1&1&1&1&-1&-1&-1&-1\\\ 0&1&1&0&0&0&0&1&1&1&1&1&1&1&1\\\
0&0&0&1&1&0&0&1&1&1&1&1&1&1&1\\\ 0&0&0&0&0&1&1&1&1&1&1&1&1&1&1\\\
0&0&0&0&0&0&0&1&-1&-1&1&1&-1&-1&1\\\ 0&0&0&0&0&0&0&1&1&-1&-1&-1&-1&1&1\\\
0&0&0&0&0&0&0&1&-1&1&-1&-1&1&-1&1\\\ 0&0&0&0&0&0&0&1&-1&1&-1&1&-1&1&-1\\\
0&0&0&0&0&0&0&1&1&-1&-1&1&1&-1&-1\\\ 0&0&0&0&0&0&0&1&1&1&1&-1&-1&-1&-1\\\
0&0&0&0&0&0&0&1&-1&-1&1&-1&1&1&-1\\\ 0&0&0&0&0&0&0&1&1&1&1&1&1&1&1\\\
\end{pmatrix},$ (35)
$\mathbf{N}_{H}=\setcounter{MaxMatrixCols}{15}\begin{pmatrix}[r]1&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&1&0&0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&1&0&0&0&0&0&0&0&0&0&0&0\\\ -c_{s}^{2}&0&0&0&1&0&0&0&0&0&0&0&0&0&0\\\
-c_{s}^{2}&0&0&0&0&1&0&0&0&0&0&0&0&0&0\\\
-c_{s}^{2}&0&0&0&0&0&1&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&1&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&1&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&1&0&0&0&0&0\\\
0&-c_{s}^{2}&0&0&0&0&0&0&0&0&1&0&0&0&0\\\
0&0&-c_{s}^{2}&0&0&0&0&0&0&0&0&1&0&0&0\\\
0&0&0&-c_{s}^{2}&0&0&0&0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&1&0\\\
2c_{s}^{4}&0&0&0&-c_{s}^{2}&-c_{s}^{2}&-c_{s}^{2}&0&0&0&0&0&0&0&1\\\
\end{pmatrix}.$ (36)

Here it should be noted that in the D3Q15 lattice structure, the Hermite
matrix $\mathbf{H}=\mathbf{N}_{H}\mathbf{M}$ cannot be directly obtained from
the commonly used Hermite polynomials, and some row transformations for the
Hermite moments are needed to satisfy the weighted orthogonality of
$\mathbf{H}$.

## References

  * [1] V. Alvarado, E. Manrique, Enhanced oil recovery: An update review, Energies 3 (9) (2010) 1529–1575.
  * [2] A. Maghzi, S. Mohammadi, M. H. Ghazanfari, R. Kharrat, M. Masihi, Monitoring wettability alteration by silica nanoparticles during water flooding to heavy oils in five-spot systems: A pore-level investigation, Experimental Thermal and Fluid Science 40 (2012) 168–176.
  * [3] C. Fan, D. Elsworth, S. Li, L. Zhou, Z. Yang, Y. Song, Thermo-hydro-mechanical-chemical couplings controlling CH4 production and CO2 sequestration in enhanced coalbed methane recovery, Energy 173 (2019) 1054–1077.
  * [4] S. Korneev, Z. Wang, V. Thiagarajan, S. Nelaturi, Fabricated shape estimation for additive manufacturing processes with uncertainty, Computer-Aided Design 127 (2020) 102852.
  * [5] B. Blakey-Milner, P. Gradl, G. Snedden, M. Brooks, J. Pitot, E. Lopez, M. Leary, F. Berto, A. du Plessis, Metal additive manufacturing in aerospace: A review, Materials & Design 209 (2021) 110008.
  * [6] L. M. Abriola, Modeling multiphase migration of organic chemicals in groundwater systems–a review and assessment, Environmental Health Perspectives 83 (1989) 117–143.
  * [7] J. W. Cahn, C. M. Elliott, A. Novick-Cohen, The Cahn–Hilliard equation with a concentration dependent mobility: Motion by minus the laplacian of the mean curvature, European Journal of Applied Mathematics 7 (3) (1996) 287–301.
  * [8] V. Badalassi, H. Ceniceros, S. Banerjee, Computation of multiphase systems with phase field models, Journal of Computational Physics 190 (2) (2003) 371–397.
  * [9] D. Jacqmin, Calculation of two-phase Navier–Stokes flows using phase-field modeling, Journal of Computational Physics 155 (1) (1999) 96–127.
  * [10] H. Abels, D. Depner, H. Garcke, On an incompressible Navier–Stokes/Cahn–Hilliard system with degenerate mobility, Annales de l’Institut Henri Poincaré C, Analyse non linéaire 30 (6) (2013) 1175–1190.
  * [11] C. S. Peskin, Numerical analysis of blood flow in the heart, Journal of Computational Physics 25 (3) (1977) 220–252.
  * [12] A. L. Fogelson, C. S. Peskin, A fast numerical method for solving the three-dimensional Stokes’ equations in the presence of suspended particles, Journal of Computational Physics 79 (1) (1988) 50–69.
  * [13] P. Angot, Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows, Mathematical Methods in the Applied Sciences 22 (16) (1999) 1395–1412.
  * [14] I. Ramière, P. Angot, M. Belliard, A fictitious domain approach with spread interface for elliptic problems with general boundary conditions, Computer Methods in Applied Mechanics and Engineering 196 (4) (2007) 766–781.
  * [15] P. Angot, C.-H. Bruneau, P. Fabrie, A penalization method to take into account obstacles in incompressible viscous flows, Numerische Mathematik 81 (1999) 497–520.
  * [16] D. Kolomenskiy, K. Schneider, A Fourier spectral method for the Navier–Stokes equations with volume penalization for moving solid obstacles, Journal of Computational Physics 228 (16) (2009) 5687–5709.
  * [17] Y. Nakayama, R. Yamamoto, Simulation method to resolve hydrodynamic interactions in colloidal dispersions, Physical Review E 71 (2005) 036707.
  * [18] Y. Nakayama, K. Kim, R. Yamamoto, Simulating (electro)hydrodynamic effects in colloidal dispersions: Smoothed profile method, The European Physical Journal E 26 (2008) 361–368.
  * [19] X. Li, J. Lowengrub, A. Rätz, A. Voigt, Solving pdes in complex geometries: A diffuse domain approach, Communications in Mathematical Sciences 7 (1) (2009) 81–107.
  * [20] K. Y. Lervag, J. Lowengrub, Analysis of the diffuse-domain method for solving PDEs in complex geometries, Communications in Mathematical Sciences 13 (6) (2015) 1473–1500.
  * [21] P. Horgue, M. Prat, M. Quintard, A penalization technique applied to the “Volume-Of-Fluid” method: Wettability condition on immersed boundaries, Computers & Fluids 100 (2014) 255–266.
  * [22] J. K. Patel, G. Natarajan, Diffuse interface immersed boundary method for multi-fluid flows with arbitrarily moving rigid bodies, Journal of Computational Physics 360 (2018) 202–228.
  * [23] A. Shahmardi, M. E. Rosti, O. Tammisola, L. Brandt, A fully Eulerian hybrid immersed boundary-phase field model for contact line dynamics on complex geometries, Journal of Computational Physics 443 (2021) 110468.
  * [24] E. L. Sharaborin, O. A. Rogozin, A. R. Kasimov, The coupled volume of fluid and Brinkman penalization methods for simulation of incompressible multiphase flows, Fluids 6 (9) (2021).
  * [25] S. Aland, J. Lowengrub, A. Voigt, Two-phase flow in complex geometries: A diffuse domain approach, Computer Modeling in Engineering & Sciences 57 (1) (2010) 77–106.
  * [26] Z. Guo, F. Yu, P. Lin, S. Wise, J. Lowengrub, A diffuse domain method for two-phase flows with large density ratio in complex geometries, Journal of Fluid Mechanics 907 (2021) A38.
  * [27] J. Yang, Z. Tan, J. Wang, J. Kim, Modified diffuse interface fluid model and its consistent energy-stable computation in arbitrary domains, Journal of Computational Physics 488 (2023) 112216.
  * [28] F. Boyer, C. Lapuerta, Study of a three component Cahn-Hilliard flow model, ESAIM: Mathematical Modelling and Numerical Analysis 40 (4) (2006) 653–687.
  * [29] J. Kim, Phase field computations for ternary fluid flows, Computer Methods in Applied Mechanics and Engineering 196 (45) (2007) 4779–4788.
  * [30] F. Boyer, S. Minjeaud, Hierarchy of consistent $n$-component Cahn–Hilliard systems, Mathematical Models and Methods in Applied Sciences 24 (14) (2014) 2885–2928.
  * [31] S. Dong, Multiphase flows of $n$ immiscible incompressible fluids: A reduction-consistent and thermodynamically-consistent formulation and associated algorithm, Journal of Computational Physics 361 (2018) 1–49.
  * [32] S. Yi, A phase field method for the numerical simulation of rigid particulate in two-phase flows, Fluid Dynamics Research 52 (1) (2020) 015512.
  * [33] C. Rohde, L. von Wolff, A ternary Cahn–Hilliard–Navier–Stokes model for two-phase flow with precipitation and dissolution, Mathematical Models and Methods in Applied Sciences 31 (01) (2021) 1–35.
  * [34] J. R. Panter, A. R. Konicek, M. A. King, A. Jusufi, M. S. Yeganeh, H. Kusumaatmaja, Rough capillary rise, Communications Physics 6 (1) (2023) 44\.
  * [35] P. Yue, C. Zhou, J. J. Feng, Spontaneous shrinkage of drops and mass conservation in phase-field simulations, Journal of Computational Physics 223 (1) (2007) 1–9.
  * [36] F. J. Higuera, S. Succi, R. Benzi, Lattice gas dynamics with enhanced collisions, Europhysics Letters 9 (4) (1989) 345.
  * [37] S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid flows, Annual Review of Fluid Mechanics 30 (1) (1998) 329–364.
  * [38] C. K. Aidun, J. R. Clausen, Lattice-Boltzmann method for complex flows, Annual Review of Fluid Mechanics 42 (1) (2010) 439–472.
  * [39] T. Young, III. An essay on the cohesion of fluids, Philosophical Transactions of the Royal Society of London 95 (1805) 65–87.
  * [40] A. Karma, Phase-field formulation for quantitative modeling of alloy solidification, Physical Review Letters 87 (2001) 115701.
  * [41] B. Echebarria, R. Folch, A. Karma, M. Plapp, Quantitative phase-field model of alloy solidification, Physical Review E 70 (2004) 061604.
  * [42] Z. Chai, B. Shi, Multiple-relaxation-time lattice Boltzmann method for the Navier-Stokes and nonlinear convection-diffusion equations: Modeling, analysis, and elements, Physical Review E 102 (2020) 023306.
  * [43] Z. Chai, X. Yuan, B. Shi, Rectangular multiple-relaxation-time lattice Boltzmann method for the Navier-Stokes and nonlinear convection-diffusion equations: General equilibrium and some important issues, Physical Review E 108 (2023) 015304.
  * [44] C. Zhan, Z. Chai, B. Shi, Consistent and conservative phase-field-based lattice Boltzmann method for incompressible two-phase flows, Physical Review E 106 (2022) 025319.
  * [45] X. Liu, Z. Chai, C. Zhan, B. Shi, W. Zhang, A diffuse-domain phase-field lattice Boltzmann method for two-phase flows in complex geometries, Multiscale Modeling & Simulation 20 (4) (2022) 1411–1436.
  * [46] D. d’Humières, Generalized Lattice-Boltzmann Equations, Vol. 159, AIAA, Washington, DC, 1992, pp. 450–458.
  * [47] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability, Physical Review E 61 (2000) 6546–6562.
  * [48] P. V. Coveney, S. Succi, D. d’Humières, I. Ginzburg, M. Krafczyk, P. Lallemand, L.-S. Luo, Multiple–relaxation–time lattice Boltzmann models in three dimensions, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 360 (1792) (2002) 437–451.
  * [49] M. Geier, A. Greiner, J. G. Korvink, Cascaded digital lattice Boltzmann automata for high reynolds number flow, Physical Review E 73 (2006) 066705.
  * [50] K. N. Premnath, S. Banerjee, Inertial frame independent forcing for discrete velocity Boltzmann equation: Implications for filtered turbulence simulation, Communications in Computational Physics 12 (3) (2012) 732–766.
  * [51] C. Coreixas, G. Wissocq, G. Puigt, J.-F. Boussuge, P. Sagaut, Recursive regularization step for high-order lattice Boltzmann methods, Physical Review E 96 (2017) 033306.
  * [52] C. Coreixas, B. Chopard, J. Latt, Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations, Physical Review E 100 (2019) 033305.
  * [53] K. K. Mattila, P. C. Philippi, J. Hegele, Luiz A., High-order regularization in lattice-Boltzmann equations, Physics of Fluids 29 (4) (2017) 046103.
  * [54] A. J. C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation, Journal of Fluid Mechanics 271 (1994) 285–309.
  * [55] Z. Guo, C. Zheng, B. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chinese Physics 11 (4) (2002) 366.
  * [56] C. M. Elliott, H. Garcke, On the Cahn–Hilliard equation with degenerate mobility, SIAM Journal on Mathematical Analysis 27 (2) (1996) 404–423.
  * [57] A. Fakhari, D. Bolster, Diffuse interface modeling of three-phase contact line dynamics on curved boundaries: A lattice Boltzmann model for large density and viscosity ratios, Journal of Computational Physics 334 (2017) 620–638.
  * [58] C. Chung, M. Lee, K. Char, K. H. Ahn, S. J. Lee, Droplet dynamics passing through obstructions in confined microchannel flow, Microfluidics and Nanofluidics 9 (2010) 1151–1163.
  * [59] N. Tesla, Valvular conduit, U.S. Patent No. 1329559 (1920).
  * [60] Y. Bao, H. Wang, Numerical study on flow and heat transfer characteristics of a novel Tesla valve with improved evaluation method, International Journal of Heat and Mass Transfer 187 (2022) 122540.
  * [61] J. E. Santos, M. J. Pyrcz, M. Prodanović, 3D Dataset of binary images: A collection of synthetically created digital rock images of complex media, Data in Brief 40 (2022) 107797.
  * [62] S. H. Kim, H. Pitsch, A generalized periodic boundary condition for lattice Boltzmann method simulation of a pressure driven flow in a periodic geometry, Physics of Fluids 19 (10) (2007) 108101.