Figure 7.28:

Dimensionless equilibrium constant and Thiele modulus versus reactor volume.

Code for Figure 7.28

Text of the GNU GPL.

main.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

global P Rg T Ncof No2f NIf rhop rhob Dco k Kco a xco

%
% units: g, mol, sec, cm, K (but atm pressure)
%

T = 838;
Rg =  82.06;
P  =  1;
ccof = 0.167*(P/(Rg*T));
co2f = 0.833*(P/(Rg*T));
cIf  = 0.0*(P/(Rg*T));
Kco = 8.099e6*exp(409/T);
% fix rate expression, 12/1/2013
%k = 1.048e12*exp(-13500/T)*co2f;
k = 1.3828e19*exp(-13500/T)*co2f;

% fix feed flowrate
%Qf = 60;
Qf = 792e3;
Ncof = Qf*ccof;
No2f = Qf*co2f;
NIf =  Qf*cIf;
rhop = 0.68;
rhob = 0.60;
Dco = 0.0487;
Rp = 0.1;
a  = Rp/3;
xco = 0.95;



npts   = 100;
vfinal = 500;
vsteps = vfinal*linspace(0,1,npts)';
x0=Ncof;
opts = odeset ('Events', @stop, 'AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[vout,x] = ode15s (@pbr, vsteps, x0, opts);
if ( length(vout) == npts )
  fprintf ('hey, did not reach final conversion, increase rstop\n');
end
vplot = vout;
VR = vplot(end)
Nco = x;
No2  = No2f + (Nco-Ncof)/2;
Nco2 = Ncof - Nco;
Nt   = Nco + No2 + Nco2 + NIf;
cco  = P/(Rg*T) * Nco./Nt;
co2  = P/(Rg*T) * No2./Nt;
cco2 = P/(Rg*T) * Nco2./Nt;
phi = Kco*cco;
Phi = phi./(1+phi) .* sqrt( k*a*a./( 2*Dco*(phi-log(1+phi)) )  );
eta = 1.0 ./ Phi;
table = [vplot cco co2 cco2 phi Phi eta];
save -ascii fbhw.dat table;
if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
subplot (3, 1, 1);
plot (table(:,1), table(:,2:4));
% TITLE fbhw

subplot (3, 1, 2);
plot (table(:,1), table(:,5));
% TITLE fbhw_phi

subplot (3, 1, 3);
plot (table(:,1), table(:,6));
% TITLE fbhw_phi
end % PLOTTING

stop.m


1
2
3
4
5
6
function [retval, isterminal, direction] = stop(t,x)
  global Ncof xco
  Nco = x(1);
  retval = Nco - (1-xco)*Ncof;
  isterminal = 1;
  direction = 0;

pbr.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
function xdot = pbr (t,x)
  global P Rg T Ncof No2f NIf rhop rhob Dco k Kco a
  Nco = x(1);
  Nt  = No2f + (Nco+Ncof)/2 + NIf;
  cco = P/(Rg*T) * Nco/Nt;
  phi = Kco*cco;
  Phi = phi/(1+phi) * sqrt( k*a*a/( 2*Dco*(phi-log(1+phi)) )  );
  if (Phi <= 1.0)
    eta = 1.0;
  else
    eta = 1./Phi;
  end
  xdot = -rhob/rhop*eta*k*cco/(1+Kco*cco);