静止中被加热的流体 — 自然对流与瑞利数1708

在下方点火，流体却会安静一阵子。底部很热、顶部很冷，液体却保持静止，只把热量往上传。然后在某一瞬间，一整片条纹状的涡旋同时开启。本文追踪这个「临界点」为什么是瑞利数1708，并一直讲到决定受热壁面冷却好坏的努塞尔特关联式。最后你将亲自越过临界值，并用 Python 计算竖直平板的散热。

浮力与粘性碰撞之处#

自然对流（无外部泵、由密度差驱动的流动）是三种力的较量。加热底部，其上方的流体膨胀变轻。较轻的团块想要上浮。这就是浮力。

若只有浮力，流动本应立即开始。但有两样东西拖后腿。粘性通过与邻域的摩擦把上浮团块拽下来。热扩散迅速抹平热团块的温度，消除作为浮力原动力的温度差本身。

所以结论是一场简单的竞赛。两个刹车（粘性加热扩散）胜过浮力，流体就静止。浮力胜出，对流就开启。用一个数回答哪一方获胜，正是本文的目标。

布辛涅斯克近似 — 把密度变化压进一项#

密度随温度变化时方程会变得繁琐。布辛涅斯克近似（仅在重力项中把密度当作变量的简化）减轻了这个负担。在其他所有地方密度是常数，只有浮力项对温度作出响应。

定义体积热膨胀系数（温度每升高一度密度减少的比例）。

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

其中 $\beta$ 是膨胀系数， $\rho$ 是密度， $T$ 是温度。理想气体则干净地落到 $\beta = 1/T$ 。

用这个定义把密度差换成温度差。

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

$\rho_\infty$ 、 $T_\infty$ 是远处静止流体的密度与温度。这一行在竖直平板边界层的动量方程里埋下一个浮力源。

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$ 、 $v$ 是速度分量， $g$ 是重力加速度， $\nu$ 是运动粘度。右端第一项是浮力，第二项是粘性阻力。能量方程靠扩散搬运温度。

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

$\alpha$ 是热扩散率（热量扩散的快慢）。由浮力捆在一起的这一对（动量与能量）就是自然对流的骨架。

格拉斯霍夫与瑞利 — 两个无量纲数的分工#

无量纲化后浮力项凝成一个数。那就是格拉斯霍夫数。

\mathrm{Gr} = \frac{g\beta\,(T_s - T_\infty)\,L^3}{\nu^2} \sim \frac{\text{浮力}}{\text{粘性力}}

$T_s$ 是壁温， $L$ 是特征长度。格拉斯霍夫承担了强制对流中雷诺数所担的角色。它衡量浮力压倒粘性的程度。

但浮力的敌人不止粘性。热扩散也抹去温度差。要把这第二个刹车计入，需乘以普朗特数（动量扩散与热扩散之比， $\mathrm{Pr} = \nu/\alpha$ ）。结果就是瑞利数。

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

分母里同时坐着粘性 $\nu$ 与热扩散 $\alpha$ 。这是关键。对流是否点燃，由 Ra 而非 Gr 单独决定。浮力必须胜过两个刹车的乘积，流动才能存活。

名为1708的临界值#

把自下方加热的流体层放进线性稳定性分析，会出现一个干净的临界值。上下都是固体壁时，临界瑞利数是 $\mathrm{Ra}_c = 1708$ 。低于它，一切扰动都被粘性和热扩散吃掉而消亡。流体静止，热量仅以纯导热传递。

Ra 越过1708的那一刻，图案被开启：规则排列、反向旋转的卷胞，即贝纳尔胞。胞的宽度约为层厚的两倍。

在下面的模拟里亲手操作一下。

Ra = 1.0e+4

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

把 Ra 降到1708以下，粒子几乎冻结，温度场上下干净分层。越过1708，反向旋转的胞开启，粒子开始环流。再把 Ra 调高，环流变得猛烈，边界层变薄。

Python — 追踪竖直平板的努塞尔特数#

受热壁面冷却得有多好？答案是努塞尔特数（对流相对导热的传热比， $\mathrm{Nu} = hL/k$ ）。竖直等温平板有一个在全区间都适用的丘吉尔–朱关联式。来算一面立在20度空气中的60度壁。

import numpy as np
 
def rayleigh_number(Ts, Tinf, L, beta, nu, alpha, g=9.81):
    # 浮力与 (粘性 x 热扩散) 之比
    return g * beta * (Ts - Tinf) * L**3 / (nu * alpha)
 
def churchill_chu_vertical(Ra, Pr):
    # 竖直等温平板 -- 全区间有效
    denom = (1.0 + (0.492 / Pr)**(9 / 16))**(8 / 27)
    return (0.825 + 0.387 * Ra**(1 / 6) / denom)**2
 
# 以膜温为基准的空气物性
nu, alpha, Pr, k = 1.85e-5, 2.60e-5, 0.71, 0.027
Ts, Tinf, L = 60.0, 20.0, 0.30          # 60C 壁, 0.3 m 高
beta = 1.0 / (0.5 * (Ts + Tinf) + 273.15)  # 理想气体 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)                     # 单位面积散热
 
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

$\mathrm{Ra} \approx 7\times10^7$ ，还在层流区。每平方米平板散失约198瓦。改变温差或高度，Ra 立刻移动。

下面的图把同一关联式沿 Ra 轴铺开。

Pr = 0.69Ra marker = 1e7.0

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

降低 Pr（朝液态金属一侧）曲线下移，升高 Pr（朝油一侧）则上移。把 Ra 标记拖过 $10^9$ ，斜率会从层流的 $\mathrm{Ra}^{1/4}$ 切换到湍流的 $\mathrm{Ra}^{1/3}$ 。

努塞尔特关联式一束#

每种几何形状都有自己的关联式。这里汇总常用的几个。

竖直平板： $\mathrm{Nu} = 0.59\,\mathrm{Ra}^{1/4}$ （层流， $10^4$ – $10^9$ ）， $\mathrm{Nu} = 0.10\,\mathrm{Ra}^{1/3}$ （湍流， $10^9$ – $10^{13}$ ）。全区间用丘吉尔–朱。
水平加热板，朝上： $\mathrm{Nu} = 0.54\,\mathrm{Ra}^{1/4}$ （ $10^4$ – $10^7$ ）， $0.15\,\mathrm{Ra}^{1/3}$ （ $10^7$ – $10^{11}$ ）。
水平加热板，朝下： $\mathrm{Nu} = 0.27\,\mathrm{Ra}^{1/4}$ 。
自下方加热的封闭腔：越过 $\mathrm{Ra}_c = 1708$ 就出现贝纳尔胞。

水平板的特征长度取面积除以周长 $L = A_s/P$ 。湍流区的指数 $1/3$ 藏着一层含义。因为 $\mathrm{Nu} \propto \mathrm{Ra}^{1/3}$ 且 $\mathrm{Ra} \propto L^3$ ，在 $h = \mathrm{Nu}\,k/L$ 中 $L$ 相消。湍流自然对流的传热系数与尺寸无关。

值得记住的#

浮力的敌人有两个：粘性与热扩散。所以对流的开始不由 Gr 单独决定，而由把两者相乘的 Ra 决定。
$\mathrm{Ra} < 1708$ 时流体静止（纯导热），越过它则贝纳尔胞开启。
$\mathrm{Nu}$ 在层流里随 $\mathrm{Ra}^{1/4}$ ，在湍流里随 $\mathrm{Ra}^{1/3}$ 。后者中传热系数不再依赖平板尺寸。