Figure 7.13 (page 380):

Effectiveness factor versus an inappropriate Thiele modulus in a slab; Hougen-Watson kinetics.

Code for Figure 7.13

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.

%%
%% Try to solve for the dimensionless pellet profile using
%% a Langmuir Hinshelwood mechanism.
%%
%% d^2 c/ d r^2 -Phi^2 c / (1 + phi c) = 0
%% c(1)      = 1
%% dc/dr (0) = 0
%%
%% goal is to not use a shooting method to avoid numerical problems
%% try it as a pure ODE
%%
%% jbr, 6/12/01
%%

%% first effort: consider dimensional problem
%%
%% x1 = c_A
%% x2 = dc_A/dr
%%
%% dx2/dr = p1 p2 c_A / (1 + p2 c_A)
%% dx1/dr = x2
%% 
%% x1(0) = c0 
%% x2(0) = 0
%%
%% solve this IVP to r=r1 and x1=c1
%%
%% p1 = k cm/De
%% p2 = KA
%% Phi^2 = p1*p2*r1^2
%% phi = p2*c1
%%
%% problem up to here is both Phi and phi will change with r1,c1
%% make p2 change with time as an ODE as well, to keep phi=constant
%% dp2/dr = -p2*x2/x1
%% p2(0)  = phi/x1(0)
%%
%% OK, don't think the above idea of changing p2 is solving the right problem.
%% Just march the ODE to desired phi, change p2 and solve ODE again
%% to same phi.  The two values of Phi should be different, however, 
%% and that allows a family of Phi at constant phi. Drawback is that
%% we are solving the ODE a lot.
%%
more off;
global p1 p2 phi
p1=1;
p2=1;
c0 = 0.1;
x0=[c0;0;0;0];
npts=200;
rstop = 10;
rsteps=linspace(0.001,rstop,npts)';


phivec = [0.1 10 100 1000];
nphi = length(phivec);
results = [];
for j= 1:nphi
  %% March ODE solver out to phi stopping condition for a range of p2
  %%
  %% next pairs go together for smooth log log plot of etavec vs. Phivec 
  %%
  phi=phivec(j);
  p2vec=[1e-40 1e-10 .1 1 10 20 30 40 50 55 60 65 70 75 80 90 95 99 99.5];
  if (phi >= 100)
    p2vec(1)=1e-20;
  end
  np2=length(p2vec);
  c0=1e-2*phi;
  Phivec = zeros(np2,1);
  etavec = zeros(np2,1);
  for i = 1:np2
    p2=p2vec(i);
    rstop = max(100,100*sqrt(1/p2));
    rsteps=[0;linspace(0.001,rstop,npts)'];
    x0 = [c0; 0; 0; 0];
    opts = odeset ('Events', @g, 'AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
    [rout,x] = ode15s (@hougensrt, rsteps, x0, opts);
    if ( length(rout) == npts + 1)
      fprintf ('hey, did not reach requested phi value, increase rstop\n');
    end
    nout = length(rout);
    cbar = x(:,1)/x(nout,1);
    rbar = rout/rout(nout);
    %%  plot (rbar,cbar);
    Phivec(i) = sqrt(p1*p2)*rout(nout);
    etavec(i) = (1+phi)*x(nout,2)*rout(nout)/x(nout,1)/Phivec(i)^2;
    %%
    %% check, should be zero if solved the dim'less  model 
    %% hmm, apparently cannot trust the last xdot.  Could make a call 
    %% to hougensrt to correct this last  one if this were important
    %% check passed, 6/13/01
    %%
    %% ode15s doesn't return xdot at all the solution times, so would
    %% need to compute it here to check.
    %%
    %% check = rout(nout)^2/x(nout,1)*xdot(:,2) - ...
    %% Phivec(i)^2*cbar./(1+phi*cbar);
%%    fprintf ('norm of error in solution: %g\n', max(abs(check(1:nout-1))));
  end
  Phiprimevec = phi/(1+phi) / sqrt( 2*(phi-log(1+phi)) ) * Phivec;
  results = [results Phivec Phiprimevec etavec];
end

subplot (2, 1, 1);
loglog (results(:,1), results(:,3), ...
	results(:,4), results(:,6), ...
	results(:,7), results(:,9), ...
	results(:,10), results(:,12));
title ('Figure 7.13')

subplot (2, 1, 2);
loglog (results(:,2), results(:,3), ...
	results(:,5), results(:,6), ...
	results(:,8), results(:,9), ...
	results(:,11), results(:,12));
title ('Figure 7.13')

hougensrt.m

function xdot = hougensrt(t,x)
  global p1 p2
  ca    = x(1);
  dcadr = x(2);
  s1    = x(3);
  s2    = x(4);
  xdot = [dcadr; ...
	  p1*p2*ca/(1+p2*ca); ...
	  s2; ...
	  p1*ca/(1+p2*ca)^2 + p1*p2*s1/(1+p2*ca)^2];

g.m

function [stp, isterminal, direction] = g(t,x)
  global p2 phi
  phipellet = p2*x(1);
  stp = phipellet-phi;

  isterminal = 1;
  direction = 0;