[Paper Review] When 30 K Breaks Boussinesq — Demou (2018) NOB Natural Convection and the Pressure Split

A 30 K water gap is 24 times the 1.25 K Gray–Giorgini limit, and the density correction stays near one percent.

Gray and Giorgini left two numbers in 1976. For room-temperature air, Oberbeck–Boussinesq stays inside a 10% error only up to 28.6 K. For water, the limit is 1.25 K. Solar-tower air swings through hundreds of kelvins. PV/T water warms by 50 K. Demou, Frantzis and Grigoriadis (2018) drop the assumption. Properties follow temperature. Only the pressure Poisson stays constant-coefficient.

Move the hot and cold walls and, at the same ΔT, the density curve stays almost flat while the viscosity curve tips hard.

T_hot: 333 K (60 °C)T_cold: 293 K (20 °C)

Dashed line marks the reference value at T₀ = 313 K. Δ% is the relative change between hot and cold. For ΔT ≈ 50 K, density moves about 2%, but viscosity changes by tens of percent — the OB approximation hides that.

Across 50 K the density moves about 2%. Viscosity moves about 30% on that same interval, and that gap builds the asymmetric boundary layer.

Paper card#

Authors: A. D. Demou, C. Frantzis, D. G. E. Grigoriadis
Journal: International Journal of Heat and Mass Transfer 121 (2018) 1148–1162
DOI: 10.1016/j.ijheatmasstransfer.2018.04.144
Contribution: in non-OB natural convection, keep temperature-dependent properties and turn the variable-coefficient pressure Poisson into a constant-coefficient Poisson plus an explicit correction, by the Dodd–Ferrante split. The FFT solver stays. The reported speedup is 2.5–5×.

Where Boussinesq actually breaks#

Standard Boussinesq holds three statements at once. Density varies only in the gravity term. Every other property is constant. Viscous dissipation is dropped. Once ΔT grows, viscosity $\mu(T)$ , conductivity $k(T)$ , and specific heat $C_p(T)$ move together. In a liquid the density shift can be small while the viscosity shift is not.

Keep every property and the momentum equation picks up a variable coefficient.

\frac{\partial u_i}{\partial t} + \frac{\partial(u_i u_j)}{\partial x_j} = -\frac{1}{\hat\rho}\frac{\partial P}{\partial x_i} + \frac{1}{\hat\rho}\frac{\Pr}{\sqrt{\mathrm{Ra}}}\frac{\partial}{\partial x_j}\!\left[\hat\mu\!\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)\right] + \frac{1}{\mathrm{Fr}^2}\delta_{i3}

Hats are properties normalized by the reference-temperature value. The same variability enters the energy equation.

The variable Poisson trap#

A plain fractional-step correction asks the projection step for a divergence-free field. That step is what comes after Chorin removes pressure from the momentum equation.

\frac{\partial}{\partial x_i}\!\left(\frac{1}{\hat\rho^{n+1}}\frac{\partial P^{n+1}}{\partial x_i}\right) = \frac{1}{\Delta t}\frac{\partial u_i^*}{\partial x_i}.

The factor $1/\hat\rho^{n+1}$ sits inside the divergence. When it varies in space and time, the matrix changes every step. An FFT direct solver accepts only a constant-coefficient Laplacian. The fallback is multigrid or an iterative solver. If 60–80% of a Navier–Stokes step is the pressure solve, that fallback is the bill.

Dodd–Ferrante splitting — one line#

The borrowed trick comes from two-fluid flow (Dodd & Ferrante 2014). The pressure gradient is cut in two.

\frac{1}{\hat\rho^{n+1}}\frac{\partial P^{n+1}}{\partial x_i} \longrightarrow \frac{1}{\hat\rho_0}\frac{\partial P^{n+1}}{\partial x_i} + \left(\frac{1}{\hat\rho^{n+1}} - \frac{1}{\hat\rho_0}\right)\frac{\partial \widehat P}{\partial x_i}

The first piece multiplies the new pressure by the constant $1/\hat\rho_0$ . The second multiplies the extrapolated pressure $\widehat P = 2P^n - P^{n-1}$ . Take the divergence and the Poisson is constant-coefficient.

\frac{\partial^2 P^{n+1}}{\partial x_i \partial x_i} = \frac{\hat\rho_0}{\Delta t}\frac{\partial u_i^*}{\partial x_i} + \frac{\partial}{\partial x_i}\!\left[\left(1 - \frac{\hat\rho_0}{\hat\rho^{n+1}}\right)\frac{\partial \widehat P}{\partial x_i}\right]

The Laplacian on the left is the same every step. FFT in two directions leaves a stack of Helmholtz problems in $z$ , solved directly. Only the right-hand correction is explicit.

The size that matters is $(1 - \hat\rho_0/\hat\rho)$ . It sits near 1 in two-fluid flow and is tiny in water NOB convection.

A NumPy walk-through#

Sections 2.2–2.3 are compressed to one 1D slab. No 2D DNS. One pressure correction.

import numpy as np
from numpy.fft import fft, ifft
 
# 1D slab: periodic, one Helmholtz stack
Nz = 64
Lz = 1.0
dz = Lz / Nz
z = (np.arange(Nz) + 0.5) * dz
 
# temperature-dependent density (water polynomial, simplified)
def rho_hat(T_field, T0=313.0):
    dT = T_field - T0
    return 1.0 * (1 - 3.736e-4 * dT - 3.98e-6 * dT**2)
 
# initial temperature profile: linear + small perturbation
T0 = 313.0
dT = 30.0
T_field = T0 - dT/2 + dT * z + 0.1 * np.sin(2*np.pi*z)
rho = rho_hat(T_field)
 
# reference density: the lightest (hottest) value
rho_ref = rho.min()
 
# fake divergence (in a full solver this comes from ∇·u* after the fractional step)
div_u_star = 0.01 * np.sin(2*np.pi*z)
dt = 1e-3
 
# extrapolated pressure (paper Eq. 15)
P_n_minus_1 = np.zeros(Nz)
P_n = 0.01 * np.cos(2*np.pi*z)
P_hat = 2 * P_n - P_n_minus_1
dP_hat = np.gradient(P_hat, dz)
 
# correction term (RHS second piece of Eq. 16)
correction = (1 - rho_ref / rho) * dP_hat
rhs = rho_ref / dt * div_u_star + np.gradient(correction, dz)
 
# constant-coefficient Laplacian solve (FFT, periodic BC)
k = 2*np.pi * np.fft.fftfreq(Nz, dz)
laplacian_eigenvalues = -(k**2)
laplacian_eigenvalues[0] = 1.0  # remove null space
P_new = np.real(ifft(fft(rhs) / laplacian_eigenvalues))
P_new -= P_new.mean()
 
# size of the correction term
correction_magnitude = np.max(np.abs(1 - rho_ref / rho))
print(f"max |1 - rho_0/rho|  = {correction_magnitude:.4f}")
print(f"P_new range           = [{P_new.min():.3e}, {P_new.max():.3e}]")

The output is:

max |1 - rho_0/rho|  = 0.0117
P_new range           = [-1.572e-06, 1.572e-06]

ρ_min / ρ_ref: 0.978ρ_max / ρ_ref: 1.000x-range:

The splitting scheme rewrites the variable Poisson coefficient as constant ρ₀ plus a residual proportional to |1 − ρ₀/ρ|. For water at ΔT = 50 K the correction stays under 3% — that is why a constant-coefficient FFT solver can replace an iterative variable solver without losing accuracy.

The default slider is water across 50 K. The correction stays inside 1%. Switch to the two-fluid range and the same term moves toward 0.5. The code's ΔT = 30 K print of 0.0117 is a narrower window than that 50 K slider.

The Gray–Giorgini water limit is 1.25 K. The code's ΔT = 30 K is $30/1.25 = 24$ times that limit. The air 10% line is 28.6 K, so 30 K just clears it. A PV/T gap of 50 K is 40 times the water limit. Over that 50 K, density moves about 2% and viscosity about 30%.

The linear coefficient in the density polynomial is $3.736\times 10^{-4}$ . Times 30 K it is $3.736\times 10^{-4}\times 30 = 0.0112$ , next to the printed $\max|1-\rho_0/\rho| = 0.0117$ . The quadratic coefficient $3.98\times 10^{-6}$ times $(30/2)^2 = 225$ is about $9.0\times 10^{-4}$ . Most of 0.0117 is the linear term. With $T_0 = 313\,\mathrm{K}$ , $dt = 10^{-3}$ and a divergence amplitude of 0.01, the pressure range is $\pm 1.572\times 10^{-6}$ .

The two-fluid example is $\rho_{\min} = 1$ , $\rho_{\max} = 1000$ . Then $1 - 1/1000 = 0.999$ . That is two orders away from the water value 0.0117. When the pressure solve is 60–80% of a step, keeping the FFT is how the paper's 2.5–5× shows up. Watching only $\mathrm{CFL}_{\max} = 0.2$ misses $\mathrm{VSL} = \nu\Delta t/\Delta x^2$ , which crosses 0.03 first. The theory allows 0.3, and the higher-order statistics are already off by 3.5%.

The sentence that the correction at ΔT = 30 K stays inside 1.2% points at the same printed 0.0117. The same code with $\rho_{\min} = 1$ and $\rho_{\max} = 1000$ jumps to 0.999. That is why a two-fluid solver using this split still asks for a small $\Delta t$ .

With $N_z = 64$ and $L_z = 1$ , $dz = 1/64$ . The extrapolated pressure is $\widehat P = 2P^n - P^{n-1}$ . In the listing the amplitude of $P^n$ is 0.01 and $P^{n-1}$ is 0, so the amplitude of $\widehat P$ is 0.02. The printed 0.0117 is the coefficient on that gradient: the explicit correction is 1.17% of it. Table 2 spans $Ra$ from $10^6$ to $10^8$ . Loosening the VSL ceiling from 0.03 to 0.3 on that range shifts the higher-order statistics by 3.5%.

Three pitfalls during verification#

Table 2 (air DHC, $Ra = 10^6$ to $10^8$ ) bites in three places.

Choosing $\hat\rho_0$ — Dodd–Ferrante uses the minimum density, so $(1-\hat\rho_0/\hat\rho)$ stays at least 0 and below 1. The paper says a reference-temperature density gives almost the same NOB answer, but a large ΔT spends the stability margin.
CFL and VSL together — explicit Adams–Bashforth means $\mathrm{CFL}_{\max} = 0.2$ is not the first limit. $\mathrm{VSL} = \nu\Delta t/\Delta x^2 < 0.03$ is the practical ceiling. Theory allows 0.3, and higher-order statistics drift by 3.5%.
Arakawa convective form — momentum and kinetic energy have to be conserved together at the discrete level. Plain central differences leak energy, and long-time averages tilt at large Ra.

Whether to pin rho_0 at the minimum density#

Dodd–Ferrante pins $\hat\rho_0$ at the minimum density so that $(1-\hat\rho_0/\hat\rho)$ stays at least 0 and below 1. The paper says a reference-temperature density is close enough for NOB. Past ΔT = 30 K that gap spends the margin. The printed 0.0117 is a 1D slab with the minimum density. It is not the bill for a 2D DNS.

The 2.5–5× figure assumes the pressure solve is 60–80% of the step, and that an FFT has a periodic direction to eat. Explaining the asymmetric boundary layer by a 2% density change misses the 30% viscosity change. The temperature where that assumption closes is not the OB window in natural convection at Rayleigh number 1708. For water it is already closed at 1.25 K.