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

只在风洞里测量模型的阻力，就能算出它的表面冷却得有多快。这听起来荒唐，却正是设计飞机机翼和涡轮叶片热负荷时实际使用的方法。原因在于摩擦与传热遵循同一形状的方程。本文推导这种相似从何而来，并在平板上用代码直接验证。读到最后，你将掌握一个仅凭表面摩擦系数就能估算传热系数的公式。

两个边界层并肩成长#

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

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

请操作下面的模拟。

log₁₀ Re_L: 5.0  (Re_L = 1.0e+5)Prandtl number Pr: 1.00  (δ_T/δ = 1.00)C_f/2 = 1.050e-3  |  St = 1.050e-3
St · Pr^(2/3) = 1.050e-3  ← Pr≈1: C_f/2 = St

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

无量纲化后只剩 Re 和 Pr#

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

u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{dP}{dx} + \nu \frac{\partial^2 u}{\partial y^2}

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

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

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

u^\* \frac{\partial u^\*}{\partial x^\*} + v^\* \frac{\partial u^\*}{\partial y^\*} = -\frac{dP^\*}{dx^\*} + \frac{1}{Re_L}\frac{\partial^2 u^\*}{\partial y^{\*2}}

u^\* \frac{\partial T^\*}{\partial x^\*} + v^\* \frac{\partial T^\*}{\partial y^\*} = \frac{1}{Re_L \, Pr}\frac{\partial^2 T^\*}{\partial y^{\*2}}

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

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

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

Pr = \frac{\nu}{\alpha} = \frac{\text{动量扩散率}}{\text{热扩散率}}

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

\frac{\delta}{\delta_T} \approx Pr^{1/3}

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

方程变成双胞胎的那一刻#

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

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

雷诺比拟: Cf/2 = St#

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

C_f = \frac{\tau_s}{\tfrac{1}{2}\rho U_\infty^2}, \qquad Nu_L = \frac{hL}{k_f} = \left.\frac{\partial T^\*}{\partial y^\*}\right|_{y^\*=0}

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

\frac{C_f}{2} Re_L = Nu_L

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

St \equiv \frac{h}{\rho U_\infty c_p} = \frac{Nu_L}{Re_L \, Pr}

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

\frac{C_f}{2} = St \qquad (Pr \approx 1,\ dP/dx \approx 0)

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

当 Pr 不为 1 时 — Chilton-Colburn#

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

\frac{C_f}{2} = St \cdot Pr^{2/3} \qquad (0.6 < Pr < 60)

这就是 Chilton-Colburn 比拟(修正雷诺比拟)。当 $Pr = 1$ 时它退回原式。请操作下面的图表。

log₁₀ Re_x: 5.0marker Pr: 0.70St(Pr=0.70) = 1.332e-3  |  C_f/2 = 1.050e-3
St · Pr^(2/3) = 1.050e-3 (= C_f/2)

橙色曲线(原始 $St$ )随 $Pr$ 上下移动。但青色曲线( $St \cdot Pr^{2/3}$ )无论 $Pr$ 取何值都紧贴虚线( $C_f/2$ )。修正指数 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$ 。

值得记住的要点#

速度边界层与温度边界层遵循同一形状的无量纲方程。唯一的差别是扩散系数， $Re$ 与 $Re \cdot Pr$ 。
当 $Pr = 1$ 且压力梯度为零时，速度场与温度场重合， $C_f/2 = St$ 。这就是雷诺比拟。
对一般流体，用 Chilton-Colburn 比拟 $C_f/2 = St \cdot Pr^{2/3}$ 从摩擦估算传热。