% LEP 14-7: Flow, Diffusion and Reaction in a Packed bed
%

Wspan = [0 100]; % Range for the independent variable
y0 = [0]; % Initial values for the dependent variables 
U=1;
DAB=1e-8;
dp=0.001;
kr=0.5;
rho=1000;
mu=0.001;
Ca0=1;
phi=0.4;
ys=1;
Ac=0.05;
rhoc=3000   ;

[w y]=ode45(@(w,y)ODEfun(w,y,U,DAB,dp,kr,rho,mu,Ca0,phi,ys,Ac,rhoc),Wspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.5,1000,800])
nexttile
plot(w,y(:,1),'LineWidth',2);
set(gca,'ytick',linspace(0,1,5))
leg=legend('X');
ylabel('Conversion');
xlabel('t(s)');
axis([0 100 0 1])
title('Conversion Plot')
leg.ItemTokenSize=[10,10];
grid on
z=size(y);
nu=mu/rho;
Re=dp*U/nu;
Reprime=Re/((1-phi)*ys);
Sc=nu/DAB;
Shprime=((Reprime)^0.5)*((Sc)^(1/3));
kc=DAB*Shprime*(1-phi)*ys/(dp*phi);
k=kr*kc/(kc+kr);
% Calculating values of selectivity
for i=1:z(1,1)
     Ca(i)=Ca0*(1-y(i,1));
     ra(i)=-k*Ca(i);
     kc(i)=k;
     
end
nexttile
plot(w,-ra(1:z(1,1)),'LineWidth',2)
leg=legend('W_A');
leg.ItemTokenSize=[10,10];
title('Molar Flux')
ylabel('Molar flux (mol/m^2.s)');
xlabel('W (kg)');
set(gca,'ytick',linspace(0,5e-3,5))
axis([0 100 0 0.005]);
grid on
nexttile
plot(w,kc,'LineWidth',2)
leg=legend('k');
leg.ItemTokenSize=[10,10];
title('Rate constant')
set(gca,'ytick',linspace(0,0.01,5))
ylabel('Rate constant (m/s)');
xlabel('W (kg)');
axis([0 100 0 0.01]);
grid on

%

function dYfuncvecdV = ODEfun(W,Yfuncvec,U,DAB,dp,kr,rho,mu,Ca0,phi,ys,Ac,rhoc); 
X = Yfuncvec(1);
% Explicit equations
nu=mu/rho;
Re=dp*U/nu;
Reprime=Re/((1-phi)*ys);
Sc=nu/DAB;
Shprime=((Reprime)^0.5)*((Sc)^(1/3));
kc=DAB*Shprime*(1-phi)*ys/(dp*phi);
k=kr*kc/(kc+kr);
Fa0=Ca0*U*Ac;
Ca=Ca0*(1-X);
ra = -k * Ca;
acprime=6*(1-phi)/(dp*rhoc);
raprime=ra*acprime;

% Differential equations
dXdt = -raprime/Fa0; 
dYfuncvecdV = [dXdt];  
end