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

把 1/3 指数套到同一个 Ra 上，热流大约降到 149。要到 1e9，高度得接近 0.73 m。

加热底板，流体仍会静一阵。往上走的只有热。条纹状的涡旋是一起打开的。那道门上的数字是瑞利数 1708。热壁冷得有多快，要看后面的平板计算。

把 Ra 滑条降到 1708 以下，粒子停住，温度上下分开。推过 1708，彼此反向的卷胞就打开。

Ra = 1.0e+4

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

刚过门槛时环流很弱。再把滑条往上推，边界层变薄，胞变强。

浮力与粘性碰撞之处

自然对流是没有泵、只靠密度差推动的流动。底板变热，上面的流体膨胀变轻。轻的一团想浮上去。这就是浮力。

只有浮力的话，流动马上就开始。粘性用摩擦拽住上升的那一团。热扩散抹掉热团的温度，连浮力所靠的温差一起抹掉。

两道刹车胜过浮力，流体就静止。浮力胜过它们，对流才打开。

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

密度在每一项里都可变，方程就重了。布辛涅斯克近似只让密度在浮力项里随温度变化。其余地方的密度是常数。

体积热膨胀系数，是温度升高 1 度时密度下降的比例。

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

β 是这个系数，ρ 是密度，T 是温度。理想气体里 β 等于 1/T。

密度差于是换成温差。

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

ρ∞ 和 T∞ 是远处静止流体的密度和温度。这一行把浮力埋进竖直平板边界层的动量方程。

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 是重力加速度，ν 是运动粘度。右边第一项是浮力，第二项是粘性。能量方程用扩散带走温度。

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

α 是热扩散系数。动量和能量被浮力绑在一起，这就是自然对流的骨架。

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

方程无量纲化之后，浮力项收成一个数。格拉斯霍夫数。

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

Ts 是壁温，L 是特征长度。格拉斯霍夫数接过强迫对流里只靠长度给流动定位的雷诺数的工作。它量的是浮力把粘性压到什么程度。

对手不只有粘性。热扩散也在擦掉温差。第二道刹车以普朗特数进入，也就是动量扩散与热扩散之比 Pr = ν/α。乘积是瑞利数。

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

分母里同时有 ν 和 α。对流是否打开，不由 Gr 单独决定，而由 Ra 决定。浮力必须胜过两道刹车的乘积。

名为1708的临界值#

把从下方加热的一层放进线性稳定性分析，会得到一个临界值。上下都是固体壁时，Ra_c 是 1708。低于它，扰动被粘性和热扩散吃掉。流体静止，热量只靠导热传递。

越过 1708，反向旋转的卷胞排开。那就是贝纳尔胞。胞宽大约是层厚的两倍。前面的 Ra 滑条跨的就是这道门。

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

壁面冷得怎么样，努塞尔特数 Nu = hL/k 来回答。它是对流相对于导热的比。竖直等温平板有一条全区间都用的丘吉尔–朱公式。情形是 60 度的壁面立在 20 度的空气里，高 0.30 m。

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

Pr = 0.69Ra marker = 1e7.0

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

把 Pr 滑条往液态金属一侧拉，曲线下降；往油一侧推，曲线上升。把 Ra 标记推过 1e9，斜率从 1/4 换成 1/3。

印出来的数是 Ra 7.03e7、Nu 54.9、h 4.94、热流 197.6。膜温度是 40 加上 273.15，即 313.15 K，膨胀系数是它的倒数，约 0.003193。粘性与热扩散的乘积是 1.85e-5 乘 2.60e-5，也就是 4.81e-10。高度 0.30 m 的立方是 0.027。这些碎片遇上代码里的重力 9.81 和温差 40，就聚成印出的 7.03e7。Nu 54.9 乘 0.027 再除以 0.30，得到 4.941，与输出的 h 4.94 一致。4.94 乘 40 等于 197.6，热流那一行也闭上了。这个 Ra 除以 Pr 0.71，Gr 约为 9.90e7。它不能和 1708 并排放。

同一个 Ra 上，层流别名是 0.59 乘 Ra 的 1/4 次方，约 54.0。输出的 54.9 比它大约高出 0.9。湍流别名是 0.10 乘 Ra 的 1/3 次方，约 41.3。用这个 Nu 重算 h，大约 3.72，热流大约 149。这是 197.6 的大约四分之三。

Ra 要到 1e9，需要 7.03e7 的大约 14.2 倍。高度以立方进入，所以 0.30 m 得长到大约 0.73 m 才碰到那道门。在此之前，这块平板留在 1/4 指数一侧。同样的空气、同样的 40 度温差下，Ra 等于 1708 的长度大约是 8.7 mm。那个长度是封闭层的厚度。把 1708 读成 0.30 m 平板的门，几何就错了。

输出的 197.6 是每平方米的热流。宽度取得和高度一样，0.30 m，面积就是 0.09 m^2。整块板失去的热是 197.6 乘 0.09，大约 17.8 W。

只跟层流别名走的话，Nu 按 Ra 的 1/4 次方，Ra 按高度的立方，于是 h 按高度的 -1/4 次方移动。高度加倍时，h 变成 2 的 -1/4 次方，大约 0.841 倍。4.94 变成大约 4.15，热流降到大约 166。在 1/3 指数下，这个下降消失，h 不再跟着高度走。这次输出的 Ra 还不是那个指数。

努塞尔特关联式一束

每种形状有自己的式子。

竖直平板： $\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}$ 正比于 $\mathrm{Ra}^{1/3}$ ， $\mathrm{Ra}$ 正比于 $L^3$ ，于是 $L$ 从 $h = \mathrm{Nu}\,k/L$ 里消掉。湍流自然对流的传热系数不挂在尺寸上。上面这块 0.30 m 的平板，还没到能用这个抵消的高度。

把湍流指数套到7.03e7上，热流会被削掉#

把 1/3 套到这次输出上，Nu 54.9 会降到约 41.3，197.6 会削到约 149。这个指数把高度从式子里擦掉，0.30 m 也一起消失。层流一支里，改高度，h 还跟着动。先选错指数，这个变动会因错误的理由死掉。

1e9 要到大约 0.73 m 才来。胞已经打开的封闭层里，只留三个模态时这份热流还剩多少，把卷胞折成三个方程的 Lorenz 模型带着这个问题。