Achieving thermodynamic consistency in a class of free-energy multiphase lattice Boltzmann models
11 Achieving thermodynamic consistency in a class of free-energy multiphase lattice Boltzmann models
Q. Li * , Y. Yu, and R. Z. Huang School of Energy Science and Engineering, Central South University, Changsha 410083, China * Corresponding author: [email protected]
Abstract
The free-energy lattice Boltzmann (LB) model is one of the major multiphase models in the LB community. The present study is focused on a class of free-energy LB models in which the divergence of thermodynamic pressure tensor or its equivalent form expressed by the chemical potential is incorporated into the LB equation via a forcing term. Although this class of free-energy LB models may be thermodynamically consistent at the continuum level, it suffers from thermodynamic inconsistency at the discrete lattice level owing to numerical errors [Guo et al ., Physical Review E , 036707 (2011)]. The numerical error term mainly includes two parts, one comes from the discrete gradient operator and the other can be identified in a high-order Chapman-Enskog analysis. In this paper, we propose two schemes to eliminate the thermodynamic inconsistency of the aforementioned class of free-energy LB models. The first scheme is devised by removing the major numerical error term that causes the thermodynamic inconsistency, while the other scheme is constructed by modifying the equation of state of the standard LB equation, through which the discretization of c is no longer involved in the force calculation and then the numerical errors can be significantly reduced. Numerical simulations are subsequently performed to validate the proposed schemes. Both schemes are shown to be capable of eliminating the thermodynamic inconsistency and the latter scheme is found to be relatively more accurate. PACS number(s): 47.11.-j. I. Introduction
The lattice Boltzmann (LB) method [1-5], which is a mesoscopic numerical approach originating from the lattice gas automaton (LGA) method [6], has been proven to be particularly suitable for studying multiphase and multicomponent systems [7,8] where the interfacial dynamics and phase transition are present. In the past three decades, significant progress has been made in this direction and a variety of multiphase LB models have been developed, such as the color-gradient LB model [9-11], the free-energy LB model [12-15], the pseudopotential LB model [16-20], and the phase-field LB model [21-23]. Although these models are devised from different points of view, they share the same feature, i.e., in these models the interface between different phases/components is represented by a diffuse interface. An important advantage of diffuse interface models is that the motion of interface does not need to be tracked explicitly [8]. Among these multiphase LB models, the free-energy model proposed by Swift et al . [12,13] was devised based on thermodynamic theory. They proposed to modify the second-order moment of the equilibrium density distribution function so as to include a non-ideal thermodynamic pressure tensor. Therefore in this model the phase separation is described by a non-ideal equation of state in thermodynamic theory such as the van der Waals equation of state. However, the original free-energy multiphase LB model was shown to break Galilean invariance due to some non-Navier-Stokes terms recovered in the macroscopic momentum equation [24], which arises from incorporating the pressure tensor via the equilibrium density distribution function. To restore the Galilean invariance, several correction terms should be added to the equilibrium density distribution function [13-15,25]. In 2006, based on the consideration that the thermodynamics of a multiphase system can be equivalently taken into account through a forcing term, Wagner and Li [26] proposed a free-energy LB model that uses a forcing term to incorporate the divergence of thermodynamic pressure tensor ( P ) into the LB equation. In the meantime, a similar LB model was devised by Lee and Fischer [27] and they found that the spurious currents can be considerably reduced when the divergence of pressure tensor is expressed by a chemical-potential form, i.e., c P , where c is the chemical potential. This class of free-energy LB models can be referred to as “forcing-based” free-energy models in comparison with the original free-energy model proposed by Swift et al . [12,13]. This forcing-based free-energy LB method has been recently extended to multiphase systems with large density ratios by Mazloomi M et al . [28] and Wen et al . [29]. Nevertheless, Guo et al . [30] found that, for the aforementioned class of free-energy LB models, the force balance condition does not hold at the discrete lattice level regardless of using the pressure-tensor form P or the chemical-potential form c . Subsequently, Lou and Guo [31] showed that the force imbalance at the discrete level leads to thermodynamic inconsistency, i.e., the coexisting liquid and gas densities given by the forcing-based free-energy LB models gradually deviate from the results of the Maxwell construction when the reduced temperature decreases. As shown in Refs. [30,31], the force imbalance or the thermodynamic inconsistency of the forcing-based free-energy LB models is caused by the numerical errors at the discrete lattice level. To reduce the effects of force imbalance, Lou and Guo proposed [31] a Lax-Wendroff scheme for the forcing-based free-energy LB models and numerically demonstrated that the thermodynamic inconsistency can be eliminated using the Lax-Wendroff scheme. In a similar way, Qiao et al . [32] developed a thermodynamic-consistent free-energy LB model based on a Beam-Warming scheme. However, these two schemes are inconsistent with the philosophy of the standard LB method (the standard streaming-collision procedure). Moreover, the implementation of Beam-Warming scheme involves not only the nearest-neighbor nodes but also the next-nearest-neighbor nodes. In this paper, we aim to propose alternative schemes that are still within the framework of the standard streaming-collision procedure for eliminating the thermodynamic inconsistency of the forcing-based free-energy LB models. Specifically, two schemes are proposed. The first scheme is devised by removing the major numerical error term that causes the thermodynamic inconsistency, while the other scheme is constructed by modifying the equation of state of the standard LB equation, through which the discretization of s c is no longer required in the force calculation and then the numerical errors can be significantly reduced. The rest of the present paper is organized as follows. In Sec. II, the forcing-based free-energy LB method is briefly introduced. The schemes are presented in Sec. III and the numerical validation is provided in Sec. IV. Finally, a summary is given in Sec. V. II. Forcing-based free-energy LB method
A. Basic formulations
For multiphase systems, the free energy functional can be expressed by [12,13,33-36] , d d2 V f V E , (1) where V is the region of space occupied by the system, , is the free-energy density, f E represents the bulk free-energy density, which leads to an equation of state that allows for the coexistence of liquid and gas phases, and denotes the interfacial free-energy density, in which is a positive constant. The chemical potential c is defined as the variation of the free energy functional with respect to the density, i.e., δδ f E , (2) where d d f f E E . Then a non-local pressure can be defined as follows [37]: c , p . (3) Substituting Eq. (2) and the expression of , into Eq. (3) yields p p , (4) where EOS f f p E E is the non-ideal equation of state. Correspondingly, the thermodynamic pressure tensor is defined as [37-39] , p p P I I , (5) where I is the unit tensor. After some standard algebra, the following equation can be obtained [8,27] c P . (6) The left-hand side of Eq. (6) is usually referred to as the pressure-tensor form, whereas the right-hand side is referred to as the chemical-potential form. For the forcing-based free-energy LB method, the thermodynamics of a multiphase system is taken into account through a forcing term [26,27]. Without loss of generality, in the present study we adopt the LB equation with a multiple-relaxation-time (MRT) collision operator [40,41] based on the consideration that the single-relaxation-time (SRT) collision operator [42] can be viewed as a special case of the MRT collision operator. Generally, the MRT-LB equation can be written as follows [43,44]: ,, , , 0.5 eqt t t tt f t f t f f G G xx x e x , (7) where f is the density distribution function, eq f is the equilibrium density distribution function, x is the spatial position, t is the time, t is the time step, e is the discrete velocity in the th direction, G is the forcing term in the discrete velocity space, and M M is the collision operator, in which M is a transformation matrix and is a diagonal matrix for relaxation times [43,44]. For the two-dimensional nine-velocity (D2Q9) lattice, the relaxation matrix is given by diag , , , , , , , , c e c q c q v v , (8) where c is the non-dimensional relaxation time related the conserved moments, and q are free parameters, while e and v are the non-dimensional relaxation times determining the bulk and shear viscosities, respectively, i.e.,
2B s e t c and v t c . Through the transformation matrix M , the right-hand side of Eq. (7), namely the collision step, can be implemented in the moment space: eq t m m m m I S , (9) where m Mf , eq eq m Mf , and S MG is the forcing term in the moment space. Then the streaming step can be implemented as follows: , , t t f t f t x e x , (10) where f M m and M is the inverse matrix of the transformation matrix. For the D2Q9 lattice, the standard equilibria eq m are given by T2 2 2 2
1, 2 3 , 1 3 , , , , , , eq x x y y x y x y u u u u u u u u u u m , (11) where x y u u u . The macroscopic density and velocity are calculated by , 2 t f f u e F , (12) where F is the force exerted on the system. For the forcing-based free-energy LB method, the force is defined as follows [26,27]: or c c F P F , (13) where c . The left one is the pressure-tensor form, whereas the right one is the chemical-potential form. The forcing term G in Eq. (7) is given by [45] G c c e ue u e F , (14) where are the weights. For the D2Q9 lattice, are given by , , and . With the aid of Eq. (14), the forcing term in the moment space can be obtained via S MG , in which
T0 1 8 , , ,
G G G G . B. Numerical error term at the discrete level
The macroscopic equations recovered from the forcing-based free-energy LB method can be derived by the Taylor expansion analysis [46] or the Chapman-Enskog analysis [47]. As argued by Wagner [48], a second-order analysis is inadequate for multiphase LB models because higher-order terms are ignored, which may be necessary to achieve thermodynamic consistency. Through a third-order Chapman-Enskog analysis [49], the following macroscopic equations can be obtained: t u , (15) t ct u uu F F , (16) where is the viscous stress tensor. The last term on the right-hand side of Eq. (16) is a high-order term that cannot be identified by a second-order Chapman-Enskog analysis. Meanwhile, in numerical implementation the gradient terms are usually evaluated by the following isotropic finite-difference scheme: discrete tt c x e e . (17) According to the Taylor series expansion, we can obtain t tt k t k k l k l k l m k l m e e e e e e x e x x x x . (18) Substituting Eq. (18) into Eq. (17) gives discrete t . (19) In the LB method we usually adopt 1 t c , where c is the lattice constant. Hence the force at the discrete level can be expressed via discrete F F F . (20) Combining Eq. (16) with Eq. (20), we can obtain the following numerical error term: F F . (21) At the liquid-gas interface, the variation of the chemical potential is usually much smaller than that of the density, i.e., c . Accordingly, the major numerical error term is given by c c . (22) Since j i i i i j , we can rewrite Eq. (22) as c . (23) Note that I . Consequently, the isotropic part of the thermodynamic pressure tensor is changed by the numerical error term, which leads to the thermodynamic inconsistency of the forcing-based free-energy LB models. III. Alternative schemes for achieving thermodynamic consistency
In this section, two schemes that are within the framework of the standard streaming-collision procedure are proposed to eliminate the thermodynamic inconsistency of the forcing-based free-energy LB models. According to the analysis in the previous section, the first scheme is devised by adding a source term to the LB equation to remove the major numerical error term that causes the thermodynamic inconsistency, i.e., the error term given by Eq. (23). On the other hand, if we modify the equation of state produced by the standard LB equation, the discretization of c will be no longer involved in the force calculation. Then the numerical errors can be significantly reduced. Following such a line, another scheme is also constructed. A. Scheme I
By introducing a source term to the LB equation, the collision step in the moment space can be written as follows: eq t t m m m m I S Q , (24) where the source term Q is given by T2 22 2s s c c Q . (25) Through the Chapman-Enskog analysis, it can be verified that the source term given by Eq. (25) adds the following term to the right-hand side of Eq. (16): c c A I . (26) With such an additional term, the major numerical error term given by Eq. (23) can be eliminated. Note that the discretization of is involved in the calculation of the chemical potential c . Hence no additional discretization is required when using this scheme. In the remaining of this paper, Eqs. (24) and (25) are referred to as the scheme I. B. Scheme II
In fact, the appearance of c in the force given by Eq. (13) is attributed to the equation of state p c produced by the standard LB equation, which yields the term c on the right-hand side of Eq. (16). Obviously, when such an equation of state is modified as m p p , the term c in Eq. (13) should be replaced by m p . In the present study we choose m c p . Then the chemical-potential form of the force is given by m c c p F . (27) For the sake of changing the equation of state produced by the standard LB equation, the second-order moment of the equilibrium density distribution function is modified as follows: m ijeq ii j j ue e f p u , (28) where ij is the Kronecker delta. However, it is well-known that, when Eq. (28) is implemented without other changes, the LB model will suffer from the lack of Galilean invariance owing to some non-Navier-Stokes terms recovered in the macroscopic momentum equation. Such an issue can be found in the original free-energy LB model and the color-gradient LB models with variable density ratios [10,50]. Several approaches have been devised in the literature for addressing this issue [13,14,25,51]. In the present study, we adopt the approach proposed by Li et al . [51], which has been recently applied to eliminate the error terms of color-gradient LB models [52,53]. Following this approach, the third-order moment of the equilibrium density distribution function is given by [51] , if ,, others. i jk j ik k ijeqi j k i jk j ik k ij c u u u i j ke e e f p u u u (29) According to Eqs. (28) and (29), the equilibria eq m in the moment space are now given by T2 2 2 2
1, 4 3 2 , 4 3 3 , , 2 , , 2 , , eq x x y y x y x y u u u u u u u u u u m , (30) where m p . From Eq. (29) it can be seen that the diagonal elements ( i j k ) of the third-order moment deviate from the require relationship. As a result, there are still a couple of error terms that should be removed by introducing a correction term into the LB equation, i.e., eq t m m m m I S C , (31) in which the correction term C is given by T1 7
0, 9 , 0, 0, 0, 0, 0, 3 , 0
C C C with x x y y C u u , . x x y y C u u (32) where
2s m c p . The dynamic shear viscosity is now given by m v t p . To sum up, Eqs. (27), (30), and (31) constitute the scheme II for eliminating the thermodynamic inconsistency of the forcing-based free-energy LB models. Compared with the scheme I, the present scheme introduces some additional computations, namely the discretization of x x u and y y u . But it should also be noted that the discretization of c is no longer needed when using the scheme II. Besides, owing to such a feature, the scheme II is found to be relatively more accurate than the scheme I, which will be illustrated in the following part. IV. Numerical validation
In the preceding section, we have proposed two schemes for eliminating the thermodynamic inconsistency of the forcing-based free-energy LB models. In this section numerical simulations are carried out to validate the proposed schemes.
A. Flat interface
Firstly, the test of one-dimensional flat interfaces is considered. The grid system is taken as 100 100 x y
L L with the periodic boundary condition being applied in the x and y directions. Initially, the flat interfaces are placed at 0.25 x x L and 0.75 x x L , namely the central region is filled with the liquid phase, while the rest is occupied by the gas phase. The van der Waals equation of state is employed [12,13,54], which corresponds to the bulk free-energy density ln 1 f E RT b a . Then the following chemical potential can be obtained according to Eq. (2): RT ab b . (33) where a is the attraction parameter, b is the repulsion parameter, and R is the gas constant. In our simulations, the parameters are chosen as follows: 9 392 a , 2 21 b , 1 R , and . The critical density and temperature are given by c b and c T a Rb . Figure 1 displays the numerical coexistence curves predicted by different forcing-based free-energy LB models in the cases of 0.15 and 0.03 , where is the kinematic shear viscosity. For comparison, the analytical coexistence curve given by the Maxwell construction is also shown there. From the figure it can be clearly observed that the numerical coexistence curve predicted by the standard forcing-based free-energy LB model gradually deviates from the analytical one as the reduced temperature decreases. In contrast, the numerical coexistence curves obtained by the schemes I and II are in good agreement with the analytical one except that some deviations are observed in the case of 0.03 for the scheme I at low temperatures, which may be attributed to the influences of higher-order error terms yielded by the discretization of c . (a) (b) FIG. 1 . Comparison of the numerical coexistence curves with the analytical coexistence curve given by the Maxwell construction. The kinematic viscosity is (a) and (b) . -1 c T / T c Maxwell construction Standard FE model FE model with scheme I FE model with scheme II -1 c T / T c Maxwell construction Standard FE model FE model with scheme I FE model with scheme II For a flat interface, the chemical potential is theoretically constant across the interface as no curvature effect exists [31]. However, the standard forcing-based free-energy LB model results in a non-constant chemical potential for flat interfaces due to the force imbalance at the discrete lattice level, as shown in Fig. 2, where the chemical potential profiles obtained by different forcing-based free-energy LB models are compared at c T T with . It is seen that the chemical potential predicted by the standard forcing-based free-energy LB model varies significantly across the interfaces. Contrarily, the free-energy LB model with the scheme II produces a rigorously uniform chemical potential ( c ) over the whole domain, and the chemical potential given by the scheme I is basically constant except for slight variations at the interfaces. c x / L x Standard FE model FE model with scheme I FE model with scheme II
FIG. 2 . Comparison of the chemical potential profiles predicted by different forcing-based free-energy LB models at c T T . Furthermore, it is noticed that the numerical error term given by Eq. (23) not only alters the coexisting liquid and gas densities but also affects the interface thickness. Figure 3 displays the density profiles obtained by different forcing-based free-energy LB models at c T T with 0.15 . From the figure we can see that the interface thickness produced by the standard forcing-based free-energy LB model is much thicker than that given by the schemes I and II. Such a phenomenon indicates that the major numerical error term described by Eq. (23) serves as a numerical dissipation term for the forcing-based free-energy LB models, which smoothes the liquid-gas interface. Here it is also worth mentioning that the interface in simulations usually becomes thinner with the decrease of the reduced temperature (see Fig. 3 in Ref. [19] for details), but it can be widened by adjusting the parameter a in the non-ideal equation of state. x / L x FE model with scheme I FE model with scheme II Standard FE model
FIG. 3 . Comparison of the density profiles across the liquid-gas interface obtained by different forcing-based free-energy LB models at c T T . B. Circular droplet
In this subsection, numerical simulations are performed for the problem of two-dimensional circular droplet. In this test, a liquid droplet with radius r is placed at the center of a square domain and the rest of the domain is filled with the gas phase. The grid system is chosen as 120 120 x y L L with the periodic boundary condition being applied in the x and y directions. The density field is initialized as follows:
2, tanh2 2 l g l g r rx y W , (34) where l and g are the densities of the liquid and gas phases, respectively, W is the initial interface thickness, and r x x y y , in which , x y is the center of the computational domain. For circular droplets, the liquid and gas densities at equilibrium state can be obtained according to the non-ideal equation of state and the Young-Laplace’s law [31]. The equilibrium liquid and gas densities of flat surfaces correspond to the case of a circular droplet with the radius r . As the droplet size decreases, both the liquid and gas densities increase above the respective flat interface values according to the Young-Laplace’s law, which gives the following expressions for the pressures inside and outside of the circular droplet, respectively [55]: satsat sat sat ll l l g p p r , satsat sat sat gg g l g p p r , (35) where l p and g p are the pressures of the liquid and gas phases, respectively, is the surface tension, and the superscript sat denotes the properties of flat interfaces given by the Maxwell construction. When the surface tension and the droplet radius r are known, the pressures l p and g p can be determined at a given reduced temperature ( c T T ). Then the liquid and gas densities can be theoretically obtained by solving
EOS l l p p and EOS g g p p , respectively. In the present test, the parameters a , b , R , and are the same as those used in the previous section. The liquid and gas densities obtained by different forcing-based free-energy LB models are displayed in Figs. 4 and 5 for the cases of c T T and 0.7 , respectively. The results of the standard forcing-based free-energy LB model are unavailable at c T T as the model is unstable in this case. The kinematic viscosity is chosen as . For comparison, the theoretical results are also presented in the figures. From Fig. 4 it can be seen that the liquid and gas densities predicted by the standard forcing-based free-energy LB model significantly deviate from those of the theoretical solution. In contrast, Figs. 4 and 5 show that the numerical results given by the scheme II are in excellent agreement with the theoretical ones in both cases. Meanwhile, it is also observed that the numerical results obtained by the scheme I basically agree well with the theoretical ones except that some visible deviations are found in the case of c T T for the gas density. Specifically, in the case of c T T the maximum relative errors yielded by the standard forcing-based free-energy LB model, the scheme I, and the scheme II are about 35.3%, 0.85% and, 0.23%, respectively. Obviously, the free-energy LB models with the present schemes perform much better than the standard forcing-based free-energy LB model. Furthermore, similar to the previous test, the present test also shows that the scheme II is relatively more accurate than the scheme I. Such a difference may arise from the fact that the discretization of c is not involved in the force calculation of the scheme II.
20 25 30 35 406.46.56.66.76.86.9 li qu i d den s i t y r Theoretical Standard FE model FE model with scheme I FE model with scheme II
20 25 30 35 400.50.60.70.80.9 ga s den s i t y r Theoretical Standard FE model FE model with scheme I FE model with scheme II
FIG. 4 . Simulations of circular droplets. Comparison of the liquid and gas densities obtained by different forcing-based free-energy LB models at c T T .
20 25 30 35 407.47.57.67.7 li qu i d den s i t y r Theoretical FE model with scheme I FE model with scheme II
20 25 30 35 400.20.30.40.50.60.7 ga s den s i t y r Theoretical FE model with scheme I FE model with scheme II
FIG. 5 . Simulations of circular droplets. Comparison of the liquid and gas densities obtained by different forcing-based free-energy LB models at c T T . In LB literature [14,25] it has been reported that a circular droplet in a uniform flow field will become an elliptic one when employing a two-phase LB model with broken Galilean invariance. As discussed in the previous section, when Eq. (28) is implemented without other changes, the LB model will suffer from the lost of Galilean invariance. To verify the Galilean invariance of the scheme II, a moving circular droplet in a two-dimensional channel is simulated. The grid system is also chosen as 120 120 x y L L and a circular droplet of 25 r is initially placed at the center of the computational domain. The parallel top and bottom plates in the y direction begin to move with a constant velocity 0.1 U at 0 t . The no-slip boundary scheme [56] is applied at the plates and the periodic boundary condition is employed in the x direction. The reduced temperature and the kinematic viscosity are taken as c T T and 0.15 , respectively. Figure 6 displays some snapshots of a moving circular droplet simulated by the scheme II. It can be seen that the circular shape of the droplet is well preserved during the simulation, confirming the Galilean invariance of the scheme II. FIG. 6 . Density contours of a moving circular droplet simulated by the free-energy LB model with the scheme II. (a) 20000 t t , (b) 40000 t , (c) 70000 t , and (d) 90000 t . V. Conclusions a b c d In this paper, we have investigated the problem of thermodynamic inconsistency of a class of free-energy LB models in which the divergence of thermodynamic pressure tensor or its equivalent form expressed by the chemical potential is incorporated into the LB equation via a forcing term. It is shown that the numerical error term that causes the thermodynamic inconsistency mainly includes two parts, one comes from the discrete gradient operator and the other can be identified in a high-order Chapman-Enskog analysis. Two schemes that are within the framework of the standard streaming-collision procedure have been proposed to eliminate the thermodynamic inconsistency of the forcing-based free-energy LB models. The scheme I is devised by removing the major numerical error term that causes the thermodynamic inconsistency, while the scheme II is constructed by modifying the equation of state of the standard LB equation, through which the discretization of s c is no longer required in the force calculation. Numerical simulations have been performed for one-dimensional flat interfaces and two-dimensional circular droplets to validate the proposed schemes. The two schemes are shown to be capable of eliminating the thermodynamic inconsistency of the forcing-based free-energy LB models and the scheme II is found to be relatively more accurate. Acknowledgments
This work was supported by the National Natural Science Foundation of China (No. 51822606).
References [1] S. Chen and G. D. Doolen, Annual Review of Fluid Mechanics , 329 (1998). [2] C. K. Aidun and J. R. Clausen, Annual Review of Fluid Mechanics , 439 (2010). [3] S. Succi, Europhysics Letters , 50001 (2015). [4] Z. Guo and C. Shu, Lattice Boltzmann Method and Its Applications in Engineering (World Scientific, 2013). [5] Q. Li, K. H. Luo, Q. J. Kang, Y. L. He, Q. Chen, and Q. Liu, Progress in Energy and Combustion Science , 62 (2016). [6] U. Frisch, B. Hasslacher, and Y. Pomeau, Physical Review Letters , 1505 (1986). [7] H. Huang, M. C. Sukop, and X. Y. Lu, Multiphase Lattice Boltzmann Methods: Theory and Application (John Wiley & Sons, 2015). [8] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen,
The Lattice Boltzmann Method - Principles and Practice (Springer Nature, 2017). [9] A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti, Physical Review A , 4320 (1991). [10] D. Grunau, S. Chen, and K. Eggert, Physics of Fluids A , 2557 (1993). [11] M. Latva-Kokko and D. H. Rothman, Phys. Rev. E , 056702 (2005). [12] M. R. Swift, W. R. Osborn, and J. M. Yeomans, Physical Review Letters , 830 (1995). [13] M. R. Swift, E. Orlandini, W. R. Osborn, and J. M. Yeomans, Phys. Rev. E , 5041 (1996). [14] T. Inamuro, N. Konishi, and F. Ogino, Computer Physics Communications , 32 (2000). [15] C. M. Pooley and K. Furtado, Phys. Rev. E , 046702 (2008). [16] X. Shan and H. Chen, Phys. Rev. E , 1815 (1993). [17] X. Shan and H. Chen, Phys. Rev. E , 2941 (1994). [18] X. Shan, Phys. Rev. E , 066702 (2008). [19] Q. Li, K. H. Luo, and X. J. Li, Phys. Rev. E , 053301 (2013). [20] M. Sbragaglia and D. Belardinelli, Phys. Rev. E , 013306 (2013). [21] T. Lee and L. Liu, Journal of Computational Physics , 8045 (2010). [22] A. Fakhari and M. H. Rahimian, Phys. Rev. E , 036707 (2010). [23] A. Fakhari and D. Bolster, , 620 (2017). [24] A. Kuzmin, A. A. Mohamad, and S. Succi, International Journal of Modern Physics C , 875 (2008). [25] A. N. Kalarakis, V. N. Burganos, and A. C. Payatakes, Phys. Rev. E , 056702 (2002). [26] A. J. Wagner and Q. Li, Physica A Statistical Mechanics and IIts Applications , 105 (2006). [27] T. Lee and P. F. Fischer, Phys. Rev. E , 046709 (2006). [28] A. Mazloomi M, S. S. Chikatamarla, and I. V. Karlin, Physical Review Letters , 174502 (2015). [29] B. Wen, L. Zhao, W. Qiu, Y. Ye, and X. Shan, Phys. Rev. E , 013303 (2020). [30] Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E , 036707 (2011). [31] Q. Lou and Z. Guo, Phys. Rev. E , 013302 (2015). [32] Z. Qiao, X. Yang, and Y. Zhang, Int. J. Heat Mass Transf. , 1216 (2019). [33] J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity (Oxford, Clarendon Press, 1982). [34] O. Penrose and P. C. Fife, Physica D: Nonlinear Phenomena , 44 (1990). [35] J. W. Cahn and J. E. Hilliard, The Journal of Chemical Physics , 258 (1958). [36] J. W. Cahn and J. E. Hilliard, The Journal of Chemical Physics , 688 (1959). [37] Evans and R., Advances in Physics , 143 (1979). [38] D. N. Siebert, P. C. Philippi, and K. K. Mattila, Phys. Rev. E , 053310 (2014). [39] D. M. Anderson, G. B. McFadden, and A. A. Wheeler, Annual Review of Fluid Mechanics , 139 (1998). [40] P. Lallemand and L.-S. Luo, Phys. Rev. E , 6546 (2000). [41] L.-S. Luo, W. Liao, X. Chen, Y. Peng, and W. Zhang, Phys. Rev. E , 056710 (2011). [42] Y. H. Qian, D. d'Humières, and P. Lallemand, Europhysics Letters , 479 (1992). [43] M. E. McCracken and J. Abraham, Phys. Rev. E , 036701 (2005). [44] K. N. Premnath and J. Abraham, Journal of Computational Physics , 539 (2007). [45] Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E , 046308 (2002). [46] D. Lycett-Brown and K. H. Luo, Phys. Rev. E , 023305 (2015). [47] S. Chapman and T. G. Cowling, The Mathematical Theory of Non-uniform Gases (Cambridge University Press, Cambridge, 1970). [48] A. J. Wagner, Phys. Rev. E , 056703 (2006). [49] R. Huang and H. Wu, Journal of Computational Physics , 121 (2016). [50] H. Liu, A. J. Valocchi, and Q. Kang, Phys. Rev. E , 046309 (2012). [51] Q. Li, K. H. Luo, Y. L. He, Y. J. Gao, and W. Q. Tao, Phys. Rev. E , 016710 (2012). [52] Y. Ba, H. Liu, Q. Li, Q. Kang, and J. Sun, Phys. Rev. E , 023310 (2016). [53] Z. X. Wen, Q. Li, Y. Yu, and K. H. Luo, Phys. Rev. E , 023301 (2019). [54] B. Wen, X. Zhou, B. He, C. Zhang, and H. Fang, Phys. Rev. E (2016). [55] E. S. Kikkinides, A. G. Yiotis, M. E. Kainourgiakis, and A. K. Stubos, Phys. Rev. E , 036702 (2008). [56] Q. Li, K. Luo, Q. Kang, and Q. Chen, Phys. Rev. E90