Figure 8.34 (page 479):

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

%% 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)];

subplot (3, 1, 1);
plot (table1(:,1), table1(:,4:5));
title ('Figure 8.34')

subplot (3, 1, 2);
plot (table1(:,1), table1(:,2:3));
title ('Figure 8.34')

subplot (3, 1, 3);
plot (table2(:,1), table2(:,2:3));
title ('Figure 8.34')

tworeac.m

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

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

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

onestep.m

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