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
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

1
2
3
4
5
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);