Discrete methods of the energy equations in the pseudo-potential lattice Boltzmann model based simulations

Computers and Fluids 179 (2019) 645–654

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Discrete methods of the energy equations in the pseudo-potential lattice Boltzmann model based simulations Anjie Hu a,∗, Rizwan Uddin b, Dong Liu a a b

School of Civil Engineering and Architecture, Southwest University of Science and Technology, Mianyang 621010, PR China Department of Nuclear, Plasma and Radiological Engineering, University of Illinois at Urbana Champaign, Urbana, IL 61801, USA

a r t i c l e

i n f o

Article history: Received 29 May 2018 Revised 21 August 2018 Accepted 5 December 2018 Available online 6 December 2018 Keywords: Lattice Boltzmann method Pseudo-potential model Energy equations

a b s t r a c t Pseudo-potential lattice Boltzmann (LB) Methods have been widely applied in the gas-liquid multiphase flow and heat transfer simulations such as the bubble growth in the boiling process. Various formats of energy equations and discrete models are deduced to solve the heat transition in the simulations. However, the accuracy of these models has rarely been discussed. In this paper, four common LB thermal models are analyzed and numerically compared with the finite difference method by simulating the phenomena of single-phase natural convection and the one-dimensional heat transfer across the vaporliquid interface. Simulation results show that all these LB thermal models lead to unwanted errors compared with the benchmark simulation results, while, with the proper discrete scheme, the finite difference method provides results with much better accuracy. With the presented finite difference model, we further compared the energy conservation of the energy equations in the phase change simulation. It is found that the accuracy and the thermodynamic consistency of the temperature-based energy equation is better than the internal-energy-based energy equation. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Heat transfer within the gas–liquid multiphase flow widely exists in nature and many engineering heat exchangers. Phenomena such as boiling and condensation play important roles in many engineering applications due to their large latent heat capacity and high–heat transfer efficiency. These phenomena are often complex hence it is difficult to study by the normal experimental methods precisely. With the development of the computer science in the recent decades, the numerical methods play a more and more important role in the study of multiphase flow and heat transfer research. One of these methods is lattice Boltzmann equation (LBE) method [1] also known as lattice Boltzmann method (LBM). Due to its mesoscopic background, simplicity and strict second-order accuracy, it has attracted much attention in recent years. Soon after its appearance, several multiphase LBE models have been developed. These models can be summarized into four categories [2]: color models [3], pseudo–potential models [4–6], free energy models [7, 8] and kinetic models [9, 10]. Among them, the pseudopotential models are the most widely applied, and with the efforts

∗

Corresponding author. E-mail address: [email protected] (A. Hu).

https://doi.org/10.1016/j.compfluid.2018.12.005 0045-7930/© 2018 Elsevier Ltd. All rights reserved.

have been made in the recent years [11-16], these models have been highly improved on the aspects of stability and accuracy. In the last two decades, efforts have been taken to apply the LB methods to vapor–liquid phase change problems such as boiling and evaporation. Several models are established by introducing thermal equations into LB multiphase models [17-20], and the most widely applied phase change models are pseudo-potential models. Among these models, two energy equations are applied to solve the heat transfer: one is the internal-energy-based equation which was adopted in the first pseudo-potential phase change model proposed by Zhang and Chen [21]; The other one is the temperature-based equation proposed by Haze and Markus [19], this equation has been widely applied in the later researches. Besides of the energy equations, the discrete models are also proposed in these simulations. Among them, the lattice Boltzmann thermal models [22-26] are most popular. However, Li et al. [27,28] analyzed these LB thermal models and pointed out that errors exist in these models, modifications are necessary when applied them in the pseudo-potential modeling of heat transfer simulations. To avoid the errors produced by the LB thermal models, the finite difference (FD) scheme [29-32] is adopted to solve the energy equations in their simulations of phase change problems. Since the discrete models and the energy equations have great influence on the pseudo-potential phase change model, to evaluate the accuracy of these methods, in the present work, four

646

A. Hu, R. Uddin and D. Liu / Computers and Fluids 179 (2019) 645–654

most common used LB thermal models are analyzed with the Chapman-Enskog expansion and numerically tested and compared with the FD scheme by simulating the single-phase natural convection problem and the heat transfer within the multiphase flow. The model that performed best in the simulations was then applied in the energy consistence verification of the pseudo-potential phase change model. The rest of this paper is arranged as follows, Section 2 briefly introduces the basic theory of pseudo-potential LB method; these LB thermal models are introduced and theoretically analyzed in Section 3; then the numerical verification and comparison are made in Section 4; the energy conservation verification is given in Section 5; finally, a brief conclusion is made in Section 6.

2. Pseudo-potential model In this work, the MRT pseudo-potential model is adopted due to its advantages of stability and easier tunable parameters. In the LBE method, the motion of the fluid is described by the density distribution function evolution which can be written in the form of MRT operator [33,34] as

˜ i j ( f j (x, t ) − f eq (x, t )) fi (x + ei δt , t + δt ) = fi (x, t ) −  j





˜ i j S j, +δt Ii j − 0.5

(1)

where fi (x, t) is the mass distribution function of particles at node x, time t; ei is the discrete velocity where i = 0, 1, 2 · · · N; eq fi (x, t ) is the equilibrium distribution. The right side of the equa˜ i j is the collision matrix equals to tion is a collision operator and  M−1 M, in which M is the orthogonal transformation matrix and  is a diagonal matrix which is consisted of a set of relaxation times (D2Q9):

 = diag(τρ−1 , τe−1 , τξ−1 , τ j−1 , τq−1 , τ j−1 , τq−1 , τυ−1 , τυ−1 , ).

(2)

where τυ−1 is related to the viscosity of the fluid υ , the relationship between the relax time and the viscosity can be written as υ = 1 1 3 (τυ − 2 )δ t, for simplicity, other the relax times are all set to 1 in the simulation. Sj is the term related to external forces. In MRT operator, the collision is calculated in the moment space. The density distribution function and its equilibrium distribution function can be transferred into the space of moments by m = Mf and meq = Mfeq . For the D2Q9 lattice, the equilibria meq are given by





meq = ρ 1, −2 + 3|u|2 , 1 − 3|u|2 , ux , −ux , uy , −uy , u2x − u2y , ux uy . (3) where u is the macroscopic velocity of the fluid. With Eqs. (2) and (3), the right-hand side of the MRT LB Eq. (1) can be rewritten as [35]



m∗ = m − (m − meq ) + I −



1  S˜ , 2

(4)

The streaming process is given by

fi (x + ei δt , t + δt ) = f ∗ (x, t ) = M−1 m∗

(5)

In the MRT LB method, the force is incorporated via the following forcing scheme [34,36]

⎡

⎤

0 ⎢6(ux Fx + uy Fy ) ⎥ ⎢−6(ux Fx + uy Fy )⎥ ⎢ ⎥ ⎢Fx ⎥ ⎢ ⎥ S = ⎢−Fx ⎥. ⎢Fy ⎥ ⎢ ⎥ ⎢−Fy ⎥ ⎣ ⎦ 2(ux Fx − uy Fy ) (ux Fx + uy Fy )

(6)

In the pseudo-potential model, phase separations are simulated by the introduction of interaction force [37,38], and its format has great influence on the stability and accuracy of the model. In this work, the hybrid model proposed by Kupershtokh et al. [15,39] in the simulation to guarantee the stability and accuracy of the two phases simulations, which is given as follows:

F(x, t ) = −A

8

ω (|eα |2 )U (x + ei )ei

i=0

− ( 1 − A ) Gψ ( x , t )

8 i=0

ω (|eα |2 )ψ (x + ei , t )ei.

(7)

where G is the interaction strength, ω(|ei |2 ) are the weights, and ψ (x, t) is the effective density. The weights ω(|eα |2 ) are given by ω (1 ) = 1/3 and ω (2 ) = 1/12, A is related to the stability of the model, which is set as −0.3 in the simulation. U(x, t) is the potential function which is equal to Gψ 2 (x, t)/2. In practice, both of effective density and potential function can be obtained by introducing a non-ideal EOS [5]:

p0 = cs2 ρ + cG[ψ (ρ )]2 /2 = cs2 ρ + cU (ρ ). where c is lattice speed which equals to δδ , t

(8)

δ is the lattice length.

For simplicity, δ and δ t are both set to 1 in the simulation. 3. Lattice Boltzmann thermal models

The energy equation of the model can be deduced with the thermodynamic theory. Neglecting the viscous dissipation, the energy conservation equation can be written as

Q ds = ρT = ∇ · (λ∇ T ), dt dt

(9)

where s is the entropy, λ is the thermal conductivity. According to the thermodynamic relations, the following equation can be obtained:

Tds = de + pdv,

Tds =

∂e ∂T

(10)



∂e dT + ∂v v

 dv + pdv,

(11)

T

where, e is the internal energy which equals to Cv T. According to Eq. (9), the energy balance equation can be obtained by expanding Eq. (10):

ρ (∂t (Cv T ) + u · ∇ (Cv T )) = ∇ · (λ∇ (T )) − pdv,

(12)

where dv = d ( ρ ) = ρ 2 (∂ (ρ ) + u · ∇ (ρ )). Considering the continuity equation of ∂t (ρ ) + ∇ (ρ u ) = 0, the above equation can be rewritten as: 1

1

ρ (∂t (Cv T ) + u · ∇ (Cv T )) = ∇ · (λ∇ (T )) − pEOS ∇ u

(13)

This equation was adopted in Zhang and Chen’s work [21], though it is well known as the first energy equation of the pseudopotential phase change model, it has been rarely applied in practice. Hazi and Markus [19] further deduced Eq. (11) with the thermodynamic relation of

∂e ∂v



∂p =T ∂ T T



− p. v

(14)

A. Hu, R. Uddin and D. Liu / Computers and Fluids 179 (2019) 645–654

Eq. (11) can then be rewritten as

Tds =

∂e ∂T



dT + T v

∂p ∂T

3.1. Chapman-Enskog expansion of LB thermal models



dv.

(15)

v

With Eq. (9), the following equation can be obtained:

∂ pEOS ρCv (∂ T + u · ∇ T ) = ∇ · (λ∇ (T )) − T ∂T



∇ u.

(16)

v

This equation is widely applied in the phase change simulations. Compared with Eq. (13), the convection term in above equation contains the gradient of T instead of the internal energy e. In essence, Eqs. (13) and (16) have no different, however, the definition of Cv in Zhang and Chen’s work does not consistent with the EOS. In this work, we will compare the performance of them after analyzing the LB thermal models. The energy equations LB thermal models are always solved by the LB thermal models. These models can all be written in the form of lattice Bhatnagar-Gross-Krook (LBGK) model as:

gi (x + ei t, t + t ) − gi (x, t ) = −(gi (x, t ) − geq (x, t ))/τg , i

(17)

eq

where gi represents thermal distribution functions and gi is the corresponding equilibrium distribution. According to the definition eq of gi , these models can be divided into two categories. One is the internal-energy-based thermal LB model, in which gi represents internal energy distribution functions and the corresponding expreseq sion of gi is given by:

geq i

= Cv T

647

fieq ,

By Chapman-Enskog expansion, Li et al. [27,28] found that errors existed when these LB thermal models are applied in the simulations with the pseud-potential model. To further investigate the accuracy of these thermal models, in this work, we first analysis these models by Chapman-Enskog expansion based on Li et al.’s work.  1 The internal-energy-based thermal LB (ITLB) model: geq = C p T f eq . The Eq. (17) can be written as follow by Tayler expansion

δt (∂t + ei · ∇ )gi +

2

(∂t + ei · ∇ )2 gi + · · · =

1

τg

(gi − geq ). i

(24)

By multi-scale expansion, the partial derivative and distribution function in the above equation can be written as

∂t = ε ∂t0 + ε 2 ∂t1 , ∂x = ε ∂x , gi = geq + ε g(i 1) + ε 2 g(i 2) . i

(25)

Eq. (24) can then be rewritten in the consecutive orders of δ t

O(ε ) : (∂t0 + ei · ∇ )geq =− i

1

δt τg

g(i 1) ,

O(ε 2 ) : ∂t1 geq + (∂t0 + ei · ∇ )g(i 1) + i

(26)

1 (2 ) 1 δt (∂t0 + ei · ∇ )2 geq =− g . i 2 δt τg i (27)

According to Eq. (26), Eq. (27) can be rewritten as

(18) eq fi

δt2



1 1 (2 ) g(i 1) = − g . 2 τg δt τg i

where T is the macroscopic temperature, is the density distribution functions of LBGK model, which is given by

O(ε 2 ) : ∂t1 geq + (∂t0 + ei · ∇ ) 1 − i

fieq (x, t ) = wi ρ (x, t )[1 + (ei · u )/cs2 + (ei · u )2 /2cs4 − u2 /(2cs2 )].

Taking the summations of distribution functions of all directions in Eqs. (26) and (28), we can get

(19) The macroscopic temperature T is related to the internal-energy distribution functions g by

 T =

i gi

Cpρ



=

geq i

i

Cpρ

.

(20)

This model is rarely been used in the pseudo-potential modeling except in the research [40] which adopts the energy equation proposed by Zhang and Chen. In 2014, Li and Luo [27] discussed the effect of forces on this LB thermal model in the framework of pseudo-potential model, hence it is also adopted in this work. The other is the temperature-based thermal LB models, in which gi represents temperature distribution functions. The corresponding equilibrium distribution is given by [18]

geq i

=

T fieq

ρ

.

(21)

In this case, the macroscopic temperatures T are related to the temperature distribution functions gi by

T =



i

gi =



geq . i i

According to Eq. (26),

i

ei g(i 1) = −τg



eq

∂t 0

 i



ei gi + ∇ ·

 i

ei ei geq i

(30)



.

(31)

eq

Noting gi = Cv T fi , the following equation can be obtained,

∂t 0





ei geq = ∂t 0 (Cv T ρ u ) = u∂t 0 (ρCv T ) + ρCv T ∂t0 u, i

i

∇·

 i

(32)



ei ei geq = ∇ · (ρCv T uu ) + ρ cs2 ∇ (Cv T ) + Cv T cs2 ∇ρ , (33) i

where

ρ∂t0 u = ∂t0 (ρ u ) − u∂t0 ρ .

(34)

Considering the macroscopic equations of LB fluid flow equation

∂t0 ρ + ∇ · (ρ u ) = 0, ∂t0 ρ u + ∇ · (ρ uu ) = −∇ cs2 ρ + F,

(35)

With Eqs. (34) and (35), Eq. (32) can be rewritten as

(23)

The above equation can be viewed as a simplified expression of eq the Eq. (21) without the terms of u2 in fi . Both of these two temperature-based LB thermal models are widely applied in the simulation of phase change problem, hence they are adopted to solve the temperature based energy equation (Eq. (16)) proposed by Hazi and Markus [19].

(29)

  1 ∂t1 (ρCv T ) + ∇ · 1 − ei g1i = 0. i 2 τg

(22)

In the early research [19], Instead of Eq. (21), the temperature equilibrium distribution is given by

geq (x, t ) = wi T (x, t )[1 + (ei · u )/cs2 ]. i

∂t0 (ρCv T ) + ∇ · (ρCv T u ) = 0,

(28)

∂t0



α

ei geq i



= −u∇ (ρCv T u ) − Cv T ∇ · (ρ uu ) − Cv T cs2 ∇ρ + Cv T F + Cv T u∇ · (ρ u ).



(36)

Substituting Eqs. (33) and (36) into Eq. (31), we have

i





ei g(i 1) = −τg cs2 ρ∇ (Cv T ) + Cv T F .

(37)

648

A. Hu, R. Uddin and D. Liu / Computers and Fluids 179 (2019) 645–654

Combining Eq. (29) with Eq. (30) and using Eq. (37), we can obtain the corresponding macroscopic equation:



 λ ∂ (ρCv T ) + ∇ · (ρCv T u )−∇ · ρα∇ (Cv T ) + 2 Cv T F = 0, (38) cs ρ

where α = ρλCv = 13 (τg − 0.5 ) is the thermal diffusivity. Note that the above equation is slightly different from literature [27], in which the thermal capacity Cv is considered as a constant and the term of ∇ (Cv T) is given by Cv ∇ (T). The thermal conductivity in the above equation is given by λ = ραCv . Tf

eq

eq  2 The temperature-based thermal LB (TTLB) model: g = ρi . i

The expansions of this model in the orders of δ t Eqs. (26)– (28) are the same as the ITLB model, taking the summations of Eqs. (26) and (28), we can obtain:

∂t0 (T ) + ∇ · (T u ) = 0,

∂t1 (T ) + ∇ · i

1 1− 2 τg

ei g(i 1) = −τg

(39)





i

ei g(i 1) = 0,

(40)

     ∂t 0 ei geq +∇ · ei ei geq . i i i i

(41)

eq

With the equilibrium distribution function gi , the terms on the right side of the equation can be written as

∂t 0





i

∇·

ei geq = ∂t 0 (T u ) = u∂t 0 (T ) + T ∂t0 u, i





i

ei geq i

∂t 0

i

ei geq i



ρ cs

∇ ( cα2 ∇ T uu) has the second order of Mach number, hence, it can s

be neglected when the flow velocities are small. However, the errors caused by ∇ ( 2λ Cv T F ) and ∇ ( α2 T ∇ PEOS ) are related to the c ρ ρc s

s

uneven forces and pressure distribution, since the interaction force and pressure are both uneven when the density changes, these errors terms cannot be neglected in the case of multiphase simulation. It also should be noted that the thermal diffusivity α is adopted in TTLB model and STTLB model instead of thermal conductivity λ, the heat diffusivity term can only be consistent with Eq. (16) when ρ Cv is constant, which cannot be satisfied for most two-phase flows. A similar problem also exists in ITLB model, which contains the term of ∇ (Cv T) instead of Cv ∇ (T), the large fluid properties changing across the interface of a two-phase flow may highly influence the accuracy of these models. To discuss the influence of the non-uniform Cv on the accuracy of the LB thermal models, error terms produced by the pseudopotential force in Eq. (38) should be eliminated, hence, the modified ITLB (MITLB) model proposed by Li and Luo is adopted in this work.





gi (x + ei t, t + t ) − gi (x, t ) = − gi (x, t ) − geq ( x , t ) / τg i + δ t i ( x , t ) ,

(43)



(49)

i =

1 1− 2 τg

 ωiCv T

( ei · F ) cs2

.

(50)

The corresponding macroscopic equation is then given by

= −u∇ (T u ) +



cs

where

ei ei geq = ∇ · (T uu ) + cs2 ∇ T . i

Considering the macroscopic momentum equations, the above equations can be written as

∂t 0

cs ρ

in Eq. (48), and ∇ ( α2 T ∇ PEOS )in Eqs. (46) and (48). The term of

(42)



i

which may cause errors in the multiphase flow heat transfer simulation. These error terms are ∇ ( 2λ Cv T F ) in Eq. (38), ∇ ( α2 ∇ T uu )

T

ρ





−∇ · (ρ uu ) − ∇ cs2 ρ + F + u∇ · (ρ u ) , (44)

 T = −∇ (T uu ) + −∇ cs2 ρ + F . ρ

∂ (ρCv T ) + ∇ · (ρCv T u ) − ∇ · [ρα∇ (Cv T )] = 0.

(51)

With this treatment, the error term of ∇ ( 2λ Cv T F ) is elimics ρ

(45)

nated in the equation, while Cv is still in the secondary partial derivative term.

Hence, the macroscopic temperature equation can be obtained by taking Eqs. (43) and (45) into Eq. (40):



 α ∂ (T ) + ∇ · (T u ) − ∇ · α∇ T − 2 T ∇ PEOS = 0, ρ cs

4. Numerical evaluation of LB thermal models

(46)

= ρλCv = 13 (τg − 0.5 ) is the thermal diffusivity in this case.

where α This equation is still different from Eq. (16) for multiphase simulations, since ρ Cv is always different for different phases.  3 The simplified temperature-based thermal LB (STTLB) model: eq gi (x, t ) = wi T (x, t )[1 + (ei · ueq )/cs2 ]. The expansion of this model is the same as TTLB model, except the Eq. (43) is given by

∇·





i

ei ei geq = cs2 ∇ T . i

(47)

The macroscopic temperature equation of this model is given as follows:



 α α ∂ (T ) + ∇ · (T u ) − ∇ · α∇ T − 2 T ∇ PEOS − 2 ∇ T uu = 0. ρ cs cs

(48) Eqs. (38), (46) and (48) are the macroscopic temperature equations for these LB thermal models, respectively. It can be seen in these equations that, compared with the standard energy equations Eqs. (13) and (16), additional terms exist in these equations,

To evaluate these LB thermal models, the problems of nature convection and one-dimensional heat transfer through the gasliquid interface are studied with the above LB thermal models. The corresponding fluid flow and density distributions are simulated with the presented pseudo-potential model. As a comparison, in this work, we also presented the simulation results by directly discretizing the energy equation with the finite difference scheme (FD scheme). The comparation between Eqs. (13) and (16) will be discussed in the next section, here we only adopted Eq. (13) in the simulation. In the FD scheme, the time step and the lattice size are given the same as the LB thermal models. The forth-order RungeKutta scheme Ref. [32] is adopted for the time discretization. These energy equations can be written as ∂t ϕ = K (ϕ ), the evolution is given by

ϕ t+δt = ϕ t +

δt 6

(h1 + h2 + h3 + h4 ),

(52)

where h1 = K (ϕ t ), h2 = K (ϕ t + δ2t h1 ),h3 = K (ϕ t + δ2t h2 ) and h4 = K (ϕ t + δt h3 ). The second order accuracy schemes proposed by Lee and Lin [41] are applied for the spatial discretization:

A. Hu, R. Uddin and D. Liu / Computers and Fluids 179 (2019) 645–654

     ∂ϕ  = ϕ(m+1, n ) − ϕ(m−1, n ) /3δ + ϕ(m+1, n+1) − ϕ(m−1, n−1) /12δ  ∂ x (m,n)   + ϕ(m+1, n−1) − ϕ(m−1, n+1) /12δ,      ∂ϕ  = ϕ(m ,n+1) − ϕ(m, n−1 ) /3δ + ϕ(m+1, n+1) − ϕ(m−1, n−1) /12δ ∂ y (m,n)   + ϕ(m−1, n+1) − ϕ(m+1, n−1) /12δ,

Table 1 Comparison of Nu numbers obtained by different models.

FD scheme MITLB model ITLB model TTLB model STTLB model Barakos et al. [35]

(53)

∂ 2ϕ ∂ 2ϕ + ∂ x2 ∂ y2

= [ϕ(m+1, n+1) + ϕ(m−1, n+1) + ϕ(m+1, n−1 ) +ϕ(m−1, n−1) + 4ϕ(m+1, n ) + 4ϕ(m−1, n )

FD scheme MITLB model ITLB model TTLB model STTLB model

+4ϕ(m, n+1) + 4ϕ(m, n−1) − 20ϕ(m, n ) ]/6δ 2 , (54) where m, n donate the lattice positions in the two coordinate directions, respectively. 4.1. Cavity natural convection The problem of natural convection in a two-dimensional square cavity is widely used to evaluate the accuracy of thermal flow simulations. In this problem, the sidewalls of the cavity are maintained at two different constant temperatures, respectively, whereas the bottom and top walls are adiabatic. This temperatures difference causes the density change and natural convection flow under the gravity force. The natural convection can be characterized by the Prandtl number and the Rayleigh number which can be written as [42]

Ra =

gβρ02 (Th − Tl )L3 P r

μ2

,

(55)

where Th and Tl are the sidewalls’ temperatures, g is the gravity acceleration, β is the thermal expansion coefficient which is given by T 2+T for ideal gas, L is the distance between the walls, Pr is h

l

the Prandtl number which is given by υ /α . To apply the pseudopotential model in the simulation, the EOS of ideal gas is adopted here, which is given by

p(ρ ) = ρ RT ,

Numax

YNumax

Nuave

3.4987 3.27697 2.3818 4.8902 4.8915 3.53

0.85333 0.84666 1 1 1 0.8531

2.1963 2.0558 1.1185 2.9028 2.9040 2.245

Table 2 Relative error of results obtained by different models.

 (m,n )

649

(56)

where R is the gas constant. The simulation domain is set as Nx × Ny = 150 × 150, and temperatures of the sidewalls are 1.03 and 0.97, respectively. The Prandtl number is fixed at Pr = 0.71. The non-equilibrium extrapolation [43] and the halfway bounce back schemes are adopted to deal with the thermal and flow boundaries, respectively. Simulation results of Ra = 10 0 0 0 are given as follows. Fig. 1 shows the temperature distributions obtained by different thermal models. It can be seen from this figure that, compared with literature [44], the temperature distribution obtained by ITLB model is obviously different from others, while the results given by other models are nearly the same. To be specific, we further compared the Nusselt numbers given by these models with the benchmark results, as shown in Tables 1 and 2, in which Numax , YNumax , Nuave donate the largest Nu number, its vertical coordinate and the corresponding average Nu number, respectively. It can be seen from these tables that, although the temperature distributions are very close to each other except ITLB model, the simulated parameters are quite different. The simulated Nu numbers given by TTLB model and STTLB model are obviously larger than the benchmark results. While simulation results given by the FD scheme and MITLB model are much closer to the benchmark

Numax

YNumax

Nuave

−1.17% −7.43% −32.72% 38.14% 38.18%

−0.43% −1.21% 16.69% 16.69% 16.69%

−2.17% −8.43% −50.18% 29.30% 29.36%

Table 3 Parameters of the EOS. Density ratio

ρ1

ρ2

θv

θm

θl

100

1.486

94.65

0.1633

−0.02

0.33333

results. Except the Nu numbers given by MITLB model are slightly deviated from the benchmark solution, which may be caused by the high order errors in the LB equation. These results show that, due to the existence of the error terms, these LB thermal models cannot get accurate results in the pseudo-potential model based simulation without modification. 4.2. Heat transition through static phase interface The errors’ influence in the above simulation also existed in the multiphase flow since the pressure gradient in the flow is much larger when the non-ideal gas EOS is applied. Besides of the error terms, the macroscopic equation of these LB thermal equations are still different from Eqs. (13) and (16). The accuracy of these models for uneven parameters must be verified in a multiphase flow. In this part, the heat transfer across a static phase interface is adopted. To eliminate the influence of the temperature on the density distribution and easily control the density ratios, the selftuning EOS proposed by Colosqui et al. [45] is adopted to simulate two-phases flow in present work. The EOS is given by



p( ρ ) =



ρθv if ρ ≤ ρ1 ρ1 θ v + ( ρ − ρ1 ) θ m if ρ1 ≤ ρ ≤ ρ2 , (57) ρ1 θv + (ρ2 − ρ1 )θm + (ρ − ρ2 )θl if ρ > ρ2 





where θv = (∂ p/∂ ρ )v and θl = (∂ p/∂ ρ )l are the speeds of sound of the vapor-phase and liquid-phase, respectively. θ m is the slope of the unstable branch (∂ p/∂ ρ < 0). The unknown variables ρ 1 and ρ 2 can be obtained after set the values of θ and equilibrium phase densities. The parameters of the EOS are given in Table 3. It should be noted that, in the simulation, the EOS is independent of the temperature T, hence the density distribution is not influenced by the heat transfer. The simulation domain is chosen as Nx × Ny = 150 × 20 as shown in Fig. 2. For the sake of contrastive analysis, two combinations of fluids’ thermal properties are considered in the simulation respectively, which are given in Table 4. The temperatures of the left side wall are set as 1.01, 1.02, 1.03, and 1.04 for each simulation, and the temperatures of the right wall are given as 0.99, 0.98, 0.97 and 0.96, respectively. The top and bottom wall are set as adi-

650

A. Hu, R. Uddin and D. Liu / Computers and Fluids 179 (2019) 645–654

(a) ITLB model

(b) TTLB model

(d) MITLB model

(c) STTLB model

(e) FD scheme

Fig. 1. Temperature distributions obtained by different models.

ten as

interface adiabatic boundary

Th

∂ (ρCv T ) + ∇ · (ρCv T u ) = λ∇ 2 T + ∇ λ∇ T

heat

liquid

vapor

adiabatic

Tl

boundary

Table 4 Thermal parameters of fluids.

(61)

Cv l = Cv v

Cv l = Cv v

λ

Cv

λ

Cv

20 1

1 1

20 1

1 10

with these treatments, the finite difference model in the present work is still in second-order accuracy.

abatic boundaries. Initially, the interface is located at the center of the domain, and the temperature within the simulation domain is given as 1. According to the energy conservation law, heat flow on the right and left wall should be equal after the heat transfer get to the equilibrium state since phase change does not exist in this case, and the analytic heat flux can be obtained with the onedimensional thermal conductivity law given as follows

q=

T , Rv + Rl

(58)

where Rv and Rl are the thermal resistance of the gas vapor and liquid phase respectively. The thermal parameters across the interface are assumed linearly changing with the fluids’ density as shown in the following equation:

X (ρ ) = Xv +

ρ −ρ  v

ρl − ρv

(Xl − Xv ),

To ensure the energy conservation of the discrete of the equation, the first order derivatives are applied on the term of ∇λ and ∇ T, which are given as follows

      ∂ϕ  ∂ϕ  = ϕ(m, j ) − ϕ(n−1, j ) /2, = ϕ(m, n ) − ϕ(m n−1 ) /2, ∂ x (m,n) ∂ y (m,n)

Fig. 2. Simulation of heat transition through static phase interface.

Liquid Vapor

(60)

(59)

where X donates the thermal factors of the fluids. One thing should be mentioned that since the thermal conductivity is not constant, in the FD discrete model Eq. (13) can be writ-

4.2.1. Simulation results for Cv l = Cv v Fig. 3 shows the simulated temperature distributions across the phase interface obtained by the presented models. It can be seen from this figure that the temperature gradients in the vapor phase are much larger than the liquid phase, this is because the heat conductivity of the liquid is much larger than the vapor. The highest temperature given by ITLB model is located at the area beside the interface, shows a strong heat source is generated in this area, which is not consistent with the physical model since phase change does not exist in the simulation. The distributions given by TTLB model and STTLB model are very close to each other but obviously different from the results given by MITLB model and FD scheme. To further evaluate the accuracy of these models, the simulated q +q average heat fluxes q¯ (given by H 2 L , where qH and qL donate the heat fluxes on the left and right boundaries) are compared with the analytical results, as shown in Fig. 4. It can be seen from this figure that the heat fluxes simulated by MITLB model and FD scheme are much closer to the analytical results, and both linearly change with the temperature difference. While the simulation results given by TTLB model and STTLB model are nearly two times as the analytical results. The heat flux differences between the left and right boundaries are given in Fig. 5 to evaluate the energy conservation of these models. It can be seen from this figure that all

A. Hu, R. Uddin and D. Liu / Computers and Fluids 179 (2019) 645–654

651

Fig. 3. Temperature distributions across the phase interface obtained by different models.

these four LB models suffer nonconservation of energy at different levels, the relative errors of MITLB model are about 7% and others are larger than 100%. However, the FD scheme in the present work is much better in the energy conservation, whose errors is no bigger than 1%. These results show that the simulations of the heat transfer across the interface are highly influenced by the error terms shown in the Chapman-Enskog expansion. By eliminating these error terms, MITLB model can give much better simulation results, however, the simulation accuracy is still worse than the finite difference model, this is perhaps because of some influences of the higher order errors due to the large fluid properties gradient across the interface.

Fig. 4. The simulated average heat flux q¯ given by different models.

Fig. 5. The difference between reduced heat fluxes from the walls obtained in simulations for different models.

4.2.2. Simulation results for Cv l = Cv v The simulation temperature distributions of the uneven heat capacity two-phase flow obtained by these models are shown in Fig. 6. Compared with Fig. 3, the temperature distributions obtained by these models are all obviously different from the case of even heat capacities except the FD scheme. This is against the heat transfer rules since the heat capacity does not influence the temperature distribution of statistic fluids under the equilibrium state. Among these results, the temperature distribution given by MITLB model changes most compared with Fig. 3, which shows a significant influence of heat capacity on the model. The simulated average heat flux and heat flux difference are also investigated in the present work as shown in Figs. 7 and 8. It can be seen from these figures that the simulated heat flux given by MITLB model is totally different from the analytic results and the heat fluxes differences between the side walls are quite large. These results show reasonable results cannot be obtained by MITLB model for the uneven heat capacity case. Compared with Figs. 4 and 5, the simulation results given by TTLB model and STTLB model are closer to the analytical heat flux. This is because the gradients of density and heat capacity are partly canceled out. However, the errors of these models are still quite large. It should also be noted that the simulation results given by presented FD scheme are still close to the analytic results and heat fluxes differences between the side walls are quite small, which shows the discretization model adopted in the present work a good act of accuracy and energy conservation.

652

A. Hu, R. Uddin and D. Liu / Computers and Fluids 179 (2019) 645–654

Fig. 6. Temperature distributions across the phase interface of uneven heat capacity.

The performance of the presented models in the multiphase flow heat transfer simulation is discussed in this section. According to these results, these LB thermal models are highly influenced by the heat capacity distribution. It should be noted that, in the latest research [28,46,47], additional terms are added to the right side of the equation of temperature based thermal LB model to correct its thermal diffusivity term and error terms. This is not recommended, as the additional terms need also to be discretized by the finite different method. 5. Comparation of energy equations of pseudo-potential phase change model

Fig. 7. The simulated average heat flux q¯ given by different models.

According to the simulation results in the last section, the presented FD scheme can accurately discrete the thermal equation. Hence, in this section, the FD scheme is adopted to further test the energy consistence of the pseudo-potential model in the phase change simulation. The coupling between the multiphase LB model for the flow field and the temperature field is established via the nonideal equation of state, and the Peng-Robinson equation of state of Ref. [16] is adopted:

p=

aα (T )ρ 2 ρ RT − 1 − bρ 1 + 2bρ − b2 ρ 2

where



Fig. 8. The simulated relate heat fluxes differences between the side walls given by different models.

)] 2 ,

(62)

α (T ) = [1 + (0.37464 + 1.54226ω0 − 0.26992ω02 ) × (1 −

T /Tc a = 0.45724R2 Tc2 / pc , b = 0.0778RTc / pc . The parameters are chosen asa = 3/49 b = 2/21, R = 1, ω0 = 0.344 [32]. To numerically compare the performance of these two energy equations in the simulation of phase change problems. In this work, the one-dimensional liquid film evaporation is simulated as shown in Fig. 9. The simulation domain is chosen as Lx × Ly = 20 × 100. Initially, a liquid film with thickness of 50 is placed at the bottom of the domain. The bottom boundary is set as a solid wall and a heat flux density of q is applied on the boundary, the half bounce back boundary condition is adopted to deal with the evolution of density distribution function. The open boundary condition of constant pressure is employed at the upper boundary where the density is set constant as ρ v and the temperature is set as the saturated temperature of 0.86Tc . The left and right boundaries are both set as periodic boundaries. The thermal conductivities are given as λl = 10 and λv = 1, and the thermal parameters in the interface are given by Eq. (59). According to energy conser-

A. Hu, R. Uddin and D. Liu / Computers and Fluids 179 (2019) 645–654

653

Fig. 9. Liquid film evaporation with constant heat flux density.

vation law, the phase change mass is given by:

m =

q t Lx , hlg

(63)

where hlg is the latent heat. Though Eqs. (13) and (16) are all deduced from Eq. (9), their expressions of the latent heat are actually quite different. For Eq. (13), the energy balance equation can be written as

Q pEOS = du + pdv = d(Cv T ) − 2 dρ . dt ρ

(64)

Since the temperature is constant during the phase change process, the latent heat is given by

hlg =

 ρv ρl

d(Cv T ) −

 ρv p ρl

ρ

dρ = (Cvv − Cvl )Tsat − 2

 ρv pEOS ρl

ρ2

dρ . (65)

It should be noted that the latent heat given by this equation is related to the heat capacities. Hence, the specific heat of the liquid cannot be set too large, or the hlg may become negative, since  ρv pEOS ρ ρ 2 dρ is constant at a fixed temperature of Tsat . The reason l

of the conflict is that the definition of u = Cv T cannot satisfy the thermodynamic consistence since Cv is artificially defined in the simulation. Similarly, we can also get the latent heat of Eq. (16), which is given by

hlg =

 vv T vl

∂p ∂T



dv = v

 ρv 1 ρl

ρ

T 2

∂p ∂T



dρ .

(66)

ρ

Compared with Eq. (65), the latent heat is independent with the heat capacities, the thermodynamic inconsistency of Eq (65) is actually avoided by using the temperature-based energy equation. To evaluate the effect of the specific heat at constant volume on the performance of the energy equations, two combinations of specific heat are adopted in the simulations: one is the uniform specific heats case of Cvl = Cvv = 1; while for the other, the specific heats are given by Cvl = 1.5, Cvv = 1, respectively. Fig. 10 shows the simulation results of these two energy equations. It can be found that the phase change mass linearly increases with the heat flux q, showing that the latent heat is independently with q. The simulations results given by the temperaturebased energy equation are independently with the latent heat, while the results given by the internal-energy-based energy equation is obviously changed with the latent heat. This is consistent with our previous analysis. The evaporation speed given by the he internal-energy-based energy equation is much larger than

Fig. 10. Comparison between LBM results and analytical results: (a) simulaiton results of the internal-energy-based energy equation ( t = 10 0 0 0δt ); (b) simulaiton results of the temperature-based energy equation ( t = 60 0 0 0δt ).

the temperature-based energy equation, since the latent heats of internal-energy-based energy equation is much smaller. It should be noted that all these simulation results are consistent with the theoretical analysis results of Eq. (63) except the ones of nonuniform Cv given by the internal-energy-based equation, which are obviously different from the theoretical one. The possible reason is that, since the parameters change dramatically across the interface, the term of ∇ Cv is difficult to be precisely discretized in just 4–6 lattice. In addition, the latent heat given by Eq. (65) is much smaller than the temperature-based energy equation, making it more easily to be influenced by the discretization errors. In this section, we compared these two most common used energy equations of the pseudo-potential phase change model, it is found that the temperature-based energy equation is more suitable in the simulation for both thermodynamic consistence and accuracy. 6. Conclusion In this work, we have theoretically analyzed four lattice Boltzmann thermal models and compared them with the FD scheme

654

A. Hu, R. Uddin and D. Liu / Computers and Fluids 179 (2019) 645–654

in the simulations of the cavity natural convection and onedimensional heat transfer across the interface of the static phases. The energy conservation of the pseudo-potential phase change model is them verified with the FD scheme which performed best in the previous simulations. Based on the simulation results, following conclusions can be made: (1) The presented LB thermal models all contain error terms which lead to computational errors in the simulation, while the FD scheme is much better in accuracy. (2) Compared with the internal-energy based energy equation, the temperature-based energy equation is more suitable in the simulation since it is better in accuracy and thermodynamic consistency. Acknowledgments This work is sponsored by National Natural Science Foundation of China under grant 51606159; and Sichuan Provincial Department of Education under 17ZA0405. References [1] McNamara GR, Zanetti G. Use of the Boltzmann Equation to Simulate Lattice– Gas Automata. Phys Rev Lett 1988;61(20):2332–5. [2] Guo Z, Shu C. Lattice Boltzmann method and its applications in engineering. World Scientific; 2013. [3] Gunstensen AK, et al. Lattice Boltzmann model of immiscible fluids. Phys Rev A 1991;43(8):4320–7. [4] Shan X, Doolen G. Multicomponent lattice-Boltzmann model with interparticle interaction. J Stat Phys 1995;81(1-2):379–93. [5] Shan X, Chen H. Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation. Phys Rev E 1994;49(4):2941–8. [6] Shan X, Chen H. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys Rev E 1993;47(3):1815–19. [7] Swift MR, et al. Lattice Boltzmann simulations of liquid-gas and binary fluid systems. Phys Rev E 1996;54(5):5041–52. [8] Swift MR, Osborn WR, Yeomans JM. Lattice Boltzmann simulation of nonideal fluids. Phys Rev Lett 1995;75(5):830–3. [9] He X, Chen S, Zhang R. A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh–Taylor instability. J Comput Phys 1999;152(2):642–63. [10] He X, Shan X, Doolen GD. Discrete Boltzmann equation model for nonideal gases. Phys Rev E 1998;57(1):R13–16. [11] Hu A, et al. Contact angle adjustment in equation-of-state-based pseudopotential model. Phys Rev E 2016;93(5):053307. [12] Li Q, Luo KH. Achieving tunable surface tension in the pseudopotential lattice Boltzmann modeling of interface dynamics. Phys Rev E 2013;88(5):053307. [13] Li Q, Luo KH, Li XJ. Lattice Boltzmann modeling of multiphase flows at large density ratio with an improved pseudopotential model. Phys Rev E 2013;87(5):053301. [14] Hu A, et al. On equations of state in pseudo-potential multiphase lattice Boltzmann model with large density ratio. Int J Heat Mass Transfer 2013;67(6):159–63. [15] Kupershtokh AL, Medvedev DA, Karpov DI. On equations of state in a lattice Boltzmann method. Comput Math Appl 2009;58(5):965–74. [16] Yuan P, Schaefer L. Equations of state in a lattice Boltzmann model. Phys Fluids 2006;18(4):42101. [17] Safari H, Rahimian MH, Krafczyk M. Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow. Phys Rev E 2013;88(1):013304. [18] Gong S, Cheng P. A lattice Boltzmann method for simulation of liquid–vapor phase-change heat transfer. Int J Heat Mass Transfer 2012;55(17–18):4923–7. [19] Hazi G, Markus A. On the bubble departure diameter and release frequency based on numerical simulation results. Int J Heat Mass Transfer 2009;52(5–6):1472–80.

[20] Dong Z, Li W, Song Y. Lattice Boltzmann simulation of growth and deformation for a rising vapor bubble through superheated liquid. Numer Heat Transfer Part A 2009;55(4):381–400. [21] Zhang R, Chen H. Lattice Boltzmann method for simulations of liquid-vapor thermal flows. Phys Rev E 2003;67(6):066711. [22] Shi Y, Zhao TS, Guo ZL. Thermal lattice Bhatnagar-Gross-Krook model for flows with viscous heat dissipation in the incompressible limit. Phys Rev E 2004;70(2):066310. [23] Shu C, Chew YT, Peng Y. Simplified thermal lattice Boltzmann model for incompressible thermal flows. Phys Rev E 2003;68(2):026701. [24] Guo Z, Shi B, Zheng C. A coupled lattice BGK model for the Boussinesq equations. Int J Numer Methods Fluids 2002;39(4):325–42. [25] Inamuro T, et al. A lattice Boltzmann method for a binary miscible fluid mixture and its application to a heat-transfer problem. J Comput Phys 2002;179(1):201–15. [26] Gillissen JJJ, et al. Lattice-Boltzmann-based two-phase thermal model for simulating phase change. Phys Rev E 2013;88(3):033302. [27] Li Q, Luo KH. Effect of the forcing term in the pseudopotential lattice Boltzmann modeling of thermal flows. Phys Rev E 2014;89(5):053022. [28] Li Q, Zhou P, Yan HJ. Improved thermal lattice Boltzmann model for simulation of liquid-vapor phase change. Phys Rev E 2017;96(6):063303. [29] Li Q, et al. Enhancement of boiling heat transfer using hydrophilic-hydrophobic mixed surfaces: A lattice Boltzmann study. Appl Thermal Eng 2018;132:490–9. [30] Yu Y, et al. Investigation of droplet evaporation on heterogeneous surfaces using a three-dimensional thermal multiphase lattice Boltzmann model.. Appl Thermal Eng, 2017;127:1346–54. [31] Li Q, Zhou P, Yan HJ. Pinning–Depinning mechanism of the contact line during evaporation on chemically patterned surfaces: A lattice Boltzmann study. Langmuir 2016;32(37):9389–96. [32] Li Q, et al. Lattice Boltzmann modeling of boiling heat transfer: The boiling curve and the effects of wettability. Int J Heat Mass Transfer 2015;85:787–96. [33] Fan L, Yu Z. Multirelaxation-time interaction-potential-based lattice Boltzmann model for two-phase flow. Phys Rev E 2010;82(4):046708. [34] Luo L, Lallemand P. Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability. Phys Rev E 20 0 0;61(6):6546–62. [35] He YL, et al. Improved axisymmetric lattice Boltzmann scheme. Phys Rev E 2010;81(5):056707. [36] Abraham J, McCracken ME. Multiple-relaxation-time lattice-Boltzmann model for multiphase flow. Phys Rev E 2005;71(3):036701. [37] Shan X, Chen H. Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation. Phys Rev E 1994;49(4):2941–8. [38] Shan X, Chen H. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys Rev E 1993;47(3):1815–19. [39] Kupershtokh AL, et al. Stochastic models of partial discharge activity in solid and liquid dielectrics. IET Sci Meas Technol 2007;1(6):303–11. [40] Wei Y, Qian Y. Lattice Boltzmann simulations for thermal vapor-liquid two-phase flows. J Comput Multiphase Flows 2012;4(1):14007. [41] Lee T, Lin C. A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio. J Comput Phys 2005;206(1):16–47. [42] Vierendeels J, Merci B, Dick E. Benchmark solutions for the natural convective heat transfer problem in a square cavity with large horizontal temperature differences. Int J Numer Methods Heat Fluid Flow 2003;13(8):1057–78. [43] Guo ZL, Zheng CG, Shi BC. Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method. Chin Phys 2002;11(4):366–74. [44] Barakos G, Mitsoulis E, Assimacopoulos D. Natural convection flow in a square cavity revisited: Laminar and turbulent models with wall functions. Int J Numer Methods Fluids 1994;18(7):695–719. [45] Colosqui CE, et al. Mesoscopic simulation of non-ideal fluids with self-tuning of the equation of state. Soft Matter 2012;8(14):3798–809. [46] Mu YT, et al. Nucleate boiling performance evaluation of cavities at mesoscale level. Int J Heat Mass Transfer 2017;106:708–19. [47] Fang WZ, et al. Lattice Boltzmann modeling of pool boiling with large liquid– gas density ratio. Int J Thermal Sci 2017;114:172–83.