a Kármán vortex street (the alternating swirls that form behind any blunt object in fast-moving flow). You’ve seen this pattern trailing behind bridge pillars, airplane wings, and chimney stacks. It’s one of the most recognizable phenomena in fluid dynamics.
I generated it without solving a single fluid equation.
Instead, this Lattice Boltzmann Method (LBM) simulation treats fluid as what it actually is at the microscopic level: a statistical distribution of particles colliding with each other. The macroscopic behavior (vortices, pressure gradients, viscous damping) falls out as a consequence, not an input.
LBM occupies a genuinely unusual position in computational physics, halfway between a molecular dynamics simulation and a Navier-Stokes solver, operating at the mesoscopic scale. In this article, I’ll derive it from first principles, implement it in ~200 lines of C++, and run it on the MareNostrum 5 supercomputer. The full source and a simplified Python/NumPy version are on GitHub.
Quick take: LBM replaces the intractable Navier-Stokes PDE with a much simpler update rule on a discrete velocity lattice. The NS equations are recovered in the correct limit, not assumed. This makes LBM theoretically elegant, practically fast, and embarrassingly parallelizable. It works best in the low-Mach, subsonic regime, and becomes numerically unstable at very low viscosities.
Why the obvious approach is painful
The standard route to fluid simulation goes through the Navier-Stokes equations, the macroscopic conservation laws that govern virtually everything we care about in fluid dynamics. They’re analytically intractable except in a handful of toy geometries, so in practice you discretize them on a mesh and solve numerically: Finite Volume, Finite Difference, or Finite Element methods.
These work well, but they carry a hidden cost that compounds for complex geometries. For a curved boundary (a cylinder, a car body, a porous medium) you either need a body-fitted mesh that’s expensive to generate and hard to automate, or specialized stencil corrections near the wall. Every time the geometry changes, the mesh has to be regenerated. It’s an engineering problem layered on top of the physics problem.
LBM sidesteps this almost entirely. The reason becomes clear once you understand where it starts.
A different starting point: statistical mechanics
Instead of asking “what is the velocity field at this point?”, LBM asks “how many particles are moving in each direction at this point?”
The central object is the single-particle distribution function f(x, v, t). It tells you the expected number of particles near position x moving with velocity v at time t. The macroscopic fields we actually care about are just moments of f:
\[
\rho(\mathbf ysrc >= Ny);
f_new[x][y][i] = (oob ,t) = m\int f\, d^3v \qquad \rho\,\mathbf
const int xsrc = x – EIX[i];
const int ysrc = y – EIY[i];
const bool oob = (xsrc < 0 (\mathbf,t) = m\int \mathbf
const int xsrc = x – EIX[i];
const int ysrc = y – EIY[i];
const bool oob = (xsrc < 0 \,f\, d^3v
\]
At global equilibrium (uniform density, no bulk flow) f takes the familiar Maxwell-Boltzmann form:
\[
f^{eq}(\mathbf{v}) = \frac{\rho}{(2\pi RT)^{D/2}} \exp\!\left(-\frac{|\mathbf{v} – \mathbf{u}|^2}{2RT}\right)
\]
Away from equilibrium, f takes some other shape. The equation governing how it evolves is the Boltzmann equation:
\[
\frac{\partial f}{\partial t} + \mathbf{v}\cdot\nabla_\mathbf{x} f = \Omega[f]
\]
The left-hand side is pure transport: particles carry their velocity distribution as they stream through space. The right-hand side, Ω[f], is the collision operator that drives f back toward f^eq by redistributing particles among velocity classes. The problem is that Ω[f] is a five-dimensional integral over all possible collision partners and scattering angles, encoding the full microscopic collision physics. For any practical simulation, it’s completely intractable.
The one-line miracle: BGK
In 1954, Bhatnagar, Gross and Krook proposed a brilliantly blunt simplification. Instead of modelling every collision in detail, just capture the net effect: collisions drive f toward equilibrium at a rate proportional to how far from equilibrium it currently sits.
\[
\Omega[f] \approx -\frac{1}{\tau}\bigl(f – f^{eq}\bigr)
\]
This is the BGK approximation, and it reduces an intractable five-dimensional integral to a scalar multiplication. The parameter τ is the relaxation time, meaning the average time the system takes to return to local equilibrium after a perturbation. Large τ means slow relaxation; small τ means fast.
What makes this more than a convenient hack is what τ means physically. Via the Chapman-Enskog expansion (a systematic perturbative analysis) you can prove that in the limit of small Knudsen number (many collisions per mean free path), the BGK equation recovers the Navier-Stokes equations exactly, with kinematic viscosity:
\[
\nu = c_s^2\!\left(\tau – \frac{1}{2}\right)
\]
Note that this relation operates in dimensionless lattice units; converting to physical units requires a separate length and timescale calibration. But within the simulation, τ is the entire physical knob: turn it up for thick, viscous flow; push it toward 0.5 for high-Reynolds-number flow (though values below ~0.55 risk numerical instability).
Crucially, this is why LBM is more than Navier-Stokes in disguise. The connection to NS is derived, not assumed: you’re solving a kinetic equation that provably converges to NS behavior in the appropriate limit. This is a very common theme in all of statistical mechanics.
From continuous to discrete: the D2Q9 lattice
The BGK equation still lives in continuous velocity space. To simulate it we need a finite set of velocities. The approach is a Taylor expansion of f^eq in powers of the bulk velocity u, valid when u is small compared to the lattice speed of sound (the low-Mach assumption):
\[
f_i^{eq} = \rho\, w_i\!\left[1 + 3(\mathbf{c}_i\cdot\mathbf{u}) + \frac{9}{2}(\mathbf{c}_i\cdot\mathbf{u})^2 – \frac{3}{2}|\mathbf{u}|^2\right]
\]
The continuous velocity c is replaced by a finite set of discrete velocities c_i, each carrying a quadrature weight w_i chosen to correctly recover the velocity moments. For 2D simulations the standard choice is D2Q9: 2 dimensions, 9 velocity directions: rest, four axis-aligned directions, and four diagonals.
The weights are fixed by isotropy: 4/9 for the rest particle, 1/9 for axis directions, 1/36 for diagonals. The lattice speed of sound is cs = 1/√3. Everything else in the algorithm follows from these nine numbers.
The algorithm: two steps, endlessly repeated
With the discrete equilibrium in hand, the full LBM algorithm reduces to four operations per timestep, run in a loop until convergence:
- Macroscopic recovery: compute ρ and u at each node from the current populations f_i
- Collision: relax each f_i toward f_i^eq at rate 1/τ (purely local, no neighbor communication)
- Streaming: propagate each population to the neighboring node in its direction of travel
- Boundary conditions: reconstruct missing populations at domain edges and solid walls
The collision step is where all the physics happens. The streaming step is purely kinematic bookkeeping. Their alternation is what produces fluid behavior.
void collide() {
const double omega = 1.0 / tau;
#pragma omp parallel for collapse(2) schedule(static)
for (int x = 0; x < Nx; x++) {
for (int y = 0; y < Ny; y++) {
if (solid[x][y]) continue;
double feq[9];
computeEquilibrium(rho[x][y], ux[x][y], uy[x][y], feq);
for (int i = 0; i < 9; i++)
f[x][y][i] -= omega * (f[x][y][i] - feq[i]);
}
}
}The streaming implementation turns out to matter enormously for performance. The naïve push pattern (where each thread at (x,y) writes outgoing populations into neighbouring nodes) causes false sharing: threads write to adjacent memory regions, hammering the same cache lines and forcing constant invalidation. The pull (gather) pattern below fixes this, because each thread reads from its upstream neighbours and only writes to its own node.
void stream() {
#pragma omp parallel for collapse(2) schedule(static)
for (int x = 0; x < Nx; x++) {
for (int y = 0; y < Ny; y++) {
if (solid[x][y]) continue;
for (int i = 0; i < 9; i++) {
const int xsrc = x - EIX[i];
const int ysrc = y - EIY[i];
const bool oob = (xsrc < 0 || xsrc >= Nx ||
ysrc < 0 || ysrc >= Ny);
f_new[x][y][i] = (oob || solid[xsrc][ysrc])
? f[x][y][OPP[i]] // bounce-back from wall
: f[xsrc][ysrc][i]; // stream from upstream neighbour
}
}
}
}Because collision is entirely local and streaming only requires nearest-neighbour reads, the whole algorithm is embarrassingly parallel: every node updates independently, with no global solves, no implicit systems, and no matrix inversions. This maps cleanly onto multi-core CPUs, GPUs, and distributed clusters; LBM implementations on CUDA routinely achieve near-linear scaling with thread count up to the memory bandwidth limit.
Boundary conditions: complex geometry for free
After each streaming step, nodes adjacent to solid walls are missing some incoming populations, as they would have arrived from inside the solid, which has no fluid to stream. How you fill those missing populations defines the boundary condition.
Bounce-back handles no-slip walls. Any population that streams into a solid node is reflected back in the opposite direction: f_ī = f*_i. The physical picture is a particle reversing its velocity at the wall, and the no-slip condition (zero velocity at the boundary) emerges automatically, with no explicit velocity constraint needed. The physical wall sits halfway between the last fluid node and the first solid node, giving the scheme second-order spatial accuracy.
The elegant consequence is that any obstacle shape is just a boolean flag. To add a cylinder, mark the relevant nodes as solid. The bounce-back rule handles the boundary physics at every flagged node automatically, with no mesh, no stencil corrections, and no code changes between geometries.
At the inlet we use the Zou-He scheme: a fixed velocity uw is imposed, and the three missing rightward-pointing populations (f_1, f_5, f_8) are reconstructed consistently from the remaining known populations and the imposed velocity, recovering the correct inlet density without imposing it independently.
At the outlet, a zero-gradient condition copies all populations from the penultimate column, allowing fluid to leave the domain freely without reflecting spurious waves back upstream.
Results: two regimes
The simulation domain is a 400×400 lattice with a circular obstacle of radius 40 lattice units, centred at one quarter of the domain length and offset slightly vertically to break symmetry. The inlet velocity is uw = 0.1 lattice units, well below cs ≈ 0.577 (low-Mach assumption satisfied throughout). Two values of τ reveal qualitatively different flow regimes.
At τ = 0.75 (Re ≈ 96) the wake is steady and symmetric. The distribution function has relaxed close to local equilibrium everywhere, as collisions are vigorous enough to damp any perturbation before it can grow. The flow is viscous and laminar.
At τ = 0.55 (Re ≈ 480) the picture changes completely.
Vortices shed alternately from the upper and lower surfaces, forming the Kármán street. From the statistical mechanics perspective this is the key insight: the vortex shedding is driven by the departure of f from its Maxwell-Boltzmann equilibrium shape. The smaller τ means collisions relax f toward equilibrium more slowly, allowing the non-equilibrium component f − f^eq to grow large enough to feed back into the macroscopic momentum field and sustain the oscillation. A perfectly equilibrium fluid would show no such structure.
The boolean-flag approach means changing geometry costs nothing. Here is the same solver, unmodified, applied to three square obstacles in a row:
The wakes couple, the vortex streets interact, and rich multi-body flow structure emerges. The only change to the code was marking different nodes as solid.
Scaling to MareNostrum 5
With #pragma omp parallel for on both the collision and streaming loops, the solver scales to multiple cores with minimal effort. The important choice is the streaming pattern: pull streaming ensures each thread exclusively owns its output memory, eliminating the false sharing that causes the naïve push implementation to plateau early.
We benchmarked on a full node of MareNostrum 5 at the Barcelona Supercomputing Center: 2× Intel Xeon Platinum 8480+ (Sapphire Rapids), 112 cores across two NUMA domains.
Three regimes emerge cleanly. From 1 to 32 threads, efficiency stays above 75%, reaching a 24.6× speedup, which is genuinely strong scaling for a memory-heavy stencil code. At 64 threads, efficiency drops to 42% as the memory controller saturates (the ~30 MB of distribution function data is being read and written faster than the bus can sustain). At 112 threads, the speedup actually falls to 23.5×, because spanning both physical sockets means half of all memory accesses cross the inter-socket interconnect, adding latency on top of already saturated bandwidth. The practical sweet spot for a 400×400 grid is 32 cores on a single socket.
One number worth pausing on: the single-thread baseline after the pull-streaming rewrite is 31.7 seconds, down from 113 seconds with the original push implementation. That 3.6× improvement came from a code change alone, on a single core, before touching any parallelism. The memory access pattern matters long before the hardware does.
When to use LBM (and when not to)
LBM is the right tool when geometry is complex or changes over time, when you need to run on massively parallel hardware, or when you’re modelling mesoscale phenomena like velocity slip near walls that continuum solvers need empirical corrections for. It’s also an excellent foundation for multiphase simulations, the subject of the next article in this series.
It is not the right tool for high-Mach flows (the low-order Taylor expansion breaks down above Mach ~0.3), very high Reynolds numbers (the τ > 0.5 stability requirement limits how low viscosity can go without more advanced collision operators), or situations where precise physical unit calibration is critical from the outset.
What’s next: multiphase flows with Shan-Chen
The solver described here handles a single-component, single-phase fluid. Many of the most physically interesting problems involve two immiscible fluids (oil and water, rising bubbles, droplet dynamics in airflow).
The Shan-Chen model extends LBM to handle these cases by adding short-range inter-particle forces that produce an effective equation of state with a phase transition. Interfacial tension and phase separation emerge from the collision statistics rather than being imposed as boundary conditions. In the next article, I’ll derive the Shan-Chen model, implement it as a modest extension of this solver, and show what spontaneous phase separation looks like when you let two fluids compete for the same lattice.
Closing thought
What I find most satisfying about LBM is what it reveals about the relationship between scales. The Navier-Stokes equations feel fundamental: they’re the equations of fluid dynamics, written at the scale we can see and touch. But they’re not fundamental at all. They’re a statistical consequence of many particles colliding, recoverable from a kinetic description in the appropriate limit.
LBM makes that relationship concrete and computational. Viscosity isn’t a parameter you set; it’s a collision timescale. The Kármán vortex street isn’t a fluid instability you program; it’s what emerges when collisions can’t keep up with inertia. Every macroscopic phenomenon in the simulation has a direct microscopic interpretation in the statistical mechanics underneath it.
That bridge between scales, made tangible in ~200 lines of C++, is the whole point.
References:
For those who want to go deeper, the following are the canonical references in roughly increasing order of technical difficulty.
Mohamad, A.A. (2019). Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes. Springer. The most accessible entry point; worked examples and code throughout.
Krüger, T. et al. (2017). The Lattice Boltzmann Method: Principles and Practice. Springer. The comprehensive reference; derivations, boundary conditions, multiphase models, and GPU implementation all covered.
Bhatnagar, P.L., Gross, E.P. & Krook, M. (1954). A Model for Collision Processes in Gases. Physical Review, 94(3), 511–525. The original BGK paper.
Chapman, S. & Cowling, T.G. (1970). The Mathematical Theory of Non-Uniform Gases. Cambridge University Press. For the Chapman-Enskog expansion in full rigour.
The full C++ source, OpenMP benchmark scripts, and Python/NumPy reference implementation for this article are available on GitHub.