"""
Nitric acid — NO oxidation in the nitrous-gas train of a dual-pressure plant
bdrworks.org · Handbook · Nitric acid, worked model 3

The gas leaving the waste-heat boiler oxidises in adiabatic runs of piping and is cooled in
the tail-gas preheater and the economizer (flow sheet of the dual-pressure process):
    WHB outlet 355 °C -> pipe run -> 390 °C -> tail-gas preheater III -> 280 °C
                      -> pipe run -> 310 °C -> economizer -> 200 °C
Each run is an adiabatic plug-flow section followed in time (Lagrangian gas parcel):
    d(xi)/dt = kp * n_tot * (pNO^2 pO2 - pNO2^2/K_ox) / P                 (xi = kmol O2 reacted)
    dT/dt    = -dH_ox(T) * d(xi)/dt / sum(n_i cp_i)                        (adiabatic)
An event stops each run at the outlet temperature of the flow sheet; the time it takes is
the residence time the flow sheet implies. Heat exchangers are treated as instantaneous
coolers (oxidation inside them is lumped into the runs). N2O4 and N2O3 are negligible above
~170 °C at 4.5 bar and are left out.

Run:  python hno3_train.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

RG = 8.314                                 # J/(mol K)
ATM = 1.01325                              # bar per atm
SP = ["NO", "O2", "NO2", "H2O", "N2"]
NU = np.array([-2.0, -1.0, 2.0, 0.0, 0.0])  # 2 NO + O2 -> 2 NO2

# NASA 7-coefficient polynomials (GRI-Mech 3.0), 200-1000 K
NASA = np.array([
    [4.21847630e+00, -4.63897600e-03, 1.10410220e-05, -9.33613540e-09, 2.80357700e-12, 9.84462300e+03, 2.28084640e+00],   # NO
    [3.78245636e+00, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12, -1.06394356e+03, 3.65767573e+00],  # O2
    [3.94403120e+00, -1.58542900e-03, 1.66578120e-05, -2.04754260e-08, 7.83505640e-12, 2.89661790e+03, 6.31199170e+00],   # NO2
    [4.19864056e+00, -2.03643410e-03, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12, -3.02937267e+04, -8.49032208e-01],  # H2O
    [3.29867700e+00, 1.40824040e-03, -3.96322200e-06, 5.64151500e-09, -2.44485400e-12, -1.02089990e+03, 3.95037200e+00],   # N2
])


def cp(T):
    a = NASA
    return RG*(a[:, 0] + a[:, 1]*T + a[:, 2]*T**2 + a[:, 3]*T**3 + a[:, 4]*T**4)


def h(T):
    a = NASA
    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 g_RT(T):
    a = NASA
    return (a[:, 0]*(1 - np.log(T)) - a[:, 1]*T/2 - a[:, 2]*T**2/6 - a[:, 3]*T**3/12 - a[:, 4]*T**4/20
            + a[:, 5]/T - a[:, 6])


def K_ox(T):
    """pNO2^2/(pNO^2 pO2), 1/atm."""
    return np.exp(-(NU @ g_RT(T))) * ATM


def kp(T):
    return 10.0**(652.1/T - 1.0366)          # atm^-2 s^-1


def burner_gas(y_nh3=0.10, S=0.965):
    """Gas after the gauze per kmol of NH3-air mixture (see hno3_burner.py)."""
    air = 1.0 - y_nh3
    NO = y_nh3*S
    O2 = 0.2095*air - y_nh3*(1.25*S + 0.75*(1 - S))
    N2 = 0.7905*air + 0.5*y_nh3*(1 - S)
    return np.array([NO, O2, 0.0, 1.5*y_nh3, N2])


def rate(n, T, P):
    """d(xi)/dt for a parcel with moles n at T (K) and P (atm)."""
    nt = n.sum()
    p = n/nt*P
    return kp(T)*nt*(p[0]**2*p[1] - p[2]**2/K_ox(T))/P


def run(n_in, T_in, P, T_out=None, t_max=60.0):
    """Adiabatic run from (n_in, T_in); stops when T reaches T_out (or at t_max)."""
    def rhs(t, y):
        n = n_in + NU*y[0]
        r = rate(n, y[1], P)
        return [r, -(NU @ h(y[1]))*r/(n @ cp(y[1]))]
    ev = (lambda t, y: y[1] - T_out) if T_out else None
    if ev:
        ev.terminal = True
    sol = solve_ivp(rhs, (0, t_max), [0.0, T_in], method="LSODA", rtol=1e-9, atol=1e-12,
                    events=ev, dense_output=True, max_step=0.05)
    return sol


def alpha_eq(T, n_in, P):
    def f(al):
        n = n_in + NU*al*n_in[0]/2
        p = n/n.sum()*P
        return p[2]**2 - K_ox(T)*p[0]**2*p[1]
    return brentq(f, 1e-12, 1 - 1e-12)


def duty(n, T1, T2):
    """Heat removed when cooling n from T1 to T2, J per kmol-basis mol."""
    return n @ (h(T1) - h(T2))


if __name__ == "__main__":
    P = 4.5/ATM                              # atm, after the burner and the waste-heat boiler
    n0 = burner_gas()                        # per mol NH3-air mixture
    N_mix = 6826/3600.0                      # kmol/s of mixture: 1000 t/d HNO3 at 279 kg NH3/t, 10 % NH3
    steps = [("run 1", 355, 390), ("tail-gas preheater III", 390, 280), ("run 2", 280, 310), ("economizer", 310, 200)]
    n, T, t_tot = n0.copy(), 355 + 273.15, 0.0
    tr, Tr, ar, path = [], [], [], []
    print(f"gas after burner: " + ", ".join(f"{s} {100*x:.2f} %" for s, x in zip(SP, n0/n0.sum())))
    for name, Ti, To in steps:
        if name.startswith("run"):
            sol = run(n, T, P, To + 273.15)
            tau = sol.t[-1]
            ts = np.linspace(0, tau, 100)
            xi = sol.sol(ts)[0]; Ts = sol.sol(ts)[1]
            tr += list(t_tot + ts); Tr += list(Ts - 273.15); ar += list(100*(n[2] + 2*xi)/n0[0])
            path.append((np.array(ar[-100:]), Ts - 273.15, "-"))
            a_in = 100*n[2]/n0[0]
            n = n + NU*xi[-1]; T = Ts[-1]; t_tot += tau
            V = N_mix*n.sum()*1e3*RG*T/(P*101325)
            print(f"{name}: {Ti} -> {To} °C in {tau:.2f} s; oxidation {a_in:.1f} -> {100*n[2]/n0[0]:.1f} %; "
                  f"{V:.1f} m3/s, {20*tau:.0f} m of pipe at 20 m/s")
        else:
            Q = duty(n, T, To + 273.15)*N_mix*1e3/1e6
            a = 100*n[2]/n0[0]
            path.append((np.array([a, a]), np.array([T - 273.15, To]), ":"))
            T = To + 273.15
            tr.append(t_tot); Tr.append(To); ar.append(a)
            print(f"{name}: {Ti} -> {To} °C, duty {Q:.1f} MW at {a:.1f} % oxidation")
    # a single long adiabatic run from 355 °C for comparison
    long = run(n0, 355 + 273.15, P, None, 60.0)
    tl = np.linspace(0, 60, 400); xl, Tl = long.sol(tl)
    print(f"one long adiabatic run: after 60 s {100*2*xl[-1]/n0[0]:.1f} % at {Tl[-1]-273.15:.0f} °C "
          f"(equilibrium at that T: {100*alpha_eq(Tl[-1], n0, P):.1f} %)")
    dT = -(NU @ h(630))*n0[0]/2/(n0 @ cp(630))
    print(f"adiabatic temperature rise for complete oxidation: {dT:.0f} K")

    fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
    Te = np.linspace(150, 600, 200)
    ax[0].plot(Te, [100*alpha_eq(t + 273.15, n0, P) for t in Te], "k--", lw=1.5, label="equilibrium, 4.5 bar")
    for i, (a, Tt, ls) in enumerate(path):
        ax[0].plot(Tt, a, ls, color="tab:red" if ls == "-" else "tab:blue", lw=2, zorder=3,
                   label=("adiabatic oxidation runs" if i == 0 else "gas cooling" if i == 1 else None))
    ax[0].plot(Tl - 273.15, 100*2*xl/n0[0], color="0.6", lw=1, zorder=2, label="one long adiabatic run")
    ax[0].set(xlabel="temperature (°C)", ylabel="degree of oxidation (%)", xlim=(150, 600), ylim=(0, 100),
              title="Nitrous-gas train, temperature–oxidation path")
    ax[0].legend(fontsize=8); ax[0].grid(alpha=0.3)
    ax[1].plot(tr, Tr, color="tab:red", lw=2); ax[1].set_ylabel("temperature (°C)", color="tab:red")
    b = ax[1].twinx(); b.plot(tr, ar, color="tab:green", lw=2); b.set_ylabel("degree of oxidation (%)", color="tab:green")
    ax[1].set(xlabel="cumulative residence time in the runs (s)", title="Temperature and oxidation along the train")
    ax[1].grid(alpha=0.3)
    fig.tight_layout(); fig.savefig("hno3_train.png", dpi=150)
    plt.show()