% LEP 8-8 : Membrane Reactor to improve Selectivity in multiple reactions
%

Vspan = [0 50]; % Range for the independent variable i.e. volume of the reactor 
y0 = [4; 0; 0; 0]; % Initial values for the dependent variables i.e. Fa,Fb,Fd,and Fu
Vspan1 = [0 50]; % Range for the independent variable i.e. volume of the reactor 
y01 = [4; 4; 0; 0]; % Initial values for the dependent variables i.e. Fa,Fb,Fd,and Fu
Ct0 = 0.8; 
k1a = 2; 
k2a = 3; 
Vt = 50; 
Fbo = 4;
[v,y]=ode45(@(v,y)ODEfun1(v,y,Ct0,k1a,k2a,Vt,Fbo),Vspan,y0);
[v1,y1]=ode45(@(v1,y1)ODEfun2(v1,y1,Ct0,k1a,k2a,Vt,Fbo),Vspan1,y01);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,800,400])
nexttile
plot(v,y(:,1),v,y(:,2),v,y(:,3),v,y(:,4),'LineWidth',1.5)
xlabel('V(dm^3)')
ylabel('F_i(mol/s)')
leg=legend('F_A','F_B','F_D','F_U');
leg.ItemTokenSize=[10,10];
title('Molar flow rates in a Membrane Reactor')
axis([0 50 0 4])
grid on
nexttile
plot(v1,y1(:,1),v1,y1(:,2),v1,y1(:,3),v1,y1(:,4),'LineWidth',1.5)
xlabel('V(dm^3)')
ylabel('F_i(mol/s)')
leg=legend('F_A','F_B','F_D','F_U');
leg.ItemTokenSize=[10,10];
title('Molar flow rates in a Plug Flow Reactor')
axis([0 50 0 4])
grid on
z=size(y);
Fao=4;
% Calculating values of selectivity
for i=1:z(1,1)
    if v(i,1)>0.0001
    Scd(i)=y(i,3)/y(i,4);
    else 
        Scd(i)=0;
    end
    X(i)=1-y(i,1)/Fao;
end
z1=size(y1);
Fao=4;
% Calculating values of selectivity
for i=1:z1(1,1)
    if v1(i,1)>0.0001
    Scd1(i)=y1(i,3)/y1(i,4);
    else 
        Scd(i)=0;
    end
    X1(i)=1-y1(i,1)/Fao;
end
nexttile
plot(v,Scd(1:z(1,1)),v1,Scd1(1:z1(1,1)),'LineWidth',1.5)
axis([0 50 0 60])
xlabel('V(dm^3)')
ylabel('S_D_/_U')
leg=legend('MR','PFR');
leg.ItemTokenSize=[10,10];
title('Selectivity as a function of volume')
grid on
nexttile
plot(v,X(1:z(1,1)),v1,X1(1:z1(1,1)),'LineWidth',1.5)
axis([0 50 0 1])
set(gca,'ytick',linspace(0,1,6))
xlabel('V(dm^3)')
ylabel('Conversion')
leg=legend('MR','PFR');
leg.ItemTokenSize=[10,10];
title('Conversion as a function of volume')
grid on

%

function dYfuncvecdV = ODEfun1(V,Yfuncvec,Ct0,k1a,k2a,Vt,Fbo); 
Fa = Yfuncvec(1); 
Fb = Yfuncvec(2); 
Fd = Yfuncvec(3); 
Fu = Yfuncvec(4); 
% Explicit equations
Ft = Fa + Fb + Fd + Fu; 
Cb = Ct0 * Fb / Ft; 
Ca = Ct0 * Fa / Ft; 
ra = 0 - (k1a * Ca ^ 2 * Cb) - (k2a * Ca * Cb ^ 2); 
rb = ra; 
Cd = Ct0 * Fd / Ft; 
Cu = Ct0 * Fu / Ft; 
rd = k1a * Ca ^ 2 * Cb; 
ru = k2a * Ca * Cb ^ 2; 
Rb = Fbo / Vt; 
Sdu = Fd / (Fu + 0.0000000000001); 
% Differential equations
dFadV = ra; 
dFbdV = rb + Rb; 
dFddV = rd; 
dFudV = ru; 
dYfuncvecdV = [dFadV; dFbdV; dFddV; dFudV];  
end

function dYfuncvecdV = ODEfun2(V,Yfuncvec,Ct0,k1a,k2a,Vt,Fbo); 
Fa = Yfuncvec(1); 
Fb = Yfuncvec(2); 
Fd = Yfuncvec(3); 
Fu = Yfuncvec(4); 
% Explicit equations
Ft = Fa + Fb + Fd + Fu; 
Cb = Ct0 * Fb / Ft; 
Ca = Ct0 * Fa / Ft; 
ra = 0 - (k1a * Ca ^ 2 * Cb) - (k2a * Ca * Cb ^ 2); 
rb = ra; 
Cd = Ct0 * Fd / Ft; 
Cu = Ct0 * Fu / Ft; 
rd = k1a * Ca ^ 2 * Cb; 
ru = k2a * Ca * Cb ^ 2; 
Sdu = Fd / (Fu + 0.0000000000001); 
% Differential equations
dFadV = ra; 
dFbdV = rb; 
dFddV = rd; 
dFudV = ru; 
dYfuncvecdV = [dFadV; dFbdV; dFddV; dFudV]; 
end