Note: The following is adapted from 'LBflow:
an extensible lattice Boltzmann framework for the simulation
of geophysical flows. Part I: theory and implementation.'
which is available for download here.
Introduction to the theory of the lattice Boltzmann
Historically, the lattice Boltzmann
method developed heuristically, from attempts to improve
the performance of the lattice gas cellular automaton
(McNamara and Zanetti, 1988; Higuera and Jimenez,
1989). Subsequent analysis has put the lattice Boltzmann
method on a sound theoretical footing by showing that,
providing that certain conditions are met, the lattice
Boltzmann equations are formally equivalent to the
Navier–Stokes equations, discretized in time
and space (He and Luo, 1997; Chen and Doolen, 1998).
The literature contains several excellent and rigorous
theoretical treatments of the lattice Boltzmann method
(Succi, 2001; Nourgaliev et al., 2003; Yu et al.,
2003) which this section does not aim to duplicate.
Rather, the relevant equations are presented and their
physical meaning is discussed, echoing the method’s
According to the kinetic theory of
gases, a control volume of gas, with spatial dimensions
larger than the molecular mean free path and smaller
than the hydrodynamic lengthscale of the fluid, contains
a large number of molecules moving in all directions
and at all velocities and exchanging momentum through
collisions with one another. The probability that
a molecule within this volume, centred on position
r, at time t, has velocity v, is described by the velocity distribution function f(r,v,t). Collisions between the molecules act rapidly to bring
the velocity distribution function to an equilibrium
described by the Maxwell-Boltzmann distribution:
where ρ is the gas density, R the universal gas constant, T the absolute temperature of the gas and u is the bulk velocity of the gas packet (u = 0 for a gas at rest). For each orthogonal component
α of velocity, the Maxwell-Boltzmann distribution
is a normal distribution about uα. Note that the bulk velocity of the gas must
be much smaller than the root mean square velocity
of the molecules. Figure 1
(a and b) shows, schematically, the Maxwell-Boltzmann
distribution of velocities for a point in a 2-dimensional
flow around a cylinder; note that figure 1
b shows the distribution in velocity space.
||Diagrammatic representation of the velocity distribution function
for a small packet of flowing
fluid. a) Two-dimensional
flow around a cylinder; b)
velocity distribution, in
velocity space, for packet
of fluid with mean velocity
u; c) the same distribution discretized onto a two-dimensional,
The lattice Boltzmann method discretizes
space into a regular lattice of nodes with spacing
Δx. The lattice must have a high degree of symmetry and,
typically, lower-dimension projections of a four-dimensional
face-centred hypercubic lattice are used. All of the
fluid molecules within a volume Δx3 centred on a node are treated as if they exist at
the node. Each node, therefore, is analogous to the
control volume discussed above, and has a velocity
distribution associated with it. Since time is also
discretized, into timesteps of duration Δt, velocity space is also discrete. It is important to
note that there are two distinct sets of units that
are used in the numerical implementation of the lattice
Boltzmann method: physical units and simulation units
(denoted by a ‘ ’ above a quantity).
For simplicity the simulation units are usually chosen
as Δ = 1 and Δ= 1. Distances in the
simulation are, therefore, represented in units of
lattice-spacings and times in units of timesteps.
The mapping between simulation and physical units,
and vice versa, is discussed later in this section.
In a three-dimensional lattice, each
node has 26 immediate neighbours. If we impose the
constraint that information may only propagate from
a node to its nearest neighbours in a single timestep,
this gives 27 ‘allowed’ velocities (including
a ‘rest velocity’ for stationary molecules).
The velocities are denoted i(i = 0…n) where n is the number of allowed velocities. On a lattice, the
continuous velocity distribution f(r,v,t) can be discretized into components fi. fi can be thought of as the fraction of the total mass
of fluid at a node that is moving with each of the
allowed velocities. Figure 1
c shows the discretization of the continuous Maxwell-Boltzmann
distribution of velocities for a two dimensional model,
with nine allowed velocities.
Although each node in a three dimensional
lattice has 26 nearest neighbours, the most commonly-used
three-dimensional lattice Boltzmann models have only
15 or 19 allowed velocities. Figure 2
shows the allowed velocities for a so-called
D3Q15 lattice (three-dimensional, fifteen velocity).
The 15 velocities are:
The lattice velocities i are expressed in simulation units, hence, a velocity
of [1,0,0] moves one lattice-spacing in the x–direction each timestep.
The lattice Boltzmann algorithm
At each node, at the beginning of
each timestep, each fraction of mass fi propagates to the neighbour node in the direction
of velocity i (or remains at the current node for f0). Post-propagation, therefore, each node has a new
set of fi which have arrived from the neighbour nodes. The
new total density at each node is given by the sum
of the fractional masses arriving at the node:
where the mass is expressed in simulation
units. The lattice is typically initialized such that
the density at each node = 1. The average
fluid velocity at a node is found
by first multiplying the mass fraction fi by the velocity i for each velocity and summing, to give the momentum
per unit volume, then dividing by the density at the
node (equation 2
for example, given that Δ ∕ Δ= 1:
For a calculated fluid density and
average velocity, there will be a corresponding equilibrium
velocity distribution with the same density and average
velocity (mass and momentum are conserved). On a lattice,
this is given by a discrete version of equation 1
where wi are ‘weights’ associated with each velocity,
typically, for a D3Q15 lattice: w0 = 2 ∕ 9 (rest weight), w1…6 = 2 ∕ 9 (orthogonal weights) and w7…14 = 1 ∕ 72 (diagonal weights). The quantity ĉs = ĉ ∕ (where ĉ = Δ ∕ Δ ) is referred to as
the lattice pseudo-sound-speed. This value of ĉs is valid for a D3Q15 lattice with the given value of w0 but may differ for other geometries. The derivation
of the discretized equilibrium velocity distribution
) from the continuous equation (1
) assumes constant temperature and low Mach number
Post-propagation, the incoming distribution
fi is adjusted towards the equilibrium distribution
fieq. Physically, this process is analogous to the redistribution
of momentum that occurs when molecules collide with
one another and which gives rise to viscosity. The
degree of this adjustment is controlled by the dimensionless
relaxation parameter τ. After collision, time is incremented by Δ , the new fractions
fi are propagated to the neighbour node in the direction
of velocity i and the procedure repeats. Both the propagation
and collision steps are encapsulated in the lattice
The collision operator (r.h.s. of
) is only valid when the difference between the
incoming distribution and the equilibrium distribution
is small. This version of the lattice Boltzmann equation,
with a single relaxation time, is called the LBGK
model and is the most commonly used lattice Boltzmann
model. Other schemes with multiple relaxation times
are more complex, but may offer greater numerical
stability (Lallemand and Luo, 2000; d’Humières
et al., 2002).
One of the major, generic advantages
of the lattice Boltzmann method is that operations
at each node are local. This means that, at each timestep,
operations on the state of each node only require
information held at that node (this information may
have propagated there from an immediate neighbour
at the beginning of the timestep). This locality of
operations eliminates the need for time-consuming
iterative convergence at each timestep, which is characteristic
of many conventional computational fluid dynamic techniques.
Boundaries and driving flow
One of the great advantages of the
lattice Boltzmann method over conventional computational
fluid dynamic techniques is its simple handling of
boundary conditions. Figure 3
illustrates the most common method for implementing
||Fluid–solid boundary handling. The solid is shown in dark
grey, bounded by a solid black
line. Lattice nodes within
the solid are shown as solid
black squares, nodes within
the fluid with at least one
solid neighbour (boundary
nodes) are shown as grey-filled
squares. Mass components propagating
from a boundary node towards
a solid node are reflected
and returned to the node of
origin within the same propagation
step. This amounts to siting
the fluid–solid boundary
half way between the nodes
(dashed line). The ‘staircased’
solid represented within the
simulation is shown in light
grey, bounded by a dashed
Lattice nodes that are within a solid
do not have a fluid distribution associated with them.
Fluid nodes that have at least one nearest neighbour
that is a solid node are designated as boundary nodes
and require special treatment. During the propagation
step, boundary nodes do not receive an incoming mass
component from solid nodes. Conversely, the mass component
propagating from a boundary node to solid node would
be lost since solid nodes do not have a fluid distribution.
The solution is to reverse the outgoing mass component
and send it back to the boundary node from which it
originated. This so-called ‘half-way bounce
back’ boundary condition amounts to siting the
fluid–solid interface half way between the boundary
and solid nodes and perpendicular to the vector joining
the nodes. The practical impact of the resulting ‘staircasing’
of the boundary is considered in Part II.
There are many ways of driving the
flow in a lattice Boltzmann simulation (Succi, 2001,
provides an overview). One of the simplest and most
flexible methods is to impose a constant and uniform
body force G on the fluid at every point, physically analogous to
a gravitational force acting on the fluid. This is
achieved by adding an extra term gi to each of the mass components fi prior to propagation. gi is formulated in such a way that the total mass
is conserved but the momentum is adjusted to account
for the force acting on the mass at the node for the
duration of the timestep:
This term is added to the r.h.s. of
. Care should be taken to ensure that gi ≪ fi so that the condition that the difference between
incoming and equilibrium distributions is small is
Geophysical flow problems involve
physical units of length, time and mass which are
distinct from the simulation units used internally
by the lattice Boltzmann algorithm. For a simulation
to be of practical use, we require a mapping between
simulation and physical units: L0 → x, T0 → t and M0 → m. The mapping allows us to convert, for example, velocities
from lattice spacings per timestep in the simulation,
to metres per second in the modelled system.
The length mapping L0 is imposed by the user’s choice of resolution
of the simulation. For example, if we model flow along
a square-sided pipe with side 1 mm and choose our
simulation lattice resolution such that there are
20 nodes across the pipe, the lattice spacing Δx = 5 × 10-5m. Since Δ ≡ 1, L0 = Δx ∕ Δ = 5 × 10-5m.
The mass mapping M0 is obtained from the density of the fluid. Suppose
that, in the current example, we are simulating the
flow of water which, at 20∘C, has a density ρ = 998.29Kg m-3. The simulation is initialized with a density = 1 at all nodes.
Since density has dimensions [M][L]-3, ρ = M0L0-3 , hence, M0 = 1.25 × 10-10Kg.
The time mapping T0 is determined from the user’s choice of fluid
kinematic viscosity (ν = 1.004 × 10-6m2s-1 for water at 20∘C). In the LBGK model, the viscosity is related
to the relaxation parameter and the lattice pseudo-sound-speed:
and has dimensions: [L]2[T]-1. Note that this definition of the lattice viscosity
implies that the relaxation parameter τ > 0.5. From the dimensions of viscosity, we can see that ν = L02 ∕ T0, hence, we can obtain T0. With a relaxation parameter τ = 1, = 1 ∕ 6, hence, in the current example, T0 = 4.15 × 10-4s.
LBflow handles dimensional mapping internally. The user
inputs parameters in physical units and LBflow outputs results in physical units. It is, however,
important that the user understands the relationship
between the mappings because they place practical
restrictions on the flows that can be accurately simulated.
In the calculation above, the physical
timestep Δt is determined by the choice of resolution and relaxation
parameter τ. For problems involving time dependent flow, it
is necessary that the timestep is much smaller than
the timescale over which the flow evolves. Suppose
that the pipe in the current example had a side of
1 m, rather than 1 mm. It is unlikely to be practical
to maintain the resolution Δx = 5 × 10-5m due to memory limitations so a new resolution of
Δx = 0.05m might be chosen. If all other parameters were kept
the same, the new physical timestep, given by equation 8
, would be T0 = 415s. Clearly this would be unsuitable for simulating
unsteady flow. Theoretically, this problem can be
overcome by varying the relaxation parameter, in this
case by making (τ - 0.5) very small. However, very small values of τ lead to very slow convergence of the simulation
resulting in unfeasibly long simulation runs. Conversely,
the accuracy of the simulation decreases markedly
for large values of τ, effectively restricting the relaxation parameter
to the range 0.5 < τ ≲ 3.
Care should also be taken to ensure
that the low Mach number assumption of the lattice
Boltzmann method is not violated (section 1.1
). The lattice pseudo-sound-speed gives the effective
maximum speed at which information can propagate through
the lattice. The error due to non-zero Mach number
is O(Ma3) (Sterling and Chen, 1996) where Ma = ∣max∣ ∕ ĉs. To keep this error to less than 1%, the maximum
flow velocity should be ∣max∣≲ĉs ∕ 5. The maximum expected flow velocity for a simulation
should be estimated in physical units, giving the
condition ∣umax∣≲ Δx ∕ 10Δt. Given equation 8
, the fact that Δx ∕ Δt = L0 ∕ T0, and assuming τ = 1, we have a condition for the minimum resolution
required for a given ∣umax∣:
This condition and consideration of
memory resources restrict the range of Reynold’s
number that can be achieved. The Reynold’s number
is given by Re = uL ∕ ν where u and L are the characteristic velocity and lengthscale respectively.
If we assume u = ∣umax∣ and the characteristic lengthscale is resolved
by n nodes (i. e. L = nL0), then, from equation 9
, we obtain:
to ensure low Mach number flow. High
Reynold’s number simulations, therefore, require
high resolutions and are computationally expensive,
particularly for flows in which L is small compared with the size of the simulation domain.
Modelling real materials
The lattice Boltzmann LBGK model developed
above is, strictly, only valid for a weakly-compressible
ideal gas; the ‘weakly compressible’ condition
results from the slow speed of sound on the lattice.
In practice, however, the LBGK model is found to be
accurate for all Newtonian fluids, including incompressible
liquids, at low Mach number. The LBGK model is encapsulated
in the LBflow module ‘lbgk’ which is tested and validated extensively in
Lattice Boltzmann models for more complex
materials now abound in the literature. The approach
that these models take to introducing more complex
material properties varies. Common approaches include
adjusting the equilibrium distribution equation to
introduce non-ideal behaviour (e. g. Swift
et al., 1995; Luo, 1998) and propagating multiple
species simultaneously on the same lattice (e. g. Shan
and Chen, 1993; Kang et al., 2002). Other methods
iteratively couple a lattice Boltzmann model for fluid
flow with some other model which describes the behaviour
of another phase, for example a lattice-spring model
for a surrounding elastic shell (Buxton et al., 2005)
or a discrete element model for suspended particles
(Cook et al., 2004).
The modular design of LBflow allows users to code modules which encapsulate the
many lattice Boltzmann models that have been developed.
It should be possible to create modules for most lattice
Boltzmann models which preserve locality of operations
(see end of section 1
). LBflow also supports iterative coupling to other models.
The design and operation of LBflow are discussed in detail in the following section.
Buxton, G.A., Verberg, R., Jasnow, D., Balazs, A.C., 2005.
Newtonian fluid meets and elastic solid: Coupling
lattice Boltzmann and lattice-spring models. Physical Review, E 71, Article No: 056707.
Chen, S., Doolen, G.D., 1998. Lattice Boltzmann method for fluid
flows. Annual Review of Fluid Mechanics, 30, 329–364.
He, X., Luo, L-S., 1997. Theory of the lattice Boltzmann method:
from the Boltzmann equation to the lattice Boltzmann
equation. Physical Review, E 55, 6811–6817.
Higuera, F., Jimenez, J., 1989. Boltzmann approach to lattice
gas simulations. Europhysics Letters, 9, 663–668.
Kang, Q., Zhang, D., Chen, S., He, X., 2002. Lattice Boltzmann
simulation of chemical dissolution in porous media.
Physical Review, E 65, Article No: 036318.
Lallemand, P., Luo, L-S., 2000. Theory of the lattice Boltzmann
method: dispersion, dissipation, isotropy, Galilean
invariance, and stability. Physical Review, E 64, 6546–6562.
Luo, L-S., 1998. Unified theory of the lattice Boltzmann models
for nonideal gases. Physical Review Letters, 81, 1618–1621.
McNamara, G., Zanetti, G., 1988. Use of the Boltzmann equation
to simulate lattice gas automata. Physical Review Letters, 61, 2332–2335.
Nourgaliev, R.R., Dinh, T.N., Theofanous, T.G., Joseph,
D., 2003. The lattice Boltzmann equation method:
theoretical interpretation, numerics and implications.
International Journal of Multiphase Flow, 29, 117–169.
Shan, X., Chen, H., 1993. Lattice Boltzmann model for simulating
flows with multiple phases and components. Physical Review, E 47, 1815–1819.
Sterling, J., Chen, S., 1996. Stability analysis of Lattice
Boltzmann methods. Journal of Computational Physics, 123, 196–206.
Succi, S., 2001. The lattice Boltzmann equation for fluid dynamics
and beyond. Oxford University Press, ISBN: 0198503989,
Swift, M.R., Osborn, W.R., Yeomans, J.M., 1995. Lattice Boltzmann
simulations of nonideal fluids. Physical Review Letters, 75, 830–833.