Theory & Methodology

A self-contained treatment of the Lattice Boltzmann Method: from the kinetic theory of gases through the D2Q9 discretization, BGK collision operator, boundary conditions, and the Chapman-Enskog connection to the Navier-Stokes equations.

1. From Boltzmann to Lattice Boltzmann

The Boltzmann equation describes the evolution of a particle distribution function \(f(\mathbf{x}, \mathbf{c}, t)\) in phase space:

$$ \frac{\partial f}{\partial t} + \mathbf{c} \cdot \nabla f = \Omega(f) $$
Eq. 1: Boltzmann Equation

In the Lattice Boltzmann Method, we discretize this equation in velocity space using a finite set of velocity vectors \(\mathbf{c}_i\), physical space using a regular Cartesian grid, and time using a simple explicit scheme. The collision operator \(\Omega(f)\) is replaced by the Bhatnagar-Gross-Krook (BGK) approximation:

$$ \Omega_i = -\frac{1}{\tau} \left( f_i - f_i^{\text{eq}} \right) $$
Eq. 2: BGK Collision Operator

where \(\tau\) is the dimensionless relaxation time and \(f_i^{\text{eq}}\) is the equilibrium distribution. This gives the discrete lattice Boltzmann equation:

$$ f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f_i(\mathbf{x}, t) - \frac{1}{\tau} \left( f_i(\mathbf{x}, t) - f_i^{\text{eq}}(\mathbf{x}, t) \right) $$
Eq. 3: Discrete Lattice Boltzmann Equation (BGK)

This equation forms the heart of the solver. The left side represents the streaming step (distributions propagate to neighboring lattice sites), and the right side represents the collision step (distributions relax toward equilibrium).

Algorithm Flow. Each time step executes: (1) Apply boundary conditions, (2) Collide: relax distributions toward equilibrium, (3) Stream: propagate distributions to neighbors, (4) Extract macroscopic properties and forces.

2. The D2Q9 Lattice

For 2D simulations, the D2Q9 lattice uses 9 discrete velocity vectors: one rest particle, four axial directions, and four diagonal directions.

$$ \mathbf{c}_i = \begin{cases} (0, 0) & i = 0 \\ (\pm 1, 0), (0, \pm 1) & i = 1,2,3,4 \\ (\pm 1, \pm 1) & i = 5,6,7,8 \end{cases} $$
Eq. 4: D2Q9 Velocity Set

Velocity vector layout on the lattice:

    6 ( -1,  1)    2 ( 0,  1)    5 ( 1,  1)
    3 ( -1,  0)    0 ( 0,  0)    1 ( 1,  0)
    7 ( -1, -1)    4 ( 0, -1)    8 ( 1, -1)

The quadrature weights \(w_i\) are chosen to satisfy the moment isotropy conditions required to recover the Navier-Stokes equations:

$$ w_i = \begin{cases} 4/9 & i = 0 \quad \text{(rest)} \\ 1/9 & i = 1,2,3,4 \quad \text{(axial)} \\ 1/36 & i = 5,6,7,8 \quad \text{(diagonal)} \end{cases} $$
Eq. 5: D2Q9 Quadrature Weights

The lattice speed of sound for D2Q9 is \(c_s = 1/\sqrt{3}\) in lattice units, which follows from the velocity moment conditions.

3. Equilibrium Distribution

The equilibrium distribution is a second-order expansion of the Maxwell-Boltzmann distribution, discretized onto the D2Q9 lattice:

$$ f_i^{\text{eq}} = w_i \rho \left[ 1 + 3 (\mathbf{c}_i \cdot \mathbf{u}) + \frac{9}{2} (\mathbf{c}_i \cdot \mathbf{u})^2 - \frac{3}{2} (\mathbf{u} \cdot \mathbf{u}) \right] $$
Eq. 6: D2Q9 Equilibrium Distribution

This expansion is valid in the low Mach number limit (\(M = |\mathbf{u}|/c_s \ll 1\)), which constrains the inlet velocity to approximately \(u_{\text{inflow}} \lesssim 0.3\) in lattice units.

The macroscopic density \(\rho\) and velocity \(\mathbf{u}\) are recovered as moments of the distribution functions:

$$ \rho = \sum_{i=0}^{8} f_i, \qquad \rho \mathbf{u} = \sum_{i=0}^{8} f_i \mathbf{c}_i $$
Eq. 7: Macroscopic Properties from Distribution Moments

4. Chapman-Enskog: Connecting LBM to Navier-Stokes

The macroscopic Navier-Stokes equations emerge from the lattice Boltzmann equation through the Chapman-Enskog multiscale expansion. The kinematic viscosity \(\nu\) is related to the relaxation time \(\tau\):

$$ \nu = c_s^2 \left( \tau - \frac{1}{2} \right) \Delta t = \frac{1}{3} \left( \tau - \frac{1}{2} \right) $$
Eq. 8: Viscosity from Relaxation Time (lattice units with \(\Delta t = 1\), \(\Delta x = 1\))

The Reynolds number is therefore controlled entirely by \(\tau\):

$$ \text{Re} = \frac{u_{\text{inflow}} N_x}{\nu} = \frac{3 u_{\text{inflow}} N_x}{\tau - 1/2} $$
Eq. 9: Reynolds Number in Lattice Units

The recovered macroscopic equations are the incompressible Navier-Stokes equations with an error of \(O(M^2)\):

$$ \nabla \cdot \mathbf{u} = 0 $$ $$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} $$
Eq. 10: Recovered Navier-Stokes Equations (low Mach limit)

The pressure is related to density through the equation of state \(p = \rho c_s^2\), making LBM a weakly compressible solver.

5. Boundary Conditions

5.1 Zou/He Velocity Inlet

At the left boundary (\(x = 0\)), the velocity is enforced to \((u, v) = (u_{\text{inflow}}, 0)\). The unknown distribution functions are computed from the known ones using the Zou/He method, which applies the bounce-back rule to the non-equilibrium part of the distributions:

$$ \rho_{\text{inlet}} = \frac{1}{1 - u_{\text{inflow}}} \left[ f_0 + f_2 + f_4 + 2(f_3 + f_6 + f_7) \right] $$
Eq. 11: Inlet Density from Zou/He

All unknown inlet distributions are then set to their equilibrium values using \(\rho_{\text{inlet}}\) and \((u_{\text{inflow}}, 0)\).

5.2 Convective Outlet

At the right boundary (\(x = N_x - 1\)), a zero-gradient (convective) condition is applied:

$$ f_i(N_x - 1, y) = f_i(N_x - 2, y) \quad \text{for all } i $$
Eq. 12: Zero-Gradient Outlet

5.3 Bounce-Back (No-Slip Walls)

Solid boundaries (cylinder surface, top/bottom walls) use the bounce-back scheme. Distributions that stream into a solid node are reflected back along their incoming direction:

$$ f_{\bar{i}}(\mathbf{x}_f, t + \Delta t) = f_i^*(\mathbf{x}_f, t) $$
Eq. 13: Bounce-Back where \(\bar{i}\) is the direction opposite to \(i\)

This enforces the no-slip condition at the wall, located halfway between the fluid and solid nodes.

5.4 Periodic Boundaries

The top and bottom boundaries (\(y = 0\) and \(y = N_y - 1\)) are periodic, meaning fluid that exits the top re-enters at the bottom and vice versa:

$$ f_i(x, 0) = f_i(x, N_y - 1), \quad f_i(x, N_y - 1) = f_i(x, 0) $$
Eq. 14: Periodic Y-Boundaries

6. Force Extraction via Momentum Exchange

The forces acting on the cylinder are computed using the momentum exchange method. For each boundary link connecting a fluid node \(\mathbf{x}_f\) to an obstacle node \(\mathbf{x}_b\), the momentum transfer is:

$$ \Delta \mathbf{m} = \mathbf{c}_i \left[ f_i(\mathbf{x}_f, t) - f_{\bar{i}}(\mathbf{x}_b, t) \right] $$
Eq. 15: Momentum Exchange on a Boundary Link

Summing over all boundary links gives the total force on the cylinder. The drag and lift coefficients follow:

$$ C_d = \frac{2 F_x}{\rho u_{\text{inflow}}^2 D}, \qquad C_l = \frac{2 F_y}{\rho u_{\text{inflow}}^2 D} $$
Eq. 16: Drag and Lift Coefficients

where \(D = 2R\) is the cylinder diameter.

7. Memory Architecture

High-performance LBM requires careful memory layout. Our solver uses a flat 1D array for the distribution functions, indexed as:

$$ \text{index} = (y \cdot N_x + x) \cdot 9 + i $$
Eq. 17: Flat 1D Array Indexing for \(f_i(x,y)\)

This layout ensures that the 9 distributions at a single node are stored contiguously in memory, maximizing cache utilization during the collision step. Nested vectors (\(\text{vector}[\text{vector}]\)) would scatter the same data across non-contiguous heap allocations, causing cache misses.

8. Parameter Summary

SymbolValueUnitsDescription
\(N_x\)400cellsGrid width (streamwise)
\(N_y\)150cellsGrid height (wall-normal)
\(c_s^2\)\(1/3\)--Lattice speed of sound squared
\(u_{\text{inflow}}\)0.1lu ts\(^{-1}\)Inlet velocity (Mach ~ 0.17)
\(D\)30cellsCylinder diameter (\(N_y/5\))
\(\tau\)\(0.5 + 3\nu\)tsRelaxation time (Re-dependent)
\(\nu\)\(u_{\text{inflow}} N_x / \text{Re}\)lu\(^2\) ts\(^{-1}\)Kinematic viscosity