Figure 7.26:

Molar flow of A versus reactor volume for second-order, isothermal reaction in a fixed-bed reactor.

Code for Figure 7.26

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
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
% 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)
%

k = 2.25e5;
Nafin = 10;
T = 550;
rhop = 0.68;
rhob = 0.60;
Da = 0.008;
Rp = 0.45;
Rg =  82.06;
P  =  4;
n  = 2;
xa = 0.75;
nvs   = 100;
vfinal = 5e5;

p       = struct();
p.k     = k;
p.Nafin = Nafin;
p.T     = T;
p.rhop  = rhop;
p.rhob  = rhob;
p.Da    = Da;
p.Rp    = Rp;
p.Rg    = Rg;
p.P     = P;
p.n     = n;
p.xa    = xa;
p.nvs   = nvs;
p.vfinal= vfinal;



% Vanal = -2/(-n+1)*(2^((n+1)/2))*( (1-xa)^((-n+1)/2) - 1 )* ...
%          Nafin * Rp/3 * sqrt( (n+1)/2/Da) /                  ...
%          ( rhob/rhop*sqrt(k)*((P/(Rg*T))^((n+1)/2)) );

% Vanal2 = 4*( (1-xa)^(-1/2) - 1 )* ...
%          Nafin * Rp/3 * sqrt(3/Da) /  ...
%          ( rhob/rhop*sqrt(k)*((P/(Rg*T))^(3/2)) );

% solve reactor with: eta = 1./Phi*( 1./tanh(3*Phi) - 1/(3*Phi) );
vtotal = p.vfinal*linspace(0,1,p.nvs)';
vsteps = vtotal;

p.eta = @(x) 1./x*( 1./tanh(3*x) - 1/(3*x) );
x0 = p.Nafin;

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 numel(vout) == p.nvs
    fprintf ('hey, did not reach final conversion, increase stopping time\n');
end%if


% solve reactor with:  eta = 1./Phi;
vtotal = p.vfinal*linspace(0,1,p.nvs)';
vsteps = vtotal;

p.eta = @(x) 1./x;
x0 = p.Nafin;

opts = odeset ('Events', @(t,x) stop(t,x,p), 'AbsTol', sqrt(eps),...
               'RelTol', sqrt (eps));
[voutasy, xasy] = ode15s(@(t, x) pbr(t, x, p), vsteps, x0, opts);

if numel(vout) == p.nvs
    fprintf ('hey, did not reach final conversion, increase stopping time\n');
end%if

vout = vout/1000.;
VR = vout(end);

voutasy = voutasy/1000.;
VRasy = voutasy(end);
table = [vout x];
tableasy =  [voutasy xasy];
Naout = (1 - p.xa)*p.Nafin;
Natop = (1 - p.xa + 0.10)*p.Nafin;
dashedlines = ...
    [0,      Naout, VRasy, 0,     VR, 0,     VR, 2, VRasy, 2;
     1.1*VR, Naout, VRasy, Natop, VR, Natop, VR, 3, VRasy, 3];

save fborder2.dat table tableasy dashedlines

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
    plot (table(:,1), table(:,2), ...
           tableasy(:,1), tableasy(:,2), ...
           dashedlines(:,1), dashedlines(:,2), ...
           dashedlines(:,3), dashedlines(:,4), ...
           dashedlines(:,5), dashedlines(:,6));
    % TITLE
end % PLOTTING

pbrsolve.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
function [vout, x] = pbrsolve (par)
    vtotal = par.vfinal*linspace(0,1,par.nvs)';
    vsteps = vtotal;
    x0=par.Nafin;

    opts = odeset ('Events', @(t,x) stop (t,x,par), 'AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
    [vout, x] = ode15s (@(t,x) pbr (t,x,par), vsteps, x0, opts);

    if ( numel(vout) == par.nvs )
        fprintf ('hey, did not reach final conversion, increase stopping time\n');
    end%if
end%function

function xdot = pbr (t, x, par)
    Na = x(1);
    ca = par.P/(par.Rg*par.T) * Na/(2*par.Nafin);
    Phi = par.Rp/3*sqrt((par.n+1)/2*par.k*ca/par.Da);
    xdot = -par.rhob/par.rhop*par.eta(Phi)*par.k*ca^par.n;
end%function

function [retval, isterminal, direction] = stop(t, x, par)
    Na = x(1);
    retval = Na - (1-par.xa)*par.Nafin;
    isterminal = 1;
    direction = 0;
end%function

pbr.m

1
2
3
4
5
function xdot = pbr(t, x, p)
    Na = x(1);
    ca = p.P/(p.Rg*p.T) * Na/(2*p.Nafin);
    Phi = p.Rp/3*sqrt((p.n + 1)/2*p.k*ca/p.Da);
    xdot = -p.rhob/p.rhop*p.eta(Phi)*p.k*ca^p.n;

stop.m

1
2
3
4
5
function [retval, isterminal, direction] = stop(t, x, p)
    Na = x(1);
    retval = Na - (1 - p.xa)*p.Nafin;
    isterminal = 1;
    direction = 0;