1 μm のアルミ粒子は連続体の式より 4.6 倍遅く温まった — Knudsen 遷移領域の Sherman 調和補間

1 μm の粒子はなぜ計算より遅く溶けたのか#

1500 K の空気に常温のアルミニウム粒子を投げ込みます。直径 50 μm と 1 μm、先に溶けるのはどちらでしょうか。当然、小さいほうです。質量は $d^3$ で減り、熱伝導は $d$ でしか減らないからです。連続体の式では、1 μm 粒子は 1.95 μs で融点に達します。

ところが同じ粒子を気体分子運動論で計算すると 8.95 μs になります。4.6 倍の差です。粒子が周囲気体の平均自由行程(分子が次の衝突までに飛ぶ平均距離)と同程度の大きさになったためです。この記事では、アルミニウム凝集粒子の燃焼モデルの最初の段階である予熱(pre-heating)で、この差がどこから来るのかを扱います。連続体側と自由分子側の式をそれぞれ立て、Sherman の調和補間でつなぎます。最後に、結果が一行の式 $t/t_c = 1 + Q_c/Q_{fm}$ にまとまることを確かめます。

粒子ひとつから見た気体の粒々感

気体が連続体に見えるかどうかは、観測者の大きさで決まります。その比が Knudsen 数です。

$\lambda$ は平均自由行程、$d$ は粒子直径、$\sigma$ は分子衝突直径、$p$ は圧力です。1500 K・1 atm の空気で $\sigma = 3.7\times10^{-10}$ m とすると $\lambda = 0.336$ μm です。常温の約 5 倍の長さです。同じ圧力で温度が上がると分子数密度が下がるためです。

一般に $\mathrm{Kn} < 0.01$ を連続体、$\mathrm{Kn} > 10$ を自由分子領域と呼びます。その間が遷移領域(transition regime)です。50 μm 粒子は $\mathrm{Kn} = 0.013$、1 μm 粒子は 0.67 です。マイクロサイズのアルミニウム粉末のかなりの部分が遷移領域にかかります。

下のシミュレーションで、粒子直径と適応係数を実際に動かしてみましょう。

直径スライダーを下げると、左側の分子の直進距離が粒子より長くなります。粒子に触れて冷えた分子(青)は、連続体側では粒子のそばに留まり、自由分子側ではまっすぐ遠くへ飛び去ります。右下のレースでは、2 本のバーの到達時間比が赤い点の縦の値の逆数と一致することも確かめてください。

連続体側:熱は拡散だけで届く

$\mathrm{Kn} \to 0$ では、粒子まわりの気体は Fourier の法則に従います。静止気体中の球の定常熱伝導は Nusselt 数(Nu、対流熱伝達を伝導で割った無次元数)2 を与えます。

$k$ は気体の熱伝導率、$T_g$ と $T_p$ は気体と粒子の温度です。熱は直径にしか比例して入りません。質量は $d^3$ なので時間は $d^2$ です。直径を 10 分の 1 にすると 100 倍速く温まります。

表の $t_c$ がその二乗をそのまま出します。50 μm は $4.883\times10^{-3}$ s、1 μm は $1.953\times10^{-6}$ s です。直径比 50 の二乗は 2500 で、時間比も 2500 です。10 μm と 1 μm、3 μm と 0.3 μm は直径比 10 に対して時間比 100 です。$1.953\times10^{-4}/1.953\times10^{-6}=100$。$1.758\times10^{-5}/1.758\times10^{-7}=100$。連続体側では直径だけで加熱時間の比が閉じます。液滴が乾く時間も同じ $d^2$ です。その釣り合いは拡散が支配する蒸発でもう一度出てきます。

この式は、粒子表面の気体温度が $T_p$ に等しいと仮定しています。実際には、表面から $\lambda$ 程度内側の気体は十分な衝突を経験していません。その層の厚さが粒子に比べて無視できるときだけ、この仮定が成り立ちます。

自由分子側:当たる分子を数える

$\mathrm{Kn} \to \infty$ では、気体分子は粒子に届くまで互いに衝突しません。熱伝達は、単位時間に表面へ当たる分子数と、分子 1 個が渡すエネルギーの積です。

$\alpha$ は熱適応係数(thermal accommodation coefficient、衝突した分子がどれだけ壁温度に揃って離れるかの割合)、$\bar c$ は平均分子速さ、$\gamma$ は比熱比、$m$ は分子質量です。今度は $Q_{fm} \propto d^2$ です。熱は表面積を通って入ります。加熱時間は $d$ に比例するだけです。

二つの式の比をとると構造が見えます。

$k \sim \rho c_v \bar c \lambda$ を代入すると、この比はおおよそ $\alpha/\mathrm{Kn}$ に比例します。粒子が小さいほど、自由分子の式は連続体の式より小さくなります。小さいほうがボトルネックです。

Sherman の調和補間:二つの熱抵抗を直列につなぐ#

遷移領域では二つの機構が順番に働きます。粒子表面から $\lambda$ 程度までは分子が飛んで熱を運び、その外側は拡散が運びます。二段階を直列の熱抵抗とみなすと、Sherman(1963)の調和形が得られます。

両端で厳密です。$\mathrm{Kn} \to 0$ では $1/Q_{fm}$ が消え、$\mathrm{Kn} \to \infty$ では $1/Q_c$ が消えます。元の発表資料は、この形のほかに Fuchs の境界球(boundary sphere)アプローチと、Fuchs–Sutugin、Loyalka の調和形もまとめています。境界球アプローチは、粒子を半径 $\sim\lambda$ の仮想球で包み、内側を自由分子、外側を連続体として解き、境界で二つの解を接続します。調和形はその結果を係数ひとつに圧縮した近似です。

コード:直径 6 通りと適応係数 4 通り#

下のコードは上の三つの式をそのまま移したものです。粒子の熱収支 $m c_p\, dT/dt = Q(T)$ を 300 K からアルミニウムの融点 933 K まで積分します。時間ではなく温度で積分すると $dt = m c_p\, dT / Q$ となり、ループひとつで済みます。放射と表面反応はあえて外しました。

import math
 
KB = 1.380649e-23          # Boltzmann constant [J/K]
AMU = 1.66054e-27          # atomic mass unit [kg]
 
# gas: air at 1500 K, 1 atm
TG, P = 1500.0, 101325.0
K_GAS = 0.095              # thermal conductivity [W/m/K]
GAMMA = 1.4
M_GAS = 28.97 * AMU
SIGMA = 3.7e-10            # collision diameter [m]
 
# particle: aluminum
RHO_P, CP_P = 2700.0, 1100.0
T0, T_MELT = 300.0, 933.0
 
def mean_free_path(T, p):
    return KB * T / (math.sqrt(2.0) * math.pi * SIGMA**2 * p)
 
def q_continuum(d, Tp):
    return 2.0 * math.pi * d * K_GAS * (TG - Tp)
 
def q_free_molecular(d, Tp, alpha):
    cbar = math.sqrt(8.0 * KB * TG / (math.pi * M_GAS))
    return alpha * math.pi * d**2 * P * cbar / 8.0 \
        * (GAMMA + 1.0) / (GAMMA - 1.0) * (TG - Tp) / TG
 
def q_sherman(d, Tp, alpha):
    qc, qf = q_continuum(d, Tp), q_free_molecular(d, Tp, alpha)
    return 1.0 / (1.0 / qc + 1.0 / qf)
 
def heat_up_time(d, law, n=4000):
    # integrate m*cp*dT/dt = Q(T) from T0 to T_MELT in T (dt = m cp dT / Q)
    m = RHO_P * math.pi * d**3 / 6.0
    dT = (T_MELT - T0) / n
    t = 0.0
    for i in range(n):
        T = T0 + (i + 0.5) * dT
        t += m * CP_P * dT / law(d, T)
    return t
 
lam = mean_free_path(TG, P)
print(f"lambda at {TG:.0f} K, 1 atm = {lam*1e6:.3f} um")
print()
print("alpha = 1.0")
print(f"{'d [um]':>8} {'Kn=2l/d':>8} {'Qfm/Qc':>8} {'Qsh/Qc':>8} {'t_c [s]':>10} {'t_sh [s]':>10} {'ratio':>6}")
for d_um in [50, 10, 3, 1, 0.3, 0.1]:
    d = d_um * 1e-6
    kn = 2 * lam / d
    r_fm = q_free_molecular(d, T0, 1.0) / q_continuum(d, T0)
    r_sh = q_sherman(d, T0, 1.0) / q_continuum(d, T0)
    tc = heat_up_time(d, q_continuum)
    ts = heat_up_time(d, lambda dd, T: q_sherman(dd, T, 1.0))
    print(f"{d_um:8.1f} {kn:8.3f} {r_fm:8.2f} {r_sh:8.3f} {tc:10.3e} {ts:10.3e} {ts/tc:6.2f}")
print()
print("d = 1 um, varying accommodation coefficient")
print(f"{'alpha':>6} {'Qsh/Qc':>8} {'t_sh/t_c':>9}")
d = 1e-6
tc = heat_up_time(d, q_continuum)
for a in [1.0, 0.5, 0.1, 0.02]:
    ts = heat_up_time(d, lambda dd, T: q_sherman(dd, T, a))
    print(f"{a:6.2f} {q_sherman(d, T0, a)/q_continuum(d, T0):8.3f} {ts/tc:9.2f}")

lambda at 1500 K, 1 atm = 0.336 um
 
alpha = 1.0
  d [um]  Kn=2l/d   Qfm/Qc   Qsh/Qc    t_c [s]   t_sh [s]  ratio
    50.0    0.013    13.96    0.933  4.883e-03  5.233e-03   1.07
    10.0    0.067     2.79    0.736  1.953e-04  2.653e-04   1.36
     3.0    0.224     0.84    0.456  1.758e-05  3.857e-05   2.19
     1.0    0.672     0.28    0.218  1.953e-06  8.949e-06   4.58
     0.3    2.240     0.08    0.077  1.758e-07  2.275e-06  12.94
     0.1    6.721     0.03    0.027  1.953e-08  7.191e-07  36.82
 
d = 1 um, varying accommodation coefficient
 alpha   Qsh/Qc  t_sh/t_c
  1.00    0.218      4.58
  0.50    0.122      8.16
  0.10    0.027     36.82
  0.02    0.006    180.09

4.58 は一行の式から出てくる#

どちらの式も $(T_g - T_p)$ に比例するので、比 $Q_{fm}/Q_c$ は粒子温度に依存しません。すると Sherman の加熱率は、連続体の加熱率に定数 $s = 1/(1 + Q_c/Q_{fm})$ を掛けたものになります。加熱時間の比はその逆数です。

1 μm では $Q_{fm}/Q_c = 0.28$ なので $1 + 1/0.28 \approx 4.6$ です。表の 4.58 と一致します。0.1 μm では $1 + 1/0.03$ で約 37 倍です。積分ループはこの算数を確かめているだけです。

調和側も同じ割り算です。$Q/Q_c=r/(1+r)$、$r=Q_{fm}/Q_c$ です。1 μm は $0.28/1.28\approx0.219$ で、表の 0.218 です。10 μm は $2.79/3.79=0.736$ で表と一致します。3 μm は $0.84/1.84\approx0.457$、表は 0.456 です。50 μm は $13.96/14.96=0.933$ です。出力は丸め 1 桁の中でこの式に付きます。

$Q_{fm}=Q_c$ なら時間比はちょうど 2 です。熱は半分、時間は 2 倍です。本文の 3.6 μm、$\mathrm{Kn}\approx0.19$ がその点です。連続体と呼ぶには早い直径で、すでに 2 倍です。

$\alpha$ も $r$ にそのまま掛かります。$\alpha=0.5$ なら $r=0.14$、$1+1/0.14\approx8.14$ です。表の 8.16 です。$\alpha=0.1$ なら $r\approx0.028$、$1+1/0.028\approx36.7$ です。表の 36.82 です。$\alpha=0.02$ は $1+1/0.0056\approx180$、表の 180.09 です。$r/(1+r)$ も 0.123、0.027、0.006 で、表の 0.122、0.027、0.006 に付きます。新しい計算ではありません。出力ブロックを割り直しただけです。

より注目すべきは 10 μm の行です。$\mathrm{Kn} = 0.067$ と、ふつう連続体と呼ぶ範囲なのに、伝わる熱量が 26% 減っています。$Q_{fm} = Q_c$ となる交差直径は約 3.6 μm、$\mathrm{Kn} \approx 0.19$ です。交差点は $\mathrm{Kn} = 1$ より 5 倍小さいところにあります。$(\gamma+1)/(\gamma-1) = 6$ と $1/16$、そして $k$ と $\lambda$ を結ぶ係数が絡んだ結果です。「Kn が 0.1 より小さいから連続体」という感覚は、熱伝達では緩すぎます。

適応係数 α は直径と同じくらい効く

表の下半分は、1 μm で $\alpha$ だけを変えた結果です。$\alpha = 0.1$ で加熱時間は 37 倍、0.02 で 180 倍に伸びます。$Q_{fm}$ には $\alpha d$ の積としてしか入らないからです。$\alpha$ を 10 分の 1 にすることは、直径を 10 分の 1 にすることと同じです。

問題は、$\alpha$ がよくわかっていないことです。金属表面と高温気体の間の値は、表面の酸化膜、温度比、気体の種類によって報告値の幅が広いです。Filippov & Rosner(2000)は温度比が大きい条件でのエネルギー伝達そのものを見直しました。ナノアルミニウムの着火遅れを合わせるために $\alpha$ を 1 よりかなり小さく取った研究もあります。遷移領域の粒子モデルでは、$\alpha$ は事実上の調整パラメータになりがちです。その値が結果を何倍も変えるという事実は、報告書に書いておくべきです。

調和補間が揺らぐ区間

Sherman の形は両端を合わせますが、中間の形は近似です。$\mathrm{Kn} \sim 0.1$–1 では、Fuchs の境界球解、Fuchs–Sutugin の式、Loyalka の補正が互いに数十 % の範囲で食い違います。上の 10 μm の 26% も、このモデル幅の中にある数字として読むべきです。

$\lambda$ の定義も同程度の不確かさを生みます。剛体球の式か粘性に基づく定義か、Kn を $\lambda/r$ と書くか $2\lambda/d$ と書くかで、同じ粒子の Kn は簡単に数十 % 変わります。他の論文から相関式を借りるときは、その論文の $\lambda$ の定義も一緒に持ってくる必要があります。そうしないと境界球の半径がおかしな場所に置かれます。

燃焼に進んでもボトルネックは残ります。粒子のそばの酸化剤拡散も、物質移動側の Knudsen 補正を受けます。凝集体になると有効直径が大きくなり、Kn は再び小さくなります。そのとき Nu=2 に戻してよいかは、その有効直径で Kn を測り直さないと分かりません。

Lagrangian 粒子に熱伝達の一行を書く前に#

Euler–Lagrange コードの粒子熱伝達は、たいてい $\mathrm{Nu} = 2 + 0.6\,\mathrm{Re}^{1/2}\mathrm{Pr}^{1/3}$ の一行から始まります。この式は連続体側です。裾が 3 μm なら表の比は 2.19、1 μm なら 4.58 です。2〜3 倍という感覚は 3 μm 付近だけで、1 μm ではすでに 4 倍を超えます。着火遅れが小さい直径で妙に短いなら、この一行です。

直す演算は軽いです。セルの $\lambda$、粒子の $Q_{fm}$、その調和平均です。代わりに $\alpha$ と $\lambda$ の定義を書かないと 4.58 は再現しません。その倍率はほぼその二つの数字です。熱流束を上げても傾きが保たれるわけでもありません。沸騰曲線が折れ曲がる限界熱流束と同じで、機構が変わると伝達量が折れます。ここでの折れは $Q_{fm}$ です。