% LEP 13-6 : Safety in Chemical Plants with Exothermic Runaway Reactions
%

tspan = [0 4.]; % Range for the independent variable time
y0 = [4.3; 5.1; 3.; 4.4; 422.]; % Initial values for the dependent variables i.e CA,CB,CS,P,and T
E1A=128000;
E2S=800000;
DHRx1A=-45400;
DHRx2S = -320000;
UA=0;
pset=28.2;
VH=5000;
Cv1=3360;
Cv2=53600;
V0 = 4000; 
Ta=373.15;
SumNCp = 12600000; 
A1A = 400000000000000; 
A2S = 1E+84;
[t y]=ode45(@(t,y)ODEfun(t,y,E1A,E2S,DHRx1A,DHRx2S,UA,pset,VH,Cv1,Cv2,V0,Ta,SumNCp,A1A,A2S),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,400])
nexttile
plot(t,y(:,5),'LineWidth',2)
xlabel('t(hours)')
ylabel('T(K)')
title('Temperature-time trajectory')
axis([0 4 400 600]);
grid on
nexttile
plot(t,y(:,4),'LineWidth',2)
xlabel('t(hours)')
ylabel('P(atm)')
title('Pressure-time trajectory')
axis([0 4 0 50]);
grid on
nexttile
plot(t,y(:,1),t,y(:,2),t,y(:,3),'LineWidth',2)
xlabel('t(hours)')
ylabel('C_i(mol/dm^3)')
leg=legend('C_A','C_B','C_S');
title('Concentration-time trajectory')
leg.ItemTokenSize=[10,10];
axis([0 4 0 6]);
grid on
z=size(y);
for i=1:z(1,1)
k1A(i) = A1A * exp(0 - (E1A / (8.31 * y(i,5)))); 
k2S(i) = A2S * exp(0 - (E2S / (8.31 * y(i,5)))); 
r1A(i) = 0 - (k1A(i) * y(i,1) * y(i,2)); 
r2S(i) = 0 - (k2S(i) * y(i,3));
Qg(i)=V0 * (r1A(i) * DHRx1A + r2S(i) * DHRx2S);
Qr(i)= UA * (y(i,5) - Ta);
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];
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,E1A,E2S,DHRx1A,DHRx2S,UA,pset,VH,Cv1,Cv2,V0,Ta,SumNCp,A1A,A2S); 
CA = Yfuncvec(1); 
CB = Yfuncvec(2); 
CS = Yfuncvec(3); 
P = Yfuncvec(4); 
T = Yfuncvec(5); 
% Explicit equations
 
k1A = A1A * exp(0 - (E1A / (8.31 * T))); 
k2S = A2S * exp(0 - (E2S / (8.31 * T))); 
if (T > 600 | P > 45) 
    SW1 = 0;
else
    SW1 = 1;
end 
r1A = 0 - (k1A * CA * CB); 
r2S = 0 - (k2S * CS); 
FD = (-0.5 * r1A - (3 * r2S)) * V0; 
Qg=V0 * (r1A * DHRx1A + r2S * DHRx2S);
Qr= UA * (T - Ta);
if (FD < 11400) 
    Fvent = FD;
else
    if (P < pset) 
        Fvent = (P - 1) * Cv1;
    else
        Fvent = (P - 1) * (Cv1 + Cv2);
    end
end 
% Differential equations
dCAdt = SW1 * r1A; 
dCBdt = SW1 * r1A; 
dCSdt = SW1 * r2S; 
dPdt = SW1 * (FD - Fvent) * 0.082 * T / VH; 
dTdt = SW1 * (Qg-Qr) / SumNCp; 
dYfuncvecdt = [dCAdt; dCBdt; dCSdt; dPdt; dTdt]; 

end