Fluid Dynamics
The incompressible Navier-Stokes equations govern everything from coffee stirring to Jupiter's Great Red Spot. This simulation solves them in real time using Jos Stam's unconditionally stable algorithm: semi-Lagrangian advection traces fluid parcels backwards in time, and a Helmholtz decomposition projects the velocity field onto the divergence-free manifold at each step. Drag your cursor through the field to inject dye and momentum. The Kelvin-Helmholtz preset shows what happens when two fluid layers slide past each other — the interface is unstable and rolls up into a row of vortices within seconds.
The equations
The incompressible Navier-Stokes equations for a viscous fluid at density ρ = 1:
∂u/∂t = −(u·∇)u − ∇p + ν∇²u + f
∇·u = 0
The second equation is the incompressibility constraint: the fluid has no sources or sinks — whatever flows into a region must flow out. This makes the pressure p a non-local quantity determined by the entire flow field simultaneously, which is the hard part.
Stam's algorithm
Each timestep performs three operations in sequence. First, self-advection: each grid point asks where a fluid parcel at that location came from Δt seconds ago (tracing the velocity field backwards), then interpolates to get the new value. This semi-Lagrangian scheme is unconditionally stable regardless of timestep size. Second, viscous diffusion: an implicit solve of the heat equation on the velocity field, also unconditionally stable. Third, projection: the velocity is made divergence-free by subtracting its gradient component. Solving ∇²p = ∇·u via Gauss-Seidel relaxation (16 iterations) and then correcting u ← u − ∇p enforces incompressibility at each step.
Kelvin-Helmholtz instability
When two fluid layers move past each other at different velocities, the interface is linearly unstable to perturbations of any wavelength. A small ripple in the boundary grows exponentially — the growth rate for wavelength λ is proportional to the shear |Δu| / λ — then rolls up into a row of vortex spirals. The preset seeds the interface with random noise to trigger the instability; without it, the shear layer would remain flat indefinitely (in a perfect solver with no round-off error).
Vorticity and enstrophy
Vorticity ω = ∇ × u measures the local rotation rate of fluid elements. In 2D it reduces to a scalar: ω = ∂v/∂x − ∂u/∂y. The vorticity view colours positive (anticlockwise) rotation red and negative (clockwise) rotation blue — cleaner than the velocity field for identifying coherent structures like vortex cores. Enstrophy Ω = ½ ∫ω² dA is the L² norm of vorticity; in inviscid 2D flow it is conserved, so the readout tracks how dissipation is removing it from the simulation.
Grid and performance
The simulation runs on a 128 × 128 grid with no-slip boundary conditions (velocity
vanishes at walls). The Gauss-Seidel pressure solver runs 16 iterations per step.
At 60 fps on a modern machine this is comfortably real-time in pure JavaScript using
Float32Array typed arrays. The RGB dye is advected by the same velocity
field as the fluid, making it a passive scalar tracer — it reveals the flow geometry
without altering the dynamics.