[论文综述] 30 K 让Boussinesq撒谎时 — Demou(2018) NOB自然对流与压力分裂

ΔT=30 K是Gray–Giorgini水极限1.25 K的24倍，那时密度修正在百分之一附近。

Gray 和 Giorgini 在 1976 年留下两个数。常温空气里，Oberbeck–Boussinesq 停在 10% 误差内的最大温差是 28.6 K。水只有 1.25 K。太阳能塔里的空气来回几百 K，PV/T 的水会升高 50 K。Demou、Frantzis 与 Grigoriadis(2018)丢掉这个假设。物性跟着温度走，只有压力 Poisson 保持常系数。

拨动热壁和冷壁，同一个 ΔT 上密度曲线几乎躺平，粘性曲线明显倾斜。

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.

50 K 上密度大约动 2%。粘性在同一段里大约动 30%，这个差别做出不对称边界层。

论文便签

作者: A. D. Demou, C. Frantzis, D. G. E. Grigoriadis
期刊: International Journal of Heat and Mass Transfer 121 (2018) 1148–1162
DOI: 10.1016/j.ijheatmasstransfer.2018.04.144
贡献: 非 OB(NOB)自然对流里物性保持随温度变化，再用 Dodd–Ferrante 分裂把变系数压力 Poisson 换成常系数 Poisson 加显式修正。FFT 直接求解器留在原处。文中的加速是 2.5–5 倍。

Boussinesq 实际崩溃的地方#

标准 Boussinesq 同时放着三句话。密度只在重力项里变。其余物性是常数。粘性耗散丢掉。ΔT 变大后，粘度 $\mu(T)$ 、热导率 $k(T)$ 、比热 $C_p(T)$ 一起动。液体的密度变化可以很小，粘度变化却不小。

物性都留下，动量方程就带上变系数。

\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}

$\hat\rho, \hat\mu, \hat C_p, \hat k$ 都是用参考温度归一化的无量纲物性。同样的变化也进能量方程。

可变系数 Poisson 的陷阱#

分数步原样用下去，压力修正要这样保住无散。这一步排在 Chorin 把压力从动量里拿掉之后。

\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}.

左端乘着 $1/\hat\rho^{n+1}$ 。系数随空间和时间变，矩阵每步都换。FFT 直接求解器只收常系数拉普拉斯。剩下的路是多重网格或迭代求解器。NS 求解器 60–80% 的时间花在压力上时，这条退路就是账单。

Dodd–Ferrante 分裂 — 一行替换#

借来的技巧来自两相流(Dodd & Ferrante 2014)。压力梯度切成两段。

\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}

第一段用常数 $1/\hat\rho_0$ 乘新压力。第二段把外推压力 $\widehat P = 2P^n - P^{n-1}$ 乘上修正。取散度之后，Poisson 变成常系数。

\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]

左端拉普拉斯每步相同。两个方向做 FFT，剩下 z 方向的 Helmholtz，直接解。只有右端修正显式更新。

要看的大小是 $(1 - \hat\rho_0/\hat\rho)$ 。两相流里它接近 1，NOB 水对流里它很小。

NumPy 走完一个循环#

论文 §2.2–2.3 压成一条 1D RB 板。没有二维 DNS，只看一步压力修正。

import numpy as np
from numpy.fft import fft, ifft
 
# 网格: x 周期, z 壁面 (为简化只看一束 1D Helmholtz)
Nz = 64
Lz = 1.0
dz = Lz / Nz
z = (np.arange(Nz) + 0.5) * dz
 
# 温度依赖密度 (水的多项式, 简化)
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)
 
# 初始温度场: 线性 + 小扰动
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)
 
# 参考密度: 最轻 (最热) 的一侧
rho_ref = rho.min()
 
# 虚拟散度场 (完整求解器中由分数步后的 ∇·u* 给出)
div_u_star = 0.01 * np.sin(2*np.pi*z)
dt = 1e-3
 
# 外推压力 (论文 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)
 
# 修正项 (Eq. 16 右端第二项)
correction = (1 - rho_ref / rho) * dP_hat
rhs = rho_ref / dt * div_u_star + np.gradient(correction, dz)
 
# 常系数 Laplace 解 (FFT, 假定周期边界)
k = 2*np.pi * np.fft.fftfreq(Nz, dz)
laplacian_eigenvalues = -(k**2)
laplacian_eigenvalues[0] = 1.0  # 去掉零空间
P_new = np.real(ifft(fft(rhs) / laplacian_eigenvalues))
P_new -= P_new.mean()
 
# 查看修正项大小
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}]")

输出如下。

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.

默认滑块是水的 50 K。修正留在 1% 以内。拨到两相范围，同一项靠近 0.5。代码里 ΔT = 30 K 印出的 0.0117，比这个 50 K 窗口更窄。

Gray–Giorgini(1976)给水的极限是 1.25 K。代码的 ΔT = 30 K 是 $30/1.25 = 24$ 倍。空气的 10% 线是 28.6 K，30 K 也刚刚越过那条线。PV/T 的 50 K 是水极限的 40 倍。这 50 K 里密度大约动 2%，粘度大约动 30%。

密度多项式的线性系数是 $3.736\times 10^{-4}$ 。乘上 30 K 得到 $3.736\times 10^{-4}\times 30 = 0.0112$ 。它和输出的 $\max|1-\rho_0/\rho| = 0.0117$ 并排。二次系数 $3.98\times 10^{-6}$ 乘上 $(30/2)^2 = 225$ ，大约是 $9.0\times 10^{-4}$ 。0.0117 的大部分是线性项。 $T_0 = 313\,\mathrm{K}$ ， $dt = 10^{-3}$ ，散度振幅 0.01 时， $P$ 的范围是 $\pm 1.572\times 10^{-6}$ 。

两相例子是 $\rho_{\min} = 1$ 、 $\rho_{\max} = 1000$ 。 $1 - 1/1000 = 0.999$ 。和水的 0.0117 隔着两位。压力求解占一步的 60–80% 时，这个分裂把 FFT 留下来，就接到论文里的 2.5–5 倍。只看 $\mathrm{CFL}_{\max} = 0.2$ ， $\mathrm{VSL} = \nu\Delta t/\Delta x^2$ 会先越过 0.03。理论上到 0.3 仍稳定，高阶统计已经偏 3.5%。

ΔT = 30 K 的水把修正圈在 1.2% 里，指的就是输出 0.0117。同一段代码放到 $\rho_{\min} = 1$ 、 $\rho_{\max} = 1000$ ，这个值涨到 0.999。两相求解器在同一种分裂上仍要很小的 $\Delta t$ ，原因在这里。

网格 $N_z = 64$ 、 $L_z = 1$ ，于是 $dz = 1/64$ 。外推压力是 $\widehat P = 2P^n - P^{n-1}$ ，代码里 $P^n$ 的振幅是 0.01， $P^{n-1}$ 是 0，所以 $\widehat P$ 的振幅是 0.02。输出 0.0117 是乘在这个梯度上的系数。显式修正是外推梯度的 1.17%。Table 2 的 $Ra$ 从 $10^6$ 到 $10^8$ ，跨两个数量级。在这段里把 VSL 上限从 0.03 放到 0.3，高阶统计会偏 3.5%。

验证时遇到的三个坑

复现 Table 2(空气 DHC， $Ra=10^6\!-\!10^8$ )时，有三处绊脚。

$\hat\rho_0$ 的选择 — Dodd–Ferrante 原文用最小密度。这样 $(1-\hat\rho_0/\hat\rho)$ 留在 0 与 1 之间、不到 1。论文写 NOB 里用参考温度密度结果也接近，但大 ΔT 会吃掉稳定裕度。
CFL 和 VSL 都要 — 显式 Adams–Bashforth 下，只看 $\mathrm{CFL}_{\max} = 0.2$ 会让粘性极限先破。 $\mathrm{VSL} = \nu\Delta t/\Delta x^2 < 0.03$ 是实用上限。理论上到 0.3 仍稳定，高阶统计会偏 3.5%。
Arakawa 型对流项 — 动量和动能要在离散层同时守恒。普通中心差分会漏能量，大 Ra 下时间平均悄悄歪掉。

是否把 ρ₀ 固定在最小密度

Dodd–Ferrante 把 $\hat\rho_0$ 放在最小密度上，为的是让 $(1-\hat\rho_0/\hat\rho)$ 处在 0 与 1 之间、不到 1。论文说 NOB 里用参考温度密度也接近。ΔT 超过 30 K，这个差别会削掉裕度。输出 0.0117 是用最小密度的一维板读到的数，不是二维 DNS 的账单。

2.5–5 倍说的是压力时间占 60–80% 的时候。只有 FFT 吃得到周期方向，这个倍数才出现。用密度 2% 解释边界层不对称，会漏掉粘度 30%。这个假设关上的温度，和 Rayleigh 数 1708 的自然对流里的 OB 窗口不是一回事。对水来说，1.25 K 就已经关上了。