a Kármán vortex avenue (the alternating swirls that kind behind any blunt object in fast-moving move). You’ve seen this sample trailing behind bridge pillars, airplane wings, and chimney stacks. It’s one of the crucial recognizable phenomena in fluid dynamics.
I generated it with out fixing a single fluid equation.
As an alternative, this Lattice Boltzmann Technique (LBM) simulation treats fluid as what it truly is on the microscopic stage: a statistical distribution of particles colliding with one another. The macroscopic habits (vortices, strain gradients, viscous damping) falls out as a consequence, not an enter.
LBM occupies a genuinely uncommon place in computational physics, midway between a molecular dynamics simulation and a Navier-Stokes solver, working on the mesoscopic scale. On this article, I’ll derive it from first ideas, implement it in ~200 strains of C++, and run it on the MareNostrum 5 supercomputer. The complete supply and a simplified Python/NumPy model are on GitHub.
Fast take: LBM replaces the intractable Navier-Stokes PDE with a a lot less complicated replace rule on a discrete velocity lattice. The NS equations are recovered within the right restrict, not assumed. This makes LBM theoretically elegant, virtually quick, and embarrassingly parallelizable. It really works finest within the low-Mach, subsonic regime, and turns into numerically unstable at very low viscosities.
Why the apparent strategy is painful
The usual path to fluid simulation goes via the Navier-Stokes equations, the macroscopic conservation legal guidelines that govern just about every thing we care about in fluid dynamics. They’re analytically intractable besides in a handful of toy geometries, so in apply you discretize them on a mesh and resolve numerically: Finite Quantity, Finite Distinction, or Finite Factor strategies.
These work properly, however they carry a hidden price that compounds for complicated geometries. For a curved boundary (a cylinder, a automotive physique, a porous medium) you both want a body-fitted mesh that’s costly to generate and onerous to automate, or specialised stencil corrections close to the wall. Each time the geometry modifications, the mesh must be regenerated. It’s an engineering downside layered on high of the physics downside.
LBM sidesteps this virtually fully. The explanation turns into clear when you perceive the place it begins.
A unique place to begin: statistical mechanics
As an alternative of asking “what’s the velocity subject at this level?”, LBM asks “what number of particles are shifting in every route at this level?”
The central object is the single-particle distribution operate f(x, v, t). It tells you the anticipated variety of particles close to place x shifting with velocity v at time t. The macroscopic fields we truly care about are simply moments of f:
[
rho(mathbf{x},t) = mint f, d^3v qquad rho,mathbf{u}(mathbf{x},t) = mint mathbf{v},f, d^3v
]
At world equilibrium (uniform density, no bulk move) f takes the acquainted Maxwell-Boltzmann kind:
[
f^{eq}(mathbf{v}) = frac{rho}{(2pi RT)^{D/2}} exp!left(-frac{|mathbf{v} – mathbf{u}|^2}{2RT}right)
]
Away from equilibrium, f takes another form. The equation governing the way it evolves is the Boltzmann equation:
[
frac{partial f}{partial t} + mathbf{v}cdotnabla_mathbf{x} f = Omega[f]
]
The left-hand aspect is pure transport: particles carry their velocity distribution as they stream via area. The correct-hand aspect, Ω[f], is the collision operator that drives f again towards f^eq by redistributing particles amongst velocity courses. The issue is that Ω[f] is a five-dimensional integral over all potential collision companions and scattering angles, encoding the total microscopic collision physics. For any sensible simulation, it’s fully intractable.
The one-line miracle: BGK
In 1954, Bhatnagar, Gross and Krook proposed a brilliantly blunt simplification. As an alternative of modelling each collision intimately, simply seize the online impact: collisions drive f towards equilibrium at a charge proportional to how removed from equilibrium it at the moment sits.
[
Omega[f] approx -frac{1}{tau}bigl(f – f^{eq}bigr)
]
That is the BGK approximation, and it reduces an intractable five-dimensional integral to a scalar multiplication. The parameter τ is the rest time, that means the common time the system takes to return to native equilibrium after a perturbation. Massive τ means gradual rest; small τ means quick.
What makes this greater than a handy hack is what τ means bodily. By way of the Chapman-Enskog enlargement (a scientific perturbative evaluation) you possibly can show that within the restrict of small Knudsen quantity (many collisions per imply free path), the BGK equation recovers the Navier-Stokes equations precisely, with kinematic viscosity:
[
nu = c_s^2!left(tau – frac{1}{2}right)
]
Observe that this relation operates in dimensionless lattice models; changing to bodily models requires a separate size and timescale calibration. However inside the simulation, τ is your entire bodily knob: flip it up for thick, viscous move; push it towards 0.5 for high-Reynolds-number move (although values beneath ~0.55 danger numerical instability).
Crucially, this is the reason LBM is greater than Navier-Stokes in disguise. The connection to NS is derived, not assumed: you’re fixing a kinetic equation that provably converges to NS habits within the applicable restrict. This can be a quite common theme in all of statistical mechanics.
From steady to discrete: the D2Q9 lattice
The BGK equation nonetheless lives in steady velocity area. To simulate it we want a finite set of velocities. The strategy is a Taylor enlargement of f^eq in powers of the majority velocity u, legitimate when u is small in comparison with the lattice velocity of sound (the low-Mach assumption):
[
f_i^{eq} = rho, w_i!left[1 + 3(mathbf{c}_icdotmathbf{u}) + frac{9}{2}(mathbf{c}_icdotmathbf{u})^2 – frac{3}{2}|mathbf{u}|^2right]
]
The continual velocity c is changed by a finite set of discrete velocities c_i, every carrying a quadrature weight w_i chosen to appropriately recuperate the rate moments. For 2D simulations the usual selection is D2Q9: 2 dimensions, 9 velocity instructions: relaxation, 4 axis-aligned instructions, and 4 diagonals.
The weights are mounted by isotropy: 4/9 for the remainder particle, 1/9 for axis instructions, 1/36 for diagonals. The lattice velocity of sound is cs = 1/√3. Every little thing else within the algorithm follows from these 9 numbers.
The algorithm: two steps, endlessly repeated
With the discrete equilibrium in hand, the total LBM algorithm reduces to 4 operations per timestep, run in a loop till convergence:
- Macroscopic restoration: compute ρ and u at every node from the present populations f_i
- Collision: chill out every f_i towards f_i^eq at charge 1/τ (purely native, no neighbor communication)
- Streaming: propagate every inhabitants to the neighboring node in its route of journey
- Boundary situations: reconstruct lacking populations at area edges and stable partitions
The collision step is the place all of the physics occurs. The streaming step is only kinematic bookkeeping. Their alternation is what produces fluid habits.
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 (stable[x][y]) proceed;
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 seems to matter enormously for efficiency. The naïve push sample (the place every thread at (x,y) writes outgoing populations into neighbouring nodes) causes false sharing: threads write to adjoining reminiscence areas, hammering the identical cache strains and forcing fixed invalidation. The pull (collect) sample beneath fixes this, as a result of every thread reads from its upstream neighbours and solely writes to its personal 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 (stable[x][y]) proceed;
for (int i = 0; i < 9; i++) stable[xsrc][ysrc])
? f[x][y][OPP[i]] // bounce-back from wall
: f[xsrc][ysrc][i]; // stream from upstream neighbour
}
}
}As a result of collision is fully native and streaming solely requires nearest-neighbour reads, the entire algorithm is embarrassingly parallel: each node updates independently, with no world solves, no implicit programs, and no matrix inversions. This maps cleanly onto multi-core CPUs, GPUs, and distributed clusters; LBM implementations on CUDA routinely obtain near-linear scaling with thread rely as much as the reminiscence bandwidth restrict.
Boundary situations: complicated geometry without cost
After every streaming step, nodes adjoining to stable partitions are lacking some incoming populations, as they might have arrived from contained in the stable, which has no fluid to stream. The way you fill these lacking populations defines the boundary situation.
Bounce-back handles no-slip partitions. Any inhabitants that streams right into a stable node is mirrored again in the wrong way: f_ī = f*_i. The bodily image is a particle reversing its velocity on the wall, and the no-slip situation (zero velocity on the boundary) emerges routinely, with no specific velocity constraint wanted. The bodily wall sits midway between the final fluid node and the primary stable node, giving the scheme second-order spatial accuracy.
The elegant consequence is that any impediment form is only a boolean flag. So as to add a cylinder, mark the related nodes as stable. The bounce-back rule handles the boundary physics at each flagged node routinely, with no mesh, no stencil corrections, and no code modifications between geometries.
On the inlet we use the Zou-He scheme: a hard and fast velocity uw is imposed, and the three lacking rightward-pointing populations (f_1, f_5, f_8) are reconstructed persistently from the remaining identified populations and the imposed velocity, recovering the right inlet density with out imposing it independently.
On the outlet, a zero-gradient situation copies all populations from the penultimate column, permitting fluid to depart the area freely with out reflecting spurious waves again upstream.
Outcomes: two regimes
The simulation area is a 400×400 lattice with a round impediment of radius 40 lattice models, centred at one quarter of the area size and offset barely vertically to interrupt symmetry. The inlet velocity is uw = 0.1 lattice models, properly beneath cs ≈ 0.577 (low-Mach assumption happy all through). Two values of τ reveal qualitatively totally different move regimes.
At τ = 0.75 (Re ≈ 96) the wake is regular and symmetric. The distribution operate has relaxed near native equilibrium in all places, as collisions are vigorous sufficient to damp any perturbation earlier than it may well develop. The move is viscous and laminar.
At τ = 0.55 (Re ≈ 480) the image modifications fully.
Vortices shed alternately from the higher and decrease surfaces, forming the Kármán avenue. From the statistical mechanics perspective that is the important thing perception: the vortex shedding is pushed by the departure of f from its Maxwell-Boltzmann equilibrium form. The smaller τ means collisions chill out f towards equilibrium extra slowly, permitting the non-equilibrium element f − f^eq to develop massive sufficient to feed again into the macroscopic momentum subject and maintain the oscillation. A superbly equilibrium fluid would present no such construction.
The boolean-flag strategy means altering geometry prices nothing. Right here is identical solver, unmodified, utilized to a few sq. obstacles in a row:
The wakes couple, the vortex streets work together, and wealthy multi-body move construction emerges. The one change to the code was marking totally different nodes as stable.
Scaling to MareNostrum 5
With #pragma omp parallel for on each the collision and streaming loops, the solver scales to a number of cores with minimal effort. The essential selection is the streaming sample: pull streaming ensures every thread solely owns its output reminiscence, eliminating the false sharing that causes the naïve push implementation to plateau early.
We benchmarked on a full node of MareNostrum 5 on the Barcelona Supercomputing Middle: 2× Intel Xeon Platinum 8480+ (Sapphire Rapids), 112 cores throughout two NUMA domains.
Three regimes emerge cleanly. From 1 to 32 threads, effectivity stays above 75%, reaching a 24.6× speedup, which is genuinely sturdy scaling for a memory-heavy stencil code. At 64 threads, effectivity drops to 42% because the reminiscence controller saturates (the ~30 MB of distribution operate knowledge is being learn and written sooner than the bus can maintain). At 112 threads, the speedup truly falls to 23.5×, as a result of spanning each bodily sockets means half of all reminiscence accesses cross the inter-socket interconnect, including latency on high of already saturated bandwidth. The sensible candy spot for a 400×400 grid is 32 cores on a single socket.
One quantity price pausing on: the single-thread baseline after the pull-streaming rewrite is 31.7 seconds, down from 113 seconds with the unique push implementation. That 3.6× enchancment got here from a code change alone, on a single core, earlier than touching any parallelism. The reminiscence entry sample issues lengthy earlier than the {hardware} does.
When to make use of LBM (and when to not)
LBM is the precise device when geometry is complicated or modifications over time, when it’s essential to run on massively parallel {hardware}, or whenever you’re modelling mesoscale phenomena like velocity slip close to partitions that continuum solvers want empirical corrections for. It’s additionally a superb basis for multiphase simulations, the topic of the subsequent article on this sequence.
It’s not the precise device for high-Mach flows (the low-order Taylor enlargement breaks down above Mach ~0.3), very excessive Reynolds numbers (the τ > 0.5 stability requirement limits how low viscosity can go with out extra superior collision operators), or conditions the place exact bodily unit calibration is important from the outset.
What’s subsequent: multiphase flows with Shan-Chen
The solver described right here handles a single-component, single-phase fluid. Lots of the most bodily attention-grabbing issues contain two immiscible fluids (oil and water, rising bubbles, droplet dynamics in airflow).
The Shan-Chen mannequin extends LBM to deal with these circumstances by including short-range inter-particle forces that produce an efficient equation of state with a section transition. Interfacial stress and section separation emerge from the collision statistics quite than being imposed as boundary situations. Within the subsequent article, I’ll derive the Shan-Chen mannequin, implement it as a modest extension of this solver, and present what spontaneous section separation seems like whenever you let two fluids compete for a similar lattice.
Closing thought
What I discover most satisfying about LBM is what it reveals in regards to the relationship between scales. The Navier-Stokes equations really feel basic: they’re the equations of fluid dynamics, written on the scale we will see and contact. However they’re not basic in any respect. They’re a statistical consequence of many particles colliding, recoverable from a kinetic description within the applicable restrict.
LBM makes that relationship concrete and computational. Viscosity isn’t a parameter you set; it’s a collision timescale. The Kármán vortex avenue isn’t a fluid instability you program; it’s what emerges when collisions can’t sustain with inertia. Each macroscopic phenomenon within the simulation has a direct microscopic interpretation within the statistical mechanics beneath it.
That bridge between scales, made tangible in ~200 strains of C++, is the entire level.
References:
For many who wish to go deeper, the next are the canonical references in roughly rising order of technical problem.
Mohamad, A.A. (2019). Lattice Boltzmann Technique: Fundamentals and Engineering Purposes with Pc Codes. Springer. Essentially the most accessible entry level; labored examples and code all through.
Krüger, T. et al. (2017). The Lattice Boltzmann Technique: Rules and Observe. Springer. The excellent reference; derivations, boundary situations, multiphase fashions, and GPU implementation all lined.
Bhatnagar, P.L., Gross, E.P. & Krook, M. (1954). A Mannequin for Collision Processes in Gases. Bodily Overview, 94(3), 511–525. The unique BGK paper.
Chapman, S. & Cowling, T.G. (1970). The Mathematical Idea of Non-Uniform Gases. Cambridge College Press. For the Chapman-Enskog enlargement in full rigour.
The complete C++ supply, OpenMP benchmark scripts, and Python/NumPy reference implementation for this text can be found on GitHub.
