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