% LEP 12-2d : Production of Acetic Anhydride (Counter-current)
%

Vspan = [0 0.001]; % Range for the independent variable i.e.volume of the reactor
FA0=0.0376;
FI0=0;
Cpa=163;
Cpb=83;
CT0=18.8;
mc=0.111;
Cpcoolant=34.5;
Cpc=71;
Ua=16500;
Ta0=995.15;
To=1035;
CpI=32.68;

y0 = [Ta0;To; 0]; % Initial values for the dependent variables Ta,T, and x
[v y]=ode45(@(v,y)ODEfun(v,y,FA0,FI0,Cpa,Cpb,CT0,mc,Cpcoolant,Cpc,To,Ua,CpI),Vspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,600,400])
nexttile
plot(v,y(:,2),v,y(:,1),'LineWidth',1.5);
leg=legend('T','T_a');
ylabel('T(K)');
xlabel('V(m^3)');
axis([0 0.001 900 1300]);
leg.ItemTokenSize=[10,10];
title('Temperature profile')
grid on
z=size(y);
delCp=Cpb+Cpc-Cpa;
FT0 = FA0 + FI0;
vo = (FT0/CT0)*1000;
Cao = (FA0/vo)*1000;
for i=1:z(1,1)
    ra(i) = -Cao*3.58*exp(34222*(1/To-1/y(i,2)))*(1-y(i,3))*(To/y(i,2))/(1+y(i,3));
    deltaH(i)=80770+delCp*(y(i,2)-298);
Qg(i)= ra(i) * deltaH(i);
Qr(i)=(Ua * (y(i,2) - y(i,1)));
end
nexttile
plot(v,y(:,3),'LineWidth',1.5)
leg=legend('X');
leg.ItemTokenSize=[10,10];
title('Conversion profile')
ylabel('Conversion');
xlabel('V(m^3)');
axis([0 0.001 0 1]);
grid on
nexttile
plot(v,-ra(1:z(1,1)),'LineWidth',1.5)
leg=legend('-r_A');
leg.ItemTokenSize=[10,10];
title('Rate profile')
ylabel('Rate (mol/m^3.s)');
xlabel('V(m^3)');
axis([0 0.001 0 100]);
grid on
nexttile
plot(v,Qg(1:z(1,1)),v,Qr(1:z(1,1)),'LineWidth',1.5)
leg=legend('Q_g','Q_r');
leg.ItemTokenSize=[10,10];
title('Heat profile')
ylabel('Q(J/m^3.s)');
xlabel('V(m^3)');
axis([0 0.001 -1e7 2e6]);
grid on

%

function dYfuncvecdV = ODEfun(V,Yfuncvec,FA0,FI0,Cpa,Cpb,CT0,mc,Cpcoolant,Cpc,To,Ua,CpI); 
Ta = Yfuncvec(1); 
T = Yfuncvec(2); 
X = Yfuncvec(3); 
% Explicit equations
delCp=Cpb+Cpc-Cpa;
thetaI=FI0/FA0;
FT0=FA0+FI0;
vo=(FT0/CT0)*1000;
Cao=(FA0/vo)*1000;
deltaH=80770+delCp*(T-298);
ra= -Cao*3.58*exp(34222* (1/To - 1/T))*(1 - X)*(To/T)/(1 + X);
Qg=ra*deltaH;
Qr=Ua*(T-Ta);
sumcp=Cpa+thetaI*CpI;

% Differential equations
dTadV = -Ua* (T - Ta) /mc/Cpcoolant; 
dTdV = (Qg-Qr) / (FA0*(sumcp+X*delCp)); 
dXdV = 0 - (ra / FA0); 
dYfuncvecdV = [dTadV; dTdV; dXdV];  
end