Figure 8.19:

Dimensionless effluent concentration versus dimensionless rate constant for second-order reaction.

Code for Figure 8.19

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.

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

% Bradley Schwab reports that Matlab fails with ode15s,
% so try implicit solver ode15i instead.
%
% jbr, 8/13/2018

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;
    dcbar = 0;

    [dummy,cmm] = ode15i(@(z,c,dcdz) mm_imp(z,c,dcdz,K), z, cbar, dcbar, opts);
    cmmfin(i) = cmm(npts);

    [dummy,cseg] = ode15s(@(z, x) seg(z, x, K), 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_imp.m

1
2
function res = mm_imp(z, c, dcdz, K)
    res = dcdz*z^2 + 4*(1-z)/(2-z)*(c-1) + K*c*c;

seg.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
function xdot = seg(z, x, 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