## Abstract

Recently, we proposed a theoretical framework to include thermal fluctuations into the Lattice Boltzmann (LB) method for non-ideal fluids. Here, we apply a variant thereof to a certain class of force-based non-ideal fluid LB models. We find that ideal-gas-like noise is an exact result of the fluctuation–dissipation theorem in the hydrodynamic regime. It is shown that satisfactory equilibration of the density and fluid momentum can be obtained in a simulation over a wide range of length scales.

## 1. Introduction

The Lattice Boltzmann (LB) method is a deterministic approach to the numerical solution of the hydrodynamical equations, hence, *a priori* it does not take into account the effect of the thermal motion of the fluid particles, which gives rise to random fluctuations of the fluid density and momentum. If such effects are to be simulated within the LB method, the algorithm must be extended to include explicit random forces with magnitudes being determined from an appropriate fluctuation–dissipation theorem (FDT). The theory of thermal fluctuations in the ideal-gas model is now well established [1–3], and has been applied to various problems in soft matter physics [4]. Until now, the implementation of thermal fluctuations into *non-ideal* LB fluids remained unsolved, although there exist many potential applications of such a simulation algorithm in nano- and microfluidics of free surface flows, such as jet instabilities, droplet spreading or dewetting processes [5,6].

Recently, we have shown how the Langevin framework originally devised by Adhikari *et al.* [2] for the ideal LB gas, can be extended to describe thermal fluctuations in a non-ideal fluid [7]. In this approach, a Fourier-space version of the FDT is derived in a straightforward manner directly from the discrete time dynamics of the LB equation (LBE). Application of the FDT to the non-ideal fluid model of Swift *et al.* [8] led to a spatially correlated form of noise that was shown to ensure complete equilibration of a homogeneous fluid at all length scales. In situations where strong inhomogeneities, such as coexisting liquid and vapour phases, play a role, the use of spatially correlated noise becomes technically cumbersome due to lost translational invariance. In these cases, it is preferable to use a spatially uncorrelated form of noise, which can be applied to inhomogeneous systems by taking into account the local values of the density and transport coefficients. In case of the Swift model, the wavenumber-dependent noise terms are found to be proportional to the ‘square-gradient’ parameter [7], and hence, by reducing this parameter sufficiently, the offending terms become negligible and spatially uncorrelated noise results. Alternatively, uncorrelated noise can be derived from the FDT in the zero wavenumber limit. This is the focus of the present work.

A natural question that arises is, whether and to what extent noise derived in the zero wavenumber limit can thermalize an LB fluid with sufficient accuracy also at finite wavenumbers. Physically, one might expect this to be true at least in the hydrodynamic regime (low wavenumbers), as the non-ideal-gas interactions are usually reversible in the hydrodynamic limit, i.e. they do not (or, if at all, only weakly) contribute to the dissipative terms in the Navier–Stokes equations. A well-known analogous example is the linear Langevin equation of a Brownian particle coupled to a spring [9], in which case, the statistical properties of the random-noise term are completely independent of the properties of the Hookean force. Indeed, as has been shown in Gross *et al.* [7], in the case of the Swift model, in a certain parameter window, spatially uncorrelated noise ensures correct equilibration of all degrees of freedom in an LB simulation, even up to the highest wavenumbers.

Based on these findings, we expect similar behaviour to also apply when noise derived in the zero wavenumber limit is employed in force-based non-ideal-fluid models. Therefore, we describe in the following how the appropriate FDTs can be constructed both at the level of the hydrodynamic equations and, independently, at the level of the LBE. By focusing on the zero wavenumber limit, technical complications associated with the derivation of the FDT from the LBE of force-based models at finite wavenumbers are avoided. While the presentation below is kept as general as possible, we apply the theory in the following to the model of Shan & Chen [10,11], which is a widely used lattice formulation of an effective mean-field theory. Results indicate that satisfactory thermalization of the model, even at the smallest scales, is indeed possible with the noise obtained by the present approach.

## 2. Theory

After introducing the Langevin extension of the LBE, taking as a particular example the Shan–Chen model, we proceed to analyse thermal fluctuations in the associated stochastic hydrodynamic equations and derive an expression for the density correlation function. Finally, an FDT in the zero wavenumber limit is derived from the LBE. The results obtained here hold in principle for any non-ideal-fluid LB model that employs an interaction force **F**^{int} which can be written as the divergence of a tensor, **F**^{int}=−∇⋅*P* or, more simply, the gradient of a scalar. Besides the Shan–Chen model, this is fulfilled in particular for force-based models of the He–Shan–Doolen type [12]. Moreover, as will become clear below, the details of how this interaction force is incorporated into the model (i.e. via an explicit forcing term or by redefining the macroscopic velocity) is not relevant concerning the FDT.

### (a) Model

We consider the general case of a D*d*Q*n* lattice, i.e. a *d*-dimensional grid with *n* different velocity vectors **c**_{i} (*i*=1,…,*n*) associated with every grid point. The distribution function in the Langevin extension of the Shan–Chen model evolves according to the fluctuating LBE
2.1
where is the fluid density, is the macroscopic flow velocity, is an equilibrium distribution (second-order polynomial in the fluid velocity) and *τ* is a relaxation time. The *ξ*_{i} are a set of local Gaussian random variables of zero mean and a covariance that will be specified below by invoking an FDT.

In the Shan–Chen model, the thermodynamic interactions are mediated through a body force **F**^{int} that is used in the evaluation of the equilibrium distribution in the LBE (2.1). Taylor expanding the interaction force up to the third order shows that **F**^{int} can be written as [10,13]
where is the ideal-gas pressure tensor and *σ*^{2}_{s}=1/3 is a numerical constant (equal to the squared speed of sound of an ideal-gas model). The non-ideal (Shan–Chen) pressure tensor reads
2.2
Here, is an interaction constant and *ψ* is a mean-field potential taken as .

### (b) Hydrodynamics

As can be shown through a Chapman–Enskog analysis [7,10], at long time and length scales, the above model approximates the full nonlinear stochastic Navier–Stokes equations,
2.3
where is the kinematic shear viscosity and *ν*_{b}=*ν*_{s} is the kinematic bulk viscosity. *R*_{αβ} is a random stress tensor that is directly related to the random variables *ξ*_{i} in the fluctuating LBE (2.1) [7].

The statistical properties of the *R*_{αβ} can, in principle, be determined independently of the *ξ*_{i} by deriving an FDT from the hydrodynamic equations [1,4,7]. To this end, we linearize equations (2.3) around a quiescent state, writing *ρ*=*ρ*_{0}+*δρ*, **u**=*δ***u**. Computing the pressure tensor to linear order in the density fluctuation gives
2.4
where *ψ*_{0}≡*ψ*(*ρ*_{0}). We identify the thermodynamic speed of sound of the model as . After Fourier transforming the linearized hydrodynamic equations, it is convenient to separate components longitudinal and transverse to the wavevector **k**, using , where , , , and analogously , where , . We finally obtain
2.5
2.6
and
2.7
where *ν*_{l}=[*ζ*+*η*(2−2/*d*)]/*ρ*_{0} denotes the longitudinal and *ν*_{t}=*η*/*ρ*_{0} the transverse kinematic viscosity. We see from equation (2.6) for the longitudinal momentum, that the only effect of the interaction force is to give rise to a wavenumber-dependent speed of sound,
2.8

For an equilibrium fluid, the static correlation function of the momentum density must be independent of wavenumber and is given by Landau & Lifshitz [14]
2.9
In the above, *T* is the temperature and Δ*V* is a volume element that arises by dimensional consistency and is henceforth set equal to 1 in lattice units. As the dissipative terms in the Navier–Stokes equations (2.6) are not affected by the interactions, it immediately follows that the covariance of the random stress tensor is given by Gross *et al.* [7],
2.10
which, in real space, translates to the well-known Landau–Lifshitz expression [14]
2.11

Despite the fact that the Shan–Chen model cannot be obtained from a free-energy functional, it is nevertheless possible to derive an expression for the structure factor in the hydrodynamic regime with the help of the FDT (2.10). For this task, we insert the continuity equation (2.5) into the equation for the longitudinal momentum (2.6), thereby obtaining a single equation for the density fluctuations,
Fourier transforming in time, we find a linear-response relation between the random force and the density fluctuation,
Hence, we get a relation between the dynamic structure factor and the longitudinal random stress correlations,
from which the desired static structure factor follows as
The integral can be computed using contour integration, resulting in
2.12
Comparing to the structure factor derived from a square-gradient free-energy functional [7,9], we see that the quantity plays the role of an effective ‘square-gradient’ parameter. We remark that, if a well-defined thermodynamic framework for the structure factor is needed, it might be advantageous to invoke a modified formulation of the Shan–Chen model based on a pseudo-free-energy functional, as has recently been proposed in Sbragaglia *et al.* [15].

### (c) Lattice Boltzmann fluctuation–dissipation theorem

Based on our recently introduced framework [7], we will describe in the following how the covariance matrix of all the *ξ*_{i} in the hydrodynamic regime can be directly obtained from the LBE. The theoretical development of the FDT at the level of the fluctuating LBE proceeds most easily in the space of moments of the distribution function *f*_{i}. For this task, a set of basis vectors *T*_{ai} orthogonal under a weighted inner product, *w*_{i}*T*_{ai}*T*_{bi}=*N*_{a}*δ*_{ab}, are constructed from the lattice velocities **c**_{i}. Here, *N*_{a} is the length of the *a*th basis vector *T*_{a}, and *w*_{i} are a set of weights that appear in the global-Maxwellian form of the equilibrium distribution function, *f*_{i}(*ρ*,**u**=0)=*ρw*_{i}. Several choices existing in the literature for the *T*_{ai} are possible [3,7,16]. The moments *m*_{a} (*a*=1,…,*n*) of the distribution function are now defined by projecting *f*_{i} onto the basis vectors,
For the present purpose, we will only assume that *m*_{1}=*ρ*, *m*_{2,…,d+1}=*ρu*_{α}, while the next *d*(*d*+1)/2 modes correspond to the stresses, and the last *n*−(*d*+1)−*d*(*d*+1)/2 modes are the so-called ghost or kinetic modes.

As equation (2.11) shows, the random stress tensor obtained for the present non-ideal fluid model is independent of the wavenumber and of the same form as for the ideal-gas model. Since the random stress tensor is directly connected to the subset of noise variables that pertain to the stress modes [7], one can expect the LB noise also to be independent of *k* in the hydrodynamic regime. This statement relies on the observation that the influence of the ghost modes on the LB dynamics is generally weak in the hydrodynamic regime. In fact, as a Chapman–Enskog analysis directly shows, ghost modes are decoupled from the hydrodynamic equations up to the second order of the multi-time scale expansion. Hence, a possible influence of the non-ideality on the ghost noise can become significant—if at all—only at larger wavenumbers. This crucial insight allows us in the following to restrict the derivation of the FDT for the fluctuating LBE to the zero wavenumber limit.

To derive the FDT, the fluctuating LBE (2.1) is linearized around a quiescent reference state described by ,
2.13
where , and is the linearized interaction force. In a translationally invariant system, the advection operator becomes equal to unity in the limit of zero wavenumber, as . Furthermore, we have , as the interaction force can be expressed fully as a derivative. Of course, this is tantamount to stating that the interaction force conserves momentum globally [11]. After projecting the linearized LBE onto the basis set *T*_{a}, we finally obtain the fluctuating LBE in the zero wavenumber limit (indicated by the superscript ‘0’),
2.14
where and we introduced a linear matrix collision operator *Ω*=diag(0,…,0,−1/*τ*,…), with the first *d*+1 entries being zero owing to mass and momentum conservation. We note that, in the zero wavenumber limit, the linearized fluctuating LBE of the Shan–Chen model reduces to the one of the ideal gas.

The FDT follows now by multiplying equation (2.14) by from the right, averaging over the noise distribution and assuming stationarity of equal-time correlators, obtaining
where we defined the equilibrium correlation matrix of the modes . This matrix is not immediately provided by the LBE, but must instead be derived from a statistical mechanical framework. A basic relation for correlations in the distribution function is provided by kinetic theory [7],
2.15
Here, *S*^{0} is the zero wavenumber limit of the structure factor (2.12), given by . The quantity *μ* is a mass parameter whose value has to be determined by requiring self-consistency of equation (2.15) with the the structure factor and the momentum correlation function (2.9), when expressed in terms of the LB quantities. In fact, computing the correlation of the LB momentum, using with equation (2.15), and comparing with equation (2.9) shows that one must set . The desired noise covariance matrix in the hydrodynamic limit is then finally obtained as
2.16
We note that this noise is identical to the one used in the ideal-gas case [3] and that, furthermore, conservation of mass and momentum is an ‘automatic’ result of the FDT.

## 3. Simulations

Results of simulations of the fluctuating LBE (2.1) using Gaussian noise with covariance defined by equation (2.16) are presented below. For brevity, we focus here only on the density and momentum fluctuations. Note that, as expressed by equation (2.16), thermal noise is put only onto the stress and ghost modes, while density and momentum fluctuations are a result of the LBE-induced coupling between the conserved and non-conserved modes. Simulations are performed in a uniform one-phase system of density *ρ*_{0}=*ρ*_{L}, using a D3Q19 implementation of the Shan–Chen model, with the interaction constant fixed to , giving rise to equilibrium liquid and vapour densities of *ρ*_{V}=0.134, *ρ*_{L}=1.82. The relaxation time is set to *τ*=1.0.

In figure 1*a*, the equal-time density correlation function obtained from simulation is compared to the theoretical structure factor, equation (2.12). We find reasonable agreement for low and intermediate wavenumbers, with a relative error generally smaller than 10 per cent. Deviations at higher wavenumbers are attributed partly to the fact that the linearized pressure tensor (2.4) used to obtain the theoretical structure factor is based on a Taylor expansion of the interaction force, while the interaction force itself is computed in the algorithm without the use of derivatives [11].

As figure 1*b* shows, the fluid momentum is almost perfectly thermalized at low and intermediate wavenumbers (error below 5%). Stronger deviations are only found in the high-*k* region of the longitudinal momentum (see inset of figure 1*b*), which point to a possible lack of reversibility of the forcing term at small scales. These results are found to also remain qualitatively similar for other values of the interaction constant and relaxation time.

A detailed investigation of the dependence of the equilibration errors on the system size and correlation length (interface width), as well as the origin of the apparent FDT violation at small scales are left for future research.

## Acknowledgements

M.G. thanks G. Zikos for useful discussions. M.E.C. acknowledges support from the Royal Society and funding from EPSRC EP/E030173. M.G. and F.V. further acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) under grant no. Va205/3-3 (within the Priority Programme SPP1164), as well as funding from the industrial and public sponsors of ICAMS.

## Footnotes

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

- This journal is © 2011 The Royal Society