function hno3_oxidation()
% Nitric acid - isothermal gas-phase oxidation 2 NO + O2 -> 2 NO2
% bdrworks.org - Handbook - Nitric acid, worked model 2
%   d(xi)/dt = kp*n_tot*(pNO^2*pO2 - pNO2^2/K_ox)/P ,  log10 kp = 652.1/T - 1.0366  (atm, s)
%   2 NO2 <=> N2O4 and NO + NO2 <=> N2O3 are fast equilibria
% Save as hno3_oxidation.m and run:  hno3_oxidation      (MATLAB or GNU Octave)

n0 = nitrous_gas(0.10, 0.965, 0.14);
fprintf('dry nitrous gas: NO %.2f %%, O2 %.2f %%, N2 %.2f %%\n', 100*n0);
figure('Position',[60 100 1300 380]);

subplot(1,3,1);
t = logspace(-3, 2, 200);
for Tc = [30 50 100 200 350]
    [ts, xi] = oxidise(Tc + 273.15, 11, n0, 100);
    semilogx(t, 100*2*interp1(ts, xi, t)/n0(1), 'LineWidth', 2); hold on;
end
semilogx([1e-3 1e2], [90 90], 'k:');
xlabel('residence time (s)'); ylabel('degree of oxidation (%)'); title('Isothermal oxidation at 11 bar');
legend('30 \circC','50 \circC','100 \circC','200 \circC','350 \circC','Location','northwest'); grid on;

subplot(1,3,2);
Ts = 20:20:400;
for P = [1 4.5 11]
    t90 = arrayfun(@(Tc) time90(Tc + 273.15, P, n0), Ts);
    semilogy(Ts, t90, 'LineWidth', 2); hold on;
end
xlabel('temperature (\circC)'); ylabel('time to 90 % oxidation (s)'); title('Negative temperature coefficient');
legend('1 bar','4.5 bar','11 bar'); grid on;

subplot(1,3,3);
Te = 200:10:900;
for P = [1 4.5 11]
    plot(Te, 100*arrayfun(@(Tc) alpha_eq(Tc + 273.15, P, n0), Te), 'LineWidth', 2); hold on;
end
xlabel('temperature (\circC)'); ylabel('equilibrium degree of oxidation (%)'); title('Equilibrium limit');
legend('1 bar','4.5 bar','11 bar'); grid on;

fprintf('K_ox(298 K) = %.3e, K_ox(600 K) = %.3e 1/atm;  kp(298 K) = %.2f atm^-2 s^-1\n', K_ox(298.15), K_ox(600), kp(298.15));
for P = [1 4.5 11]
    fprintf('%5.1f bar: t90 = %.3g, %.3g, %.3g, %.3g s  (30, 50, 100, 200 C)\n', P, ...
            time90(303.15, P, n0), time90(323.15, P, n0), time90(373.15, P, n0), time90(473.15, P, n0));
end
for P = [1 4.5 11]
    fprintf('%5.1f bar, 30 C: closed-form t90 = %.3g s (constant moles, no dimers)\n', P, t_analytic(303.15, P, n0, 0.9));
end
for P = [1 4.5 11]
    fprintf('%5.1f bar: equilibrium oxidation 400 C %.1f %%, 600 C %.1f %%, 890 C %.2f %%\n', P, ...
            100*alpha_eq(673.15, P, n0), 100*alpha_eq(873.15, P, n0), 100*alpha_eq(1163.15, P, n0));
end
[a, b, c, d] = speciate(0.1*n0(1), 0.9*n0(1), n0(2) - 0.45*n0(1) + n0(3), 323.15, 11/1.01325);
fprintf('11 bar, 50 C, 90 %% oxidised: NO %.4f, NO2 %.4f, N2O4 %.4f, N2O3 %.4f; N2O4 carries %.0f %% of the NO2\n', ...
        a, b, c, d, 100*2*c/(0.9*n0(1)));
end

% ---------------------------------------------------------------- model
function n = nitrous_gas(y, S, sec)
% dry gas after the burner plus secondary air (share sec of the total air): NO, O2, N2
prim = 1 - y;  s = prim*sec/(1 - sec);
n = [y*S; 0.2095*(prim + s) - y*(1.25*S + 0.75*(1 - S)); 0.7905*(prim + s) + 0.5*y*(1 - S)];
n = n/sum(n);
end

function [ts, xi] = oxidise(T, Pbar, n0, tend)
P = Pbar/1.01325;
opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-12);
[ts, xi] = ode45(@(t, x) rhs(x, T, P, n0), [0 tend], 0, opts);
end

function t90 = time90(T, Pbar, n0)
% time to 90 % oxidation, NaN if the equilibrium lies below 90 %
if alpha_eq(T, Pbar, n0) <= 0.9, t90 = NaN; return; end
[ts, xi] = oxidise(T, Pbar, n0, 1e4);
al = 2*xi/n0(1);
k = find(al >= 0.9, 1);
t90 = interp1(al(k-1:k), ts(k-1:k), 0.9);
end

function dx = rhs(xi, T, P, n0)
NOt = max(n0(1) - 2*xi, 0);  O2 = n0(2) - xi;  NO2t = 2*xi;
[a, b, ~, ~, nt] = speciate(NOt, NO2t, n0(3) + O2, T, P);
pNO = a/nt*P;  pO2 = O2/nt*P;  pNO2 = b/nt*P;
dx = kp(T)*nt*(pNO^2*pO2 - pNO2^2/K_ox(T))/P;
end

function [a, b, c, d, nt] = speciate(NOt, NO2t, nrest, T, P)
% free NO (a), NO2 (b), N2O4 (c), N2O3 (d) from the totals; P in atm
% outer fixed point on the total moles, inner Newton iteration on the free NO2
k4 = K4(T)*P;  k3 = K3(T)*P;  nt = nrest + NOt + NO2t;  b = NO2t;
for it = 1:50
    for jt = 1:60
        f  = b + 2*k4*b^2/nt + k3*NOt*b/(nt + k3*b) - NO2t;
        df = 1 + 4*k4*b/nt + k3*NOt*nt/(nt + k3*b)^2;
        bn = min(max(b - f/df, 0), NO2t);
        if abs(bn - b) <= 1e-15*max(NO2t, 1e-30), b = bn; break; end
        b = bn;
    end
    a = NOt/(1 + k3*b/nt);  c = k4*b^2/nt;  d = k3*a*b/nt;
    new = nrest + a + b + c + d;
    if abs(new - nt) < 1e-13, nt = new; break; end
    nt = new;
end
end

function t = t_analytic(T, Pbar, n0, x)
% closed-form isothermal time for -dpNO/dt = 2 kp pNO^2 pO2 at constant total moles
P = Pbar/1.01325;  a = n0(1)*P;  b = n0(2)*P;  k = 2*kp(T);  m = 2*b - a;
t = 2/(k*m)*(x/(a*(1 - x)) - log((2*b - a*x)/(2*b*(1 - x)))/m);
end

function al = alpha_eq(T, Pbar, n0)
P = Pbar/1.01325;
f = @(al) resid(al, T, P, n0);
al = fzero(f, [1e-12 1 - 1e-12]);
end

function r = resid(al, T, P, n0)
xi = al*n0(1)/2;
n = [n0(1) - 2*xi; n0(2) - xi; 2*xi];  nt = sum(n) + n0(3);
p = n/nt*P;
r = p(3)^2 - K_ox(T)*p(1)^2*p(2);
end

function k = kp(T),  k = 10^(652.1/T - 1.0366);  end        % atm^-2 s^-1
function K = K4(T),  K = 0.698e-9*exp(6866/T);   end        % N2O4/NO2^2, 1/atm
function K = K3(T),  K = 65.3e-9*exp(4740/T);    end        % N2O3/(NO NO2), 1/atm

function K = K_ox(T)
% 2 NO + O2 = 2 NO2 from NASA polynomials (GRI-Mech 3.0): NO, O2, NO2; 1/atm
lo = [4.21847630e+00 -4.63897600e-03 1.10410220e-05 -9.33613540e-09 2.80357700e-12 9.84462300e+03 2.28084640e+00
      3.78245636e+00 -2.99673416e-03 9.84730201e-06 -9.68129509e-09 3.24372837e-12 -1.06394356e+03 3.65767573e+00
      3.94403120e+00 -1.58542900e-03 1.66578120e-05 -2.04754260e-08 7.83505640e-12 2.89661790e+03 6.31199170e+00];
hi = [3.26060560e+00 1.19110430e-03 -4.29170480e-07 6.94576690e-11 -4.03360990e-15 9.92097460e+03 6.36930270e+00
      3.28253784e+00 1.48308754e-03 -7.57966669e-07 2.09470555e-10 -2.16717794e-14 -1.08845772e+03 5.45323129e+00
      4.88475420e+00 2.17239560e-03 -8.28069060e-07 1.57475100e-10 -1.05108950e-14 2.31649830e+03 -1.17416950e-01];
if T < 1000, a = lo; else, a = hi; end
g = a(:,1)*(1 - log(T)) - a(:,2)*T/2 - a(:,3)*T^2/6 - a(:,4)*T^3/12 - a(:,5)*T^4/20 + a(:,6)/T - a(:,7);
K = exp(-(2*g(3) - 2*g(1) - g(2)))*1.01325;
end