% LEP 13-5 : Multiple Reactions in a Semibatch Reactor
%

tspan = [0 1.5]; % Range for the independent variable time
y0 = [1.; 0; 0; 290.]; % Initial values for the dependent variables i.e. Ca,Cb,Cc,and T
Cao = 4; 
vo = 240; 
dHrx1=-6500;
dHrx2=8000;
UA=35000;
Ta=298;
T0=305;
E1A=9500;
E2B=7000;
[t y]=ode45(@(t,y)ODEfun(t,y,Cao,vo,dHrx1,dHrx2,UA,Ta,T0,E1A,E2B),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,400])
nexttile
plot(t,y(:,1),t,y(:,2),t,y(:,3),'LineWidth',2)
xlabel('t(h)')
ylabel('C_i(mol/dm^3)')
leg=legend('C_A','C_B','C_C');
leg.ItemTokenSize=[10,10];
title('Concentration-time curve')
axis([0 1.5 0 10]);
grid on
nexttile
plot(t,y(:,4),'LineWidth',2)
xlabel('t(h)')
ylabel('T(K)')
title('Temperature-time curve')
axis([0 1.5 250 600]);
grid on
z=size(y);
for i=1:z(1,1)
V(i) = 100 + vo * t(i);
k1a(i) = 1.25 * exp((E1A / 1.987) * (1 / 320 - (1 / y(i,4)))); 
k2b(i) = 0.08 * exp((E2B / 1.987) * (1 / 290 - (1 / y(i,4))));
Qg(i)=((dHrx1) * (-k1a(i) * y(i,1)) + (dHrx2)*(- k2b(i) * y(i,2))) * V(i);
Qr(i)=UA*(y(i,4)-Ta)+Cao * vo * 30 * (y(i,4) - T0);
end
nexttile
plot(t,Qg(1:z(1,1)),t,Qr(1:z(1,1)),'LineWidth',2)
xlabel('t(hours)')
ylabel('Q(J/hr)')
leg=legend('Q_g','Q_r');
title('Heat-time trajectory')
leg.ItemTokenSize=[10,10];
axis([0 1.5 -2*10^7 8*10^7]);
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,Cao,vo,dHrx1,dHrx2,UA,Ta,T0,E1A,E2B); 
Ca = Yfuncvec(1); 
Cb = Yfuncvec(2); 
Cc = Yfuncvec(3); 
T = Yfuncvec(4); 
% Explicit equations
k1a = 1.25 * exp((E1A / 1.987) * (1 / 320 - (1 / T))); 
k2b = 0.08 * exp((E2B / 1.987) * (1 / 290 - (1 / T))); 
ra = 0 - (k1a * Ca); 
V = 100 + vo * t; 
rc = 3 * k2b * Cb; 
rb = k1a * Ca / 2 - (k2b * Cb); 
Qg=((dHrx1) * (-k1a * Ca) + (dHrx2)*(- k2b * Cb)) * V;
Qr=UA*(T-Ta)+Cao * vo * 30 * (T - T0);
% Differential equations
dCadt = ra + (Cao - Ca) * vo / V; 
dCbdt = rb - (Cb * vo / V); 
dCcdt = rc - (Cc * vo / V); 
dTdt = (Qg - Qr)/((Ca * 30 + Cb * 60 + Cc * 20) * V + 100 * 35);
dYfuncvecdt = [dCadt; dCbdt; dCcdt; dTdt]; 
end