Code for Figure 8.19

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

%
% Compute the solution to maximum mixedness, segregated and ideal
% mixing models for 2 cstrs in series.
%
% jbr, 3/19/00.
%
% maximum mixedness case
% dc/dt --> 0 as t --> \infty
%
%  assume f(t)     = p(t)           is finite as t --> \infty, f(\infty)=P
%                   ------
%           \int_t^\infty p(t) dt
%
% Let cbar be solution to P(c-c0) + R(c) = 0
%
% Using transformation z=1/(1+t) or t=(1-z)/z
% puts t \in [0,\infty) onto z \in (0,1]
% transformed ode is
% dc/dz = - 1/z^2 [ f((1-z)/z)(c-c0) + R(c)]
%
% example in Table 1 from Zwietering, 2 cstrs
% f(t) = 4t/(2t+1) (so f(\infty)=2), R(c)=Kc^2, c0=1,
% cbar = -1+sqrt(1+K)), choose positive sign.
%
global K
Kmin = 1e-2;
Kmax = 1e3;
nks=40;
Kseq=logspace(log10(Kmin),log10(Kmax),nks)';
%
% define mm and seg models
%

npts=50;
z=linspace(0,1,npts);
tcrit=1.0;
c0=[1;0];
opts = odeset ('AbsTol', sqrt (eps));

cmmfin = zeros (nks, 1);
csegfin = zeros (nks, 1);

for i = 1:nks
K=Kseq(i);
cbar = ( -1+sqrt(1+2*K) )/K;
[dummy,cmm] = ode15s(@mm,z,cbar,opts);
cmmfin(i) = cmm(npts);
[dummy,cseg] = ode15s(@seg,z,c0,opts);
csegfin(i) = cseg(npts,2);
end
%
%compute result for 2 ideal CSTRs
%
cbar1 = (-1+sqrt(1 + 2*Kseq)) ./ (Kseq);
cbar2 = (-1+sqrt(1 + 2*Kseq.*cbar1)) ./ (Kseq);
table = [Kseq cmmfin csegfin cbar2];

save -ascii mmseg2cstr.dat table;
if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
semilogx (table(:,1), table(:,2:4));
% TITLE
end % PLOTTING



mm.m


1
2
3
4
5
6
7
8function dcdz = mm(z,c)
global K
if (z == 0.)
dcdz = 0;
else
t = (1-z)/z;
dcdz = - (4*t/(2*t+1)*(c-1) + K*c*c)/(z*z);
end



seg.m


1
2
3
4
5
6
7
8
9
10function xdot = seg(z,x)
global K
c    = x(1);
ctot = x(2);
if (z < 1)
t = z/(1-z);
xdot = [-K*c*c/(1-z)^2; 4*t*exp(-2*t)/(1-z)^2*c];
else
xdot = [0; 0];
end