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(:);