A thin fluid layer is held hot at the bottom and cold at the top. When the temperature difference is small, heat simply diffuses upward by conduction and nothing moves. Once the Rayleigh number exceeds the critical value Rac ≈ 1708, buoyancy overwhelms viscous damping and the fluid spontaneously organises into rotating convection rolls that carry heat far more efficiently.
Ra = g·β·ΔT·H³ / (ν·κ) — Rayleigh number
∂T/∂t + (u·∇)T = κ ∇²T — temperature advection-diffusion
∂ω/∂t + (u·∇)ω = ν ∇²ω + g·β·∂T/∂x — vorticity with buoyancy
∇²ψ = -ω, u = ∂ψ/∂y, v = -∂ψ/∂x — stream function
The hexagonal and roll patterns you see here also paint the surface of the Sun: each bright "granule" is the top of a convection cell roughly the size of a continent, lasting only a few minutes before sinking back down.
Rayleigh-Bénard convection is the buoyancy-driven motion that develops in a thin horizontal fluid layer heated uniformly from below and cooled from above. When the heating is strong enough, the still fluid becomes unstable and breaks into regular, rotating convection cells called rolls.
The Rayleigh number Ra is a dimensionless ratio comparing the strength of buoyancy (which drives convection) to the damping effects of viscosity and thermal diffusion. Ra = (g·β·ΔT·H³)/(ν·κ). Higher Ra means buoyancy dominates and convection is more vigorous.
For a fluid layer between two rigid, no-slip plates, linear stability theory predicts that convection first sets in at Ra_c ≈ 1708. Below this value, heat is carried only by conduction and the fluid stays still; above it, infinitesimal disturbances grow into convection rolls.
The Boussinesq approximation treats the fluid density as constant everywhere except in the buoyancy (gravity) term, where small temperature-induced density changes are kept. This greatly simplifies the equations while capturing the essential physics of thermal convection.
Among all possible disturbances, the one with a particular horizontal wavelength (about twice the layer depth) grows fastest. This mode wins the competition and organises the flow into evenly spaced counter-rotating rolls, the lowest-energy way to transport the imposed heat.
Conduction moves heat molecule by molecule, slowly. Convection physically carries hot parcels of fluid upward and cold parcels downward, advecting thermal energy in bulk. Above Ra_c the convective heat flux quickly exceeds the conductive flux, measured by the Nusselt number.
The Nusselt number Nu is the ratio of total heat transport to the heat that pure conduction alone would carry. Nu = 1 means no convection; Nu > 1 means convection is enhancing heat transfer. It rises with the Rayleigh number.
It appears in the Sun's surface granulation, in the Earth's mantle and outer core, in the atmosphere as cloud streets and hexagonal cells, in the oceans, and even in a pot of heated soup or oil. It is one of the most studied examples of pattern formation.
The Prandtl number Pr = ν/κ compares momentum diffusion to thermal diffusion. Low Pr fluids (like liquid metals) develop large-scale, inertia-dominated flows, while high Pr fluids (like oils) form smoother, more viscous rolls. It shapes the texture of the convection.
It solves a 2D Boussinesq model on a lattice grid: temperature is advected and diffused, the temperature difference produces a buoyancy force, and a stream-function / vorticity update enforces an incompressible velocity field. Colour shows temperature and arrows show velocity.
This simulation runs a 2D Boussinesq model of Rayleigh-Bénard convection on a 96×48 lattice grid. Each step advects and diffuses a temperature field, adds a buoyancy force proportional to the Rayleigh number (Ra) and the local temperature gradient, then solves the Poisson equation ∇²ψ = -ω (where ψ is the stream function and ω is vorticity) with 20 Gauss-Seidel iterations to recover a divergence-free velocity field from u = ∂ψ/∂y, v = -∂ψ/∂x. The bottom row is fixed hot (T=1), the top row fixed cold (T=0), and the side walls use Thom's no-slip vorticity formula. Raise Ra above the critical value 1708 and stationary conduction gives way to rotating convection rolls.
A cross-section of fluid trapped between a heated floor and a cooled ceiling, coloured from blue (cold) to orange-red (hot). Below the critical Rayleigh number the field stays smoothly stratified — pure conduction; above it, the buoyancy term destabilises the layer and the colour map reorganises into alternating rising and sinking convection rolls, with white arrows showing the local velocity field.
Drag the Rayleigh number slider (500-50,000) to push the system past Rac ≈ 1708 and trigger convection; the Prandtl number slider (0.1-10) changes how viscous vs. thermal the fluid behaves, reshaping roll texture. Sim speed scales the number of physics substeps per frame, and the velocity-arrows toggle hides or shows the flow overlay. Click or drag directly on the canvas to stir the fluid by hand, injecting vorticity and heat at that point. Pause/Play and Reset control playback, and the stat box reports the live state, Ra/Rac ratio and Nusselt number.
The Nusselt number shown in the HUD (Nu = 1 + a normalised convective heat flux) is exactly the quantity physicists measure in real Rayleigh-Bénard experiments with tanks of water or gas to compare against theory: Nu stays at 1 for pure conduction and climbs as convection strengthens, following well-studied scaling laws with Ra.
The stream-function/vorticity method requires solving a Poisson equation with Gauss-Seidel relaxation every single frame (20 iterations here), and that cost grows quickly with grid size. A 96×48 lattice is coarse enough to update in real time in the browser while still being fine enough to resolve the roll pattern, whose natural wavelength is roughly twice the layer depth.
The stirAt() function converts your pointer position into the nearest interior grid cell and adds +0.6 to its vorticity and +0.3 to its temperature (capped at 1). This lets you manually kick-start convection cells below Rac, or disrupt an existing roll pattern to watch it reorganise.
The bottom row is pinned to T = 1 (hot) and the top row to T = 0 (cold) every step, which is what drives the whole system. The side walls are insulated (zero temperature gradient) so no heat escapes sideways. All four walls are rigid and no-slip, enforced on vorticity with Thom's formula ωwall = -2·ψadjacent, which reproduces zero velocity at the boundary without needing a ghost layer.
Instead of tracking pressure, the code uses the vorticity-stream function formulation: vorticity ω is transported and diffused directly, then ∇²ψ = -ω is solved (via 20 Gauss-Seidel iterations per frame) to recover a stream function ψ whose derivatives, u = ∂ψ/∂y and v = -∂ψ/∂x, are automatically divergence-free by construction — incompressibility is guaranteed by the math, not enforced with an extra pressure-correction step.
It checks two things: whether Ra/Rac is above 1 and whether the average kinetic energy of the velocity field has grown past a small threshold (1e-4). Near the threshold the ratio alone can be misleading — the perturbation added at seed time takes a moment to grow into visible rolls — so the badge only switches to "Convecting" once real motion is detected.