Intro to the Lattice Boltzmann Equation


The Lattice Boltzmann Equation (LBE) is a mesoscale modelling method with a statistical mechanics approach. Instead of tracking the motion of individual particles, distribution functions can be defined to describe the probability of finding particles at a given point in time and space with a particular momentum. The distribution functions evolve over time by considering how the particles collide and subsequently move, while summing their moments over all possible particle momenta gives macroscopic properties such as fluid density and momentum. Confining the particles to a grid and to preset links between grid points simplifies the calculations but still provides enough information to correctly calculate hydrodynamic behaviour. Macroscopic properties such as fluid density and momentum can be obtained by summing up moments of the distribution functions over all available particle momenta.

The particles in LBE can be considered flexibly, thus allowing larger length and time scales than atomistic modelling methods while still incorporating some detail of molecular interactions. Boundary conditions can be treated in intuitive and simple ways, allowing systems with complex geometries to be modelled nearly as efficiently as simpler ones. Methods to allow multiple fluids and phases to be modelled in LBE without fundamentally changing the algorithm have also been devised, as have ways to represent fluids with non-standard rheological behaviour and to incorporate diffusion and heat transfer effects.


LBE can be considered as a variant form of lattice gas cellular automata (LGCA), a modelling method that directly models particles colliding and moving around a regular grid. By substituting binary particle occupation numbers with probability distribution functions [Frisch1987], several of LGCA’s shortcomings as a method to model fluid flows (e.g. noise, lack of Galilean invariance) are overcome.

Given a lattice link \(i\) with vector \(\mathbf{e}_i\), we can define an associated distribution function \(f_i\) for each grid point. The distribution functions evolve via the following equation:

(1)\[f_i \left(\mathbf{x}+\mathbf{e}_i \Delta t, t + \Delta t \right) - f_i \left(\mathbf{x}, t \right) = C_i,\]

where \(C_i\) is a collision operator for link \(i\) that may depend upon all distribution functions for a given grid point \(\mathbf{x}\). Equation (1) can be also expressed by separating out the collision and propagation stages as equations (2) and (3) respectively:

(2)\[f_i \left(\mathbf{x}, t^{+}\right) = f_i \left(\mathbf{x}, t\right) + C_i,\]
(3)\[f_i \left(\mathbf{x}+\mathbf{e}_i \Delta t, t + \Delta t \right) = f_i \left(\mathbf{x}, t^{+}\right).\]

The collision operator \(C_i\) can take several different forms, but the simplest and most common form (without considering forces acting on the fluid) is based on the Bhatnagar-Gross-Krook approximation that uses a relaxation time \(\tau\) as the parameter defining the timescale of collisions [Qian1992]:

(4)\[C_i = -\frac{f_i \left(\mathbf{x}, t\right) - f_i^{eq} \left(\mathbf{x}, t\right)}{\tau},\]

where \(f_i^{eq}\) is a distribution function corresponding to a local equilibrium state. (An interpretation of equation (4) is the collision pushes the system towards an equilibrium state, with \(\tau\) determining how quickly this occurs.) This local equilibrium distribution function can be determined as a function of macroscopic fluid density and velocity, usually of the form for mildly compressible fluids at speeds significantly lower than the speed of sound \(c_s\):

(5)\[f_i^{eq} \left(\rho, \mathbf{u}\right) = \rho w_i \left[1 + 3 \left(\mathbf{e}_i \cdot \mathbf{u} \right) + \frac{9}{2} \left(\mathbf{e}_i \cdot \mathbf{u} \right)^2 - \frac{3}{2} u^2 \right],\]

where \(w_i\) is a weighting parameter dependent on link and lattice scheme. To use equation (5) in the collision operator given in equation (4), we can calculate the macroscopic density and velocity by finding the zeroth and first moments of the distribution functions:

(6)\[\rho \left(\mathbf{x}, t\right) = \sum_i f_i \left(\mathbf{x}, t\right),\]
(7)\[\mathbf{u} \left(\mathbf{x}, t\right) = \frac{1}{\rho \left(\mathbf{x}, t\right)} \sum_i \mathbf{e}_i f_i \left(\mathbf{x}, t\right).\]

It is possible to show that the above equations can accurately represent fluid flow by assuming distribution functions are predominately local equilibrium values and their non-equilibrium parts scale with the Knudsen number (ratio of molecular mean free path to molecular length scale). This allows equation (1) to be expanded about the Knudsen number in time and space and the application of a Chapman-Enskog expansion to various time and length scales. By subsequently separating out equations based on different orders of Knudsen number and summing over all lattice links, conservation equations for mass and momentum can be obtained, leading to the Navier-Stokes equation for systems with small variations in density [Chen1998].

The above procedure allows us to define the speed of sound for the LBE fluid:

\[c_s = \frac{1}{\sqrt{3}} \frac{\Delta x}{\Delta t},\]

the equation of state:

\[p = \rho c_s^2,\]

and a relationship between kinematic viscosity and the relaxation time:

\[\nu = \frac{\mu}{\rho} = \frac{1}{3}\left(\tau - \frac{1}{2}\right) \frac{\Delta x^2}{\Delta t}.\]

The value of \(\tau\) can be selected to avoid numerical instabilities in LBE simulations: along with the kinematic viscosity and speed of sound for a given fluid, this sets the length scale (lattice spacing \(\Delta x\)) and time scale (timestep \(\Delta t\)).

On top of the above, it is possible to include additional fluids by either modelling them in separate lattices and calculating interaction forces between them, which can then be applied as an additional part of the collision operator, or by modifying equation (5) to obtain the correct free energy density function and apply interfacial tensions.

Alternative forms of equation (5) can be used to represent completely incompressible fluids and diffusive processes: the latter can represent diffusion of solutes or heat transfers and can be coupled to bulk fluid flows.

Boundary conditions can be defined to specify missing distribution functions ‘entering’ the system and satisfy required fluid densities and/or velocities. The simplest form is bounce back, where a boundary grid point reflects distribution functions entering it: this has the effect of applying a no-slip condition at any arbitrary boundary point without any extensive calculations, making it straightforward and computationally inexpensive to represent flows in porous media or past complex shapes.

[Frisch1987]U Frisch, B Hasslacher, P Lallemand, Y Pomeau and JP Rivet, Lattice gas hydrodynamics in two and three dimensions, Complex Systems, 1, 649-707, 1987.
[Qian1992]YH Qian, D d’Humières and P Lallemand, Lattice BGK models for Navier-Stokes equation, EPL, 17, 479-484, 1992, doi: 10.1209/0295-5075/17/6/001.
[Chen1998]S Chen and GD Doolen, Lattice Boltzmann Method for fluid flows, Annual Review of Fluid Mechanics, 30, 329-364, 1998, dos:10.1146/annurev.fluid.30.1.329.