Achieving thermodynamic consistency in a class of free-energy multiphase lattice Boltzmann models

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 83, 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 $\nabla \left( {\rho c_{\rm{s}}^{\rm{2}}} \right)$ 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.


I. Introduction
The lattice Boltzmann (LB) method [1][2][3][4][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][10][11], the free-energy LB model [12][13][14][15], the pseudopotential LB model [16][17][18][19][20], and the phase-field LB model [21][22][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][14][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.,  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 where I is the unit tensor. After some standard algebra, the following equation can be obtained [8,27] 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]: 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 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   Through the transformation matrix M , the right-hand side of Eq. (7), namely the collision step, can be implemented in the moment space: 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: The macroscopic density and velocity are calculated by 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]: 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]   where   are the weights. For the D2Q9 lattice,   are given by 0 4 9  

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: 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: According to the Taylor series expansion, we can obtain Substituting Eq. (18) into Eq. (17) gives 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 Combining Eq. (16) with Eq. (20), we can obtain the following numerical error term: 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 . 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.

A. Scheme I
By introducing a source term to the LB equation, the collision step in the moment space can be written as follows: where the source term Q is given by 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): With such an additional term, the major numerical error term given by Eq. (23) can be eliminated. Note that the discretization of 2   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   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: 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]   in which the correction term C is given by   T 1 7 0, 9 , 0, 0, 0, 0, 0, 3 , 0 But it should also be noted that the discretization of   2 s 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 , 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 Then the following chemical potential can be obtained according to Eq. (2): 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:   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 . 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 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 . 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.

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: t a n h 2 2  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 0.7 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.

V. Conclusions
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 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.