역반응 상수를 Kp로 나눴더니 2500 K에서 AlO가 99.6% 쪼개졌다 — 평형상수의 단위 (p°/RT)^Δn

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$은 역반응 상수, $K_c$는 농도 기반 평형상수, $\nu_i$는 생성물 양수·반응물 음수인 화학량론 계수다.

그런데 열역학이 바로 주는 값은 $K_c$가 아니다. Gibbs 자유에너지로 계산하는 값은 압력 기반 $K_p$다.

$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$에 넣으면 두 상수의 관계가 나온다.

$\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만 배 커진다.

아래 상자에서 직접 조작해보자.

$\Delta n = 0$ 버튼을 누르면 세 모드의 상자가 똑같아진다. $\Delta n = +1$에서 CGS 모드를 고르면 결합선이 거의 다 끊긴다. 압력 슬라이더를 움직이면 참 평형은 변하지만 $K_p$ 눈금은 그대로다.

R2~R5만 틀린 이유#

발표 자료는 원인을 적지 않았다. 그러나 "일부 반응의 역반응만 틀렸다"는 증상은 이 인자와 정확히 맞는다. 모든 반응에 같은 실수를 해도 $\Delta n = 0$ 반응은 영향을 받지 않는다. 틀리는 것은 해리·재결합처럼 몰수가 바뀌는 반응뿐이다.

알루미늄 기상 메커니즘에는 그런 반응이 섞여 있다. $\mathrm{AlO + M}$ 해리, $\mathrm{Al_2O}$ 생성 같은 3체(third body, 충돌 에너지만 주고받는 제3의 분자) 반응이 그렇다. 슬라이드에 "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]

위 코드를 브라우저에서 직접 돌려보자.

기본값에서 세 곡선은 출력 블록의 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$다. 반응기는 잘못된 평형에 정확히 수렴했다. 잔차는 0이고 경고도 없다. 수렴 판정으로는 이 실수를 잡을 수 없다는 뜻이다.

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}$을 곱한다. 3체 반응은 $[M]$이 하나 더 붙으니 $10^{-12}$이다. 차수마다 인자가 달라진다.

그래서 단위 실수는 늘 반응 차수나 $\Delta n$과 함께 나타난다. 모든 반응이 같이 틀리면 오히려 찾기 쉽다. 일부만 틀리면 메커니즘 자체를 의심하게 된다. 발표 자료가 "역반응상수와 Kp의 관계를 다시 확인할 필요"라고 쓴 것은 정확한 방향이다.

역반응 하나가 의심스러울 때 돌려볼 두 줄#

평형 계산기 하나만 있으면 된다. 같은 T·p에서 CEA나 Cantera로 평형 조성을 구한다. 내 코드를 등온 반응기로 오래 돌린 뒤 반응마다 $Q/K_c$를 찍는다.

모든 반응이 1이면 통과다. 어떤 반응이 $RT/p^\circ$나 그 역수로 멈추면 $K_p$를 그대로 쓴 것이다. $10^5$~$10^6$ 근처라면 CGS 변환까지 겹친 것이다. 숫자 자체가 원인을 가리킨다. 반응 열두 개를 하나씩 들여다보기 전에 이 표 하나를 먼저 만들어보자.