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

tspan = [0 125]; % Range for the independent variable
y0 = [448;0]; % Initial values for the dependent variables T and X
Nao=9.0448;
Nbo=33;
Nco=103.7;
UA=35.85;
DeltaHrx=-590000;
t1=45;
t2=55;
Ea=11273;
[t y]=ode23(@(t,y)ODEfun(t,y,Nao,Nbo,Nco,UA,DeltaHrx,t1,t2,Ea),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,400])
nexttile
plot(t,y(:,1),'LineWidth',2);
leg=legend('T');
leg.ItemTokenSize=[10,10];
title('Temperature-time trajectory')
ylabel('T(K)');
xlabel('t(mins)');
axis([0 125 400 1000]);
grid on
nexttile
plot(t,y(:,2),'LineWidth',2);
leg=legend('X');
leg.ItemTokenSize=[10,10];
title('Conversion-time trajectory')
ylabel('Conversion');
xlabel('t(mins)');
axis([0 125 0 1]);
grid on
z=size(y);
Vaqammonia=3.9;
VONCB= Nao*157.55/1199;
V = Vaqammonia + VONCB;
Theata=Nbo/Nao;
for i=1:z(1,1)
    k(i) = 0.00017 * exp((Ea / 1.987 )* (1 / 461 - (1 / y(i,1)))); 
    ra(i) = 0 - (k(i) * Nao ^ 2 * (1 - y(i,2)) * (Theata - (2 * y(i,2))) / (V ^ 2)); 
    Qg(i) = ra(i) * V * DeltaHrx;
 if (t(i) > t1 & t(i) < t2) 
    Qr(i) = 0;
else
    Qr(i) = UA * (y(i,1) - 298);
end 
end
nexttile
plot(t,Qg(1:z(1,1)),t,Qr(1:z(1,1)),'LineWidth',2);
leg=legend('Q_g','Q_r');
leg.ItemTokenSize=[10,10];
title('Heat-time trajectory')
ylabel('Q(kcal/min)');
xlabel('t(mins)');
axis([0 125 0 30000]);
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,Nao,Nbo,Nco,UA,DeltaHrx,t1,t2,Ea); 
T = Yfuncvec(1); 
X = Yfuncvec(2); 
% Explicit equations
k = 0.00017 * exp((Ea / 1.987) * (1 / 461 - (1 / T))); 
Theata=Nbo/Nao;
Vaqammonia=3.9;
VONCB= Nao*157.55/1199;
V = Vaqammonia + VONCB;
ra = 0 - (k * Nao ^ 2 * (1 - X) * (Theata - (2 * X)) / (V ^ 2)); 
NCp = Nao*40+Nbo*8.38+Nco*18;
Qg = ra * V * DeltaHrx;
 if (t > t1 & t < t2) 
    Qr = 0;
else
    Qr = UA * (T - 298);
end 

% Differential equations
if (t < t1) 
    dTdt = 0;
else
    dTdt = (Qg - Qr) / NCp;
end 
dXdt = (0 - ra) * V / Nao; 
dYfuncvecdt = [dTdt; dXdt];  
end