%% --- Highly exothermic reaction in PFR
%% --- eqquiped with external cooling
%% --- Reaction: XY + 3 O2 --> PAH + 3 H2O
%%
%% --- Global variables (shared with routine/function
%%     dfun)
global kr Ac F0 Tr Ea Tm cpg kapa HR p0 rd L Q T0 S Rg Mw yi0

%% --- Variables declaration
kr = 2.0822;        % rate constant [1/s]
HR = -1.284e6;      % Heat of reaction [kJ/kmol]
L  = 1.5;           % reactor lenght [m]
Q  = 2.6371e-3;     % Mass flow [kg/s]
T0 = 615;           % Temperature [K]
p0 = 101.325;       % pressure [kPa]
rd = 0.0125;        % reactor radius [m]
S  = pi().*rd.^2;   % S [m2] 
Ac = rd .* pi();    % [m] 
Ea = 13636.0;       % activation eng. factor [K] 
Tm = 625;           % cool. medium temp. [K]
Tr = 625;           % ref. temperature [K]
cpg   = 0.992;      % spec. heat capacity [kJ/kg/K]
kapa  = 0.373;      % Overall heat trans. coef [kW/m2/K] 
Rg = 8.314;         % universal gas constant [kJ/kmol/K]
% Molar mass [kg/mol]
%      XY      O2         PAH  H2O N2
Mw  = [106.17, 32,        148, 18, 28]; 
yi0 = [0.02,   0.98*0.21, 0,   0,  0.98*0.79];

MM = Mw * yi0';     % mean molar mass [kg/kmol] 
F0 = Q ./ MM;       % Inlet molar flow [kmol/s]

%% -- CALL of ODE solver
[z X] = ode45(@dfun,[0 L], [0 T0]);

%% --- Display results (figs)
subplot(2,1,1)
plot(z,X(:,1),'-b', 'LineWidth',2)
grid on
ylim([0 1])
xlabel('z / m')
ylabel('XY Conversion')
legend('X')

subplot(2,1,2)
plot(z,X(:,2),'-r', 'LineWidth',2)
grid on
xlabel('z / m')
ylabel('Temperature / K')
legend('T(reactor)')

