Figure 8.33:

Conversion of reactant for single, ideal CSTR, and as a function of internal flowrate in a 2-CSTR mixing model.

Code for Figure 8.33

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

global k1 k2 theta1 theta2 theta caf cbf rho alpha n ...
       a B

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


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


npts = 100;
rho0 = 0;
rhofin = 1;
rhovec = linspace(rhofin,rho0,npts)';
x0=[0; cbf-caf; 0; cbf-caf];
y0=[0; cbf-caf];
[y, fval, info] = fsolve('onereac',y0);
fsolve_failed = info <= 0;
if  (fsolve_failed)
  fprintf ('warning, fsolve failed, info = %d\n', info);
  rho, y0
end
conv1  = (alpha*caf-(1+alpha)*y(1))/(alpha*caf)*ones(npts,1);
yield1 = (cbf-(1+alpha)*y(2))/(alpha*caf-(1+alpha)*y(1))*ones(npts,1);
conv2 = zeros (npts, 1);
yield2 = zeros (npts, 1);
for i = 1:npts
  rho = rhovec(i);
  [x, fval, info] = fsolve('tworeac',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)  = (alpha*caf-(1+alpha)*x(3))/(alpha*caf);
  yield2(i) = (cbf-(1+alpha)*x(4))/(alpha*caf-(1+alpha)*x(3));

end
table1 =  [rhovec yield1 yield2 conv1 conv2];

%
% step test results
%



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

x0 = [0;0];
rho = 0;
a = [alpha/theta1; 0];
B = [-(alpha+rho)/theta1, rho/theta1;
      (alpha+rho)/theta2, -(1+alpha+rho)/theta2 ];
opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tsolver, x1] = ode15s(@twostep, times, x0, opts);
impx1 = [diag(a)*ones(2,nts)+B*x1']';
rho = 1;
a = [alpha/theta1; 0];
B = [-(alpha+rho)/theta1, rho/theta1;
      (alpha+rho)/theta2, -(1+alpha+rho)/theta2 ];
opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tsolver, x2]  = ode15s(@twostep, 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:3));
% 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
function res = tworeac (x)
  global k1 k2 theta1 theta2 caf cbf rho alpha n
  ca1 = x(1);
  cb1 = x(2);
  ca2 = x(3);
  cb2 = x(4);
  r11 = k1*ca1*cb1;
  r21 = k2*ca1^n;
  r12 = k1*ca2*cb2;
  r22 = k2*ca2^n;
  R = [ -(r11+r21)*theta1;
        -r11*theta1;
        -(r12+r22)*theta2;
        -r12*theta2];

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

onereac.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
function res = onereac (x)
  global k1 k2 theta caf cbf n alpha
  ca = x(1);
  cb = x(2);
  r1 = k1*ca*cb;
  r2 = k2*ca^n;
  R = [ -(r1+r2)*theta;
        -r1*theta];

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

twostep.m


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

onestep.m


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