Dividir la constante inversa entre Kp disoció el 99.6% del AlO a 2500 K — el (p°/RT)^Δn oculto en la constante de equilibrio

La constante inversa sale de kf/Kc. Como Kc = Kp(p°/RT)^Δn, usar Kp no hace daño si Δn=0, pero cuando cambia el número de moles el error va de 4 a 250 000 veces según el sistema de unidades.

Chemkin y el código propio difirieron en 134 K#

Existe una presentación de laboratorio sobre un modelo propio de combustión de partículas de aluminio. Doce reacciones en fase gas entraron como expresiones de Arrhenius. El mismo mecanismo se corrió también en Chemkin. La temperatura máxima fue 3800 K en un caso y 3666 K en el otro. Al comparar las velocidades una por una, solo las constantes inversas de R2–R5 estaban mal. La siguiente diapositiva se titulaba "Constante de equilibrio".

Este artículo completa esa diapositiva. La fórmula que construye la constante inversa esconde una unidad. Si se omite, algunas reacciones siguen correctas y otras fallan. Un único reactor de disociación de AlO reproduce el error y deja una regla para identificar qué reacciones están en riesgo.

La constante inversa no viene por separado#

Un archivo de mecanismo suele traer solo los coeficientes directos: $A$ , $b$ y $E_a$ en $k_f = A T^b \exp(-E_a/RT)$ . La reacción inversa se obtiene del balance detallado (detailed balance: en equilibrio, las velocidades directa e inversa son iguales).

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

$k_r$ es la constante inversa, $K_c$ la constante de equilibrio en concentraciones y $\nu_i$ el coeficiente estequiométrico, positivo en productos y negativo en reactivos.

Sin embargo, la termodinámica no entrega $K_c$ directamente. Lo que se calcula con la energía libre de Gibbs es la $K_p$ basada en presiones.

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$ es la presión estándar (1 atm), y $\Delta H^\circ$ , $\Delta S^\circ$ son la entalpía y la entropía de reacción. $K_p$ es adimensional. $K_c$ tiene unidades de (concentración) $^{\Delta n}$ .

Entre Kp y Kc está (p°/RT)^Δn#

En un gas ideal, $p_i = [X_i]RT$ . Al sustituirlo en $K_p$ aparece la relación entre ambas constantes.

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

$\Delta n$ es el cambio neto de moles de gas. $p^\circ/RT$ es la concentración molar de un gas en el estado estándar.

Este factor actúa distinto en cada reacción. $\mathrm{Al + O_2 \rightleftharpoons AlO + O}$ tiene dos moles a cada lado. Con $\Delta n = 0$ resulta $K_c = K_p$ , así que usar $K_p$ tal cual sigue siendo correcto. $\mathrm{AlO + M \rightleftharpoons Al + O + M}$ convierte una molécula en dos, de modo que $\Delta n = +1$ . A 3000 K, $p^\circ/RT = 4.06\ \mathrm{mol/m^3}$ , por lo que un código en SI subestima $K_c$ unas 4 veces. En un código estilo Chemkin con concentraciones en $\mathrm{mol/cm^3}$ , el mismo factor vale $4.06\times10^{-6}$ . Ahora $K_c$ queda 250 000 veces más grande.

La caja de abajo permite manipularlo directamente.

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

Con el botón $\Delta n = 0$ , los tres modos producen la misma caja. Con $\Delta n = +1$ y el modo CGS, casi todas las líneas de enlace se rompen. Al mover la presión, el equilibrio verdadero cambia, pero la marca de $K_p$ no se mueve.

Por qué solo fallaban R2–R5#

La presentación no indica la causa. Pero el síntoma, "solo algunas constantes inversas estaban mal", encaja exactamente con este factor. Aunque se cometa el mismo error en todas las reacciones, las de $\Delta n = 0$ no se ven afectadas. Solo fallan las que cambian el número de moles, como la disociación y la recombinación.

El mecanismo de aluminio en fase gas incluye justamente reacciones así. La disociación de AlO + M y la formación de $\mathrm{Al_2O}$ son reacciones de tercer cuerpo (third body: una tercera molécula que solo intercambia energía de colisión). La nota "Third body=?" de la diapositiva quizá no sea casualidad.

La segunda trampa es la dirección. Usar $K_p$ en SI deja $K_c$ demasiado pequeña: se disocia de menos. Usarla en CGS la deja demasiado grande: se disocia de más. El mismo error se desvía en sentidos opuestos según el sistema de unidades. Por eso tampoco se puede anticipar si la temperatura máxima saldrá alta o baja.

Código: 1% de AlO en un reactor a 3000 K#

Un reactor isotérmico a volumen constante contiene 1% de AlO diluido en argón. Se resuelve solo la disociación, con la constante inversa construida de tres maneras. Los datos termodinámicos se fijan en valores aproximados a 298 K. El factor preexponencial $A$ es ilustrativo. El objetivo no es un cálculo de combustión exacto, sino medir el tamaño del error de unidades.

El avance temporal es Euler implícito (backward Euler: el lado derecho se evalúa en el instante nuevo). La actualización es cuadrática en $\xi$ , así que la fórmula cuadrática resuelve cada paso de forma exacta. Nada diverge aunque la reacción inversa se acelere 250 000 veces. Cómo el avance implícito domina los términos rígidos también aparece en el artículo sobre el esquema Gear semi-implícito.

import math
 
R = 8.314462618        # J/(mol K)
P0 = 101325.0          # presión estándar, Pa
 
# entalpía/entropía de reacción a 298 K fijas (valores JANAF aproximados)
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))   # adimensional
 
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 leído como mol/m^3
    if mode == "Kp as Kc, CGS":
        return kp_of(rx, T) * 1e6 ** rx["dn"]            # Kp leído como mol/cm^3
 
def k_forward(T):
    return 1.0e11 * math.exp(-DISSOC["DH"] / (R * T))    # m^3/(mol s), valor ilustrativo
 
def run_reactor(T, p, x_alo, mode, t_end=0.2, nsteps=500):
    M = p / (R * T)                                  # diluido en Ar: [M] ~ constante
    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 implícito, raíz exacta de la cuadrática
        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]

El mismo código puede correrse directamente en el navegador.

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

Con los valores por defecto, las tres curvas llegan a 0.2997, 0.1626 y 1.0000, como en el bloque de salida. Al bajar T a 2500 K, la curva verde se queda cerca de 0.05 mientras solo la roja se dispara hasta 1.

Lo que dicen 0.2997, 0.1626 y 1.0000#

El modo correcto coincide con el equilibrio exacto de 0.2997 en cuatro cifras. $Q/K_c$ (el cociente de reacción dividido entre la constante verdadera) también vale 1.000. El reactor alcanzó el equilibrio real.

El error en SI se detiene en $Q/K_c = 0.2462$ . Ese número es $1/4.062$ , es decir, $RT/p^\circ$ . El reactor convergió con exactitud al equilibrio equivocado. El residuo es cero y no hay ninguna advertencia. Ningún criterio de convergencia detecta este error.

El $2.462\times10^5$ del error en CGS es $10^6/4.062$ . El factor de conversión $10^6$ empuja el equilibrio sin más. A 2500 K debería disociarse solo el 5% del AlO. Aquí se disoció el 99.6%. La disociación es endotérmica, así que en un cálculo de combustión real el gas se enfriaría mucho. Son posibles errores muy superiores a los 134 K de la presentación.

Al cambiar la presión, el error se ve distinto#

El equilibrio verdadero de una reacción con $\Delta n = +1$ depende de la presión: 0.66 a 0.1 atm y 0.11 a 10 atm. Como predice Le Chatelier, la presión baja favorece el lado con más moles.

El error en SI sigue esa tendencia, solo que más abajo. A cualquier presión su $K_c$ es $RT/p^\circ$ veces menor. Como la tendencia es correcta, la gráfica parece verosímil. Frente a experimentos, es el tipo de error que termina absorbido en un "ajuste de constantes del modelo".

El error en CGS casi borra la dependencia con la presión. Incluso a 10 atm marca 0.9997. Si una simulación de cámara de combustión a alta presión muestra esa insensibilidad, conviene revisar primero las unidades. El flujo de calor hacia una partícula de aluminio también se separa según presión y tamaño; eso se trata en el artículo sobre el régimen de transición de Knudsen.

Dónde se esconden las unidades al migrar una tabla de constantes#

No es solo la conversión de $K_c$ . La $A$ en formato Chemkin está en $\mathrm{cm^3/(mol\,s)}$ . Para una reacción bimolecular, pasar a SI exige multiplicar por $10^{-6}$ . Una reacción de tercer cuerpo lleva un $[M]$ más, así que el factor es $10^{-12}$ . El factor cambia con el orden de reacción.

Por eso los errores de unidades siempre aparecen ligados al orden de reacción o a $\Delta n$ . Si todas las reacciones fallaran a la vez, sería fácil encontrarlo. Cuando solo fallan algunas, se termina dudando del mecanismo mismo. La nota de la presentación, "hay que revisar la relación entre la constante inversa y Kp", apuntaba justo en la dirección correcta.

Dos líneas para correr cuando una constante inversa resulta sospechosa#

Basta con una herramienta de equilibrio. Se obtiene la composición de equilibrio con CEA o Cantera a la misma T y p. Luego se corre el código propio como reactor isotérmico durante largo tiempo y se imprime $Q/K_c$ para cada reacción.

Si todas valen 1, el código pasa. Si alguna se detiene en $RT/p^\circ$ o su inverso, usó $K_p$ tal cual. Si ronda $10^5$ – $10^6$ , además se le sumó la conversión CGS. El número mismo señala la causa. Antes de revisar doce reacciones una por una, conviene armar primero esa tabla.