How the flow is calculated
No OpenFOAM inside, sorry. It's just a fun toy of an OpenFOAM user. Everything runs in your browser in about 300 lines of JavaScript. The core is Jos Stam’s Stable Fluids (1999), the classic real-time method, with one upgrade to the advection. It was created by AI to have some fun.
Equations and mesh
The solver integrates the 2D incompressible Navier–Stokes equations with a Boussinesq buoyancy term:
∂u/∂t + (u·∇)u = −∇p + ν∇²u + βgT ĵ + f, ∇·u = 0
The mesh is a uniform Cartesian grid of 176 × 88 cells with all variables collocated at cell centres. Walls, the step, the cylinder and the tree are just a mask of blocked cells, so every curved surface is a staircase. Everything is in grid units, Δx = Δt = 1, and 100 time steps are shown as one second.
One time step, one term at a time
Instead of solving the equation as a whole, each step is split into sub-steps, and each sub-step applies one term. This is exactly why the terms in the game can be switched off: a switch simply skips its sub-step.
- Body forces, explicitly: buoyancy βgT in hotRoom and your wind f in the last case.
- Viscous diffusion, implicit Euler: (I − νΔt∇²)u = u*, a few Jacobi sweeps. No-slip comes for free because blocked cells hold zero velocity.
- Convection, semi-Lagrangian: each cell centre is traced backwards along the local velocity and the old field is interpolated there bilinearly. Plain semi-Lagrangian is very diffusive, so it is wrapped in BFECC (back-and-forth error compensation), with a limiter against overshoots. This is what keeps the Kármán street alive.
- Pressure projection: solve ∇²p = ∇·u* with 20 Jacobi sweeps and subtract ∇p. Zero gradient at walls, p = 0 at the outlet. In the closed rooms this is the only job of the pressure: keeping ∇·u = 0.
- Temperature is a passive scalar, advected semi-Lagrangian, with the heater cells held at fixed T.
Compared with pimpleFoam this is a non-iterative fractional-step projection (Chorin’s method): no PISO correctors, no outer loops, first-order splitting error, Jacobi instead of GAMG and no Rhie–Chow, so the collocated checkerboard is just smoothed away. The “residual” in the log is the mean divergence before projection, which is at least an honest number.
What is faked
Δp is not a pressure boundary condition. The bulk inlet velocity comes from a lumped law, Δp = aνU + bU², a Hagen–Poiseuille-like part plus a quadratic loss, relaxed in time so the flow spins up. It is imposed as a plug profile with smoothed edges.
The crash is a game rule. Semi-Lagrangian advection is unconditionally stable, so this solver would happily run at Courant 10. The fatal error at Co > 2 is there for the drama.
The Reynolds number is optimistic. At low ν the numerical diffusion of the advection is larger than the physical viscosity, so the effective Re is lower than the one in the corner. The shedding, recirculation and plumes are nevertheless produced by the flow itself, not scripted.
Units are decorative. The ν shown is the grid value rescaled so that the lowest setting looks like water; the fluid names, Pa, K and planets are for fun.
Particles are one-way-coupled tracers: xn+1 = xn + u(xn)Δt − wsĵ, with a settling velocity ws for the heavy ones and a small random walk for the dye. A leaf stays on the tree until the local speed exceeds its grip for a few steps; then it falls with a little flutter.
The term bars show the RMS change of velocity produced by each sub-step, normalised by the largest one. When the flow becomes steady, ∂u/∂t drops and the other terms balance each other, which is the equation doing exactly what it says.
Further reading
A. J. Chorin, Numerical solution of the Navier–Stokes equations, Math. Comp. 1968.
J. Stam, Stable fluids, SIGGRAPH 1999.
B. Kim, Y. Liu, I. Llamas, J. Rossignac, FlowFixer: using BFECC for fluid simulation, 2005.