Chaos That Begins at the Bottom of a Pot — Rayleigh–Bénard Convection and the Lorenz Equations

In 1900 Henri Bénard heated a thin oil layer from the bottom and watched it split itself into a hexagonal honeycomb. Nobody stirred it. Sixteen years later Lord Rayleigh wrote the threshold as one number. Forty-seven years after that, Lorenz pulled chaos out of the same system.

The threshold is Rayleigh number 1708. Below it, conduction. Above it, rolls. Cut the roll down to three equations and the heating parameter ρ is what remains. Past 24.74 the fixed points no longer hold the trajectory.

Drop the Ra slider under 1708 and the particles nearly freeze. Push it above 1708 and counter-rotating rolls switch on.

Below 1708 the temperature only stacks in horizontal layers. Past the threshold, pushing the slider farther makes the same rolls spin faster.

Why a Still Liquid Starts to Move#

Hot on the bottom, cold on top. The warmed fluid expands and gets lighter. Light fluid below, heavy fluid above. Buoyancy wants the stack flipped.

It does not jump at once. Viscosity holds the motion. Thermal diffusion leaks heat out of a rising blob and kills the buoyancy that lifted it. Both are brakes.

If buoyancy wins, a circulation appears: a convection roll. If the brakes win, the liquid stays put and heat rides conduction only. Buoyancy can also be spent in the path before a roll ever forms. The draft that loses 7 Pa when the smoke path drops by 1 m is that budget.

The Rayleigh Number — Overturning Force over Restraining Force#

Pressed into one dimensionless number, the tug of war is the Rayleigh number.

$g$ is gravitational acceleration, $\beta$ the thermal expansion coefficient (how much volume grows as temperature rises), $\Delta T$ the top-to-bottom temperature difference, $L$ the layer thickness, $\nu$ the kinematic viscosity, and $\alpha$ the thermal diffusivity. The numerator overturns. The denominator holds.

The heavy term is $L^3$. Thicken the layer and Ra grows as the cube. A thin oil film rarely convects. The deep ocean and the mantle move on a tiny temperature difference. A slightly deeper pot walks the threshold much closer.

Critical Rayleigh Number 1708 — the Threshold Where Convection Turns On#

Rayleigh's linear stability analysis changes the number with the boundary. Free top and bottom: $\mathrm{Ra}_c = 27\pi^4/4 \approx 657.5$. The rigid walls common in experiments: $\mathrm{Ra}_c \approx 1708$.

When $\mathrm{Ra} < \mathrm{Ra}_c$, a small wobble is eaten by viscosity and diffusion. The fluid returns to rest. The moment $\mathrm{Ra} > \mathrm{Ra}_c$, a wobble of one wavelength grows on its own. That is the seed of a roll.

Divide 1708 by 657.5 and the ratio is about 2.6. Same fluid, same depth: rigid walls ask for about 2.6 times the temperature difference. Paid in thickness instead, the cube root is $2.6^{1/3} \approx 1.37$. That is the cost of swapping a free-slip lid for a rigid one.

The Shape of a Roll and the Nusselt Number#

Just past threshold the roll picks a wavelength. The critical wavenumber is $k_c \approx 3.117/L$. One wavelength is $2\pi/k$, and dividing by $3.117$ gives about $2.02L$. Two rolls side by side fill that wavelength, so one roll is about $L$ wide. "Square" is that division. The mode that switches on at the lowest Ra has that wavelength.

How much extra heat leaves is the Nusselt number (total heat transfer over pure conduction).

$q_\text{total}$ is the actual heat flux and $q_\text{conduction}$ the flux with no convection. Conduction alone means $\mathrm{Nu}=1$. Once convection is strong, experiments sit near $\mathrm{Nu} \sim \mathrm{Ra}^{1/3}$. Read only that exponent and an eightfold Ra doubles Nu. Right at onset it does not apply. Nu is still near 1 there.

The Convection Lorenz Reduced to Three Equations#

In 1963 Edward Lorenz kept only the coarsest modes of temperature and flow. A Galerkin truncation. The partial differential equations collapse to three ordinary ones.

$x$ is the roll's circulation strength, $y$ the temperature contrast between rising and sinking columns, $z$ the distortion of the vertical temperature profile. $\sigma$ is the Prandtl number (momentum diffusivity over thermal diffusivity), $\beta$ a geometric constant, and $\rho$ is $\mathrm{Ra}/\mathrm{Ra}_c$.

Raising $\rho$ is heating the pot harder. For $\rho<1$ every trajectory goes to the origin. Conduction. For $1<\rho<24.74$ it settles on an off-origin fixed point. A steady roll. Past $\rho>24.74$ it never sits on a fixed point and wanders between two wings. Chaos.

A steady state is where all three derivatives are zero. If $\sigma \neq 0$, then $y=x$. If $x \neq 0$, then $z=\rho-1$. The remaining equation is $x^2 = \beta(\rho-1)$. With the code's $\beta = 8/3$ and $\rho = 15$, $z=14$ and $x^2 = 112/3 \approx 37.3$. The middle printed row, steady roll, means the last 4000 samples of $x$ span less than 0.1. At $\rho=28$, $z=27$ and $x^2=72$, so the two wings sit there. The printed chaos means that tail failed the 0.1 gate.

Python — Integrating the Three Faces of Convection#

Fourth-order Runge–Kutta. $\sigma=10$, $\beta=8/3$, step 0.01, 20000 steps. The tail is the last 4000 samples.

import numpy as np
 
SIGMA, BETA = 10.0, 8.0 / 3.0
 
def lorenz_deriv(state, rho):
    x, y, z = state
    return np.array([
        SIGMA * (y - x),        # change in roll circulation strength
        x * (rho - z) - y,      # change in rising/sinking temperature contrast
        x * y - BETA * z,       # change in vertical temperature distortion
    ])
 
def integrate_rk4(rho, dt=0.01, n=20000, s0=(0.6, 0.4, 0.2)):
    s = np.array(s0, dtype=float)
    traj = np.empty((n, 3))
    for i in range(n):
        k1 = lorenz_deriv(s, rho)
        k2 = lorenz_deriv(s + 0.5 * dt * k1, rho)
        k3 = lorenz_deriv(s + 0.5 * dt * k2, rho)
        k4 = lorenz_deriv(s + dt * k3, rho)
        s = s + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
        traj[i] = s
    return traj
 
def convection_regime(rho):
    tail = integrate_rk4(rho)[-4000:]          # tail after dropping the transient
    spread = tail[:, 0].max() - tail[:, 0].min()
    if abs(tail[-1, 0]) < 1e-3 and rho < 1.0:
        return "conduction"                     # converges to the origin
    if spread < 0.1:
        return "steady roll"                    # converges to a fixed point
    return "chaos"                              # oscillates between two wings
 
for rho in (0.8, 15.0, 28.0):
    print(f"rho={rho:5.1f} -> {convection_regime(rho)}")

Three lines come back.

rho=  0.8 -> conduction
rho= 15.0 -> steady roll
rho= 28.0 -> chaos

Set the ρ slider below 24.74 and the trajectory falls into one yellow fixed point. Cross it and the butterfly wings open. Reset shifts the start a little. Two paths that overlap at first split soon after.

Stopping at 1708 and calling the job finished is a failure#

Double the depth and $L^3$ multiplies Ra by eight. A layer that had just hit 1708 then sits at $\rho = \mathrm{Ra}/\mathrm{Ra}_c = 8$. The Lorenz steady-roll window runs from 1 to 24.74. Eight is inside it. Chaos wants ρ at 24.74, which at the same temperature difference means a depth of $24.74^{1/3} \approx 2.9$ times. Labeling the roll just past onset as chaos misses that gap.

The gap is built by raising $\Delta T$ while holding $\beta$ fixed. Once the temperature difference passes 30 K, does that constant fail first? Where Boussinesq breaks past that temperature shows which number, 1708 or 24.74, loses its meaning sooner.