Figure 9.35:
Predictions of LC measurement for reduced model.
Code for Figure 9.35
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 | % 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/29/2018
%
% vrdiclo initial moles of dichloro
% k1k2ratio k1/k2
% vro initial reactor volume cu dm
% teef TEA feed concentration kgmol/cu dm
% teaden TEA density kg/cu dm
% tflow, flowrate measured flowrate at time tflow
% read from "flow.dat" containing time
% (min) and flowrate (kg/min) converted to
% (cu dm/min)
% tlc, lcmeas lc measurement at time tlc (min)
%
% optimal values
% This m-file loads data file 'flow.dat'.
% This m-file loads data file 'lc.dat'.
p.vrdiclo = 2.3497;
p.k1k2ratio = 2;
p.vro = 2370;
p.teaf = 0.00721;
p.teaden = 0.728;
flow = load('flow.dat');
lc = load('lc.dat');
p.tflow = [0.; flow(:,1)];
p.flowrate = [0.; flow(:,2)./p.teaden];
tlc = lc(:,1);
lcmeas = lc(:,2);
ntimes = 80;
tlin = linspace(0, p.tflow(end), ntimes)';
p.timeout = sort(unique([tlin; p.tflow; tlc]));
% I can't think of a way to do this without a loop.
idx = zeros(size(tlc));
for i = 1:length (tlc)
% This assignment will fail if find returns more than one element,
% but that's OK, because each element of tlc should only appear once
% in timeout.
idx(i) = find (p.timeout == tlc(i));
end
% sanity check
if (any (p.timeout (idx) ~= tlc))
error ('extracting index vector for tlc');
end
x0 = [p.vro; 0];
opts = odeset('AbsTol', sqrt(eps));
ks = [0.5; 1; 2; 8; 32];
nks = length(ks);
for i = 1:nks
p.k1k2ratio = ks(i);
p.theta(1) = p.vrdiclo;
p.theta(2) = p.k1k2ratio;
[dummy, x] = ode15s(@(time, x) res_reduced(time, x, p), p.timeout, x0, opts);
% express the other states in terms of vr and eps2
vr = x(:,1);
eps2 = x(:,2);
badded = (vr - p.vro)*p.teaf;
eps1 = badded - eps2;
ca = (p.vrdiclo - eps1)./vr;
cb = badded - (eps1 + eps2); % sanity check, should be zero
cc = (eps1 - eps2)./vr;
ccd = eps2./vr;
lcpred(:,i) = cc./(cc + 2*ccd);
% avoid 0/0 at time = 0 and next few times if flowrate is zero and
% therefore x doesn't change;
firstposflow = min(find(p.flowrate > 0));
lcpred(1:firstposflow,i) = 1;
end
table = [p.timeout lcpred];
data = struct('lc', lc, 'table', table);
gnuplotsave('lc_reduced.dat', data);
if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
plot (table(:,1), table(:,2:6), lc(:,1), lc(:,2), '+');
% TITLE
end % PLOTTING
|
flow.dat
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 | 9.000000 9.7817838E-02 4.786000
19.00000 1.290601 14.21800
29.00000 1.262779 13.99800
39.00000 1.240016 13.81800
49.00000 1.215736 13.62600
59.00000 1.187156 13.40000
69.00000 1.160599 13.19000
79.00000 0.5545961 8.398001
89.00000 1.109508 12.78600
99.00000 1.361925 14.78200
109.0000 1.332333 14.54800
119.0000 1.294648 14.25000
129.0000 1.271885 14.07000
139.0000 1.244569 13.85400
149.0000 1.220794 13.66600
159.0000 1.197778 13.48400
169.0000 1.170463 13.26800
179.0000 1.150988 13.11400
189.0000 1.128983 12.94000
199.0000 1.105968 12.75800
209.0000 1.092563 12.65200
219.0000 1.360155 14.76800
229.0000 1.344473 14.64400
239.0000 1.329045 14.52200
249.0000 1.305524 14.33600
259.0000 1.275426 14.09800
269.0000 1.259239 13.97000
279.0000 1.244063 13.85000
289.0000 1.219783 13.65800
299.0000 1.191961 13.43800
309.0000 1.173751 13.29400
319.0000 1.158322 13.17200
329.0000 0.3163429 6.514000
339.0000 0.5495375 8.358001
349.0000 1.224841 13.69800
359.0000 0.3228684 6.565600
369.0000 9.4845297E-04 4.020000
379.0000 4.6794413E-04 4.016201
389.0000 5.4383278E-04 4.016800
399.0000 4.9326418E-04 4.016400
409.0000 3.6680698E-04 4.015400
419.0000 0.2701340 6.148601
429.0000 1.185891 13.39000
439.0000 0.9299334 11.36600
449.0000 0.5259147 8.171202
459.0000 8.9793204E-04 4.019600
469.0000 8.7256433E-04 4.019400
479.0000 7.9672335E-04 4.018800
489.0000 4.9319270E-04 4.016400
499.0000 4.6799183E-04 4.016200
509.0000 3.6678315E-04 4.015400
519.0000 0.2330809 5.855600
529.0000 0.7784327 10.16800
539.0000 7.5244429E-03 4.072000
549.0000 7.5244186E-03 4.072000
559.0000 5.7539940E-03 4.058000
569.0000 8.0302954E-03 4.076001
579.0000 0.5019882 7.982001
589.0000 -3.0983209E-03 3.988000
599.0000 9.4850064E-04 4.020000
609.0000 -6.3252446E-05 4.012000
619.0000 -1.3278484E-03 4.002000
629.0000 -4.6158792E-03 3.976000
639.0000 -6.3204767E-05 4.012000
649.0000 -2.0865917E-03 3.996000
659.0000 0.1752122 5.398000
669.0000 0.3024322 6.404000
679.0000 4.9952744E-03 4.052001
689.0000 3.9835931E-03 4.044001
699.0000 0.5114223 8.056601
709.0000 3.4151078E-04 4.015200
719.0000 3.1619071E-04 4.015000
729.0000 1.1386871E-04 4.013400
739.0000 6.3300133E-05 4.013000
749.0000 3.4151078E-04 4.015200
759.0000 0.1059619 4.850400
769.0000 6.7019463E-04 4.017800
779.0000 6.8099439E-02 4.551000
789.0000 8.7261200E-04 4.019400
799.0000 9.7372534E-04 4.020200
809.0000 8.7256433E-04 4.019400
819.0000 7.9669955E-04 4.018800
829.0000 7.2083471E-04 4.018200
839.0000 5.6912901E-04 4.017000
849.0000 4.6801567E-04 4.016200
859.0000 5.1856041E-04 4.016600
869.0000 4.4267176E-04 4.016000
|
lc.dat
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 | 414 0.17050
424 0.16040
434 0.13250
444 0.10840
493 0.10140
503 0.10420
513 0.10450
523 0.09700
533 0.08240
582 0.06780
592 0.06970
602 0.07020
612 0.06920
621 0.06840
631 0.06750
641 0.06730
651 0.06730
661 0.06540
671 0.05770
681 0.05710
690 0.05630
700 0.04660
710 0.04550
720 0.04490
730 0.04560
740 0.04500
750 0.04340
759 0.04280
769 0.04260
779 0.04090
789 0.03990
799 0.03970
848 0.03940
858 0.03740
868 0.03710
|
find_time_index.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 | % 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.
% Assumes T is sorted and contains no duplicate values.
function [index, is_switching_time] = find_time_index (t, now)
if (nargin ~= 2)
error ('usage: [index, is_switching_time] = find_time_index (t, now)');
end
if (~ isvector (t) || ~ isscalar (now))
error ('find_time_index: T must be a vector and NOW must be a scalar');
end
% use last flowrate past final time
if (now < t(1))
error ('find_time_index: NOW = %f is out of range [%f, %f]', ...
now, t(1), t(end));
end
if (now > t(end))
index = length(t);
else
index = max (find (t <= now));
end
is_switching_time = (t(index) == now);
end%function
|
res_reduced.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 | function xdot = res_reduced(time, x, p)
% Passed variable Local equivalent
%
% x(1) reactor contents volume
% x(2) extent of second reaction
%
vrdiclo = p.theta(1);
k1k2ratio = p.theta(2);
vr = x(1);
if ( vr <= 0.0 )
error ('(res_reduced).... vr is zero, time= %d',time);
end
eps2 = x(2);
badded = (vr - p.vro)*p.teaf;
eps1 = badded - eps2;
% look up flowrate at current time
[index, is_switching_time] = find_time_index(p.tflow, time);
fteaf = p.flowrate(index) * p.teaf;
% calculate xdots
xdot(1) = p.flowrate(index);
if (badded == 0.)
xdot(2) = 0.;
else
xdot(2) = fteaf/( 1 + p.k1k2ratio*(vrdiclo - (badded - eps2) )/...
(badded - 2*eps2) );
end
xdot = xdot(:);
|