Figure 7.10 (page 374):

Effectiveness factor versus Thiele modulus in a spherical pellet; reaction orders less than unity.

Code for Figure 7.10

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.

%%
%% jbr modifications to run on telemark; seems to have some parameters
%% set awfully high. Probably works on hazel with some better precision
%% than telemark.
%%
%% 6/4/00
%%

%% The pellet problem - Explanation of the structure of this program
%% This program computes the effectiveness factors for different thiele
%% moduli , for any order of the reaction and for the cases of
%% slab,cylinder and sphere(q=0,1,2).

%% The orders can be can be categorised into order=1 ,order<1 and
%% order>1. 

%% For order>1, the problem is set up as an initial value problem. We
%% start from the centre of the pellet (t=0) with the inner boundary condition
%% (dc/dr=0) and some c0=constant, and start marching forward in the
%% radial direction upto t=tsteps. The profile is then rescaled to satisfy
%% the outer boundary condition (c((q+1)*tsteps)=1). For different
%% values of tsteps we have different thiele modulus. We use the ODE
%% solver with stopping condition 
%% (ex.. thiele modulus=100). The solver marches forward until the
%% stopping condition is satisfied.(The stopping condition helps to
%% avoid the problem of blowing up)

%% For order=1, q=0 or 2, the analytical solutions have been used
%% For order=1, q=1 the solution is computed just like the case of
%% order>1 (Analytical solution for this is actually available, but it
%% involves Modified Bessels function which limits the value of thiele
%% modulus that can be attained.

%% For order<1 , the problem is solved in two stages
%% Stage 1. 
%% The same intitial value problem approach as explained abouve is used
%% , but  here we march upto t=large no.This would yield all values of
%% thiele  modulus for which the concentration at the centre is positive
%% Stage 2
%% Here the discontinuity in the concentration profile is
%% accomodated. We slightly modify the above initial value problem
%% approach. In this we start at a radius t=start, with the initial
%% conditions as c0=0 and a very small derivative, then we march forward
%% upto t=start+tsteps. The profile is rescaled to satisfy the outer
%% boundary condition c((q+1)*(start+tsteps)=1). Unlike the previous
%% case, the values of thiele modulus decreases for increase in
%% tsteps.The required solution is [Stage1;Stage2]

%% Limitations:

%% At extremely high thiele values (i.e >10000) for  order >1, the k and
%% c0 values chosen here will have to be changed suitably so that the
%% sensitivity of concentration w.r.t tsteps, is within the limit of
%% the ODE solver tolerances.
%% If one needs extremely high thiele modulus for order < 1, then ,the
%% initial derivative needs to be reduced and also the tsteps in stage 2
%% need to be very small.
%% Order =1,q=1 will not be able to yield very high thiele values. This
%% is because the thiele modulus is only proportional to tsteps, but
%% overflows occur at high tsteps.
  
global order k q thiele

%% Choose the orders , geometry (q=0,1,2 for slab,cylinder,sphere) and
%% the thiele modulus upto which you require eta values.(If you do not
%% have an order=1, q=1 case then you can choose values of thiele
%% modulus to be as high as 1000, else you will have to have thiele to
%% be not greater than 100.

ordervec=[0;0.25;0.5;0.75;1];
q=2;
thiele=10;

%% These are constants chosen arbitrarily
%% k     rate constant/diffusivity
%% c0    initial concentration at the centre
%% der   initial derivative for the case of dry regions in order <1
%% start the radius from which we start to march in the case of order<1 dry

k=2;
c0=0.01;
%der=1e-30;
der=1e-8;
start=1;

%% 
axis([0.1,thiele,1/thiele,1])
more off;
tab = [];
%%

%%
%% ODE's for order<1


%% Stopping Critierion for order>2


%% Main program

for j=1:size(ordervec,1)
  order=ordervec(j);

  if (order < 1)

    %% Stage 1

    tsteps=[0;logspace(-2,5,200)'];
    x0=[c0;0];
    ode_opts = odeset ('AbsTol', 1e-18, 'RelTol', sqrt (eps));
    [tsolver,x] = ode15s(@ivp,tsteps,x0,ode_opts);
    phi = sqrt((order+1)/2.*x(:,1).^(order-1)*k).*(tsteps)/(q+1);
    eta = x(:,2)./ (k*(tsteps).*x(:,1).^order) *(q+1);
    phi(1) =0;
    eta(1) = 1.0;

    %% Stage 2

%%   avoid duplicating time points in ODE solver call
    epsi = 1e-3;
    tsteps=[linspace(start,start+2,30)';logspace(log10(start+2+epsi),log10(40),20)'];
    x0=[0;der];
    ode_opts = odeset ('AbsTol', 1e-18, 'RelTol', sqrt (eps));
    [tsolver,x] = ode15s(@ivp,tsteps,x0,ode_opts);
    phidry = sqrt((order+1)/2.*x(:,1).^(order-1)*k).*(tsteps+start)/(q+1);
    etadry = x(:,2)./ (k*(tsteps+start).*x(:,1).^order) *(q+1);
    tab =[ tab  [phi , eta ;flipud(phidry(2:end)), ...
		 flipud(etadry(2:end))] ];
    
  elseif (order==1)
    phi=logspace(-2,2,250)';
    if q==0
      eta=tanh(phi)./phi;
    elseif q==1
      eta = besseli(2*phi,1)./phi./besseli(2*phi,1);
    elseif q==2
      eta=(1./phi).*((1./tanh(3*phi))-(1./(3*phi)));
    end
    
    tab=[ tab phi eta];
    
  elseif (order>1)
    x_0=[c0;0];
    y=k*c0^order;

    t_out=[0;logspace(-1,10,2000)'];

    ode_opts = odeset ('Events', @g, 'AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
    [t,x] = ode15s (@ivp, t_out, x_0, ode_opts);
    
    phi = sqrt((order+1)/2.*x(:,1).^(order-1)*k).*t./(q+1);
    eta = x(:,2)./ (k*t.*x(:,1).^order) *(q+1);
    tab=[ tab phi eta];
  end
 
end

loglog (tab(:,1), tab(:,2), ...
	tab(:,3), tab(:,4), ...
	tab(:,5), tab(:,6), ...
	tab(:,7), tab(:,8), ...
	tab(:,9), tab(:,10));
title ('Figure 7.10')

ivp.m

function xdot = ivp(t,x)
  global order q k
  xdot = zeros (2, 1);
  xdot(1) = x(2);
  if (t == 0)
    xdot(2) = k*x(1)^order;
  else
    xdot(2) = -q/t*x(2) + k*x(1)^order;
  end

g.m

function [stp, isterminal, direction] = g(t,x)
  global order q k thiele
  stp=sqrt((order+1)/2*x(1)^(order-1)*k)*t/(q+1)-thiele;

  isterminal = 1;
  direction = 0;