function meoh_tubular()
% Methanol synthesis - boiling-water-cooled multitubular reactor (Lurgi type)
% bdrworks.org - Handbook - Methanol, worked model 2
%   dF/dz = rho_b A_t nu' r ,  (sum F cp) dT/dz = rho_b A_t sum(-dH r) - U pi D (T - Tc)
% Save as meoh_tubular.m and run:  meoh_tubular

p = params();
Tc = 523;
[z, u] = tube(p, Tc);
[meoh, tpd, xc, y] = summary(p, u);
[Tmax, i] = max(u(:,8));
fprintf('hot spot %.1f K at z = %.2f m, outlet %.1f K\n', Tmax, z(i), u(end,8));
fprintf('MeOH out %.2f mol-%%, CO+CO2 conversion %.1f %%, production %.0f kmol/h = %.0f t/d\n', 100*y, 100*xc, meoh, tpd);

% coolant temperature sweep
tc = 470:2.5:550;  prod = zeros(size(tc));  hs = prod;
for k = 1:numel(tc)
    [~, uk] = tube(p, tc(k));  [~, prod(k)] = summary(p, uk);  hs(k) = max(uk(:,8));
end
[~, ib] = max(prod);
fprintf('coolant temperature for maximum production = %.0f K\n', tc(ib));
figure('Position', [100 100 1000 380]);
subplot(1,2,1); plot(z, u(:,8), 'LineWidth', 2); xlabel('z (m)'); ylabel('T (K)'); grid on;
subplot(1,2,2); plot(tc, prod, 'LineWidth', 2); xlabel('T_c (K)'); ylabel('methanol (t/d)'); grid on;
end

function p = params()
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.Tin = 498;  p.Nt = 5000;  p.D = 0.040;  p.L = 7;  p.rho = 1140;  p.U = 250;  p.At = pi*p.D^2/4;
end

function [z, u] = tube(p, Tc)
f0 = p.F0*1000/3600/p.Nt*p.y0;                         % mol/s per tube
opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[z, u] = ode45(@(z,u) rhs(u, p, Tc), linspace(0, p.L, 400), [f0; p.Tin], opts);
end

function du = rhs(u, p, Tc)
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 - Tc))/(F'*cp_all(T))];
end

function [meoh, tpd, xc, y] = summary(p, u)
Fi = u(1,1:7);  Fo = u(end,1:7);
meoh = (Fo(4) - Fi(4))*p.Nt*3.6;  tpd = meoh*32.04*24/1000;
xc = 1 - (Fo(1) + Fo(2))/(Fi(1) + Fi(2));  y = Fo(4)/sum(Fo);
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