I am often asked questions about a lattice Boltzmann hydrodynamic code that I published some time ago. I usually try to answer, but email isn't often the best way of describing the required theory. I've written that below, as clearly and concisely as I can.

It's a lightly-edited version of a dissertation introduction. If you have any suggestions that would improve this page for other readers, please get in touch. If you would like a pdf of the entire dissertation, drop me an email.

But first, a motivating animation:

Classical Boltzmann Equation

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 is far 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 enough to understand the fluid's properties. Any introduction to statistical mechanics covers this - for example, Chandler.

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$ \begin{equation} \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} \end{equation} could be derived from the Boltzmann equation Cercignani, 1987.

Lattice Boltzmann Equation

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 (Chen) to solve the incompressible Navier-Stokes equations.

Hydrodynamic Model

The Boltzmann equation, tracking the evolution of particle distributions $f_c$, with microscopic velocities $\xi_c$ is \begin{equation} \label{eq:boltzmann_equation} \frac{\partial f_c}{\partial t_c} + \xi_c \partial f_c = \Omega_c. \end{equation}

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 considers the collision as a linear relaxation from the non-equilibrium (propagated) particle distribution towards an equilibrium state $f^{eq}_c$. The time taken for $f_c$ to reach the equilibrium state $f_c^{eq}$ is termed the relaxation time $\lambda_c$. The BGK collision is therefore defined as: \begin{equation} \Omega_c = -\frac{1}{\lambda_c}(f_c-f^{eq}_c). \end{equation} 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$: \begin{equation} f^{eq}_c = \frac{\rho_c}{2\pi RT} \exp\left(\frac{-({{\bf \xi}_c - {\bf u}_c})^2}{2RT}\right). \label{eq:maxwell-boltzmann} \end{equation}

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 the Boltzmann equation above (\ref{eq:boltzmann_equation}) with the BGK collision operator in time, and neglecting terms of order ${\Delta t}^2$ or smaller, leads to the equation solved in the lattice Boltzmann method, \begin{equation} \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), \end{equation} where $\tau$ is a rescaled equivalent of the relaxation time $\lambda_c$. 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.

The equilibrium particle density $f^{eq}$, derived from the Maxwell-Boltzmann distribution above and retaining only terms up to order ${\bf u}^2$ during discretisation, is given by

\begin{equation} \label{eq:feq} 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), \end{equation}

where $w_i$, weightings from partition of unity, and the value of the constant $RT$ is determined by the lattice geometry. Because of the truncation, this equation is valid only for low Mach values $\frac{u}{RT} = Ma \ll 1\frac{\Delta x}{\Delta t}$, with a value of $u=0.15\frac{\Delta x}{\Delta t}$ regarded as a reasonable maximum; above this value, significant numerical inaccuracies can occur.

The above equations define the single phase hydrodynamic LBGK model. As with the continuous Boltzmann equation, hydrodynamic variables are calculated from the (microscopic velocity) moments of the particle distributions. The density $\rho$ and momentum $\rho{\bf u}$ are respectively: \begin{equation} \rho = \sum_{i=1}^{q}{f_i({\bf r},t)} \label{eq:lbm_density} \end{equation}

\begin{equation} \rho{\bf u} = \sum_{i=1}^{q} {\bf c}_i {f_i({\bf r},t)}. \label{eq:lbm_momentum} \end{equation}

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 \begin{equation} P=\frac{1}{3} \rho. \label{eq:lbm_pressure} \end{equation} The liquid's kinematic viscosity is given by \begin{equation} \nu = \frac{1}{3}\left( \frac{1}{\omega} - \frac{1}{2} \right) \frac{\Delta x^2}{\Delta t}, \end{equation} imposing a condition of $\omega<2$ for positive viscosity.

Boundary Conditions

Auxiliary conditions, including initial and boundary conditions, are required when simulating physical systems. It is then necessary to solve the update of $f_i$ in two steps; propagation \begin{equation} \tilde{f}(r, t + \Delta t) = f_i({\bf r} - {\bf \xi}_i \Delta t, t), \label{eq:propagation} \end{equation} and collision

\begin{equation} \label{eq:collision} f_i({\bf r} + {\bf \xi}_i, t + \Delta t) = \left(1 - \omega \right) \tilde{f}({\bf r}, t) + \omega f_i^{eq}(\rho({\bf r},t), {\bf u}({\bf r}, t)). \end{equation}

The new parameter $\tilde{f}$ 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 steps.

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 equation \ref{eq:feq}.
• Calculate equilibrium. Set $f^{eq}(\rho,{\bf u})$, also calculated using equation \ref{eq:feq}.
• Collision. Update particle distributions on each node, as a function of those from the previous state and the equivalent equilibrium state, using equation \ref{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 (equation \ref{eq:propagation}).
• Completed timestep. Set $t \leftarrow t+\Delta t$.
• State parameters. Calculate $\rho(t,{\bf r})$ and $\rho{\bf u}(t,{\bf r})$ using equations \ref{eq:lbm_density} and \ref{eq:lbm_momentum}.
• Simulation complete? If $t<t_{max}$ go back to **Calculate equilibrium**; 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.

Looking at the algorithm and the equations used in it, you can see that the method is computationally simple, containing only a relatively small number of inexpensive floating-point arithmetic operations to update. This is in contrast to other numerical methods which solve the Navier-Stokes equations using a system of linear equations whose solution is expensive (for example, by Gaussian elimination and back-substitution). In addition, locality makes the method relatively simple to parallelize.

It might help with understanding to notice 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. Under this constraint, the diffusive LBGK model corresponds to a Lax-Wendroff scheme. But no-one uses LBGK like this because it is more efficient to retain the flexibility in $\Delta t$ and $\Delta x$ allowed by $\tau\neq 1$.

Code

Complete self-contained 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 available at exolete.com.

(to be continued … I'm converting from the original LaTeX and simplifying as I have time)

Bibliography

• P. Bhatnagar, E. Gross and M. Krook. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94 (1954), 511–525. A breakthrough simplification that was found to be very useful some decades later. If you're going to study this field, you should certainly read it, but it makes a poor starting point if you're being practical. ↩︎
• C. Cercignani. The Boltzmann Equation and its Applications. Springer-Verlag, 1987. See also S. Chapman and T.G. Cowling The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases. Springer, 1990. The concluding chapter gives a historical summary of the progress towards an accepted microscopic model, due to Enskog, for the transport of dense gases and liquids. ↩︎
• C. Cercignani. Ludwig Boltzmann: The man who Trusted Atoms. Oxford University Press, 1998. It's a really interesting book, and certainly the most accessible in this list. ↩︎
• D. Chandler. Introduction to Modern Statistical Mechanics. Oxford University Press, 1987. There are many like it, but this one is good. ↩︎
• H. Chen, S. Chen and W.H. Matthaeus. Recovery of the Navier–Stokes equations using a lattice–gas Boltzmann method. Phys. Rev. A 45, 8 (1992). ↩︎

⁂