Fluid Ink Simulator
Paint into a real incompressible fluid. Drag to inject coloured ink and velocity and watch it swirl, curl and diffuse — a live Navier-Stokes solver using Jos Stam's Stable Fluids method with vorticity confinement.
Grid resolution
Ink colour
The incompressible Navier-Stokes equations, term by term
Everything on the canvas above is driven by two equations. The first is momentum conservation for a fluid of constant density; the second is the incompressibility constraint. Written for a velocity field u(x, t):
∂u/∂t = −(u · ∇)u − (1/ρ)∇p + ν∇²u + f
∇ · u = 0
- −(u · ∇)u — advection: the fluid carries itself along. Every parcel of ink and every bit of momentum is transported by the velocity field it sits in. This is the only nonlinear term, and it is the reason fluids are hard.
- −(1/ρ)∇p — pressure gradient: fluid accelerates from high pressure toward low pressure. In an incompressible fluid the pressure is not a material property you look up; it is whatever field is needed at each instant to stop the flow from compressing.
- ν∇²u — viscous diffusion: internal friction. Fast-moving fluid drags its slower neighbours along and is slowed by them, smoothing velocity differences. The coefficient ν is the kinematic viscosity, in m²/s — the "Viscosity" slider above.
- f — external force: anything pushing from outside — gravity, a fan, a stirring spoon. Here it is your pointer: dragging adds a force proportional to how fast you moved, plus a blob of dye at the cursor.
The dye you see obeys a companion equation, ∂d/∂t = −(u · ∇)d + κ∇²d − αd: it is advected by the same velocity field, spreads by molecular diffusion κ ("Dye diffusion"), and fades at rate α ("Dye fade"). The dye is a passive tracer — it is pushed by the fluid but never pushes back, exactly like a drop of food colouring in water.
What "incompressible" means, and why there is a pressure solve
∇ · u = 0 — the divergence of the velocity field is zero everywhere. Divergence at a point measures net outflow from an infinitesimal box around it. Zero divergence means whatever flows into a cell flows out of it: no fluid is created, destroyed, squeezed or expanded. Water at everyday pressures is incompressible to within a fraction of a percent, and so is air below roughly Mach 0.3, which is why aerodynamicists use the incompressible equations for cars and cyclists but not for jets.
The trouble is that if you advance the velocity field using only advection, viscosity and external forces, the result is not divergence-free. Sources and sinks appear, ink piles up in corners, and the picture stops looking like a fluid. So each step ends with a correction called projection.
It rests on the Helmholtz-Hodge decomposition: any vector field w splits uniquely into a divergence-free part and the gradient of a scalar, w = u + ∇p. Take the divergence of both sides and, since ∇ · u = 0, you get the pressure Poisson equation:
∇²p = ∇ · w then u = w − ∇p
Discretised on the grid, that is a large sparse linear system — one unknown per cell. This page solves it with Gauss-Seidel relaxation, about 20 sweeps per projection, which is the same iterative scheme Stam used and is cheap enough to run twice per frame. Note what this implies physically: pressure propagates instantaneously across the whole domain. Poke the fluid on the left and the constraint is felt on the right in the same timestep. That is genuinely how the incompressible limit behaves — it is the limit of infinite speed of sound.
Semi-Lagrangian advection: unconditionally stable, unavoidably blurry
Before 1999, real-time fluid animation kept exploding. Explicit finite-difference advection is bound by the CFL condition: the timestep must be small enough that fluid moves less than one cell per step,Δt < Δx / |u|. Violate it and the simulation diverges to infinity within a handful of frames.
Jos Stam's Stable Fluids (SIGGRAPH 1999) sidestepped it entirely. Instead of asking "where does the stuff in this cell go?", ask "where did the stuff now in this cell come from?" Trace the cell centre backwards along the velocity field for one timestep, land at some arbitrary point between grid nodes, and bilinearly interpolate the previous field there. That is the new value.
Why it can never blow up: a bilinear interpolation of four values is a weighted average with non-negative weights that sum to 1. The new value is therefore always bounded by the maximum of the old field. Nothing can grow. The scheme is unconditionally stable for any timestep and any velocity — which is exactly what a 60 fps interactive toy needs, because you cannot control how hard the user drags.
The price is numerical dissipation. Every step replaces sharp structure with an average of neighbours, so the method behaves as if the fluid had extra artificial viscosity it was never given. Fine vortices smear out and the flow drifts toward a bland, syrupy laminar look — the well-known first-order error of the scheme. Set viscosity to exactly zero in the controls and the ink still goes soft: that residual smoothing is the advection step, not the physics.
Vorticity confinement: putting the eddies back
Vorticity ω = ∇ × u measures local spin. In 2D it is a single scalar per cell: ω = ∂v/∂x − ∂u/∂y. Numerical dissipation destroys vorticity fastest at the smallest scales, which is precisely where the visually interesting curling lives.
Vorticity confinement, introduced to graphics by Fedkiw, Stam and Jensen (2001) and borrowed from Steinhoff's work on helicopter rotor wakes, adds an artificial force that pushes energy back into those eddies. Compute the normalised gradient of vorticity magnitude, which points toward the centre of each existing vortex, then apply a force perpendicular to it:
N = ∇|ω| / |∇|ω||
f_conf = ε · h · (N × ω)
The effect is to spin up whatever swirl the solver still has rather than to invent new ones — it amplifies surviving vortices instead of creating them from nothing. The strength ε is a purely artistic dial with no physical justification; it is a correction for a numerical error, not a term in Navier-Stokes. Drag the vorticity slider to 0 and watch the flow go visibly flatter and creamier within a few seconds; push it to 40 and the ink breaks into tight, persistent curls.
Reynolds number: the dial between laminar and turbulent
The single most important number in fluid dynamics is dimensionless: Re = ρUL/μ = UL/ν, where U is a characteristic speed, L a characteristic length and ν the kinematic viscosity. It is the ratio of inertial forces to viscous forces. Low Re means viscosity wins and the flow is smooth and layered (laminar); high Re means inertia wins, small disturbances grow instead of damping out, and the flow becomes turbulent. Pipe flow transitions around Re ≈ 2,300.
| Situation | Approx. Re | Regime |
|---|---|---|
| Bacterium swimming | 10⁻⁵ – 10⁻⁴ | Creeping (Stokes) flow |
| Blood in a capillary | ~0.001 | Laminar |
| Honey off a spoon | ~0.1 | Laminar |
| Smoke rising from a candle | ~10² – 10³ | Transitional — the visible break-up |
| Water in a garden hose | ~10⁴ | Turbulent |
| Airliner wing at cruise | ~10⁷ | Fully turbulent |
| Blue whale cruising | ~3 × 10⁸ | Fully turbulent |
You can walk that axis with the controls above. Raise viscosity and drop force strength and the ink settles into thick, honey-like sheets that refuse to break up. Drop viscosity toward zero, raise force and turn vorticity confinement up and the flow shreds into fine filaments and curls. The simulation cannot reach genuinely turbulent Reynolds numbers — resolving turbulence needs roughly Re9/4 grid cells in 3D, so a real Re = 10⁴ flow wants millions of cells — but on a 128×128 grid it captures the visual character of the transition well enough that your eye reads it as fluid.
A note on resolution: the cost of each frame scales roughly with the number of cells times the iteration count, so 192² is about 2.25× the work of 128². Higher resolution also means less numerical diffusion per unit distance travelled, so the ink stays sharper — the trade is straightforward, and the fps readout under the canvas tells you where your machine lands.