## Code for Figure 7.26

### 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% 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.

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

par.k     = k;
par.Nafin   = Nafin;
par.T     = T;
par.rhop  = rhop;
par.rhob  = rhob;
par.Da    = Da;
par.Rp    = Rp;
par.Rg    = Rg;
par.P     = P;
par.n     = n;
par.xa    = xa;
par.nvs   = nvs;
par.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) );
par.eta = (@(x) 1./x*( 1./tanh(3*x) - 1/(3*x) ));
[vout, x] = pbrsolve(par);

% solve reactor with:  eta = 1./Phi;
par.eta = (@(x) 1./x);
[voutasy, xasy] = pbrsolve(par);

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

voutasy = voutasy/1000.;
VRasy = voutasy(end);
table = [vout x];
tableasy =  [voutasy xasy];
Naout = (1-par.xa)*par.Nafin;
Natop = (1-par.xa+0.10)*par.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
24function [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

```