% Lecture 9, Isothermal Reactors PFRs
%
% Example 1 A+2B --> C
home
format compact
close all
clear

nu =[-1 -2 1];
k = 0.50 ; % l/mol s

CA_in=1.0; % mol/l
vA0 = 2; % l/s
CB_in=2.5; % mol/l
vB0 = 2; % l/s
v0=vA0+vB0;

F0=[CA_in*vA0 CB_in*vB0 0];
C0=F0/v0;
sprintf('Initial concentrations %f %f %f M',C0)

lim=-F0./nu;
% A limiting in forward direction

C=@(X) C0-(nu/nu(1))*X*C0(1);
CA=@(X) C0(1)*(1-X);
CB=@(X) C0(2)-(nu(2)/nu(1))*X*C0(1);
CC=@(X) C0(3)-(nu(3)/nu(1))*X*C0(1);

nrA=@(X) 1./(k*CA(X).*CB(X));

CA0=C0(1);
tau0=CA0*quad(nrA,0,.9);
Volume=v0*tau0;

disp(sprintf('Residence time %f s, Volume %f l',tau0,Volume))

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')

subplot(3,1,1);
fplot(C,[0 1])
title('A + 2 B \rightarrow C at constant density')
ylabel('Concentrations (M)')
legend('C_A','C_B','C_C')
%
%  Reciprocal of rate:
subplot(3,1,2);
fplot(nrA,[0 .95])
title('Levenspiel plot for r=kC_AC_B')
ylabel('-1/rate')
axis([0 1 0 250])

%

dXdtau=@(t,X) (k/CA0)*CA(X).*CB(X);
tau_range=[0 20];
init=0;
% concentrations vs. tau
[tau X]=ode45(dXdtau,tau_range,init);

V=tau*v0;
subplot(3,1,3)
plot(X,V);
axis([0 0.95 0 50])
title('Volume vs. Desired Conversion')
xlabel('Conversion'); ylabel('Volume (l)');
axis([0 1 0 50])

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
plot(V,CA(X),V,CB(X),V,CC(X))
title('Constant density A + 2 B \rightarrow C in PFR')
xlabel('Reactor volume (l)'); ylabel('Concentration (M)')
legend('C_A','C_B','C_C')

%
% Levenspiel plots for various orders
%
% h=figure('color',[1 1 1]);
% set(h,'WindowStyle','Docked')
% Xl=[0:0.02:0.98];
% first=1./(1-Xl); second=first.^2; third=first.^3;
% plot(Xl,first,Xl,second,Xl,third)
% axis([0 1 0 100])
% xlabel('Conversion');ylabel('1/rate');title('Comparison of reaction orders in CSTR')
% legend('1st order','2nd order','3rd order','Location','Northwest')

%
%  Variable density reactor
%
F=@(X) F0-(nu/nu(1))*X*F0(1);
FA=@(X) F0(1)*(1-X);
FB=@(X) F0(2)-(nu(2)/nu(1))*X*F0(1);
FC=@(X) F0(3)-(nu(3)/nu(1))*X*F0(1);
yA0=FA(0)/sum(F(0));

R0=0.0821; T=400;
P=sum(F(0))*R0*T/v0

v0v=@(X) 1-sum(nu)*yA0*X/nu(1);

Cv=@(X) F(X)./v0v(X);
CAv = @(X) CA(X)./v0v(X);
CBv = @(X) CB(X)./v0v(X);
CCv = @(X) CC(X)./v0v(X);

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
subplot(3,1,1);
fplot(Cv,[0 1])
title('A + 2 B \rightarrow C at constant pressure')
ylabel('Concentrations (M)')
hold all
fplot(v0v,[0 1])
legend('C_A','C_B','C_C','v/v_0')

subplot(3,1,2)
%title('Volumetric flow rate change vs. conversion')

nrAv =@(X) F0(1)./(k*CAv(X).*CBv(X));
fplot(nrAv,[0 0.99])
ylabel('1/rate')
title('Levenspiel plot for r=kC_AC_B')
axis([0 1 0 250])

Volume=quad(nrAv,0,.9);
space_time=Volume*v0;
disp(sprintf('Space time %f s, Volume %f l',space_time,Volume))

dXdV=@(t,X) (k/F0(1))*CAv(X).*CBv(X);

V_range=[0 40];
init=0;
% concentrations vs. volume
[Vd Xd]=ode45(dXdV,V_range,init);

subplot(3,1,3)
plot(Xd,Vd)
xlabel('Conversion');ylabel('Volume (l)')
title('Volume vs. Desired Conversion')
axis([0 1 0 50])

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
subplot(2,1,1)
plot(V,CA(X),V,CB(X),V,CC(X))
title('Constant density A + 2 B \rightarrow C')
ylabel('Concentration (M)')
legend('C_A','C_B','C_C')
axis([0 40 0 1.5])
subplot(2,1,2)
plot(Vd,CAv(Xd),Vd,CBv(Xd),Vd,CCv(Xd))
title('Constant pressure A + 2 B \rightarrow C in PFR')
xlabel('Reactor volume (l)'); ylabel('Concentration (M)')
axis([0 40 0 1.5])

h=figure('color',[1 1 1]);
set(h,'WindowStyle','Docked')
subplot(2,1,1)
plot(V,v0*CA(X),V,v0*CB(X),V,v0*CC(X))
title('Constant density A + 2 B \rightarrow C')
legend('F_A','F_B','F_C')
ylabel('Molar flow rate (mol/s)')
axis([0 40 0 5])
subplot(2,1,2)
plot(Vd,FA(Xd),Vd,FB(Xd),Vd,FC(Xd))
title('Constant pressure A + 2 B \rightarrow C in PFR')
xlabel('Reactor volume (l)'); ylabel('Molar flow rate (mol/s))')
axis([0 40 0 5])
ans =
Initial concentrations 0.500000 1.250000 0.000000 M
Residence time 8.236955 s, Volume 32.947821 l
P =
   57.4700
Space time 50.237135 s, Volume 12.559284 l