function meoh_loop()
% Methanol synthesis loop - reactor, separator, recycle and purge
% bdrworks.org - Handbook - Methanol, worked model 3
%   tear stream R = (1 - phi) * gas( reactor(M + R) ), solved by Wegstein iteration
% Save as meoh_loop.m and run:  meoh_loop

p.P = 75;  p.Tin = 498;  p.Tc = 523;  p.Nt = 5000;  p.D = 0.040;  p.L = 7;  p.rho = 1140;  p.U = 250;  p.At = pi*p.D^2/4;
My = [0.17 0.08 0.72 0 0 0.02 0.01]';  M = 2500*My;
gas = [1 1 1 0 0 1 1]';
phis = [0.02 0.03 0.04 0.05 0.07 0.10 0.14 0.20];
R = 3*2500*[0.05 0.04 0.70 0 0 0.15 0.06]';
fprintf('SN = %.2f\n', (My(3) - My(2))/(My(1) + My(2)));
for ph = phis
    g = @(R) (1 - ph)*gas.*reactor(M + R, p);
    R = wegstein(g, R);
    out = reactor(M + R, p);
    fprintf('purge %4.0f %%: recycle %.2f, inerts %.1f %%, %.0f t/d, carbon eff. %.1f %%\n', 100*ph, ...
        sum(R)/2500, 100*(R(6) + R(7))/sum(R), out(4)*32.04*24/1000, 100*out(4)/(M(1) + M(2)));
end
end

function x = wegstein(g, x)
% bounded Wegstein acceleration of the fixed-point iteration x = g(x)
gx = g(x);  x1 = gx;
for it = 1:200
    gx1 = g(x1);
    s = (gx1 - gx)./(x1 - x + 1e-12);  q = s./(s - 1);  q = min(max(q, -5), 0);
    x = x1;  gx = gx1;  x1 = q.*x + (1 - q).*gx;
    if norm(x1 - x) < 1e-6*norm(x), break; end
end
x = x1;
end

function out = reactor(Fin, p)
f = Fin*1000/3600/p.Nt;
opts = odeset('RelTol', 1e-7, 'AbsTol', 1e-10);
[~, u] = ode45(@(z,u) rhs(u, p), [0 p.L], [f; p.Tin], opts);
out = u(end,1:7)'*p.Nt*3.6;
end

function du = rhs(u, p)
F = u(1:7);  T = u(8);  nu = stoich();
r = rates(T, p.P*F/sum(F));
du = [p.rho*p.At*(nu'*r);
      (p.rho*p.At*(-dH_rxn(T))'*r - p.U*pi*p.D*(T - p.Tc))/(F'*cp_all(T))];
end

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

function c = cp_all(T)
% ideal-gas cp, J/(mol K) (Reid, Prausnitz & Poling)
C = cpdata();  c = C*[1; T; T^2; T^3];
end

function H = h_all(T)
% H_f(298) + integral of cp, J/mol
C = cpdata();  T0 = 298.15;
Hf = [-110.53; -393.51; 0; -200.94; -241.83; -74.87; 0]*1e3;
H = Hf + C*[T-T0; (T^2-T0^2)/2; (T^3-T0^3)/3; (T^4-T0^4)/4];
end

function C = cpdata()
C = [30.87 -1.285e-2  2.789e-5 -1.272e-8;    % CO
     19.80  7.344e-2 -5.602e-5  1.715e-8;    % CO2
     27.14  9.274e-3 -1.381e-5  7.645e-9;    % H2
     21.15  7.092e-2  2.587e-5 -2.852e-8;    % CH3OH
     32.24  1.924e-3  1.055e-5 -3.596e-9;    % H2O
     19.25  5.213e-2  1.197e-5 -1.132e-8;    % CH4
     31.15 -1.357e-2  2.680e-5 -1.168e-8];   % N2
end

function dH = dH_rxn(T)
dH = stoich()*h_all(T);                      % J/mol, reactions 1 and 2
end

function [K1, K3] = K_eq(T)
% Graaf et al.: K1 (bar^-2) for CO2 + 3H2 = CH3OH + H2O, K3 for CO + H2O = CO2 + H2
K1 = 10^(3066/T - 10.592);  K3 = 10^(2073/T - 2.029);
end

function r = rates(T, p)
% Vanden Bussche & Froment (1996), mol/(kg_cat s), p in bar [CO CO2 H2 MeOH H2O ...]
p = max(p, 1e-10);  RG = 8.314;
k = @(A, B) A*exp(B/(RG*T));
k1 = k(1.07, 36696); k2 = k(3453.38, 0); k3 = k(0.499, 17197); k4 = k(6.62e-11, 124119); k5 = k(1.22e10, -94765);
[K1, K3] = K_eq(T);
den = 1 + k2*p(5)/p(3) + k3*sqrt(p(3)) + k4*p(5);
r = [k1*p(2)*p(3)*(1 - p(5)*p(4)/(K1*p(3)^3*p(2)))/den^3;
     k5*p(2)*(1 - K3*p(5)*p(1)/(p(2)*p(3)))/den];
end