% LEP 13-4 : Heat Effects in a Semibatch Reactor
%

tspan = [0 360.]; % Range for the independent variable time
y0 = [5.; 0; 0; 300.; 6.14]; % Initial values for the dependent variables i.e. Ca,Cb,Cc,T,and Nw
v0=0.004;
Cb0=1;
UA=3000;
cp=75240;
T0=300;
dH=-79076000;
Cw0=55;
Vi=0.2;
cpa=170700;
Cpc=18;
Ta1=285;
mc=280;
[t y]=ode45(@(t,y)ODEfun(t,y,v0,Cb0,UA,cp,T0,dH,Cw0,Vi,cpa,Cpc,Ta1,mc),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(s)')
ylabel('C_i(kmol/m^3)')
leg=legend('C_A','C_B','C_C');
leg.ItemTokenSize=[10,10];
title('Concentration-time curve')
axis([0 360 0 5])
grid on
z=size(y);
Fw=Cw0*v0;
Fb0 = Cb0 * v0;
for i=1:z(1,1)
    % Calculating the values of Ta2 from the following expression
    Ta2(i) = y(i,4) - ((y(i,4) - Ta1) * exp(0 - (UA /mc /Cpc))); 
    Qrs2(i) = -mc * Cpc * (Ta1 - y(i,4)) * (1 - exp(0 - (UA / mc / Cpc))); 
Qrs1(i)=(Fb0+Fw)*cp*(y(i,4)-T0);
Qrs(i)=Qrs1(i)+Qrs2(i);
k(i) = 0.39175 * exp(5472.7 * (1 / 273 - (1 / y(i,4))));
Cd(i)=y(i,3);
Kc(i) = 10 ^ (3885.44 / y(i,4));
ra(i) = 0 - (k(i) * (y(i,1) * y(i,2) - (y(i,3) * Cd(i) / Kc(i)))); 
V(i) = Vi + v0 * t(i);
Qgs(i)=ra(i) * V(i) * dH;
end
nexttile
plot(t,y(:,4),t,Ta2(1:z(1,1)),'LineWidth',2)
xlabel('t(s)')
ylabel('T(K)')
axis([0 360 260 380])
leg=legend('T','T_{a2}');
title('Temperature-time curve')
leg.ItemTokenSize=[10,10];
grid on

nexttile
plot(t,Qgs(1:z(1,1)),t,Qrs(1:z(1,1)),'LineWidth',2)
xlabel('t(s)')
ylabel('Q(kJ/s)')
axis([0 360 0 800000])
leg=legend('Q_g','Q_{r}');
leg.ItemTokenSize=[10,10];
title('Heat-time trajectory')
grid on

nexttile
plot(t,-ra(1:z(1,1)),'LineWidth',2)
xlabel('t(s)')
ylabel('Rate(kmol/m^3.s)')
axis([0 360 0 0.05])
leg=legend('-r_A');
leg.ItemTokenSize=[10,10];
title('Rate-time trajectory')
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,v0,Cb0,UA,cp,T0,dH,Cw0,Vi,cpa,Cpc,Ta1,mc); 
Ca = Yfuncvec(1); 
Cb = Yfuncvec(2); 
Cc = Yfuncvec(3); 
T = Yfuncvec(4); 
Nw = Yfuncvec(5); 
% Explicit equations
k = 0.39175 * exp(5472.7 * (1 / 273 - (1 / T))); 
Cd = Cc; 
Kc = 10 ^ (3885.44 / T); 
V = Vi + v0 * t; 
Fb0 = Cb0 * v0; 
ra = 0 - (k * (Ca * Cb - (Cc * Cd / Kc))); 
Na = V * Ca; 
Nb = V * Cb; 
Nc = V * Cc; 
rb = ra; 
rc = 0 - ra; 
Nd = V * Cd; 
NCp = cp * (Nb + Nc + Nd + Nw) + cpa * Na; 
Fw=Cw0*v0;
Qrs2 = -mc * Cpc * (Ta1 - T) * (1 - exp(0 - (UA / mc / Cpc))); 
Qrs1=(Fb0+Fw)*cp*(T-T0);
Qrs=Qrs1+Qrs2;
Qgs=ra * V * dH;
% Differential equations
dCadt = ra - (v0 * Ca / V); 
dCbdt = rb + v0 * (Cb0 - Cb) / V; 
dCcdt = rc - (Cc * v0 / V); 
dTdt = (Qgs-Qrs) / NCp; 
dNwdt = v0 * Cw0; 
dYfuncvecdt = [dCadt; dCbdt; dCcdt; dTdt; dNwdt]; 
end