Lattice Boltzmann

Ludwig Boltzmann showed that a fluid’s macroscopic behaviour could be understood as the manifestation of vast numbers of ‘billiard-ball’ type collisions between individual molecules. Each molecule may possess a near-infinite variety of positions and velocities, and tracking such many-body collisions for a significant fluid volume remains significantly out of computational reach.

Boltzmann discovered a statistical solution to this problem, simultaneously removing the need to consider individual molecules and providing a convincing argument for the existence of atoms before they were universally accepted @Cercignani:1998. Instead of tracking the paths of individual molecules from collision to collision, Boltzmann saw that within a certain volume of space there would be a distribution of molecular velocities and that the averaged effect of collisions would be the redistribution of energy. The statistical mechanics of the molecules were sufficient to understand the fluid’s properties @Chandler:1987.

Development of statistical mechanics (notably by Maxwell, Chapman and Enskog) extended its reach beyond gases into liquids, and showed that the Navier-Stokes equation, describing the behaviour of a fluid with velocity $\bf u$, density $\rho$, pressure $p$ and viscosity $\mu$ subject to body force $\bf F$

$$ \rho \left( \frac{d{\bf u} }{dt} + {\bf u \cdot \nabla {\bf u} } \right) = - \nabla p + \mu \nabla^2{\bf u} + {\bf F}, \label{eqn:navier-stokes} $$

could be derived from the Boltzmann equation @Cercignani:1987. The concluding chapter of @Chapman:1990 gives a historical summary of the progress towards an accepted microscopic model, due to Enskog, for the transport of dense gases and liquids, and for a general introduction, see for example @Bird:1960.

The lattice Boltzmann (LB) method further restricts the distribution of particles to discrete velocities and locations, which are more amenable to numerical solution. Using lattice Boltzmann, fluid dynamics are simulated through a cycle of propagation and collision of particle densities on a regular lattice. This has been shown @Qian:1992 [@Chen:1992] to solve the incompressible Navier-Stokes equations for the class of problems considered here, and this is used as the hydrodynamic solver in the present work.

This chapter presents the theoretical foundations of the lattice Boltzmann models used for the present study, along with validation and verification.

Hydrodynamic Model

The Boltzmann equation, tracking the evolution of particle distributions $f_c$, with microscopic velocities $\xi_c$ is

$$ \label{eq:boltzmann_equation} \frac{\partial f_c}{\partial t_c} + \xi_c \partial f_c = \Omega_c. $$

The $c$ subscript denotes the continuous, classical parameters, to distinguish from their discretised counterparts, which are introduced in the following section. This continuous Boltzmann equation updates $f_c$ according to propagation along microscopic velocities $\xi_c$, and a collision behaviour defined by the collision operator $\Omega_c$.

The widely used Bhatnagar-Gross-Krook (BGK) model @Bhatnagar:1954 considers the collision as a linear relaxation from the non-equilibrium (propagated) particle distribution towards an equilibrium state $f^{eq}_c$, with a relaxation time $\lambda_c$:

$$ \Omega_c = -\frac{1}{\lambda_c}(f^{eq}_c-g_c). $$

The equilibrium distribution state is defined by the Maxwell-Boltzmann distribution function, which depends on the mean velocity ${\bf u}_c$, temperature $T$ and density $\rho_c$, as well as the gas constant $R$:

$$ \label{eq:maxwell-boltzmann} f^{eq}_c = \frac{\rho_c}{2\pi RT} \exp\left(\frac{-({ {\bf \xi}_c - {\bf u}_c})^2}{2RT}\right). $$

The lattice Boltzmann equation restricts the particle densities $f_i$ to discrete lattice nodes $\bf r$, connected by $q$ vectors ${\bf \xi}_i$ ($i=[1..q]$ and note the switch to discretised parameters), along which the densities propagate with velocity ${\bf c}_i$. The values of ${\bf c}_i$ are chosen so that the particle densities $f_i$ propagate to nearby lattice nodes ${\bf r} + {\bf\xi}_i$ in precisely one timestep $\Delta t$. That is, ${\bf\xi}_i = {\bf c}_i\Delta t$.

So rather than tracking the movement and interaction of vast numbers of individual molecules in free space, the lattice Boltzmann method tracks fictitious particle groupings $f_i$ moving in unison between points on a lattice $\bf r$. This spatio-temporal discretisation restricts the Boltzmann equation enough to simulate it numerically, while retaining enough of the underlying physics to produce faithful models of reality. Discretising ([eq:boltzmann~e~quation]) with the BGK collision operator in time (see @He:1997 for details of the procedure), and neglecting terms of order ${\Delta t}^2$ or smaller, leads to the equation solved in the lattice Boltzmann method,

$$ \label{eq:lbm} f_i\left({\bf r}+{\bf\xi}_i,t+\Delta t\right) - f_i({\bf r},t) = -\frac{1}{\tau}\left(f_i({\bf r},t) -f^{eq}_i(\rho, {\bf u} ) \right), $$

where $\tau$ is a rescaled equivalent of the relaxation time $\lambda$. For later convenience we also define a relaxation frequency, $\omega = \tau^{-1}$, and also the spacing of neighbouring nodes $\Delta x = \min{(|{\bf\xi}_i|)}$.

This equation simulates fluid flow through propagation ($f_i\left(\bf r+{\bf\xi}_i,t+\Delta t\right) - f_i(\bf r,t)$) and collision ($(f_i-f^{eq}_i)/\tau$) of particle densities on a lattice, which is generally regular. When used for modelling physical systems, these operations are handled separately, since phase-boundaries are handled between the propagation and collision steps.

image [fig:sandstone2d]

The equilibrium particle density, $f^{eq}$, derived from the Maxwell-Boltzmann distribution ([eq:maxwell-boltzmann]) retaining only terms up to order $\bf u^2$ during discretisation, is given by

$$ f_i^{eq}(\rho,\bf u) = w_{i}\rho \left(1 + \frac{\bf c_i\cdot\bf u}{RT} + \frac{(\bf c_i\cdot\bf u)^2}{2(RT)^2} -\frac{\bf u^2}{2RT} \right), \label{eq:feq} $$

where $w_i$, weightings from partition of unity, and the value of $RT$ are determined by the lattice geometry (and is $\frac{1}{3}$ for the D2Q9 lattice used here). Due to the truncation, this equation is valid only when $\frac{u}{RT} = Ma \ll 1$, with a value of $u=0.15$ regarded as a reasonable maximum @He:1997. Above this value, significant numerical inaccuracies can occur.

The above equations define the single phase hydrodynamic LBGK model @Qian:1992 [@Chopard:1998; @Kehrwald:2002]. As with the continuous Boltzmann equation, hydrodynamic variables are calculated from the (microscopic velocity) moments of the particle distributions:

$$ \rho = \sum_{i=1}^{q}{f_i(\bf r,t)} \label{eq:lbm_density} $$
$$ \label{eq:lbm_momentum} \rho\bf r = \sum_{i=1}^{q} \ci {f_i(\bf r,t)}. $$

Unlike in other numerical methods for solving flow problems, the local pressure $P$ is not an independent variable, being linearly related to the local density by the relationship

$$ P=\frac{1}{3} \rho. \label{eq:lbm_pressure} $$

The liquid’s kinematic viscosity is given by

$$ \nu = \frac{1}{3}\left( \frac{1}{\omega} - \frac{1}{2} \right) \frac{\Delta x^2}{\Delta t}, $$

imposing a condition of $\omega<2$ for positive viscosity.

Auxiliary conditions, including initial and boundary conditions (see Section [sec:boundary~c~onditions]), are required when simulating physical systems. It is then necessary to solve ([eq:lbm]) in two steps; propagation

$$ \label{eq:propagation} \ftil(r, t + \Delta t) = f_{i}(\bf r - {\bf \xi}_{i}\Delta t, t), $$

and collision

$$ \label{eq:collision} f_{i}(\bf r + {\bf \xi}_{i}, t + \Delta t) = \left(1 - \omega \right) \ftil(\bf r, t) + \omega \feq(\rho(\bf r,t), \bf r(\bf r, t)). $$

The new parameter $\ftil$ represents the non-equilibrium state of the system following propagation, before collisions which result in relaxation to equilibrium. This is only an intermediate state, and not a solution to the system at any time.

Lattice Boltzmann Algorithm

The method is summarised, excluding auxiliary conditions, in the following algorithm

Given an initial state defined in terms of density $\rho({\bf r},t)$ and velocity $\bf u({\bf r},t)$, compute the evolution of those parameters in a hydrodynamic system for time $t=[0,t_{max} ]$ using the lattice Boltzmann equation with BGK collision operator.

  • [Initialise.] Set $t\leftarrow 0$. Initialise $\rho(\bf r,0)$, $\bf u(\bf r,0)$ and $t_{max}$. Set $f_i(\bf r,0) \leftarrow f^{eq}_i(\rho,\bf u)$ using ([eq:feq]).
  • [Calculate equilibrium.] Set $f^{eq}(\rho,\bf u)$, calculated using ([eq:feq]).
  • [Collision.] Update particle distributions on each node, as a function of those from the previous state and the equivalent equilibrium state, using ([eq:collision]). This is often referred to as a relaxation, since the BGK collision process is a relaxation to the equilibrium state.
  • [Propagation.] Move particle densities along their propagation vectors $\xi_i$ to adjacent nodes ([eq:propagation]).
  • [Completed timestep.] Set $t \leftarrow t+\Delta t$.
  • [State parameters.] Calculate $\rho(t,\bf r)$ and $\rho\bf r(t,\bf r)$ using ([eq:lbm~d~ensity],[eq:lbm~m~omentum]).
  • [Simulation complete?] If $t go back to A2; otherwise the algorithm terminates. Depending on the system being simulated, other convergence checks may be more appropriate (such as change in $\Delta \bf u/\Delta t$ or $\Delta \rho / \Delta t$ being within some tolerance). [alg:lbm]

It will be noticed that this algorithm is computationally simple, containing only a relatively small number of inexpensive floating-point arithmetic operations. This is in contrast to numerical methods which solve the Navier-Stokes equations using a system of linear equations, whose solution (for example, by Gaussian elimination and back-substitution @Smith:2004) is comparatively expensive. It is this simplicity and the locality which make the scheme amenable to parallelisation. Short, complete code examples for a 2D porous domain flow model and a 1D distributed-memory diffusion model, on a single page of Matlab and C++ respectively, are contained in Appendix [sec:code].

It is interesting to note that, with $\tau=\omega=1$, $f_i$ correspond precisely to $f_i^{eq}$ following collision, allowing for a much simplified implementation, since the non-equilbrium densities do not need to be retained. Boundary conditions may also be implemented solely as macroscopic variables (i.e. equilibrium macroscopic quantities are imposed directly, rather than calculating $f_i$), also allowing greater efficiency. Under this constraint, the diffusive LBGK model corresponds to a Lax-Wendroff scheme @vanderSman:2000. In common with the majority of published studies, the current work retains the flexibility in $\Delta t$ and $\Delta x$ afforded by $\tau\neq 1$, which results in higher computational efficiency.

Lattice Geometries

Having described the lattice Boltzmann method, all that remains before it may be implemented is a valid lattice. Lattices, whose nodes are generally equally-spaced, are conventionally named D$d$Q$q$ according to their dimensionality $d$ and number of distribution vectors at each node $q$. Several acceptable lattices exist in one, two and three dimensions, which use ([eq:feq]) with different weightings and values of $RT$.

The lattices considered here, D1Q3, D2Q9 and D3Q19 are those most commonly used in engineering applications of LBGK. Each is athermal and has a constant value of $RT=\frac{1}{3}$. Their distribution vectors are defined as:

$$ [\xi_0,\xi_1,\xi_2] = \left[\begin{array}{rcc} -1 & 0 & 1 \end{array}\right], $$
$$ [\xi_0,\xi_1,\xi_2,\xi_3,\xi_4,\xi_5,\xi_6,\xi_7,\xi_8] \\ = \left[ \begin{array}{cccrrrrrr} 0 & 1 & 1 & 0 & -1 & -1 & -1 & 0 & 1\\ 0 & 0 & 1 & 1 & 1 & 0 & -1 & -1 & -1 \end{array} \right] $$
$$ \begin{aligned} & [\xi_0,\xi_1,\xi_2, \xi_3,\ldots,\xi_{17}, \xi_{18}] = \nonumber \\ & { \small \begin{bmatrix} 0 & 1 & 0 & 0 & -1 & 0 & 0 & 1 & -1 & -1 & 1 & 0 & 0 & 0 & 0 & 1 & -1 & -1 & 1\\ 0 & 0 & 1 & 0 & 0 & -1 & 0 & 1 & 1 & -1 & -1 & 1 & -1 & -1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1 & 1 & 1 & -1 & -1 \end{bmatrix} }\end{aligned} $$

The corresponding weightings $w_i$ are listed in Table [table:lb~w~eightings], and the 2D and 3D model vectors are displayed in Figure [fig:lbmvectors].


|c c|c c|c c| &&\ $|\xi_i|$ & $w_i$ & $|\xi_i|$ & $w_i$ & $|\xi_i|$ & $w_i$\ $0$ & ${2}/{3}$ & $0$ & ${4}/{9}$ & $0$ & ${1}/{3}$\ $\Delta x$ & ${1}/{6}$ & $\Delta x$ & ${1}/{9}$ & $\Delta x$ & ${1}/{18}$\ & & $\sqrt{2}\Delta x$ & ${1}/{36}$ & $\sqrt{2}\Delta x$ & ${1}/{36}$\



Other Lattice Boltzmann Formulations

Multiple Relaxation Times

The most popular LB method remains the lattice BGK, which uses a single relaxation time. Advantages of using multiple relaxation time (MRT) models have been studied analytically @DHumieres:2002 and by simulating flow through porous media @Pan:2005. Enhanced stability at higher flow speeds means that MRT models are capable of greater efficiency than LBGK formulations, and do not suffer from inaccuracies at fluid-solid boundaries when under-relaxed (i.e. with $\tau < 1.5$, see Section [sec:review~p~orous~d~omain~h~ydrodynamics]), but BGK remains popular due to its simplicity.

In the BGK model’s update equation, the particle densities $f_i$ depend only on the previous $f_i$ and the parameters $\rho$ and $\bf r$. There is only an indirect dependence, through $\rho$ and $\bf r$, on the other local particle densities. This simplifies the model, but limits numerical stability compared to other formulations. Multiple-relaxation-time (MRT) LB models @DHumieres:2002 retain the simple propagation described for LBGK, but replace the single relaxation parameter with a collision matrix $\mathcal{S}$: $|f_i\left(\bf r+\ci,t+\Delta t\right)\rangle - |f_i(\bf r,t)\rangle = -\mathcal{S}\left[|f_i(\bf r,t)\rangle -|f^{eq}_i(\rho, \bf r )\rangle \right],$ where $|f_i\rangle$ is a column matrix of $f_i$ and $\mathcal{S}$ is a $q\times q$ matrix representing the relative contributions of each $f_i$ at $t$ to the calculated $f_i$ at $t+\Delta t$. Note that LBGK is a special case of this MRT model.

The added complexity of allowing separate relaxation rates for each $f_i$ increases the computational time per $\Delta t$ per node by $\simeq 15\%$, but this can be counterbalanced by the increased stability @DHumieres:2002. For example, the speed attainable in a LBGK model is limited to $\bf r\leq0.1{\Delta x}/{\Delta t}$, above which value conservation of mass and momentum ceases. This restriction is due to assumption of low speed, made during derivation of the LB equation. The flexibility of MRT allows parameters to be tuned for maximum stability, and the restriction to be raised to $\leq0.19$. This can allow an overall reduction in solution time. Note that numerical simulations are needed to determine these limits, since no analytical solutions have been described.

Entropic Lattice Boltzmann

A separate approach to increasing LBM stability is via Boltzmann’s H-Theorem, which is a statement of the second law of thermodynamics from the perspective of statistical mechanics. It shows that the Boltzmann equation predicts systems whose entropy can only increase. By retaining higher orders during the discretisation of the Boltzmann equation, it is possible to construct a LB equation which locally obeys the H-Theorem. A global H-Theorem in a lattice Boltzmann model guarantees stability, by constraining the system to a long-term homogeneous distribution of particles. Such entropic LB (ELB) models @Ansumali:2005 depend upon an increased set of particle distribution vectors, for example a D2Q107 model with $|\bfc|\leq4\Delta x$. Extending the update rule to depend on more than near-neighbours clearly increases computational demands, but there are additional drawbacks. Boundary conditions are a problem with these models, both in terms of implementation (due to $f_i$ from many nodes crossing each domain boundary node), and since none have been shown to increase entropy unconditionally.

The use of the entropic lattice Boltzmann method has advantages, but is still not common, perhaps due to misunderstanding about its implementation. For example, an experienced and respected LB researcher recently claimed that ELB had no significant stability advantage over ‘traditional’ lattice Boltzmann models @Luo:2011. However, the study had implemented the ELB method incorrectly @Karlin:2011. Whatever the reason, no studies to date have reported using ELB for engineering applications.

Reactive Mass Transport LB

The advection-dispersion equation (ADE), defining the evolution of the concentrations $C^s$ of multiple species with chemical reactions may be expressed as

$$ \frac{\partial C^s}{\partial t} + ({\bf u \cdot \nabla}) C^s = {\bf \nabla}\cdot(D_m^s\nabla C^s) +R^s, \label{eq:RADE} $$

where $D_m^s$ and $R^s$ are respectively the molecular diffusion and the rate of generation of species $s$.

Modelling ([eq:RADE]) using lattice Boltzmann involves introducing a reaction term $r^s$ as a modifier to the collision operator in ([eq:lbm]), to define the chemical reaction source or sink. By linking the generation of one chemical species to the depletion of another, via rate constants introduced below, the species are coupled and chemical reactions are modelled.

The lattice Boltzmann equation solved in chemically-reactive mass transport problems is therefore

$$ g_i^s\left(\bf r+{\bf\xi}_i,t+\Delta t\right) - g_i^s(\bf r,t) = -\frac{1}{\tau^s}\left(g_i^s(\bf r,t) -g^{s,eq}_i(C^s, \bf r) \right) + \omega_ir^s, \label{eq:lbm_rd} $$

with $g_i^s$ the distribution function of the concentration of the $s^{th}$ solute and $\tau^s$ is its relaxation time. The LB reaction term is related to the macroscopic equivalent by

$$ r^s=\frac{R^s}{m_s}\Delta x^2 \Delta t, $$

where $m_s$ is the mass, taken here as 1, and $\Delta x$ and $\Delta t$ are spatial and temporal discretisations. For single-species (e.g. tracer) dispersion, $r= R^s = 0$. Modelling complex kinetic reactions between multiple reactants and products using this coupling term is not complicated (see e.g. @PonceDawson:1993 [@Yu:2002], or Section [sec:selkov]).

It can be seen that the non-reactive terms in ([eq:lbm~r~d]) correspond to those in the hydrodynamic model (i.e. ([eq:lbm])). Their equilibrium distribution is (compare with ([eq:feq]) for $f_i^{eq}$, noting that all common terms retain their meanings, and in particular that the values of the weightings $w_i$ are identical):

$$ g_i^{s,eq}(C^s,\bf r) = w_{i}C^s \left(1 + \frac{\bfc_i\cdot\bf r}{RT} \right), \label{eq:lbm_equilibrium_rd} $$

and the concentration of species $s$ is calculated from the zeroth moment of the particle densities

$$ C^s = \sum_{i}g_i^s. \label{eq:lb_rd_concentration} $$

This is identical to the calculation in the hydrodynamic model, except that no higher moments are required in calculating $g^{eq}$, as momentum is not a conserved quantity in this LB model. The advecting species are treated as passive tracers, with negligible mass compared to the liquid by which they are transported. This method is commonly known as an ‘overlay’ model, since the mass-transport LB simulation ([eq:lbm~r~d])–([eq:lb~r~d~c~oncentration]) is performed on top of the hydrodynamic simulation ([eq:lbm])–([eq:lbm~d~ensity]). In practice the similarities of the hydrodynamic and mass transfer models mean that only a small proportion of the implementation is exclusive to each. This overlay method has been used to model mineral dissolution @Kang:2002 and heat/mass transfer @Yoshino:2003.

It has been shown using the Chapman-Enskog expansion method @PonceDawson:1993 that these LB equations correspond to the reactive ADE ([eq:RADE]), with the diffusivity of each species given by (in parallel with the kinematic viscosity of the hydrodynamic model)

$$ D_m^s = \frac{1}{3} \left( \tau^s - \frac{1}{2} \right) \frac{\Delta x^2}{\Delta t}. \label{eq:lb_species_diffusion_constant} $$

Although not a requirement of the method, for all of the simulations reported here, the hydrodynamic and mass transport simulations operated sequentially. In common with related physical experiments @Jones:2005, analysis of the reactive dynamics was simplified by using a fully developed (steady-state) flow field to advect the chemicals. So the hydrodynamic simulation ran until convergence was reached, and that steady-state flow field (i.e. $\bf r$ in ([eq:lbm])) was used as the only flow field in the mass transport model (i.e. $\bf r$ in ([eq:lbm~r~d]) and ([eq:lbm~e~quilibrium~r~d])).

Boundary Conditions

The choice of appropriate boundary conditions, including phase and simulation-domain boundaries, is essential when solving partial differential equations for physical systems, regardless of the solution method. One of the main advantages of lattice Boltzmann over other discretisation methods, and the primary reason for its use in porous domain simulations, is the ease with which fluid-structure interfaces may be handled.

Having described in previous sections the update of an LB simulation from $t\rightarrow t+\Delta t$, auxiliary conditions are now detailed which are subsequently used to model reactions between chemical species transported through void networks. Boundary conditions for hydrodynamic and mass-transport LB models, which are physically quite different, are in several cases numerically similar and are therefore described with reference to one another in adjacent sections below. In some highlighted cases the handling of boundaries is identical and the particle densities of the hydrodynamic model, $f_i$, can be read interchangeably with the mass-transport equivalent $g_i$.

Initial conditions are typically defined by a simulation-specific macroscopic equilibrium state, from which starting $f^{eq}$ and $g^{eq}$ are calculated. Following propagation of $f_i$ (see algorithm A) adjacent to a phase or domain boundary defined by surface normal $\bf n$, there are unknowns for which the the distribution vectors $\xi_{i (i\cdot {\bf n} \geq 0)}$ originate outside the computational domain. Figure [fig:boundary~d~istributions] demonstrates this for a boundary ${\bf n}={\bf x}$, where the computational domain exists only in positive ${\bf n}$, with a D2Q9 lattice. This may represent for example be the boundary of an impermeable solid or the inlet of a flowing channel. Three particle densities ($f_{1,2,8}$) are unknown. For consistency the following sections each consider this configuration, although the descriptions are equally valid for other edges due to symmetry.

image [fig:boundary~d~istributions]

Boundary conditions are defined by the method of determining $f_{i (i\cdot{\bf n}\geq 0)}$, often based on the known densities $f_{i (i\cdot{\bf n}< 0)}$ in conjunction with pre-determined macroscopic properties. The definition of appropriate boundary conditions for LB remains under active discussion in the literature @Zhang:2002; the following sections describe the methods used in this research.

Periodic Boundaries

When the modelled domain is periodic in some direction, the unknown $f_i$ emerging from one edge of the domain are complemented by equivalent $f_i$ leaving the opposite side of the domain:

$$ \begin{aligned} f_1(t+\Delta t,0,y) &= f_1(t,x_{max},y)\\ f_2(t+\Delta t,0,y) &= f_2(t,x_{max},y-\Delta x)\\ f_8(t+\Delta t,0,y) &= f_8(t,x_{max},y+\Delta x).\end{aligned} $$

That is, the nodes at $x_{max}$ are separated by $\Delta x$ from the nodes at $x=0$, the $f_i$ travel between them in $\Delta t$, and nodal connectivity is identical to that elsewhere in the domain.

In the simulations of mass-transport past regularly-spaced obstructions presented below, the obstructions are periodic in the direction perpendicular to the initial pressure and concentration gradients. Periodic boundaries are applied on the other two domain edges (those without specified $P$ or $C$), both for $f_i$ and $g_i$, so that any mass leaving the domain re-enters at the opposite side.

Impermeable Boundaries

Impermeable boundaries defining the fluid-solid interface are treated in the present research as no-slip boundaries and are handled by reflection of particle densities $f_i$ that encounter a ‘solid’ node back to the fluid node from which they originated.

So, for Figure [fig:boundary~d~istributions], when interface $\bf n$ represents a solid boundary, the reflecting particle densities are defined as follows

$$ \begin{aligned} f_8(t+\Delta t) &= f_4(t),\\ f_1(t+\Delta t) &= f_5(t), \text{~and}\\ f_2(t+\Delta t) &= f_6(t).\end{aligned} $$

Since the $f_i$ re-emerge at their original positions after $\Delta t$, the location of the solid interface is precisely half-way between the two nodes. This is probably the most common method of defining impermeable boundaries with lattice Boltzmann, but as noted by another researcher @Humieres:2006, many studies reporting the use of these ‘bounceback’ impermeable boundaries do not state when the $f_i$ re-enter the fluid domain. The alternatives range from immediate re-entry (defining $f_8(t) = f_4(t)$ with the interface at the fluid node) to delayed (defining $f_8(t + 2\Delta t) = f_4(t)$ with the interface at the solid node). The accuracies of each method are studied in detail in @Humieres:2006.

More accurate boundary conditions exist @Noble:1996, but they cannot be directly used for more complex morphologies. Bounceback boundaries are useful because tomographic data from real sampled rocks can be used as direct inputs to LB simulations.

The accuracy of bounceback is lower (with a $1^{st}$ order rate of convergence) when flow is not parallel to the boundary @Gallivan:1997. This effect is accounted for in the simulations described here by ensuring that the discretisation is small enough that obstructions occupy a significant number of lattice spacings.

A significant source of the error in LB simulations of flow around a circular obstruction (or any other non-linear body) can be attributed to the inability of a square grid to represent a curved surface, for example in @Gallivan:1997 results for octagonal obstructions, which are better modelled by a square grid, exhibited $1\%$ errors at the same discretisation as circular obstructions showing $2\%$ error.

Curved boundaries are possible in lattice Boltzmann, for example by combining the simple bounceback method with spatial interpolations @Bouzidi:2001b, which also permits the modelling of moving boundaries. However, the simulations described in this research seek to compare non-reactive to reactive mass transport past regular obstructions, where the modelling of approximations to curved surfaces alters the simulations equally. That is, if the simulations are considered as mass transport past approximations to circles and other shapes, the errors when comparing results (for example spatial moments of mass distribution) will be lower than those reported in the above studies.

Pressure Boundaries

In the lattice Boltzmann scheme, pressure and density are linearly related by ([eq:lbm~p~ressure]), so imposing a nodal pressure may be achieved by specifying the equivalent density, with the restriction that the density fluctuation $\frac{\Delta P}{\Delta x}\ll \rho$ @He:1997c. To define a density at the domain boundary, $\rho_0$, some assumptions are required. We assume the $u_{\bf n}=0$, where

$$ \sum_{i=2,3,4}f_i - \sum_{i=6,7,8}f_i = 0, \mbox{ and } \label{eq:pbc_fi_zero} $$
$$ \sum_{i=1,2,8}f_i - \sum_{i=4,5,6}f_i = \rho_0 u_x \label{eq:pbc_rho_ux} $$

equating $f_1+f_2+f_8$ in ([eq:lbm~d~ensity]) and ([eq:pbc~r~ho~u~x]) gives

$$ u_x = 1 - \frac{\sum_{i=0,3,7}f_i + 2 \sum_{i=4,5,6}f_i}{\rho_0} \label{eq:pbc_ux} $$

Zou and He proposed that bounceback holds @Zou:1997 in the direction normal to the boundary, that is $f_5 = f_5^{eq} = f_1 - f_1^{eq}$, and since most of the terms of $f_1^{eq}-f_5^{eq}$ cancel (q.v. ([eq:feq])),

$$ f_1 = f_5 + \frac{2}{3}\rho_0 u_x. \label{eq:pbc_f1} $$

Using ([eq:pbc~f~1]) with ([eq:pbc~f~i~z~ero]) and ([eq:pbc~r~ho~u~x]) enables definition of the final two unknown particle densities,

$$ \begin{aligned} f_8 &= \frac{1}{6}\rho_0 u_x + f_4 + \frac{f_3-f_7}{2}, \text{~and}\\ f_2 &= \frac{1}{6}\rho_0 u_x + f_6 + \frac{f_7-f_3}{2}.\end{aligned} $$

To calculate flow fields $\bf u(r)$ for a given geometry, Dirichlet boundary conditions can be imposed at the inlet and outlet to produce a pressure drop $\Delta P$ along the channel of length $L$. The hydrodynamic simulation continues until a constant velocity field is reached @Shan:1993. The initial velocity field influences only the number of timesteps required to reach steady state, not the calculated field.

Gravitational Force

In addition to constant pressure or velocity boundary conditions at the boundary of a simulated domain, it is possible to consider a body force addition to the calculation of $f^{eq}$ acting over the entire mass of the fluid @Buick:2000. All of the models of fluid flow in the present document used a pressure gradient to drive the flow.

Constant Concentration Boundaries

In many applications of mass transport models, concentrations are fixed at various positions. For example, the inlet concentration of a channel may be held at a constant value in order to measure the resulting concentration profile at some downstream location. As discussed in Section [sec:constant~p~ressure], constant pressure boundaries are applied to a hydrodynamic LB model by altering the unknown $f_i$ to maintain a constant nodal density. The nodal densities of a mass-transport LB model may be similarly maintained altering $g_i$;

$$ \begin{aligned} g_1 &= g_5 + \frac{2}{3}C_0 u_x,\\ g_8 &= \frac{1}{6}C_0 u_x + g_4 + \frac{g_3-g_7}{2}, \text{~and}\\ g_2 &= \frac{1}{6}C_0 u_x + g_6 + \frac{g_7-g_3}{2}.\end{aligned} $$

Previous Lattice Boltzmann Studies

Having described the specific lattice Boltzmann models used in the current research, this section reviews previous work relevant to the aims of the present study.

In general, it can be said that the majority of the published reports of engineering applications of lattice Boltzmann models relate to the calculation of flow profiles of liquids, with fewer concerning mass transport and only a handful concerning the transport of reacting chemical species.

The outcomes of these previous studies demonstrate that the choice of lattice Boltzmann for the current work is appropriate, and that the work represents something absent from the existing literature.

Porous Domain Hydrodynamics

A study of flow around regular arrays of circular and hexagonal obstructions in 2D found agreement to within $0.5\%$ with benchmark finite difference simulations @Noble:1996, with a large range of velocities ($0.01<Re<100$). Using the bounceback method to model no-slip boundaries (see Section [sec:boundary~c~onditions]), errors at the boundary are known to depend on the value of the LB relaxation parameter @Succi:2001. In a subsequent study of flow around circles and hexagons, an upper limit ($\tau<1.5$) on the value was suggested @Gallivan:1997, and the reason for the limit is described in a study of a lattice Boltzmann model of gas mixtures @Bennett:2010.

Considering flow in more realstic domains, The hydraulic conductivity of artificially-generated soil samples was measured using lattice Boltzmann hydrodynamics in @Zhang:2005, finding that high-porosity soils were well-modelled, but low porosities were less so. This was suggested to be due to insufficient resolution (i.e. $\Delta x$ was too high), which is not suprising, and essentially reiterates that higher resolutions increase the accuracy of these numerical models @Martys:1999, which is true of any discretisation-based approach to solving partial differential equations. The above studies used a pressure-gradient to drive the flow of an incompressible LBM formulation @Zou:1997, and due to the stochastic generation of soil, only statistically homogeneous morphologies were considered. A similar study @Okabe:2004 conducted permeability measurements of reconstructed sandstone using lattice Boltzmann, and used two-point correlations to categorise and reconstruct the porous media.

Similar studies were conducted for flow through bead packs (with spheres of equal radius) @Maier:1998 and static mixer elements of fixed geometry @Kandhai:1998, finding good agreement with physical experiments, and local imhomogeneities in flow speed that were greater than expected. No investigation of the inhomogeneities was reported, and the studies were not concerned with the effect this localised behaviour would have upon mass transport through the domain.

Mass Transport

In an early use of lattice Boltzmann to model mass transport problems, models of diffusion rates in partially saturated porous media (overlapping randomly-placed spheres) were found to follow the semi-empirical Archie’s second law relationship to saturation @Martys:1999, using a $100 \Delta x^3$ domain. One investigated model successfully used monodisperse spheres of radius $9 \Delta x$, which is not dissimilar to some of the 2D domains used in the present studies.

Non-reactive dispersion of mass in a straight walled channel, known as Taylor-Aris dispersion @Aris:1956, has been approached as a benchmark problem independently for Poiseuille flow in a pipe @Lowe:2002, and in a 2D slit @Sukop:2005. Neither study described deviations from the analytical solution, nor modelled high $Pe$ systems. Both studies exclusively reported on the steady-state behaviour of the model.

Mass transport through porous domains has been investigated in many studies using lattice Boltmzann, including a study with a wide variation in Reynolds numbers @Yoshino:2003. Several schemes for incorporating elements of the porous domain at the level of individual computational nodes have been introduced in attempts to speed up calculations. Porosity and drag terms were added at each node in @Guo:2002 with good comparisons against existing finite difference results for benchmark problems in porous domains. Such methods are not used in the current studies, partly because no benchmarks or studies exist of their use for reactive mass transport problems. However, for more complicated domains than those used here, they could produce significant gains in computational speed.

Reactive Mass Transport

The early lattice Boltzmann studies of reactive flows are reviewed in @Chen:1998. Mixing of species at simple ‘crossroads’ fracture junctions, found a higher mixing rate predicted by LB simulation than would be predicted by streamline tracking @Stockman:1997 [@Stockman:2001], due to deviations from the simple laminar flow pattern. The importance of boundary conditions for modelling the fluid-solid interface was again emphasised.

Several studies have used the lattice Boltzmann method to simulate gas packed-bed (adsorption) reactors @Freund:2003 [@Manjhi:2006], finding quantitative agreement with observed physical simulation. Pockets of extremely fast flow (and consequent lower reaction product) at the edges of the reactor, due to a looser bead packing were found in both studies. As was noted above for related hydrodynamic simulations, although the local inhomogeneity was measured and reported, no study of the effect it had upon reaction rates was reported, and the distances before such inhomogeneities would dissappear were not investigated.

Simulating models where chemical reactions alter the morphology of the porous domain were made possible by a simple alteration to the lattice Boltzmann formulation allowing sub-grid boundary conditions @Verberg:2000. These have been used in a number of published studies, including a study of chemical erosion in rough fractures @Verberg:2002 and mineral dissolution and precipitation @OBrien:2002, each of which were validated against laboratory experiments. In each of these studies, LB was used for the hydrodynamic flow simulation only and different methods were used to model chemical transport. In the latter report, an analytical reaction model was used, while the former used a stochastic tracer diffusion algorithm. The reason for this was stated as a good handling of flows at high Peclet numbers, and that ‘grid-based schemes suffer significantly from numerical diffusion, leading to serious errors in the concentration distribution at Peclet numbers that are typical of flow in fractures’.

As described in Section [sec:lgca], similar considerations meant that this current research did not use a lattice gas method for the mass transport modelling. Although similar problems at relatively low $Pe$ have also been modelled using lattice Boltzmann for the chemical reaction @Kang:2002, the following chapter carefully validates the model used in this research at high Peclet numbers.


The implementation of the lattice Boltzmann model used in this research was written specifically for it by the author, and therefore needed to be validated. Some of the validation tests are described in this section in order to illustrate the method, with further examples of validations described in the following chapters.

Channel Flow

A liquid flowing through an unobstructed tube travels faster in the middle of the tube than near the edges. The speed profile is parabolic and described in 2D by

$$ u_{y} = \frac{Q}{2\mu}(L^{2} - y^{2}) \label{eq:poiseuille} $$

where $Q$ is the flow, $\mu$ the liquid’s viscosity, and $2L$ the width of the channel.

Using a lattice Boltzmann model with an initial velocity ${\bf u} = 0$ and uniform density, a flow was induced using a constant pressure gradient along a channel using the method in Section [sec:constant~p~ressure]. The resulting flow when a steady state had been reached is plotted in Figure [fig:poiseuille], and the characterisitic Poiseuille parabola was seen. For reference, the source code for running this simulation is listed in Figure [fig:lbm~l~ibrary~c~lient~c~ode].

This is a standard test for any numerical flow model, and detailed analysis of the behaviour of lattice Boltzmann when modelling this system is available in the literature, for example Chapter 6 of @Succi:2001 and Section 5.1 of @Sukop:2005.

image [fig:poiseuille]

Diffusion in Unbounded Domains

In the absence of reaction or advection, the advection-dispersion equation ([eq:RADE]) in one-dimension becomes

$$ \frac{dC}{dt}=D\frac{d^2 C}{dx^2} $$

where $C$, $t$ and $x$ have their usual meanings of concentration, time and location.

The analytical solution for diffusion from an inital Dirac delta source located at $x=0$ in an infinite empty domain in 1D is

$$ C(x,t)=\frac{M_0}{2\sqrt{\pi Dt} }e^{-\frac{x^2}{4Dt} }. \label{eq:diffusion_1d} $$

This result is compared with the result of a lattice Boltzmann simulation with $\omega =1$ at $t=3000\Delta t$ in Figure [fig:test1].

Validation of lattice Boltzmann model for diffusion in 1-D, comparison with ([eq:diffusion~1~d]) at $t=3000\Delta t$. [fig:test1]

Diffusion away from a Dirac delta pulse in two dimensions can be expressed by the equation

$$ C=\frac{M_0}{4\pi Dt} $$

A lattice Boltzmann simulation plotted against the analytical result at $y=0$ in Figure [fig:test2]. This is a relevant test because in grid-based discretisation methods for solving partial differential equations, it is possible for the underlying lattice (i.e. the spatial points at which the equation is solved) to influence the solution. The equivalent 2D LGCA diffusion model, used in the conference papers in the Appendix, exhibits from grid dependence which can require finer discretisations, and longer simulation times, to overcome @Chopard:1998.

Comparison of analytical and numerical results for 2-D unbounded diffusion at $t=200\Delta t$ [fig:test2]

Lattice Boltzmann Models of Chemical Reactions

In addition to hydrodynamic and mass transport modelling, lattice Boltzmann is used for mass transport modelling involving chemical reactions using ([eq:lbm~r~d]). The results of some simple lattice Boltzmann models of chemical reactions are presented in the following sections. These reactive simulations take place homogenously, between fully mixed reactants so that the simulation’s spatial extent is irrelevant, and the rate of reaction is entirely determined by the kinetic constants. Note that the D2Q9 model was used for all of these results. Validation of the first and third order reactions repeat published results. The second-order reaction validation does not appear in the published literature.

When chemical species are initially separated, their reaction will be further limited by the extent of mixing, which in a hydrodynamic scenario occurs by diffusion and advection. Validation of advection-diffusion models is included in Section [sec:reactive~m~ass~t~ransport]

First Order Chemical Reactions

A simple first-order chemical reaction involving a single reactant and product can be described as

$$ A\rightarrow B. $$

Assuming the reaction to be elementary, the rate equations are

$$ r_a \equiv \frac{dC_a}{dt} = -kC_a = -r_b \equiv -\frac{dC_b}{dt} = -kC_b \label{eq:rate_AB} $$

where $k$ is the kinetic constant, and $r$ and $C$ represent the reaction rates and concentrations of the species denoted by the subscript. To model this system with the lattice Boltzmann model described above, ([eq:lbm~r~d]) is solved for species $A$ and $B$, using the reactive rates defined by ([eq:rate~A~B]). With initial concentrations of $C_a=C_0$ and $C_b=0$, solving the rate equations gives

$$ C_a = C_0e^{-kt} \quad {\rm and} \quad C_b = C_0(1-e^{-kt}). \label{eq:concentration_AB} $$

The lattice Boltzmann model explicitly conserves mass during the reaction, which is essential for application to reactive flows, and Figure [fig:atob~r~eaction] compares the output of the numerical simulation with the analytical solution ([eq:concentration~A~B]). The combined concentration of reactant and product is seen in Figure [fig:atob~r~eaction:b] to deviate by no more than $3\times10^{-14}$ from $C_0$ during the simulation, which is constant within the limits of double-precision floating-point accuracy.


Second Order Reactions

Consider the simple second-order chemical reaction $A + B \rightarrow C,$ where the reaction rates for each species are defined by

$$ r_A = r_B = -C_aC_b k, \quad r_C=C_aC_b k. \label{eq:rate_ABC} $$

Here the reaction rate is limited by the species with the least concentration. The evolution of $C_a$ can be shown @Clark:1996 to follow

$$ \frac{C_a}{C_{a0} } = \frac{1-C_{b0}/C_{a0} }{(1-C_{b0}/C_{a0})\exp{[-(C_{a0}-C_{b0})kt}]}. \label{eqn:second_order_benchmark} $$

The system is solved using ([eq:lbm~r~d]) for each of the three species, using the reactive rates defined by ([eq:rate~A~BC]). The results of simulations with various initial ratios $C_a:C_b$, are shown in Figure [fig:reactions~s~econd~o~rder] to follow the analytical solution ([eqn:second~o~rder~b~enchmark]).


Third Order Oscillatory Reactions

To demonstrate the capability of the LB reaction model, and to validate the current implementation, an oscillatory Sel’kov reaction @selkov:1968 was simulated, following @PonceDawson:1993 - the nomenclature here follows that of the original. The reactions within this system are defined by

$$ {A \rightleftharpoons X, X}+2{Y} \rightleftharpoons 3 {Y, Y\rightleftharpoons B} $$

where the concentration of species, $C_X$ and $C_Y$

$$ \begin{aligned} R_1 &= k_1\rho_a - k_{-1}\rho_1 - k_2\rho_1\rho_2^2 + k_{-2}\rho_2^3\\ R_2 &= -k_3\rho_2 + k_{-3}\rho_b + k_2\rho_1\rho_2^2 - k_{-2}\rho_2^3\end{aligned} $$


 Rate     Value ($\Delta t^{-1}$)   ---------- -------------------------
$k_1$      $1.9006\times10^{-4}$    $k_{-1}$      $5\times10^{-4}$
$k_2$        $1\times10^{-2}$    $k_{-2}$      $1\times10^{-2}$    $k_{3}$       $5\times10^{-3}$    $k_{-3}$     $5.59\times10^{-6}$


image [fig:selkov]

The evolution of this system tends towards an oscillatory reaction where the concentrations of species $X$ and $Y$ reach a limit cycle, demonstrated in Figure [fig:selkov]. This simulation demonstrates that complex chemical reactions are easily modelled, and importantly that mass is conserved.

Parallelising the LB Algorithm

Observing the rate of chemical reaction during the approach to ergodic conditions requires domains of many obstruction lengths. Performing detailed computational pore-scale hydrodynamic and reactive advection-dispersion simulations of such systems requires high spatio-temporal resolution over the domains. The high number of degrees of freedom means that solving the system is computationally expensive, regardless of the numerical method chosen to solve it.

The lattice Boltzmann method inherits from lattice gas cellular automata a natural capacity for parallel computation which allows faster simulation of larger systems. During update from $t\rightarrow t+\Delta t$ only information from nearest-neighbour lattice nodes (i.e. $f_i$) is used for each of the lattice geometries used during the work reported here. This allows simple decomposition of the global computational domain into local domains which can be modelled concurrently with minimal intra-nodal communication. During the main temporal update loop, the particle densities, $f_i$, crossing nodal interfaces are the only data passed between computational nodes, and because of this, lattice Boltzmann belongs to the class of computer models described as ‘embarrasingly parallelizable’. Several reports on results from lattice Boltzmann simulations also describe their parallelisation strategy, for example @Kim:2003.

This study uses the Message Passing Interface, MPI @Snir:1998 [@Gropp:1999] to solve LB systems on distributed memory systems. Using this method, many individual computer processors operate in parallel on separate memory stores.

Multiple-data simulations are typically more difficult to implement than multiple processor-single data simulations, for which modern compilers require only certain lines of code to be marked as capable of being run in parallel, for example using OpenMP @Chapman:2007. However, the distributed memory systems on which MPI codes can run are significantly cheaper for the same computational ability than shared-memory equivalents. Because of this, on the computing cluster used for some of the simulations presented here, only 32 of the 256 processors were capable of operating with shared memory. It should be noted that significant changes in computer power and relative costs are certain to change any assumptions that are currently valid regarding optimum computer configurations.

The major steps required to parallelise the lattice Boltzmann algorithms are now defined. Due to the similarities between the hydrodynamic and reactive mass-transport lattice Boltzmann methods described elsewhere, the changes required for parallelisation are identical for all models, so no distinction is made here. Unless stated otherwise, ‘node’ in this section refers to a computational node rather than a lattice node as in other sections.

Domain Decomposition

A principal aim when developing distributed-memory solvers is to minimise the communication between processors, which normally limits the maximum rate of calculation. The maximum attainable speed-up is approximated by Amdahl’s law (originally suggested in @Amdahl:1967, but superceded for large data sets by Gustafson’s law @Gustafson:1988, this approximation is acceptable for the indication of expected speedup for the current work) as

$$ \mbox{ Max Speed Gain } = \frac{1}{s+(1-s)/n}, $$

where $s$ is the proportion of the calculation that is sequential, and $n$ the number of processors. Minimising the proportion of the simulation that must be performed sequentially is clearly of great importance to the efficiency of the overall parallel program.


The only data passed between processors in a parallel LB algorithm are the particle densities $f_i$ propagating to nodes stored on a remote computer. So in order to maximise the potential speedup running on a computing cluster, the proportion of nodes adjacent to computational boundaries should be minimised. The larger 2D domains considered in this work are channels much longer in the prevailing direction of flow than perpendicular to it, and for these, 1D domain decomposition is appropriate and efficient.

A typical decomposition is sketched in Figure fig:domain~d~ecomposition. Figures fig:domain~d~ecomposition and (c) highlight what happens at the interface between computational nodes, where propagation between the nodes takes $\Delta t$. The distribution vectors are those of the D2Q9 lattice - see Figures [fig:lbmvectors] and [fig:boundary~d~istributions].

1D decomposition of a 2D domain. [fig:domain~d~ecomposition]

Result Output and Convergence Checking

Results from all of the distributed-memory simulations were analysed on a single computer, and collecting data from each of the processors is a step that, with some methods, can limit the scalability of the algorithm. Lattice Boltzmann requires no global matrix assembly, but care is still required when writing results files.

Since all processors generally have access to a single filesystem, the results can be collected by having each process write results from its portion of the domain into a single file on that filesystem. In contrast to collecting all data on a single ‘master node’, this method scales better since domains are not limited in size to the largest that can be stored in the memory of a single node.

For large distributed problems, complete mass information for the entire domain can be too large a dataset to store and analyse at the temporal resolution required for concentration ‘outflow’ data. Therefore, in addition to saving distributed $u$ and $\rho$ to a single file at a single time, $\rho$ at a few locations (typically averaged over channel width to provide concentrations past certain obstacles) are stored at shorter intervals to a single file. The same ‘semaphore’ strategy is used to ensure that data is written in the correct order, with the result file appended to by each processor at each reporting timestep.

While planning the implementation, a very simple solver for a 1D lattice Boltzmann model of diffusion was written, using MPI to distribute the work between two computational nodes. The entire C++ source listing for that model can be contained in a single page and is included as Figure [fig:mpi~d~iffusion~c~ode] in Appendix [sec:code].

In summary, running the models required for this research in a distributed-memory cluster computer proved to be relatively straightforward.


This chapter has described the lattice Boltzmann models used to simulate physical systems throughout the current work, and has described where similar models have been used to perform related simulations in previous studies. Using some simple benchmarks, the models have been demonstrated to be accurate for hydrodynamic, mass transport and chemically reactive simulations.