function hno3_burner()
% Nitric acid - ammonia burner (Pt-Rh gauze), adiabatic energy balance
% bdrworks.org - Handbook - Nitric acid, worked model 1
%   4 NH3 + 5 O2 -> 4 NO + 6 H2O  (selectivity S),  4 NH3 + 3 O2 -> 2 N2 + 6 H2O  (rest)
%   H_in(T_mix) = H_out(T_gauze)
% Save as hno3_burner.m and run:  hno3_burner      (MATLAB or GNU Octave)

ys = linspace(0.06, 0.14, 41);
figure('Position',[100 100 1100 400]);
subplot(1,2,1); hold on;
cols = {'b','m','r'};  tm = [100 200 300];
for k = 1:3
    T = arrayfun(@(y) gauze_temperature(y, tm(k) + 273.15, 0.96), ys) - 273.15;
    plot(100*ys, T, cols{k}, 'LineWidth', 2);
    plot(100*ys, tm(k) + 68*100*ys, [cols{k} '--']);
end
plot([10 10.7], [890 900], 'ks', 'MarkerFaceColor', 'k');
plot([11 11], [500 1250], 'k:'); plot([13.5 13.5], [500 1250], 'k:');
xlabel('NH_3 in NH_3-air mixture (mol-%)'); ylabel('gauze temperature (\circC)');
title('Adiabatic gauze temperature, S = 96 %'); grid on;
legend('mix 100 \circC','rule t_{mix}+68(%NH_3)','mix 200 \circC','','mix 300 \circC','','plant data','Location','northwest');

subplot(1,2,2); hold on;
for S = [0.93 0.96 0.98]
    sur = arrayfun(@(y) o2_surplus(y, S), ys);
    plot(100*ys, sur, 'LineWidth', 2);
end
plot([6 14], [0 0], 'k');
xlabel('NH_3 in NH_3-air mixture (mol-%)'); ylabel('O_2 surplus after burner (mol per mol NH_3)');
title('Oxygen left for NO oxidation and absorption'); legend('S = 93 %','S = 96 %','S = 98 %'); grid on;

for t = tm
    [T, n1] = gauze_temperature(0.10, t + 273.15, 0.96);  y = n1/sum(n1);
    fprintf('10 %% NH3, mixture %d C: gauze %.0f C (rule %.0f C); NO %.2f %%, O2 %.2f %%, H2O %.2f %%\n', ...
            t, T - 273.15, t + 680, 100*y(4), 100*y(2), 100*y(5));
end
[T1] = gauze_temperature(0.10, 220 + 273.15, 0.965);
[T2] = gauze_temperature(0.107, 175 + 273.15, 0.96);
fprintf('plant check: 10 %% at 220 C -> %.0f C (plant ca. 890); 10.7 %% at 175 C -> %.0f C (plant ca. 900)\n', T1 - 273.15, T2 - 273.15);
for S = [0.93 0.96 0.98]
    fprintf('S = %.2f: O2 surplus after burner = %+.3f mol O2 per mol NH3\n', S, o2_surplus(0.10, S));
end
end

% ---------------------------------------------------------------- model
function n = feed(y)
% 1 mol NH3-air mixture: NH3, O2, N2, NO, H2O
air = 1 - y;
n = [y; 0.2095*air; 0.7905*air; 0; 0];
end

function n1 = burn(n0, S)
a = n0(1);
n1 = n0 + [-a; -a*(1.25*S + 0.75*(1 - S)); 0.5*a*(1 - S); a*S; 1.5*a];
end

function [T, n1] = gauze_temperature(y, Tmix, S)
n0 = feed(y);  n1 = burn(n0, S);
Hin = n0'*h_all(Tmix);
T = fzero(@(T) n1'*h_all(T) - Hin, [Tmix 2500]);
end

function s = o2_surplus(y, S)
n1 = burn(feed(y), S);
s = (n1(2) - 0.75*n1(4))/y;      % O2 left minus O2 needed for 4 NO + 3 O2 + 2 H2O -> 4 HNO3
end

% ---------------------------------------------------------------- NASA 7-coefficient polynomials (GRI-Mech 3.0)
function hh = h_all(T)
% molar enthalpies incl. heats of formation, J/mol: NH3, O2, N2, NO, H2O
lo = [4.28602740e+00 -4.66052300e-03 2.17185130e-05 -2.28088870e-08 8.26380460e-12 -6.74172850e+03 -6.25372770e-01
      3.78245636e+00 -2.99673416e-03 9.84730201e-06 -9.68129509e-09 3.24372837e-12 -1.06394356e+03 3.65767573e+00
      3.29867700e+00 1.40824040e-03 -3.96322200e-06 5.64151500e-09 -2.44485400e-12 -1.02089990e+03 3.95037200e+00
      4.21847630e+00 -4.63897600e-03 1.10410220e-05 -9.33613540e-09 2.80357700e-12 9.84462300e+03 2.28084640e+00
      4.19864056e+00 -2.03643410e-03 6.52040211e-06 -5.48797062e-09 1.77197817e-12 -3.02937267e+04 -8.49032208e-01];
hi = [2.63445210e+00 5.66625600e-03 -1.72786760e-06 2.38671610e-10 -1.25787860e-14 -6.54469580e+03 6.56629280e+00
      3.28253784e+00 1.48308754e-03 -7.57966669e-07 2.09470555e-10 -2.16717794e-14 -1.08845772e+03 5.45323129e+00
      2.92664000e+00 1.48797680e-03 -5.68476000e-07 1.00970380e-10 -6.75335100e-15 -9.22797700e+02 5.98052800e+00
      3.26060560e+00 1.19110430e-03 -4.29170480e-07 6.94576690e-11 -4.03360990e-15 9.92097460e+03 6.36930270e+00
      3.03399249e+00 2.17691804e-03 -1.64072518e-07 -9.70419870e-11 1.68200992e-14 -3.00042971e+04 4.96677010e+00];
if T < 1000, a = lo; else, a = hi; end
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