流体力学基础

浮力与自然对流：原理与诊断验证

从 Boussinesq 近似的密度误差出发，说明浮力在动量方程与湍流输运方程中的落点、封闭腔与开口边界的处理方式，并给出 Rayleigh-Bénard 临界值与热边界层网格的量化判据。

林工程师更新于 2026/9/22
CFD结果诊断与可信度验证物理建模与适用边界Boussinesq 浮力近似流体力学基础浮力与自然对流瑞利数能量收支浮力生成项进阶

浮力与自然对流：原理与诊断验证

浮力项写错位置或密度模型选得过粗，会让整个算例在正确的网格上解一个错误的方程。本文只讨论建模层面的三个决定：用 Boussinesq 还是变密度、浮力项落在动量方程与 kk 方程的哪一处、以及边界与域尺寸如何决定解的物理性。每条结论都给出可核算的数字，便于判断当前工况是否越界。自然对流算例最典型的失败方式不是发散，而是收敛到一个几乎不流动的解：残差漂亮下降、温度场光滑，但速度只有 10−4 m/s10^{-4}\,\mathrm{m/s} 量级，浮力被压力梯度完全抵消。诊断的第一步是确认解里确实存在浮力驱动的流动，再用相似准则、能量收支和独立基准三路证据交叉验证。

浮力项在动量方程中的写法

不可压 Boussinesq 求解器把浮力作为体积力加入动量方程，并把参考密度并入修正压力：

∂u∂t+(u⋅∇)u=−∇prgh+∇⋅(ν∇u)+[1−β(T−T0)]g\frac{\partial \mathbf{u}}{\partial t} + \left(\mathbf{u}\cdot\nabla\right)\mathbf{u} = -\nabla p_{rgh} + \nabla\cdot\left(\nu\nabla\mathbf{u}\right) + \left[1-\beta\left(T-T_0\right)\right]\mathbf{g}

其中 prgh=p−ρ0 g⋅xp_{rgh} = p - \rho_0\,\mathbf{g}\cdot\mathbf{x} 是修正压力，它把静水压从压力场中扣除。这一步是浮力算例最容易出错的地方：若求解器使用的是 prghp_{rgh}，那么压力边界必须给修正压力（开口处通常为 0），直接给绝对压力会让整个压力场平移一个 ρ0gH\rho_0 g H。对 H=2 mH = 2\,\mathrm{m} 的空气柱，这个偏差约 1.18×9.81×2=23 Pa1.18\times9.81\times2 = 23\,\mathrm{Pa}，看似不大，但对以 Pa 为量级的自然对流压差已经是致命量级。

湍流方程里的浮力生成项

标准 kk-ε\varepsilon 需要额外提供浮力生成项，否则热分层对湍流强度的影响完全丢失：

Gk=−μtρ σT∂ρ∂xigi=β μtσT g ∂T∂yG_k = -\frac{\mu_t}{\rho\,\sigma_T}\frac{\partial \rho}{\partial x_i}g_i = \beta\,\frac{\mu_t}{\sigma_T}\,g\,\frac{\partial T}{\partial y}

（取重力沿 −y-y 方向，σT\sigma_T 为湍流普朗特数，空气常取 0.85。）符号决定效果：热面向上时 ∂T/∂y<0\partial T/\partial y < 0，Gk>0G_k > 0，湍流被增强；稳定分层时 Gk<0G_k < 0，湍流被抑制。稳定分层下 GkG_k 可能超过耗散项，若不施加限制会导致 kk 出现负值，因此多数实现会对 GkG_k 做上限约束或改用 μt/μ\mu_t/\mu 的 Richardson 数修正。

判断是否需要打开浮力项，用梯度理查森数：

Rig=gβ ∂T/∂y(∂U/∂y)2Ri_g = \frac{g\beta\,\partial T/\partial y}{\left(\partial U/\partial y\right)^{2}}

∣Rig∣≫1|Ri_g| \gg 1 表示浮力主导湍流生成，必须打开 GkG_k。

流动区制与临界值

封闭腔的传热强弱由瑞利数决定，封闭腔换热常用 Globe-Dropkin 关联式（3×105<Ra<7×1093\times10^{5} < Ra < 7\times10^{9}）：

Nu=0.069 Ra1/3 Pr0.074Nu = 0.069\,Ra^{1/3}\,Pr^{0.074}

核算一个 H=0.5 mH = 0.5\,\mathrm{m}、ΔT=20 K\Delta T = 20\,\mathrm{K} 的空气腔体：Ra=9.81×3.33×10−3×20×0.53/(1.57×10−5×2.2×10−5)=2.37×108Ra = 9.81 \times 3.33\times10^{-3} \times 20 \times 0.5^{3} / (1.57\times10^{-5}\times2.2\times10^{-5}) = 2.37\times10^{8}，Ra1/3=618.6Ra^{1/3} = 618.6，Pr0.074=0.975Pr^{0.074} = 0.975，得 Nu=0.069×618.6×0.975=41.6Nu = 0.069\times618.6\times0.975 = 41.6，壁面热流 q=Nu kΔT/H=41.6×0.0262×20/0.5=43.6 W/m2q = Nu\,k\Delta T/H = 41.6\times0.0262\times20/0.5 = 43.6\,\mathrm{W/m^2}。

底部加热的水平流体层另有硬阈值：Rayleigh-Bénard 对流在 Rac=1708Ra_c = 1708 处起始（刚性上下边界）。低于该值只有纯导热，任何"解出对流胞"的结果都说明边界或初始条件引入了扰动。这条判据在验证求解器时非常有用。

建模选择清单

  1. 由 βΔT\beta\Delta T 判断 Boussinesq 是否成立，必要时切换变密度。
  2. 明确求解器的压力口径是绝对压力还是修正压力，边界条件与之匹配。
  3. 热分层显著时打开湍流方程的浮力生成项，并检查其符号与幅值。
  4. 封闭腔传热先核对 RaRa 落在哪条关联式的有效区间。
  5. 用 δT=H/Nu\delta_T = H/Nu 估算热边界层厚度，据此定首层厚度与增长率。
  6. 开口算例的域尺寸要能容纳卷吸流，且外边界给温度与总压。

网格与边界：热边界层才是限制

自然对流不需要解析黏性底层，但必须解析热边界层。其厚度可由 NuNu 估计：

δT≈HNu\delta_T \approx \frac{H}{Nu}

上例中 δT≈0.5/41.6=12 mm\delta_T \approx 0.5/41.6 = 12\,\mathrm{mm}，要求壁面法向至少 5 层单元落在 δT\delta_T 内，即首层厚度约 1～2 mm，增长率不超过 1.2。开口边界则应给总压与温度（entrainment 方向），并把计算域向外扩到 5H5H 以上，否则卷吸流被壁面人为阻断，NuNu 会偏高 10% 以上。

OpenFOAM 中一个可直接套用的 Boussinesq 物性字典：

// constant/thermophysicalProperties
thermophysicalProperties
{
    type            Boussinesq;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState Boussinesq;
    specie          pureMixture;

mixture
    {
        beta        3.33e-03;   // 1/K
        T0          300;        // K
        rho0        1.18;       // kg/m^3
        mu          1.85e-05;   // Pa.s
        Pr          0.71;
        Cp          1005;       // J/(kg.K)
        Hf          0;
    }
}

// constant/g
dimensions      [0 1 -2 0 0 0 0];
value           (0 -9.81 0);

时间步与统计窗口

自然对流中的羽流与 Bénard 胞以低频摆动，瞬时场永远不会"收敛"。判定办法是把采样窗口对半切开，比较前后半段的时均努塞尔数，相对差小于 2% 才认为统计量稳定。时间步要同时满足对流与浮力两个限制：

Δt≤min⁡(ΔxUmax⁡, ΔxgβΔT)\Delta t \le \min\left(\frac{\Delta x}{U_{\max}},\ \sqrt{\frac{\Delta x}{g\beta\Delta T}}\right)

取 Δx=2 mm\Delta x = 2\,\mathrm{mm}、ΔT=30 K\Delta T = 30\,\mathrm{K}，浮力限制给出 Δt≤0.002/0.98≈0.045 s\Delta t \le \sqrt{0.002/0.98} \approx 0.045\,\mathrm{s}，比 Umax⁡=0.5 m/sU_{\max} = 0.5\,\mathrm{m/s} 给出的 0.004 s0.004\,\mathrm{s} 宽松，因此实际由对流项主导。

故障模式与判定试验

现象 根因 判定试验
热羽流上升速度比经验值低 30% 以上 βΔT>0.1\beta\Delta T > 0.1 却仍用 Boussinesq 计算 βΔT\beta\Delta T，超过 0.1 换变密度求解器复算
压力场整体偏移约 ρ0gH\rho_0 gH 把绝对压力赋给了修正压力边界 检查压力边界类型与量纲，改回 prghp_{rgh} 口径
稳定分层区 kk 出现负值或爆掉 GkG_k 未做上限约束 监控 Gk/εG_k/\varepsilon，必要时限制或改用代数热通量模型
底加热腔在 Ra=800Ra = 800 时出现对流胞 初始扰动或侧壁边界引入非物理对流 用纯导热初值重算，确认 Ra<1708Ra < 1708 时 Nu→1Nu \to 1
NuNu 比关联式高 15% 且随域尺寸变化 开口域太小，卷吸被阻断 把域扩到 5H5H 以上重算，看 NuNu 是否回落
网格加密后 NuNu 持续漂移 热边界层内单元数不足 统计 δT\delta_T 内层数，补到 5 层以上
残差 10−610^{-6} 但 NuNu 比关联式低 40% 首层太厚，热边界层被数值导热抹平 统计壁面 y+y^{+} 与热边界层内单元数，要求不少于 5 层
速度量级只有 10−4 m/s10^{-4}\,\mathrm{m/s}、温度分层光滑 常密度求解器未加浮力源项，或参考密度填错 关闭重力重算，结果几乎不变说明浮力项未生效
封闭腔两壁热流相差 5% 以上 温度场未进入统计稳态 比较前后半窗时均 NuNu，或把时长延长 3 倍
对称双涡结构与实验单涡不符 忽略辐射，或侧壁绝热假设不成立 打开辐射模型并复核侧壁边界条件
加密网格后 NuNu 单调上升不收敛 解仍在网格相关区，或存在未解析的角点奇异 用三套网格做 Richardson 外推，看观测阶是否接近 2
减小时间步后时均量跳变 显式推进越过对流 CFL 限制 记录最大 CFL 数，压到 0.5 以下复算

Boussinesq 近似的误差到底有多大

Boussinesq 近似只保留密度随温度变化的一阶项：

ρ(T)=ρ0[1−β(T−T0)],β=−1ρ(∂ρ∂T)p\rho(T) = \rho_0\left[1-\beta\left(T-T_0\right)\right], \qquad \beta = -\frac{1}{\rho}\left(\frac{\partial \rho}{\partial T}\right)_p

对理想气体 β=1/T0\beta = 1/T_0。近似成立的常用判据是 βΔT≲0.1\beta\Delta T \lesssim 0.1，即密度变化不超过 10%。核算两档温差：

  • ΔT=20 K\Delta T = 20\,\mathrm{K}、$T_0 = 300,\mathrm{K}：：\beta\Delta T = 20/300 = 0.067$，密度变化 6.7%，Boussinesq 可用。
  • ΔT=100 K\Delta T = 100\,\mathrm{K}：βΔT=0.333\beta\Delta T = 0.333，密度变化 28%，此时必须改用理想气体状态方程的变密度求解器，否则热羽流上升速度会被系统性低估。

注意除密度外的物性（μ\mu、kk、cpc_p）在 Boussinesq 框架下仍被当作常数。若 ΔT\Delta T 跨越大范围，空气黏度在 300～500 K 间从 1.85×10−51.85\times10^{-5} 升到 2.67×10−5 Pa⋅s2.67\times10^{-5}\,\mathrm{Pa\cdot s}（+44%），这个变化与浮力项同阶，也需要一并评估。

用瑞利数确定该出现什么流动

竖直等温平板、高度 LL、温差 ΔT\Delta T 的空气自然对流由瑞利数控制：

Ra=gβΔTL3να=Gr Pr,Gr=gβΔTL3ν2Ra = \frac{g\beta \Delta T L^{3}}{\nu \alpha} = Gr\,Pr, \qquad Gr = \frac{g\beta \Delta T L^{3}}{\nu^{2}}

取 L=0.3 mL = 0.3\,\mathrm{m}、ΔT=30 K\Delta T = 30\,\mathrm{K}、膜温 Tf=300 KT_f = 300\,\mathrm{K} 的空气物性 ν=1.57×10−5 m2/s\nu = 1.57\times10^{-5}\,\mathrm{m^2/s}、α=2.2×10−5 m2/s\alpha = 2.2\times10^{-5}\,\mathrm{m^2/s}、β=1/Tf=3.33×10−3 K−1\beta = 1/T_f = 3.33\times10^{-3}\,\mathrm{K^{-1}}：

  • Gr=9.81×3.33×10−3×30×0.33/(1.57×10−5)2=1.07×108Gr = 9.81 \times 3.33\times10^{-3} \times 30 \times 0.3^{3} / (1.57\times10^{-5})^{2} = 1.07\times10^{8}
  • Ra=1.07×108×0.71=7.6×107Ra = 1.07\times10^{8} \times 0.71 = 7.6\times10^{7}

Ra=7.6×107Ra = 7.6\times10^{7} 位于竖板层流边界层区（转捩约在 Ra≈109Ra \approx 10^{9}），因此壁面平均努塞尔数应接近 Churchill-Chu 层流关联式：

Nu=0.68+0.670 Ra1/4[1+(0.492/Pr)9/16]4/9Nu = 0.68 + \frac{0.670\,Ra^{1/4}}{\left[1+\left(0.492/Pr\right)^{9/16}\right]^{4/9}}

代入 Ra1/4=93.5Ra^{1/4} = 93.5、Pr=0.71Pr = 0.71，分母为 [1+0.6939/16]4/9=1.303\left[1+0.693^{9/16}\right]^{4/9} = 1.303，得 Nu≈48.7Nu \approx 48.7。换算成换热系数 h=Nu k/L=48.7×0.0262/0.3=4.25 W/(m2⋅K)h = Nu\,k/L = 48.7 \times 0.0262/0.3 = 4.25\,\mathrm{W/(m^2\cdot K)}，热流密度 q=hΔT≈128 W/m2q = h\Delta T \approx 128\,\mathrm{W/m^2}。这就是写进验收表的硬数字：若 CFD 给出 q<60 W/m2q < 60\,\mathrm{W/m^2} 或 >250 W/m2> 250\,\mathrm{W/m^2}，先怀疑解被抹平或网格不足，而不是先改关联式。

有来流时先分离出浮力份额

存在强迫对流时，浮力是否必须保留由理查森数决定：

Ri=GrRe2=gβΔTLU2Ri = \frac{Gr}{Re^{2}} = \frac{g\beta\Delta T L}{U^{2}}

同一平板置于 U=1 m/sU = 1\,\mathrm{m/s} 的来流，Re=1×0.3/1.57×10−5=1.91×104Re = 1\times0.3/1.57\times10^{-5} = 1.91\times10^{4}，Ri=1.07×108/(1.91×104)2=0.29Ri = 1.07\times10^{8}/(1.91\times10^{4})^{2} = 0.29。经验分区为：Ri<0.1Ri < 0.1 可忽略浮力，0.1<Ri<100.1 < Ri < 10 属混合对流、惯性项与浮力项都要保留，Ri>10Ri > 10 以自然对流为主。Ri=0.29Ri = 0.29 落在混合对流区，此时关掉重力会带来 10% 量级的热流偏差，这是"结果与实验对不上"最常见的根因。

能量收支比残差更可靠

稳态下进出控制体的净焓流加壁面换热应闭合：

Q˙wall−m˙cp(Tˉout−Tˉin)−Q˙rad≈0\dot{Q}_{\text{wall}} - \dot{m}c_p\left(\bar{T}_{\text{out}}-\bar{T}_{\text{in}}\right) - \dot{Q}_{\text{rad}} \approx 0

封闭腔体内则简化为 Q˙热壁+Q˙冷壁≈0\dot{Q}_{\text{热壁}} + \dot{Q}_{\text{冷壁}} \approx 0（其余壁面绝热），工程上要求相对不平衡量低于 1%～2%。常见情形是残差已降到 10−610^{-6}，两壁热流却相差 8%，原因通常是温度场未进入统计稳态（羽流仍在低频摆动），或近壁热边界层内只有 1～2 个单元。

与独立证据对照

  • 关联式：竖板用 Churchill-Chu；水平热面朝上用 Nu=0.54Ra1/4Nu = 0.54Ra^{1/4}（104<Ra<10710^{4}<Ra<10^{7}）；封闭腔用 Globe-Dropkin。允许偏差 10%～15%，超出即视为模型或网格问题。
  • 网格无关性：三套网格的 NuNu 与最大速度，观测阶应接近 2；若低于 1，说明首层网格仍控制结果。
  • 能量收支：不平衡量低于 2%。
  • 物性敏感性：把 β\beta 从 1/Tf1/T_f 换成查表值，NuNu 变化应小于 5%；否则说明温度场对物性过于敏感，需改用变物性求解器。

后处理用独立脚本重算面积平均 NuNu，避免依赖求解器内部的归一化口径：

# 提取壁面热流，再用独立公式换算 Nu
postProcess -func "wallHeatFlux(patches=(hotWall))" -time 200:400
python3 - <<'PY'
k, L, dT = 0.0262, 0.3, 30.0     # W/(m.K), m, K
q = 128.4                        # 面积平均壁面热流, W/m^2
print("Nu =", q * L / (k * dT))  # 期望 ~48.7
PY

结论可信前必须回答的问题

  1. RaRa 与 RiRi 是否落在预期区间，模型选择与流动区制是否匹配？
  2. 壁面热流与关联式是否在 15% 以内一致？
  3. 能量收支不平衡量是否低于 2%？
  4. 前后半窗时均量差是否低于 2%？
  5. 三套网格的观测阶是否接近 2？
  6. 常物性与变物性两版计算的 NuNu 差异是否在容差内？

参考资料

  1. Gray D.D., Giorgini A., “The Validity of the Boussinesq Approximation for Liquids and Gases,” International Journal of Heat and Mass Transfer, 1976.
  2. de Vahl Davis G., “Natural Convection of Air in a Square Cavity: A Bench Mark Numerical Solution,” International Journal for Numerical Methods in Fluids, 1983.
  3. Henkes R.A.W.M., van der Vlugt F.F., Hoogendoorn C.J., “Natural-Convection Flow in a Square Cavity Calculated with Low-Reynolds-Number Turbulence Models,” International Journal of Heat and Mass Transfer, 1991.
  4. Ferziger J.H., Perić M., Street R.L., Computational Methods for Fluid Dynamics, 4th ed., Springer, 2020.
  5. Churchill S.W., Chu H.H.S., “Correlating Equations for Laminar and Turbulent Free Convection from a Vertical Plate,” International Journal of Heat and Mass Transfer, 1975.
  6. Globe S., Dropkin D., “Natural-Convection Heat Transfer in Liquids Confined by Two Horizontal Plates and Heated from Below,” Journal of Heat Transfer, 1959.
  7. Bejan A., Convection Heat Transfer, 4th ed., Wiley, 2013.
  8. ASME, Standard for Verification and Validation in Computational Fluid Dynamics and Heat Transfer, V&V 20-2009.