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.


%% 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
%% 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.


%% 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


more off;
tab = [];

%% ODE's for order<1

%% Stopping Critierion for order>2

%% Main program

for j=1:size(ordervec,1)

  if (order < 1)

    %% Stage 1

    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;
    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)
    if q==0
    elseif q==1
      eta = besseli(2*phi,1)./phi./besseli(2*phi,1);
    elseif q==2
    tab=[ tab phi eta];
  elseif (order>1)


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

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')


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;
    xdot(2) = -q/t*x(2) + k*x(1)^order;


function [stp, isterminal, direction] = g(t,x)
  global order q k thiele

  isterminal = 1;
  direction = 0;