## Abstract

The thermodynamic consistency of kinetic models for non-ideal mixtures in non-isothermal conditions is investigated. A kinetic model is proposed that is suitable for deriving high-order lattice Boltzmann equations by an appropriate discretization of the velocity space, satisfying the Galilean invariance condition and free of spurious terms in the first moment equations.

## 1. Introduction

From a molecular standpoint, the interface between two immiscible fluids is a transition layer, where the molecules of each fluid are subjected to two opposing processes: the random thermal motion of the molecules, responsible for mixing, and the intermolecular attractive forces on each molecule from its own phase, promoting segregation.

The relative importance of these two competing effects is mainly dependent on the temperature. Indeed, miscibility is not an intrinsic property of a pair of liquids, but is dependent on the thermodynamic state of the mixture. Two liquids that can mix in all proportions forming homogeneous solutions at a giving temperature can partially segregate when the temperature is decreased. The amount of solute that remains in the solution at this lower temperature is called the solubility and this last property is dependent on the temperature and pressure.

In non-equilibrium conditions, the transport equations for non-ideal fluids and their mixtures, in the miscibility range of these fluids, are usually written with the same form as in the ideal case, with an appropriate equation of state for the thermodynamic pressure. These macroscopic equations have not, nevertheless, been able to describe the segregation process that happens when two immiscible or partially miscible fluids are mixed, and this subject has been treated in a completely distinct framework.

In the lattice Boltzmann (LB) framework, immiscible flows have usually been treated by using ad hoc assumptions, e.g. Gunstensen *et al.* [1] and Gunstensen & Rothman [2] with their colour model, Shan & Chen [3,4] and the effective-mass model, and Santos *et al.* [5] with a model based on field mediators. A model based on a temperature shift besides Shan & Chen’s velocity shift was proposed by Sbrabaglia *et al.* [6] for non-isothermal two-phase flows. Free-energy approaches were proposed for multi-component and more frequently for multi-phase systems [7,8], and ad hoc assumptions were proposed [9–12] for correcting the energy transport, for reducing the spurious currents and for achieving Galilean invariance. Lee & Fischer [13] proposed a higher-order finite-difference treatment of the kinetic terms arising from molecular interactions, and Shan [14] showed that spurious currents are due to the insufficient isotropy of the discrete gradient operator.

Owing to the microscopic nature of intermolecular forces, LB models appear to be very suitable for producing mesoscopic methods that can improve the understanding of complex physical phenomena that are very difficult to describe at the hydrodynamic scale.

The description of real fluids requires consideration of the attractive forces among particles. He & Doolen [15] considered that the integral term in the Boltzmann kinetic equation was composed of two terms, in accordance with the range of integration in physical space, split by the radius *σ* of the sphere of influence around the target particle. The first term within this sphere is a short-range repulsion term corrected by considering Enskog’s volume exclusion effect. The second one ranging outside *σ* is a long-range attraction term.

These short-range repulsion and long-range attraction terms are the subject of the present work, when two or more kinds of particles are involved. The purpose of this paper is to derive a lattice Boltzmann equation (LBE), for non-ideal mixtures in non-isothermal problems, aiming to give a unified approach that will enable the influence of temperature on segregation to be understood. In this attempt, the thermodynamic consistency of the derived kinetic equations is thoroughly investigated in this paper. Owing to space limitations, only the kinetic equation in continuous space is presented in this paper, the LBEs themselves being considered as discrete forms of the continuous kinetic model, after using an appropriate velocity discretizaton scheme [16–18]. In this respect, it is important to distinguish errors due to the kinetic model used for non-isothermal segregation from errors due to the discretization procedure. The great majority of presently known multi-component, multi-phase models were planned for isothermal conditions and do not satisfy the condition of invariance of the total energy and, consequently, are not thermodynamically consistent. The conceptual problems related to the continuous models appear to be major questions to be solved for this class of problems.

## 2. Kinetic equation derived from the Liouville equation

Considering the mixture as a mechanical system of material points belonging to *r* different species, a kinetic model was derived for non-ideal multi-component mixtures starting from the Liouville equation. Defining *f*^{(p)} as the probable *amount* of *p* particles inside ** x**±

*d*

**with velocity**

*x***±**

*ξ**d*

**, it can be concluded that the evolution equation for**

*ξ**f*

^{(p)}can be written as 2.1Here

*g*^{(p,e)}is the acceleration due to an external conservative field acting on the

*p*particles and 2.2is the

*p*–

*s*interaction term related to the force

*g*^{(ps)}

_{η}acting on particles

*p*inside

**±**

*x**d*

**due to the interaction potential with particles**

*x**s*located in an elementary volume around

**+**

*x***, 2.3where**

*η**Φ*

^{(sk)}is supposed to be a central potential depending solely on the distance between particles

*p*and

*s*, and

*f*

^{(2,ps)}is related to the joint probability of finding a particle of species

*p*in position

**with velocity**

*x***and a particle**

*ξ**s*in position

**+**

*x***with velocity**

*η*

*ξ*_{1}.

Since the spatial variable *η* in the integral equation (2.2) has an infinite range, each *Ω*^{(ps)} integral can be split into two integrals,
2.4the first one in the range (0,*σ*_{ps}) and the second one in the range . When *σ*_{ps} is chosen as the position where the two-particle *p*–*s* interaction potential changes from repulsion to attraction [15], the first integral *Ω*^{(ps)}_{rep} can be considered as a classical Boltzmann collision term, whereas the second integral *Ω*^{(ps)}_{att} is to be related to the attractive forces on the *p* particle from the *s* particles that are outside the *σ*_{ps} sphere.

The repulsion integral term *Ω*^{(ps)}_{rep} can be further updated, using the Enskog theory for dense gases [19], which considered the molecules to be smooth rigid spheres. In this paper, this term will be written considering the shielding and amplifying effects on each *p* particle to be given as the overall effect of all *s* particles around it,
2.5with
2.6Here , *b*^{(s)} is the molecular volume that is attributed to molecules *s*, *x*^{(s)}=*n*^{(s)}/*n* is the number fraction of molecules *s*, *χ* is a dimensionless volume parameter dependent on *bn*, with *χ*=1 when *b*=0, *T* is the temperature,
2.7is the centre-of-mass velocity, with , and the dimensionless peculiar velocity is given by
2.8

The collision term *Ω*^{(ps)}_{rep}(*b*=0) corresponds to the repulsion term for molecules without volume satisfying the conservation of mass, momentum and energy in *p*–*s* collisions. A model for this term was discussed by Santos *et al.* [5]. Finally, *f*^{(p)}_{MB}(** u**) is the Maxwell–Boltzmann equilibrium distribution centred around

*n*

^{(p)},

**and**

*u**T*, 2.9the parameter

*D*representing the Euclidean dimension of physical space.

Using the same reasoning as in He & Doolen [15], considering that molecular chaos prevails for |** η**|>

*σ*

_{ps}, enables the attraction term

*Ω*

^{(ps)}

_{att}to be written as 2.10where 2.11is a mean potential that takes account of the

*p*–

*s*interaction with the

*s*molecules that are outside the

*σ*

_{ps}sphere centred on the

*p*molecule.

Using a Taylor expansion of *n*^{(s)}(** x**+

**) around**

*η***, the two first non-null terms in this expansion are 2.12where**

*x**a*

^{(ps)}and

*κ*

^{(ps)}are molecular parameters related to the potential energy

*Φ*

^{(ps)}.

With these assumptions, the full kinetic model for non-ideal mixtures can be written as
2.13This model retrieves the advection–diffusion equations and the momentum balance equations when the condition
2.14is satisfied, with the pressure tensor given as
2.15where is the surface tension tensor
2.16and *P*_{s} is the scalar pressure
2.17where *P* is the thermodynamic pressure and
2.18

The second-order terms in the scalar pressure and in the interfacial tension tensor, although not explicitly appearing in the Shan & Chen model [3,4], result from higher-order Taylor expansions in the calculation of the spatial derivatives. Therefore, segregation and interfacial tension are predicted to be related to effects with the same order of magnitude as the errors in the discrete approximation, when the isotropy of the lattice tensors used in the discrete calculation of the spatial derivatives is not assured up to a given rank. The same comments hold for the Santos *et al*. [5] model, with the particularity that their field mediators enable the local calculation of spatial derivatives. The presence of these second-order terms in the pressure tensor is consistent with the Helmholtz free energy approaches for fluid interfaces and the terms are due to the strong number-density gradients in these layers.

## 3. Thermodynamic consistency

The kinetic level of description, leading to discrete LB equations, may be thought as a two-way bridge linking the molecular and the macroscopic levels of description. Therefore, although derived from the molecular domain and, in consequence, featuring some properties of fluid interfaces that are inaccessible to the macroscopic level of description, the molecular parameters that appear in the kinetic model, equation (2.12), may be related to macroscopic measurable quantities, and the kinetic model itself can be rewritten in such a manner as to recover the correct macroscopic equations.

When is given by equation (2.12) and the volume parameter *χ* in equation (2.5) is written as 1/(1−*bn*), the model, equation (2.13), predicts a van der Waals equation of state for the thermodynamic pressure *P*. Nevertheless, adopting a reverse standpoint, equation (2.14) can be used to find the potential for any given cubic equation of state *P*=*P*(*n*,*T*). Indeed, the kinetic model, equation (2.13), was derived considering a number of assumptions: molecules were first considered as material points and, although Enskog’s volume exclusion theory led to the correct repulsion part of the thermodynamic pressure, molecular shapes were not taken into account and the molecular chaos hypothesis for the attraction term was used in the derivation. On the other hand, considering their great technological importance, the equilibrium properties of real, non-ideal liquids and their mixtures have been exhaustively studied in the last century, since van der Waals first attempted to understand the thermodynamic properties of vapours and phase transitions [20]. Equations of state and mixing rules have been heuristically proposed and successfully validated against pressure–temperature–volume experiments and their use in the LB method (LBM) was shown to be beneficial, improving the stability of numerical schemes [21].

The main drawback of the kinetic model, equation (2.13), is, nevertheless, its inability to correctly predict the balance equation for the thermodynamic internal energy. Since the molecules are modelled by particles without any other form of kinetic energy, except for that related to their translational degrees of freedom, the thermodynamic internal energy per unit volume of the mixture, *ρe*, is the sum of the peculiar kinetic energy of the molecules,
3.1and the volumetric intermolecular potential energy, . Therefore, the balance equation for *ρe* can be derived as a moment of equation (2.13), as
3.2where ** q** is the heat flow vector and is the viscous stress tensor. The result is that only the repulsive part of the thermodynamic pressure appears to be responsible for the compression work in this equation, and this lack of consistency is to be attributed to the assumptions mentioned above in deriving equation (2.13). As can be concluded from this equation, the mean-field theory leads to a kinetic model where the intermolecular forces, ∇

*Φ*

_{m}, mimic the effects of the external body forces,

*g*^{(p)}. In their 2002 paper for pure compounds, He & Doolen [15] found the compression work to be given by the full pressure tensor, considering the intermolecular potential energy equation

*with a spurious source term*that cannot be reduced to the divergence of a flux. This is wrong, since the source of the total energy in a unit volume must be zero, 3.3where

*S*(

*E*

_{c}) and

*S*(

*E*

_{g}) are, respectively, the sources of kinetic and body force potential energies. Accordingly, the source of the intermolecular potential energy

*Φ*in this volume must be equal to −[

*S*(

*E*

_{c})+

*S*(

*E*

_{g})] to give the correct balance for the total energy. In fact, in the absence of gravity the only sources of kinetic energy are the intermolecular forces among the particles, and this source term can be calculated, in the present case, to be 3.4

This lack of consistency may be corrected by rewriting the kinetic model, equation (2.13), as
3.5The second term on the right-hand side of the above equation was rewritten by replacing *f*^{(p)} by its *zeroth* approximation , since this form is more appropriate for LBEs and this change does not have any influence on the macroscopic equations. In the same manner, since is orthogonal to 1 and to *ξ*_{α}, the correction term, related to *F*, does not affect the advection–diffusion and momentum balance equations, giving the correct energy balance when *F* satisfies
3.6

The use of heuristic correction terms is very common in LBM practice. Nevertheless, these corrections are usually made on the discrete LB equation itself and lead to spurious terms in the transport equations that are difficult to control and, frequently, to a lack of Galilean invariance [11]. Indeed, there are several sources of errors in deriving LBEs. As in all numerical schemes, some errors are only dependent on the lattice spacing Δ and become meaningless when . However, the most serious sources of errors in the LBM framework are due to the velocity discretization, usually leading to spurious terms in the macroscopic equations that are dependent on the local velocity, ** u**. High-accuracy models such as the one shown in equation (3.5) demand a velocity discretization with an improved set of discrete velocities. In fact, in its early development, the LBM was considered to be applicable only to flows with small Mach number and constant temperature. A key breakthrough in this area was the identification of the correspondence between the LBE and the kinetic continuous equation, and the invention of high-order LB schemes [16–18], leading to velocity sets that, when used in a discrete velocity kinetic scheme, ensure accurate recovery of the high-order hydrodynamic moments.

## 4. Conclusion

In this paper, the thermodynamic consistency in deriving kinetic models from the Liouville equation was investigated in relation to the study of segregation in non-ideal mixtures in non-isothermal conditions. The kinetic model, equation (3.5), satisfies all the requirements imposed by the macroscopic level of description for describing non-ideal mixtures and, simultaneously, is accurate enough for describing the molecular-induced physical phenomena that happen inside fluid interfaces. The presence in the model of first-, second- and even third-order spatial derivatives requires the use of improved neighbouring and, consequently, the use of increased velocity sets, in the calculation of these derivatives. Assuring increasingly higher order of isotropy of the lattice tensors is a necessary condition for taming spurious currents [14]. The discretization conditions imposed by the kinetic model, equation (3.5), require the norm preservation of sixth-order moments and will be discussed in a forthcoming paper.

## Acknowledgements

The authors are greatly indebted to CNPq Brazilian Research Council, Finep Brazilian Agency for Research and Projects, and Petrobras Brazilian Petroleum Company for financial support. The authors thank Xiaowen Shan, from Exa Corporation, for very fruitful discussions.

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society