% LEP 8-5 : Multiple Gas Phase Reactions in a PBR
%

Wspan = [0 1000]; % Range for the weight of the catalyst 
y0 = [10; 10; 0; 0; 1]; % Initial values for the dependent variables i.e. Fa,Fb,Fc,Fd,and y
k1a = 100; 
k2c = 1500; 
Cto = 0.2; 
alpha = 0.0019; 
[w,y]=ode45(@(w,y)ODEfun(w,y,k1a,k2c,Cto,alpha),Wspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,400])
nexttile
plot(w,y(:,1),w,y(:,2),w,y(:,3),w,y(:,4),'LineWidth',2)
leg=legend('F_A','F_B','F_C','F_D');
leg.ItemTokenSize=[10,10];
set(gca,'ytick',linspace(0,10,6))
xlabel('W (kg)')
ylabel('F_i (mol/min)')
title('Molar Flow Rate Profiles')
grid on
Fao=10;
Fbo=10;
z=size(y);
for i=1:z(1,1)
    if w(i,1)>0.0001
    Scd(i)=y(i,3)/y(i,4);
    else 
        Scd(i)=0;
    end
   Xa(i)=(Fao-y(i,1))/Fao;
   Xb(i)=(Fbo-y(i,2))/Fbo; 
   Yc(i)=y(i,3)/(Fao-y(i,1));
   Yd(i)=y(i,4)/(Fao-y(i,1));
end
nexttile
plot(w,Scd(1:z(1,1)),'LineWidth',2)
xlabel('W (kg)')
ylabel('Selectivity')
leg=legend('S_C_/_D');
leg.ItemTokenSize=[10,10];
axis([0 1000 0 10000])
set(gca,'ytick',linspace(0,10000,6))
title('Selectivity Profile')
grid on
nexttile
plot(w,Xa(1:z(1,1)),w,Xb(1:z(1,1)),w,y(:,5),'LineWidth',2)
xlabel('W (kg)')
ylabel('X,p')
leg=legend('X_A','X_B','p');
leg.ItemTokenSize=[10,10];
axis([0 1000 0 1])
set(gca,'ytick',linspace(0,1,6))
title('Conversion/pressure drop Profile')
grid on
nexttile
plot(w,Yc(1:z(1,1)),w,Yd(1:z(1,1)),'LineWidth',2)
xlabel('W (kg)')
ylabel('Yield')
leg=legend('Y_C','Y_D');
set(gca,'ytick',linspace(0,1,6))
leg.ItemTokenSize=[10,10];
axis([0 1000 0 1])
title('Yield Profile')
grid on

%

function dYfuncvecdW = ODEfun(W,Yfuncvec,k1a,k2c,Cto,alpha); 
Fa = Yfuncvec(1); 
Fb = Yfuncvec(2); 
Fc = Yfuncvec(3); 
Fd = Yfuncvec(4); 
p = Yfuncvec(5); 
% Explicit equations
Ft = Fa + Fb + Fc + Fd; 
Ca = Cto * Fa / Ft * p; 
Cb = Cto * Fb / Ft * p; 
Cc = Cto * Fc / Ft * p; 
r1a = 0 - (k1a * Ca * Cb ^ 2); 
r1b = 2 * r1a; 
rb = r1b; 
r2c = 0 - (k2c * Ca ^ 2 * Cc ^ 3); 
r2a = 2 / 3 * r2c; 
r2d = -1 / 3 * r2c; 
r1c = 0 - r1a; 
rd = r2d; 
ra = r1a + r2a; 
rc = r1c + r2c; 
Cd = Cto * Fd / Ft * p; 
Fto = 20; 
if (W > 0.0001) 
    Scd = Fc / Fd;
else
    Scd = 0;
end 
% Differential equations
dFadW = ra; 
dFbdW = rb; 
dFcdW = rc; 
dFddW = rd; 
dpdW = 0 - (alpha / 2 / p * Ft / Fto); 
dYfuncvecdW = [dFadW; dFbdW; dFcdW; dFddW; dpdW]; 
end

%