function smr_equilibrium()
% Steam methane reforming - chemical equilibrium (the thermodynamic limit)
% bdrworks.org - Handbook - Steam reforming, worked model 1
% Save as smr_equilibrium.m and run:  smr_equilibrium

Ts = linspace(800, 1300, 51);
figure('Position', [100 100 1000 380]); hold on;
for sc = [2 3 4]
    n0 = [1; sc; 0; 0; 0.05; 0; 0];  X = zeros(size(Ts));
    for k = 1:numel(Ts), n = equilibrium(Ts(k), 29, n0); X(k) = 1 - n(1)/n0(1); end
    plot(Ts, 100*X, 'LineWidth', 2);
    n = equilibrium(1123, 29, n0);  dry = sum(n) - n(2);
    fprintf('1123 K, 29 bar, S/C %d: CH4 conversion %.1f %%, slip %.2f %% dry, H2/CO %.2f\n', sc, 100*(1 - n(1)), 100*n(1)/dry, n(5)/n(3));
end
xlabel('T (K)'); ylabel('CH_4 conversion (%)'); legend('S/C = 2', 'S/C = 3', 'S/C = 4'); grid on;
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