用 Kp 去除逆反应速率常数，2500 K 下 AlO 竟解离了 99.6%

逆反应速率常数由 kf/Kc 得到。由于 Kc = Kp(p°/RT)^Δn，Δn=0 时用 Kp 也对；只有摩尔数变化的反应，会因单位制不同错 4 倍到 25 万倍。

Chemkin 和自编代码差了 134 K#

有一份实验室报告，讲的是自己编写的铝颗粒燃烧模型。十二个气相反应用 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 mol。 $\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_p$ ， $K_c$ 偏大，解离过度。同一个错误，随单位制不同而朝相反方向偏。这也是为什么无法预先判断最高温度会偏高还是偏低。

代码：把 1% 的 AlO 放进 3000 K 反应器#

在等温、等容反应器中放入用氩气稀释的 1% AlO。只求解解离反应，并用三种方式构造逆反应速率常数。热力学数据取 298 K 近似值并固定为常数。正反应指前因子 $A$ 是示例值。目的不是精确的燃烧计算，而是看清单位错误的量级。

时间推进采用后向 Euler(backward Euler，用新时刻的值计算右端项的隐式方法)。方程是 $\xi$ 的二次式，所以每一步都可以用求根公式精确求解。即使逆反应快了 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, 二次方程精确根
        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 前四位一致。 $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 本来只有 5% 的 AlO 解离，这里却解离了 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)}$ 。双分子反应移植到 SI 时要乘 $10^{-6}$ 。三体反应多一个 $[M]$ ，因子是 $10^{-12}$ 。反应级数不同，因子也不同。

所以单位错误总是和反应级数或 $\Delta n$ 一起出现。如果所有反应一起错，反而容易发现。只错一部分时，人们会去怀疑机理本身。报告里写下"需要重新确认逆反应速率常数与 Kp 的关系"，方向完全正确。

怀疑某个逆反应时可以跑的两行

只需要一个平衡计算工具。在相同 T、p 下用 CEA 或 Cantera 求平衡组成。把自己的代码当作等温反应器长时间运行，然后逐个反应输出 $Q/K_c$ 。

所有反应都是 1，就算通过。某个反应停在 $RT/p^\circ$ 或其倒数，说明直接用了 $K_p$ 。若在 $10^5$ ~ $10^6$ 附近，说明还叠加了 CGS 换算。数字本身就指向原因。在逐个检查十二个反应之前，先做出这张表。

用 Kp 去除逆反应速率常数，2500 K 下 AlO 竟解离了 99.6% — 平衡常数里藏着的 (p°/RT)^Δn