%  Startup of a CSTR, Adiabatic kinetics, multiple steady states
clc
clear all; close all;
k10 =17e11; %h-1
EbyR = 9000; %K

Q= 9.5; %m3/h
Volume = 1.25; %m3
CpA = 120;%J/mol/K
CpB = 60;
CpC = 60;
CpI = 100;

Fain = 20000;% mol/h
Fbin = 200000;
Fiin = 33000;

Ca0 = Fain/Q;
Na0 = Ca0*Volume;

Tin = 300;%K
Tc = 290;

U = 1e6; %J/s/m2/h
A = 0; %m2
delH = -85000; %J/molA reacted
Temperature = 280:0.1:400;

% mass balance equations
Tau = Volume/Q;
k = k10*exp(-EbyR./Temperature);
Xmb = k*Tau./(1+k*Tau);

% energy balance equations
SumFjCpJ = Fain * CpA + Fbin * CpB + Fiin * CpI;

Xeb = (SumFjCpJ * (Temperature-Tin)+U*A*(Temperature-Tc))/Fain/(-delH);

figure;
plot(Temperature,Xmb,'b');
hold on;
plot(Temperature,Xeb,'r');
xlabel('Temperature / K');
ylabel('Conversion');
legend('Mass Balance','Energy Balance');
set(gca,'ylim',[0 1]);


% Heat generated
G=Xmb*Fain*(-delH);
R=SumFjCpJ*(Temperature-Tin)+U*A*(Temperature-Tc);

figure;
plot(Temperature,G,'b');
hold on;
plot(Temperature,R,'r');
xlabel('Temperature / K');
ylabel('Heat / J');
legend('Generated','Removed');

clear Temperature;
clear k;

% Unsteady operation of the CSTR
t=[0 1.4];  % time horizon

X0=[2000 0 0 0 300.0]; % [Cai Cbi Cci Cii T0]

[t,X] = ode15s(@test,t,X0);

figure
plot (t(:,1),X(:,5));
xlabel('Time / h');
ylabel('Temperature / K');

% Temperature =302.5
% k = k10*exp(-EbyR./Temperature);
% Xmb = k*Tau./(1+k*Tau)
% Xeb = (SumFjCpJ * (Temperature-Tin)+U*A*(Temperature-Tc))/Fain/(-delH)
% 
% Temperature =340
% k = k10*exp(-EbyR./Temperature);
% Xmb = k*Tau./(1+k*Tau)
% Xeb = (SumFjCpJ * (Temperature-Tin)+U*A*(Temperature-Tc))/Fain/(-delH)
% 
% 
% Temperature =392.25
% k = k10*exp(-EbyR./Temperature);
% Xmb = k*Tau./(1+k*Tau)
% Xeb = (SumFjCpJ * (Temperature-Tin)+U*A*(Temperature-Tc))/Fain/(-delH)





