function hno3_train()
% Nitric acid - NO oxidation in the nitrous-gas train of a dual-pressure plant
% bdrworks.org - Handbook - Nitric acid, worked model 3
%   adiabatic runs:  d(xi)/dt = kp*n_tot*(pNO^2*pO2 - pNO2^2/K_ox)/P
%                    dT/dt    = -dH_ox(T)*d(xi)/dt / sum(n_i cp_i)
%   each run stops at the outlet temperature of the flow sheet (event); coolers are instantaneous
% Save as hno3_train.m and run:  hno3_train      (MATLAB or GNU Octave)

P = 4.5/1.01325;                     % atm
nu = [-2; -1; 2; 0; 0];              % NO, O2, NO2, H2O, N2
n0 = burner_gas(0.10, 0.965);
Nmix = 6826/3600;                    % kmol/s of NH3-air mixture (1000 t/d HNO3)
fprintf('gas after burner: NO %.2f %%, O2 %.2f %%, H2O %.2f %%, N2 %.2f %%\n', 100*n0([1 2 4 5])/sum(n0));
steps = {'run 1', 355, 390; 'tail-gas preheater III', 390, 280; 'run 2', 280, 310; 'economizer', 310, 200};
n = n0;  T = 355 + 273.15;  ttot = 0;
figure('Position',[100 100 1100 400]);
subplot(1,2,1); hold on;
Te = linspace(150, 600, 120);
plot(Te, 100*arrayfun(@(t) alpha_eq(t + 273.15, n0, P, nu), Te), 'k--');
for k = 1:size(steps, 1)
    To = steps{k,3} + 273.15;
    if strncmp(steps{k,1}, 'run', 3)
        ev = @(t, u) event_T(u, To);
        opts = odeset('RelTol', 1e-9, 'AbsTol', 1e-12, 'Events', ev, 'MaxStep', 0.05);
        [ts, u] = ode45(@(t, u) rhs(u, n, P, nu), [0 60], [0; T], opts);
        a_in = 100*n(3)/n0(1);
        al = 100*(n(3) + 2*u(:,1))/n0(1);
        subplot(1,2,1); plot(u(:,2) - 273.15, al, 'r', 'LineWidth', 2);
        subplot(1,2,2); plot(ttot + ts, u(:,2) - 273.15, 'r', 'LineWidth', 2); hold on;
        n = n + nu*u(end,1);  T = u(end,2);  ttot = ttot + ts(end);
        fprintf('%s: %d -> %d C in %.2f s; oxidation %.1f -> %.1f %%\n', steps{k,1}, steps{k,2}, steps{k,3}, ts(end), a_in, 100*n(3)/n0(1));
    else
        Q = n'*(h_all(T) - h_all(To))*Nmix*1e3/1e6;
        a = 100*n(3)/n0(1);
        subplot(1,2,1); plot([T To] - 273.15, [a a], 'b:', 'LineWidth', 2);
        T = To;
        fprintf('%s: %d -> %d C, duty %.1f MW at %.1f %% oxidation\n', steps{k,1}, steps{k,2}, steps{k,3}, Q, a);
    end
end
opts = odeset('RelTol', 1e-9, 'AbsTol', 1e-12, 'MaxStep', 0.05);
[tl, ul] = ode45(@(t, u) rhs(u, n0, P, nu), [0 60], [0; 355 + 273.15], opts);
fprintf('one long adiabatic run: after 60 s %.1f %% at %.0f C\n', 100*2*ul(end,1)/n0(1), ul(end,2) - 273.15);
subplot(1,2,1); plot(ul(:,2) - 273.15, 100*2*ul(:,1)/n0(1), 'Color', [0.6 0.6 0.6]);
xlabel('temperature (\circC)'); ylabel('degree of oxidation (%)'); xlim([150 600]); ylim([0 100]); grid on;
title('Temperature-oxidation path');
subplot(1,2,2); xlabel('cumulative residence time in the runs (s)'); ylabel('temperature (\circC)'); grid on;
end

function [value, isterminal, direction] = event_T(u, To)
value = u(2) - To;  isterminal = 1;  direction = 0;
end

function du = rhs(u, nin, P, nu)
n = nin + nu*u(1);  T = u(2);
r = rate(n, T, P);
du = [r; -(nu'*h_all(T))*r/(n'*cp_all(T))];
end

function r = rate(n, T, P)
nt = sum(n);  p = n/nt*P;
r = kp(T)*nt*(p(1)^2*p(2) - p(3)^2/K_ox(T))/P;
end

function n = burner_gas(y, S)
air = 1 - y;
n = [y*S; 0.2095*air - y*(1.25*S + 0.75*(1 - S)); 0; 1.5*y; 0.7905*air + 0.5*y*(1 - S)];
end

function al = alpha_eq(T, nin, P, nu)
f = @(al) eqres(al, T, nin, P, nu);
al = fzero(f, [1e-12 1 - 1e-12]);
end

function r = eqres(al, T, nin, P, nu)
n = nin + nu*al*nin(1)/2;  p = n/sum(n)*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

% ---------------------------------------------------------------- NASA polynomials (GRI-Mech 3.0), 200-1000 K
function a = nasa()
a = [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
end

function c = cp_all(T)
a = nasa();  c = 8.314*(a(:,1) + a(:,2)*T + a(:,3)*T^2 + a(:,4)*T^3 + a(:,5)*T^4);
end

function hh = h_all(T)
a = nasa();  hh = 8.314*T*(a(:,1) + a(:,2)*T/2 + a(:,3)*T^2/3 + a(:,4)*T^3/4 + a(:,5)*T^4/5 + a(:,6)/T);
end

function K = K_ox(T)
a = nasa();
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