Figure 8.22:

Total concentration of A in the reactor effluent versus particle size.

Code for Figure 8.22

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

global thetabar caf cbf alpha kma kmb k1 ncolpt A B Qz z r p cabar cbbar
k1  = 1e3;  %cm^3/mol min
kma = 0.01/60;  %cm/min
kmb = 0.01/60;  %cm/min
alpha = 1;
r = 5; %cm
caf = 1e-3; %mol/cm^3
cbf = 1e-3; %mol/cm^3
thetabar = 10; %min

ncolpt = 20;
[z, A, B, Q] = colloc(ncolpt-1,'left');
theta = z ./ (1 -z);
p = exp(-theta/thetabar) / thetabar;
Qz = Q./((1 -z) .^ 2);

r0=1;
rfin = 1e-5;
nrs = 50;
rvec = logspace(log10(r0), log10(rfin), nrs)';
cabar = zeros (nrs, 1);
cbbar = zeros (nrs, 1);
catot = zeros (nrs, 1);
cbtot = zeros (nrs, 1);
for i=1:nrs
  r = rvec(i);
  x0 = [caf*ones(ncolpt,1); cbf*ones(ncolpt,1); caf; cbf];
  [x, fval, info] = fsolve('particles',x0);
  fsolve_failed = info <= 0;
  if (fsolve_failed)
    fprintf ('warning, fsolve failed\n');
    info, x0
  end
  ca = x(1:ncolpt);
  cb = x(ncolpt+1:2*ncolpt);
  cabar(i) = x(2*ncolpt+1);
  cbbar(i) = x(2*ncolpt+2);
  inta = Qz'*(ca .* p);
  intb = Qz'*(cb .* p);
  catot(i) = (inta*alpha + cabar(i))/(1+alpha);
  cbtot(i) = (intb*alpha + cbbar(i))/(1+alpha);
end

y0 = [caf; cbf];
opts = optimset ('MaxFunEvals', 2000*numel (x0), ...
                 'MaxIter', 500*numel (x0));
[y, fval, info] = fsolve('cstr',y0,opts);
fsolve_failed = info <= 0;
if (fsolve_failed)
  fprintf ('warning, fsolve failed\n');
  info, y0
end
cacstr = y(1);
cbcstr = y(2);

table1 = [rvec*1e4 1000*catot 1000*cacstr*ones(nrs,1) 1000*caf*alpha/(1+alpha)*ones(nrs,1) ];


rvec = [1e-4; 1e-3; 1e-2];
table = [theta];
for i=1:length(rvec)
  r = rvec(i);
  x0 = [caf*ones(ncolpt,1); cbf*ones(ncolpt,1); caf; cbf];
  [x, fval, info] = fsolve('particles',x0);
  fsolve_failed = info <= 0;
  if (fsolve_failed)
    fprintf ('warning, fsolve failed\n');
    info, x0
  end
  ca = x(1:ncolpt);
  cb = x(ncolpt+1:2*ncolpt);
  cabar(i) = x(2*ncolpt+1);
  cbbar(i) = x(2*ncolpt+2);
  inta = Qz'*(ca .* p);
  intb = Qz'*(cb .* p);
  catot(i) = (inta*alpha + cabar(i))/(1+alpha);
  cbtot(i) = (intb*alpha + cbbar(i))/(1+alpha);
  table = [table 1000*ca 1000*cb];
end
table2 = [table z p];

save suspension.dat table1 table2;

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
subplot (2, 1, 1);
semilogx (table1(:,1), table1(:,2:4));
% TITLE suspension_1

subplot (2, 1, 2);
plot (table2(:,1), table2(:,2:7));
% TITLE suspension_2
end % PLOTTING

particles.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
function rhs = particles(x)
  global thetabar caf cbf alpha kma kmb k1 ncolpt A B Qz z r p cabar cbbar
  ca = x(1:ncolpt);
  cb = x(ncolpt+1:2*ncolpt);
  cabar = x(2*ncolpt+1);
  cbbar = x(2*ncolpt+2);
  inta = Qz'*(ca .* p);
  intb = Qz'*(cb .* p);
  rx =  k1 * ca.*cb;
  intrx = thetabar*Qz'*(rx .* p);
  first = [A*ca; A*cb];
  rhs = first - [kma*3/r*(cabar-ca)+rx; kmb*3/r*(cbbar-cb)+rx];
  % initial condition
  rhs(1) = ca(1)-caf;
  rhs(ncolpt+1) = cb(1);
  % total balances
  rhs(2*ncolpt+1) = caf - inta - intrx -k1*cabar*cbbar*thetabar/alpha ...
                    - cabar/alpha;
  rhs(2*ncolpt+2) = cbf - alpha*intb - alpha*intrx -k1*cabar*cbbar*thetabar ...
                    - cbbar;

cstr.m


1
2
3
4
5
6
function res = cstr(x)
  global caf cbf thetabar k1  alpha
  ca = x(1);
  cb = x(2);
  res= [caf*alpha/(1+alpha) - ca - k1*ca*cb*thetabar;
        cbf/(1+alpha) - cb - k1*ca*cb*thetabar];