% LEP 13-1b : Adiabatic Batch Reactor (Part b)
%

tspan = [0 4000]; % Range for the independent variable
y0 = [286;0]; % Initial values for the dependent variables T and X
UA =10;
DeltaH =-20202; 
Ta1 = 290; 
Cpc = 4.18; 
mc =10;
NA0=54.8;
NB0=555;
NM=98.8;
Ea=18000.23;
R=1.987;
[t y]=ode45(@(t,y)ODEfun(t,y,UA,DeltaH,Ea,Ta1,Cpc,mc,NA0,NB0,NM,R),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,800,400])
nexttile
plot(t,y(:,2),'LineWidth',2);
leg=legend('X');
ylabel('Conversion');
xlabel('t(s)');
axis([0 4000 0 1])
title('Heat Exchange: Batch Reactor')
leg.ItemTokenSize=[10,10];
grid on
z=size(y);
thetaB=NB0/NA0;
thetaM=NM/NA0;
Cps=35+18*thetaB+19.5*thetaM;
% Calculating values of selectivity
for i=1:z(1,1)
     Ta2(i) = y(i,1)-(y(i,1)-Ta1)*exp(-UA/mc/Cpc);
     k(i) = .000273*exp((Ea/R)*(1/297-1/y(i,1)));
     Qr(i) = mc*Cpc*(y(i,1)-Ta1)*(1-exp(-UA/mc/Cpc));
   Qg(i) =NA0* k(i)*(1-y(i,2))*(-DeltaH); 
end
nexttile
plot(t,y(:,1),t,Ta2(1:z(1,1)),'LineWidth',2);
leg=legend('T','T_{a2}');
leg.ItemTokenSize=[10,10];
title('Heat Exchange: Batch Reactor')
ylabel('T(K)');
xlabel('t(s)');
axis([0 4000 280 360]);
grid on
nexttile
plot(t,k(1:z(1,1)),'LineWidth',2);
leg=legend('k');
leg.ItemTokenSize=[10,10];
title('Heat Exchange: Batch Reactor')
ylabel('k(s)');
xlabel('t(s)');
axis([0 4000 0 0.03]);
grid on
nexttile
plot(t,Qr(1:z(1,1)),t,Qg(1:z(1,1)),'LineWidth',2);
leg=legend('Q_r ','Q_g');
leg.ItemTokenSize=[10,10];
title('Heat Exchange: Batch Reactor')
ylabel('Q (cal/s)');
xlabel('k(s)');
axis([0 4000 0 2700]);
grid on

%

function dYfuncvecdV = ODEfun(V,Yfuncvec,UA,DeltaH,Ea,Ta1,Cpc,mc,NA0,NB0,NM,R); 
T = Yfuncvec(1);
X = Yfuncvec(2);
% Explicit equations
thetaB=NB0/NA0;
thetaM=NM/NA0;
Cps=35+18*thetaB+19.5*thetaM;
k = .000273*exp((Ea/R)*(1/297-1/T)); 
Ta2 = T-(T-Ta1)*exp(-UA/mc/Cpc);
Qr = mc*Cpc*(T-Ta1)*(1-exp(-UA/mc/Cpc)); 
Qg =NA0* k*(1-X)*(-DeltaH); 
DeltaQ = Qr-Qg;
% Differential equations
dTdt = (Qg-Qr)/Cps/NA0;
dXdt = k*(1-X); 
dYfuncvecdV = [dTdt;dXdt]; 
end