function meoh_quench()
% Methanol synthesis - adiabatic multi-bed quench converter (ICI type)
% bdrworks.org - Handbook - Methanol, worked model 1
%   dF/dW = nu' r ,  (sum F cp) dT/dW = sum(-dH r) ;  cold-shot quench between beds
% Save as meoh_quench.m and run:  meoh_quench

p.P = 75;  p.F0 = 10000;  y0 = [0.10 0.05 0.70 0.005 0.001 0.12 0.024]';  p.y0 = y0/sum(y0);
p.beds = [8e3 10e3 14e3 18e3];  p.Tbed = 503;  p.Tq = 400;

f1 = fzero(@(f) run(f, p) - 1, [0.05 0.95]);           % fraction of feed to bed 1
[~, F, prof] = run(f1, p);
fprintf('feed to bed 1: %.1f %%\n', 100*f1);
figure; hold on;
for k = 1:numel(prof)
    u = prof{k}{2};
    fprintf('bed %d: T out %.1f K, MeOH out %.2f %%\n', k, u(end,8), 100*u(end,4)/sum(u(end,1:7)));
    plot(u(:,8), 100*u(:,4)./sum(u(:,1:7),2), 'LineWidth', 2);
end
fprintf('production %.0f kmol/h = %.0f t/d\n', F(4) - p.F0*p.y0(4), (F(4) - p.F0*p.y0(4))*32.04*24/1000);
xlabel('T (K)'); ylabel('CH_3OH (mol-%)'); grid on;
end

function [used, F, prof] = run(f1, p)
F = f1*p.F0*p.y0;  used = f1;  prof = {};  Tout = p.Tbed;
opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
for k = 1:numel(p.beds)
    if k > 1                                            % quench to Tbed (enthalpy balance)
        q = (F'*(h_all(Tout) - h_all(p.Tbed)))/(p.y0'*(h_all(p.Tbed) - h_all(p.Tq)));
        F = F + q*p.y0;  used = used + q/p.F0;
    end
    [W, u] = ode45(@(W,u) bed_rhs(u, p), [0 p.beds(k)], [F; p.Tbed], opts);
    prof{end+1} = {W, u};  %#ok<AGROW>
    F = u(end,1:7)';  Tout = u(end,8);
end
end

function du = bed_rhs(u, p)
F = u(1:7);  T = u(8);  nu = stoich();
r = rates(T, p.P*F/sum(F))*3.6;                         % kmol/(kg h)
du = [nu'*r; (-dH_rxn(T))'*r/(F'*cp_all(T))];           % kJ/kmol cancels J/mol
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