Figure 7.28 (page 407):

Dimensionless equilibrium constant and Thiele modulus versus reactor volume.

Code for Figure 7.28

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 P Rg T Ncof No2f NIf rhop rhob Dco k Kco a xco

%%
%% units: g, mol, sec, cm, K (but atm pressure)
%%

T = 838;
Rg =  82.06;
P  =  1;
ccof = 0.167*(P/(Rg*T));
co2f = 0.833*(P/(Rg*T));
cIf  = 0.0*(P/(Rg*T));
Kco = 8.099e6*exp(409/T);
% fix rate expression, 12/1/2013
%k = 1.048e12*exp(-13500/T)*co2f;
k = 1.3828e19*exp(-13500/T)*co2f;

%% fix feed flowrate
%Qf = 60;
Qf = 792e3;
Ncof = Qf*ccof;
No2f = Qf*co2f;
NIf =  Qf*cIf;
rhop = 0.68;
rhob = 0.60;
Dco = 0.0487;
Rp = 0.1;
a  = Rp/3;
xco = 0.95;



npts   = 100;
vfinal = 500;
vsteps = vfinal*linspace(0,1,npts)';
x0=Ncof;
opts = odeset ('Events', @stop, 'AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[vout,x] = ode15s (@pbr, vsteps, x0, opts);
if ( length(vout) == npts )
  fprintf ('hey, did not reach final conversion, increase rstop\n');
end
vplot = vout;
VR = vplot(end)
Nco = x;
No2  = No2f + (Nco-Ncof)/2;
Nco2 = Ncof - Nco;
Nt   = Nco + No2 + Nco2 + NIf;
cco  = P/(Rg*T) * Nco./Nt;
co2  = P/(Rg*T) * No2./Nt;
cco2 = P/(Rg*T) * Nco2./Nt;
phi = Kco*cco;
Phi = phi./(1+phi) .* sqrt( k*a*a./( 2*Dco*(phi-log(1+phi)) )  );
eta = 1.0 ./ Phi;
table = [vplot cco co2 cco2 phi Phi eta];
subplot (3, 1, 1);
plot (table(:,1), table(:,2:4));
title ('Figure 7.28')

subplot (3, 1, 2);
plot (table(:,1), table(:,5));
title ('Figure 7.28')

subplot (3, 1, 3);
plot (table(:,1), table(:,6));
title ('Figure 7.28')

stop.m

function [retval, isterminal, direction] = stop(t,x)
  global Ncof xco
  Nco = x(1);
  retval = Nco - (1-xco)*Ncof;
  isterminal = 1;
  direction = 0;

pbr.m

function xdot = pbr (t,x)
  global P Rg T Ncof No2f NIf rhop rhob Dco k Kco a
  Nco = x(1);
  Nt  = No2f + (Nco+Ncof)/2 + NIf;
  cco = P/(Rg*T) * Nco/Nt;
  phi = Kco*cco;
  Phi = phi/(1+phi) * sqrt( k*a*a/( 2*Dco*(phi-log(1+phi)) )  );
  if (Phi <= 1.0)
    eta = 1.0;
  else 
    eta = 1./Phi;
  end
  xdot = -rhob/rhop*eta*k*cco/(1+Kco*cco);