% ... Exercise so2so3-maxrate (4.1)
% ... Adiabatic microstructured PFR reactor for SO2 oxidation
% ... SO2 + 1/2 O2 --> SO3
% ... 
% ... 9.4.2013:
% ... Error found in the value of kr !!! 
% ... instead of 0.5154 we have to use 1.05e7!!!
% ...
% ... Milan Bernauer, ICT, Prague 2013

% Declarations of variables
global n kr EA R A1 A2 HR cp Feed ni
y0(1) = 0.10; % SO2
y0(2) = 0.11; % O2
y0(3) = 0.00; % SO3
y0(4) = 0.79; % N2
Feed  = 1.5e-6 ; % Feed rate [mol/s]
F0 = Feed .* y0; % Inlet feed rate of components [mol/s]
T0 = 873;     % Inlet temperature [K]
dr = 250e-6;  % reactor diameter [m]
Lr = 0.15;    % Reactor length [m]
Vr = pi .* dr.^2 ./ 4 .* Lr; % Reactor volume [m3]
n = 4;        % Number of compound          
kr = 1.05e7;  % Preexponential factor [mol3/m3/s]
EA = 60000;   % Activation energy
A1 = -11.12;  % Constant
A2 = 11750;   % Constant
cp = 30 ;     % Reaction mixture cp [J/mol/K]
HR = -98000;  % Reaction enthalpy [J/mol]
R = 8.314;    % Universal gas constant [J/mol/K]
ni = [ -1 -0.5 1 0]; % Stoichiometric matrix
% Initial conditions (feed ant T0)
u0      = F0;
u0(n+1) = T0;

% ... Adiabatic case
[z_adiab u] = ode15s(@dfdv,[0 Vr],u0);
% conversion
conv_adiab = 1 - u(:,1)/F0(1);
Tadb = u(:,n+1);
% ... Isothermal case
HR = 0;
[z_isoT u] = ode15s(@dfdv,[0 Vr],u0);
% conversion
conv_isoT = 1 - u(:,1)/F0(1);


subplot(2,1,1)
plot(z_adiab,conv_adiab,'-',z_isoT,conv_isoT,'-')
title('Conversion profile')
grid on
xlabel('V_R / m')
ylabel('X(SO_2)')
legend('Adiabatic','Isothermal','Location','Best')

subplot(2,1,2)
plot(z_adiab,Tadb,'-',z_isoT,u(:,n+1),'-')
title('Temperature profile')
xlabel('V_R / m')
ylabel('T / K')
legend('Adiabatic','Isothermal','Location','Best')
grid on