% LEP 12-5 : Parallel Reactions in a PFR with Heat Effects
%

Vspan = [0 1]; % Range for the independent variable i.e. volume of the reactor
y0 = [423; 0; 0; 100]; % Initial values for the dependent variables i.e. T,Fc,Fb,and Fa
Cto=0.1;
Ua=4000;
dHrx1=-20000;
dHrx2=-60000;
k1a0=10;
k2a0=0.09;
Ta=373;
E1=33256;
E2=74826;
To=423;
[v y]=ode45(@(v,y)ODEfun(v,y,Cto,Ua,dHrx1,dHrx2,k1a0,k2a0,Ta,E1,E2,To),Vspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,600,400])
nexttile
plot(v,y(:,1),'LineWidth',1.5);
leg=legend('T');
ylabel('T(K)');
xlabel('V(dm^3)');
axis([0 1 300 1100]);
leg.ItemTokenSize=[10,10];
title('Temperature profile')
grid on
nexttile
plot(v,y(:,4),v,y(:,3),v,y(:,2),'LineWidth',1.5);
leg=legend('F_A','F_B','F_C');
ylabel('F_i(mol/s)');
xlabel('V(dm^3)');
axis([0 1 0 100]);
leg.ItemTokenSize=[10,10];
title('Molar flow rate profile')
grid on
z=size(y);
for i=1:z(1,1)
k1a(i) = k1a0 * exp((E1/8.314) * (1 / 300 - (1 / y(i,1)))); 
k2a(i) = k2a0 * exp((E2/8.314) * (1 / 300 - (1 /y(i,1)))); 
Ft(i) = y(i,4) + y(i,3) + y(i,2); 
Ca(i) = Cto * y(i,4) ./ Ft(i) * To /y(i,1); 
r1a(i) = 0 - (k1a(i) * Ca(i)); 
r2a(i) = 0 - (k2a(i) * Ca(i).^ 2); 
Qg(i)=(-r1a(i)) * (-dHrx1) + (-r2a(i)) * (-dHrx2) ;
Qr(i)=Ua*(y(i,1)-Ta);
Soverall(i)= y(i,3)/(y(i,2)+ 0.00001);
end
nexttile
plot(v,Qg(1:z(1,1)),v,Qr(1:z(1,1)),'LineWidth',1.5)
leg=legend('Q_g','Q_r');
leg.ItemTokenSize=[10,10];
title('Heat profile')
ylabel('Q(J/m^3.s)');
xlabel('V(m^3)');
axis([0 1 0 9e7]);
grid on
nexttile
plot(v,Soverall(1:z(1,1)),'LineWidth',1.5)
leg=legend('S_{B/C}');
leg.ItemTokenSize=[10,10];
title('Overall Selectivity profile')
ylabel('Selectivity');
xlabel('V(m^3)');
axis([0 1 0 40]);
grid on

%

function dYfuncvecdV = ODEfun(V,Yfuncvec,Cto,Ua,dHrx1,dHrx2,k1a0,k2a0,Ta,E1,E2,To); 
T = Yfuncvec(1); 
Fc = Yfuncvec(2); 
Fb = Yfuncvec(3); 
Fa = Yfuncvec(4); 
% Explicit equations
k1a = k1a0 * exp((E1/8.314) * (1 / 300 - (1 / T))); 
k2a = k2a0 * exp((E2/8.314) * (1 / 300 - (1 / T))); 
Ft = Fa + Fb + Fc; 
Ca = Cto * Fa / Ft * To / T; 
Cb = Cto * Fb / Ft * To / T; 
Cc = Cto * Fc / Ft * To / T; 
r1a = 0 - (k1a * Ca); 
r2a = 0 - (k2a * Ca ^ 2); 
Qg=(-r1a) * (-dHrx1) + (-r2a) * (-dHrx2) ;
Qr=Ua*(T-Ta);
% Differential equations
dTdV = (Qg-Qr) / (90 * Fa + 90 * Fb + 180 * Fc); 
dFcdV = 0 - (r2a / 2); 
dFbdV = 0 - r1a; 
dFadV = r1a + r2a; 
dYfuncvecdV = [dTdV; dFcdV; dFbdV; dFadV]; 

end