Figure 7.27:
Molar concentrations versus reactor volume.
Code for Figure 7.27
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
86
87
88
89
90
91 | % 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.
%
% Revised 8/15/2018
%
% units: g, mol, sec, cm, K (but atm pressure)
%
p.T = 838;
p.Rg = 82.06;
p.P = 1;
p.ccof = 0.167*(p.P/(p.Rg*p.T));
p.co2f = 0.833*(p.P/(p.Rg*p.T));
p.cIf = 0.0*(p.P/(p.Rg*p.T));
p.Kco = 8.099e6*exp(409/p.T);
% fix rate expression, 12/1/2013
%p.k = 1.048e12*exp(-13500/T)*co2f;
p.k = 1.3828e19*exp(-13500/p.T)*p.co2f;
% fix feed flowrate
%Qf = 60;
p.Qf = 792e3;
p.Ncof = p.Qf*p.ccof;
p.No2f = p.Qf*p.co2f;
p.NIf = p.Qf*p.cIf;
p.rhop = 0.68;
p.rhob = 0.60;
p.Dco = 0.0487;
p.Rp = 0.1;
p.a = p.Rp/3;
p.xco = 0.95;
npts = 100;
vfinal = 500;
vsteps = vfinal*linspace(0,1,npts)';
x0 = p.Ncof;
opts = odeset('Events', @(t, x) stop(t, x, p), 'AbsTol', sqrt(eps), ...
'RelTol', sqrt(eps));
[vout,x] = ode15s(@(t, x) pbr(t, x, p), 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 = p.No2f + (Nco - p.Ncof)/2;
Nco2 = p.Ncof - Nco;
Nt = Nco + No2 + Nco2 + p.NIf;
cco = p.P/(p.Rg*p.T) * Nco./Nt;
co2 = p.P/(p.Rg*p.T) * No2./Nt;
cco2 = p.P/(p.Rg*p.T) * Nco2./Nt;
phi = p.Kco*cco;
Phi = phi./(1 + phi) .* sqrt(p.k*p.a^2./...
(2*p.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
| function [retval, isterminal, direction] = stop(t, x, p)
Nco = x(1);
retval = Nco - (1 - p.xco)*p.Ncof;
isterminal = 1;
direction = 0;
|
pbr.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 | function xdot = pbr (t, x, p)
Nco = x(1);
Nt = p.No2f + (Nco + p.Ncof)/2 + p.NIf;
cco = p.P/(p.Rg*p.T) * Nco/Nt;
phi = p.Kco*cco;
Phi = phi/(1 + phi) * sqrt(p.k*p.a^2/...
(2*p.Dco* (phi - log(1 + phi)) ));
if (Phi <= 1.0)
eta = 1.0;
else
eta = 1./Phi;
end
xdot = -p.rhob/p.rhop*eta*p.k*cco/(1 + p.Kco*cco);
|