% LEP 12-7b : Complex Reactions with Heat Effects in a PFR (Counter-current)
%

Vspan = [0 10.]; % Range for the independent variable i.e length of the reactor
E1=8000;
E2=12000;
Cto=0.2;
Cpco=10;
m=50;
DH1b=-15000;
DH2a=-10000;
Ua=80;
Cpa=10;
Cpb=12;
Cpc=14;
Cpd=16;
Tao=507;
p=1;
[v y]=ode45(@(v,y)ODEfun(v,y,E1,E2,p,Cto,Cpco,m,DH1b,DH2a,Ua,Cpa,Cpb,Cpc,Cpd,Tao),Vspan,y0);
y0 = [5.; 10.; 0; 0; 300.; Tao]; % Initial values for the dependent variables i.e. Fa,Fb,Fc,Fd,T,and Ta
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,600,400])
nexttile
plot(v,y(:,5),v,y(:,6),'LineWidth',1.5);
leg=legend('T','T_a');
leg.ItemTokenSize=[10,10];
ylabel('T(K)');
xlabel('V(dm^3)');
axis([0 10 200 1400]);
title('Temperature profile')
grid on
nexttile
plot(v,y(:,1),v,y(:,2),v,y(:,3),v,y(:,4),'LineWidth',1.5);
leg=legend('F_A','F_B','F_C','F_D');
ylabel('F_i(mol/min)');
xlabel('V(dm^3)');
axis([0 10 0 10]);
leg.ItemTokenSize=[10,10];
title('Flow rate profile')
grid on
z=size(y);
R = 1.987;
To = 300; 
for i=1:z(1,1)
k1a(i) = 40 * exp(E1 / R * (1 / 300 - (1 / y(i,5)))); 
k2c(i) = 2 * exp(E2 / R * (1 / 300 - (1 / y(i,5)))); 
Ft(i) = y(i,1) + y(i,2) + y(i,3) + y(i,4); 
Ca(i) = Cto * y(i,1) / Ft(i) * To / y(i,5) * p; 
Cb(i) = Cto * y(i,2) / Ft(i) * To / y(i,5) * p; 
Cc(i) = Cto *  y(i,3) / Ft(i) * To /y(i,5)  * p; 
r1a(i) = 0 - (k1a(i) * Ca(i) * Cb(i).^ 2); 
r1b(i) = 2 * r1a(i);
r2c(i) = 0 - (k2c(i) * Ca(i).^ 2 * Cc(i).^ 3); 
r2a(i) = 2 / 3 * r2c(i); 
Qg(i) = r1b(i) * DH1b + r2a(i) * DH2a; 
Qr(i) = Ua * (y(i,5) - y(i,6)); 
Scd(i)=y(i,3)/(y(i,4)+0.00000001);
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(cal/m^3.min)');
xlabel('V(dm^3)');
axis([0 10 -1e5 4e5]);
grid on
nexttile
plot(v,Scd(1:z(1,1)),'LineWidth',1.5)
leg=legend('S_{C/D}');
leg.ItemTokenSize=[10,10];
title('Selectivity profile')
ylabel('Overall Selectivity');
xlabel('V(dm^3)');
axis([0 10 0 1000000]);
grid on

%

function dYfuncvecdV = ODEfun(V,Yfuncvec,E1,E2,p,Cto,Cpco,m,DH1b,DH2a,Ua,Cpa,Cpb,Cpc,Cpd,Tao); 
Fa = Yfuncvec(1); 
Fb = Yfuncvec(2); 
Fc = Yfuncvec(3); 
Fd = Yfuncvec(4); 
T = Yfuncvec(5); 
Ta = Yfuncvec(6); 
% Explicit equations
R = 1.987;
To = 300; 
Ft = Fa + Fb + Fc + Fd; 
k2c = 2 * exp(E2 / R * (1 / 300 - (1 / T))); 
Ca = Cto * Fa / Ft * To / T * p; 
Cb = Cto * Fb / Ft * To / T * p; 
Cc = Cto * Fc / Ft * To / T * p; 
r2c = 0 - (k2c * Ca ^ 2 * Cc ^ 3); 
k1a = 40 * exp(E1 / R * (1 / 300 - (1 / T))); 
r1a = 0 - (k1a * Ca * Cb ^ 2); 
r1b = 2 * r1a; 
rb = r1b; 
r2a = 2 / 3 * r2c; 
r1c = 0 - r1a; 
sumFiCpi = Cpa * Fa + Cpb * Fb + Cpc * Fc + Cpd * Fd; 
rc = r1c + r2c; 
r2d = -1 / 3 * r2c; 
ra = r1a + r2a; 
rd = r2d; 
Qg = r1b * DH1b + r2a * DH2a; 
Qr = Ua * (T - Ta); 
% Differential equations
dFadV = ra; 
dFbdV = rb; 
dFcdV = rc; 
dFddV = rd; 
dTdV = (Qg - Qr) / sumFiCpi; 
dTadV = -Ua * (T - Ta) / m / Cpco; 
dYfuncvecdV = [dFadV; dFbdV; dFcdV; dFddV; dTdV; dTadV]; 
end