% LEP 10-2 : Catalytic Reactor Design
%

Wspan = [0 10000]; % Range for the independent variable W
y0 = [0]; % Initial values for the dependent variables i.e X
FTo=50;
k=0.00087;
KT=1.038;
KB=1.39;
alpha = 0.000098;   
Po = 40; 
yH2=0.45;
[w,y]=ode45(@(w,y)ODEfun(w,y,FTo,k,KT,KB,alpha,Po,yH2),Wspan,y0);
tiledlayout(1,2)
set (gcf,'Position',[0,0.1,800,300])
z=size(w);
PTo = 0.3*Po;
thetaH2=yH2/0.3;
for i=1:z(1,1)
p(i) = (1-alpha*w(i))^0.5;
PH2(i)= PTo*(thetaH2-y(i))*p(i); 
PB(i)= PTo*y(i)*p(i); 
PT(i)= PTo*(1-y(i))*p(i);
end
nexttile
plot(w,y,w,p(1:z(1,1)),'LineWidth',2)
ylabel('X,p');
xlabel('W(kg)');
leg=legend('X','p');
leg.ItemTokenSize=[10,10];
axis([0 10000 0 1])
title('Conversion and pressure ratio profiles')
grid on
nexttile
plot(w,PH2(1:z(1,1)),w,PB(1:z(1,1)),w,PT(1:z(1,1)),'LineWidth',2)
ylabel('P_i(atm)');
xlabel('W(kg)');
leg=legend('P_{H2}','P_B','P_T');
leg.ItemTokenSize=[10,10];
axis([0 10000 0 20])
title('Partial pressure profiles')
grid on

%

function dYfuncvecdt = ODEfun(W,Y,FTo,k,KT,KB,alpha,Po,yH2); 
X = Y(1); 
% Explicit equations

PTo = 0.3*Po; 
p = (1-alpha*W)^0.5;
P = p*Po; 
thetaH2=yH2/0.3;
PH2 = PTo*(thetaH2-X)*p; 
PB = PTo*X*p; 
PT = PTo*(1-X)*p; 
rt = -k*PT*PH2/(1+KB*PB+KT*PT); 
RATE = -rt; 
% Differential equations
dXdW= -rt/FTo; 
dYfuncvecdt = [dXdW]; 

end
