Figure 8.37:

Yield of desired product B versus reactor length for different dispersion numbers.

Code for Figure 8.37

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
% 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 k1 k2 v D A B caf cbf ncolpt order

k1 = 1;
k2 = 1;
v =  1;
%D = 1000;
%D = 1;
caf = 1;
cbf = 0;
length = 0.5;
ncolpt = 50;
[z, A, B, Q] = colloc(ncolpt-2,'left','right');
z = z*length;
A = A/length;
B = B/(length*length);

Dvec = [1000; 1; 0.1; 1e-3];
nD = size(Dvec,1);
table(1:nD)  = {[]};
for i = 1:nD
  D=Dvec(i);
  x0=[caf*ones(ncolpt,1); caf*ones(ncolpt,1)];
  opts = optimset ('TolFun', 1e-16, 'TolX', 1e-10);
  [x, fval, info] = fsolve('pfrcol',x0,opts);
  fsolve_failed = info <= 0;
  if  (fsolve_failed)
    fprintf ('warning, fsolve failed, info = %d\n', info);
    i, D, x0
    fprintf('i= %g \n',i);
  end
  ca = x(1:ncolpt);
  cb = x(ncolpt+1:2*ncolpt);
  yield = cb ./ (caf - ca);
  conv  = (caf-ca) ./ (caf);
  table(i) = {[z conv yield]};
end

save dispersedpfr.dat table;

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
subplot (2, 1, 1);
hold on
for i = 1:nD
  plot (table{i}(:,1),table{i}(:,2))
end
% TITLE dispersedpfr_1

subplot (2, 1, 2);
hold on
for i = 1:nD
  plot (table{i}(:,1),table{i}(:,3))
end
hold off
% TITLE dispersedpfr_2
end % PLOTTING

pfrcol.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
function res = pfrcol(x)
  global v D A B k1 k2 caf cbf ncolpt order
  %
  % express the diff eq at every collocation point
  %
  ca = x(1:ncolpt);
  cb = x(ncolpt+1:2*ncolpt);
  Rab  = [-k1*ca; k1*ca - 2*k2*cb.^2];
  first  = [A*ca;A*cb];
  second = [B*ca;B*cb];
  res = -v*first + D*second + Rab;
  % write over the first and last points with the
  % boundary conditions
  res(1) = v*ca(1)-D*first(1)-v*caf;
  res(ncolpt) = first(ncolpt);
  res(ncolpt+1)= v*cb(1)-D*first(ncolpt+1)-v*cbf;
  res(2*ncolpt) = first(2*ncolpt);