始于锅底的混沌 — 瑞利-贝纳对流与洛伦兹方程

1900年，亨利·贝纳从底部加热薄油层，看见液体自己裂成六角形蜂窝。没有人搅拌。十六年之后，瑞利勋爵把这道门槛写成一个数。又过了四十七年，洛伦兹从同一组方程里取出了混沌。

门槛是瑞利数1708。下面是导热，上面是对流卷。把对流卷切成三个方程，剩下加热参数ρ。一过24.74，定点就抓不住轨迹。

把Ra滑条拉到1708以下，粒子几乎停住；推过1708，彼此反向的对流卷就点着。

1708以下，温度只堆成水平层。过了门槛再把滑条往上推，同一对卷转得更快。

静止的液体为何开始运动

下热上冷。被加热的下层膨胀变轻。轻的在下，重的在上。浮力想把这叠翻过来。

它不会马上涌起。黏性拽住流动，热扩散把上升热团的热量漏给四周，抬它的浮力就死了。两样都是刹车。

浮力赢了，就出现环流，也就是对流卷。刹车赢了，液体不动，热只走导热。浮力若在半路上被花掉，卷根本长不出来。烟道降低1 m、抽力消失7 Pa的情形就是这笔预算。

瑞利数 — 翻转之力与拽住之力之比

把这场拉锯压成一个无量纲数，就是瑞利数。

$g$ 是重力加速度，$\beta$ 是热膨胀系数（温度升高时体积增大的程度），$\Delta T$ 是上下温差，$L$ 是液层厚度，$\nu$ 是运动黏度，$\alpha$ 是热扩散率。分子是要翻转的浮力，分母是拽住的黏性与扩散。

分量大的是 $L^3$。层一变厚，瑞利数按立方增大。薄油膜很难对流。深海和地幔在很小的温差下就会动。锅只加深一点，门槛就近了一大截，原因在这个立方。

临界瑞利数1708 — 对流点燃的门槛#

瑞利的线性稳定性分析，边界一换数字就换。上下自由的理想边界，门槛是 $\mathrm{Ra}_c = 27\pi^4/4 \approx 657.5$。实验里常见的刚性壁，$\mathrm{Ra}_c \approx 1708$。

$\mathrm{Ra} < \mathrm{Ra}_c$ 时，小扰动被黏性和扩散吃掉。液体回到静止。一过 $\mathrm{Ra} > \mathrm{Ra}_c$，某个波长的扰动自己生长。那就是对流卷的种子。

1708除以657.5，大约是2.6。同样的流体、同样的厚度，刚性壁要求温差大约再大2.6倍。改用厚度来补，是立方根，$2.6^{1/3} \approx 1.37$ 倍就够。把自由滑移的盖子换成刚性壁，代价就是这个倍数。

对流卷的形状与努塞尔数

刚过门槛的卷会挑一个波长。临界波数是 $k_c \approx 3.117/L$。一个波长的长度是 $2\pi/k$，除以 $3.117$ 大约是 $2.02L$。并排的两个卷占满这个波长，单个卷的宽度大约是一半，靠近 $L$。说它接近正方形，就是这次除法。在最低 Ra 点燃的模态，波长正好是它。

热多送出去多少，用**努塞尔数（总传热与纯导热之比）**量。

$q_\text{total}$ 是实际热流密度，$q_\text{conduction}$ 是没有对流时的导热热流。只有导热时 $\mathrm{Nu}=1$。对流变强以后，实验大约是 $\mathrm{Nu} \sim \mathrm{Ra}^{1/3}$。只读这个指数，Ra 变成8倍时 Nu 变成2倍。门槛刚过的地方用不上。那里 Nu 还靠近1。

洛伦兹缩成三个方程的对流

1963年，爱德华·洛伦兹只留下温度场和流场最粗的几个模态。这是伽辽金截断。偏微分方程塌成三个常微分方程。

$x$ 是对流卷的环流强度，$y$ 是上升柱和下降柱的温差，$z$ 是竖直温度分布的扭曲。$\sigma$ 是普朗特数（动量扩散与热扩散之比），$\beta$ 是几何常数，$\rho$ 就是 $\mathrm{Ra}/\mathrm{Ra}_c$。

把 $\rho$ 调高，等于把锅烧得更猛。$\rho<1$ 时轨迹奔向原点，那是导热。$1<\rho<24.74$ 时收到偏离原点的定点，那是定常对流卷。一过 $\rho>24.74$，轨迹不再落在任何定点上，只在两片翼之间走。这就是混沌。

定常点是把三个导数都放成0的地方。$\sigma \neq 0$ 则 $y=x$。$x \neq 0$ 则 $z=\rho-1$。剩下的式子是 $x^2 = \beta(\rho-1)$。代码里 $\beta = 8/3$、$\rho = 15$ 时，$z=14$，$x^2 = 112/3 \approx 37.3$。中间那行输出 steady roll，意思是尾段4000个点里 $x$ 的跨度窄于0.1。$\rho=28$ 时 $z=27$、$x^2=72$，两片翼就停在那里。输出 chaos 的意思是，那段尾巴没能通过0.1这道门。

Python — 积分对流的三副面孔#

四阶龙格–库塔。$\sigma=10$，$\beta=8/3$，时间步0.01，共20000步。尾巴是最后4000个点。

import numpy as np
 
SIGMA, BETA = 10.0, 8.0 / 3.0
 
def lorenz_deriv(state, rho):
    x, y, z = state
    return np.array([
        SIGMA * (y - x),        # 卷环流强度的变化
        x * (rho - z) - y,      # 上升/下降温差的变化
        x * y - BETA * z,       # 竖直温度扭曲的变化
    ])
 
def integrate_rk4(rho, dt=0.01, n=20000, s0=(0.6, 0.4, 0.2)):
    s = np.array(s0, dtype=float)
    traj = np.empty((n, 3))
    for i in range(n):
        k1 = lorenz_deriv(s, rho)
        k2 = lorenz_deriv(s + 0.5 * dt * k1, rho)
        k3 = lorenz_deriv(s + 0.5 * dt * k2, rho)
        k4 = lorenz_deriv(s + dt * k3, rho)
        s = s + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
        traj[i] = s
    return traj
 
def convection_regime(rho):
    tail = integrate_rk4(rho)[-4000:]          # 丢弃暂态后的尾段
    spread = tail[:, 0].max() - tail[:, 0].min()
    if abs(tail[-1, 0]) < 1e-3 and rho < 1.0:
        return "conduction"                     # 收敛到原点
    if spread < 0.1:
        return "steady roll"                    # 收敛到固定点
    return "chaos"                              # 在两片翼之间振荡
 
for rho in (0.8, 15.0, 28.0):
    print(f"rho={rho:5.1f} -> {convection_regime(rho)}")

打出来是三行。

rho=  0.8 -> conduction
rho= 15.0 -> steady roll
rho= 28.0 -> chaos

把ρ滑条放在24.74以下，轨迹会掉进一个黄色定点。推过去，蝴蝶翅膀就张开。reset 把起点挪开一点。开头叠在一起的两条轨迹，很快就分开。

过了1708就当计算结束，会失败#

厚度加倍，$L^3$ 让 Ra 变成8倍。一层刚刚贴在1708上的流体，加深之后 $\rho = \mathrm{Ra}/\mathrm{Ra}_c = 8$。洛伦兹的定常卷区间在1和24.74之间。8在里面。要到混沌，ρ得升到24.74；温差不变的话，厚度要到 $24.74^{1/3} \approx 2.9$ 倍。把刚过门槛的卷写成混沌，错的就是这段间隔。

这段间隔是把 $\beta$ 当成常数、只加大 $\Delta T$ 做出来的。温差一过30 K，这个常数会不会先变成假的？打开这个温度上布辛涅斯克失效的那一边，就能看到1708和24.74哪一个先把数字丢掉。