Lattice Boltzmann Method — fluid dynamics from kinetic theory

1. Kinetic theory — the Boltzmann equation

Classical fluid dynamics (Navier-Stokes) works with macroscopic fields: density ρ, velocity u and pressure p. Kinetic theory takes a mesoscopic view: the fluid is described by a distribution function f(x, ξ, t) — the density of particles at position x with microscopic velocity ξ at time t.

The continuous Boltzmann transport equation governs the evolution of f:

The key insight of LBM is to discretise both space and velocity. Space becomes a regular lattice; the continuous velocity space collapses to a small finite set of allowed velocities e_i. The collision term is approximated by the simple BGK model. Despite these drastic simplifications, the Chapman-Enskog expansion shows that the resulting equations reproduce the incompressible Navier-Stokes equations to second order in the Mach number.

2. D2Q9 lattice — velocity set and weights

The name "D2Q9" encodes the dimensionality and velocity count: 2 dimensions, 9 velocity directions. Every node on the lattice propagates density along 9 directions (including rest).

The weights w_i are chosen so that the lattice correctly recovers the first and second moments (mass and momentum) of the equilibrium distribution. The constraint Σ w_i = 1 ensures mass conservation.

3. BGK collision operator

The Bhatnagar-Gross-Krook (BGK) approximation replaces the true collision integral with a relaxation toward local equilibrium. Each distribution function f_i relaxes at rate 1/τ:

The relaxation time τ controls viscosity: τ = 1 gives ν = 1/6 (moderately viscous). Stability requires τ > 0.5. As τ → 0.5, viscosity → 0 and the simulation becomes unstable for high-velocity flows.

4. Streaming step

After collision, each post-collision distribution f_i* propagates to the neighbouring lattice node in direction e_i:

The full LBM update cycle per time step is therefore:

1. Compute macroscopic ρ and u from the current f
2. Compute equilibrium distributions f^eq
3. Apply BGK collision: f* = f − (1/τ)(f − f^eq)
4. Stream: shift each f_i* one step in direction e_i
5. Apply boundary conditions

5. Macroscopic recovery — ρ, u from f

The macroscopic fluid quantities are moments of the distribution function. Mass density ρ is the zeroth moment; momentum ρu is the first:

This local computation — requiring only the 9 f_i at a single node — is what makes LBM highly parallelisable. Each node is independent during the collision step; only the streaming step requires neighbouring data.

6. Boundary conditions

No-slip walls — Bounce-Back

The simplest no-slip boundary: any distribution arriving at a solid wall node is reflected back in the opposite direction — as if the particle bounced off the wall.

Moving walls — Zou-He

For inflow/outflow with known velocity u or known pressure ρ, the Zou-He scheme extrapolates the unknown distributions analytically. For a left inlet with velocity u = (u_x, 0):

7. Full JavaScript D2Q9 implementation

// ── D2Q9 constants ────────────────────────────────────────────────
const NX = 200, NY = 100;  // grid size
const TAU = 0.6;            // relaxation time  (ν = (τ−0.5)/3 ≈ 0.033)
const OMEGA = 1 / TAU;      // relaxation rate
const N = NX * NY;

// Velocity vectors: 0=rest, 1–4 axis, 5–8 diagonal
const EX = [0,  1,  0, -1,  0,  1, -1, -1,  1];
const EY = [0,  0,  1,  0, -1,  1,  1, -1, -1];
const W  = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36];
const OPP = [0, 3, 4, 1, 2, 7, 8, 5, 6]; // opposite direction index

// Distribution functions: f[i * N + x + y*NX]
const f    = new Float32Array(9 * N).fill(0);
const fTmp = new Float32Array(9 * N);
const solid = new Uint8Array(N); // 1 = obstacle

function idx(x, y) { return x + y * NX; }
function fi(i, x, y) { return i * N + idx(x, y); }

// ── Initialise: uniform flow u_x = u0 ─────────────────────────────
function init(u0 = 0.1) {
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      for (let i = 0; i < 9; i++) {
        f[fi(i, x, y)] = feq(i, 1.0, u0, 0.0);
      }
    }
  }
  // Add a circular obstacle near the left
  const cx = NX / 4, cy = NY / 2, r = NY / 8;
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      if ((x-cx)**2 + (y-cy)**2 < r*r) solid[idx(x,y)] = 1;
    }
  }
}

// ── Equilibrium distribution ──────────────────────────────────────
function feq(i, rho, ux, uy) {
  const eu  = EX[i] * ux + EY[i] * uy;
  const u2  = ux * ux + uy * uy;
  return W[i] * rho * (1 + 3*eu + 4.5*eu*eu - 1.5*u2);
}

// ── Single LBM time step ──────────────────────────────────────────
function step(u0 = 0.1) {
  // 1. Collision + stream (pull scheme)
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      if (solid[idx(x, y)]) continue;

      // Compute macroscopic quantities
      let rho = 0, ux = 0, uy = 0;
      for (let i = 0; i < 9; i++) {
        const fi_val = f[fi(i, x, y)];
        rho += fi_val;
        ux  += fi_val * EX[i];
        uy  += fi_val * EY[i];
      }
      ux /= rho; uy /= rho;

      // BGK collision, pull stream
      for (let i = 0; i < 9; i++) {
        // Source node for streaming
        const sx = x - EX[i], sy = y - EY[i];
        if (sx < 0 || sx >= NX || sy < 0 || sy >= NY) {
          fTmp[fi(i, x, y)] = feq(i, 1.0, u0, 0.0); // inflow
          continue;
        }
        if (solid[idx(sx, sy)]) {
          // Bounce-back from obstacle
          fTmp[fi(i, x, y)] = f[fi(OPP[i], x, y)];
          continue;
        }
        const f_src = f[fi(i, sx, sy)];
        fTmp[fi(i, x, y)] = f_src - OMEGA * (f_src - feq(i, rho, ux, uy));
      }
    }
  }
  fTmp.copyTo ? fTmp.copyTo(f) : f.set(fTmp);
}

// ── Render: colour by velocity magnitude ─────────────────────────
function render(ctx, canvas) {
  const img = ctx.createImageData(NX, NY);
  const d   = img.data;
  for (let y = 0; y < NY; y++) {
    for (let x = 0; x < NX; x++) {
      const p = idx(x, y);
      if (solid[p]) { d[p*4]   = 200; d[p*4+1] = 200; d[p*4+2] = 200; d[p*4+3] = 255; continue; }
      let ux = 0, uy = 0, rho = 0;
      for (let i = 0; i < 9; i++) {
        const v = f[fi(i, x, y)]; rho += v; ux += v*EX[i]; uy += v*EY[i];
      }
      ux /= rho; uy /= rho;
      const spd = Math.hypot(ux, uy) / 0.25; // normalise to [0,1]
      const t   = Math.min(spd, 1);
      d[p*4]   = (t * 56)  | 0;  // r
      d[p*4+1] = (t * 189) | 0;  // g (→ cyan)
      d[p*4+2] = (t * 248) | 0;  // b
      d[p*4+3] = 255;
    }
  }
  ctx.putImageData(img, 0, 0);
}

8. Applications and extensions

• Kármán vortex street: the classic flow-past-cylinder benchmark — at Re ≈ 100 a periodic vortex shedding pattern develops behind the obstacle. LBM reproduces it naturally.
• D3Q19 / D3Q27: 3D lattice variants with 19 or 27 velocities. Identical structure — just more velocity directions and correspondingly higher memory cost.
• Multi-relaxation time (MRT): replaces the single τ with a relaxation matrix in moment space. Improves stability at low viscosity (high Re) and reduces non-physical shear artefacts.
• Multiphase LBM: the Shan-Chen model adds an inter-particle force term to the collision step, producing phase separation, surface tension and droplet dynamics — entirely within the LBM framework.
• Thermal LBM: a second distribution function tracks energy, enabling buoyancy-driven convection (Rayleigh-Bénard instability).