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

Gray 和 Giorgini 在 1976 年算出两个数字。常温空气下,Oberbeck–Boussinesq 近似保持 10% 误差的最大温差是 28.6 K。水的极限只有 1.25 K。太阳能塔里空气往返几百度,光伏热水板里水中午就升温 50 K,这些场景一旦走出 OB 适用窗口,误差就开始累积。Demou、Frantzis 和 Grigoriadis(2018)直接抛弃该近似:每一项都用温度相关物性,而压力 Poisson 方程依然保持常系数。本文沿着那条让一切串起来的一行替换走一遍。

一行总结#

作者: 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}

帽子号表示用参考温度值归一化后的无量纲物性。同样的可变性也出现在能量方程中。

直接看物性曲线#

下图是水的四项关键物性随温度的变化。拖动热端、冷端滑块,看看每一条曲线摆动的幅度有多不同。

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%。粘度同样区间却变化几十个百分点。只用平均密度走 OB 近似,会完全错过可变粘度造成的不对称边界层。Demou 第 §3.2 节展示了水的 RB 对流中,正是这种不对称把 Nusselt 数分布扭歪。

可变系数 Poisson 的陷阱#

分数步法直接用,压力修正阶段的无散度条件要求:

\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 直接求解器(FDS)只接受常系数 Laplace 算子,无法直接用,只能退回多重网格或迭代法。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}$ ,在第 $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]

左侧 Laplace 算子每步相同。对两个方向做 FFT,问题降为沿 $z$ 方向的 Helmholtz 方程组,直接解。只有右侧修正项需要显式更新。

整套方案的成败取决于 $(1 - \hat\rho_0/\hat\rho)$ 的大小。两流体情况下该项接近 1,水的 NOB 对流中却极小。下图把两个情形并列。

ρ_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 — 这就是为什么分裂方案在两流体下需要小 $\Delta t$ ,而在 NOB 下几乎免费。

NumPy 走完一个循环#

把论文 §2.2–2.3 的流程压缩成 1D RB 板的一步。不做真正的 2D 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]

水在 ΔT = 30 K 下修正项被限制在 1.2% 以内。把同一份代码用在两流体( $\rho_{\min} = 1, \rho_{\max} = 1000$ )上,这个数字会跳到 0.999 — 也是两流体求解器在同样分裂方案下必须缩小 $\Delta t$ 的根本原因。

验证时遇到的三个坑#

复现论文 Table 2(空气 DHC, $Ra = 10^6\!-\!10^8$ )时,有三处容易绊倒。

$\hat\rho_0$ 的选取 — Dodd–Ferrante 原作用最小密度,这样 $(1-\hat\rho_0/\hat\rho)$ 锁在 $[0, 1)$ 内。论文说在 NOB 自然对流里用参考温度密度结果几乎一样,但大 ΔT 下两者差会蚕食稳定裕度。
CFL 与 VSL 都要看 — 显式 Adams–Bashforth 下只盯 $\mathrm{CFL}_{\max} = 0.2$ 不够,粘性限制 $\mathrm{VSL} = \nu\Delta t/\Delta x^2 < 0.03$ 通常先破。理论稳定区间可到 0.3,但 0.03 以外高阶统计会漂移 3.5%。
Arakawa 型对流项 — 必须让动量与动能在离散层面同时守恒。普通中心差分会偷偷漏能量,大 Ra 下时间平均统计会悄悄倾斜。

收尾#

NOB 自然对流的意义在于让可变粘度与可变热导率塑造出不对称边界层。Demou 等把这种可变性留在动量与能量方程中,只在压力 Poisson 处用分裂技巧抽掉。结果是一行替换加一项显式修正,而这一行让 FFT 直接求解器仍能上场 — 正是这点把 3D NOB DNS 从论文里搬到机器上。