""" Steam methane reforming — autothermal reformer (ATR) bdrworks.org · Handbook · Steam reforming, worked model 3 Heat is generated inside the reactor: oxygen burns part of the methane in a burner zone, and the hot gas then reforms over an adiabatic nickel catalyst bed. burner (assumed complete, adiabatic): CH4 + 3/2 O2 -> CO + 2 H2O catalyst bed (adiabatic, Xu & Froment kinetics): dF_i/dW = sum_j nu_ij eta_j r_j, (sum F_i cp_i) dT/dW = sum_j (-dH_j) eta_j r_j Run: python atr.py (needs numpy, scipy, matplotlib) """ import numpy as np from scipy.integrate import solve_ivp from scipy.optimize import brentq import matplotlib.pyplot as plt # ---------------------------------------------------------------- data SPECIES = ["CH4", "H2O", "CO", "CO2", "H2", "N2", "O2"] # CH4 H2O CO CO2 H2 N2 O2 NU = np.array([[-1.0, -1.0, 1.0, 0.0, 3.0, 0.0, 0.0], # 1: CH4 + H2O <=> CO + 3 H2 [0.0, -1.0, -1.0, 1.0, 1.0, 0.0, 0.0], # 2: CO + H2O <=> CO2 + H2 [-1.0, -2.0, 0.0, 1.0, 4.0, 0.0, 0.0]]) # 3: CH4 + 2 H2O <=> CO2 + 4 H2 RG = 8.314 # J/(mol K) MW = np.array([16.043, 18.015, 28.010, 44.009, 2.016, 28.013, 31.999]) # kg/kmol # NASA 7-coefficient polynomials (GRI-Mech 3.0 thermodynamic data): [low T < 1000 K], [high T] NASA_LO = np.array([ [5.14987613e+00, -1.36709788e-02, 4.91800599e-05, -4.84743026e-08, 1.66693956e-11, -1.02466476e+04, -4.64130376e+00], # CH4 [4.19864056e+00, -2.03643410e-03, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12, -3.02937267e+04, -8.49032208e-01], # H2O [3.57953347e+00, -6.10353680e-04, 1.01681433e-06, 9.07005884e-10, -9.04424499e-13, -1.43440860e+04, 3.50840928e+00], # CO [2.35677352e+00, 8.98459677e-03, -7.12356269e-06, 2.45919022e-09, -1.43699548e-13, -4.83719697e+04, 9.90105222e+00], # CO2 [2.34433112e+00, 7.98052075e-03, -1.94781510e-05, 2.01572094e-08, -7.37611761e-12, -9.17935173e+02, 6.83010238e-01], # H2 [3.29867700e+00, 1.40824040e-03, -3.96322200e-06, 5.64151500e-09, -2.44485400e-12, -1.02089990e+03, 3.95037200e+00], # N2 [3.78245636e+00, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12, -1.06394356e+03, 3.65767573e+00], # O2 ]) NASA_HI = np.array([ [7.48514950e-02, 1.33909467e-02, -5.73285809e-06, 1.22292535e-09, -1.01815230e-13, -9.46834459e+03, 1.84373180e+01], # CH4 [3.03399249e+00, 2.17691804e-03, -1.64072518e-07, -9.70419870e-11, 1.68200992e-14, -3.00042971e+04, 4.96677010e+00], # H2O [2.71518561e+00, 2.06252743e-03, -9.98825771e-07, 2.30053008e-10, -2.03647716e-14, -1.41518724e+04, 7.81868772e+00], # CO [3.85746029e+00, 4.41437026e-03, -2.21481404e-06, 5.23490188e-10, -4.72084164e-14, -4.87591660e+04, 2.27163806e+00], # CO2 [3.33727920e+00, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14, -9.50158922e+02, -3.20502331e+00], # H2 [2.92664000e+00, 1.48797680e-03, -5.68476000e-07, 1.00970380e-10, -6.75335100e-15, -9.22797700e+02, 5.98052800e+00], # N2 [3.28253784e+00, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, -1.08845772e+03, 5.45323129e+00], # O2 ]) def _a(T): return NASA_LO if T < 1000.0 else NASA_HI def cp(T): """Ideal-gas heat capacities, J/(mol K).""" a = _a(T) return RG * (a[:, 0] + a[:, 1]*T + a[:, 2]*T**2 + a[:, 3]*T**3 + a[:, 4]*T**4) def h(T): """Molar enthalpies including heats of formation, J/mol.""" a = _a(T) return RG * T * (a[:, 0] + a[:, 1]*T/2 + a[:, 2]*T**2/3 + a[:, 3]*T**3/4 + a[:, 4]*T**4/5 + a[:, 5]/T) def dH_rxn(T): """Heats of reactions 1-3, J/mol.""" return NU @ h(T) def K_eq(T): """Equilibrium constants of Xu & Froment: K1, K3 in bar^2, K2 dimensionless.""" K1 = np.exp(-26830.0 / T + 30.114) K2 = np.exp(4400.0 / T - 4.036) return K1, K2, K1 * K2 def rates(T, p, eta=(1.0, 1.0, 1.0)): """Xu & Froment (1989) intrinsic rates times effectiveness factors, kmol/(kg_cat h); p in bar.""" pCH4, pH2O, pCO, pCO2, pH2 = [max(x, 1e-10) for x in p[:5]] RT = RG * T / 1000.0 # kJ/mol k1 = 4.225e15 * np.exp(-240.1 / RT) k2 = 1.955e6 * np.exp(-67.13 / RT) k3 = 1.020e15 * np.exp(-243.9 / RT) KCO = 8.23e-5 * np.exp(70.65 / RT) KH2 = 6.12e-9 * np.exp(82.90 / RT) KCH4 = 6.65e-4 * np.exp(38.28 / RT) KH2O = 1.77e5 * np.exp(-88.68 / RT) K1, K2, K3 = K_eq(T) den = 1.0 + KCO*pCO + KH2*pH2 + KCH4*pCH4 + KH2O*pH2O/pH2 r1 = k1 / pH2**2.5 * (pCH4*pH2O - pH2**3*pCO/K1) / den**2 r2 = k2 / pH2 * (pCO*pH2O - pH2*pCO2/K2) / den**2 r3 = k3 / pH2**3.5 * (pCH4*pH2O**2 - pH2**4*pCO2/K3) / den**2 return np.array([eta[0]*r1, eta[1]*r2, eta[2]*r3]) def equilibrium(T, P, n0): """Equilibrium of reactions 1 and 2 (3 is their sum) from flows n0 at T (K), P (bar).""" from scipy.optimize import fsolve K1, K2, _ = K_eq(T) def res(x): n = np.maximum(n0 + NU[0]*x[0] + NU[1]*x[1], 1e-14) p = P * n / n.sum() return [np.log(p[2]*p[4]**3 / (p[0]*p[1])) - np.log(K1), np.log(p[3]*p[4] / (p[2]*p[1])) - np.log(K2)] x = fsolve(res, [0.7*n0[0], 0.2*n0[0]], xtol=1e-12) return n0 + NU[0]*x[0] + NU[1]*x[1] # ---------------------------------------------------------------- ATR P = 30.0 # bar F_CH4 = 1000.0 # natural gas (as CH4), kmol/h S_C = 0.6 # steam-to-carbon ratio (low, as in modern ATRs) T_GAS, T_O2 = 900.0, 500.0 # preheat of CH4 + steam, and of the oxygen, K W_CAT = 20e3 # catalyst, kg ETA = (0.03, 0.03, 0.03) BURN = np.array([-1.0, 2.0, 1.0, 0.0, 0.0, 0.0, -1.5]) # CH4 + 1.5 O2 -> CO + 2 H2O def burner(o2_c): """Adiabatic flame: all O2 burns CH4 to CO + H2O. Returns flows (kmol/h) and temperature (K).""" gas = np.array([F_CH4, S_C * F_CH4, 0, 0, 0.02 * F_CH4, 0.01 * F_CH4, 0]) o2 = np.zeros(7); o2[6] = o2_c * F_CH4 H_in = np.dot(gas, h(T_GAS)) + np.dot(o2, h(T_O2)) n = gas + o2 + BURN * (o2_c * F_CH4 / 1.5) T = brentq(lambda T: np.dot(n, h(T)) - H_in, 600, 4000) return n, T def rhs(W, u): F, T = u[:7], u[7] r = rates(T, P * F / F.sum(), ETA) # kmol/(kg h) return np.append(NU.T @ r, np.dot(-dH_rxn(T), r) * 1e3 / np.dot(F * 1e3, cp(T))) def atr(o2_c, n=300): n0, Tb = burner(o2_c) sol = solve_ivp(rhs, (0, W_CAT), np.append(n0, Tb), method="LSODA", rtol=1e-8, atol=1e-8, dense_output=True) w = np.linspace(0, W_CAT, n) return w, sol.sol(w), Tb def summary(u): F = u[:7, -1] return dict(T=u[7, -1], X=1 - F[0] / F_CH4, slip=100 * F[0] / (F.sum() - F[1]), h2co=F[4] / F[2], syngas=F[4] + F[2]) if __name__ == "__main__": o2c = 0.57 w, u, Tb = atr(o2c) s = summary(u) ne = equilibrium(s["T"], P, u[:7, 0]) print(f"ATR {P:.0f} bar, S/C {S_C}, O2/C {o2c}: burner exit {Tb:.0f} K, bed outlet {s['T']:.0f} K") print(f"CH4 conversion {100*s['X']:.1f} % (equilibrium at outlet {100*(1-ne[0]/F_CH4):.1f} %), " f"slip {s['slip']:.2f} % dry, H2/CO {s['h2co']:.2f}, H2+CO {s['syngas']:.0f} kmol/h") o2s = np.arange(0.50, 0.681, 0.02) rows = [summary(atr(x, 60)[1]) for x in o2s] for x, r in zip(o2s, rows): print(f" O2/C {x:.2f}: outlet {r['T']:.0f} K, slip {r['slip']:.2f} %, H2/CO {r['h2co']:.2f}") fig, ax = plt.subplots(1, 3, figsize=(14, 4.2)) ax[0].plot(w / 1e3, u[7], lw=2); ax[0].set(xlabel="catalyst (t)", ylabel="temperature (K)", title=f"Bed temperature, O$_2$/C = {o2c}") dry = u[:7].sum(axis=0) - u[1] for i, lab in [(0, "CH$_4$"), (4, "H$_2$"), (2, "CO"), (3, "CO$_2$")]: ax[1].plot(w / 1e3, 100 * u[i] / dry, lw=2, label=lab) ax[1].set(xlabel="catalyst (t)", ylabel="mol-% (dry)", title="Composition"); ax[1].legend() ax[2].plot(o2s, [r["slip"] for r in rows], "o-", lw=2, label="CH$_4$ slip (% dry)") ax[2].plot(o2s, [r["h2co"] for r in rows], "s-", lw=2, label="H$_2$/CO") ax[2].set(xlabel="O$_2$/C (mol/mol)", title="Oxygen demand"); ax[2].legend() a2 = ax[2].twinx(); a2.plot(o2s, [r["T"] for r in rows], "--", c="C3"); a2.set_ylabel("outlet T (K)", color="C3") for a in ax[:2]: a.set_xlim(0, 2) # equilibrium is reached in the first tonnes ax[0].set_xlabel("catalyst (t) — first 2 t of the bed"); ax[1].set_xlabel("catalyst (t) — first 2 t of the bed") for a in ax: a.grid(alpha=0.3) fig.tight_layout(); fig.savefig("atr.png", dpi=150) plt.show()