The Dance of Fluids

Building a Navier-Stokes solver from scratch

2026-03-18 simulation

Watch smoke rise from a candle. Pour cream into coffee. Look at clouds shearing against the wind. You're watching the Navier-Stokes equations in action — the same equations that govern weather systems, ocean currents, the airflow over a wing, and the blood flowing through your veins.

These equations are among the most important in physics, and among the most difficult to solve. Whether smooth solutions always exist in 3D is one of the seven Millennium Prize Problems — a million dollars for a proof or disproof.

But we don't need to solve them analytically. We can simulate them. In 1999, Jos Stam published "Stable Fluids," a beautifully simple algorithm that simulates incompressible flow with unconditional stability. The core idea fits in a few hundred lines of code. I built one from scratch — first in Python for analysis, then in JavaScript for the interactive demo below.

The Equations

The Navier-Stokes equations for an incompressible fluid describe two things: how the velocity field changes over time, and the constraint that fluid cannot be created or destroyed.

∂u/∂t = -(u·∇)u + ν∇²u - ∇p/ρ + f

∇·u = 0

The first equation is Newton's second law applied to a fluid element. Each term has a physical meaning:

-(u·∇)u

Advection. The fluid carries itself. Each point's velocity is transported by the flow — this is why smoke curls and eddies form.

ν∇²u

Diffusion. Viscosity smooths out velocity differences. Honey has high ν; air has low ν. This is why honey flows smoothly and air turbulates.

-∇p/ρ

Pressure. Fluid flows from high pressure to low. This term enforces incompressibility — it's the "push back" that prevents fluid from piling up.

f

External forces. Gravity, your mouse dragging through the fluid, a fan blowing air — anything acting on the fluid from outside.

The second equation, ∇·u = 0, is the incompressibility constraint. It says the divergence of the velocity field is zero — fluid neither compresses nor expands. What flows into a region must flow out. This is the equation that makes fluids hard to simulate: you can't just integrate velocities forward, because the result might violate this constraint.

Stam's Insight: Operator Splitting

The genius of Stable Fluids is to solve each physical effect separately, in sequence. Instead of tackling the full equation at once, we split it into four steps, each handling one term:

Add Forces

Apply external forces to the velocity field. Simple: u += dt * f

Diffuse

Spread velocity via viscosity. Solved implicitly for stability.

Advect

Transport velocity by itself. Semi-Lagrangian: trace backwards, interpolate.

Project

Enforce ∇·u = 0 by subtracting the gradient of pressure.

This isn't physically exact — splitting introduces a small error proportional to dt. But it's unconditionally stable: no matter how large the time step or how fast the flow, the simulation never blows up. That's Stam's key contribution. Previous methods required tiny time steps (the CFL condition) or they'd explode. Stable Fluids trades some accuracy for guaranteed stability.

Try It Yourself

Click and drag on the canvas below to inject fluid and velocity. The simulation runs the complete Stable Fluids algorithm in your browser at 60fps.

Color: Viscosity:

Click and drag to paint fluid. Right-click to add density without velocity.

Semi-Lagrangian Advection

The advection step is where the magic happens. The naive approach — moving fluid forward along velocity vectors — is numerically unstable. Push too much fluid into a cell and values explode.

Stam's trick is to work backwards. For each grid cell, ask: "Where did the fluid that's now here come from?" Trace a virtual particle backwards along the velocity field by one time step, find where it lands, and interpolate the value from the surrounding grid cells.

// For each grid cell (i, j):
// 1. Trace backwards along velocity
x = i - dt * u[i][j]
y = j - dt * v[i][j]

// 2. Clamp to grid bounds
x = clamp(x, 0.5, N + 0.5)
y = clamp(y, 0.5, N + 0.5)

// 3. Bilinear interpolation from source position
d_new[i][j] = bilerp(d_old, x, y)

This is unconditionally stable because we're always interpolating between existing values. The result is bounded by the input — no matter how fast the flow, no cell can get a value larger than what already exists in the field. The trade-off is numerical dissipation: interpolation slightly blurs the field each step, which is why simulated fluids gradually smooth out over time.

The Pressure Projection

After advection, the velocity field generally has non-zero divergence — fluid is "piling up" in some places and leaving voids in others. The projection step fixes this using a beautiful mathematical result: the Helmholtz-Hodge decomposition.

Any vector field w can be uniquely decomposed into a divergence-free part u and the gradient of a scalar field p:

w = u + ∇p where ∇·u = 0

Taking the divergence of both sides and using ∇·u = 0:

∇·w = ∇²p

This is Poisson's equation — we can solve it iteratively using Gauss-Seidel relaxation. Once we have p, we subtract its gradient from w to get the divergence-free velocity field u. This is the pressure projection, and it's what enforces incompressibility.

// 1. Compute divergence of velocity
div[i][j] = -0.5 * h * (u[i+1][j] - u[i-1][j]
                       + v[i][j+1] - v[i][j-1])

// 2. Solve ∇²p = div (Gauss-Seidel, ~20 iterations)
p[i][j] = (div[i][j] + p[i-1][j] + p[i+1][j]
                      + p[i][j-1] + p[i][j+1]) / 4

// 3. Subtract gradient of pressure
u[i][j] -= 0.5 * N * (p[i+1][j] - p[i-1][j])
v[i][j] -= 0.5 * N * (p[i][j+1] - p[i][j-1])

The physical interpretation: when fluid converges (positive divergence), pressure builds up, pushing fluid outward. When fluid diverges (negative divergence), low pressure pulls fluid inward. Pressure is the mechanism by which incompressibility is enforced.

The Gallery

I built the simulator in Python first, using NumPy for vectorized operations. Here are scenes from six different scenarios, each revealing different aspects of fluid dynamics.

Rising Smoke

Smoke plume rising with turbulent eddies

Left: A smoke plume rising from a heat source, with turbulent eddies forming at the boundary between the fast-moving plume and the still air. Right: Four jets colliding at the center, creating a violent turbulent mixing zone.

Vortices and Flow

Left: A counter-rotating vortex pair visualized by vorticity — blue is clockwise rotation, red is counterclockwise. The two vortices orbit each other as they propagate. Right: Kelvin-Helmholtz instability where two fluid layers moving in opposite directions develop characteristic rolling waves.

Ink and Wind

Left: A drop of ink falling through water, forming a mushroom-shaped vortex ring. Right: Wind tunnel flow around a circular obstacle, showing the wake region downstream.

Why It Matters

Stam's algorithm (and its descendants) power the fluid effects in almost every film, game, and visualization you've seen. The key insight — that unconditional stability matters more than accuracy for visual simulation — was controversial when published but proved exactly right for interactive applications.

Real-world computational fluid dynamics (CFD) uses more sophisticated methods: finite volume methods, spectral methods, adaptive mesh refinement, turbulence models. But the core mathematical structure is the same: advect, diffuse, project. The Navier-Stokes equations are the same whether you're simulating a candle flame or designing a jet engine — the difference is resolution, boundary conditions, and how much accuracy you need.

What I find most beautiful about fluid simulation is the gap between the simplicity of the equations and the complexity of the behavior. Four terms. Two constraints. The mathematics fits on a napkin. But from those four terms emerge turbulence, vortex streets, Kelvin-Helmholtz instabilities, the fractal cascade of energy from large eddies to small ones — all the endless, beautiful complexity of flowing water.

The equations are simple. The solutions are not. And that, perhaps, is why smooth solutions in 3D remain a million-dollar question.