"""
Steam methane reforming — chemical equilibrium (the thermodynamic limit)
bdrworks.org · Handbook · Steam reforming, worked model 1

Equilibrium of CH4 + H2O = CO + 3H2 and CO + H2O = CO2 + H2 from a methane/steam feed
as a function of temperature, pressure and steam-to-carbon ratio (ideal gas).
Every reformer, whatever its design, is bounded by these curves.

Run:  python smr_equilibrium.py      (needs numpy, scipy, matplotlib)
"""
import numpy as np
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]


def feed(sc):
    return np.array([1.0, sc, 0.0, 0.0, 0.05, 0.0, 0.0])


def summary(n, n0):
    dry = n.sum() - n[1]
    return 1 - n[0] / n0[0], 100 * n[0] / dry, n[4] / n[2], n[4] / n0[0]


if __name__ == "__main__":
    Ts = np.linspace(800, 1300, 51)
    fig, ax = plt.subplots(1, 3, figsize=(14, 4.2))
    for sc in [2.0, 3.0, 4.0]:
        n0 = feed(sc)
        res = [summary(equilibrium(t, 29.0, n0), n0) for t in Ts]
        ax[0].plot(Ts, [100 * r[0] for r in res], lw=2, label=f"S/C = {sc:.0f}")
        ax[1].plot(Ts, [r[1] for r in res], lw=2, label=f"S/C = {sc:.0f}")
    for P in [1.0, 10.0, 29.0, 40.0]:
        n0 = feed(3.0)
        ax[2].plot(Ts, [100 * summary(equilibrium(t, P, n0), n0)[0] for t in Ts], lw=2, label=f"{P:.0f} bar")
    ax[0].set(xlabel="temperature (K)", ylabel="CH$_4$ conversion (%)", title="Equilibrium conversion, 29 bar")
    ax[1].set(xlabel="temperature (K)", ylabel="CH$_4$ slip (mol-%, dry)", title="Methane slip, 29 bar", yscale="log")
    ax[2].set(xlabel="temperature (K)", ylabel="CH$_4$ conversion (%)", title="Effect of pressure, S/C = 3")
    for a in ax:
        a.grid(alpha=0.3, which="both"); a.legend()
    fig.tight_layout(); fig.savefig("smr_equilibrium.png", dpi=150)
    for sc in [2.0, 3.0, 4.0]:
        n0 = feed(sc); X, slip, h2co, h2 = summary(equilibrium(1123.0, 29.0, n0), n0)
        print(f"1123 K, 29 bar, S/C {sc:.0f}: CH4 conversion {100*X:.1f} %, slip {slip:.2f} % dry, H2/CO {h2co:.2f}, H2 per CH4 {h2:.2f}")
    for P in [1.0, 10.0, 40.0]:
        n0 = feed(3.0); print(f"1123 K, {P:.0f} bar, S/C 3: CH4 conversion {100*summary(equilibrium(1123.0, P, n0), n0)[0]:.1f} %")
    plt.show()