用同一个式子导出μ和k，Pr偏了1.98倍 — 决定123.5 K壁温的系数9/4

从同一平均自由程导出μ和k，Pr就被钉在γ上；一旦把转动能量单独计数，它就落到4γ/(9γ−5)。

用同一套导出给出μ和k，Pr必然等于γ#

把空气压到10个大气压，导热系数几乎不变。密度涨了10倍，热却以同样的速度流动。分子运动论用两行解释了这个反直觉的事实。可是从同样这两行里再取出黏性系数，普朗特数(Prandtl number，动量扩散与热扩散之比)就被钉在 $\gamma$ 上。对空气而言是1.400。实测值是0.707。偏大1.98倍的数写进CFD会把壁温算偏，1913年Eucken补上的系数 $9/4$ 负责把它拉回来。数密度滑块 number density 从1推到8，crossings 和 λ 反向移动，k̂ 几乎不动。

number density n × 4n=0 · λ=68px · crossings=0 · k̂=0.00

两条柱对冲的时候，k̂ 停在原地。只有 n = 1 时 $\lambda$ 超过壁间距，k̂ 才掉下去。

穿过统计面的分子带走两样东西

在流场中竖起一个虚拟平面。单位时间、单位面积穿过它的分子数为

\dot{N}'' = \tfrac{1}{4} n \bar{c}

$n$ 是数密度， $\bar{c}$ 是平均分子速率。向下穿过的分子，上一次碰撞发生在平面上方一个平均自由程(mean free path，分子在两次碰撞之间飞行的平均距离) $\lambda$ 处，所以它带着那个位置的值降下来。向上穿过的分子带着下方 $\lambda$ 处的值升上去。带走动量就成为黏性，带走动能就成为导热。

\mu = \tfrac{1}{3}\rho\,\bar{c}\,\lambda, \qquad k = \tfrac{1}{3}\rho\,\bar{c}\,\lambda\,c_v

其中 $\rho = nm$ ， $c_v$ 是单位质量定容比热。两个式子前半部分完全相同。只有末尾的 $c_v$ 不同。

密度提高8倍而导热系数不动的原因#

$\lambda$ 与密度成反比。对碰撞截面 $\sigma$ 有

\lambda = \frac{1}{\sqrt{2}\, n \sigma}

于是在 $\rho \bar{c} \lambda$ 里， $n$ 与 $1/n$ 互相抵消。

k = \frac{m c_v \bar{c}}{3\sqrt{2}\,\sigma}

密度消失了。分子数增多则穿面的分子变多，但每个分子出发的距离更短，带来的温差更小。 k̂ 是对真正穿过黄色平面的粒子所携带温度计数得到的值。n = 1 就是连续介质假设失效的地方。

把三条路线放进同一张表

现在把两个式子相除。普朗特数的定义是 $\Pr = \mu c_p / k$ 。

\Pr = \frac{\mu c_p}{k} = \frac{\tfrac{1}{3}\rho\bar{c}\lambda\, c_p}{\tfrac{1}{3}\rho\bar{c}\lambda\, c_v} = \frac{c_p}{c_v} = \gamma

$\rho$ 、 $\bar{c}$ 、 $\lambda$ 全部约掉。分子直径、温度、压力都不留在答案里。只剩下比热比。这不是导出过程的意外，而是它的结构：动量和能量被以同样的方式装在同一个载体上。黏性能用分子的一行写下来的地方，和黏性的分子起源是同一处。并不是某一个系数算错了，而是把两种不同的物理绑到了一个尺度上。

路线	Pr 公式	空气取值
平均自由程单项	$\gamma$	1.400
Eucken	$4\gamma/(9\gamma-5)$	0.737
实测	—	0.707

转动自由度不搬运动量

$\Pr = \gamma$ 失手的原因在于把 $c_v$ 当成一整块。双原子分子的 $c_v$ 分为平动3个自由度和转动2个自由度，而只有平动参与动量输运。一个旋转的哑铃即使穿过平面，它的转动能量对剪应力也毫无贡献。

Eucken把两个部分分开，各自给了不同的输运效率。平动能量与速度相关，快分子搬得更多，因此获得 $5/2$ 的权重。内能与速度不相关，权重为1。结果只有一行。

\frac{k}{\mu c_v} = 1 + \frac{9}{4}(\gamma - 1)

对 $\gamma \to 5/3$ 的单原子气体，这个值为 $1 + 9/4 \cdot 2/3 = 2.5$ ，正好是 $5/2$ ：没有转动，只剩平动权重。两边整理后

\Pr = \frac{\gamma}{1 + \tfrac{9}{4}(\gamma-1)} = \frac{4\gamma}{9\gamma - 5}

把三条路线全部套在六种气体上

脚本并排算出三条路线，再把错误的Pr一直带到真实计算里：带剪切加热的一维Couette层。一侧壁面以 $U$ 运动，另一侧壁面绝热。

\frac{d}{dy}\!\left(k\frac{dT}{dy}\right) + \mu\!\left(\frac{du}{dy}\right)^{2} = 0

$u$ 为线性时源项是常数，绝热壁的温升只与 $\Pr$ 成正比。

T(0) - T_w = \Pr\,\frac{U^{2}}{2 c_p}

import math
 
GASES = [
    ("He",  1.667, 0.68),
    ("Ar",  1.667, 0.67),
    ("N2",  1.404, 0.72),
    ("Air", 1.400, 0.707),
    ("CO2", 1.289, 0.77),
    ("H2O", 1.329, 0.98),
]
 
def pr_mean_free_path(gamma):
    """mu = (1/3) rho c_bar lam,  k = (1/3) rho c_bar lam c_v  ->  Pr = c_p/c_v."""
    return gamma
 
def pr_eucken(gamma):
    """k/(mu c_v) = 1 + (9/4)(gamma-1)  ->  Pr = 4 gamma / (9 gamma - 5)."""
    return 4.0 * gamma / (9.0 * gamma - 5.0)
 
print("gas   gamma   Pr_mfp  Pr_euck  Pr_meas   mfp/meas  euck/meas")
for name, g, meas in GASES:
    a, b = pr_mean_free_path(g), pr_eucken(g)
    print("%-5s %5.3f  %6.3f  %7.3f  %7.3f   %8.2f  %9.2f" % (name, g, a, b, meas, a / meas, b / meas))
 
def couette_wall_rise(Pr, U, cp, N=201):
    """d2T/dy2 = -Pr U^2/cp on y in [0,1]; adiabatic at y=0, T=0 at y=1. Thomas sweep."""
    h = 1.0 / (N - 1)
    src = Pr * U * U / cp
    lo = [0.0] * N; di = [0.0] * N; up = [0.0] * N; rhs = [0.0] * N
    di[0], up[0], rhs[0] = -1.0, 1.0, 0.0
    for i in range(1, N - 1):
        lo[i], di[i], up[i] = 1.0, -2.0, 1.0
        rhs[i] = -src * h * h
    lo[N - 1], di[N - 1], rhs[N - 1] = 0.0, 1.0, 0.0
    for i in range(1, N):
        w = lo[i] / di[i - 1]
        di[i] -= w * up[i - 1]
        rhs[i] -= w * rhs[i - 1]
    T = [0.0] * N
    T[N - 1] = rhs[N - 1] / di[N - 1]
    for i in range(N - 2, -1, -1):
        T[i] = (rhs[i] - up[i] * T[i + 1]) / di[i]
    return T
 
U, cp = 600.0, 1005.0
print("\nCouette shear layer, U=%.0f m/s, cp=%.0f J/kg/K, adiabatic wall at y=0" % (U, cp))
print("route            Pr    dT_wall[K]   Pr*U^2/(2cp)[K]")
rise = {}
for label, Pr in (("mean-free-path", pr_mean_free_path(1.4)),
                  ("Eucken", pr_eucken(1.4)),
                  ("measured", 0.707)):
    T = couette_wall_rise(Pr, U, cp)
    rise[label] = T[0]
    print("%-15s %5.3f  %10.1f   %14.1f" % (label, Pr, T[0], Pr * U * U / (2 * cp)))
 
print("\nwall-temperature gap = %.1f K   ratio = %.2f"
      % (rise["mean-free-path"] - rise["measured"], rise["mean-free-path"] / rise["measured"]))
 
print("\nk ~ sqrt(T) extrapolated from 300 K vs measured air conductivity")
k300 = 0.0263
for T_, k_meas in ((400, 0.0338), (600, 0.0469), (800, 0.0573), (1000, 0.0672)):
    k_law = k300 * math.sqrt(T_ / 300.0)
    print("T=%5d K   sqrtT=%.4f   measured=%.4f   error=%+6.1f%%"
          % (T_, k_law, k_meas, 100.0 * (k_law - k_meas) / k_meas))

gas   gamma   Pr_mfp  Pr_euck  Pr_meas   mfp/meas  euck/meas
He    1.667   1.667    0.667    0.680       2.45       0.98
Ar    1.667   1.667    0.667    0.670       2.49       0.99
N2    1.404   1.404    0.735    0.720       1.95       1.02
Air   1.400   1.400    0.737    0.707       1.98       1.04
CO2   1.289   1.289    0.781    0.770       1.67       1.01
H2O   1.329   1.329    0.764    0.980       1.36       0.78
 
Couette shear layer, U=600 m/s, cp=1005 J/kg/K, adiabatic wall at y=0
route            Pr    dT_wall[K]   Pr*U^2/(2cp)[K]
mean-free-path  1.400       249.5            250.7
Eucken          0.737       131.3            132.0
measured        0.707       126.0            126.6
 
wall-temperature gap = 123.5 K   ratio = 1.98
 
k ~ sqrt(T) extrapolated from 300 K vs measured air conductivity
T=  400 K   sqrtT=0.0304   measured=0.0338   error= -10.2%
T=  600 K   sqrtT=0.0372   measured=0.0469   error= -20.7%
T=  800 K   sqrtT=0.0429   measured=0.0573   error= -25.0%
T= 1000 K   sqrtT=0.0480   measured=0.0672   error= -28.5%

气体按钮和 Eucken 开关把表里的三列画成壁温曲线。关掉开关，只剩 Pr = $\gamma$ 的曲线挨着实测曲线。

γ 1.400wall speed 600 m/s

选 He 时，两条抛物线在壁端的间距张得最大。 $\gamma$ 越大，漏算转动自由度的代价越高。

决定123.5 K壁温的是系数9/4#

数值解与解析解在1 K以内吻合，这是201个网格点上的 $O(h^2)$ 误差。第三列是被拉开的温度。在 $U = 600$ m/s时，绝热壁温度用 $\gamma$ 当Pr算出249.5 K，用实测值算出126.0 K。差123.5 K。在预测超声速壁面热流的计算里，这个差距会改变选材结论。

Eucken列对表中所有单原子和双原子气体都落在实测值的2~4%以内。仅仅把一根哑铃拆成平动和转动，就把2.45倍的误差压到2%。

把空气的 $\gamma=1.400$ 代进 $4\gamma/(9\gamma-5)$ ，分子是 $4\times1.400=5.600$ ，分母是 $9\times1.400-5=7.600$ 。相除得 $0.7368$ ，和表里的Eucken值 $0.737$ 在小数第三位对齐。同一关系写成 $1+(9/4)(\gamma-1)$ ，则 $9/4\times0.400=0.900$ ， $1.400/1.900=0.737$ 。He的 $\gamma=1.667$ ，平均自由程路线是 $1.667$ ，Eucken是 $0.667$ ，实测是 $0.680$ 。 $1.667/0.680=2.45$ ，就是表里的 mfp/meas。Eucken除以实测是 $0.667/0.680=0.98$ 。

Couette里 $U=600$ 、 $c_p=1005$ ，于是 $U^2/(2c_p)=360000/2010=179.1$ K。乘 $1.400$ 得 $250.7$ K，乘 $0.707$ 得 $126.6$ K，乘 $0.737$ 得 $132.0$ K。201个网格点的数值解是 $249.5$ 、 $126.0$ 、 $131.3$ 。与解析解分别差 $1.2$ K、 $0.6$ K、 $0.7$ K，落在1 K以内。 $249.5-126.0=123.5$ ， $249.5/126.0=1.98$ 。把 $300$ K的 $k=0.0263$ 乘上 $\sqrt{1000/300}$ ， $\sqrt{3.333}\times0.0263$ 落在 $0.0480$ 附近。相对实测 $0.0672$ ， $(0.0480-0.0672)/0.0672$ 约是 $-28.6\%$ ，和输出的 $-28.5\%$ 只差一个舍入。H₂O是 $0.764/0.980=0.780$ ，即表里的 euck/meas = 0.78，大约低22%。

√T定律在1000 K处低28.5%的原因#

同一套导出也给出温度依赖。若 $\bar{c} \propto \sqrt{T}$ 且 $\sigma$ 为常数，则

k \propto \sqrt{T}

最后一段输出就是把这条定律从300 K外推的结果。400 K处低10.2%，1000 K处低28.5%。符号始终不翻转是线索： $\sigma$ 为常数的假设破了。快分子彼此擦过得更近，有效截面随之缩小。这正是 Chapman–Enskog理论用 $\Omega^{(2,2)*}$ 碰撞积分填上的位置，而工程代码则交给Sutherland 公式里的常数 $S$ 去吸收。在壁面上把摩擦和热捆成同一个比，正是雷诺比拟只靠阻力就把传热对上做的事。放进这个比的Pr偏了1.98倍，热流会跟着偏。

Eucken在H₂O上偏离22%的地方#

表中只有H₂O是 euck/meas = 0.78，低22%。水分子是极性分子，偶极相互作用支配其碰撞。转动松弛远慢于平动，所以"内能与速度不相关地被搬运"这个假设不成立。修正Eucken式给内部自由度单独乘上 $\rho D_{\text{int}}/\mu$ 来纳入这一效应。在水蒸气比例高的燃烧计算里，把Pr当常数就会把这个偏差原封不动地压到壁面热流上。

在代码里写下Pr = 0.71时真正签下的东西#

Pr = 0.71 这一行不是常数声明。它主张动量和能量由同样的分子搬运却以不同速率扩散，并且这个比值来自两个转动自由度。破坏这一主张的条件有三个：单原子气体(降到0.67)、极性多原子气体 (升到0.98)，以及 $\lambda$ 大于网格的稀薄气体(第一个模拟中的 n = 1)。如果计算落在这三者之一，改一个常数并不是活干完了。