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

###### Abstract

In this paper, a generalized lattice Boltzmann (LB) model with a mass source 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 mass source term. 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.

###### keywords:

generalized lattice Boltzmann model, mass source term, incompressible and nearly incompressible N-S equations, fluid flow system, two-phase flow^{†}

^{†}journal: Computers & Mathematics with Applications

## 1 Introduction

Fluid flow problems are ubiquitous in nature and engineering applications, such as single and multiphase flows Tryggvason (); He0 (); Aursjo0 (); Ramstad (); Jiang (); Liang1 (), thermal flows Vafai (); Peng (), and porous media flows Pan (); Jiang (). For such problems, the phenomenon of fluid flowing in and out of the system is very frequent. In order to describe this phenomenon, an alternative way is to introduce a mass source or sink term in the fluid dynamical equations (i.e., N-S equations) Aursjo (); Kuzmin (); Crowe (); Migdal (). Besides, in the presence of electro-chemical reactions or multiphase situations, mass sources in the N-S equations could be applied to describe the reaction between the components in the system Aursjo (); Dutta (). 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 mass source term.

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 Kruger (); Chai0 (). This method has gained great success in the study of fluid flow system He1 (); Guo1 (); Aidun (); Liang2019 (). 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 , and the pressure can be obtained from the equation of state (). While for the incompressible model, the macroscopic quantities we need to calculate are fluid velocity and the pressure , 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 mass source term. For instance, Halliday et al. Halliday () proposed a single-relaxation-time (SRT) LB model including a mass source term, while they employed a non-local scheme to calculate the spatial derivatives which appear in the source term. Cheng et al. Cheng1 () presented another LB model with a general mass source term and adopted a non-local scheme for temporal and spatial derivatives. Aursjø et al. Aursjo () also developed an SRT model with a mass 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. While, up to now, there is no available work on incompressible N-S equations with a mass source term. 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 mass source term, 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 mass source term, and we present a modified scheme to calculate the pressure , 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. The present model 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 mass source term. Actually, there are also some phase-field-based LB models for incompressible two-phase flows. He et al. He0 () 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 Zu (). To recover the Cahn-Hilliard (CH) equation correctly, Zheng et al. Zheng1 () and Zu et al. Zu () 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 equals to Zu (). To overcome these problems, Liang et al. Liang () proposed a new multi-relaxation-time (MRT) LB model through introducing a time-dependent source term in the evolution equation of phase field. While all the above models cannot conserve mass locally when the densities of the two fluids are different. To solve the problem, Yang et al. Yang () 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, ( and are characteristic time and length, respectively). Based on the work of Yang et al. Yang (), Zhang et al. Zhangchunhua () proposed a discrete unified gas-kinetic scheme (DUGKS) for two-phase flows which can exactly guarantee the mass conservation, while this model also give 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 Sec. II, the generalized LB model for fluid flow system with a mass source is introduced, and a phase-field-based LB model for two-phase flows is given in Sec. III. Numerical experiments to validate the present model are carried out in Sec. IV. Finally, some conclusions are given in Sec. V.

## 2 Generalized LB model for fluid flow system with a mass source

### 2.1 Governing equations

The governing equations for nearly incompressible fluid flows with a mass source term can be written as Aursjo (); Unverdi (); Liang2018 ()

(1a) | |||

(1b) |

where is the density, u denotes the fluid velocity, is a mass source or sink term, p is the hydrodynamic pressure, F is the external force, is the deviatoric stress tensor, and for Newtonian fluids,

(2) |

where denotes the dynamic viscosity by , is the kinematic viscosity, is the bulk (or volume) viscosity and is the number of spatial dimensions.

If we take in Eq. (1), the governing equations for nearly incompressible fluid flow system can be obtained Kruger (),

(3a) | |||

(3b) |

For fluid flows with small temperature changes, the flow can be regarded as incompressible under the condition of , so that the above governing equations will reduce to Kruger (); Guo1 ()

(4a) | |||

(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 mass source,

(5a) | |||

(5b) |

where is 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 mass source term. When and take different values, Eqs. (5) represents different governing equations. For example, if we take , Eqs. (5) will reduce to Eqs. (1), and further reduce to Eqs. (3) with . While if we take , and , 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 mass 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 mass source term

To obtain the evolution equation of Eq. (5), we integrate the following discrete velocity Boltzmann equation

(6) |

along a characteristic line over a time interval Luo (); Du (), and we have

(7) |

where denotes particle distribution function with velocity at position and time , represents the force term, is the collision operator approximated by

(8) |

where is the relaxation time and is the equilibrium distribution function.

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

(9) |

Let , we have

(10) |

where is the dimensionless relaxation time, and satisfies , and .

Through Taylor expansion of and neglecting the terms of order , the last term in the right hand of Eq. (10) becomes

(11) |

where .

The LB evolution equation with the BGK collision operator for the N-S equations can be expressed as Du ()

(12) |

If the up-wind scheme is used to Eq. (12), the evolution equation for the generalized N-S equations reads

(13) |

To remove the implicitness, we introduce a modified particle distribution function,

(14) |

With some simple manipulations, the explicit evolution equation can be derived,

(15) |

where is the new equilibrium distribution function and satisfies .

To derive Eq. (5) through Chapman-Enskog analysis, the equilibrium distribution function is defined as (see A for the details)

(16) |

with

(17) |

where denotes the weighting coefficient, is the discrete velocity, and is the speed of sound. Note that our model is based on the DQ lattice with velocity directions in -dimensional space. and depend on the choice of the lattice model, and in DQ model, , , , , where , with and representing the spacing and time step, respectively; while in the DQ model, is given by , , , , and is given by

in the DQ model, , , , , and is given by

Different from some previous LB models He2 (); Yang (); Liang (); Lee1 (); Lee2 (); Zu (); Ren (); Fakhari2017 (), in the present model, the force distribution function is given by

(18) |

where is a modified total force

(19) |

and can be expressed as

(20) |

(21) |

Under the Stokes hypothesis, the bulk viscosity is usually assumed to be zero George (). In addition, can also take , so that the derived macroscopic equations will not contain the term . Please note that is a complete form and doesn’t contain any unnecessary assumptions, and 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 A for the details), and simultaneously, the fluid kinematic viscosity can be determined by

(22) |

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

(23) |

where Eq. (14) is used. By taking the first-order moment of , the fluid velocity can be obtained Liang (); Liang1 (); Liang2018 (),

(24) |

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

(25) |

or

(26) |

We would also like to point out that the present pressure expression is different from that in Refs. Liang (); Liang2018 (); Ren () which can be expressed as

(27) |

where the term of has been neglected. However, the last term may be significant since is usually nonzero under the condition of , and the effect of this item cannot be ignored.

Thus, our generalized LB model for fluid flow system are 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 (), which is necessary for the construction of most LB models. Under this assumption, the following equations are established for nearly incompressible or incompressible flows, , , and . Further, we have and . 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 mass source term.

Remark 1. When taking and , the generalized model can reduce to the nearly incompressible form. Correspondingly, the macroscopic equations becomes Eqs. (1), while the equilibrium distribution function can be written as

(28) |

and the force distribution function is given by

(29) |

where the term is the order of and can be neglected. This model is almost identical to the one in Ref. Aursjo (), 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,

(30) |

Multiplying by on both sides of Eq. (30) and ignoring the term of , one can get

(31) |

where we still use to represent the new equilibrium distribution function. The force distribution function is given by

(32) |

where , , , and the terms of are abandoned in the incompressible limit. The fluid velocity and pressure can be obtained by and . The corresponding incompressible N-S equations with a source term can be obtained from Eqs. (5),

(33a) | |||

(33b) |

where . In the limit of a low Mach number (), the dynamic pressure is assumed to be and the left end term of Eq. (33a) can be ignored. Note that if we take , this model can reduce to the incompressible LB model by He et al. He1 ().

Remark 3. When taking , (e.g., ), where is the average velocity of the fluid, we can derive a new model for incompressible N-S equations with a mass source term. In this case, the macroscopic equations can be written as [see Eqs. (5)]

(34a) | |||

(34b) |

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

(35) |

and the force distribution function is given by

(36) |

where the term and are neglect in the incompressible limit.

We note that if we take , this model can reduce to the incompressible LB model by Guo et al. Guo1 () with a force term, i.e., Eqs. (4).

Remark 4. Eq. (24) may be implicit if the force is a nonlinear function of . To remove this implicitness, we can discretize by Du ()

(37) |

Then the evolution equation can be rewritten as

(38) |

This scheme is explicit and can be iterated if is known, while the results of at time must be saved additionally.

## 3 Phase-field-based LB model for two-phase flows

In Sec. II, we have given the generalized LB model for single-phase flows with a mass 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 Cahn-Hilliard (C-H) equation is usually adopted to capture the phase interface which can be given by Jacqmin (); Kendon (); Badalassi ()

(39) |

where 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,

(40) |

In this work, denotes the phase B while represents the phase A, and is the chemical potential, which is derived by the variation of the free-energy function with respect to the order parameter Jacqmin (); Jacqmin1 (); Yan (),

(41) |

where can be taken as the following form Jacqmin (); Kendon (); Badalassi (); Chai3 (),

(42) |

where is the physical domain of the system, denotes the bulk free-energy density, and accounts for the surface energy with a positive coefficient . If the system considered is a van der Waals fluid, the bulk free energy has a double-well form Jacqmin (),

(43) |

where is a constant dependent on the interfacial thickness and the surface tension Jacqmin (),

(44) |

and

(45) |

The equilibrium interface profile can be obtained by minimizing with respect to the function , i.e., . In a plane interface at the equilibrium condition, the order parameter profile across the interface (along the direction) is represented by Yan ()

(46) |

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

(47a) | |||

(47b) | |||

(47c) |

where represents the total force which is defined as , and is the body force, is the surface tension with the potential form Jacqmin (); Zheng () if not specified. Eq. (47b) can be derived from Eq. (5a), if , and .

### 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 [Eq. (47b) and Eq. (47c)] when substituting , and into Eq. (16) and Eq. (18),

(48) |

(49) |

where we take the bulk viscosity equal to , and the dynamic pressure satisfies so that the term of is neglected in the limit of a low Mach number. Actually, the relationships of and are also true. Specially, the force distribution function can be simplified as

(50) |

From the above equation, one can find that is related to since the term is the order of , and is nonzero once , thus the term in Eq. (25) cannot be ignored in the pressure expression.

We would like to point out that Eq. (47) can reduce to the incompressible phase-field model if . However, based on Eqs. (47a), (47b), (40) and , one can obtain Yang ()

(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 . Therefore, in the incompressible phase-field model, the mass is not locally conserved.

To conserve mass locally, Shen et al. Shen () proposed a quasi-incompressible phase-field model. Subsequently, based on the quasi-incompressible phase-field model, Yang et al. Yang () designed a lattice Boltzmann model for binary fluids. Actually, Eq. (47) can also reduce to quasi-incompressible phase-field model in Ref. Shen () when with . It should also be noted that the momentum equation (47c) is different from those used in Refs. Liang (); Yang (), where the terms of and are not included.

Remark 1. Here we also give a compasion between the generalized LB model and the one in Ref. Yang (). If we take , and , the equilibrium distribution function for the N-S equations can be written as

(52) |

Based on the following fact,

(53) |

the pressure can be calculated as

(54) |

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 takes the same form.

#### 3.2.2 LB model for the C-H equation

The evolution equation for the C-H equation can be given as

(55) |

where is the distribution function of order parameter , is the non-dimensional relaxation time related to the mobility, is the local equilibrium distribution function, which is defined as Fakhari2010 (); Huang (); Chai ()

(56) |

where is an adjustable parameter that controls the mobility. is the source term and is given by Liang ()

(57) |

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

(58) |

Through the Chapman-Enskog analysis, the order parameter is calculated by

(59) |

the C-H equation can be recovered with second-order accuracy and the mobility can be determined by

(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 a inverse linear form Lee2 (); Fakhari2017 () is adopted to ensure a smooth viscosity across the interface if not specified,

(61) |

In addition, to determine the gradient terms in the source term , surface tension and chemical potential , the following isotropic schemes are adopted to discretize the first-order spatial derivative and the Laplacian operator Liang ():

(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 Liang1 ().

## 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 () and quasi-incompressible [] two-phase flows. In our simulations, the DQ lattice structure is adopted. We attempt to conduct a detailed comparison between the present and the analytical solutions or some available results in each test.

### 4.1 Spinodal decomposition

Spinodal decomposition Cahn1 () 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., Eq. (25) and Eq. (27), and this example is mainly used to demonstrate the differences in calculation results obtained by present incompressible and quasi-incompressible models. The computational domain is , the initial distribution of order parameter with a small fluctuation is given by

(63) |

where is a random function with the maximum amplitude of . 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 . 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-compressible models are significantly different, which may be caused by the fact that the term has been neglected in incompressible model, and thus mass conservation is not satisfied locally.

To illustrate the difference between two pressure expressions, we presented some results in Fig. 2. As seen in Fig. 2(a), the distributions of order parameter within the red circles are distinctly different, which means that the term 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 can be ignored when , and is the order of . Based on the above results, one can conclude that the pressure expression Eq. (25) is more general, and would be adopted in the following simulations.

### 4.2 Static droplet

A static droplet is popular problem to verify LB models for two-phase flow Zu (); Ba (); Liang (); Lee1 (); Fakhari2010 (). In this subsection, we will consider this problem with different density ratios. Initially, a circular droplet with the radius ranging from to is placed in the middle of the computational domain with . The initial order parameter is given by

(64) |

where is the droplet radius, and surface tension is expressed as . In numerical simulations, we set the density ratio to be , and . The other physical parameters are fixed as , and the periodic boundary conditions are applied at all boundaries. We first verify the LB model with the well-known Laplace’s law,

(65) |

where is the pressure jump across the interface, is the radius of the droplet. is calculated by with the equation of state Zu (); ZhengL (). We performed some simulations and presented the results in Fig. 3. From Fig. 3(a), 3(c), 3(e) one can find that the results of present models and those in Refs. Liang (); Yang () 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 in Table 1. In general, the present quasi-incompressible LB model and quasi-incompressible model in Ref. Yang () are more accurate than incompressible LB models. This is because that quasi-incompressible model is physically more reasonable than incompressible model. Besides, it is also found that usually present QIM is more accurate than the model in Ref. Yang ().

In addition, we also presented the density profiles along the horizontal center line at in Figs. 3(b), 3(d), and 3(f). From these figures, one can observe that all the numerical results are very close to the analytical solutions given by Eq. (64). To give a quantitatively estimation on the accuracy of numerical results, we also measured the relative errors of density, i.e., [] with in Fig. 4(a), where 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. Liang () produce smaller errors than quasi-incompressible models, which means that the quasi-incompressible and incompressible models each have their own advantages.

Present QIM | Present IM | Yang et al. | Liang et al. | |
---|---|---|---|---|

20 | ||||

25 | ||||

30 | ||||

35 | ||||

40 |