Figure 7.32 (page 411):

Pellet CO profiles at several reactor positions.

Code for Figure 7.32

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.

%%
%% added log transformation in pellet
%% jbr
%% 12/13/03 

global  npts A Aint Bint Rint caf cbf ccf ...
        k10 E1 k20 E2 Ka0 Ea Kc0 Ec Da Db Dc T ...
        Ra Rb Rc Da Db Dc kma kmb kmc dcadr dcbdr dccdr ...
        epsb a T Ntfin Nafin Ncfin Pin Tin Qin Q Cptot U ...
        Ta Rt delH1 delH2 km Rp vis massf Ac xc Pext
        
Rg  = 8.314;  % J/K mol
Mair= 28.96; % g/mol, mol weight air
%Pin   = 1.5*1.013e5;  % N/m^2
Pin   = 2.0*1.013e5;  % N/m^2
Pext  = 1.013e5; % minimum allowed exit pressure
Tin   = 550; % K
Rt  = 10/2; %cm, tube radius
At  = pi*Rt*Rt; % cm^2, tube cross-sectional area
Ta  = 325; % K, ambient temperature
U    = 5.5e-3; % cal/(cm^2 K s), heat transfer coefficient
delH1 = (-67.63e3); % cal/mol CO;
delH2 = (-460.4e3); % cal/mol C_3H_6;
Cp  = 0.25; % cal/(g K), heat capacity of air
vis = 0.028e-2; % g/(cm s), viscosity of air
bt  = 1; % bed to tube area ratio
%u   = 500/bt;  % cm/sec, feed gas velocity entering bed
u  = 75/bt;  % cm/sec, feed gas velocity entering bed
Ac  = bt*At;  % cm^2 area of bed, 4 x area of tube
Qin = u*Ac;  % cm^3/sec, feed volumetric flowrate
P   = Pin;
T   = Tin;
Q   = Qin;
cfin= P/(Rg*T)*1e-6; % mol/cm^3
Rp  = 0.175; % cm radius of catalyst particle
a   = Rp/3;
rhob = 0.51; % g/cm^3 bed density
rhop = 0.68; % g/cm^3 particle density
epsb = 1 - rhob/rhop; % bed porosity, dimensionless
%epsb = 0.4;
xc  = 0.996;

cafin = cfin*0.02;
cbfin = cfin*0.03;
ccfin = cfin*5e-4;

Nafin = Q*cafin;
Nbfin = Q*cbfin;
Ncfin = Q*ccfin;
Nfin  = [Nafin; Nbfin; Ncfin];
Ntfin = Q*cfin;
massf = Ntfin*Mair; % mass flowrate, g/s, remains constant
Cptot = massf*Cp;

alpha = 1;
k10  = alpha*6.802e16*2.6e6*80/100*0.05/100;  %mol/cm^3 s
k20  = alpha*1.416e18*2.6e6*80/100*0.05/100;  %mol/cm^3 s
E1   = 13108; %K
E2   = 15109; %K
Ka0  = 8.099e6; % cm^3/mol
Kc0  = 2.579e8; % cm^3/mol
Ea   = - 409; %K
Ec   = 191; %K
Da   = 0.0487; % cm^2/s
Db   = 0.0469; % cm^2/s
Dc   = 0.0487; % cm^2/s
kma  = 0.4*9.76;   % cm/s
kmb  = 0.4*10.18;  % cm/s
kmc  = 0.4*9.76;   % cm/s

%%
%% calculate initial reaction rates
%%
k1      = k10*exp(-E1/Tin);
k2      = k20*exp(-E2/Tin);
Ka      = Ka0*exp(-Ea/Tin);
Kc      = Kc0*exp(-Ec/Tin);
den     = (1+Ka*cafin+Kc*ccfin)*(1+Ka*cafin+Kc*ccfin);
r1      = k1*cafin*cbfin/den;
r2      = k2*ccfin*cbfin/den;
%%
%% adiabatic temperature rise and 
%% Ucrit to achieve a zero intial temperature derivative
%%
adrise = -(cafin*delH1+ccfin*delH2)/(Cp*cfin*Mair)
Ucrit  = (r1*delH1 + r2*delH2)/(2/Rt*(Ta-T))

%%
%% global collocation 
%%
npts = 40;
[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 tube inlet
%%
caf = cafin; cbf=cbfin; ccf=ccfin;
zain = log(1e-9*caf); zaout = log(0.75*caf); 
%zain = log(1e-12*caf); zaout = log(0.9*caf); %P=1.5
za0 = zain + R/Rp*(zaout-zain);
zbin = log(0.75*cbf); zbout = log(cbf);
zb0 = zbin + R/Rp*(zbout-zbin);
zcin = log(1e-6*ccf); zcout = log(ccf);
%zcin = log(1e-9*ccf); zcout = log(0.9*ccf); %P=1.5
zc0 = zcin + R/Rp*(zcout-zcin);

z0 = [za0;zb0;zc0];
fsolve_options = optimset ('TolFun', 1e-16, 'TolX', 1e-10);
[z,fval,info] = fsolve('pellet',z0,fsolve_options);
fsolve_failed = info <= 0;

if (fsolve_failed) 
  plot(R,[z0(1:npts),z(1:npts),z0(npts+1:2*npts),z(npts+1:2*npts), ...
          z0(2*npts+1:3*npts),z(2*npts+1:3*npts)])
  error ('cannot find pellet profile at tube inlet') 
end

za = z(1:npts);
zb = z(npts+1:2*npts);
zc = z(2*npts+1:3*npts);

ca = exp(za);
cb = exp(zb);
cc = exp(zc);

%%
%% march down the bed
%%
nvs    = 201;
%vfinal = 4737;
vfinal = 2000;
vsteps = linspace (0,vfinal,nvs)';
vout   = vsteps;
y0     = [Nfin; Tin; Pin; z];
ydot0    = zeros(length(y0),1);
res      = bed(0,y0,ydot0);
ydot0(1:5)    = -res(1:5);

ymin = min(y0);
opts = odeset ('RelTol', 1e-10, 'AbsTol', 1e-10, 'Events', @stop);
t0=cputime();
[vout,y] = ode15i (@bed, vsteps, y0, ydot0, opts);
cpu = cputime()-t0
nout = length(vout);
if ( nout == nvs )
  fprintf ('hey, did not reach final conversion, increase vfinal\n');
end
xa = (Nafin-y(end,1))/Nafin
xb = (Nbfin-y(end,2))/Nbfin
xc = (Ncfin-y(end,3))/Ncfin
P=y(:,5);
T=y(:,4);
ctot = P./(Rg*T)*1e-6;
Nt = Ntfin - 1/2*(Nafin-y(:,1))+1/2*(Ncfin-y(:,3));
Q = Nt./ctot;
caf = y(:,1)./Q;
cbf = y(:,2)./Q;
ccf = y(:,3)./Q;
Patm = P/1.013e5;
table1 = [vout, caf, cbf, ccf, T, Patm];
%%
%% pick out some good length locations for pellet profiles
%%
nrows = 5;
%rowsp = ceil(linspace(1,nout,nrows));
rowsp = [1,5,10,50,90,100,nout]
vout(rowsp)
cols = [6:npts+5];
ca   = exp(y(rowsp,cols)');
cols = [npts+6:2*npts+5];
cb   = exp(y(rowsp,cols)');
cols = [2*npts+6:3*npts+5];
cc   = exp(y(rowsp,cols)');
table2 = [R, ca, cb, cc];

subplot (2, 2, 1);
semilogy (table1(:,1), table1(:,2:4));
title ('Figure 7.32')

subplot (2, 2, 2);
plot (table1(:,1), table1(:,5));
title ('Figure 7.32')

subplot (2, 2, 4);
plot (table1(:,1), table1(:,6));
title ('Figure 7.32')

subplot (2, 2, 3);
semilogy (table2(:,1), table2(:,2:8));
title ('Figure 7.32')

pellet.m

function retval = pellet(x)
  global  npts A Aint Bint Rint  caf cbf ccf ...
          k10 E1 k20 E2 Ka0 Ea Kc0 Ec Da Db Dc T ...
          Ra Rb Rc Da Db Dc kma kmb kmc dcadr dcbdr dccdr
  %%
  %% component A
  %%
  za = x(1:npts);
  zb = x(npts+1:2*npts);
  zc = x(2*npts+1:3*npts);
  ca = exp(za);
  cb = exp(zb);
  cc = exp(zc);

  k1      = k10*exp(-E1/T);
  k2      = k20*exp(-E2/T);
  Ka      = Ka0*exp(-Ea/T);
  Kc      = Kc0*exp(-Ec/T);
  den     = (1+Ka*ca+Kc*cc).^2;
  r1      = k1.*ca.*cb./den;
  r2      = k2*cc.*cb./den;
  Ra      = - r1;
  Rb      = - 1/2*r1  - 9/2*r2;
  Rc      = - r2;
  ip = 1;
  retval(ip,1)    = A(1,:)*za;
  caint = ca(2:npts-1);
  retval(ip+1:ip+npts-2,1) = Bint*za + Aint*za.*(Aint*za + 2./Rint) + ...
      Ra(2:npts-1)./(Da*caint);
  dzadr = A(npts,:)*za;
  retval(ip+npts-1,1) = Da*dzadr - kma*(caf/ca(npts) - 1);
  %%
  %% component B
  %%
  ip = npts+1;
  cbint = cb(2:npts-1);
  retval(ip,1)    = A(1,:)*zb;
  retval(ip+1:ip+npts-2,1) = Bint*zb + Aint*zb.*(Aint*zb + 2./Rint) + ...
      Rb(2:npts-1)./(Db*cbint);
  dzbdr = A(npts,:)*zb;
  retval(ip+npts-1,1) =  Db*dzbdr - kmb*(cbf/cb(npts) - 1);
  %%
  %% component C
  %%
  ip = 2*npts+1;
  ccint = cc(2:npts-1);
  retval(ip,1)    = A(1,:)*zc;
  retval(ip+1:ip+npts-2,1) = Bint*zc + Aint*zc.*(Aint*zc + 2./Rint) + ...
      Rc(2:npts-1)./(Dc*ccint);
  dzcdr = A(npts,:)*zc;
  retval(ip+npts-1,1) = Dc*dzcdr - kmc*(ccf/cc(npts) - 1);
  dcadr = A(npts,:)*ca;
  dcbdr = A(npts,:)*cb;
  dccdr = A(npts,:)*cc;

bed.m

function res = bed(t, y, ydot)
  global epsb Da Db Dc a T ...
      Ntfin Nafin Ncfin Pin Tin Qin Q Cptot U Ta Rt delH1 delH2 ...
      km Rp  vis massf Ac dcadr dcbdr dccdr caf cbf ccf
%  fprintf('time= %g \n',t);
  Naf = y(1);
  Nbf = y(2);
  Ncf = y(3);
  T   = y(4);
  P   = y(5);
  zpellet = y(6:length(y));

  Ntf  = Ntfin -1/2*(Nafin - Naf) + 1/2*(Ncfin - Ncf);
  Q  = Qin * Pin/P*T/Tin*Ntf/Ntfin;
  caf = Naf/Q;
  cbf = Nbf/Q;
  ccf = Ncf/Q;
  %%
  %% calculate pellet residual and update 
  %% total pellet reaction rate through dcadr dcbdr dccdr
  %%
  pelletres = pellet(zpellet);
  r1p   = Da/a*dcadr;
  r2p   = Dc/a*dccdr;

  Rj    = -(1-epsb)/a*[ Da*dcadr; Db*dcbdr; Dc*dccdr];
  dT    = ( -(delH1*r1p + delH2*r2p) + 2/Rt*U*(Ta-T) ) / Cptot;
%  dT    = 0;
  dP    = - (1-epsb)*Q / ( 2*Rp*epsb^3.*Ac^2. ) * ...
          ( 150*(1-epsb)*vis/(2*Rp) + 7/4*massf/Ac );
%  dP = 0;
  rhs = [Rj; dT; dP];
  res(1:5,1) = ydot(1:5) - rhs;
  res(6:length(y),1) = pelletres;

stop.m

function [retval, isterminal, direction] = stop(t,y,ydot)
  global xc Ncfin Pext
  Ncf = y(3);
  P   = y(5);
  convtest = xc - (1-Ncf/Ncfin);
  presstest = P - Pext;
  retval = min(convtest,presstest);
  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