""" Methanol synthesis — boiling-water-cooled multitubular reactor (Lurgi type) bdrworks.org · Handbook · Methanol, worked model 2 Reactions on Cu/ZnO/Al2O3 (Vanden Bussche & Froment, 1996): 1: CO2 + 3 H2 <=> CH3OH + H2O 2: CO2 + H2 <=> CO + H2O (reverse water-gas shift) One representative tube, 1-D pseudo-homogeneous plug flow, shell side at the boiling temperature of the cooling water (constant Tc): dF_i/dz = rho_b A_t * sum_j nu_ij r_j (sum F_i cp_i) dT/dz = rho_b A_t * sum_j (-dH_j) r_j - U pi D (T - Tc) Run: python meoh_tubular.py (needs numpy, scipy, matplotlib) """ import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # ---------------------------------------------------------------- data SPECIES = ["CO", "CO2", "H2", "CH3OH", "H2O", "CH4", "N2"] # CO CO2 H2 MeOH H2O CH4 N2 NU = np.array([[0.0, -1.0, -3.0, 1.0, 1.0, 0.0, 0.0], # 1: CO2 + 3H2 <=> CH3OH + H2O [1.0, -1.0, -1.0, 0.0, 1.0, 0.0, 0.0]]) # 2: CO2 + H2 <=> CO + H2O (RWGS) RG = 8.314 # J/(mol K) # ideal-gas cp = A + B T + C T^2 + D T^3, J/(mol K) (Reid, Prausnitz & Poling) CP = np.array([[30.87, -1.285e-2, 2.789e-5, -1.272e-8], # CO [19.80, 7.344e-2, -5.602e-5, 1.715e-8], # CO2 [27.14, 9.274e-3, -1.381e-5, 7.645e-9], # H2 [21.15, 7.092e-2, 2.587e-5, -2.852e-8], # CH3OH [32.24, 1.924e-3, 1.055e-5, -3.596e-9], # H2O [19.25, 5.213e-2, 1.197e-5, -1.132e-8], # CH4 [31.15, -1.357e-2, 2.680e-5, -1.168e-8]]) # N2 HF = np.array([-110.53, -393.51, 0.0, -200.94, -241.83, -74.87, 0.0]) * 1e3 # J/mol, 298.15 K def cp(T): """Heat capacities of all species, J/(mol K).""" return CP @ np.array([1.0, T, T**2, T**3]) def h(T): """Molar enthalpies H_f(298) + integral of cp, J/mol.""" T0 = 298.15 dT = np.array([T - T0, (T**2 - T0**2) / 2, (T**3 - T0**3) / 3, (T**4 - T0**4) / 4]) return HF + CP @ dT def dH_rxn(T): """Heats of reaction 1 and 2, J/mol.""" return NU @ h(T) def K_eq(T): """Equilibrium constants (Graaf et al.), partial pressures in bar. K1 for CO2 + 3H2 = CH3OH + H2O (bar^-2); K3 for CO + H2O = CO2 + H2 (-).""" K1 = 10.0 ** (3066.0 / T - 10.592) K3 = 10.0 ** (2073.0 / T - 2.029) return K1, K3 def rates(T, p): """Vanden Bussche & Froment (1996) rates, mol/(kg_cat s); p = partial pressures in bar.""" pCO, pCO2, pH2, pM, pW = [max(x, 1e-10) for x in p[:5]] k = lambda A, B: A * np.exp(B / (RG * T)) k1, k2, k3, k4, k5 = k(1.07, 36696), k(3453.38, 0), k(0.499, 17197), k(6.62e-11, 124119), k(1.22e10, -94765) K1, K3 = K_eq(T) den = 1.0 + k2 * pW / pH2 + k3 * np.sqrt(pH2) + k4 * pW r1 = k1 * pCO2 * pH2 * (1.0 - pW * pM / (K1 * pH2**3 * pCO2)) / den**3 r2 = k5 * pCO2 * (1.0 - K3 * pW * pCO / (pCO2 * pH2)) / den return np.array([r1, r2]) def equilibrium(T, P, F): """Chemical equilibrium of reactions 1 and 2 from flows F (any unit); returns outlet flows.""" from scipy.optimize import fsolve K1, K3 = K_eq(T) def res(x): n = F + NU.T @ x n = np.maximum(n, 1e-12) p = P * n / n.sum() return [np.log(p[3] * p[4] / (p[1] * p[2]**3)) - np.log(K1), np.log(p[0] * p[4] / (p[1] * p[2])) + np.log(K3)] x = fsolve(res, [0.1 * F.sum() * 0.05, 0.0], xtol=1e-12) return F + NU.T @ x # ---------------------------------------------------------------- reactor P = 75.0 # bar F0 = 10000.0 # total feed, kmol/h Y0 = np.array([0.10, 0.05, 0.70, 0.005, 0.001, 0.12, 0.024]) Y0 = Y0 / Y0.sum() T_IN = 498.0 # gas inlet, K N_T = 5000 # number of tubes D = 0.040 # tube inside diameter, m L = 7.0 # tube length, m RHO_B = 1140.0 # catalyst bulk density, kg/m3 U = 250.0 # overall heat-transfer coefficient, W/(m2 K) A_T = np.pi * D**2 / 4 def rhs(z, u, Tc): F, T = u[:7], u[7] # F in mol/s per tube p = P * F / F.sum() r = rates(T, p) # mol/(kg s) dF = RHO_B * A_T * (NU.T @ r) q = U * np.pi * D * (T - Tc) # W per m of tube dT = (RHO_B * A_T * np.dot(-dH_rxn(T), r) - q) / np.dot(F, cp(T)) return np.append(dF, dT) def solve(Tc, T_in=T_IN, n=400): F_tube = F0 * 1000 / 3600 / N_T * Y0 # mol/s per tube sol = solve_ivp(rhs, (0, L), np.append(F_tube, T_in), args=(Tc,), method="LSODA", rtol=1e-8, atol=1e-10, dense_output=True) z = np.linspace(0, L, n) return z, sol.sol(z) def summary(u): F_in, F_out = u[:7, 0], u[:7, -1] meoh = (F_out[3] - F_in[3]) * N_T * 3600 / 1000 # kmol/h xc = 1 - (F_out[0] + F_out[1]) / (F_in[0] + F_in[1]) # carbon (CO+CO2) conversion return meoh, meoh * 32.04 * 24 / 1000, xc, F_out[3] / F_out.sum() if __name__ == "__main__": Tc = 523.0 # boiling water at ~40 bar z, u = solve(Tc) meoh, tpd, xc, y = summary(u) print(f"P = {P:.0f} bar, feed {F0:.0f} kmol/h at {T_IN:.0f} K, coolant {Tc:.0f} K, " f"{N_T} tubes x {L} m, catalyst {RHO_B*A_T*L*N_T/1000:.1f} t") print(f"hot spot {u[7].max():.1f} K at z = {z[u[7].argmax()]:.2f} m, outlet {u[7,-1]:.1f} K") print(f"MeOH out {100*y:.2f} mol-%, CO+CO2 conversion {100*xc:.1f} %, " f"production {meoh:.0f} kmol/h = {tpd:.0f} t/d") # parametric sensitivity: coolant temperature Tcs = np.arange(490, 551, 10) fig, ax = plt.subplots(1, 3, figsize=(14, 4.2)) for t in Tcs: zz, uu = solve(t) ax[0].plot(zz, uu[7], lw=1.8, label=f"T$_c$ = {t:.0f} K") ax[0].set(xlabel="tube length z (m)", ylabel="gas temperature (K)", title="Hot spot vs coolant temperature") ax[0].legend(fontsize=8) ax[1].plot(z, 100 * u[3] / u[:7].sum(axis=0), lw=2, label="CH$_3$OH") ax[1].plot(z, 100 * u[0] / u[:7].sum(axis=0), lw=2, label="CO") ax[1].plot(z, 100 * u[1] / u[:7].sum(axis=0), lw=2, label="CO$_2$") ax[1].plot(z, 100 * u[4] / u[:7].sum(axis=0), lw=2, label="H$_2$O") ax[1].set(xlabel="tube length z (m)", ylabel="mol-%", title=f"Composition, T$_c$ = {Tc:.0f} K") ax[1].legend() tc = np.arange(470, 551, 2.5) res = [] for t in tc: uu = solve(t)[1] res.append(summary(uu) + (uu[7].max(),)) best = tc[int(np.argmax([r[1] for r in res]))] print(f"coolant temperature for maximum production ≈ {best:.0f} K") ax[2].plot(tc, [r[1] for r in res], lw=2, c="C0") ax[2].set(xlabel="coolant temperature T$_c$ (K)", ylabel="methanol (t/d)", title="Production and hot spot") a2 = ax[2].twinx(); a2.plot(tc, [r[4] for r in res], lw=2, c="C3", ls="--"); a2.set_ylabel("hot-spot T (K)", color="C3") for a in ax: a.grid(alpha=0.3) fig.tight_layout() fig.savefig("meoh_tubular.png", dpi=150) plt.show()