% LEP 13-3 : Startup of a CSTR (English units)
%


UA=16000;
Ta1=60;
T0=75;
dH=-36000;
mc=1000;
Fa0=80;
Fb0=1000;
Fm0=100;
V=66.809;
Ea=32400;
CAi=0;
Ti=75;
tspan = [0 4.]; % Range for the independent variable time
y0 = [CAi; 3.45; 0; 0; Ti]; % Initial values for the dependent variables i.e. Ca,Cb,Cc,Cm,and T
[t y]=ode45(@(t,y)ODEfun(t,y,UA,Ta1,T0,dH,mc,Fa0,Fb0,Fm0,V,Ea),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,400])
nexttile
plot(t,y(:,1),'LineWidth',2)
xlabel('t(h)')
ylabel('C_A(lb-mol/ft^3)')
leg=legend('C_A');
leg.ItemTokenSize=[10,10];
title('Concentration-time curve')
axis([0 4 0 0.25])
grid on
nexttile
plot(t,y(:,5),'LineWidth',2)
xlabel('t(h)')
ylabel('T(^{\circ}F)')
axis([0 4 60 200])
leg=legend('T');
title('Temperature-time curve')
leg.ItemTokenSize=[10,10];
grid on
nexttile
plot(y(:,5),y(:,1),'LineWidth',2)
xlabel('T(^{\circ}F)')
ylabel('C_A(lb-mol/ft^3)')
axis([65 250 0 0.2])
title('Conc-temp trajectory')
grid on
z=size(y);
ThetaCp = 35 + Fb0 / Fa0 * 18 + Fm0 / Fa0 * 19.5; 
for i=1:z(1,1)
    % Calculating the values of Ta2 from the following expression
    k(i) = 16960000000000 * exp(-Ea / 1.987 / (y(i,5)+460)); 
ra(i) = 0 - (k(i) * y(i,1)); 
Qg(i)=ra(i) *V * dH;
Ta2(i) = y(i,5) - ((y(i,5) - Ta1) * exp(0 - (UA / (18 * mc)))); 
Qr2(i) = mc * 18 * (Ta2(i) - Ta1);
Qr1(i)=Fa0 * ThetaCp * (y(i,5) - T0);
Qr(i)=Qr1(i)+Qr2(i);
end
nexttile
plot(t,Qg(1:z(1,1)),t,Qr(1:z(1,1)),'LineWidth',2)
xlabel('t(hr)')
ylabel('Q(Btu/hr)')
axis([0 4 0 8e6])
leg=legend('Q_g','Q_{r}');
leg.ItemTokenSize=[10,10];
title('Heat-time trajectory')
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,UA,Ta1,T0,dH,mc,Fa0,Fb0,Fm0,V,Ea); 
Ca = Yfuncvec(1); 
Cb = Yfuncvec(2); 
Cc = Yfuncvec(3); 
Cm = Yfuncvec(4); 
T = Yfuncvec(5); 
% Explicit equations
k = 16960000000000 * exp(-Ea / 1.987 / (T + 460)); 
ra = 0 - (k * Ca); 
rb = 0 - (k * Ca); 
rc = k * Ca; 
Nm = Cm * V; 
Na = Ca * V; 
Nb = Cb * V; 
Nc = Cc * V; 
ThetaCp = 35 + Fb0 / Fa0 * 18 + Fm0 / Fa0 * 19.5; 
v0 = Fa0 / 0.923 + Fb0 / 3.45 + Fm0 / 1.54; 
Ta2 = T - ((T - Ta1) * exp(0 - (UA / (18 * mc)))); 
Ca0 = Fa0 / v0; 
Cb0 = Fb0 / v0; 
Cm0 = Fm0 / v0; 
Qr2 = mc * 18 * (Ta2 - Ta1);
Qr1=Fa0 * ThetaCp * (T - T0);
Qr=Qr1+Qr2;
Qg=ra * V*dH;
tau = V / v0; 
NCp = Na * 35 + Nb * 18 + Nc * 46 + Nm * 19.5; 
% Differential equations
dCadt = 1 / tau * (Ca0 - Ca) + ra; 
dCbdt = 1 / tau * (Cb0 - Cb) + rb; 
dCcdt = 1 / tau * (0 - Cc) + rc; 
dCmdt = 1 / tau * (Cm0 - Cm); 
dTdt = (Qg -Qr) / NCp; 
dYfuncvecdt = [dCadt; dCbdt; dCcdt; dCmdt; dTdt]; 
end