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 | % 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 % % streamlined code, jbr, 12/29/2017 % 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 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 % 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. % This derivative program cvsrnsmall.m calculates the concentration % profile for n<1 to show that an inner sphere of the pellet has % exactly zero concentration of A. That is a result of the order of the % reaction being less than one. global order k q thiele % Choose the orders n, 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. order=0; q=2; % % The thiele values are chosen to so that for given order and geometry, % first one has nonzero center concentration, second one has zero % concentration right at center, and remaining ones have zero % conentration at finite radius. % thielevec=[0.4; 0.57735; 0.8; 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-8; start=1; tab = []; % % Stopping Critierion % % % ODE's % nt = length(thielevec); table(1:nt) = {[]}; opts = odeset ('Events', @g, 'AbsTol', sqrt (eps), 'RelTol', sqrt (eps)); tout=[0;logspace(-2,5,400)']; t0=start; nfine = 200; ncoarse = 100; fine = logspace(log10(start),log10(start+1),nfine)'; coarse = logspace(log10(start+2),log10(40),ncoarse)'; tout2 = [fine; coarse(2:ncoarse)]; % Main program for i = 1:nt thiele = thielevec(i); if (i <= 2) % Stage 1 x0=[c0; 0]; tpts = tout; else % Stage 2 x0=[0; der]; tpts = tout2; end [tsteps,x] = ode15s (@f, tpts, x0, opts); tabletmp = [ tsteps./tsteps(end)*(q+1), x(:,1)./x(end,1) ]; table(i) = {tabletmp}; end save cvsrnsmall.dat table if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING hold on; for i = 1:nt plot (table{i}(:,1), table{i}(:,2)); end hold off; % TITLE end % PLOTTING |

1 2 3 4 5 6 7 8 9 10 11 12 | function [stp, isterminal, direction] = g(t,x) global order k q thiele if (x(1) == 0) pelthiele = Inf; else pelthiele = sqrt((order+1)/2*x(1)^(order-1)*k)*t/(q+1); end % stop at each desired thiele modulus stp = pelthiele-thiele; isterminal = 1; direction = 0; |

1 2 3 4 5 6 7 8 9 | function xdot=f(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 |