Figure 8.34:

Yield of desired product C for single, ideal CSTR, and as a function of internal flowrate, \rho =Q_r/Q_2, in a 2-CSTR mixing model.

Code for Figure 8.34

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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
% 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/16/2018

p.k1 = 1;
p.k2 = 2;
p.theta1 = 1;
p.theta2 = 2;
p.caf = 1;
p.cbf = 1;
p.rhotest = 1;
p.n   = 2;
p.alpha = 0.1;


%
% one well-mixed reactor
%
p.theta = p.theta1 + p.theta2;


npts = 100;
rho0 = 0;
rhofin = 1;
rhovec = linspace(rhofin,rho0,npts)';
x0 = [0; p.cbf - p.caf; 0; p.cbf - p.caf];
y0 = [0; p.cbf - p.caf];
[y, fval, info] = fsolve(@(x) onereac(x, p), y0);
fsolve_failed = info <= 0;

if  (fsolve_failed)
    fprintf ('warning, fsolve failed, info = %d\n', info);
    p.rho, y0
end

conv1  = (p.alpha*p.caf - (1 + p.alpha)*y(1))/...
            (p.alpha*p.caf)*ones(npts,1);
yield1 = (p.cbf - (1 + p.alpha)*y(2))/...
            (p.alpha*p.caf - (1 + p.alpha)*y(1))*ones(npts,1);
conv2 = zeros(npts, 1);
yield2 = zeros(npts, 1);

for i = 1:npts
    p.rho = rhovec(i);
    [x, fval, info] = fsolve(@(x) tworeac(x, p), x0);
    fsolve_failed = info <= 0;

    if (fsolve_failed)
        fprintf ('warning, fsolve failed, info = %d\n', info);
        fprintf('i= %g \n',i);
        %    rho, x0
    end

    x0 = x;
    conv2(i)  = (p.alpha*p.caf - (1 + p.alpha)*x(3))/...
                    (p.alpha*p.caf);
    yield2(i) = (p.cbf - (1 + p.alpha)*x(4))/...
                    (p.alpha*p.caf - (1 + p.alpha)*x(3));
end

table1 =  [rhovec yield1 yield2 conv1 conv2];

%
% step test results
%



nts = 100;
tfin = 10*p.theta;
times = linspace(0,tfin,nts)';
y0 = 0;
opts = odeset('AbsTol', sqrt(eps), 'RelTol', sqrt(eps));
[tsolver, y] = ode15s(@(t, x) onestep(t, x, p), times, y0, opts);
impy = onestep(0, y, p);

x0 = [0; 0];
p.rho = 0;
a = [p.alpha/p.theta1; 0];
B = [-(p.alpha + p.rho)/p.theta1, p.rho/p.theta1;
      (p.alpha + p.rho)/p.theta2, ...
            -(1 + p.alpha + p.rho)/p.theta2];
opts = odeset('AbsTol', sqrt(eps), 'RelTol', sqrt(eps));
[tsolver, x1] = ode15s(@(t, x) twostep(t, x, a, B), times, x0, opts);
impx1 = [diag(a)*ones(2,nts) + B*x1']';

p.rho = 1;
a = [p.alpha/p.theta1; 0];
B = [-(p.alpha + p.rho)/p.theta1, p.rho/p.theta1;
      (p.alpha + p.rho)/p.theta2, ...
            -(1 + p.alpha + p.rho)/p.theta2];
opts = odeset('AbsTol', sqrt(eps), 'RelTol', sqrt(eps));
[tsolver, x2] = ode15s(@(t, x) twostep(t, x, a, B), times, x0, opts);
impx2 = [diag(a)*ones(2,nts) + B*x2']';

table2 = [times y x1(:,2) x2(:,2) impy impx1(:,2) impx2(:,2)];
save selectivity.dat table1 table2

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
    subplot (3, 1, 1);
    plot (table1(:,1), table1(:,4:5));
    % TITLE selectivity_1

    subplot (3, 1, 2);
    plot (table1(:,1), table1(:,2:3));
    % TITLE selectivity_2

    subplot (3, 1, 3);
    plot (table2(:,1), table2(:,2:4));
    % TITLE selectivity_3
end % PLOTTING

tworeac.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
function res = tworeac(x, p)
    ca1 = x(1);
    cb1 = x(2);
    ca2 = x(3);
    cb2 = x(4);

    r11 = p.k1*ca1*cb1;
    r21 = p.k2*ca1^p.n;
    r12 = p.k1*ca2*cb2;
    r22 = p.k2*ca2^p.n;

    R = [ -(r11 + r21)*p.theta1;
          -r11*p.theta1;
          -(r12 + r22)*p.theta2;
          -r12*p.theta2];

    flow = [p.alpha*p.caf - (p.alpha + p.rho)*ca1 + p.rho*ca2;
          -(p.alpha + p.rho)*cb1 + p.rho*cb2;
          -(1 + p.alpha)*ca2 + (p.alpha + p.rho)*ca1 - p.rho*ca2;
           p.cbf + (p.alpha + p.rho)*cb1 - ...
                p.rho*cb2 - (1 + p.alpha)*cb2 ];
    res = flow + R;

onereac.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
function res = onereac(x, p)
    ca = x(1);
    cb = x(2);
    r1 = p.k1*ca*cb;
    r2 = p.k2*ca^p.n;

    R = [-(r1 + r2)*p.theta;
         -r1*p.theta];

    flow = [p.alpha*p.caf - (1 + p.alpha)*ca;
            p.cbf - (1 + p.alpha)*cb];
    res = flow + R;

twostep.m

1
2
function rhs = twostep(t, x, a, B)
    rhs = a + B*x;

onestep.m

1
2
function rhs = onestep(t, x, p)
    rhs = (p.alpha - (1 + p.alpha)*x)  / (p.theta1 + p.theta2);