Fluid That Warms Up While Standing Still — Natural Convection and Rayleigh 1708

The one-third exponent on this same Ra drops the heat flux to about 149. Reaching 1e9 takes a height near 0.73 m, not 0.30 m.

Heat the floor and the fluid still sits there for a while. Only heat moves upward. Then a field of striped vortices switches on at once. The number on that door is the Rayleigh number 1708.

Drop the Ra slider under 1708 and the particles stop while the temperature splits top from bottom. Raise it past 1708 and counter-rotating rolls switch on.

Ra = 1.0e+4

critical Ra_c = 1708 · state: steady convection rolls (Benard cells)

Just above the threshold the circulation is weak. Push the slider further and the boundary layers thin while the cells strengthen.

Where buoyancy and viscosity collide#

Natural convection is flow from a density difference, with no pump. Heat the floor and the fluid above it expands and lightens. A light parcel wants to rise. That is buoyancy.

Buoyancy alone would start the motion at once. Viscosity holds the rising parcel by friction. Thermal diffusion wipes the parcel's temperature, and with it the gap buoyancy needs.

If those two brakes beat buoyancy, the fluid stays still. If buoyancy wins, convection turns on.

The Boussinesq approximation — density change in a single term#

Letting density vary in every term makes the equations heavy. The Boussinesq approximation keeps density as a function of temperature only inside the buoyancy term.

The volumetric expansion coefficient is the fraction by which density drops when temperature rises by one degree.

\beta \equiv -\frac{1}{\rho}\left(\frac{\partial \rho}{\partial T}\right)_P

β is that coefficient, ρ density, T temperature. For an ideal gas, β is 1/T.

The density gap then becomes a temperature gap.

\frac{\rho_\infty - \rho}{\rho} \approx \beta\,(T - T_\infty)

ρ∞ and T∞ are the density and temperature of the still fluid far away. That line plants buoyancy in the momentum equation of a vertical-plate boundary layer.

u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = g\beta\,(T - T_\infty) + \nu\frac{\partial^2 u}{\partial y^2}

u and v are velocity components, g gravity, ν kinematic viscosity. The first term on the right is buoyancy, the second is viscous drag. The energy equation carries temperature by diffusion.

u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y} = \alpha\frac{\partial^2 T}{\partial y^2}

α is the thermal diffusivity.

Grashof and Rayleigh — a division of labor#

Nondimensionalize and the buoyancy term collapses into one number. The Grashof number.

\mathrm{Gr} = \frac{g\beta\,(T_s - T_\infty)\,L^3}{\nu^2} \sim \frac{\text{buoyancy}}{\text{viscous force}}

Ts is the wall temperature, L the characteristic length. Grashof takes the job that the Reynolds number placing a flow by length alone does in forced convection.

Viscosity is not the only opponent. Thermal diffusion also erases the temperature gap. That second brake enters as the Prandtl number, momentum diffusion over heat diffusion, Pr = ν/α. The product is the Rayleigh number.

\mathrm{Ra} = \mathrm{Gr}\cdot\mathrm{Pr} = \frac{g\beta\,(T_s - T_\infty)\,L^3}{\nu\,\alpha}

Both ν and α sit in the denominator. Whether convection turns on is set by Ra, not by Gr alone.

The threshold called 1708#

A linear stability analysis of a layer heated from below returns a critical value. With solid walls above and below, Ra_c is 1708. Below it, disturbances are eaten by viscosity and thermal diffusion. The fluid stays still and heat moves by conduction only.

Past 1708, counter-rotating rolls line up. Those are Bénard cells. The cell width is about twice the layer depth. The Ra slider above is what crosses that door.

Python — tracking the Nusselt number of a vertical plate#

How well a wall cools is the Nusselt number, Nu = hL/k, convection against conduction. A vertical isothermal plate has a Churchill–Chu fit for the whole range. The case is a 60 degree wall in 20 degree air, 0.30 m tall.

import numpy as np
 
def rayleigh_number(Ts, Tinf, L, beta, nu, alpha, g=9.81):
    # ratio of buoyancy to (viscosity x thermal diffusion)
    return g * beta * (Ts - Tinf) * L**3 / (nu * alpha)
 
def churchill_chu_vertical(Ra, Pr):
    # vertical isothermal plate -- valid over the whole range
    denom = (1.0 + (0.492 / Pr)**(9 / 16))**(8 / 27)
    return (0.825 + 0.387 * Ra**(1 / 6) / denom)**2
 
# air properties at the film temperature
nu, alpha, Pr, k = 1.85e-5, 2.60e-5, 0.71, 0.027
Ts, Tinf, L = 60.0, 20.0, 0.30          # 60C wall, 0.3 m tall
beta = 1.0 / (0.5 * (Ts + Tinf) + 273.15)  # ideal gas, 1/Tf
 
Ra = rayleigh_number(Ts, Tinf, L, beta, nu, alpha)
Nu = churchill_chu_vertical(Ra, Pr)
h = Nu * k / L
q = h * (Ts - Tinf)                     # heat loss per unit area
 
print(f"Ra = {Ra:.2e}")
print(f"Nu = {Nu:.1f},  h = {h:.2f} W/m^2K")
print(f"q  = {q:.1f} W/m^2")
 
# Ra = 7.03e+07
# Nu = 54.9,  h = 4.94 W/m^2K
# q  = 197.6 W/m^2

Pr = 0.69Ra marker = 1e7.0

cyan: Churchill–Chu · green: laminar 0.59 Ra^1/4 · pink: turbulent 0.10 Ra^1/3

Lower the Pr slider toward liquid metals and the curve drops. Raise it toward oils and the curve climbs. Push the Ra marker past 1e9 and the slope switches from 1/4 to 1/3.

The printed numbers are Ra 7.03e7, Nu 54.9, h 4.94, and a heat flux of 197.6. The film temperature is 40 plus 273.15, so 313.15 K, and the expansion coefficient is the reciprocal, about 0.003193. Viscosity times thermal diffusivity is 1.85e-5 times 2.60e-5, or 4.81e-10. The cube of the 0.30 m height is 0.027. Those pieces, with the code's gravity 9.81 and a 40 degree gap, meet the printed 7.03e7. Multiply Nu 54.9 by 0.027 and divide by 0.30 and the result is 4.941, the printed h of 4.94. Then 4.94 times 40 is 197.6, so the heat-flux line closes too. Divide this Ra by Pr 0.71 and Gr is about 9.90e7. That number does not belong next to 1708.

On the same Ra the laminar alias, 0.59 times Ra to the 1/4, is about 54.0. The printed 54.9 sits about 0.9 above it. The turbulent alias, 0.10 times Ra to the 1/3, is about 41.3. Rebuilt into h that Nu gives about 3.72, and a heat flux near 149. That is about three quarters of 197.6.

Reaching Ra 1e9 takes about 14.2 times 7.03e7. Height enters as a cube, so 0.30 m has to become about 0.73 m before that door. Short of it, this plate stays on the 1/4 branch. The length that makes Ra 1708 in the same air and the same 40 degree gap is about 8.7 mm. That length is a layer thickness. Read 1708 as the door for the 0.30 m plate and the geometry is the wrong one.

The printed 197.6 is a flux per square meter. Give the plate a width equal to its height, 0.30 m, and the area is 0.09 m^2. The whole plate then loses 197.6 times 0.09, about 17.8 W.

Follow only the laminar alias and Nu scales as Ra to the 1/4 while Ra scales as height cubed, so h moves as height to the -1/4. Doubling the height multiplies h by 2 to the -1/4, about 0.841. Then 4.94 becomes about 4.15, and the flux falls to about 166. Under the 1/3 exponent that drop disappears and h no longer depends on height. The Ra of this print is not on that exponent yet.

A bundle of Nusselt correlations#

Each shape has its own fit.

Vertical plate: $\mathrm{Nu} = 0.59\,\mathrm{Ra}^{1/4}$ (laminar, $10^4$ – $10^9$ ), $\mathrm{Nu} = 0.10\,\mathrm{Ra}^{1/3}$ (turbulent, $10^9$ – $10^{13}$ ). Whole range: Churchill–Chu.
Horizontal heated plate, facing up: $\mathrm{Nu} = 0.54\,\mathrm{Ra}^{1/4}$ ( $10^4$ – $10^7$ ), $0.15\,\mathrm{Ra}^{1/3}$ ( $10^7$ – $10^{11}$ ).
Horizontal heated plate, facing down: $\mathrm{Nu} = 0.27\,\mathrm{Ra}^{1/4}$ .
Enclosure heated from below: past $\mathrm{Ra}_c = 1708$ , Bénard cells.

For a horizontal plate the length is area over perimeter, $L = A_s/P$ . At the turbulent exponent 1/3, $\mathrm{Nu}$ scales as $\mathrm{Ra}^{1/3}$ and $\mathrm{Ra}$ scales as $L^3$ , so $L$ drops out of $h = \mathrm{Nu}\,k/L$ . The heat-transfer coefficient of turbulent natural convection does not hang on size. The 0.30 m plate above is not yet tall enough for that cancellation.

The turbulent exponent on 7.03e7 cuts the heat flux#

Put 1/3 on this print and Nu 54.9 falls to about 41.3, while 197.6 falls to about 149. That exponent erases height, so the 0.30 m disappears with it. On the laminar branch a change in height still moves h. Pick the exponent first and that movement dies for the wrong reason.

1e9 arrives only near 0.73 m. Once the cells are on, how much of this heat flux survives if only three modes remain is the question carried by the Lorenz model that folds the rolls into three equations.