Improved thermal lattice Boltzmann model for simulation of liquid-vapor phase change

In this paper, an improved thermal lattice Boltzmann (LB) model is proposed for simulating liquid-vapor phase change, which is aimed at improving an existing thermal LB model for liquid-vapor phase change [S. Gong and P. Cheng, Int. J. Heat Mass Transfer 55, 4923 (2012)]. First, we emphasize that the replacement of \[{\left( {\rho {c_V}} \right)^{ - 1}}\nabla \cdot \left( {\lambda \nabla T} \right)\] with \[\nabla \cdot \left( {\chi \nabla T} \right)\] is an inappropriate treatment for diffuse interface modeling of liquid-vapor phase change. Furthermore, the error terms ${\partial_{t0}}\left( {Tv} \right) + \nabla \cdot \left( {Tvv} \right)$, which exist in the macroscopic temperature equation recovered from the standard thermal LB equation, are eliminated in the present model through a way that is consistent with the philosophy of the LB method. In addition, the discrete effect of the source term is also eliminated in the present model. Numerical simulations are performed for droplet evaporation and bubble nucleation to validate the capability of the improved model for simulating liquid-vapor phase change. Numerical comparisons show that the aforementioned replacement leads to significant numerical errors and the error terms in the recovered macroscopic temperature equation also result in considerable errors.

is an inappropriate treatment for diffuse interface modeling of liquid-vapor phase change. Furthermore, the error terms

I. INTRODUCTION
The lattice Boltzmann (LB) method, which originates from the lattice gas automata method [1], has been developed into an efficient numerical approach for a wide range of phenomena and processes in the past three decades [2][3][4][5][6][7][8][9]. The LB equation can be viewed as a special discrete solver for the kinetic Boltzmann equation with certain collision operator, such as the Bhatnagar-Gross-Krook (BGK) collision operator [3,10] and the multiple-relaxation-time (MRT) collision operator [11][12][13][14][15][16]. The fluid flow is simulated by tracking the evolution of the particle distribution function and then the distribution function is accumulated to obtain the macroscopic properties. The LB method is easy to parallelize and is far less costly in terms of data exchange owing to its explicit scheme and the local interactions.
In recent years, the LB simulations of liquid-vapor phase change have attracted much attention and three categories of thermal LB models have been developed for simulating liquid-vapor phase change. The first category is based on the phase-field multiphase LB method, such as the models developed by Dong et al. [17], Safari et al. [18,19], and Sun et al. [20]. In these models, the liquid-vapor interface is captured by solving an interface-capturing equation (e.g., the Cahn-Hilliard equation) and a source term is incorporated into the continuity equation or the interface-capturing equation to mimic the phase change. Hence the rate of the liquid-vapor phase change in these models is an artificial input.
The second category is based on the pseudopotential multiphase LB method, which is a very popular multiphase approach in the LB community [7]. In the pseudopotential multiphase LB method, the phase separation between different phases is achieved via an interparticle potential [21,22].
Therefore the liquid-vapor interface can naturally arise, deform, and migrate without using any interface-tracking or interface-capturing technique. The thermal multiphase LB models proposed by Zhang and Chen [23], Házi and Márkus [24,25], Biferale et al. [26], Gong and Cheng [27], Kamali et al. [28], and Li et al. [29] can be classified into this category. The third category is the multi-speed thermal LB method, which employs a single set of distribution functions like the standard isothermal LB method but utilizes more discrete velocities [30,31]. The equilibrium distribution function usually includes higher-order velocity terms so as to recover the energy equation. The thermal LB models presented by Gonnella et al. [32] and Gan et al. [33] for thermal liquid-vapor flows fall into this category.
In many of the aforementioned thermal multiphase LB models, a thermal LB equation is employed to recover a target temperature equation at the Navier-Stokes level. The target temperature equation is usually a convection-diffusion equation with a source term. Therefore a thermal LB equation with a source term was devised in these models. However, it has been widely found [34][35][36][37] that there exist error terms in the macroscopic equation recovered from the standard thermal LB equation, which should be treated using appropriate correction techniques. In addition, the temperature field can also be simulated by traditional numerical methods such as the finite-difference method. In Ref. [

A. The target temperature equation
Historically, the first thermal LB model for liquid-vapor phase change was proposed by Zhang and Chen [23]. In their work, the macroscopic energy equation was given by where e is the internal energy, V c is the specific heat at constant volume, and λ is the thermal conductivity. In 2009, Házi and Márkus [24] derived a target temperature equation from the local balance law for entropy [38] ( ) where s is the entropy and ( ) ( ) ( ) is the material derivative. The viscous heat dissipation has been neglected in Eq. (2). According to the thermodynamic relations of non-ideal gases, the following equation can be obtained: This equation can be found in Table 11.4-1 in Ref. [39]. For ideal gases ( EOS p RT ρ = ), the last term on the right-hand side of Eq. (4) reduces to EOS p ⋅ v ∇ . The above equation can also be written as follows: ( ) In the literature, some other forms of the energy equation for non-ideal fluids can also be found [40][41][42]. For example, Onuki [41,42] established a general equation for the total energy density of non-ideal fluids (see Eq. (9) in Ref. [41]), which can be transformed to the following equation for the internal energy density (see Eq. (2.40) in Ref. [42]): where ê e ρ = is the internal energy density ( e is the internal energy of non-ideal fluids), σ is the According to thermodynamics, the relationship between the internal energy and the entropy is given by Using Eq. (8), the internal energy equation (7) can be transformed to Substituting the continuity equation D Dt The term v σ : ∇ represents the viscous heat dissipation. Comparing Eq. (2) with Eq. (10), we can see that these two equations are basically consistent except that ( ) − v Τ σ : ∇ is neglected in Eq. (2). 6 For simplicity, Gong and Cheng [27] replaced

B. The Chapman-Enskog analysis of the Gong-Cheng model
is the thermal diffusivity. Then they established the following temperature equation: The corresponding thermal LB equation for Eq. (11) was given by [27] ( where g α is the temperature distribution function, α e is the discrete velocity in the α th direction, g τ is non-dimensional relaxation time for the temperature field, and the source term G α α ω φ = , in which φ represents the second term on the right-hand side of Eq. (11), namely [27] The equilibrium temperature distribution function eq g α was given by ( ) vv e e I e v (14) where I is the unit tensor, The macroscopic equation recovered from Eq. (12) can be derived through the Chapman-Enskog analysis, which can be implemented by introducing the following multi-scale expansions: where 0 t and 1 t are two different time scales, and t δ serves as the expansion parameter [43]. The Taylor series expansion of Eq. (12) yields With the help of Eq. (15), Eq. (16) can be rewritten in the consecutive orders of t δ as follows: 1 : Taking the summations of Eqs. (17) and (19), the following equations can be obtained, respectively: In the above derivations, the relations With the aid of Eq. (14), we have Substituting Eq. (23) into Eq. (21) gives Combining Eq. (20) with Eq. (24) through  e , while the last term on the right-hand side of Eq. (25) is caused by the discrete effect of the source term, which can be seen from Eqs. (19) and (21).
is an inappropriate treatment for multiphase flows. In fact, such a treatment requires that the following term can be neglected: For single-phase incompressible flows, the aforementioned replacement is applicable since the density variation is very small. For multiphase flows, the density varies significantly within the liquid-vapor interface, which usually has a thickness of 4−5 lattices in the LB simulations of multiphase flows.
Therefore the term given by Eq. (26) cannot be neglected at the liquid-vapor interface. Some researchers [44] found that under certain conditions ϕ is small in comparison with the source term φ given by Eq. (13). In fact, not only ϕ but also the thermal conductivity term can be small as compared with the term φ in Eq. (13), but it does not mean that ϕ or the thermal conductivity term in the temperature equation can be dropped. The comparison should be made between ϕ and the thermal conductivity term instead of comparing ϕ with the source term φ in Eq. (13), because it arises from the replacement of ∂ v is given as follows according to the Chapman-Enskog analysis of the LB equation for the flow field [36]: where F is the force acting on the system. Obviously,

III. Improved thermal LB model A. Theoretical analysis based on the BGK collision operator
The improved thermal LB model will be constructed based on the MRT collision operator. Before presenting the improved model, we would like to provide some analyses about removing the error terms in Eq. (25) within the framework of the BGK collision operator, which may be useful for general readers to better understand the improved thermal LB model in the next subsection. The target temperature equation given by Eq. (5) can be rewritten as follows [25]: The source term φ is now given by the underlined terms in Eq. (28).
According to Eqs. (22) and (23), the error term in Eq. (25) can be removed by dropping the second-order velocity terms in eq g α , and then eq g α becomes Meanwhile, the error term (25) can be eliminated by adding a correction term to the thermal LB equation where the correction term C α is given by which satisfies Theoretically, to remove the discrete effect of the source term, namely the error term 0 t φ ∂ in Eq.
(25), the source term G α in Eq. (30) should also contain the coefficient ( ) in the correction term given by Eq. (31), which has been extensively demonstrated in the literature when a forcing or source term is incorporated into the LB equation [45,46]. However, when this coefficient is placed in front of the source term, the temperature should be calculated by 0 the calculation of the temperature will become implicit and iterations will be required.
Hence another treatment is considered. If we retain the definition of the temperature T g α α =  , the source term should take the form of  (19). Therefore the source term can take the following form: Then the thermal LB equation becomes The Chapman-Enskog analysis can also be performed for Eq. (33). Using the multi-scale expansions given by Eq. (15), the correction term C α should be expanded as Then the following equations can be obtained: (1) (2) 0 1 Note that the last term on the right-hand side of Eq. (35) has been used to eliminate the same term generated on the left-hand side of Eq. (36). The summations of Eqs. (34) and (36) lead to, respectively In the above derivations, the relations Using Eq. (29), we have  It can be found that the following treatments have been employed in the above analyses. First, a correct target temperature equation is adopted. Second, the error terms removed by dropping the second-order velocity terms in the equilibrium temperature distribution function and adding a correction term to the thermal LB equation. Furthermore, the discrete effect of the source term is eliminated by incorporating an additional term into the thermal LB equation.

B. The improved thermal MRT-LB model
In this subsection, the improved thermal LB model is presented based on the MRT collision operator. Using the MRT collision operator, the thermal LB equation can be written as follows: where S α ′ is the source term in the discrete velocity space and  .
The second-order velocity terms in the equilibrium distribution function eq g α should be dropped to remove the error term The right-hand side of Eq. (42) can be implemented in the moment space as follows: where 0 0.5 t t S φ δ φ = + ∂ . As discussed in the previous subsection, the additional term 0.5 t t δ φ ∂ is used to eliminate the discrete effect of the source term. With Eq. (45), the streaming process is given by The above equilibria can also be found in Ref. [37]. Using Eq. (52), the modifications given by Eqs. (48) and (49) can be rewritten as follows: where the non-equilibrium parts and Huang and Wu [37].

A. Droplet evaporation
First, the well-known D2 law for droplet evaporation is considered, which predicts that the square of the droplet diameter changes linearly over time [18,49]. This law is established based on the  Fig. 2. For comparison, the available data in Ref. [50], which were obtained using a finite-difference scheme to solve the temperature equation (5) [50], which were obtained by a finite-difference scheme [29].
As mentioned earlier, the influence of the replacement of To further illustrate the above points, the droplet evaporation on a solid surface is also considered.
In the above test, the thermal conductivity is chosen to be constant according to the requirement of the D2 law. In the present test, the thermal conductivity is taken as  [29].
The variation of the droplet volume with time is shown in Fig. 4, where l.u. represent lattice units.
For comparison, the numerical results obtained by a finite-difference scheme for solving the temperature equation [29] are also shown in Fig. 4.

V. Conclusions
In this paper, we have presented an improved thermal LB model for simulating liquid-vapor phase change. The Chapman-Enskog analysis has been performed for the Gong-Cheng model, which shows that the term found that the numerical errors caused by the error terms are non-negligible. We believe that the theoretical analyses as well as the numerical results in the present paper are useful for clarifying some critical issues and the present study would be helpful for general readers to better understand the thermal LB models for liquid-vapor phase change. With the modifications given by Eqs. (54) and (55), the following equations can be obtained: