μとkを同じ式から導いたらPrが1.98倍ずれた — 壁温123.5 Kを分けた係数9/4

μとkを同じ導出で出すとPrは必ずγになる#

空気を10気圧まで圧縮しても熱伝導率はほとんど変わりません。密度が10倍になったのに、熱は同じ 速さで流れます。分子運動論はこの反直観を2行で説明します。ところが同じ2行から粘性係数まで 取り出すと、プラントル数(Prandtl number、運動量拡散と熱拡散の比)が $\gamma$ に釘付けに なります。空気なら1.400です。実測値は0.707です。1.98倍ずれた数をCFDコードに入れると壁温がずれ、1913年にEuckenが埋めた係数 $9/4$ がそのずれを引き戻します。 数密度スライダー number density を1から8へ上げると、crossings と λ は逆向きに動き、k̂ はほとんど動きません。

二つの棒が逆に動くあいだ、k̂ はほぼその場に残ります。n = 1 のときだけ $\lambda$ が壁間隔を超えて k̂ が落ちます。

集計面を横切る分子が運ぶ二つのもの

流れの中に仮想的な平面を1枚立てます。この面を単位時間・単位面積あたり横切る分子数は

$n$ は数密度、$\bar{c}$ は平均分子速さです。下向きに横切る分子は、平均自由行程 (mean free path、衝突と衝突の間に分子が飛ぶ平均距離)$\lambda$ だけ上で最後に衝突しました。 ですからその位置の値を持って降りてきます。上向きに横切る分子は $\lambda$ だけ下の値を持って 上がります。運動量を運べば粘性になり、運動エネルギーを運べば熱伝導になります。

$\rho = nm$、$c_v$ は単位質量あたりの定積比熱です。二つの式は前半が完全に同じです。後ろに 付いた $c_v$ だけが違います。

密度を8倍にしても熱伝導率が変わらない理由#

$\lambda$ は密度に反比例します。衝突断面積 $\sigma$ に対して

したがって $\rho \bar{c} \lambda$ の中で $n$ と $1/n$ が互いを消します。

密度が消えました。分子数が増えれば面を横切る分子は多くなりますが、各分子はより短い距離から 出発するので運んでくる温度差が小さくなります。k̂ は、実際に黄色い平面を横切った粒子が運んだ 温度を数えて作った値です。n = 1 は連続体仮定が破れる場所です。

三つの経路を一つの表に並べる

二つの式を割ります。プラントル数の定義は $\Pr = \mu c_p / k$ です。

$\rho$、$\bar{c}$、$\lambda$ がすべて約分されます。分子直径も温度も圧力も答えに残りません。 残るのは比熱比だけです。これは導出の事故ではなく導出の構造です。運動量とエネルギーを同じ 運び手に同じやり方で載せたからです。粘性が分子の一行で書ける場所は、粘性の分子的な起源と同じです。どこか一つの係数を間違えたのではなく、二つの異なる物理を一つの スケールに束ねたのです。

回転自由度は運動量を運ばない

$\Pr = \gamma$ が外れる理由は $c_v$ を一塊で書いたことにあります。二原子分子の $c_v$ は 並進3自由度と回転2自由度に分かれます。運動量輸送に関与するのは並進だけです。回転する ダンベルが平面を横切っても、その回転エネルギーはせん断応力に寄与しません。

Euckenは二つの成分を分離し、それぞれに違う輸送効率を与えました。並進エネルギーは速度と 相関しているので、速い分子がより多く運び、$5/2$ 倍の重みを受けます。内部エネルギーは速度と 無相関なので重みは1です。結果は1行です。

$\gamma \to 5/3$ の単原子気体ではこの値が $1 + 9/4 \cdot 2/3 = 2.5$、つまりちょうど $5/2$ に なります。回転がないので並進の重みだけが残ります。両辺を整理すると

気体6種に三つの経路をすべて掛けてみた#

三つの経路を並べて計算し、間違ったPrが実際の計算で何を変えるかまでつなげました。せん断加熱 のある1次元Couette層です。片方の壁が $U$ で動き、反対側の壁は断熱です。

$u$ が線形なら源項は定数になり、断熱壁の温度上昇は $\Pr$ にのみ比例します。

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$ の曲線だけが実測曲線の横に残ります。

Heを選ぶと、二つの放物線の壁端の間隔が最も大きく開きます。$\gamma$ が大きいほど、回転自由度を落とした代償は大きくなります。

壁温123.5 Kを分けたのは係数9/4だった#

数値解と解析解が1 K以内で一致します。格子201点での $O(h^2)$ 誤差です。3列目が分けた温度です。 $U = 600$ m/sで断熱壁の温度は、Prを $\gamma$ で取ると249.5 K、実測値で取ると126.0 Kです。 123.5 Kの差です。超音速の壁熱流束を予測する計算なら、この差で材料選定が変わります。

Eucken列はすべての単原子・二原子気体で実測値の2〜4%以内に入ります。ダンベル1本を並進と回転 に分けただけで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$ と小数第3位で合います。同じ関係を $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$ が定数なら

最後の出力ブロックがこの法則を300 Kから外挿した結果です。400 Kで10.2%、1000 Kで28.5%低く 出ます。符号が片側だけに寄るのがヒントです。$\sigma$ を定数と見た仮定が破れたのです。速い 分子は互いをより近くかすめるので有効断面積が縮みます。Chapman–Enskog理論が $\Omega^{(2,2)*}$ 衝突積分で埋める場所であり、実務ではSutherland式の定数 $S$ が代わりに 吸収します。壁で摩擦と熱を一つの比に束ねる仕事は、レイノルズの相似が抗力だけで熱伝達を当てる話と同じです。その比に入れるPrが1.98倍ずれれば、熱流束も一緒にずれます。

H₂OでEuckenが22%外れる場所#

表でH₂Oだけが euck/meas = 0.78 です。22%低い。水分子は極性で、双極子相互作用が衝突を 支配します。回転緩和が並進よりはるかに遅いため、内部エネルギーが「速度と無相関に運ばれる」 という仮定が成り立ちません。修正Eucken式は内部自由度に $\rho D_{\text{int}}/\mu$ を別に 掛けてこの効果を入れます。水蒸気比率の高い燃焼計算でPrを定数に置くと、この差がそのまま壁 熱流束に乗ります。

コードにPr = 0.71と書くとき実際に署名するもの#

Pr = 0.71 の1行は定数宣言ではありません。同じ分子が両方を運ぶのに運動量とエネルギーが 違う速さで拡散するという主張であり、その比が回転自由度2個から来るという主張です。その主張が 破れる条件は3つです。単原子(0.67まで下がる)、極性多原子(0.98まで上がる)、そして $\lambda$ が格子より大きくなる希薄気体(最初のシミュレーションの n = 1)です。この3つの どれかに入る計算なら、定数を一つ変えて終わりにはなりません。