## Abstract

The lattice Boltzmann equation was introduced about 20 years ago as a new paradigm for computational fluid dynamics. In this paper, we revisit the main formulation of the lattice Boltzmann collision integral (matrix model) and introduce a new two-parametric family of collision operators, which permits us to combine enhanced stability and accuracy of matrix models with the outstanding simplicity of the most popular single-relaxation time schemes. The option of the revised lattice Boltzmann equation is demonstrated through numerical simulations of a three-dimensional lid-driven cavity.

## 1. Introduction

Twenty years ago, the first computationally competitive lattice Boltzmann (LB) scheme was introduced in a form whereby interparticle collisions are represented through a scattering matrix between the various discrete-velocity Boltzmann distributions living on the same node of the spatial lattice [1]. By prescribing the spectral structure of the scattering matrix based on the desired macroscopic target equations (Navier–Stokes) rather than on an underlying microscopic dynamics, this top-down formulation enabled values of the Reynolds number to be achieved as high as permitted by the grid resolution, thereby making the first competitive contact with computational fluid dynamics at large (e.g. the earliest review [2], as well as the most recent review [3] covering LB applications in many fields of fluid dynamics, from turbulence to micro-flows, and references therein). The original matrix formulation was quickly superseded by the so-called lattice Bhatnagar–Gross–Krook (LBGK) equation, in which the scattering matrix is reduced to a diagonal form, with a single parameter controlling the relaxation rate of fluid momentum (viscosity), as well as any other non-conserved kinetic moment [4,5]. Despite its obvious limitations, namely all kinetic modes relaxing at the same pace to local equilibrium, the LBGK has rapidly assumed a dominant role in the field, mostly on account of its simplicity. At around the same time that the LBGK was developed, the original matrix version was revisited, optimized and turned into what has now come to be known as the multiple relaxation time (MRT) LB [6–8], recently improved with Galilean invariance by the so-called cascaded formulation [9]. Although the MRT provides a useful outgrowth of the LBGK, the mainstream of LB research remains with the LBGK because of its unsurpassed simplicity.

As a result, it appears that an optimum compromise between the enhanced stability and flexibility of the MRT and the computational handiness of the LBGK is much needed. In this work, we shall present precisely such an optimum compromise. We introduce a new scheme (the revised matrix LB equation; RM hereafter) which aims to combine the best of the two options, namely enhanced stability with virtually the same simplicity as the LBGK. Our option is based on a *two-step* relaxation BGK-like collision operator, in which the two time scales describe relaxation to two *distinct* equilibria, one of which is the (usual) fixed point of the relaxation and the other is an intermediate (quasi-equilibrium) state. Models of this type are well known in kinetic theory (e.g. [10] and references therein) and are used to achieve a description of fluids with realistic Prandtl or Schmidt numbers not possible with the BGK. At variance with previous formulations, the present RM method is aimed at boosting the stability and accuracy of athermal LB models in the incompressible limit. In this paper, we provide an example of the RM construction on the commonly used 19-velocity lattice, which retains much of the factorization property. On the practical side, the RM numerical algorithm is pretty simple: all it takes is just a few line changes in existing LBGK codes. Supporting numerical evidence is provided through the simulation of a lid-driven cavity, which shows the viability of the revised matrix model.

## 2. Lattice Boltzmann method

### (a) Matrix lattice Boltzmann model

Let *f*(** x**,

*t*) be the vector of populations

*f*

_{i}corresponding to the

*D*-dimensional discrete velocities

*v*_{i},

*i*=1,…,

*b*, at the site

**at time**

*x**t*. Throughout the paper, we denote a sum over the discrete-velocity index. The matrix LB equation for the incompressible flow simulation has the form [1] 2.1 where

*f*

^{eq}is a local equilibrium which depends on

*ρ*=〈

*f*〉 (local density) and

*ρ*

**=〈**

*u*

*v**f*〉 (local momentum with

**the local flow velocity). The spectral properties and computer implementation of the scattering matrix**

*u**A*are extensively discussed in [6,7]. If

*A*=−

*ω*

_{1}Id, where Id is the identity matrix, then the matrix LB (2.1) reduces to a one-parametric LBGK model [4,5], 2.2 With the equilibrium

*f*

^{eq}satisfying the isotropy conditions for the pressure tensor

*ρΠ*

_{αβ}=〈

*v*

_{α}

*v*

_{β}

*f*〉 and the third-order moment

*Q*

_{αβγ}=〈

*v*

_{α}

*v*

_{β}

*v*

_{γ}

*f*〉, 2.3 and 2.4 where

*c*

_{s}is the lattice speed of sound, the LBGK model (2.2) remains the most popular version to date of the LB for the simulation of incompressible flow, mostly because of its outstanding simplicity.

### (b) Matrix lattice Boltzmann reloaded model

A characteristic property of the matrix LB equation (2.1) is that only one state is specified (the local equilibrium *f*^{eq}). We here suggest a different two-parametric family of LB equations based on an intermediate (quasi-equilibrium) state of relaxation *f*^{C},
2.5
Let us introduce a local temperature variable *Θ* (*D* is the space dimension),
2.6
With this, we consider an intermediate state *f*^{C} which satisfies the following properties:
2.7
2.8
2.9
Condition (2.7) maintains the local density and momentum conservation, while condition (2.8) guarantees that the zero of the collision integral is at the local equilibrium. Finally, condition (2.9) regulates the pressure tensor *ρΠ*^{C}_{αβ} at the intermediate state of relaxation. Two limiting cases are most interesting: at *λ*=0, this tensor is kept at the equilibrium, , while the choice *Λ*=*θ* (*λ*=1) corresponds to the familiar ideal gas equation of state. When the two relaxation parameters are equal (*ω*_{1}=*ω*_{2}), equation (2.5) reduces to the LBGK (2.2). On the other hand, equation (2.5) can be written in the LBGK-like form,
2.10
where a generalized equilibrium *f*^{GE},
2.11
is a convex (because of the realizability condition, *ω*_{2}/*ω*_{1}≤1) combination between *f*^{C} and *f*^{eq}. As we shall see below, the proper choice of *f*^{C} endows the LBGK-like relaxation process (2.10) with some valuable properties not accessible to the single-time LBGK scheme (2.2).

With an intermediate state *f*^{C} satisfying the conditions (2.7)–(2.9) the two-parametric family of *revised matrix* (RM) equations (2.5) recovers the athermal Navier–Stokes equations with the kinematic (shear) viscosity *ν* and the second (bulk) viscosity *ξ* as follows:
2.12
and
2.13
Shear viscosity *ν* (2.12) is the only transport coefficient which remains pertinent in the low-Mach number (incompressible flow) limit. In that respect, the two-parametric RM family defines an equivalence class of LB models for the incompressible flow simulation, with the LBGK being a particular element of this family. On the contrary, the second (bulk) viscosity coefficient depends on the choice of *f*^{C}. Note that, as it is pertinent to the kinetics without the energy conservation, the slope of the pressure *ρΛ* with respect to the ‘temperature’ *Θ* in equation (2.9) defines the ratio between the bulk and the shear viscosity. In the two limiting cases mentioned above, we have, at *λ*=0, the shear and bulk viscosities are equal, whereas, at *λ*=1, the bulk viscosity is decoupled from the shear, and is defined by the second relaxation parameter *ω*_{2} alone (and thus *ξ* can be regarded as a free tunable parameter to enhance numerical stability [11,12]). In certain cases, when the intermediate state *f*^{C} can be described as the minimum of the entropy function, one can prove that the bulk viscosity in the present model should be larger than or equal to the shear (e.g. [11]). Note that conditions (2.7)–(2.9) do not yet define the intermediate state *f*^{C} uniquely, and the actual performance of the scheme may be strongly affected by a particular choice of *f*^{C}. For example, even for *λ*=0, the RM is not necessarily equivalent to the LBGK, even though they share the same shear and bulk viscosities.

### (c) Unidirectional quasi-equilibrium on the Maxwell lattice

The essential ingredient of the RM model (2.5) is the intermediate state of relaxation *f*^{C}, which needs to be explicitly designed so as to accomplish the aforementioned optimization. This task is facilitated by considering a specific family of populations, providing a universal template for constructing various values of *f*^{C} with tailored properties [11,13].

In the dimension *D*=2,3, the fundamental Maxwell lattices are generated by a tensor product of *D* copies of the one-dimensional velocities {−1,0,1} (the D2Q9 and D3Q27 lattices), which we consider first. For *D*=3, the discrete velocities are denoted as *v*=(*i*,*j*,*k*), *i*,*j*,*k*∈{−1,0,1} in a fixed Cartesian coordinate system. A weight *W*_{(i,j,k)}=*W*_{(i)}*W*_{(j)}*W*_{(k)} corresponds to each velocity where *W*_{(0)}=2/3 and *W*_{(−1)}=*W*_{(1)}=1/6 are one-dimensional weights, while the speed of sound is independent of the dimension. With this, the entropy function of the D3Q27 model is defined as [14].

Following Asinari & Karlin [11,15] and Karlin & Asinari [13], for D3Q27, we define a special *unidirectional quasi-equilibrium* population (UniQuE) *f** as a minimizer of the entropy function *H* under fixed density, momentum, and three diagonal components of the pressure tensor,
2.14
Minimization is achieved by a fully factorized UniQuE population,
2.15
where
2.16
We further refer to the three-tuple (2.16) as a *unidirectional particle* (in a direction *α*). The UniQuE is the main source for derivations of the intermediate states of relaxation *f*^{C}, with tailored properties for the RM.

### (d) Unidirectional quasi-equilibrium on sub-Maxwell lattices

The two sub-lattices of the D3Q27 which also satisfy the isotropy relation (2.4) are known as the D3Q19 and D3Q15 lattices. In this case, the analogue of the UniQuE (2.15) is readily constructed by a projection pruning algorithm [13]. Let us introduce a unidirectional *anti*-particle that is a three-tuple,
2.17
The (analogue of) UniQuE on D3Q19 is then represented as the difference between the factorized distributions of unidirectional particles and those of the anti-particles:
2.18

## 3. Example of the revised matrix lattice Boltzmann model on the D3Q19 lattice

Recent studies have indicated that kinetic modes stemming from the dynamics of the third-order moments can severely affect the stability of the LBGK [7,9,16]. Decoupling of the relaxation of the third-order moments *ρQ*_{αββ} from the relaxation of the rest of the moments represents, therefore, our objective here in constructing the intermediate state *f*^{C} for the RM. In order to construct the corresponding intermediate state *f*^{C} from the UniQuE of the D3Q19 lattice (2.18), the following two steps are taken. (i) Unbias the third-order UniQuE moments : replace *u*_{α}*Π*_{ββ}→*Q*_{αββ} in equation (2.18), where *ρQ*_{αββ} is the corresponding third-order moment evaluated at the current population *f*. (ii) Equilibrate the rest of the moments: replace . Furthermore, the local equilibrium *f*^{eq} is constructed by equilibration of all the moments *Π*_{αα}→*Π*^{eq}_{αα} in *f** (2.18). The intermediate state *f*^{C} and the equilibrium *f*^{eq} can be combined into the generalized equilibrium *f*^{GE} (2.11),
3.1
where *σ*,*μ*,*δ*∈{−1,1}, and
3.2
Substituting the generalized equilibrium (3.1) and (3.2) into the LBGK-like form (2.10), we set up the RM model on the D3Q19 lattice.

We shall now proceed with a numerical validation of the present scheme. For that, we have chosen the diagonally lid-driven three-dimensional cavity flow, the standard test for LB models suggested in d’Humières *et al*. [7]. The cavity is a cubic box with a unit edge and it fits into a Cartesian coordinate system, namely (*x*,*y*,*z*). The boundary condition at the top plane (*x*,*y*,1) is so that *u*_{L}=∥*u*_{L}∥=1/10. The other five planes are subject to no-slip boundary conditions. The kinematic (shear) viscosity is equal to *ν*=3/1000. The computational domain is discretized by a uniform collocated grid with *N*^{3} points with *N*=60. The boundaries are located a half-cell away from the computational nodes. Hence, the Reynolds number is Re=2000, and consequently the relaxation parameter controlling the kinematic viscosity is equal to *ω*_{1}=1000/509≈1.9646. Let us denote *x*_{w} the generic boundary computational node. In all inner computational nodes (** x**≠

*x*_{w}), equation (2.5) holds for any lattice velocity

*v*_{i}. Following d’Humières

*et al*. [7], in the generic boundary computational node

*x*_{w}, the following condition holds:

*f*(

*x*_{w},

*t*+

*δt*)=

*f*

^{eq}(

*ρ*

_{w},

*u*_{w}), where

*ρ*

_{w}=(1−

*η*

_{w})

*ρ*(

*x*_{w}+

*n*_{w},

*t*)+

*η*

_{w}

*ρ*(

*x*_{w},0), , where

*η*

_{w}is a tunable parameter,

*n*_{w}is the vector normal to the wall and pointing towards the fluid (e.g.

*n*_{w}=(0,0,−1) for the sliding wall) and

*u*_{d}is the desired velocity at the wall (located a half-cell away from the computational node). In particular,

*η*

_{w}=0 and

*u*_{d}=

**0**for all the walls, with the exception of the sliding wall for which

*η*

_{w}=1 and . Finally, the number of time steps (roughly 18 000) has been selected in order to reach

*E*(

*t*+1)/

*E*(

*t*)−1≤5×10

^{−5}, where

*E*(

*t*) is the total kinetic energy in the cavity. In summary, the boundary conditions and the Reynolds number are taken from the MRT set-up studied in d’Humières

*et al*. [7] in order to make a comparison of the performance of various LB models.

Three different models are considered for the bulk equation, namely (i) the LBGK model: equation (2.1) with *A*=−*ω*_{1}Id; (ii) the MRT model: equation (2.1) with the scattering matrix *A* as given in d’Humières *et al*. [7], appendix A; and (iii) the RM model: equation (2.10) with *ω*_{2}=1.2<*ω*_{1} and *f*^{GE} given by equation (3.1). Note that the present choice of the second relaxation parameter *ω*_{2}=1.2 is made in order to match the corresponding rate of the third-order moment relaxation of the MRT model of d’Humières *et al*. [7].

For the present RM model, the velocity vectors of the diagonally driven cavity flow for Re=2000 at *z*=0.5 (mid-plane parallel to the sliding plane) and the pressure field are reported in figure 1.

It has already been reported in d’Humières *et al*. [7] that the pressure field of the LBGK model in the present set-up is severely contaminated by the chessboard mode. Our simulation confirmed this observation. In figure 2, we report the velocity and the pressure fields as obtained by the MRT model of d’Humières *et al*. [7] (the realization of the MRT model as suggested in d’Humières *et al*. [7] was followed). Results presented in figure 2 are consistent with those of d’Humières *et al*. [7]. Note that, while the velocity fields are practically identical for the MRT and the present RM models, they both significantly improve the pressure, when compared with the LBGK case. Moreover, the present RM model does so essentially at a computational cost of the standard LBGK.

Summarizing, the advantage of this revised matrix LB option is to considerably enlarge the stability domain of the most popular LBGK scheme, while retaining its simplicity and computational efficiency. Three-dimensional simulations on the D3Q19 lattice indicate that the RM scheme enhances the LBGK model, basically at the same computational cost. Remarkably, the implementation of the RM is very straightforward, requiring, as it does, just the change of a few lines in existing LBGK codes. Based on the above, it is hoped that the RM can be brought to broad fruition in LB research at large.

## Acknowledgements

P.A. acknowledges support of EnerGRID project. S.S. wishes to acknowledge ETHZ for financial support.

## Footnotes

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

- This journal is © 2011 The Royal Society