逆反応定数を Kp で割ったら 2500 K で AlO が 99.6% 解離した

逆反応定数は kf/Kc で作ります。Kc = Kp(p°/RT)^Δn なので、Δn=0 なら Kp でも正しく、モル数が変わる反応だけが単位系次第で 4 倍から 25 万倍ずれます。

Chemkin と自作コードが 134 K ずれた#

アルミニウム粒子の燃焼モデルを自作した研究室の発表資料があります。気相反応 12 本を Arrhenius 式で組み込み、同じ機構を Chemkin でも計算しました。最高温度は一方が 3800 K、もう一方が 3666 K でした。反応速度を一つずつ比べると、R2〜R5 の逆反応定数だけが間違っていました。次のスライドの題は「平衡定数」でした。

この記事はその次のスライドを埋めます。逆反応定数を作る式には単位が隠れています。それを落とすと、正しいままの反応と間違う反応に分かれます。AlO 解離反応器ひとつでこの誤りを再現し、どの反応が危ないかを見分ける規則を得ます。

逆反応定数は別に与えられない

反応機構ファイルは通常、正反応の係数しか持ちません。 $k_f = A T^b \exp(-E_a/RT)$ の $A$ 、 $b$ 、 $E_a$ の三つです。逆反応は詳細釣り合い(detailed balance、平衡で正反応と逆反応の速度が等しいという条件)から作ります。

k_r = \frac{k_f}{K_c}, \qquad K_c = \prod_i [X_i]_{eq}^{\nu_i}

$k_r$ は逆反応定数、 $K_c$ は濃度基準の平衡定数、 $\nu_i$ は生成物で正・反応物で負の化学量論係数です。

ところが熱力学が直接くれるのは $K_c$ ではありません。Gibbs 自由エネルギーから計算されるのは圧力基準の $K_p$ です。

K_p = \exp\left(\frac{\Delta S^\circ}{R} - \frac{\Delta H^\circ}{RT}\right) = \prod_i \left(\frac{p_i}{p^\circ}\right)^{\nu_i}

$p^\circ$ は標準圧力(1 atm)、 $\Delta H^\circ$ ・ $\Delta S^\circ$ は反応エンタルピーとエントロピーです。 $K_p$ は無次元で、 $K_c$ は(濃度) $^{\Delta n}$ の単位を持ちます。

Kp と Kc のあいだに (p°/RT)^Δn がある#

理想気体では $p_i = [X_i]RT$ です。これを $K_p$ に代入すると二つの定数の関係が出ます。

K_c = K_p \left(\frac{p^\circ}{RT}\right)^{\Delta n}, \qquad \Delta n = \sum_i \nu_i

$\Delta n$ は反応で増える気体のモル数です。 $p^\circ/RT$ は標準状態の気体のモル濃度です。

この因子は反応ごとに違う働きをします。 $\mathrm{Al + O_2 \rightleftharpoons AlO + O}$ は両辺とも 2 モルです。 $\Delta n = 0$ なので $K_c = K_p$ となり、 $K_p$ をそのまま使っても答えは合います。 $\mathrm{AlO + M \rightleftharpoons Al + O + M}$ は一つが二つになり、 $\Delta n = +1$ です。3000 K では $p^\circ/RT = 4.06\ \mathrm{mol/m^3}$ なので、SI のコードは $K_c$ を 4 倍小さく見積もります。濃度を $\mathrm{mol/cm^3}$ で扱う Chemkin 流のコードなら、同じ値は $4.06\times10^{-6}$ です。今度は $K_c$ が 25 万倍大きくなります。

下の箱で実際に操作してみましょう。

T = 3000 K p = 1.00 atm

box alpha = 0.000 — with dn = 0 all three buttons give the same box; with dn = +1 the CGS reading floods it with free atoms

$\Delta n = 0$ ボタンを押すと、三つのモードの箱は同じになります。 $\Delta n = +1$ で CGS モードを選ぶと、結合線がほぼすべて切れます。圧力スライダーを動かすと真の平衡は変わりますが、 $K_p$ の目盛りは動きません。

R2〜R5 だけが間違っていた理由#

発表資料は原因を書いていません。しかし「一部の反応の逆反応だけが違う」という症状は、この因子とぴったり合います。すべての反応に同じ誤りをしても、 $\Delta n = 0$ の反応は影響を受けません。間違うのは解離・再結合のようにモル数が変わる反応だけです。

アルミニウムの気相機構にはそうした反応が混ざっています。AlO + M の解離や $\mathrm{Al_2O}$ の生成などの三体反応(third body、衝突エネルギーだけをやり取りする第三の分子)がそうです。スライドに「Third body=?」というメモがあるのも偶然ではないかもしれません。

二つ目の落とし穴は向きです。SI で $K_p$ を使うと $K_c$ が小さすぎ、解離が足りません。CGS で使うと $K_c$ が大きすぎ、解離しすぎます。同じ誤りが単位系によって逆向きに外れます。最高温度が高く出るか低く出るかを前もって言えないのもそのためです。

コード: 1% の AlO を 3000 K の反応器に入れた#

等温・等積の反応器に、アルゴンで希釈した AlO を 1% 入れます。解離反応ひとつだけを解き、逆反応定数を三通りに作ります。熱力学値は 298 K の近似値を定数として固定しました。正反応の前指数 $A$ は例示値です。目的は正確な燃焼計算ではなく、単位の誤りの大きさを見ることです。

時間進行は後退 Euler(backward Euler、次の時刻の値で右辺を評価する陰解法)です。式が $\xi$ の 2 次式なので、毎ステップ解の公式で厳密に解けます。逆反応が 25 万倍速くなっても発散しません。陰的時間進行が硬い項を扱う方法は semi-implicit Gear スキームの記事でも扱いました。

import math
 
R = 8.314462618        # J/(mol K)
P0 = 101325.0          # 標準圧力, Pa
 
# 298 K の反応エンタルピー・エントロピーを定数で固定 (JANAF 近似値)
DISSOC = dict(name="AlO + M <=> Al + O + M", DH=487.7e3, DS=107.2, dn=1)
SHUFFLE = dict(name="Al + O2 <=> AlO + O", DH=10.7e3, DS=9.8, dn=0)
MODES = ["Kc (correct)", "Kp as Kc, SI", "Kp as Kc, CGS"]
 
def kp_of(rx, T):
    return math.exp(rx["DS"] / R - rx["DH"] / (R * T))   # 無次元
 
def kc_true(T, rx=DISSOC):
    return kp_of(rx, T) * (P0 / (R * T)) ** rx["dn"]      # (mol/m^3)^dn
 
def kc_used(T, mode, rx=DISSOC):
    if mode == "Kc (correct)":
        return kc_true(T, rx)
    if mode == "Kp as Kc, SI":
        return kp_of(rx, T)                              # Kp を mol/m^3 として読む
    if mode == "Kp as Kc, CGS":
        return kp_of(rx, T) * 1e6 ** rx["dn"]            # Kp を mol/cm^3 として読む
 
def k_forward(T):
    return 1.0e11 * math.exp(-DISSOC["DH"] / (R * T))    # m^3/(mol s), 例示値
 
def run_reactor(T, p, x_alo, mode, t_end=0.2, nsteps=500):
    M = p / (R * T)                                  # Ar 希釈: [M] ~ 一定
    c0 = x_alo * M
    a = k_forward(T) * M                             # 1/s
    b = a / kc_used(T, mode)                         # kr*[M], m^3/(mol s)
    dt = t_end / nsteps
    xi = 0.0
    for _ in range(nsteps):                          # 後退 Euler, 2次方程式の厳密解
        A, B, C = b * dt, 1 + a * dt, -(xi + a * dt * c0)
        xi = (-B + math.sqrt(B * B - 4 * A * C)) / (2 * A)
    return xi / c0, xi * xi / (c0 - xi) / kc_true(T)
 
def alpha_exact(T, p, x_alo):
    M = p / (R * T); c0 = x_alo * M; K = kc_true(T)
    xi = (-K + math.sqrt(K * K + 4 * K * c0)) / 2
    return xi / c0
 
T, p, x = 3000.0, P0, 0.01
print(f"T={T:.0f} K  p=1 atm  x_AlO={x}")
print(f"Kp={kp_of(DISSOC, T):.3e}   P0/RT={P0/(R*T):.3f} mol/m^3   Kc={kc_true(T):.3e} mol/m^3")
print(f"exact equilibrium alpha = {alpha_exact(T, p, x):.4f}")
for mode in MODES:
    al, q = run_reactor(T, p, x, mode)
    print(f"{mode:15s}  alpha={al:.4f}   Q/Kc_true={q:10.3e}")
 
print("\nT scan (alpha after 200 ms)")
print("   T  correct  SI-bug  CGS-bug  exact")
for T in [2500, 3000, 3500, 4000]:
    r = [run_reactor(T, p, x, m)[0] for m in MODES]
    print(f"{T:5d}  {r[0]:.4f}  {r[1]:.4f}  {r[2]:.4f}  {alpha_exact(T, p, x):.4f}")
 
print("\np scan at 3000 K (alpha)")
for pa in [0.1, 1, 10]:
    r = [run_reactor(3000, pa*P0, x, m)[0] for m in MODES]
    print(f"{pa:5.1f} atm  {r[0]:.4f}  {r[1]:.4f}  {r[2]:.4f}")
 
print("\ndn=0 check:", SHUFFLE["name"])
print("  Kc used / Kc true:", [round(kc_used(3000, m, SHUFFLE) / kc_true(3000, SHUFFLE), 6) for m in MODES])

T=3000 K  p=1 atm  x_AlO=0.01
Kp=1.282e-03   P0/RT=4.062 mol/m^3   Kc=5.209e-03 mol/m^3
exact equilibrium alpha = 0.2997
Kc (correct)     alpha=0.2997   Q/Kc_true= 1.000e+00
Kp as Kc, SI     alpha=0.1626   Q/Kc_true= 2.462e-01
Kp as Kc, CGS    alpha=1.0000   Q/Kc_true= 2.462e+05
 
T scan (alpha after 200 ms)
   T  correct  SI-bug  CGS-bug  exact
 2500  0.0494  0.0227  0.9962  0.0494
 3000  0.2997  0.1626  1.0000  0.2997
 3500  0.7391  0.5311  1.0000  0.7391
 4000  0.9473  0.8658  1.0000  0.9473
 
p scan at 3000 K (alpha)
  0.1 atm  0.6602  0.4258  1.0000
  1.0 atm  0.2997  0.1626  1.0000
 10.0 atm  0.1070  0.0546  0.9997
 
dn=0 check: Al + O2 <=> AlO + O
  Kc used / Kc true: [1.0, 1.0, 1.0]

上のコードをブラウザで直接動かしてみましょう。

T = 3000 K p = 1.00 atm x_AlO = 0.0100 steps/frame = 2

step 0/500 (t = 0.0 ms, backward Euler dt = 0.4 ms) | Kp = 1.282e-3 | Kc = 5.209e-3 mol/m^3
Kc (correct)   : alpha = 0.0000 | Q/Kc_true = 0.000e+0
Kp as Kc, SI   : alpha = 0.0000 | Q/Kc_true = 0.000e+0
Kp as Kc, CGS  : alpha = 0.0000 | Q/Kc_true = 0.000e+0

既定値では三本の曲線が出力ブロックの 0.2997、0.1626、1.0000 に向かいます。T スライダーを 2500 K まで下げると、緑の曲線は 0.05 付近にとどまり、赤い曲線だけが 1 まで跳ね上がります。

0.2997、0.1626、1.0000 が語ること#

正しいモードは厳密な平衡解 0.2997 と 4 桁まで一致します。 $Q/K_c$ (反応商を真の平衡定数で割った値)も 1.000 です。反応器が本当の平衡に達したということです。

SI のバグは $Q/K_c = 0.2462$ で止まります。この数は $1/4.062$ 、つまり $RT/p^\circ$ です。反応器は間違った平衡に正確に収束しました。残差はゼロで、警告も出ません。収束判定ではこの誤りを捕まえられないということです。

CGS のバグの $2.462\times10^5$ は $10^6/4.062$ です。単位換算の因子 $10^6$ がそのまま平衡を押しています。2500 K は本来 AlO が 5% しか解離しない温度です。それが 99.6% 解離しました。解離は吸熱なので、実際の燃焼計算なら気体は大きく冷えます。発表資料の 134 K よりはるかに大きな誤差もありえます。

圧力を変えるとバグの見え方が変わる

$\Delta n = +1$ の反応の真の平衡は圧力に敏感です。0.1 atm で 0.66、10 atm で 0.11 です。ルシャトリエの原理どおり、低圧ではモル数が増える側へ進みます。

SI のバグはこの傾向に従いつつ値だけが低くなります。どの圧力でも $K_c$ が $RT/p^\circ$ 倍小さいからです。傾向が合っているので、グラフだけ見るともっともらしく見えます。実験と比べるとき「モデル定数の調整」で覆い隠されやすい誤差です。

CGS のバグでは圧力依存性がほとんど消えます。10 atm でも 0.9997 です。高圧燃焼室の解析でこのような無感度が見えたら、まず単位を疑う価値があります。アルミニウム粒子が受ける熱の流れも圧力や粒径で分かれます。その話は Knudsen 遷移領域の記事にあります。

係数表を移すときに単位が隠れる場所

$K_c$ の換算だけではありません。Chemkin 形式の $A$ は $\mathrm{cm^3/(mol\,s)}$ 単位です。2 分子反応なら SI に移すとき $10^{-6}$ を掛けます。三体反応は $[M]$ がもう一つ付くので $10^{-12}$ です。次数ごとに因子が変わります。

そのため単位の誤りは、いつも反応次数や $\Delta n$ と一緒に現れます。すべての反応が一緒に間違っていれば、むしろ見つけやすいのです。一部だけが間違うと、機構そのものを疑ってしまいます。発表資料が「逆反応定数と Kp の関係を確認し直す必要」と書いたのは正しい方向でした。

逆反応が一つ怪しいときに回す 2 行#

平衡計算ツールが一つあれば十分です。同じ T・p で CEA や Cantera から平衡組成を求めます。自作コードを等温反応器として長く回し、反応ごとに $Q/K_c$ を出力します。

すべての反応が 1 なら合格です。ある反応が $RT/p^\circ$ やその逆数で止まるなら、 $K_p$ をそのまま使っています。 $10^5$ 〜 $10^6$ 付近なら CGS 換算まで重なっています。数字そのものが原因を指しています。12 本の反応を一つずつ眺める前に、まずこの表を一枚作ってみてください。

逆反応定数を Kp で割ったら 2500 K で AlO が 99.6% 解離した — 平衡定数に潜む (p°/RT)^Δn