"""
Steam methane reforming — one fired reformer tube
bdrworks.org · Handbook · Steam reforming, worked model 2

Catalyst tube in a top-fired furnace. Heat enters through the tube wall with a
prescribed flux profile; the gas is heated and reformed on Ni/Al2O3 catalyst.
    dF_i/dz = rho_b A sum_j nu_ij eta_j r_j                        (Xu & Froment, 1989)
    (sum F_i cp_i) dT/dz = pi d_o q(z) + rho_b A sum_j (-dH_j) eta_j r_j
    dP/dz  = -(Ergun)
Tube-wall temperatures follow from the inside film coefficient (Leva) and wall conduction.

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

# ---------------------------------------------------------------- tube
L, D_I, D_O = 12.0, 0.1016, 0.1320     # heated length, inside and outside diameter, m
A = np.pi * D_I**2 / 4                 # flow area, m2
RHO_B = 1360.0                         # catalyst bulk density, kg/m3
D_P, EPS = 0.0173, 0.50                # equivalent pellet diameter (m), bed voidage
ETA = (0.03, 0.03, 0.03)               # effectiveness factors of reactions 1-3
MU, LAM_G, K_W = 3.0e-5, 0.12, 25.0    # gas viscosity (Pa s), gas conductivity and wall conductivity (W/m K)
P0, T0 = 29.0, 793.0                   # inlet pressure (bar) and temperature (K)
F_CH4, S_C = 6.0, 3.0                  # methane per tube (kmol/h), steam-to-carbon ratio
F0 = np.array([F_CH4, S_C * F_CH4, 0.0, 0.1, 0.3, 0.1, 0.0]) / 3.6   # mol/s per tube (H2 keeps Ni reduced)


def q_flux(z, q_avg):
    """Heat flux on the outer tube surface, W/m2: highest near the burners at the top."""
    return q_avg * (1.3 - 0.6 * z / L)


def rhs(z, u, q_avg):
    F, T, P = u[:7], u[7], u[8]
    y = F / F.sum()
    r = rates(T, P * y, ETA) / 3.6                                # mol/(kg s)
    dF = RHO_B * A * (NU.T @ r)
    Q = np.pi * D_O * q_flux(z, q_avg)                            # W per m of tube
    dT = (Q + RHO_B * A * np.dot(-dH_rxn(T), r)) / np.dot(F, cp(T))
    G = np.dot(F, MW) / 1000 / A                                  # kg/(m2 s)
    rho = P * 1e5 * np.dot(y, MW) / 1000 / (RG * T)
    us = G / rho
    dP = -(150 * MU * (1 - EPS)**2 / (EPS**3 * D_P**2) * us + 1.75 * rho * (1 - EPS) / (EPS**3 * D_P) * us**2) / 1e5
    return np.concatenate([dF, [dT, dP]])


def wall_T(z, T, q_avg):
    """Inner and outer tube-wall temperatures, K (Leva correlation for heating packed tubes)."""
    G = np.dot(F0, MW) / 1000 / A
    h_i = 0.813 * LAM_G / D_I * (D_P * G / MU)**0.9 * np.exp(-6 * D_P / D_I)
    q = q_flux(z, q_avg)
    Twi = T + q * D_O / (h_i * D_I)
    Two = Twi + q * D_O / 2 * np.log(D_O / D_I) / K_W
    return Twi, Two, h_i


def solve(q_avg, n=300):
    sol = solve_ivp(rhs, (0, L), np.concatenate([F0, [T0, P0]]), args=(q_avg,), method="LSODA",
                    rtol=1e-8, atol=1e-10, dense_output=True)
    z = np.linspace(0, L, n)
    return z, sol.sol(z)


def outlet(u):
    F = u[:7, -1]
    X = 1 - F[0] / F0[0]
    dry = F.sum() - F[1]
    return X, 100 * F[0] / dry, u[7, -1], u[8, -1], F[4] / F[2]


if __name__ == "__main__":
    q = 75e3                                                      # average heat flux, W/m2
    z, u = solve(q)
    Twi, Two, h_i = wall_T(z, u[7], q)
    X, slip, Tout, Pout, h2co = outlet(u)
    ne = equilibrium(Tout, Pout, F0)
    Xeq = 1 - ne[0] / F0[0]
    print(f"tube {L} m x {1000*D_I:.1f} mm, S/C {S_C}, {F_CH4} kmol/h CH4, average flux {q/1e3:.0f} kW/m2, h_i = {h_i:.0f} W/m2K")
    print(f"outlet {Tout:.1f} K, {Pout:.2f} bar, CH4 conversion {100*X:.1f} % (equilibrium at outlet {100*Xeq:.1f} %), "
          f"CH4 slip {slip:.2f} % dry, H2/CO {h2co:.2f}")
    print(f"max outer wall temperature {Two.max():.0f} K at z = {z[Two.argmax()]:.1f} m")

    fig, ax = plt.subplots(1, 3, figsize=(14, 4.2))
    ax[0].plot(z, u[7], lw=2, label="process gas")
    ax[0].plot(z, Twi, lw=2, label="inner wall")
    ax[0].plot(z, Two, lw=2, label="outer wall")
    ax[0].set(xlabel="tube length z (m)", ylabel="temperature (K)", title="Temperatures, 75 kW/m$^2$"); ax[0].legend()
    Fz = u[:7]
    dry = Fz.sum(axis=0) - Fz[1]
    for i, s in [(0, "CH$_4$"), (4, "H$_2$"), (2, "CO"), (3, "CO$_2$")]:
        ax[1].plot(z, 100 * Fz[i] / dry, lw=2, label=s)
    ax[1].set(xlabel="tube length z (m)", ylabel="mol-% (dry)", title="Composition"); ax[1].legend()
    qs = np.arange(50e3, 101e3, 5e3)
    res = []
    for qq in qs:
        zz, uu = solve(qq)
        res.append(outlet(uu) + (wall_T(zz, uu[7], qq)[1].max(),))
    ax[2].plot(qs / 1e3, [100 * r[0] for r in res], lw=2, c="C0", label="CH$_4$ conversion")
    ax[2].set(xlabel="average heat flux (kW/m$^2$)", ylabel="CH$_4$ conversion (%)", title="Firing rate")
    a2 = ax[2].twinx(); a2.plot(qs / 1e3, [r[5] for r in res], lw=2, c="C3", ls="--"); a2.set_ylabel("max outer wall T (K)", color="C3")
    for a in ax:
        a.grid(alpha=0.3)
    fig.tight_layout()
    fig.savefig("smr_tube.png", dpi=150)
    print("flux sweep (kW/m2: conversion %, outlet K, max wall K):",
          ", ".join(f"{qq/1e3:.0f}: {100*r[0]:.1f}/{r[2]:.0f}/{r[5]:.0f}" for qq, r in zip(qs, res)))
    plt.show()