Figure 8.23:

Particle concentrations of A and B versus particle age for three different-sized particles.

Code for Figure 8.23

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
% 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.
%
% Revised 8/15/2018

p.k1  = 1e3;  %cm^3/mol min
p.kma = 0.01/60;  %cm/min
p.kmb = 0.01/60;  %cm/min
p.alpha = 1;
% p.r = 5; %cm
p.caf = 1e-3; %mol/cm^3
p.cbf = 1e-3; %mol/cm^3
p.thetabar = 10; %min

ncolpt = 20;
[col.z, col.A, col.B, col.Q] = colloc(ncolpt - 1,'left');
theta = col.z ./ (1 - col.z);
col.p = exp(-theta/p.thetabar) / p.thetabar;
col.Qz = col.Q./((1 - col.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

    p.r = rvec(i);
    x0 = [p.caf*ones(ncolpt,1); p.cbf*ones(ncolpt,1); p.caf; p.cbf];
    [x, fval, info] = fsolve(@(x) particles(x, ncolpt, p, col), 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 = col.Qz'*(ca .* col.p);
    intb = col.Qz'*(cb .* col.p);
    catot(i) = (inta*p.alpha + cabar(i))/(1 + p.alpha);
    cbtot(i) = (intb*p.alpha + cbbar(i))/(1 + p.alpha);

end


y0 = [p.caf; p.cbf];
opts = optimset ('MaxFunEvals', 2000*numel (x0), ...
                 'MaxIter', 500*numel (x0));
[y, fval, info] = fsolve(@(x) cstr(x, p), 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*p.caf*p.alpha/(1 + p.alpha)*ones(nrs,1) ];
rvec = [1e-4; 1e-3; 1e-2];
table = [theta];

for i = 1:length(rvec)

    p.r = rvec(i);
    x0 = [p.caf*ones(ncolpt,1); p.cbf*ones(ncolpt,1); p.caf; p.cbf];
    [x, fval, info] = fsolve(@(x) particles(x, ncolpt, p, col), 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 = col.Qz'*(ca .* col.p);
    intb = col.Qz'*(cb .* col.p);
    catot(i) = (inta*p.alpha + cabar(i))/(1 + p.alpha);
    cbtot(i) = (intb*p.alpha + cbbar(i))/(1 + p.alpha);
    table = [table 1000*ca 1000*cb];

end

table2 = [table col.z col.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
21
22
23
24
25
function rhs = particles(x, ncolpt, p, col)

    ca = x(1:ncolpt);
    cb = x(ncolpt+1:2*ncolpt);
    cabar = x(2*ncolpt + 1);
    cbbar = x(2*ncolpt + 2);

    inta = col.Qz'*(ca .* col.p);
    intb = col.Qz'*(cb .* col.p);
    rx =  p.k1 * ca.*cb;
    intrx = p.thetabar*col.Qz'*(rx .* col.p);
    first = [col.A*ca; col.A*cb];
    rhs = first - [p.kma*3/p.r*(cabar - ca) + rx;
                   p.kmb*3/p.r*(cbbar - cb) + rx];

    % initial condition
    rhs(1) = ca(1) - p.caf;
    rhs(ncolpt + 1) = cb(1);

    % total balances
    rhs(2*ncolpt+1) = p.caf - inta - intrx - ...
                      p.k1*cabar*cbbar*p.thetabar/p.alpha - ...
                      cabar/p.alpha;
    rhs(2*ncolpt+2) = p.cbf - p.alpha*intb - p.alpha*intrx -...
                      p.k1*cabar*cbbar*p.thetabar - cbbar;

cstr.m

1
2
3
4
5
6
7
function res = cstr(x, p)

    ca = x(1);
    cb = x(2);
    res= [p.caf*p.alpha/(1 + p.alpha) - ca - ...
                p.k1*ca*cb*p.thetabar;
          p.cbf/(1 + p.alpha) - cb - p.k1*ca*cb*p.thetabar];