A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows Xiaolei Yuan a,b , Zhenhua Chai a,b , Huili Wang c , Baochang Shi a,b ,

∗

a

School of Mathematics and Statistics, 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 c School of Mathematics and Computer Science, Wuhan Textile University, Wuhan, 430073, China b

article

info

Article history: Received 1 April 2019 Received in revised form 19 July 2019 Accepted 5 October 2019 Available online xxxx Keywords: Generalized lattice Boltzmann model Source term Incompressible and nearly incompressible N–S equations Fluid flow system Two-phase flow

a b s t r a c t In this paper, a generalized lattice Boltzmann (LB) model with a source term in the continuity equation is proposed to solve both incompressible and nearly incompressible Navier–Stokes (N–S) equations. This model can be used to deal with single-phase and two-phase flows problems with a source term in the continuity equation. From this generalized model, we can not only get some existing models, but also derive new models. Moreover, for the incompressible model derived, a modified pressure scheme is introduced to calculate the pressure, and then to ensure the accuracy of the model. In this work, we will focus on a two-phase flow system, and in the frame work of our generalized LB model, a new phase-field-based LB model is developed for incompressible and quasi-incompressible two-phase flows. A series of numerical simulations of some classic physical problems, including a spinodal decomposition, a static droplet, a layered Poiseuille flow, and a bubble rising flow under buoyancy, are performed to validate the developed model. Besides, some comparisons with previous quasi-incompressible and incompressible LB models are also carried out, and the results show that the present model is accurate in the study of two-phase flows. Finally, we also conduct a comparison between quasi-incompressible and incompressible LB models for two-phase flow problems, and find that in some cases, the proposed quasi-incompressible LB model performs better than incompressible LB models. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Fluid flow problems are ubiquitous in nature and engineering applications, such as single and multiphase flows [1–8], thermal flows [9,10], porous media flows [5,11,12], turbulent flows [13,14], external aerodynamics [15] and free-surface flow problems [16]. Usually, people can use Navier–Stokes equations to describe these problems. For some fluid flow problems, such as multiphase flows and axisymmetric fluid system, the right-hand side term of the continuity equation may not be zero [17–21]. In other words, the continuity equation may contain a source term. It is very meaningful to study the problems of containing a source term in the continuity equation. Due to limitations of theoretical research and experimental methods, it is especially necessary to develop a numerical method to solve the N–S equations with a source term in the continuity equation. ∗ Corresponding author at: School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, China. E-mail address: [email protected] (B. Shi). https://doi.org/10.1016/j.camwa.2019.10.007 0898-1221/© 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

2

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

As an efficient numerical method, the lattice Boltzmann method has made rapid progress since its appearance in the late 1980s due to its simplicity, scalability on parallel computers, and ease to handle complex geometries [22–24]. This method has gained great success in the study of fluid flow system [19,25–29]. In general, LB methods for dealing with such flow problems can be divided into two categories. One is the nearly incompressible model and the other is the so-called incompressible model. For the nearly incompressible model, the macroscopic quantities those need to be solved are the density ρ and fluid velocity u, and the pressure can be obtained from the equation of state (p = ρ cs2 ). While for the incompressible model, the macroscopic quantities we need to calculate are fluid velocity u and the pressure p, where density ρ is viewed as a constant. In the past, people always considered these two types of models separately. Actually, from the perspective of model construction, one can design a generalized LB model to deal with both incompressible and nearly incompressible problems, which is the main motivation of this paper. Moreover, some previous models have more or less assumptions in the derivation process, and often cannot recover the macroscopic equations completely. In addition, under the framework of LB methods, people have less research on N–S equations with a source term in the continuity equation. For instance, Halliday et al. [30] proposed a single-relaxation-time (SRT) LB model including a source term, while they employed a non-local scheme to calculate the spatial derivatives which appear in the source term. Cheng et al. [31] presented another LB model with a general source term and adopted a non-local scheme for temporal and spatial derivatives. Aursjø et al. [21] also developed an SRT model with a source term, which does not include temporal and spatial derivatives in the source term, and it preserves the Galilean invariance. We note that the above models are limited to nearly incompressible situations. Up to now, there is no available work on incompressible N–S equations with a source term in the continuity equation. Considering the above points, in this work, we will develop a generalized SRT LB model for both incompressible and nearly incompressible N–S equations with a source term in the continuity equation, and this model is also an extension of existing models. From the generalized model, we can not only get some existing models, but also derive new models. Among these new models, we can obtain an incompressible model for N–S equations with a source term, and we present a modified scheme to calculate the pressure p, which is more accurate than the previous one. Simultaneously, our generalized model can recover the macroscopic equations without any unnecessary assumptions. Finally, we would like to point that this model can be used not only to solve single-phase fluid flows, but also to design models for two-phase flows. Based on our generalized model, we also present a phase-field-based LB model for two-phase flows. In phase-field theory, there are two types of classical phase-field model, namely Model H [32] (referred to as the incompressible model in this paper) and the quasi-incompressible model [33]. In these two classic models, the source term may appear at the right end of the continuity equation. Therefore, we use the N–S equations with a source term to solve the two models uniformly. In fact, our main purpose of this paper is to propose a generalized phase-field LB model that contains these two types of phase-field models. The proposed model can be used to handle single-phase or two-phase flow problems related to N–S equations. For two-phase flow problem, it contains both quasi-incompressible and incompressible situations, in which the quasi-incompressible model can guarantee the mass conservation, and its governing equation of the flow field can be regarded as a kind of incompressible N–S equation with a source term in the continuity equation. Actually, there are also some phase-field-based LB models for incompressible two-phase flows. He et al. [2] proposed a phase-field LB model and adopted an order parameter to track the interface of two incompressible fluids. However, there are some differences between the derived governing equations and the phase-field theory for incompressible two-phase flows [34]. To recover the Cahn–Hilliard (C–H) equation correctly, Zheng et al. [35] and Zu et al. [34] developed two different LB models, while the extra terms in the recovered macroscopic equations from their models will produce large errors in the interface capturing, and numerical instability will occur when the dimensionless relaxation time is equal to 1 [34]. To overcome these problems, Liang et al. [36] proposed a new multi-relaxation-time (MRT) LB model through introducing a time-dependent source term in the evolution equation of phase field. However, all the above models cannot conserve mass locally when the densities of the two fluids are different. To solve the problem, Yang et al. [37] presented a modified LB model based on the quasi-incompressible theory. From his model, the quasi-incompressible N–S equations in artificial compressible form can be derived. To neglect the artificial compressible effect, the model requires an additional condition, T ≫ L/cs (T and L are characteristic time and length, respectively). Based on the work of Yang et al. [37], Zhang et al. [38] proposed a discrete unified gas-kinetic scheme (DUGKS) for two-phase flows which can exactly guarantee the mass conservation, while this model also gives rise to the generation of artificial compressible effect. On the contrary, the present phase-field-based LB model can overcome the above defects. The rest of present paper is organized as follows. In Section 2 , the generalized LB model for fluid flow system with a source term is introduced, and a phase-field-based LB model for two-phase flows is given in Section 3 . Numerical experiments to validate the present model are carried out in Section 4 . Finally, some conclusions are given in Section 5. 2. Generalized LB model for fluid flow system with a source term 2.1. Governing equations The governing equations for compressible fluid flows with a source term in the continuity equation can be written as [21,39]

∂ρ + ∇ · (ρ u) = S , ∂t

(1a)

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

3

∂ (ρ u) (1b) + ∇ · (ρ uu) = − ∇ p + ∇ · τ + F, ∂t where ρ is the density, u denotes the fluid velocity, S is a source term, p is the hydrodynamic pressure, F is the external force, τ is the deviatoric stress tensor, and for Newtonian fluids, ( ) 2 τ =µ(∇ u + ∇ uT ) + ξ − µ ∇ · uI (2) d

where µ denotes the dynamic viscosity by µ = ρν , ν is the kinematic viscosity, ξ is the bulk (or volume) viscosity and d is the number of spatial dimensions. If we take S = 0 in Eq. (1), the governing equations for compressible fluid flow system can be obtained [22],

∂ρ + ∇ · (ρ u) = 0, ∂t ∂ (ρ u) + ∇ · (ρ uu) = − ∇ p + ∇ · τ + F, ∂t

(3a) (3b)

Note that the classic LB method can only be used to solve fluid flow system with low Mach numbers. Under the assumption of low Mach number, the governing equations for nearly incompressible fluid flow system can also be written as Eqs. (1) or (3). We would like to point out that our discussions in this paper are all under assumption of low Mach number, and we only consider the incompressible or nearly incompressible fluid flows. For fluid flows with small temperature changes, the flow can be regarded as incompressible under the condition of ρ = const, so that the above governing equations will become [22,26]

∇ · u = 0,

(4a)

∂u + ∇ · (uu) = − ∇ p + ∇ · [ν (∇ u)] + F. ∂t

(4b)

In this work, to simplify the following discussion and facilitate the design of a generalized model, here we consider the following generalized governing equations with a source term S,

∂ ρ˜ + ∇ · (ρ u) = S , (5a) ∂t ∂ (ρ u) (5b) + ∇ · (ρ uu) = − ∇ p + ∇ · τ + F, ∂t where ρ˜ is a physical quantity related to ρ or a constant. Note that Eqs. (5) are a kind of generalized N–S equations which contain two basic forms of incompressible and nearly incompressible governing equations, and they can be used to describe single-phase or multi-phase flows with a source term. When ρ˜ and S take different values, Eqs. (5) represent different governing equations. For example, if we take ρ˜ = ρ , Eqs. (5) will become Eqs. (1), and further reduce to Eqs. (3) with S = 0. While if we take ρ˜ = ρ = const, and S = 0, Eqs. (4) can be derived. Next, we will design the corresponding LB model for Eqs. (5), and the designed model must be able to deal with the incompressible and near incompressible single-phase and multi-phase flow problems with a source term. From this perspective, the LB model we designed is also a generalized model. 2.2. LB model for generalized N–S equations with a source term To obtain the evolution equation of Eq. (5), we integrate the following discrete velocity Boltzmann equation

∂ fi + ci · ∇ fi = Ωi + Gi ∂t

(6)

along a characteristic line over a time interval ∆t [40,41], and we have fi (x + ci ∆t , t + ∆t) − fi (x, t) =

∆t

∫ 0

Ωi (x + ci t ′ , t + t ′ )dt ′ +

∆t

∫

Gi (x + ci t ′ , t + t ′ )dt ′ .

(7)

0

where fi (x, t) denotes particle distribution function with velocity ci at position x and time t, Gi (x, t) represents the force term, Ωi (x, t) is the collision operator approximated by

Ωi = −

1

τ′

eq

(fi − fi ),

(8) eq

where τ ′ is the relaxation time and fi (x, t) is the equilibrium distribution function. Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

4

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

The integral of the collision term adopts the trapezoidal rule, then Eq. (7) becomes

∆t

fi (x + ci ∆t , t + ∆t) − fi (x, t) = ∆t

Let f¯i = fi −

2

[Ωi (x + ci ∆t , t + ∆t) + Ωi (x, t)] +

2

∆t

∫

Gi (x + ci t ′ , t + t ′ )dt ′ .

(9)

0

Ωi , we have

f¯i (x + ci ∆t , t + ∆t) − f¯i (x, t) = −

] 1 [ eq f¯i (x, t) − fi (x, t) +

∆t

∫

τg

Gi (x + ci t ′ , t + t ′ )dt ′ ,

(10)

0

+∆t ¯ ¯ where τg = 2τ2∆ is the dimensionless relaxation time, and f¯i satisfies i fi = i fi , and i ci fi = i ci fi . t Through Taylor expansion of Gi and neglecting the terms of order O(∆t 2 ), the last term on the right hand of Eq. (10) becomes ′

∆t

∫

∑

Gi (x + ci ∆t ′ , t + ∆t ′ )dt ′ =

∆t

∫

0

(

∑

∑

∑

Gi (x, t) + t ′ Di Gi (x, t) dt ′

)

0

= ∆tGi (x, t) +

∆t 2 2

(11) Di Gi (x, t),

where Di = ∂ t + ci · ∇ . The LB evolution equation with the BGK collision operator for the N–S equations can be expressed as [41]

] [ [ ] ¯fi (x + ci ∆t , t + ∆t) − f¯i (x, t) = − 1 f¯i (x, t) − f eq (x, t) + ∆t Gi (x, t) + ∆t Di Gi , i τg 2

(12)

If the up-wind scheme is used to Eq. (12), the evolution equation for the generalized N–S equations reads f¯i (x + ci ∆t , t + ∆t) − f¯i (x, t) = −

] 1 [ Gi (x + ci ∆t , t + ∆t) − Gi (x, t) eq f¯i (x, t) − fi (x, t) + ∆t Gi (x, t) + . τg 2 [

]

(13)

To remove the implicitness, we introduce a modified particle distribution function, gi (x, t) = f¯i (x, t) −

∆t 2

Gi .

(14)

With some simple manipulations, the explicit evolution equation can be derived, gi (x + ci ∆t , t + ∆t) − gi (x, t) = −

1 (

τg

eq

gi (x, t) − gi (x, t) + ∆t(1 −

)

eq

eq

1 2τg

)Gi ,

(15)

eq

where gi is the new equilibrium distribution function and satisfies gi = fi . eq To derive Eq. (5) through Chapman–Enskog analysis, the equilibrium distribution function gi is defined as (see Appendix A for the details)

{ eq gi

=

ρ˜ + cp2 (ωi − 1) + ρ si (u), i = 0 s p ωi + ρ si (u), i ̸ = 0, c2

(16)

s

with si (u) = ωi

[

ci · u cs2

+

(ci · u)2 2cs4

−

u·u 2cs2

]

,

(17)

where ωi denotes the weighting coefficient, ci is the discrete velocity, and cs is the speed of sound. Note that our model is based on the DdQq lattice with q velocity directions in d-dimensional space. ci and √ ωi depend on the choice of the lattice model, and in D1Q3 model, {ci } = (0, 1, −1)c, ω0 = 2/3, ω1,2 = 1/6, cs = c / 3, where c = ∆x/∆t, with ∆x and ∆t representing the spacing √ and time step, respectively; while in the D2Q9 model, ωi is given by ω0 = 4/9, ω1−4 = 1/9, ω5−8 = 1/36, cs = c / 3, and ci is given by

( {ci } =

0 0

1 0

0 1

−1 0

0 −1

1 1

−1 1

−1 −1

)

1 c; −1

√

in the D3Q15 model, ω0 = 2/9, ω1−6 = 1/9, ω7−14 = 1/72, cs = c / 3, and ci is given by

( {ci } =

0 0 0

1 0 0

0 1 0

0 0 1

−1 0 0

0 −1 0

0 0 −1

1 1 1

1 1 −1

1 −1 1

−1 1 1

−1 −1 −1

−1 −1

−1

1

−1

1

1 −1 c . −1

)

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

5

Different from some previous LB models [34,36,37,42–46], in the present model, the force distribution function is given by Gi = ωi

⎧ ⎨

[

ci · F

S+

cs2

⎩

+

]⎫

˜ − uuS˜ + Q (∇ · uI) ⎬ (ci ci − cs2 I) : ∂t (p − ρ˜ cs2 )I + uF˜ + Fu 2cs4

,

(18)

⎭

where F˜ is a modified total force F˜ = F − ∇ (p − ρ cs2 ),

(19)

S˜ and Q can be expressed as S˜ = S + ∂t (ρ − ρ˜ ), Q =

2 d

ρ cs2 −

(20)

ξ . ∆t(τg − 0.5)

(21)

Under the Stokes hypothesis, the bulk viscosity ξ is usually assumed to ( be zero ) [47]. In addition, ξ can also take 2ρν/d, so that the derived macroscopic equations will not contain the term ξ − 2d µ ∇ · uI. Please note that Gi is a complete form and does not contain any unnecessary assumptions, and Gi can be reasonably simplified for specific problems. We would like to point out that through the Chapman–Enskog analysis, the present LB model with the force term Eq. (18) can correctly recover Eq. (5) (see Appendix A for the details), and simultaneously, the fluid kinematic viscosity can be determined by

ν = cs2 (τg − 0.5)∆t .

(22)

In the implementation of the present model, for nearly incompressible model, the macroscopic quantities ρ˜ can be calculated as

ρ˜ =

∑

gi +

∆t 2

i

S,

(23)

where Eq. (14) is used. By taking the first-order moment of gi , the fluid velocity can be obtained [6,36], u=

1

(

) ∑

ρ

ci gi + 0.5∆tF .

(24)

i

While for the incompressible model, the macroscopic quantities we need to calculate are fluid velocity u and the pressure p. The fluid velocity is given by Eq. (24), and the pressure can be calculated as (see Appendix B for the details) p=

cs2 1 − ω0

⎡ ∑ ⎣

gi +

i̸ =0

∆t 2

⎤ S + ρ s0 (u) + ∆t(τg −

1

S + ρ s0 (u) + ∆t(τg −

1

2

)G0 ⎦ ,

(25)

or p=

cs2 1 − ω0

⎡ ∑ ⎣

gi +

i̸ =0

∆t 2

⎤ 2

)(S −

∑

Gi )⎦ .

(26)

i̸ =0

We would also like to point out that the present pressure expression is different from that in Refs. [36,45] which can be expressed as p=

cs2 1 − ω0

⎡ ∑ ⎣ i̸ =0

gi +

∆t 2

⎤ S + ρ s0 (u)⎦ ,

(27)

where the term of ∆t(τg − 12 )G0 has been neglected. However, the last term ∆t(τ − 21 )G0 may be significant since G0 is usually nonzero under the condition of S ̸ = 0, and the effect of this item cannot be ignored. Thus, our generalized LB model for fluid flow system is made up of Eqs. (15), (16), (18) and (23), (24) or (24), (25). In the derivation of the model, we only used the assumption of low Mach number (Ma ≪ 1), which is necessary for the construction of most LB models. Under this assumption, the following equations are established for nearly incompressible or incompressible flows, u = O(Ma), δρ = O(Ma2 ), and δ p = O(Ma2 ). Further, we have F = O(Ma) and S = O(Ma). Now, we will use some remarks to show that our model not only contains some existing models, but also deduces new models for N–S equations with a source term. Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

6

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Remark 1. When taking ρ˜ = ρ and p = ρ cs2 , the generalized model can reduce to the nearly incompressible form. Correspondingly, the macroscopic equations become Eqs. (1), while the equilibrium distribution function can be written as eq

gi = ρ [ωi + si (u)] ,

(28)

and the force distribution function is given by

{

Gi = ωi S +

ci · F cs2

+

(ci ci − cs2 I) : [uF + Fu − uuS + Q (∇ · uI)]

}

2cs4

,

(29)

where the term uuS is the order of O(Ma3 ) and can be neglected. This model is almost identical to the one in Ref. [21], except that the expression of the momentum equation and the value of bulk viscosity are different. Remark 2. If we rewrite the equilibrium distribution function Eq. (28) to the following form,

{

eq

gi = ωi ρ + (ρ0 + δρ )

[

ci · u cs2

+

(ci · u)2 2cs4

−

u·u

]}

2cs2

.

(30)

Multiplying by cs2 on both sides of Eq. (30) and ignoring the term of O(Ma3 ), one can get eq gi

{

= ωi p + p0

[

ci · u cs2

+

(ci · u)2 2cs4

−

u·u 2cs2

]}

,

(31)

eq

where we still use gi to represent the new equilibrium distribution function. The force distribution function is given by Gi = ωi

⎧ ⎨

[

p0 S1 +

p0 ci · F1 cs2

⎩ where p0 = const, S1 =

S , p0

F1 =

(ci ci − cs2 I) : p0 (uF1 + F1 u) +

+ F , p0

cs2 ∂ t

2 p c2 d 0 s

−

ξ ∆t(τg −0.5)

2cs4

)

]⎫

(∇ · uI) ⎬

,

(32)

⎭

and the terms of O(Ma3 ) are abandoned in the incompressible limit. The fluid velocity

∆t and pressure can be obtained by p = i gi + 2 p0 S1 and p0 u = equations with a source term S1 can be obtained from Eqs. (5),

1 ∂P

(

∑

∑

i ci gi

+

∆t 2

p0 F1 . The corresponding incompressible N–S

+ ∇ · u = S1 ,

(33a)

∂u 2 + ∇ · (uu) = − ∇ P + ∇ · [ν (∇ u + ∇ uT )] + ∇ · [(− ν )∇ · uI] + F1 , ∂t 3

(33b)

where P = ρ . In the limit of a low Mach number (Ma = |u|/cs2 ), the dynamic pressure is assumed to be δ p ∼ Ma2 and 0 the left end term of Eqs. (33a) can be ignored. Note that if we take S1 = 0, this model will reduce to the incompressible LB model by He et al. [25]. p

Remark 3. When taking ρ˜ = ρ0 , ρ = const (e.g., ρ = 1), where ρ0 is the average velocity of the fluid, we can derive a new model for incompressible N–S equations with a source term in the continuity equation. In this case, the macroscopic equations can be written as [see Eqs. (5)]

∇ · u = S,

(34a)

∂u 2 + ∇ · (uu) = − ∇ p + ∇ · [ν (∇ u + ∇ uT )] + ∇ · [(− ν )∇ · uI] + F. ∂t 3

(34b)

The corresponding equilibrium distribution function can be derived from Eq. (16)

{ eq gi

=

ρ0 + cp2 (ωi − 1) + si (u), i = 0 s p ωi + si (u), i ̸ = 0, c2

(35)

s

and the force distribution function is given by

{

Gi = ωi S +

ci · F cs2

+

(ci ci − cs2 I) : [uF + Fu + QSI] 2cs4

}

,

(36)

where the terms uuS˜ and ∂t p are neglected in the incompressible limit. We note that if we take S = 0, this model can reduce to the incompressible LB model by Guo et al. [26] with a force term, i.e., Eqs. (4). Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

7

Remark 4. Eq. (24) may be implicit if the force F is a nonlinear function of u. To remove this implicitness, we can discretize Di Gi by [41] Di Gi =

Gi (x + ci ∆t , t) − Gi (x, t − ∆t)

∆t

.

(37)

Then the evolution equation can be rewritten as gi (x + ci ∆t , t + ∆t) − gi (x, t) = −

1 (

τg

eq

[

gi (x, t) − gi (x, t) + ∆t Gi (x, t) +

)

Gi (x + ci ∆t , t) − Gi (x, t − ∆t) 2

]

.

(38)

This scheme is explicit and can be iterated if Gi is known, while the results of Gi at time t − ∆t must be saved additionally. 3. Phase-field-based LB model for two-phase flows In Section 2, we have given the generalized LB model for single-phase flows with a source term. Actually, one of the motivations of the generalized model is to deal with the two-phase flow problems based on phase-field theory. 3.1. Governing equations In the classic theory of phase field model for two-phase flows, the N–S equations are employed to describe the flow field, while the C–H equation is usually adopted to capture the phase interface which can be given by [48–50]

∂φ + ∇ · φ u = ∇ · Mφ (∇µ) , ∂t

(39)

where Mφ is the mobility coefficient, φ is the order parameter defined as the volume fraction of one of the two phases, and φ and ρ satisfy the linear relationship,

ρ − ρB φ − φB = . (40) ρA − ρB φA − φB In this work, φ = 0 denotes the phase B while φ = 1 represents the phase A, and µ is the chemical potential, which is derived by the variation of the free-energy function F (φ ) with respect to the order parameter [48,51,52], δ F (φ ) dψ (φ ) µ= = − κ∇ 2 φ δφ dφ (41) = 4βφ (φ − 1)(φ − 0.5) − κ∇ 2 φ. where F (φ ) can be taken as the following form [48–50,53], F (φ ) =

∫ [ ] κ ψ (φ) + |∇φ|2 dΩ , 2

Ω

(42)

where Ω is the physical domain of the system, ψ (φ) denotes the bulk free-energy density, and k |∇φ|2 /2 accounts for the surface energy with a positive coefficient k. If the system considered is a van der Waals fluid, the bulk free energy has a double-well form [48],

ψ (φ ) = βφ 2 (φ − 1)2 ,

(43)

where β is a constant dependent on the interfacial thickness W and the surface tension σ [48],

√ W = and

8κ

(44)

β √

σ =

2κβ

. (45) 6 The equilibrium interface profile can be obtained by minimizing F (φ ) with respect to the function φ , i.e., µ = 0. In a plane interface at the equilibrium condition, the order parameter profile across the interface (along the z direction) is represented by [52] φ (z) =

1 2

+

1 2

( tanh

2z W

)

.

(46)

In the present work, we will focus on the following quasi-incompressible phase-field system, and the governing equations are described as

∂φ + ∇ · φ u = ∇ · Mφ (∇µ) , ∂t

(47a)

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

8

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

∇ · u = S1 ,

(47b)

∂ (ρ u) + ∇ · (ρ uu) = − ∇ p + ∇ · τ + F, ∂t

(47c)

where F represents the total force which is defined as Fs + G, and G is the body force, Fs is the surface tension with the potential form Fs = µ∇φ [48,54] if not specified. Eq. (47b) can be derived from Eq. (5a), if S = u · ∇ρ + ρ S1 , and ρ˜ = ρ0 . 3.2. The LB model for quasi-incompressible phase-field system 3.2.1. LB model for the N–S equations Based on the generalized LB model for fluid flow system, one can get the equilibrium and force distribution function for the N–S equations [Eqs. (47b) and (47c)] when substituting ρ˜ = ρ0 , and S = u · ∇ρ + ρ S1 into Eqs. (16) and (18),

{ ρ0 + cp2 (ωi − 1) + ρ si (u), i = 0, s = p ωi + ρ si (u), i ̸ = 0, cs2 ( )⎤ ⎡ 2 ˜ + Fu ˜ + 2 ρ cs2 S1 I (c c − c I) : u F i i s d ci · F ⎦, Gi = ωi ⎣u · ∇ρ + ρ S1 + 2 + 4 eq gi

cs

2cs

(48)

(49)

where we take the bulk viscosity ξ equal to 0 and the dynamic pressure satisfies δ p = O(Ma2 ), so that the term of ∂t p is neglected in the limit of a low Mach number. Actually, the relationships of u = O(Ma) and F = O(Ma) are also true. Specially, the force distribution function G0 can be simplified as

( u · ∇ρ + ρ S1 −

G0 = ω0

2u · F˜ + 2d ρ cs2 S1

)

2cs2

( )] [ 1 2 1 = ω0 ρ S1 − 2 u · F − u · ∇ p + ρ cs S1 cs d { } 1 1 = ω0 (1 − )ρ S1 − 2 [u · (F − ∇ p)] . d

(50)

cs

From the above equation, one can find that G0 is related to S1 since the term u · (F − ∇ p) is the order of O(Ma2 ), and is nonzero once S1 ̸ = 0, thus the term ∆t(τg − 21 )G0 in Eq. (25) cannot be ignored in the pressure expression. We would like to point out that Eq. (47) will reduce to the incompressible phase-field model if S1 = 0. However, based on Eqs. (47a), (47b), (40) and S1 = 0, one can obtain [37]

ρA − ρB ∂ρ + ∇ · (ρ u) = ∇ · [Mφ ∇µ]. ∂t φA − φB

(51)

It is clear that the mass conservation is constrained by the term on the right hand of Eq. (51), which is nonzero in the interfacial region as long as ρA ̸ = ρB . Therefore, in the incompressible phase-field model, the mass is not locally conserved. To conserve mass locally, Shen et al. [33] proposed a quasi-incompressible phase-field model. Subsequently, based on this quasi-incompressible phase-field model, Yang et al. [37] designed a lattice Boltzmann model for binary fluids. Actually, Eq. (47) can also reduce to quasi-incompressible phase-field model in Ref. [33] when S1 = −γ ∇ · (Mφ ∇µ) with −ρB γ = φ ρρA −φ , where the quasi-incompressible model can conserve the mass. It should also be noted that the momentum ρ A B

B A

equation (47c) is different from those used in Refs. [36,37], where the terms of uS and

[( ) ] ξ − 2d µ∇ · u I are not included.

Remark between the generalized LB model and the one in Ref. [37]. If we take [ 5. Here we also give a] comparison p S = cs2 u · ∇ρ − ργ ∇ · (Mφ ∇µ) , and ρ˜ = 2 , the equilibrium distribution function for the N–S equations can be written cs as eq

gi = ωi

[

p cs2

] + ρ si (u)

(52)

Based on the following fact, p cs2

=

∑

eq

gi =

∑

i

gi +

i

∆t ∑ 2

i

Gi ,

∑

Gi = S ,

(53)

i

the pressure can be calculated as p=

∑ i

cs2 gi +

∆t 2

cs2 u · ∇ρ − ργ ∇ · (Mφ ∇µ) ,

[

]

(54)

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

9

and the fluid velocity can be obtained from Eq. (24). The above derivation shows that our generalized flow field LB model can reduce to Yang’s model once the force distribution function Gi takes the same form, and Yang’s model is in fact a nearly incompressible model. 3.2.2. LB model for the C-H equation The evolution equation for the C-H equation can be given as hi (x + ci ∆t , t + ∆t) − hi (x, t) = −

] 1 [ eq hi (x, t) − hi (x, t) + ∆tRi (x, t),

τh

(55)

where hi (x, t) is the distribution function of order parameter φ , τh is the non-dimensional relaxation time related to the eq mobility, hi (x, t) is the local equilibrium distribution function, which is defined as [55–57]

⎧ ⎨φ + (ω − 1)ηµ, i = 0 i eq hi (x, t) = ⎩ωi ηµ + ωi ci ·φ2 u , i ̸= 0, c

(56)

s

where η is an adjustable parameter that controls the mobility. Ri (x, t) is the source term and is given by [36]

( Ri =

1−

1

)

2τ

ωi ci · ∂t φ u cs2

.

(57)

In our simulations, the first-order explicit difference scheme is used to compute the time derivative in Eq. (57),

∂t φ u|(x,t) = [φ u|(x,t) −φ u|(x,t −∆t) ]/∆t .

(58)

Through the Chapman–Enskog analysis, the order parameter φ is calculated by

φ=

∑

hi ,

(59)

i

the C–H equation can be recovered with second-order accuracy and the mobility can be determined by Mφ = ηcs2 (τh − 0.5)∆t .

(60)

In a two-phase system, the viscosity is no longer a uniform value due to its jump at the interface. In this work, the following viscosity with an inverse linear form [44,46] is adopted to ensure a smooth viscosity across the interface if not specified, 1

ν

=

φ − φB 1 1 1 ( − )+ . φA − φB νA νB νB

(61)

In addition, to determine the gradient terms in the source term Gi , surface tension Fs and chemical potential µ, the following isotropic schemes are adopted to discretize the first-order spatial derivative and the Laplacian operator [36]:

∇ζ (x, t) =

∑ ωi ci ζ (x + ci ∆t , t) i̸ =0

∇ 2 ζ (x, t) =

cs2 ∆t

,

∑ 2ωi [ζ (x + ci ∆t , t) − ζ (x, t)] i̸ =0

cs2 ∆t 2

(62a)

,

(62b)

where ζ is an arbitrary function. It should be noted that the schemes (62) not only have a secondary-order accuracy in space, but also can ensure the global mass conservation of a multiphase system [6]. 4. Model validation In this section, several examples including a spinodal decomposition, a static droplet, layered Poiseuille flows and a bubble rising flow, are performed to test our LB model for incompressible (S1 = 0) and quasi-incompressible [S1 = −γ ∇ · (Mφ ∇µ)] two-phase flows. In our simulations, the D2Q9 lattice structure is adopted. We attempt to conduct a detailed comparison between the present and analytical solutions or some available results in each test. 4.1. Spinodal decomposition Spinodal decomposition [58] is a mechanism for the rapid unmixing of a fluid mixture with two different species. The spinodal decomposition phenomenon will take place when the initial homogeneous mixture is unstable in the presence of small fluctuations. In this section, the separation of two emulsified fluids is simulated with different pressure expressions, i.e., Eqs. (25) and (27), and this example is mainly used to demonstrate the differences in calculation results obtained by Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

10

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 1. (Color online) Separation of binary fluid: distribution of order parameter with modified pressure expression equation (25), (a) our quasi-incompressible model, (b) our incompressible model.

present incompressible and quasi-incompressible models. The computational domain is NY × NX = 200 × 200, the initial distribution of order parameter with a small fluctuation is given by

φ (x, y) =

1

+ rand(x, y), (63) 3 where rand(x, y) is a random function with the maximum amplitude of 1%. The initial velocity is zero in the whole domain and the periodic boundary conditions are applied at all boundaries. The model parameters are fixed as σ = 0.03, W = 4, Mφ = 0.1, φA = 1, φB = 0, ρA /ρB = 5, νA /νB = 0.1, τh = 1. Fig. 1 shows the time evolution of the order-parameter distribution by our incompressible model (IM) and quasi-incompressible model (QIM) with Eq. (25). It can be found from this figure that small fluctuations of order-parameter gradually become larger, then some small droplets are formed, and the diameters of the droplets gradually become larger. Finally, the phenomenon of fluid–fluid separation can be observed as expected. However, the results of the incompressible and quasi-incompressible models are significantly different. To quantitatively show the difference between the two models, we calculated the difference between the order parameters obtained by the two models at different times (see Fig. 2). One can observe that the difference (φIM − φQIM ) between the order parameters gradually increases with time, where φIM and φQIM represent the order parameters calculated by the incompressible and quasi-incompressible models, respectively. This difference may be caused by the fact that the term γ ∇ · (Mφ ∇µ) has been neglected in incompressible model, and thus mass conservation is not satisfied strictly. To illustrate the difference between two pressure expressions, we presented some results in Fig. 3. As seen in Fig. 3(a), the distributions of order parameter within the red circles are distinctly different, which means that the term ∆t(τg − 0.5)G0 in Eq. (25) has a significant influence on the numerical results. However, for incompressible LB model, there are no apparent differences, which is due to the fact that the effect of ∆t(τg − 0.5)G0 can be ignored when S1 = 0, and G0 is the order of O(Ma2 ). Based on the above results, one can conclude that the pressure expression equation (25) is more general, and would be adopted in the following simulations. 4.2. Static droplet 2D static droplet test is a popular problem to verify LB models for two-phase flow [34,36,43,55,59]. In this subsection, we will consider this problem with different density ratios. Initially, a circular droplet with the radius ranging from 20 to 40 is placed in the middle of the computational domain with NX × NY = 100 × 100. The initial order parameter is given by

[

2 R−

√

(x − 50)2 + (y − 50)2

]

, (64) W where R is the droplet radius, and surface tension is expressed as Fs = −φ∇µ. In numerical simulations, we set the density ratio to be ρA /ρB = 2, 10, and 50. The other physical parameters are fixed as ρB = 1, τh = 1, τg = 0.8, σ = 0.001, Mφ = 0.02, and the periodic boundary conditions are applied at all boundaries. We first verify the LB model with the well-known Laplace’s law, φ (x, y) = 0.5 + 0.5 tanh

∆P =

σ R

,

(65)

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

11

Fig. 2. (Color online) The difference between the order parameters obtained by the two models at different times: (a) t = 0, (b) t = 2000, (c) t = 6000, (d) t = 20000, where the z-axis represents (φIM − φQIM ). Table 1 Relative error of pressure jump with the density ratio ρA /ρB = 50. R

Present QIM

Present IM

Yang et al.

Liang et al.

20 25 30 35 40

0.02% 0.45% 0.42% 0.53% 0.68%

0.36% 0.73% 0.84% 0.95% 1.00%

0.50% 0.63% 0.63% 0.60% 0.48%

0.36% 0.73% 0.84% 0.95% 1.00%

where ∆P is the pressure jump across the interface, R is the radius of the droplet. P is calculated by P = p0 − κφ∇ 2 φ + κ|∇φ|2 /2 + p with the equation of state p0 = φ∂φ ψ − ψ [34,60]. We performed some simulations and presented the results in Fig. 4. From Fig. 4(a), (c), (e) one can find that the results of present models and those in Refs. [36,37] agree well with the Laplace law. In order to show the difference between these models more clearly, we also give the relative errors of pressure jump with the density ratio ρA /ρB = 50 in Table 1. Here the relative error of pressure jump is given by Error =

|∆Pn − ∆Pa | , |∆Pa |

(66)

where ∆Pa is the analytical pressure jump computed by Laplace’s law, ∆Pn = PA − PB is the numerical pressure jump with PA and PB the pressures of fluids A and B. In general, the present quasi-incompressible LB model and quasi-incompressible model in Ref. [37] are more accurate than incompressible LB models. This is because that quasi-incompressible model is physically more reasonable (it can satisfy the conservation of mass) than incompressible model. Besides, it is also found that usually present QIM is more accurate than the model in Ref. [37]. In addition, we also presented the density profiles along the horizontal center line at t = 5 × 105 ∆t in Fig. 4 (b), (d), and (f). From these figures, one can observe that all the numerical results are very close to the analytical solutions given Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

12

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 3. (Color online) Distribution of order parameter with Eq. (25) on the left, and Eq. (27) on the right, (a) the quasi-incompressible model at t = 20 000, (b) the incompressible model at t = 20 000. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

by Eq. (64). To give a quantitatively estimation on the accuracy of numerical results, we also measured the relative errors of density, i.e., [(ρ − ρ0 )/ρ0 ] with ρA /ρB = 10 in Fig. 5(a), where ρ0 is the analytical solution. Different from the results of Laplace’s law, we can find that the present incompressible LB models and the one in Ref. [36] produce smaller errors than quasi-incompressible models, which means that the quasi-incompressible and incompressible models each has its own advantages. 4.3. Layered Poiseuille flow The layered Poiseuille flow between two parallel plates is also a classical two-phase problem which provides a good benchmark for validating the LB models [34,46,59,61,62]. Considering a channel flow of two immiscible fluids driven by a const pressure gradient G in the flowing direction (x-direction). Initially, fluid A is located in the upper region of the channel (0 < y ≤ h), while fluid B is at the bottom region (−h ≤ y ≤ 0). When the flow is sufficiently slow and no instabilities occur at the interface, an steady analytical solution of velocity can be obtained,

( ) ⎧ 2 [ ( )2 −µB ⎨ 2Ghµ − hy − hy µµA +µ + A A B ( [ ( ) ) ux,a (y) = 2 ⎩ Gh − y 2 − y µA −µB + 2µB h h µA +µB

2µA

]

,

0
2µB

]

.

−h≤y≤0

µA +µB µA +µB

(67)

In this work, G is given as G = uc (µA + µB )/h , and to ensure the stability of the interface, uc is fixed as uc = 5 × 10−5 . To quantitatively describe the accuracy of the present models and also convenient compare with the existing LB models, the following relative error is adopted, 2

∑ Error =

y

|unx (y, t) − uax (y)| ∑ a , y |ux (y)|

(68)

where the subscripts a and n denote the analytical and numerical solutions. In our simulations, the computational domain is chosen as NY × NX = 100 × 10. Periodic boundary conditions are applied in the x-direction, and the halfway bounce-back boundary conditions are enforced on the top and bottom walls. We first investigated the effects of viscosity ratio. To this end, four different cases with viscosity ratios of µA /µB = 3, 10, 100, 1000 are considered. The other parameters are given as W = 4, σ = 0.001, Mφ = 0.1, ρl /ρg = 1. Based on the results in Fig. 6, one can find that the numerical results of present LB models and some others [36,37] agree well with the analytical solutions for different viscosity ratios. We also calculated the relative errors of velocity under different viscosity ratios, and the results are given in Table 2. From this table, one can observe that the relative error increases as the viscosity ratio becomes larger, and for a fixed viscosity ratio, the relative errors of all these models are almost the same. This is Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

13

Fig. 4. (Color online) Tests of Laplace’s law [(a), (c), (e)] and density profiles [(b), (d), (f)] at different density ratios (ρA /ρB = 2, 10, 50), where the droplet radius is 25.

because if the density ratios are equal to 1, which leads to γ = 0, the quasi-incompressible and incompressible models are equivalent except for some terms with the order of O(Ma2 ). We then simulated the layered Poiseuille flow with another density ratio ρA /ρB = 3, and viscosity ratio is equal to 3. Here the dynamic viscosity is given by [45]

µ=

{ µA , µB ,

φ ≥ 0.5, φ < 0.5.

(69)

The other parameters are the same as those stated above. From Fig. 7(a), one can find that there is a good agreement between the numerical solutions of the four LB models and the analytical solutions except in the interface region. Fig. 7(b) is an enlarged view of the interface region in Fig. 7(a). From this figure, it can be observed that the present QIM and IM produce smaller errors than the models of Yang et al. [37] and Liang et al. [36] in the interface area, while their models Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

14

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 5. (Color online) Tests of diffusion errors of density profiles with ρA /ρB = 10 and R = 25.

Fig. 6. (Color online) Comparison of the velocity distributions obtained by present QIM, IM, model of Yang et al. [37], and model of Liang et al. [36] with the corresponding analytical solutions (solid line): (a) µA /µB = 3, (b) µA /µB = 10, (c) µA /µB = 100, (d) µA /µB = 1000. ux is normalized by the maximum speed of analytical solution umax .

are more accurate in the bulk region. In other words, our models have a better performance in the interface region, in contrast, the models of Yang et al. [37] and Liang et al. [36] are more accurate in the bulk region. 4.4. Bubble rising under buoyancy To further demonstrate the accuracy of the present models for more complex flows, the problem of single bubble rising flow driven by the buoyancy is also considered here. Initially, a circular bubble (fluid B) without initial velocity is immersed in the bottom center of another fluid (fluid A). The radius of the initial bubble, R, occupies 32 lattice spaces. To generate the buoyancy effects, a body force, Fb,y = −(ρ − ρA )g, is added to the momentum equation, where g is the gravitational acceleration. We conducted some simulations on a uniform computational mesh with the size of Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

15

Table 2 Relative errors of velocity for different viscosity ratios. Models

µA µB

Present QIM Present IM Yang et al. [37] Liang et al. [36]

1.04% 1.04% 1.04% 1.04%

=3

µA µB

= 10

1.30% 1.30% 1.30% 1.30%

µA µB

= 100

1.90% 1.90% 1.90% 1.90%

µA µB

= 1000

2.16% 2.16% 2.16% 2.16%

Fig. 7. (Color online) Comparison of the velocity distributions obtained by present QIM, IM, model of Yang et al. [37], and model of Liang et al. [36] with the analytical solution (solid line), where ρA /ρB = 3, and µA /µB = 3: (a) velocity distribution in the whole region, (b) velocity distribution in the interface region.

NX × NY = 160 × 480, and the periodic boundary conditions are applied at all boundaries. The other related parameters are given as W = 4, g = 10−5 , σ = 0.001, τf = τg = 1, Mφ = 6.667, ρA /ρB = 2. Fig. 8 shows the density distributions of the rising bubble at different times based on our models and models of Yang et al. [37] and Liang et al. [36]. From this figure, it can be observed that the results of these models are quite similar to each other. However, a zoom-in-view of Fig. 8 at t = 40 000 shows that the incompressible models (our IM and Liang’s model) will produce numerical oscillation near the interface region (see Fig. 9). To see the differences among these LB models more clearly, we show the errors along the vertical centerline at t = 40 000 in Fig. 10, where the error between our QIM and IM increases sharply to 0.096 near the interface. This phenomenon is reasonable because the compressible term γ ∇ · (Mφ ∇µ) has an influence on the numerical results. However, the maximum error between our QIM and the one in Ref. [37] is approximately 0.028, which is caused by the difference between the two models. In the model of Yang et al. [37], the N–S equations in artificial compressible form can be derived. To neglect the artificial compressible effect, the condition T ≫ L/cs should be satisfied. From present QIM, the quasi-incompressible N–S equations can be exactly recovered in the limit of a low Mach number. Fig. 11 depicts the pressure distributions of the rising bubble at t = 40 000 based on different models. From this figure, one can see that there is a significant difference between quasi-incompressible models and incompressible models, and the pressure interface of quasi-incompressible models is more clear. Based on above observations, one can conclude that quasi-incompressible models (our QIM and Liang’s model) are more superior for complex two-phase flows. 5. Conclusions In this study, to solve single-phase flow problems with a source term in the governing equations and other problems coupled with the flow field, such as two-phase flow problems, we developed a generalized LB model for incompressible and nearly incompressible N–S equations with a source term in the continuity equation. The model we propose has the following characteristics. First, it can apply to single-phase flow and two-phase flow problems. Second, it is a generalized model. For single-phase flow, it contains incompressible and nearly incompressible situations; for two-phase flows, it includes common incompressible and quasi-incompressible situations. Third, it modifies the existing models and one can derive new models from it. From the derived models, we can get an incompressible model for N–S equations with a source term, and we present a modified scheme to calculate the pressure p, which is more accurate than the previous one. Simultaneously, our generalized model can recover the macroscopic equations without any unnecessary assumptions. Based on the generalized LB model, a new LB model is proposed for the quasi-incompressible and incompressible phase-field systems. To validate the accuracy of the proposed model, a series of numerical tests were performed. First, we conducted a detailed comparison between present scheme and the previous one [36,45]. The result shows that there is a significant difference between the two pressure schemes when S1 ̸ = 0, and theoretically, the present pressure Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

16

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 8. (Color online) Density distributions of the rising bubble at t = 10 000, 20 000, 30 000, 35 000, 40 000, (a) present QIM, (b) present IM, (c) model of Yang et al. [37], and (d) model of Liang et al. [36].

Fig. 9. (Color online) The enlarged images of density distributions of the rising bubble at t = 40 000 for present QIM, present IM, model of Yang et al. [37], and model of Liang et al. [36].

Fig. 10. (Color online) Errors along the vertical centerline at t = 40 000, where ρ0 is the result of our QIM.

scheme is more accurate. Then we investigated two basic steady problems of a static droplet and layered Poiseuille flows. The results of the former case show that present quasi-incompressible model usually performs better than the quasiincompressible model of Yang et al. [37] in terms of accuracy, and quasi-incompressible models are more accurate than incompressible models. In the latter case, we simulated the layered Poiseuille flows with ρA /ρB = 1 and ρA /ρB ̸ = 1, and Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

17

Fig. 11. (Color online) Pressure distributions of the rising bubble at t = 40 000 for (a) present QIM, (b) present IM, (c) model of Yang et al. [37], and (d) model of Liang et al. [36].

found that our models can obtain satisfactory results in the velocity under different viscosity ratios, and our models have a better performance in the interface region when ρA /ρB ̸ = 1. Finally, we carried out some simulations of single bubble rising flow driven by the buoyancy to further demonstrate the accuracy of the present models. The results indicate that the incompressible LB models will produce numerical oscillation near the interface region. While, the quasi-incompressible LB models seem more reasonable from physical point of view, and should be considered in the study of complex two-phase flows. Acknowledgments The authors are grateful to referees for their valuable comments and suggestions. This work is supported by the National Natural Science Foundation of China (Grant Nos. 51576079, 51836003), and the National Key Research and Development Program of China (Grant No. 2017YFE0100100). Appendix A. Chapman–Enskog analysis of the present model In Appendix A, we would present the details on how to obtain the proposed LB model for hydrodynamic equations [Eq. (5)] through the Chapman–Enskog (C–E) expansion. Before performing C–E expansion, we first define the zeroth to second moments of the equilibrium distribution function,

∑ i

eq

gi = M0 ,

∑ i

eq

ci gi = M1 ,

∑ i

eq

ci ci gi = M2 ,

∑

eq

ci ci ci gi = M3 .

(A.1)

i

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

18

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

In the C–E analysis, the time and space derivatives, the force and source term can be expanded as, (0)

+ ϵ gi(1) + ϵ 2 gi(2) + · · · ,

gi = gi

(1)

(A.2a)

+ ϵ 2 G(2) i ,

(A.2b)

∂t = ϵ∂t1 + ϵ 2 ∂t2 , ∂α = ϵ∂1α ,

(A.2c)

Fα = ϵ Fα(1) + ϵ 2 Fα(2) ,

(A.2d)

S = ϵ S (1) + ϵ 2 S (2) ,

(A.2e)

Gi = ϵ Gi

where ϵ is a small expansion parameter and Greek indices denote Cartesian spatial components. Using the Taylor expansion to Eq. (15), we have

∆t 2

∆tDi gi (x, t) +

D2i gi (x, t) + · · · = −

2

) 1 1 ( eq gi (x, t) − gi (x, t) + ∆t(1 − )Gi , 2τg

τg

(A.3)

where Di = ∂t + ciα ∂α , and substituting Eqs. (A.2) into (A.3), one can obtain the following multi-scale equations, (0)

O(ϵ 0 ) : gi

= gieq , (0)

O(ϵ 1 ) : D1i gi

(0)

O(ϵ 2 ) : ∂t2 gi

(A.4a)

=−

1

(1)

gi

τg ∆t

+ (1 −

∆t

+ D1i gi(1) +

2

1

(1)

)Gi ,

2τg

(0)

D21i gi

=−

1

τg ∆t

(A.4b) (2)

+ (1 −

gi

1

(2)

2τg

)Gi ,

(A.4c)

where D1i = ∂t1 + ciα ∂1α . Then, the substitution of Eqs. (A.4b) into (A.4c) yields

∂t2 gi(0) + (1 −

1

(1)

)D1i gi

2τg

+

∆t 2

1

(1 −

(1)

2τg

)D1i Gi

=−

1

(2)

τg ∆t

gi

+ (1 −

1

(2)

2τg

)Gi .

(A.5)

By summing Eqs. (A.4b) and (A.4b)×ciβ over i, the recovered equations at ϵ scale can be obtained,

∂t1 M0 + ∂1α M1α = −

1

∑

τg ∆t

∂t1 M1β + ∂1α M2αβ = −

(1)

gi

+ (1 −

i

1

∑

τg ∆t

(1)

ciβ gi

1 2τg

)

(1)

Gi ,

(A.6a)

i

+ (1 −

i

∑ 1 2τg

)

∑

(1)

ciβ Gi .

(A.6b)

i

Similarly, we can also obtain the recovered equations at ϵ 2 scale from Eq. (A.5)

∂t2 M0 + (1 −

[

1

]

) ∂t1 (

2τg

∑

+∂1α (

(1) ciα Gi )

=−

i

∂t2 M1β + (1 −

1 2τg

+∂1α (

∑

(1) ciα gi )

1

τg ∆t

i

(2) gi

+ (1 −

[

1 2τg

)

] ∑

) ∂t1 (

(1) ciβ gi )

+ ∂1α Λ(1) +

i

(1) ciα ciβ Gi )

=−

i

∆t

+

2

i

∑

] ∑

+ ∂1α (

i

] ∑

(1) gi )

1

∑

τg ∆t

i

(2) ciβ gi

+ (1 −

∑

(1 −

(2) Gi

[

1

∑

2τg

) ∂t1 (

(1)

Gi )

i

(A.7a)

,

i

∆t 2

(1 −

1 2τg

)

1 2τg

∑

[ ) ∂t1 (

(2) ciβ Gi

∑ i

(1)

ciβ Gi ) (A.7b)

,

i

(1)

where Λ(1) = is the first-order momentum flux tensor. i ciα ciβ gi Summing Eqs. (14) and (14)×ciα over i, one can obtain

∑

M0 =

∑

gi +

i

M1α =

∑ i

∆t ∑ 2

ciα gi +

Gi ,

(A.8a)

i

∆t ∑ 2

ciα Gi ,

(A.8b)

i

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

which can be further recast as ∑ (1) ∆t ∑ (1) ∑ (2) ∆t ∑ (2) Gi , Gi , gi = − gi = − 2 2 i i i i ∑ ∑ ∆t ∑ ∆t ∑ (1) (1) (2) (2) ciα gi = − ciα Gi , ciα Gi . ciα gi = − 2 2 i

i

19

(A.9a) (A.9b)

i

i

Substituting Eqs. (A.9) into (A.6) and (A.7), we have

∑

∂t1 M0 + ∂1α M1α =

(1)

Gi ,

(A.10a)

i

∑

∂t1 M1β + ∂1α M2αβ =

(1)

ciβ Gi ,

(A.10b)

i

∂t2 M0 =

∑

(2)

Gi ,

(A.11a)

i

∂t2 M1β + (1 −

1 2τg

)∂1α Λ(1) +

∆t 2

(1 −

1 2τg

)∂1α (

∑

(1)

ciα ciβ Gi ) =

∑

i

(2)

ciβ Gi .

(A.11b)

i

Combining Eqs. (A.10a) and (A.11a) at ϵ and ϵ 2 scales yields

∂t M0 + ∂1α M1α =

∑

Gi .

(A.12)

i

To recover the continuity equation with a source term [Eq. (5a)], the following conditions should be satisfied M0 = ρ, ˜ M1 = ρ u,

∑

Gi = S .

(A.13)

i

In addition, to recover the momentum equation [Eq. (5b)], M2 = pI +ρ uu is also needed. Thus, the equilibrium distribution function can be given as eq

gi =

⎧ ⎨ρ˜ +

p cs2

(ωi − 1) + ρ si (u),

i = 0,

⎩ ωi + ρ si (u), c2

(A.14)

i ̸ = 0,

p s

and Eq. (A.10b) can be rewritten as

∂t1 (ρ uβ ) + ∂1β p + ∂1α (ρ uα uβ ) =

∑

(1)

ciβ Gi .

(A.15)

i

To derive the equation at ϵ 2 scale, we first express Λ(1) as

[ Λ

(1)

= −τg ∆t ∂t1 M2αβ + ∂1γ M3αβγ − (1 −

]

1 2τg

)

∑

(1) ciα ciβ Gi

i

{ = −τg ∆t ∂t1 pδαβ + ∂t1 (ρ uα uβ ) + cs2 ∂1γ (ρ uγ δαβ ) + cs2 ∂1γ [ρ (uαδβγ + uβ δαγ )] } 1 ∑ (1) −(1 − ) ciα ciβ Gi , 2τg

(A.16)

i

where Eq. (A.4b) has been used, and the term ∂t1 (ρ uα uβ ) can be given by

∂t1 (ρ uα uβ ) = uα (

∑

(1)

ciβ Gi ) + uβ (

i

∑

(1)

ciα Gi ) − (uα ∂1β p + uβ ∂1α p) − uα uβ S˜ (1) ,

(A.17)

i

where Eqs. (A.15) and (A.10a) are used, S˜ (1) = S (1) + ∂t1 (ρ − ρ˜ ) and the term of O(Ma3 ) has been neglected. Combining Eqs. (A.15) with (A.11b) at ϵ and ϵ 2 scales, together with Eqs. (A.16) and (A.17), we now obtain

[ ] [ ∂t (ρ uβ ) + ∂α (ρ uα uβ ) = − ∂β p + ∂α ρν (∂α uβ + ∂β uα ) + ϵ ∆t(τg − 0.5)∂α ∂t1 pδαβ ∑ ∑ (1) (1) ciβ Gi ) + uβ ( ciα Gi )+ + ∂1γ (cs2 ρ uγ δαβ ) + uα ( i

i

cs2 (uα ∂1β ρ + uβ ∂1α ρ ) − (uα ∂1β p + uβ ∂1α p) − uα uβ S˜ (1)

(A.18)

] −

∑ i

(1)

ciα ciβ Gi

+

∑

ciβ Gi .

i

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

20

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

where the kinetic viscosity is determined by

ν = cs2 (τg − 0.5)∆t .

(A.19)

Compared to Eq. (5b), we can recover the momentum equation as long as the following equations hold,

∑

ciβ Gi = Fβ ,

(A.20)

i

∑

ciα ciβ Gi =∂t pδαβ + ∂γ (cs2 ρ uγ δαβ ) + uα Fβ + uβ Fα + cs2 (uα ∂β ρ + uβ ∂α ρ )−

i

(uα ∂β p + uβ ∂α p) − uα uβ S˜ +

(

2 d

ρ

cs2

(A.21)

) ξ − ∂γ uγ δαβ . ∆t(τg − 0.5)

Based on Eqs. (A.20), (A.21), and the last equation in Eq. (A.13), the force distribution function can be given by Gi = ωi

⎧ ⎨

[

S+

⎩

ci · F cs2

+

]⎫

˜ − uuS˜ + Q (∇ · uI) ⎬ (ci ci − cs2 I) : ∂t (p − ρ˜ cs2 )I + uF˜ + Fu 2cs4

,

(A.22)

⎭

where F˜ is a modified total force F˜ = F − ∇ (p − ρ cs2 ),

(A.23)

S˜ and Q can be expressed as S˜ = S + ∂t (ρ − ρ˜ ), Q =

2 d

ρ cs2 −

(A.24)

ξ . ∆t(τg − 0.5)

(A.25)

Appendix B. The computation of the pressure Now we will focus on how to calculate the pressure from the distribution function gi . According to the expression of eq g0 , we have p cs2

eq

(1 − ω0 ) = ρ˜ + ρ s0 (u) − g0 .

(B.1)

eq

Actually, once g0 in Eq. (B.1) is replaced by the distribution function gi , one can present a scheme to calculate pressure. Firstly, from Eq. (A.4b) we can derive

[ ] 1 (1) )Gi + O(∆t 2 ), ϵ gi(1) = −ϵτg ∆t D1i gi(0) − (1 − 2τg

(B.2)

or equivalently,

[

eq

(0)

gi − gi = −ϵτg ∆t D1i gi

− (1 −

1 2τg

(1)

]

)Gi

+ O(∆t 2 ).

(B.3)

Taking the zeroth-direction of Eq. (B.3), we have

[

eq

(0)

g0 − g0 = −ϵτg ∆t ∂t1 g0 − (1 −

1 2τg

(1)

)G0

]

+ O(∆t 2 )

[ ] 1 (0) = −τg ∆t ∂t g0 − (1 − )G0 + O(∆t 2 ). 2τg (0)

Note that the term ∂t g0

(B.4)

is the order of O(Ma2 ), thus, we get

−g0eq = −g0 + τg ∆t(1 −

1 2τg

)G0 + O(∆t 2 + ∆tMa2 ).

(B.5)

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

21

Neglecting the terms of O(∆t 2 + ∆tMa2 ), and substituting Eqs. (B.5) into (B.1), we can obtain p cs2

(1 − ω0 ) = ρ˜ + ρ s0 (u) − g0 + ∆t(τ −

= ρ˜ + ρ s0 (u) − (

∑

gi −

i

= ρ˜ + ρ s0 (u) − (

∑

∑

= ρ˜ + ρ s0 (u) − ρ˜ + ∆t 2

S+

2

)G0

gi ) + ∆t(τ −

i̸ =0 eq

gi −

∆t ∑

i

= ρ s0 (u) +

1

∆t 2

∑

S+

2

∑

Gi −

i

1 2

∑

gi ) + ∆t(τ −

i̸ =0

gi + ∆t(τ −

i ̸ =0

gi + ∆t(τ −

i ̸ =0

)G0

1 2

1 2

1 2

)G0

(B.6)

)G0

)G0 .

As a result, the pressure can be calculated as p=

cs2 1 − ω0

⎡ ∑ ⎣ i̸ =0

gi +

∆t 2

⎤ S + ρ s0 (u) + ∆t(τ −

1 2

)G0 ⎦ ,

(B.7)

which has an accuracy of O(∆t 2 + ∆tMa2 ). References [1] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the computations of multiphase flow, J. Comput. Phys. 169 (2) (2001) 708–759. [2] X. He, S. Chen, R. Zhang, A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh–Taylor instability, J. Comput. Phys. 152 (2) (1999) 642–663. [3] T. Ramstad, P.E. Øren, S. Bakke, Simulation of two-phase flow in reservoir rocks using a lattice Boltzmann method, SPE J. 15 (04) (2010) 917–927. [4] O. Aursjø, S.R. Pride, Lattice Boltzmann method for diffusion-limited partial dissolution of fluids, Phys. Rev. E 92 (1) (2015) 013306. [5] F. Jiang, T. Tsuji, Estimation of three-phase relative permeability by simulating fluid dynamics directly on rock-microstructure images, Water Resour. Res. 53 (1) (2017) 11–32. [6] H. Liang, B.C. Shi, Z.H. Chai, Lattice Boltzmann modeling of three-phase incompressible flows, Phys. Rev. E 93 (1) (2016) 013308. [7] G. Falcucci, S. Ubertini, D. Chiappini, S. Succi, Modern lattice Boltzmann methods for multiphase microflows, IMA J. Appl. Math. 76 (2011) 712–725. [8] G. Falcucci, G. Bella, G. Chiatti, S. Chibbaro, M. Sbragaglia, S. Succi, Lattice Boltzmann models with mid-range interactions, Commun. Comput. Phys. 2 (6) (2007) 1071–1084. [9] K. Vafai, C.L. Tien, Boundary and inertia effects on flow and heat transfer in porous media, Int. J. Heat Mass Transfer 24 (2) (1981) 195–203. [10] Y. Peng, C. Shu, Y.T. Chew, Simplified thermal lattice Boltzmann model for incompressible thermal flows, Phys. Rev. E 68 (2) (2003) 026701. [11] C. Pan, M. Hilpert, C. Miller, Lattice-Boltzmann simulation of two-phase flow in porous media, Water Resour. Res. 40 (1) (2004) W01501. [12] Z.H. Chai, H. Liang, R. Du, B.C. Shi, A lattice Boltzmann model for two-phase flow in porous media, SIAM J. Sci. Comput. 41 (4) (2019) B746–B772. [13] G.D. Ilio, B. Dorschner, G. Bella, S. Succi, I.V. Karlin, Simulation of turbulent flows with the entropic multirelaxation time lattice Boltzmann method on body-fitted meshes, J. Fluid Mech. 849 (2018) 35–56. [14] B. Dorschner, S.S. Chikatamarla, I.V. Karlin, Transitional flows with the entropic lattice Boltzmann method, J. Fluid Mech. 824 (2017) 388–412. [15] G.D. Ilio, D. Chiappini, S. Ubertini, G. Bella, S. Succi, Fluid flow around NACA 0012 airfoil at low-Reynolds numbers with hybrid lattice Boltzmann method, Comput. Fluids 166 (2018) 200C208. [16] A. Zarghami, G. Falcucci, E. Jannelli, S. Succi, M. Porfiri, S. Ubertini, Lattice Boltzmann modeling of water entry problems, Internat. J. Modern Phys. C 25 (2014) 1441012. [17] H. Abels, H. Garcke, G. Grün, Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities, Math. Models Methods Appl. Sci. 22 (2012) 1150013. [18] S. Dong, An efficient algorithm for incompressible N-phase flows, J. Comput. Phys. 276 (2014) 691–728. [19] H. Liang, Y. Li, J.X. Chen, J.R. Xu, Axisymmetric lattice Boltzmann model for multiphase flows with large density ratio, Int. J. Heat Mass Transfer 130 (2019) 1189–1205. [20] I. Halliday, L.A. Hammond, C.M. Care, K. Good, A. Stevens, Lattice Boltzmann equation hydrodynamics, Phys. Rev. E 64 (2001) 011208. [21] O. Aursjø, E. Jettestuen, J.L. Vinningland, A. Hiorth, On the inclusion of mass source terms in a single-relaxation-time lattice Boltzmann method, Phys. Fluids 30 (5) (2018) 057104. [22] T. Krüger, H. Kusumaatmaja, G. Silva, O. Shardt, A. Kuzmin, E.M. Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer International Publishing, Switzerland, 2017. [23] Z.H. Chai, B.C. Shi, Z.L. Guo, A multiple-relaxation-time lattice Boltzmann model for general nonlinear anisotropic convection–diffusion equations, J. Sci. Comput. 69 (1) (2016) 355–390. [24] H.L. Wang, Z.H. Chai, B.C. Shi, H. Liang, Comparative study of the lattice Boltzmann models for Allen-Cahn and Cahn-Hilliard equations, Phys. Rev. E 94 (2016) 033304. [25] X. He, G. Doolen, Lattice Boltzmann method on curvilinear coordinates system: flow around a circular cylinder, J. Comput. Phys. 134 (2) (1997) 306–315. [26] Z.L. Guo, B.C. Shi, N.C. Wang, Lattice BGK model for incompressible Navier–Stokes equation, J. Comput. Phys. 165 (1) (2000) 288–306. [27] C.K. Aidun, J.R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid. Mech. 42 (2010) 439–472. [28] Y. Shi, G.H. Tang, Investigation of coalesced droplet vertical jumping and horizontal moving on textured surface using the lattice Boltzmann method, Comput. Math. Appl. 75 (2018) 1213–1225.

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.

22

X. Yuan, Z. Chai, H. Wang et al. / Computers and Mathematics with Applications xxx (xxxx) xxx

[29] L.-P. Wang, H. Min, C. Peng, N. Geneva, Z.L. Guo, A lattice-Boltzmann scheme of the Navier–Stokes equation on a three-dimensional cuboid lattice, Comput. Math. Appl. (2016) http://dx.doi.org/10.1016/j.camwa.2016.06.017. [30] I. Halliday, L. Hammond, C. Care, K. Good, A. Stevens, Lattice Boltzmann equation hydrodynamics, Phys. Rev. E 64 (1) (2001) 011208. [31] Y. Cheng, J. Li, Introducing unsteady non-uniform source terms into the lattice Boltzmann model, Internat. J. Numer. Methods Fluids 56 (6) (2008) 629–641. [32] P.C. Hohenberg, B.I. Halperin, Theory of dynamic critical phenomena, Rev. Mod. Phys. 49 (1977) 435–479. [33] J. Shen, X. Yang, Q. Wang, Mass and volume conservation in phase field models for binary fluids, Commun. Comput. Phys. 13 (4) (2013) 1045–1065. [34] Y.Q. Zu, S. He, Phase-field-based lattice Boltzmann model for incompressible binary fluid systems with density and viscosity contrasts, Phys. Rev. E 87 (4) (2013) 043301. [35] H.W. Zheng, C. Shu, Y.T. Chew, Lattice Boltzmann interface capturing method for incompressible flows, Phys. Rev. E 72 (5) (2005) 056705. [36] H. Liang, B.C. Shi, Z.L. Guo, Z.H. Chai, Phase-field-based multiple-relaxation-time lattice Boltzmann model for incompressible multiphase flows, Phys. Rev. E 89 (5) (2014) 053320. [37] K. Yang, Z.L. Guo, Lattice Boltzmann method for binary fluids based on mass-conserving quasi-incompressible phase-field theory, Phys. Rev. E 93 (4) (2016) 043303. [38] C.H. Zhang, K. Yang, Z.L. Guo, A discrete unified gas-kinetic scheme for immiscible two-phase flows, Int. J. Heat Mass Transfer 126 (2018) 1326–1336. [39] S.O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, J. Comput. Phys. 100 (1) (1992) 25–37. [40] L.S. Luo, Unified theory of lattice Boltzmann models for nonideal gases, Phys. Rev. Lett. 81 (8) (1998) 1618. [41] R. Du, B.C. Shi, A novel scheme for force term in the lattice BGK model, Internat. J. Modern Phys. C 17 (07) (2006) 945–958. [42] X. He, Q. Zou, L.-S. Luo, Analytic solutions of simple flows and analysis of nonslip boundary conditions for the lattice Boltzmann BGK model, J. Stat. Phys. 87 (1–2) (1997) 115–136. [43] T. Lee, C.L. Lin, A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio, J. Comput. Phys. 206 (1) (2005) 16–47. [44] T. Lee, L. Liu, Lattice Boltzmann simulations of micron–scale drop impact on dry surfaces, J. Comput. Phys. 229 (20) (2010) 8045–8063. [45] F. Ren, B.W. Song, M.C. Sukop, H.B. Hu, Improved lattice Boltzmann modeling of binary flow based on the conservative Allen-Cahn equation, Phys. Rev. E 94 (2) (2016) 023311. [46] 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, J. Comput. Phys. 334 (2017) 620–638. [47] G. Emanuel, Bulk viscosity in the Navier–Stokes equations, Internat. J. Engrg. Sci. 36 (11) (1998) 1313–1323. [48] D. Jacqmin, Calculation of two-phase Navier–Stokes flows using phase-field modeling, J. Comput. Phys. 155 (1) (1999) 96–127. [49] V.M. Kendon, M.E. Cates, I. Pagonabarraga, J.C. Desplat, P. Bladon, Inertial effects in three-dimensional spinodal decomposition of a symmetric binary fluid mixture: a lattice Boltzmann study, J. Fluid Mech. 440 (2001) 147–203. [50] V.E. Bandalassi, H.D. Ceniceros, S. Banerjee, Computation of multiphase systems with phase field models, J. Comput. Phys. 190 (2) (2003) 371–397. [51] D. Jacqmin, Contact-line dynamics of a diffuse fluid interface, J. Fluid Mech. 402 (2000) 57–88. [52] Y.Y. Yan, Y.Q. Zu, A lattice Boltzmann method for incompressible two-phase flows on partial wetting surface with large density ratio, J. Comput. Phys. 227 (1) (2007) 763–775. [53] Z. Chai, D. Sun, H. Wang, B. Shi, A comparative study of local and nonlocal Allen-Cahn equations with mass conservation, Int. J. Heat Mass Transfer 122 (2018) 631–642. [54] H.W. Zheng, C. Shu, Y.T. Chew, A lattice Boltzmann model for multiphase flows with large density ratio, J. Comput. Phys. 218 (1) (2006) 353–371. [55] A. Fakhari, M.H. Rahimian, Phase-field modeling by the method of lattice Boltzmann equations, Phys. Rev. E 81 (3) (2010) 036707. [56] J.J. Huang, C. Shu, Y.T. Chew, Mobility-dependent bifurcations in capillarity-driven two-phase fluid systems by using a lattice Boltzmann phase-field model, Internat. J. Numer. Methods Fluids 60 (2009) 203–225. [57] Z.H. Chai, T.S. Zhao, Lattice Boltzmann model for the convection–diffusion equation, Phys. Rev. E 87 (6) (2013) 063309. [58] J.W. Cahn, Phase separation by spinodal decomposition in isotropic systems, J. Chem. Phys. 42 (1) (1965) 93–99. [59] Y. Ba, H. Liu, Q. Li, Q. Kang, J. Sun, Multiple-relaxation-time color-gradient lattice Boltzmann model for simulating two-phase flows with high density ratio, Phys. Rev. E 94 (2) (2016) 023310. [60] L. Zheng, S. Zheng, Q. Zhai, Lattice Boltzmann equation method for the Cahn-Hilliard equation, Phys. Rev. E 91 (1) (2015) 013309. [61] Y. Wang, C. Shu, H.B. Huang, C.T. Teo, Multiphase lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio, J. Comput. Phys. 280 (2015) 404–423. [62] Y.B. Gan, A.G. Xu, G.C. Zhang, S. Succi, Discrete Boltzmann modeling of multiphase flows: hydrodynamic and thermodynamic non-equilibrium effects, Soft Matter 11 (26) (2015) 5336–5345.

Please cite this article as: X. Yuan, Z. Chai, H. Wang et al., A generalized lattice Boltzmann model for fluid flow system and its application in two-phase flows, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.10.007.