Figure A.7 (page 656):

Molar flow of A versus reactor volume for second-order, isothermal reaction in a fixed-bed reactor; two approximations and exact solution.

Code for Figure A.7

Text of the GNU GPL.

main.m

%% 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 epsb a Nafin Ntfin Qf Da k1 A Bint Aint Rint npts Da caf xa
        
xa  = 0.75;
Rg  = 8.314;  % J/K mol
P   = 4*1.013e5;  % N/m^2
T   = 550; % K
ctf = P/(Rg*T)*1e-6; % mol/cm^3
Rp  = 0.45; % cm radius of catalyst particle
a   = Rp/3;
Nafin = 10; % mol/sec
NIf   = 10;
Ntf   = Nafin+NIf;
caf = ctf*0.5;
Qf  = Nafin/caf;
k1  = 2.25e5;  %cm^3/mol s
Da   = 0.008; % cm^2/s
rhob = 0.60; % g/cm^3 bed density
rhop = 0.68; % g/cm^3 particle density
epsb = 1 - rhob/rhop; % bed porosity, dimensionless
%%
%% global collocation 
%%
npts = 25;
[R A B Q] = colloc(npts-2, 'left', 'right');
R = R*Rp;
A = A/Rp;
B = B/(Rp*Rp);
Q = Q*Rp;
Aint = A(2:npts-1,:);
Bint = B(2:npts-1,:);
Rint = R(2:npts-1);


%%
%% find the pellet profile at bed inlet
%%

ca0 = logspace(log10(caf)-2,log10(caf),npts)';
tol = 1e-12;
opts = optimset ('TolFun', tol);
[ca,fval,info] = fsolve('pellet',ca0,opts);

info



%%
%% march down the bed
%%
nvs    = 100;
vfinal = 4e5;
vsteps = linspace (0,vfinal,nvs)';
vout   = vsteps;
y0     = [Nafin; ca];
ydot0    = zeros(length(y0),1);
res      = bed(y0,ydot0,0);
ydot0(1) = -res(1);

ymin = min(y0);
opts = odeset ('AbsTol', sqrt (eps), 'RelTol', 1e-10, 'Events', @stop);
[vout,y] = ode15i (@bed, vsteps, y0, ydot0, opts);
nout = length(vout);
if ( nout == nvs )
  fprintf ('hey, did not reach final conversion, increase vfinal\n');
end
xa = (Nafin-y(end,1))/Nafin;
Naf = y(:,1);
vplot = vout/1000.; %lit
VR = vplot(length(vplot));
table1 = [vplot, Naf];
%%
%% pick out some good length locations for pellet profiles
%%
rowsc = [1,nout];
colsc = [2:npts+1];
ca   = y(rowsc,colsc)';
table2 = [R, ca];
%save -ascii fb2colloc_2.dat table2;
Naout = (1-xa)*Nafin;
Natop = (1-xa+0.10)*Nafin;
dashedlines = [0,      Naout, VR, 0,     300,  Naout, VR, 2  ;
               1.1*VR, Naout, VR, Natop, 400,  Naout, VR, 3  ];

%%plot the molar flow versus reactor volume and 75% conversion line
plot(vplot, Naf, dashedlines(:,1), dashedlines(:,2), '-')
axis ([0, 400, 2, 10])
title ('Figure A.7')

pellet.m

function retval = pellet(x)
  global  k1 A Bint Aint Rint Da caf npts dcadr
  %%
  %% component A
  %%
  ca = x(1:npts);
  r1      = k1.*ca.*ca;
  Ra      = - r1;
  ip = 1;
  retval(ip,1)    = A(1,:)*ca;
  caint = ca(ip+1:ip+npts-2);
  retval(ip+1:ip+npts-2,1) = Bint*ca + 2*Aint*ca./Rint + ...
      Ra(2:npts-1)/Da;
  dcadr = A(npts,:)*ca;
%  retval(ip+npts-1) = Da*dcadr - kma*(caf - ca(npts));
  retval(ip+npts-1,1) = caf - ca(npts);

bed.m

function res = bed(t, y, ydot)
  global epsb a Ntfin Qf Da dcadr caf
%  fprintf('time= %g \n',t);
  Naf = y(1);
  cpellet = y(2:length(y));
  Q  = Qf;
  caf = Naf/Q;
  %%
  %% calculate pellet residual and update 
  %% total pellet reaction rate through dcadr
  %%
  pelletres = pellet(cpellet);
  r1p   = Da/a*dcadr;
  RA    = -(1-epsb)/a*[ Da*dcadr];
  res(1,1) = ydot(1) - RA;
  res(2:length(y),1) = pelletres;

stop.m

function [retval, isterminal, direction] = stop(t,y,ydot)
  global xa Nafin
  Naf = y(1);
  convtest = xa - (1-Naf/Nafin);
  retval = convtest;
  isterminal = 1;
  direction = 0;

../../util/matlab/colloc.m

function [r, a, b, q] = colloc (n, varargin)

  nargs = nargin;

  if (nargs < 1 || nargs > 3)
    error ('usage: [r, a, b, q] = colloc (n, ''left'', ''right'')');
  end

  if (~ (isnumeric (n) && numel (n) == 1 && round (n) == n))
    error ('colloc: first argument must be an integer scalar');
  end

  if (isnan (n) || isinf (n))
    error ('colloc: NaN is invalid as N');
  end

  if (n < 0)
    error ('colloc: first argument must be non-negative');
  end

  n0 = 0;
  n1 = 0;

  for i = 2:nargs
    s = varargin{i-1};
    if (~ ischar (s))
      error ('colloc: expecting character string argument');
    end

    s = lower (s);
    if (strcmp (s, 'r') || strcmp (s, 'right'))
      n1 = 1;
    elseif (strcmp (s, 'l') || strcmp (s, 'left'))
      n0 = 1;
    else
      error ('colloc: unrecognized argument');
    end
  end

  nt = n + n0 + n1;

  if (nt < 1)
    error ('colloc: the total number of roots must be positive');
  end

  alpha = 0;
  beta = 0;

  %% Compute roots.

  [dif1, dif2, dif3, r] = jcobi (n, n0, n1, alpha, beta);

  %% First derivative weights.

  a = zeros (nt, nt);
  for i = 1:nt
    a(i,:) = dfopr (n, n0, n1, i, 1, dif1, dif2, dif3, r)';
  end

  %% Second derivative weights.

  b = zeros (nt, nt);
  for i = 1:nt
    b(i,:) = dfopr (n, n0, n1, i, 2, dif1, dif2, dif3, r)';
  end

  %% Gaussian quadrature weights.

  id = 3;
  q = dfopr (n, n0, n1, 0, id, dif1, dif2, dif3, r);

end


%% The following routines (JCOBI, DIF, DFOPR, INTRP, AND RADAU)
%% are the same as found in Villadsen, J. and M.L. Michelsen,
%% Solution of Differential Equation Models by Polynomial
%% Approximation, Prentice-Hall (1978) pages 418-420.
%%
%% Cosmetic changes (elimination of arithmetic IF statements, most
%% GO TO statements, and indentation of program blocks) made by:
%%
%% John W. Eaton
%% Department of Chemical Engineering
%% The University of Texas at Austin
%% Austin, Texas 78712
%%
%% June 6, 1987
%%
%% Some error checking additions also made on June 7, 1987
%%
%% Further cosmetic changes made August 20, 1987
%%
%% Translated from Fortran December 14, 2006

function vect = dfopr (n, n0, n1, i, id, dif1, dif2, dif3, root)

%%   Villadsen and Michelsen, pages 133-134, 419
%%
%%   Input parameters:
%%
%%     N      : The degree of the Jacobi polynomial, (i.e. the number
%%              of interior interpolation points)
%%
%%     N0     : Determines whether x = 0 is included as an
%%              interpolation point
%%
%%                n0 = 0  ==>  x = 0 is not included
%%                n0 = 1  ==>  x = 0 is included
%%
%%     N1     : Determines whether x = 1 is included as an
%%              interpolation point
%%
%%                n1 = 0  ==>  x = 1 is not included
%%                n1 = 1  ==>  x = 1 is included
%%
%%     I      : The index of the node for which the weights are to be
%%              calculated
%%
%%     ID     : Indicator
%%
%%                id = 1  ==>  first derivative weights are computed
%%                id = 2  ==>  second derivative weights are computed
%%                id = 3  ==>  Gaussian weights are computed (in this
%%                             case, I is not used).
%%
%%   Output parameters:
%%
%%     DIF1   : vector containing the first derivative
%%              of the node polynomial at the zeros
%%
%%     DIF2   : vector containing the second derivative
%%              of the node polynomial at the zeros
%%
%%     DIF3   : vector containing the third derivative
%%              of the node polynomial at the zeros
%%
%%     VECT   : vector of computed weights

  if (n0 ~= 0 && n0 ~= 1)
    error ('dfopr: n0 not equal to 0 or 1');
  end

  if (n1 ~= 0 && n1 ~= 1)
    error ('dfopr: n1 not equal to 0 or 1');
  end

  if (n < 0)
    error ('dfopr: n less than 0');
  end

  nt = n + n0 + n1;

  if (id ~= 1 && id ~= 2 && id ~= 3)
    error ('dfopr: id not equal to 1, 2, or 3');
  end

  if (id ~= 3)
    if (i < 1)
      error ('dfopr: index less than zero');
    end

    if (i > nt)
      error ('dfopr: index greater than number of interpolation points');
    end
  end

  if (nt < 1)
    error ('dfopr: number of interpolation points less than 1');
  end

%% Evaluate discretization matrices and Gaussian quadrature
%% weights.  Quadrature weights are normalized to sum to one.

  vect = zeros (nt, 1);

  if (id ~= 3)
    for j = 1:nt
      if (j == i)
	if (id == 1)
	  vect(i) = dif2(i)/dif1(i)/2.0;
	else
	  vect(i) = dif3(i)/dif1(i)/3.0;
	end
      else
	y = root(i) - root(j);
	vect(j) = dif1(i)/dif1(j)/y;
	if (id == 2)
	  vect(j) = vect(j)*(dif2(i)/dif1(i) - 2.0/y);
	end
      end
    end
  else
    y = 0.0;

    for j = 1:nt

      x = root(j);
      ax = x*(1.0 - x);

      if (n0 == 0)
	ax = ax/x/x;
      end

      if (n1 == 0)
	ax = ax/(1.0 - x)/(1.0 - x);
      end

      vect(j) = ax/dif1(j)^2;
      y = y + vect(j);

    end

    vect = vect/y;

  end

end

function [dif1, dif2, dif3, root] = jcobi (n, n0, n1, alpha, beta)

%%   Villadsen and Michelsen, pages 131-132, 418
%%
%%   This subroutine computes the zeros of the Jacobi polynomial
%%
%%      (ALPHA,BETA)
%%     P  (X)
%%      N
%%
%%   Input parameters:
%%
%%     N      : The degree of the Jacobi polynomial, (i.e. the number
%%              of interior interpolation points)
%%
%%     N0     : Determines whether x = 0 is included as an
%%              interpolation point
%%
%%                N0 = 0  ==>  x = 0 is not included
%%                N0 = 1  ==>  x = 0 is included
%%
%%     N1     : Determines whether x = 1 is included as an
%%              interpolation point
%%
%%                N1 = 0  ==>  x = 1 is not included
%%                N1 = 1  ==>  x = 1 is included
%%
%%     ALPHA  : The value of alpha in the description of the jacobi
%%              polynomial
%%
%%     BETA   : The value of beta in the description of the Jacobi
%%              polynomial
%%
%%     For a more complete explanation of alpha and beta, see Villadsen
%%     and Michelsen, pages 57 to 59
%%
%%   Output parameters:
%%
%%     ROOT   : vector containing the n + n0 + n1 zeros of the node
%%              polynomial used in the interpolation routine
%%
%%     DIF1   : vector containing the first derivative
%%              of the node polynomial at the zeros
%%
%%     DIF2   : vector containing the second derivative
%%              of the node polynomial at the zeros
%%
%%     DIF3   : vector containing the third derivative
%%              of the node polynomial at the zeros

  if (n0 ~= 0 && n0 ~= 1)
    error ('jcobi: n0 not equal to 0 or 1');
  end

  if (n1 ~= 0 && n1 ~= 1)
    error ('jcobi: n1 not equal to 0 or 1');
  end

  if (n < 0)
    error ('jcobi: n less than 0');
  end

  nt = n + n0 + n1;

  if (nt < 1)
    error ('jcobi: number of interpolation points less than 1');
  end

  dif1 = zeros (nt, 1);
  dif2 = zeros (nt, 1);
  dif3 = zeros (nt, 1);
  root = zeros (nt, 1);

%% First evaluation of coefficients in recursion formulas.
%% recursion coefficients are stored in dif1 and dif2.

  ab = alpha + beta;
  ad = beta - alpha;
  ap = beta*alpha;
  dif1(1) = (ad/(ab + 2.0) + 1.0)/2.0;
  dif2(1) = 0.0;

  if (n >= 2)
    for i = 2:n

      z1 = i - 1.0;
      z = ab + 2*z1;
      dif1(i) = (ab*ad/z/(z + 2.0) + 1.0)/2.0;

      if (i == 2)
	dif2(i) = (ab + ap + z1)/z/z/(z + 1.0);
      else
	z = z*z;
	y = z1*(ab + z1);
	y = y*(ap + y);
	dif2(i) = y/z/(z - 1.0);
      end

    end
  end

%% Root determination by newton method with suppression of
%% previously determined roots.

  x = 0.0;

  for i = 1:n

    done = false;

    while (~ done)

      xd = 0.0;
      xn = 1.0;
      xd1 = 0.0;
      xn1 = 0.0;

      for j = 1:n
	xp = (dif1(j) - x)*xn  - dif2(j)*xd;
	xp1 = (dif1(j) - x)*xn1 - dif2(j)*xd1 - xn;
	xd = xn;
	xd1 = xn1;
	xn = xp;
	xn1 = xp1;
      end

      zc = 1.0;
      z = xn/xn1;

      if (i ~= 1)
	for j = 2:i
	  zc = zc - z/(x - root(j-1));
	end
      end

      z = z/zc;
      x = x - z;

      if (abs(z) <= 1.0e-09)
	done = true;
      end

    end

    root(i) = x;
    x = x + 0.0001;

  end

%% Add interpolation points at x = 0 and/or x = 1.

  nt = n + n0 + n1;

  if (n0 ~= 0)
    for i = 1:n
      j = n + 1 - i;
      root(j+1) = root(j);
    end
    root(1) = 0.0;
  end

  if (n1 == 1)
    root(nt) = 1.0;
  end


%% Use recursion formulas to evaluate derivatives of node polynomial
%%
%%                   N0     (ALPHA,BETA)           N1
%%     P  (X)  =  (X)   *  P (X)         *  (1 - X)
%%      NT                   N
%%
%% at the interpolation points.

  for i = 1:nt
    x = root(i);
    dif1(i) = 1.0;
    dif2(i) = 0.0;
    dif3(i) = 0.0;
    for j = 1:nt
      if (j ~= i)
	y = x - root(j);
	dif3(i) = y*dif3(i) + 3.0*dif2(i);
	dif2(i) = y*dif2(i) + 2.0*dif1(i);
	dif1(i) = y*dif1(i);
      end
    end
  end

end

Data files

../../ch7/figures/fborder2.dat

 0.00000000e+00 1.00000000e+01
 5.05050505e+00 9.71822838e+00
 1.01010101e+01 9.44840658e+00
 1.51515152e+01 9.18986553e+00
 2.02020202e+01 8.94198244e+00
 2.52525253e+01 8.70417701e+00
 3.03030303e+01 8.47590800e+00
 3.53535354e+01 8.25667005e+00
 4.04040404e+01 8.04599093e+00
 4.54545455e+01 7.84342883e+00
 5.05050505e+01 7.64857015e+00
 5.55555556e+01 7.46102741e+00
 6.06060606e+01 7.28043723e+00
 6.56565657e+01 7.10645855e+00
 7.07070707e+01 6.93877097e+00
 7.57575758e+01 6.77707337e+00
 8.08080808e+01 6.62108242e+00
 8.58585859e+01 6.47053140e+00
 9.09090909e+01 6.32516900e+00
 9.59595960e+01 6.18475829e+00
 1.01010101e+02 6.04907579e+00
 1.06060606e+02 5.91791043e+00
 1.11111111e+02 5.79106286e+00
 1.16161616e+02 5.66834463e+00
 1.21212121e+02 5.54957749e+00
 1.26262626e+02 5.43459269e+00
 1.31313131e+02 5.32323043e+00
 1.36363636e+02 5.21533930e+00
 1.41414141e+02 5.11077572e+00
 1.46464646e+02 5.00940346e+00
 1.51515152e+02 4.91109322e+00
 1.56565657e+02 4.81572222e+00
 1.61616162e+02 4.72317377e+00
 1.66666667e+02 4.63333689e+00
 1.71717172e+02 4.54610606e+00
 1.76767677e+02 4.46138086e+00
 1.81818182e+02 4.37906566e+00
 1.86868687e+02 4.29906932e+00
 1.91919192e+02 4.22130504e+00
 1.96969697e+02 4.14569003e+00
 2.02020202e+02 4.07214533e+00
 2.07070707e+02 4.00059554e+00
 2.12121212e+02 3.93096869e+00
 2.17171717e+02 3.86319603e+00
 2.22222222e+02 3.79721189e+00
 2.27272727e+02 3.73295347e+00
 2.32323232e+02 3.67036073e+00
 2.37373737e+02 3.60937617e+00
 2.42424242e+02 3.54994481e+00
 2.47474747e+02 3.49201398e+00
 2.52525253e+02 3.43553327e+00
 2.57575758e+02 3.38045434e+00
 2.62626263e+02 3.32673084e+00
 2.67676768e+02 3.27431834e+00
 2.72727273e+02 3.22317423e+00
 2.77777778e+02 3.17325763e+00
 2.82828283e+02 3.12452928e+00
 2.87878788e+02 3.07695149e+00
 2.92929293e+02 3.03048805e+00
 2.97979798e+02 2.98510417e+00
 3.03030303e+02 2.94076641e+00
 3.08080808e+02 2.89744262e+00
 3.13131313e+02 2.85510189e+00
 3.18181818e+02 2.81371446e+00
 3.23232323e+02 2.77325172e+00
 3.28282828e+02 2.73368611e+00
 3.33333333e+02 2.69499109e+00
 3.38383838e+02 2.65714111e+00
 3.43434343e+02 2.62011154e+00
 3.48484848e+02 2.58387865e+00
 3.53535354e+02 2.54841957e+00
 3.58585859e+02 2.51371224e+00
 3.60611253e+02 2.50000000e+00

../../ch7/figures/fborder2_asy.dat

 0.00000000e+00 1.00000000e+01
 5.05050505e+00 9.70318967e+00
 1.01010101e+01 9.41940067e+00
 1.51515152e+01 9.14788197e+00
 2.02020202e+01 8.88793635e+00
 2.52525253e+01 8.63891517e+00
 3.03030303e+01 8.40021482e+00
 3.53535354e+01 8.17127274e+00
 4.04040404e+01 7.95156421e+00
 4.54545455e+01 7.74059928e+00
 5.05050505e+01 7.53792001e+00
 5.55555556e+01 7.34309817e+00
 6.06060606e+01 7.15573273e+00
 6.56565657e+01 6.97544809e+00
 7.07070707e+01 6.80189187e+00
 7.57575758e+01 6.63473337e+00
 8.08080808e+01 6.47366194e+00
 8.58585859e+01 6.31838563e+00
 9.09090909e+01 6.16862972e+00
 9.59595960e+01 6.02413559e+00
 1.01010101e+02 5.88465960e+00
 1.06060606e+02 5.74997198e+00
 1.11111111e+02 5.61985604e+00
 1.16161616e+02 5.49410725e+00
 1.21212121e+02 5.37253232e+00
 1.26262626e+02 5.25494854e+00
 1.31313131e+02 5.14118313e+00
 1.36363636e+02 5.03107250e+00
 1.41414141e+02 4.92446176e+00
 1.46464646e+02 4.82120417e+00
 1.51515152e+02 4.72116056e+00
 1.56565657e+02 4.62419893e+00
 1.61616162e+02 4.53019396e+00
 1.66666667e+02 4.43902665e+00
 1.71717172e+02 4.35058393e+00
 1.76767677e+02 4.26475830e+00
 1.81818182e+02 4.18144751e+00
 1.86868687e+02 4.10055427e+00
 1.91919192e+02 4.02198593e+00
 1.96969697e+02 3.94565426e+00
 2.02020202e+02 3.87147514e+00
 2.07070707e+02 3.79936840e+00
 2.12121212e+02 3.72925755e+00
 2.17171717e+02 3.66106962e+00
 2.22222222e+02 3.59473490e+00
 2.27272727e+02 3.53018685e+00
 2.32323232e+02 3.46736188e+00
 2.37373737e+02 3.40619921e+00
 2.42424242e+02 3.34664071e+00
 2.47474747e+02 3.28863076e+00
 2.52525253e+02 3.23211615e+00
 2.57575758e+02 3.17704590e+00
 2.62626263e+02 3.12337123e+00
 2.67676768e+02 3.07104536e+00
 2.72727273e+02 3.02002349e+00
 2.77777778e+02 2.97026264e+00
 2.82828283e+02 2.92172160e+00
 2.87878788e+02 2.87436082e+00
 2.92929293e+02 2.82814234e+00
 2.97979798e+02 2.78302972e+00
 3.03030303e+02 2.73898797e+00
 3.08080808e+02 2.69598345e+00
 3.13131313e+02 2.65398385e+00
 3.18181818e+02 2.61295811e+00
 3.23232323e+02 2.57287633e+00
 3.28282828e+02 2.53370979e+00
 3.32724389e+02 2.50000000e+00

../../ch7/figures/fborder2_dashed.dat

 0.00000000e+00 2.50000000e+00 3.32724389e+02 0.00000000e+00 3.60611253e+02 0.00000000e+00 3.60611253e+02 2.00000000e+00 3.32724389e+02 2.00000000e+00
 3.96672378e+02 2.50000000e+00 3.32724389e+02 3.50000000e+00 3.60611253e+02 3.50000000e+00 3.60611253e+02 3.00000000e+00 3.32724389e+02 3.00000000e+00