function smr_tube()
% Steam methane reforming - one fired reformer tube (Xu & Froment kinetics, Ergun pressure drop)
% bdrworks.org - Handbook - Steam reforming, worked model 2
% Save as smr_tube.m and run:  smr_tube

p = params();  q = 75e3;                                   % average heat flux, W/m2
[z, u] = tube(p, q);
[Twi, Two, hi] = wall_T(p, z, u(:,8), q);
F = u(end,1:7)';  X = 1 - F(1)/p.F0(1);  dry = sum(F) - F(2);
fprintf('h_i = %.0f W/m2K; outlet %.1f K, %.2f bar, CH4 conversion %.1f %%, slip %.2f %% dry, H2/CO %.2f\n', ...
    hi, u(end,8), u(end,9), 100*X, 100*F(1)/dry, F(5)/F(3));
fprintf('max outer wall temperature %.0f K\n', max(Two));
figure; plot(z, u(:,8), z, Twi, z, Two, 'LineWidth', 2); legend('gas', 'inner wall', 'outer wall'); xlabel('z (m)'); ylabel('T (K)'); grid on;
end

function p = params()
p.L = 12; p.Di = 0.1016; p.Do = 0.1320; p.A = pi*p.Di^2/4; p.rho = 1360; p.dp = 0.0173; p.eps = 0.5;
p.eta = [0.03 0.03 0.03]; p.mu = 3e-5; p.lam = 0.12; p.kw = 25; p.P0 = 29; p.T0 = 793;
p.MW = [16.043 18.015 28.010 44.009 2.016 28.013 31.999]';
p.F0 = [6; 18; 0; 0.1; 0.3; 0.1; 0]/3.6;                   % mol/s per tube, S/C = 3
end

function q = qflux(p, z, qavg)
q = qavg*(1.3 - 0.6*z/p.L);
end

function [z, u] = tube(p, qavg)
opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[z, u] = ode45(@(z,u) rhs(z, u, p, qavg), linspace(0, p.L, 300), [p.F0; p.T0; p.P0], opts);
end

function du = rhs(z, u, p, qavg)
F = u(1:7); T = u(8); P = u(9); y = F/sum(F); nu = stoich();
r = rates(T, P*y, p.eta)/3.6;                              % mol/(kg s)
Q = pi*p.Do*qflux(p, z, qavg);
G = (F'*p.MW)/1000/p.A;  rho = P*1e5*(y'*p.MW)/1000/(8.314*T);  us = G/rho;
dP = -(150*p.mu*(1-p.eps)^2/(p.eps^3*p.dp^2)*us + 1.75*rho*(1-p.eps)/(p.eps^3*p.dp)*us^2)/1e5;
du = [p.rho*p.A*(nu'*r); (Q + p.rho*p.A*(-dH_rxn(T))'*r)/(F'*cp_all(T)); dP];
end

function [Twi, Two, hi] = wall_T(p, z, T, qavg)
G = (p.F0'*p.MW)/1000/p.A;
hi = 0.813*p.lam/p.Di*(p.dp*G/p.mu)^0.9*exp(-6*p.dp/p.Di);   % Leva
q = qflux(p, z, qavg);
Twi = T + q*p.Do/(hi*p.Di);  Two = Twi + q*p.Do/2*log(p.Do/p.Di)/p.kw;
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