Figure 7.9:

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

Code for Figure 7.9

Text of the GNU GPL.

main.m


  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
% 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=[1; 2; 5; 10];
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;
%

%
% ODE's for order<1


% Stopping Critierion for order>2


% Main program

no = length(ordervec);

for j=1:no
  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

    tsteps=[linspace(start,start+2,30)';logspace(log10(start+2),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 =[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=[phi eta];
%    save ('-ascii', sprintf ('etaordern_large_%d.dat', j), 'tab');

  elseif (order>1)
    x_0=[c0;0];
    y=k*c0^order;
    xdot_0=[0;y];
    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=[ phi(2:end,:) eta(2:end,:)];
%    save ('-ascii', sprintf ('etaordern_large_%d.dat', j), 'tab');
  end
  table(j) = {tab};

  % if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
  % hold ('on');
  % loglog (tab(:,1), tab(:,2));
  % % TITLE
  % end % PLOTTING

end

%hold ('off');

save etaordern_large.dat table

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
  hold on;
  for i = 1:no
    loglog (table{i}(:,1), table{i}(:,2));
  end
  hold off;
  % TITLE
end % PLOTTING

ivp.m


1
2
3
4
5
6
7
8
9
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


1
2
3
4
5
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;