## LBflow > Lattice Boltzmann theory

 Home Research Publications LBflow LB theory Downloads Examples Reference Collaborators PhD/postdoc opportunities PubVolc Contact Links Site last updated30th June 2011
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 method

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 heuristic origins.

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:
 (1)

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.

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 = 0n) 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.

## 1 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:
 (2)

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 ):
 (3)

for example, given that ΔΔ= 1:
 (4)

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 :
 (5)

where wi are ‘weights’ associated with each velocity, typically, for a D3Q15 lattice: w0 = 29 (rest weight), w16 = 29 (orthogonal weights) and w714 = 172 (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 equation (5 ) from the continuous equation (1 ) assumes constant temperature and low Mach number flow (∣≪ĉs).

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 Boltzmann equation:
 (6)

The collision operator (r.h.s. of equation 6 ) 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.

## 2 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 boundaries.

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:
 (7)

This term is added to the r.h.s. of equation 6 . Care should be taken to ensure that gi fi so that the condition that the difference between incoming and equilibrium distributions is small is not violated.

## 3 Dimensional considerations

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 20C, 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 20C). In the LBGK model, the viscosity is related to the relaxation parameter and the lattice pseudo-sound-speed:
 (8)

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, = 16, hence, in the current example, T0 = 4.15 × 10-4s.

### 3.1 Practical limitations

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ĉs5. 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:
 (9)

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 = umaxand the characteristic lengthscale is resolved by n nodes (i. e. L = nL0), then, from equation 9 , we obtain:
 (10)

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.

## 4 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 Part II.

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.

### References

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, 288pp.

Swift, M.R., Osborn, W.R., Yeomans, J.M., 1995. Lattice Boltzmann simulations of nonideal fluids. Physical Review Letters, 75, 830–833.