% LEP 12-1d : Butane Isomerization (Adiabatic)
%

Vspan = [0 5]; % Range for the independent variable i.e.volume of the reactor
y0 = [315; 305.; 0]; % Initial values for the dependent variables Ta,T, and x
Cpc=28;
mc=500;
Ua=0;
Cao=1.86;
deltaH=-34500;
FT0=163;
ya0=0.9;
CpI=161;
Ea=65700;
CpA=141;
[v y]=ode45(@(v,y)ODEfun(v,y,Cpc,mc,Ua,Cao,deltaH,FT0,ya0,CpI,Ea,CpA),Vspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,600,400])
nexttile
plot(v,y(:,2),'LineWidth',1.5);
leg=legend('T');
ylabel('T(K)');
xlabel('V(m^3)');
axis([0 5 300 420]);
leg.ItemTokenSize=[10,10];
title('Temperature profile')
grid on
z=size(y);
for i=1:z(1,1)
    k(i) = 31.1 * exp((Ea/8.314) * (y(i,2) - 360) / (y(i,2) * 360)); 
Kc(i) = 3.03 * exp((deltaH/8.314) * (y(i,2) - 333) / (y(i,2) * 333)); 
ra(i) = 0 - (k(i) * Cao * (1 - ((1 + 1 / Kc(i)) * y(i,3)))); 
Qg(i)= ra(i) * deltaH;
Qr(i)=(Ua * (y(i,2) - y(i,1)));
Xe(i)=Kc(i)/(1+Kc(i));
end
nexttile
plot(v,y(:,3),v,Xe(1:z(1,1)),'LineWidth',1.5)
leg=legend('X','X_e');
leg.ItemTokenSize=[10,10];
title('Conversion profile')
ylabel('Conversion');
xlabel('V(m^3)');
axis([0 5 0 1]);
grid on
nexttile
plot(v,-ra(1:z(1,1)),'LineWidth',1.5)
leg=legend('-r_A');
leg.ItemTokenSize=[10,10];
title('Rate profile')
ylabel('Rate (kmol/m^3.h)');
xlabel('V(m^3)');
axis([0 5 0 60]);
grid on
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(kJ/m^3.h)');
xlabel('V(m^3)');
axis([0 5 0 1500000]);
grid on

%

function dYfuncvecdV = ODEfun(V,Yfuncvec,Cpc,mc,Ua,Cao,deltaH,FT0,ya0,CpI,Ea,CpA); 
Ta = Yfuncvec(1); 
T = Yfuncvec(2); 
X = Yfuncvec(3);
% Explicit equations
k = 31.1 * exp((Ea/8.314) * (T - 360) / (T * 360)); 
Kc = 3.03 * exp((deltaH/8.314) * (T - 333) / (T * 333)); 
ra = 0 - (k * Cao * (1 - ((1 + 1 / Kc) * X))); 
Qg= ra * deltaH;
Qr=(Ua * (T - Ta));
thetaI=(1-ya0)/ya0;
sumcp=CpA+thetaI*CpI;
Fao=ya0*FT0/10;
% Differential equations
dTadV = Ua * (T - Ta) / (mc*Cpc); 
dTdV = (Qg -Qr ) / (sumcp * Fao); 
dXdV = 0 - (ra / Fao); 
dYfuncvecdV = [dTadV; dTdV; dXdV];  
end