% Acetone Cracking to Ketene and Methane
% PFR with Counter-current Heat Exchange
% Simplified Model with Constant Heat Capacities
% Example from Fogler, 4th Ed., p. 504.
global R
R = 8.31446; % kJ/kgmol/K  also  m3*kPa/kgmol/K
Tref = 298.15;
% Basic data
MWA = 58.08;
MWAir = 28.96;

% Reactor
NoTubes = 1000;
TotalVolume = 2; % m3
VolPerTube = TotalVolume/NoTubes; % m3
TubeID = 26.7e-3; % m
TubeXC = pi*TubeID^2/4; % m2
TubeLength = VolPerTube/TubeXC; % m

% Heat Transfer
U = 400; % kJ/m2/hr/K
A = 4/TubeID; % m2/m3

% Acetone Feed
ReactorMassFeed = 7850; % kg/hr
FAM = ReactorMassFeed/NoTubes; % kg/hr per tube
FAF = FAM/MWA; % kgmol/hr
TF = 1035; % K
P = 162; % kPa
CAF = P/R/TF; % kgmol/m3

% Air Feed
AirFeed = 11088*8; % kg/hr
FaM = AirFeed/NoTubes; % kg/hr per tube
Fa = FaM/MWAir; % kgmol/hr
TaF = 1250; % K
Ta0 = 1112.9; % K -- estimate

% Solution span
vspan = linspace(0,VolPerTube,200);

% Initial conditions

y0 = [ FAF ; TF ; Ta0 ];

% anonymous function
PFRanon = @(v,y) PFRAcetoneSimplified(v,y,FAF,CAF,TF,Fa,U,A);

% solve ODEs
[ v, yout ] = ode45(PFRanon,vspan,y0);

% unpack output
FAout = yout(:,1);
FKout = FAF - FAout;
FMout = FAF - FAout;
Tout = yout(:,2);
Taout = yout(:,3);

% conversion
Conv = (FAF - FAout)/FAF;

% plot results
figure(1)
plot(v,FAout,'k-',v,FKout,'g-')
grid
axis([ 0 2e-3 0 0.14 ]);
xlabel('Reactor Volume - m3')
ylabel('Molar Flow Rate - kgmol/hr')
title('Molar Flow Rates')
legend('Acetone','Ketene','Location','northeastoutside')

figure(2)
plot(v,Tout,'b-',v,Taout,'r-')
grid
xlabel('Reactor Volume - m3')
ylabel('Temperature - K')
title('Temperature Profiles')
legend('Reaction','Air','Location','northwest')

figure(3)
plot(v,Conv,'g-')
grid
axis([ 0 2e-3 0 1 ]);
xlabel('Reactor Volume - m3')
ylabel('Conversion')
title('Conversion Profile')

disp('Air Entry Temperature')
disp(Taout(end))

