function atr()
% Steam methane reforming - autothermal reformer (burner + adiabatic Ni bed)
% bdrworks.org - Handbook - Steam reforming, worked model 3
% Save as atr.m and run:  atr

p.P = 30; p.F = 1000; p.SC = 0.6; p.Tg = 900; p.To2 = 500; p.W = 20e3; p.eta = [0.03 0.03 0.03];
for o2c = [0.50 0.54 0.57 0.60 0.64]
    [n0, Tb] = burner(o2c, p);
    opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
    [~, u] = ode45(@(W,u) bed(u, p), [0 p.W], [n0; Tb], opts);
    F = u(end,1:7)';  dry = sum(F) - F(2);
    fprintf('O2/C %.2f: burner %.0f K, outlet %.0f K, CH4 slip %.2f %% dry, H2/CO %.2f\n', o2c, Tb, u(end,8), 100*F(1)/dry, F(5)/F(3));
end
end

function [n, T] = burner(o2c, p)
gas = [p.F; p.SC*p.F; 0; 0; 0.02*p.F; 0.01*p.F; 0];  o2 = [0; 0; 0; 0; 0; 0; o2c*p.F];
Hin = gas'*h_all(p.Tg) + o2'*h_all(p.To2);
n = gas + o2 + [-1; 2; 1; 0; 0; 0; -1.5]*(o2c*p.F/1.5);    % CH4 + 1.5 O2 -> CO + 2 H2O
T = fzero(@(T) n'*h_all(T) - Hin, [600 4000]);
end

function du = bed(u, p)
F = u(1:7);  T = u(8);  nu = stoich();
r = rates(T, p.P*F/sum(F), p.eta);                        % kmol/(kg h)
du = [nu'*r; (-dH_rxn(T))'*r/(F'*cp_all(T))];
end

% ================================================================ properties & kinetics
function nu = stoich()
%       CH4 H2O  CO CO2  H2 N2  O2
nu = [  -1  -1    1   0   3  0   0;      % CH4 + H2O   <=> CO + 3H2
         0  -1   -1   1   1  0   0;      % CO + H2O    <=> CO2 + H2
        -1  -2    0   1   4  0   0];     % CH4 + 2H2O  <=> CO2 + 4H2
end

function a = nasa(T)
% NASA 7-coefficient polynomials (GRI-Mech 3.0) for CH4 H2O CO CO2 H2 N2 O2
if T < 1000
  a = [
     5.14987613e+00, -1.36709788e-02, 4.91800599e-05, -4.84743026e-08, 1.66693956e-11, -1.02466476e+04, -4.64130376e+00;   % CH4
     4.19864056e+00, -2.03643410e-03, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12, -3.02937267e+04, -8.49032208e-01;   % H2O
     3.57953347e+00, -6.10353680e-04, 1.01681433e-06, 9.07005884e-10, -9.04424499e-13, -1.43440860e+04, 3.50840928e+00;   % CO
     2.35677352e+00, 8.98459677e-03, -7.12356269e-06, 2.45919022e-09, -1.43699548e-13, -4.83719697e+04, 9.90105222e+00;   % CO2
     2.34433112e+00, 7.98052075e-03, -1.94781510e-05, 2.01572094e-08, -7.37611761e-12, -9.17935173e+02, 6.83010238e-01;   % H2
     3.29867700e+00, 1.40824040e-03, -3.96322200e-06, 5.64151500e-09, -2.44485400e-12, -1.02089990e+03, 3.95037200e+00;   % N2
     3.78245636e+00, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12, -1.06394356e+03, 3.65767573e+00;   % O2
  ];
else
  a = [
     7.48514950e-02, 1.33909467e-02, -5.73285809e-06, 1.22292535e-09, -1.01815230e-13, -9.46834459e+03, 1.84373180e+01;   % CH4
     3.03399249e+00, 2.17691804e-03, -1.64072518e-07, -9.70419870e-11, 1.68200992e-14, -3.00042971e+04, 4.96677010e+00;   % H2O
     2.71518561e+00, 2.06252743e-03, -9.98825771e-07, 2.30053008e-10, -2.03647716e-14, -1.41518724e+04, 7.81868772e+00;   % CO
     3.85746029e+00, 4.41437026e-03, -2.21481404e-06, 5.23490188e-10, -4.72084164e-14, -4.87591660e+04, 2.27163806e+00;   % CO2
     3.33727920e+00, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14, -9.50158922e+02, -3.20502331e+00;   % H2
     2.92664000e+00, 1.48797680e-03, -5.68476000e-07, 1.00970380e-10, -6.75335100e-15, -9.22797700e+02, 5.98052800e+00;   % N2
     3.28253784e+00, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, -1.08845772e+03, 5.45323129e+00;   % O2
  ];
end
end

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

function H = h_all(T)
a = nasa(T);  H = 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);   % J/mol
end

function dH = dH_rxn(T)
dH = stoich()*h_all(T);
end

function [K1, K2, K3] = K_eq(T)
K1 = exp(-26830/T + 30.114);  K2 = exp(4400/T - 4.036);  K3 = K1*K2;      % Xu & Froment
end

function r = rates(T, p, eta)
% Xu & Froment (1989) rates times effectiveness factors, kmol/(kg_cat h); p in bar [CH4 H2O CO CO2 H2 ...]
p = max(p, 1e-10);  RT = 8.314e-3*T;
k1 = 4.225e15*exp(-240.1/RT);  k2 = 1.955e6*exp(-67.13/RT);  k3 = 1.020e15*exp(-243.9/RT);
KCO = 8.23e-5*exp(70.65/RT);  KH2 = 6.12e-9*exp(82.90/RT);  KCH4 = 6.65e-4*exp(38.28/RT);  KH2O = 1.77e5*exp(-88.68/RT);
[K1, K2, K3] = K_eq(T);
den = 1 + KCO*p(3) + KH2*p(5) + KCH4*p(1) + KH2O*p(2)/p(5);
r = eta(:).*[k1/p(5)^2.5*(p(1)*p(2) - p(5)^3*p(3)/K1)/den^2;
             k2/p(5)*(p(3)*p(2) - p(5)*p(4)/K2)/den^2;
             k3/p(5)^3.5*(p(1)*p(2)^2 - p(5)^4*p(4)/K3)/den^2];
end

function n = equilibrium(T, P, n0)
% equilibrium of reactions 1 and 2 by Newton's method on the log equations
nu = stoich();  [K1, K2] = K_eq(T);  x = [0.7*n0(1); 0.2*n0(1)];
f = @(x) eqres(x, n0, nu, P, K1, K2);
for it = 1:100
    F = f(x);  J = zeros(2);  h = 1e-7*max(abs(x), 1e-6);
    for j = 1:2, e = zeros(2,1); e(j) = h(j); J(:,j) = (f(x + e) - F)/h(j); end
    dx = -J\F;  lam = 1;
    while any(n0 + nu(1,:)'*(x(1)+lam*dx(1)) + nu(2,:)'*(x(2)+lam*dx(2)) < 0 & [1 1 1 1 1 0 0]') && lam > 1e-4, lam = lam/2; end
    x = x + lam*dx;
    if norm(dx) < 1e-12*max(1, norm(x)), break; end
end
n = n0 + nu(1,:)'*x(1) + nu(2,:)'*x(2);
end

function F = eqres(x, n0, nu, P, K1, K2)
n = max(n0 + nu(1,:)'*x(1) + nu(2,:)'*x(2), 1e-14);  p = P*n/sum(n);
F = [log(p(3)*p(5)^3/(p(1)*p(2))) - log(K1);  log(p(4)*p(5)/(p(3)*p(2))) - log(K2)];
end