Rayleigh-Bénard Convection

The Rayleigh-Bénard problem couples incompressible fluid flow with heat transport: a horizontal cavity heated from below, cooled from above, develops convective rolls once the buoyancy force exceeds the dissipative effects of viscosity and thermal diffusion. The script examples/fluid/rayleigh_benard/rayleigh_benard.py solves the steady-state Boussinesq approximation in a 2:1 aspect-ratio rectangular cavity at a configurable Rayleigh number.

This is the only example in Fluid Mechanics that involves a multi-physics coupling — momentum, continuity, and energy in one monolithic system.

Problem

The Boussinesq equations,

\[ \begin{align}\begin{aligned}\rho\, (\mathbf{u} \cdot \nabla)\mathbf{u} \;=\; -\nabla p + \mu\, \Delta \mathbf{u} - \rho\, g\, \beta\, T\, \hat{\mathbf{e}}_y, \qquad \nabla \cdot \mathbf{u} = 0,\\\rho\, c_p\, (\mathbf{u} \cdot \nabla T) \;=\; \kappa\, \Delta T,\end{aligned}\end{align} \]

on \(\Omega = [0, 2] \times [0, 1]\), with

  • velocity: no-slip on all walls,

  • temperature: \(T = 1\) on the bottom (\(y = 0\)), \(T = 0\) on the top (\(y = 1\)), no-flux on the side walls.

The Rayleigh number

\[\mathrm{Ra} \;=\; \frac{g\, \beta\, \Delta T\, L^3}{\nu\, \alpha}\]

controls the regime: below \(\mathrm{Ra}_c \approx 1708\) heat transfer is purely conductive; above it convective rolls appear; at higher \(\mathrm{Ra}\) the rolls become unsteady and eventually chaotic. The script defaults to \(\mathrm{Ra} = 10^4\), comfortably in the steady-convective regime.

Coupled assembler

The unknowns at each node are \((u, v, p, T)\) — four DOFs per node — and the assembler returns a \(4 \times 4\) block per quadrature point:

Listing 18 examples/fluid/rayleigh_benard/rayleigh_benard.py (essence)
class RayleighBenardAssembler(ElementAssembler):
    def __post_init__(self, rho=1.0, mu=0.01, kappa=0.01,
                      g=9.81, beta=0.1, tau=0.1):
        self.rho, self.mu, self.kappa = rho, mu, kappa
        self.g, self.beta, self.tau   = g, beta, tau

    def forward(self, u, v, gradu, gradv, w_prev, T_prev):
        # row 0: x-momentum    [convection+diffusion, 0,    -dxp,   0   ]
        # row 1: y-momentum    [0,    convection+diffusion, -dyp,   buoyancy ]
        # row 2: continuity    [dxv,    dyv,                 PSPG,    0   ]
        # row 3: energy        [0,      0,                     0,    convection+diffusion]
        ...

The buoyancy block at row 1 / col 3 is what couples temperature back into the momentum equation:

\[K_{1,3} \;=\; -\rho\, g\, \beta\, v\, u,\]

where here \(u\) and \(v\) are the trial- and test-shape function values (note the unfortunate name collision with the velocity \(v\) field — see Forms for the dispatch convention). The temperature equation (row 3) re-uses the convection / diffusion block but with thermal conductivity \(\kappa\) in place of viscosity \(\mu\).

The previous-iterate velocity w_prev and temperature T_prev arrive as point_data kwargs, exactly the same pattern as the single-physics cavity example.

Stabilization

The script uses PSPG only (no SUPG), with a constant \(\tau = 0.05 / n_\text{grid}\). At \(\mathrm{Ra} = 10^4\) the velocities are mild enough that SUPG is not strictly required; for higher Rayleigh numbers, swap in the adaptive \(\tau\) from Cylinder Flow (Vortex Shedding).

Picard iteration

The same fixed-point loop as Lid-Driven Cavity. Each iteration:

  1. Compute \(\mathbf{w}^n, T^n\) from the current solution vector.

  2. Reassemble \(K\) with the new point_data.

  3. Apply the Condenser for the velocity and temperature Dirichlet BCs (different masks per component!) and solve.

  4. Check convergence against the previous iterate.

Typical convergence: ~30 iterations at \(\mathrm{Ra} = 10^4\).

Rayleigh-Bénard temperature and velocity magnitude

Fig. 48 Output of rayleigh_benard.py. Left: temperature field — the warm bottom wall sends a rising plume up the centerline that splits at the cold top into two cool downward plumes along the side walls. Right: velocity magnitude — two counter-rotating convection rolls, with peak speeds along the rising / sinking columns and stagnation at the roll centers.

Output

rayleigh_benard.png shows two panels: temperature (with the characteristic warm rising plumes / cool sinking plumes) and velocity magnitude (with the two-roll structure). For a 2:1 cavity at moderate \(\mathrm{Ra}\) the steady solution features two counter-rotating rolls.

Running it

cd examples/fluid/rayleigh_benard
python rayleigh_benard.py     # writes rayleigh_benard.png

Edit the ra= argument inside the script to sweep the Rayleigh number; values much above \(10^5\) will need a finer mesh and likely a transient driver.

What’s next