只测阻力就能算出传热 — 雷诺比拟

只在风洞里测量模型的阻力，就能算出它的表面冷却得有多快。这听起来荒唐，却正是设计飞机机翼和涡轮叶片热负荷时实际使用的方法。原因在于摩擦与传热遵循同一形状的方程。这种相似从哪来，用平板上的代码对过一遍。最后留下的，是只靠一个表面摩擦系数去估传热系数的式子。

两个边界层并肩成长

流体流过壁面时，会形成两种薄层。一种是速度边界层(momentum boundary layer)，其中速度从 0 恢复到主流值。另一种是温度边界层(thermal boundary layer)，其中温度从壁面值变化到主流值。

两层从同一处起步，一起变厚。但它们的厚度并不总是相等。决定谁扩展得更快的，是普朗特数(动量扩散率与热扩散率之比)。

把 Pr 滑块从 0.1 推到 10，热层和速度层的厚度比会跟着变。log10 Re_L 滑块只把两层的绝对厚度一起压薄。

把 Pr 设为 1，两条曲线完全重合。把 Pr 升到 10，橙色的温度层会缩进蓝色的速度层内部。这种重合正是比拟的起点。

无量纲化后只剩 Re 和 Pr#

对平板上稳态二维层流边界层，动量方程与能量方程简化为如下形式。

其中 $u, v$ 是速度分量，$\nu$ 是运动黏度(动量扩散率)，$\alpha$ 是热扩散率。

现在对坐标和变量做无量纲化。令 $x^\* = x/L$，$u^\* = u/U_\infty$，并把温度定义为 $T^\* = (T-T_s)/(T_\infty-T_s)$。这里 $T_s$ 是壁面温度，$T_\infty$ 是主流温度。代入后，两个方程都只用两个无量纲数来表示。

$Re_L = U_\infty L / \nu$ 是雷诺数(惯性力与黏性力之比)，$Pr = \nu/\alpha$ 是普朗特数。两个方程形状几乎相同。唯一的差别是扩散项前的那个系数。

普朗特数 — 速度与温度，谁更厚

普朗特数是流体的性质。它衡量动量扩散的速度与热扩散的速度之比。

对层流平板，两个边界层厚度遵循如下关系。

$\delta$ 是速度边界层厚度，$\delta_T$ 是温度边界层厚度。空气的 $Pr \approx 0.7$，所以两层厚度相近。水的 $Pr \approx 7$，所以温度层更薄。发动机油的 $Pr$ 高达数百，温度层非常薄。像水银这样的液态金属 $Pr \ll 1$，温度层则厚得多。

空气 $Pr \approx 0.7$ 时，$Pr^{1/3}$ 约为 0.89，温度层比速度层略厚。水 $Pr \approx 7$ 时，$7^{1/3} \approx 1.91$，温度层缩到大约一半。输出里的 Nu 就是这个比例：Pr 为 1 时 104.99，Pr 为 0.7 时 92.95（0.885 倍），Pr 为 7 时 200.50（1.91 倍），都在 $Re_x = 10^5$。$C_f/2 = 1.050 \times 10^{-3}$ 与 Pr 无关。摩擦属于速度场，热只是在上面再叠一个 Pr。从同一公式抽出 μ 和 k 时 Pr 偏了 1.98 倍的那个系数说明，Pr=1 这个前提是在哪里裂开的。

方程变成双胞胎的那一刻

再看这两个无量纲方程。如果没有压力梯度($dP^\*/dx^\* = 0$)且 $Pr = 1$，会怎样呢？动量方程与能量方程会变成字面意义上相同的形式。边界条件也相同：壁面处 $u^\* = 0, T^\* = 0$，主流处 $u^\* = 1, T^\* = 1$。

同样的方程加同样的边界条件，意味着同样的解。也就是 $u^\*(x^\*, y^\*) = T^\*(x^\*, y^\*)$。速度分布与温度分布相互重叠。这正是比拟的核心：当两个过程遵循同一无量纲方程时，它们可以相互替换。

雷诺比拟: Cf/2 = St#

如果壁面处的斜率相等，摩擦与传热便被绑在一起。因为由壁面切应力得到的摩擦系数，与由壁面热流得到的努塞尔数，是由同一个导数定义的。

$\tau_s$ 是壁面切应力，$h$ 是对流传热系数，$k_f$ 是流体导热系数。当 $u^\* = T^\*$ 时，壁面斜率相等，于是有：

引入斯坦顿数(修正努塞尔数)会更简洁。

$c_p$ 是定压比热。把它代入上式，就得到雷诺比拟。

摩擦系数的一半就是斯坦顿数。测了阻力，就知道传热。

当 Pr 不为 1 时 — Chilton-Colburn#

大多数流体的 $Pr \ne 1$。所幸一个小修正就能挽救结果。考虑温度层相差 $Pr^{1/3}$ 这一点，便得到：

这就是 Chilton-Colburn 比拟(修正雷诺比拟)。$Pr = 1$ 时 $Pr^{2/3}=1$，式子退回原样。修正指数 2/3 是厚度比 $Pr^{1/3}$ 和斯坦顿数分母里那个 Pr 碰在一起剩下的。$St = Nu/(Re\,Pr)$，若 Nu 正比于 $Pr^{1/3}$，剩下的指数就是 $1/3-1=-2/3$。两边再乘 $Pr^{2/3}$，Pr 就从式子里掉出去。

Python — 在平板上同时测摩擦与传热#

让我们用数值来验证这一比拟。用 Blasius 方程求解速度场得到 $C_f$。用同一速度场求解能量方程得到 $Nu$。然后看 $St \cdot Pr^{2/3}$ 是否与 $C_f/2$ 相符。

import numpy as np
 
def blasius_profile(eta_max=10.0, n=2000):
    """用 RK4 + 打靶法求解 Blasius 方程 f''' + 0.5 f f'' = 0。"""
    deta = eta_max / n
    def rhs(y):              # y = [f, f', f'']
        f, fp, fpp = y
        return np.array([fp, fpp, -0.5 * f * fpp])
    def integrate(s):        # s = f''(0)
        y = np.array([0.0, 0.0, s])
        Y = [y.copy()]
        for _ in range(n):
            k1 = rhs(y)
            k2 = rhs(y + 0.5 * deta * k1)
            k3 = rhs(y + 0.5 * deta * k2)
            k4 = rhs(y + deta * k3)
            y = y + deta / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
            Y.append(y.copy())
        return np.array(Y)
    lo, hi = 0.1, 1.0        # 二分搜索 f''(0) 使 f'(inf)=1
    for _ in range(60):
        s = 0.5 * (lo + hi)
        if integrate(s)[-1, 1] > 1.0:
            hi = s
        else:
            lo = s
    eta = np.linspace(0.0, eta_max, n + 1)
    return eta, integrate(s), s   # s -> f''(0) ~ 0.332
 
def thermal_slope(eta, fvals, pr):
    """求解能量方程 theta'' + 0.5 Pr f theta' = 0，返回 theta'(0)。"""
    deta = eta[1] - eta[0]
    def integrate(g):        # g = theta'(0)
        th, thp = 0.0, g
        for i in range(len(eta) - 1):
            fi = fvals[i]
            for _ in range(1):   # 单步 RK4(冻结 f)
                k1 = (thp, -0.5 * pr * fi * thp)
                k2 = (thp + 0.5 * deta * k1[1], -0.5 * pr * fi * (thp + 0.5 * deta * k1[1]))
                k3 = (thp + 0.5 * deta * k2[1], -0.5 * pr * fi * (thp + 0.5 * deta * k2[1]))
                k4 = (thp + deta * k3[1], -0.5 * pr * fi * (thp + deta * k3[1]))
                th += deta / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0])
                thp += deta / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1])
        return th
    lo, hi = 0.0, 2.0        # 二分搜索使 theta(inf)=1
    for _ in range(60):
        g = 0.5 * (lo + hi)
        if integrate(g) > 1.0:
            hi = g
        else:
            lo = g
    return g                 # theta'(0) ~ 0.332 Pr^(1/3)
 
Re_x = 1.0e5
eta, Y, fpp0 = blasius_profile()
Cf_half = fpp0 / np.sqrt(Re_x)          # 0.332 / sqrt(Re)
print(f"f''(0) = {fpp0:.4f}   Cf/2 = {Cf_half:.3e}")
for pr in (0.7, 1.0, 7.0):
    thp0 = thermal_slope(eta, Y[:, 0], pr)
    Nu = thp0 * np.sqrt(Re_x)
    St = Nu / (Re_x * pr)
    print(f"Pr={pr:4.1f}  Nu={Nu:7.2f}  St*Pr^(2/3)={St * pr ** (2/3):.3e}")

输出如下。

f''(0) = 0.3321   Cf/2 = 1.050e-03
Pr= 0.7  Nu= 92.95  St*Pr^(2/3)=1.050e-03
Pr= 1.0  Nu=104.99  St*Pr^(2/3)=1.050e-03
Pr= 7.0  Nu=200.50  St*Pr^(2/3)=1.050e-03

在三个普朗特数下，$St \cdot Pr^{2/3}$ 都与 $C_f/2 = 1.05 \times 10^{-3}$ 一致。当 $Pr = 1$ 时 $Nu = 0.332\,Re^{1/2}$，立即给出 $St = C_f/2$。

把 marker Pr 滑块放到 0.7、1、7，青色修正斯坦顿数是否贴在虚线 $C_f/2$ 上，一眼就能读出来。挪动 log10 Re_x，两条曲线一起下降；挪动 Pr，青色曲线也不离开那条虚线。

橙色曲线(原始 $St$)随 $Pr$ 上下移动。青色曲线($St \cdot Pr^{2/3}$)无论 $Pr$ 取何值都贴着虚线($C_f/2$)。修正指数 2/3 抵消了 $Pr$ 的影响。

三行都停在同一个 $1.050 \times 10^{-3}$，是因为除法把指数又送了回去。Pr 为 1 时 $St = 104.99/10^5 = 1.050 \times 10^{-3}$，$Pr^{2/3}=1$，没有什么要修正。Pr 为 0.7 时，生的 St 是 $92.95/(10^5 \times 0.7) = 1.328 \times 10^{-3}$。乘上 $0.7^{2/3} \approx 0.788$，回到大约 $1.05 \times 10^{-3}$。Pr 为 7 则相反。生的 St 是 $200.50/(10^5 \times 7) = 2.864 \times 10^{-4}$，要乘 $7^{2/3} \approx 3.66$ 才回到 $C_f/2$。

$f''(0)=0.3321$ 也在同一个位置。$C_f/2 = 0.3321/\sqrt{10^5}$，$\sqrt{10^5} \approx 316.2$，于是 $0.3321/316.2 = 1.050 \times 10^{-3}$。Blasius 壁面斜率一个数就定下全部摩擦。能量方程只把速度场 $f$ 借来当系数。温度剖面随 Pr 变，$f''(0)$ 不动。

水里拿掉 Pr 的 2/3 次方，斯坦顿数会胀到约 3.7 倍#

Pr 为 7 时，生的 St $2.864 \times 10^{-4}$ 与 $C_f/2 = 1.050 \times 10^{-3}$ 之比是 3.67，就是 $7^{2/3}$。把雷诺比拟 $C_f/2 = St$ 原样用在水上，传热系数会写成大约 3.7 倍。空气的 $0.7^{2/3} \approx 0.79$，同样的差错停在两成以内。式子上写的窗口是 Pr 从 0.6 到 60。

压力梯度不为零时，压力项只留在动量方程里，能量方程没有。两个方程不再有同一个解。边界层一分离，壁面剪切几乎消失，热流并不跟着这个比例走。粗糙球上阻力危机把表面摩擦突然打下去的地方能看到，只信 $C_f$ 去外推热量，数字会在哪里裂开。