LBflow > Lattice Boltzmann theory

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.

Figure 1: 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, nine-velocity lattice.

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 Δx³ 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 f_i. f_i 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 D₃Q₁₅ lattice (three-dimensional, fifteen velocity).

Figure 2: D3Q15 lattice. Model has one rest direction (0), six orthogonal directions (16) and eight diagonal directions (714).

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 f_i propagates to the neighbour node in the direction of velocity _i (or remains at the current node for f₀). Post-propagation, therefore, each node has a new set of f_i 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 f_i 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 w_i are ‘weights’ associated with each velocity, typically, for a D₃Q₁₅ lattice: w₀ = 2 ∕ 9 (rest weight), w_1…6 = 2 ∕ 9 (orthogonal weights) and w_7…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 D₃Q₁₅ lattice with the given value of w₀ 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 f_i is adjusted towards the equilibrium distribution f_i^eq. 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 f_i 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.

Figure 3: 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 line.

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 g_i to each of the mass components f_i prior to propagation. g_i 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 g_i ≪ f_i 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: L₀ → x, T₀ → t and M₀ → 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 L₀ 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, L₀ = Δx ∕ Δ = 5 × 10^-5m.

The mass mapping M₀ 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, ρ = M₀L₀^-3 , hence, M₀ = 1.25 × 10^-10Kg.

The time mapping T₀ is determined from the user’s choice of fluid kinematic viscosity (ν = 1.004 × 10^-6m²s^-1 for water at 20^∘C). In the LBGK model, the viscosity is related to the relaxation parameter and the lattice pseudo-sound-speed:

(8)

and has dimensions: [L]²[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 ν = L₀² ∕ T₀, hence, we can obtain T₀. With a relaxation parameter τ = 1, = 1 ∕ 6, hence, in the current example, T₀ = 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 T₀ = 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(Ma³) (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 ∣u_max∣≲ Δx ∕ 10Δt. Given equation 8 , the fact that Δx ∕ Δt = L₀ ∕ T₀, and assuming τ = 1, we have a condition for the minimum resolution required for a given ∣u_max∣:

(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 = ∣u_max∣ and the characteristic lengthscale is resolved by n nodes (i. e. L = nL₀), 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