Figure 8.19 (page 459):

Dimensionless effluent concentration versus dimensionless rate constant for second-order reaction.

Code for Figure 8.19

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.

%%
%% Compute the solution to maximum mixedness, segregated and ideal
%% mixing models for 2 cstrs in series. 
%%
%% jbr, 3/19/00.
%%
%% maximum mixedness case
%% dc/dt --> 0 as t --> \infty
%%
%%  assume f(t)     = p(t)           is finite as t --> \infty, f(\infty)=P
%%                   ------
%%           \int_t^\infty p(t) dt
%%
%% Let cbar be solution to P(c-c0) + R(c) = 0
%%
%% Using transformation z=1/(1+t) or t=(1-z)/z
%% puts t \in [0,\infty) onto z \in (0,1]
%% transformed ode is
%% dc/dz = - 1/z^2 [ f((1-z)/z)(c-c0) + R(c)]
%%
%% example in Table 1 from Zwietering, 2 cstrs
%% f(t) = 4t/(2t+1) (so f(\infty)=2), R(c)=Kc^2, c0=1, 
%% cbar = -1+sqrt(1+K)), choose positive sign.
%%
global K
Kmin = 1e-2;
Kmax = 1e3;
nks=40;
Kseq=logspace(log10(Kmin),log10(Kmax),nks)';
%%
%% define mm and seg models 
%%

npts=50;
z=linspace(0,1,npts);
tcrit=1.0;
c0=[1;0];
opts = odeset ('AbsTol', sqrt (eps));

cmmfin = zeros (nks, 1);
csegfin = zeros (nks, 1);

for i = 1:nks
  K=Kseq(i);
  cbar = ( -1+sqrt(1+2*K) )/K;
  [dummy,cmm] = ode15s(@mm,z,cbar,opts);
  cmmfin(i) = cmm(npts);
  [dummy,cseg] = ode15s(@seg,z,c0,opts);
  csegfin(i) = cseg(npts,2);
end
%%
%%compute result for 2 ideal CSTRs
%%
cbar1 = (-1+sqrt(1 + 2*Kseq)) ./ (Kseq);
cbar2 = (-1+sqrt(1 + 2*Kseq.*cbar1)) ./ (Kseq);
table = [Kseq cmmfin csegfin cbar2];

semilogx (table(:,1), table(:,2:4));
title ('Figure 8.19')

mm.m

function dcdz = mm(z,c)
global K
  if (z == 0.)
    dcdz = 0;
  else
    t = (1-z)/z;
    dcdz = - (4*t/(2*t+1)*(c-1) + K*c*c)/(z*z);
  end

seg.m

function xdot = seg(z,x)
  global K 
  c    = x(1);
  ctot = x(2);
  if (z < 1) 
    t = z/(1-z);
    xdot = [-K*c*c/(1-z)^2; 4*t*exp(-2*t)/(1-z)^2*c];
  else
    xdot = [0; 0];
  end