Computational Methods

Copyright c 2015 by Nob Hill Publishing, LLC

• We summarize many of the fundamental computational procedures required

in reactor analysis and design.

• We illustrate these procedures using the numerical computing language Oc-

tave [4, 5] or, equivalently, MATLAB.

• Octave is freely available for a variety of hardware platforms and can be down-
loaded from [Link].

• The MATLAB package is commercially available from The MathWorks, Inc., and
is becoming a commonly available tool of industrial engineering practice.

1
 Software

• Many of the figures and tables in this text required numerical solutions of
various types of chemical reactor problems.

• All computations required to produce every figure and table in this text were
performed with Octave.

• The Octave files required to produce these figures and tables are freely avail-
able at these two websites:

[Link]
[Link]/~jbraw/chemreacfun

2
 Linear Algebra and Least Squares

• After starting an Octave session, type help -i introduction at the Octave

prompt for information on how to enter matrices and perform matrix multi-
plication.

• Consider again the water gas shift reaction chemistry

 
H
 H2
   
 
0 1 0 −1 −1 1    0 


 −1
 OH
1 1 −1 0 0  = 0 
  
  H2 O  
1 0 −1 0 −1 1   0

 CO 

CO2

3
 Linear Indepdendence

• Say we would like to find the number of linearly independent reactions,

and then compute the species production rates from given reaction rates,
[r1 r2 r3] = [1 2 3].

4
 Linear Indepdendence

• A typical set of Octave commands to perform these operations would appear

as follows
octave:1> stoi=[0 1 0 -1 -1 1; -1 1 1 -1 0 0; \
1 0 -1 0 -1 1];
octave:2> rank(stoi)
ans = 2
octave:3> r=[1;2;3];
octave:4> R=stoi’*r
R =
1
3
-1
-3
-4
4

5
 Least Squares Estimation of r

• As we explored in Chapters 2 and 9, from measured R, we can solve for r as

a least-squares problem.

• Given the production rates and reaction rates are related by

R = νT r (1)

• The least-squares solution is given by

r = (ν νT )−1νR

in which the superscript −1 indicates a matrix inverse.

6
• The command in Octave and MATLAB for the least-squares solution is

r = stoi’ \ R

• This formula can be recalled with the mnemonic of a “left division” of Equa-
tion A.1 by νT .

7
 Estimating reaction rates

Example 0.1: Estimating reaction rates

• Consider the first and second water gas shift reaction as an independent set
and let [r1 r2] = [1 2].

• Create 2000 production rate measurements by adding noise to R = νT r and

estimate the reaction rates from these data.

• Plot the distribution of estimated reaction rates.

8
The Distribution of Estimated Reaction Rates

2.1

2.05

r2 2

1.95

1.9
0.9 0.95 1 1.05 1.1
r1

9
 The Octave Code

Solution

• The Octave commands to generate these results are

stoi=[0 1 0 -1 -1 1; -1 1 1 -1 0 0];
r=[1;2];
R=stoi’*r;
npoints = 2000;
for i=1:npoints R_meas(:,i)=0.05*randn(6,1)+R; endfor
r_est=stoi’ \ R_meas;

• The figure shows the estimates.

• We know from Chapter 9 that the distribution of estimates is a multivariate

normal.

10
• We also know how to calculate the α-level confidence intervals.

11
 Nonlinear Algebraic Equations and Optimization

• Determining equilibria for reacting systems with multiple phases and reac-
tions can require significant computational effort.

• As we saw in Chapter 3, the phase and reaction equilibrium conditions gener-

ally lead to mathematical problems of two types: solving nonlinear algebraic
equations and minimizing a nonlinear function subject to constraints.

• In this section it is assumed that the required thermochemical data are avail-
able, but finding or measuring these data is often another significant chal-
lenge in computing equilibria for systems of industrial interest.

• See also Section 3.2.1 for a brief discussion of thermochemical databases.

12
 Defining Functions (function)

• Octave and MATLAB provide many convenient built-in functions.

• Some that we have used in this text include the matrix exponential, expm; in-
complete gamma function, gammai; pseudorandom number generators, rand
and randm; and several others.

• Type help -i and select the menu item Function Index for a complete
list of Octave’s built-in functions.

13
 Defining Functions (function)

• But often we need to define our own functions.

• To solve algebraic equations or optimization problems, for example, the user

needs to provide a function that can be called by the solver to evaluate the
nonlinear equations or the objective function.

• To define a function in Octave, f (x) = x 2 for example, the command struc-

ture is

function y = f(x)
y = x*x;
endfunction

14
 Defining Functions (function)

• The first line states that the value of the vari-

able y is returned by invoking the function
function y = f(x)
named f with the input value x.
y = x*x;
endfunction
• The body of the function, then, is the calcu-
lation of the value of y to be returned.

15
 The m-files

• As the computational complexity of the function increases, we would like to

store our work so that it does not need to be reentered at the beginning of
each subsequent Octave session.

• Octave commands and function definitions can be stored in text files, called
m-files

• The naming convention is filename.m.

• Editing these files with a text editor during an interactive Octave session pro-
vides a convenient means for debugging new calculations.

16
 Solving Nonlinear Algebraic Equations (fsolve)

• We consider again the problem in Example 3.5. The equilibrium condition

requires the solution of two equations in two unknowns,

P K1yI yB − yP1 = 0
P K2yI yB − yP2 = 0

• yI , yB , yP1 , yP2 are defined in terms of the two reaction extents in Equa-
tions 3.66.

17
 Function Definition

• An Octave function defining these two equations is

function residual = dgdx(x)

K1 = 108; K2 = 284; P = 2.5;
yI0 = 0.5; yB0 = 0.5; yP10 = 0; yP20 = 0;
d = 1 - x(1) - x(2);
yI = (yI0 - x(1) - x(2)) / d;
yB = (yB0 - x(1) - x(2)) / d;
yP1 = (yP10 + x(1)) / d;
yP2 = (yP20 + x(2)) / d;
residual(1) = P*K1*yI*yB - yP1;
residual(2) = P*K2*yI*yB - yP2;
endfunction

18
 Using the Function

• Notice that x is a vector containing the two reaction extents provided to the
function named dgdx when it is invoked, and residual is a vector containing
the values of the two equilibrium conditions to be returned by the function.

• We seek values of x such that residual is zero. With this function defined, we
can use the built-in nonlinear equation solver to find its solution, by entering

[x,info] = fsolve("dgdx", x0)

19
 fsolve

• The fsolve function requires the name of the function, inside quotes, and an
initial guess for the unknown extents, provided in the variable x0, and returns
the solution in x and a flag info indicating if the calculation was successful.

• It is highly recommended to examine these information flags after all calcu-

lations.

• Reactor equilibrium problems can be numerically challenging, and even the

best software can run into problems.

20
 Solving the Problem

• After the function dgdx is defined, the following is a typical session to com-
pute the solution given in Example 3.5

octave:1> x0=[0.2;0.2];
octave:2> [x,info] = fsolve("dgdx",x0)
x =
0.13317
0.35087
info = 1

• The value of info = 1 indicates that the solution has converged.

21
 Solving Nonlinear Optimization Problems (npsol)

• The other main approach to finding the reaction equilibrium is to minimize

the appropriate energy function, in this case the Gibbs energy.

• This optimization-based formulation of the problem, as shown in Exam-

ple 3.3, can be more informative than the algebraic approach discussed
above.

22
 Gibbs Energy Function

• In Chapter 3, we derived the following Gibbs energy function for an ideal-gas

mixture
 
X X
G̃ = − εi0 ln Ki + 1 + ν̄iεi0 ln P
i i
 
0
" P #
X
yj0 +
X
0 yj0 + i νij εi
+ νij εi ln
1 + i ν̄iεi0
P
j i

• The minimization of this function of εi0 then determines the two equilibrium
extents.

23
 Minimization of G̃

• The figure shows the lines of constant Gibbs energy as a function of the two
reaction extents.

24
 0.5
-2.557
0.45 -2.55
-2.53
0.4 -2.5
-2
0.35 -1
0.3 0
0 0.25
ε2
0.2
0.15
0.1
0.05
0
0 0.1 0.2 0.3 0.4 0.5
0
ε1

25
 Minimization of G̃

• We see immediately that the minimum is unique.

• Notice that G̃ is not defined for all reaction extents, and the following con-
straints are required to ensure nonnegative concentrations

0 ≤ ε10 0 ≤ ε20 ε10 + ε20 ≤ 0.5

26
 Octave Code

• First we define a function that evaluates G̃ given the two reaction extents.
function retval=gibbs(x)
K1 = 108; K2 = 284; P = 2.5;
yI0 = 0.5; yB0 = 0.5; yP10 = 0; yP20 = 0;
d = 1 - x(1) - x(2);
yI = (yI0 - x(1) - x(2)) / d;
yB = (yB0 - x(1) - x(2)) / d;
yP1 = (yP10 + x(1)) / d;
yP2 = (yP20 + x(2)) / d;
retval = - (x(1)*log(K1)+x(2)*log(K2)) + \
(1-x(1)-x(2))*log(P) + yI*d*log(yI) + \
yB*d*log(yB) + yP1*d*log(yP1) + yP2*d*log(yP2);
endfunction

27
 Minimizing

• To solve the optimization we use the Octave function npsol [6] or MATLAB
function constr.

[x,obj,info]=npsol(x0,"gibbs",lb,ub,alb,a,aub)

• x0 is the initial guess as before and the new variables are used to define the
constraints (if any).

• The variable lb provides lower bounds for the extents and ub provides upper
bounds for the extents. These are required to prevent the optimizer from
“guessing” reaction extent values for which the Gibbs energy function is not
defined.

28
 Constraints

• For example, consider what the gibbs function returns given negative reac-
tion extents.

• Finally, we specify the sum of extents is less than 0.5 with the last three argu-
ments, which define linear inequality constraints other than simple bounds,

alb ≤ Aε0 ≤ aub

ε10 + ε20 ≤ 1/2

h i
A= 1 1 aub = 1/2

• We can set alb = 0 because both extents are already specified to be positive.

29
 Solving the Problem

• A typical session to solve for equilibrium composition by minimizing Gibbs

energy is
octave:1> a=[1 1]; alb=0.001; aub=0.499;
octave:2> lb=[0.001;0.001]; ub=[0.499;0.499];
octave:3> x0=[0.2;0.2];
octave:4> [x,obj,info]=npsol(x0,"gibbs",lb,ub,alb,a,aub)
x =
0.13317
0.35087
obj = -2.5593
info = 0

• The value of info = 0 indicates that the solution has converged

30
 Comparing Solutions

0.5
-2.557
0.45 -2.55
-2.53
0.4 -2.5
-2
0.35
0.3
-1
0
• The optimization results are in good agree-
0 0.25
ε2 ment with those computed using the alge-
0.2
0.15 braic approach, and the Gibbs energy con-
0.1
tours
0.05
0
0 0.1 0.2 0.3 0.4 0.5
0
ε1

• Optimization is a powerful tool for solving many types of engineering mod-

eling and design problems.

• We also rely heavily on optimization tools in Chapter 9 on parameter estima-

tion.

31
 Differential Equations (lsode)

The real workhorse for simulating chemical reactors is the ordinary differen-
tial equation (ODE) solver. A high-quality ODE solver is an indispensable tool of
the reactor designer because chemical reactions can take place at widely differ-
ent rates and time scales, presenting difficult numerical challenges. Type

help -i differential equations

for an overview of Octave’s ODE solver, lsode [7]. Type help ode23s for an
overview of MATLAB’s ODE solvers. We recommend always using the so-called
“stiff” solvers for chemical reactor modeling. The time scales for different reac-
tions may be widely separated, leading to stiff ODEs.

The ODEs required to describe transient-batch, CSTR, and steady-state PFR

32
reactors are of the form

dx
= f (x)
dt
x(0) = x0

in which x(t) is a vector of species concentrations to be found, t is time in the

batch reactor and CSTR, but volume or length in the steady-state PFR, and f is
the function defining the material balances for the species. If the differential
equations are nonlinear (f is not a linear function) and one must keep track of
several species simultaneously, then analytical solution is not possible and a
numerical procedure is required.

The first task is to formulate a function which defines f . If we examine the

benzene pyrolysis example, Example 4.5, a suitable function is

function xdot = benzene_pyrol_rhs(ext,volume)

33
 global k1 K1 k2 K2 R T P NBf
NB = NBf - 2*ext(1) - ext(2);
ND = ext(1) - ext(2);
NH = ext(1) + ext(2);
NT = ext(2);
Q = NBf*R*T/P;
cB = NB/Q; cD = ND/Q; cT = NT/Q; cH = NH/Q;
xdot(1) = k1*(cB*cB - cD*cH/K1);
xdot(2) = k2*(cB*cD - cT*cH/K2);
endfunction

in which the meanings of the symbol names should be apparent from the exam-
ple.

We then need to specify at which t values, volumes in this case, the solution
should be reported, the initial condition (in this case, feed condition), x0, and
call the ODE solver. A typical Octave session would be

34
 NBf = 60e3; R = 0.08205; T = 1033; P = 1;
k1 = 7e5; k2 = 4e5; K1 = 0.31; K2 = 0.48;
x0=[0;0];
vout = linspace(0,1600,50)’;
x = lsode("benzene_pyrol_rhs", x0, vout);
conv = ( 2*x(:,1) + x(:,2) ) /NBf;
yB = ( NBf - 2*x(:,1) - x(:,2) ) /NBf;
yD = ( x(:,1) - x(:,2) ) /NBf;
yT = ( x(:,2) ) /NBf;
yH = ( x(:,1) + x(:,2) ) /NBf;

in which the function linspace computes 50 linearly spaced points between

0 and 1600 L. Finally, any quantities such as conversions or yields, which are
computable from x, are calculated for plotting. The results of this calculation
are shown in Figures 4.20 and 4.21.

35
 Differential-Algebraic Equations (dassl)

Some reactor models require a more general structure than the ODE,
dx/dt = f (x, t). The nonconstant density, nonconstant volume semi-batch
and CSTR reactors in Chapter 4 are more conveniently expressed as differential-
algebraic equations (DAEs). To address these models, consider the more general
form of implicit ODEs
0 = f (dx/dt, x, t)
Both DAEs and ODEs can be considered special cases of this structure. Brenan et
al. [1] provide further reading on existence and uniqueness of solutions to these
models, which are considerably more complex issues than in the case of simple
ODEs. Initial conditions are required for dx/dt as well as x in this model,

dx
(t) = ẋ 0 x(t) = x 0, at t = 0
dt

36
Petzold has provided a numerical package, dassl, to compute solutions to im-
plicit differential, and differential-algebraic equations. The main difference be-
tween using dassl and lsode is the form of the user-supplied function defining
the model. A second difference is that the user must supply ẋ 0 as well as x 0.
The reader can consult the software for Example 8.1 at [Link]/
~jbraw/chemreacfun to see the details of how to use dassl.

37
 Automatic Stopping Times (dasrt)

We often need to stop ODE solvers when certain conditions are met. Two
examples are when we have reached a certain conversion, and when we have
created a new phase and need to change the ODEs governing the system. The
program dasrt enables us to stop when specified conditions are met exactly.
The following code illustrates how to find the PFR volume at which 50% conver-
sion of benzene is achieved. The user provides one new function, stop, which
the ODE solver checks while marching forward until the function reaches the
value zero. The PFR volume is found to be VR = 403.25 L, in good agreement
with Figure 4.20.

function val = stop(ext, volume)

global NBf xBstop
NB = NBf - 2*ext(1) - ext(2);

38
 xB = 1 - NB/NBf;
val = xBstop - xB;
endfunction

octave:1> x0 = [0;0]; xdot0 = [0;0];

octave:2> vout = linspace(0,1600,50)’; xBstop = 0.5;
octave:3> [x, xdot, vstop] = dasrt("benzene_pyrol_rhs", \
"stop", x0, xdot0, vout);
octave:4> vstop(length(vstop))
ans = 403.25

39
 Parametric Sensitivities of Differential Equations (dasac)

As discussed in Chapter 9, it is often valuable to know not just the solution

to a set of ODEs, but how the solution changes with the model parameters.
Consider the following vector of ODEs with model parameters θ.

dx
= f (x; θ) (2)
dt
x(0) = g(x 0; θ) (3)

The sensitivities, S, are defined as the change in the solution with respect to the

40
model parameters
!
∂xi
Sij = (4)
∂θj θk≠j

∂x
S= T (5)
∂θ

We can differentiate the Sij to obtain

" #
dSij d ∂xi ∂ dxi ∂
= = = fi(x; θ) (6)
dt dt ∂θj ∂θj dt ∂θj

Using the chain rule to perform the final differentiation gives

! ! !
∂fi X  ∂fi  ∂xk ∂fi
= + (7)
∂θj θk≠j k
∂xk xl≠k ,θj ∂θj θk≠j
∂θj xl ,θk≠j

41
Substituting this back into Equation A.4 and writing the sum as a matrix multiply
yields
dS
= fx S + fθ (8)
dt
in which the subscript denotes differentiation with respect to that variable. No-
tice the sensitivity differential equation is linear in the unknown S. The initial
conditions for this matrix differential equation is easily derived from the defini-
tion of S
∂xi(0) ∂gi
Sij (0) = =
∂θj ∂θj

dS
= fx S + fθ
dt
S(0) = gθ

Notice that if none of the initial conditions are parameters for which one is
computing sensitivities, then S(0) = 0.

42
Example 0.2: Two simple sensitivity calculations

Consider the differential equation describing a first-order, isothermal, irre-

versible reaction in a constant-volume batch reactor

dcA
= −kcA
dt
cA(0) = cA0 (9)

1. First write out the solution to the differential equation.

2. Next consider the rate constant to be the parameter in the problem, θ = k.

Take the partial derivative of the solution directly to produce the sensitivity
to this parameter, S = ∂cA/∂k.

3. Next write the differential equation describing dS/dt. What is S(0)? Solve
this differential equation and compare to your previous result. Plot cA(t), S(t)

43
 versus t.

4. Repeat steps 2 and 3 using the initial concentration as the parameter in the
problem, θ = cA0.

Solution

1. Equation A.5 is our standard, first-order linear equation, whose solution is

cA = cA0e−kt

2. Taking the partial derivative of cA(t) with respect to θ = k produces directly

S1 = −tcA0e−kt (10)

44
Evaluating the various partial derivatives gives

fx = −k fθ = −cA gθ = 0

so the sensitivity differential equation is

dS1
= −kS1 − cA
dt
S1(0) = 0

This equation is linear and can be solved with the method used in the series
reactions in Chapter 4, or with Laplace transforms as in Exercise 4.6. The
result is
S1(t) = −tkcA0e−kt (11)
which agrees with Equation A.6. The functions cA(t), S1(t) are shown in Fig-
ure A.3.

45
 2

1.5

1
cA

0.5
S2
0

-0.5 S1

-1
0 1 2 3 4 5
t

Figure 1: Solution to first-order differential equation dcA/dt = −kcA, and sen-

sitivities S1 = ∂cA/∂k and S2 = ∂cA/∂cA0.

46
3. Taking the partial derivative of cA(t) with respect to θ = cA0 produces directly

S2(t) = e−kt (12)

Evaluating the various partial derivatives for this parameter choice gives

fx = −k fθ = 0 gθ = 1

so the sensitivity differential equation is

dS2
= −kS2
dt
S2(0) = 1

This equation is our standard, first-order linear equation whose solution is

S2(t) = e−kt

47
 which agrees with Equation A.7. This sensitivity also is plotted in Figure A.3.
Notice that only when the initial condition is the parameter is the sensitivity
initially nonzero.

Stewart and Caracotsios produced the code dasac to perform sensitivity cal-
culations conveniently [3, 2]. The following Octave commands show how to use
dasac to solve Example A.2.

function xdot = firstorder(t, x, theta)

global ca0 k
ca = x(1);
k = theta(1);
xdot = -k*ca;
endfunction

x0 = ca0; theta = [k, ca0]; sx0 = [0, 1];

48
[x, xdot, sx, sxdot] = \
dasac("firstorder", x0, theta, sx0, tout);
ca = x;
dcadk = nth(sx, 1);
dcadca0 = nth(sx, 2);

49
 Boundary-Value Problems and Collocation (colloc)

In the collocation method, we approximate a function by passing a polyno-

mial through values of the function at selected points. The selected points are
known as collocation points. The locations of the collocation points have a large
impact on how well the method works. Evenly spaced points, which seems a nat-
ural first choice, turns out to have mediocre properties. Choosing the points as
zeros of a member of a family of orthogonal polynomials turns out to have much
better properties. This choice is referred to as orthogonal collocation.

50
 c
dc X
= Aij cj
dr ri
j

r1 r2 r3 r4 r5
r

Figure 2: Function c(r ) and its values at five collocation points. Derivatives and
integrals of the polynomial interpolant are linear combinations of the function
values at the points.

Figure A.4 shows values of a function c(r ) at five collocation points. The

51
function is approximated by passing a polynomial through these points. To
solve differential equations and boundary-value problems (BVP), we first com-
pute the required derivatives of the polynomial approximation. Then we then
find the values of the function such that the differential equation is satisfied at
the collocation points. If we increase the number of collocation points, nc , we
require the differential equation to be satisfied at more locations, and we obtain
a more accurate solution.

The derivatives and integrals of the polynomial interpolant can be computed

52
as linear combinations of the values at the collocation points

nc
dc X
= Aij cj
dr ri
j=1
nc
d2c X
= Bij cj
dr 2 ri j=1
Z1 nc
X
f (c)dc = Qj f (cj )
0
j=1

To obtain the locations of the collocation points and derivatives and integral
weighting matrices and vectors, we use the function colloc, based on the meth-
ods described by Villadsen and Michelsen [8].

[R A B Q] = colloc(npts-2, "left", "right");

53
The strings "left" and "right" specify that we would like to have collocation
points at the endpoints of the interval in addition to the zeros of the orthogonal
polynomial. We solve for the concentration profile for the reaction-diffusion
problem in a catalyst pellet to illustrate the collocation method.

Example 0.3: Single-pellet profile

Consider the isothermal, first-order reaction-diffusion problem in a spherical

pellet
1 d dc
 
2 2
r − Φ c=0 (13)
r 2 dr dr

c = 1, r =3
dc
= 0, r =0
dr

54
The effectiveness factor is given by

1 dc
η=
Φ2 dr r =3

1. Compute the concentration profile for a first-order reaction in a spherical

pellet. Solve the problem for φ = 10.

2. Plot the concentration profile for nc = 5, 10, 30 and 50. How many collocation
points are required to reach accuracy in the concentration profile for this value
of Φ.

3. How many collocation points are required to achieve a relative error in the
effectiveness factor of less than 10−5?

55
 1
nc = 50
nc = 30
0.8 nc = 10
nc =5
0.6

0.4
c
0.2

-0.2

-0.4
0 0.5 1 1.5 2 2.5 3
r

Figure 3: Dimensionless concentration versus dimensionless radial position for

different numbers of collocation points.

56
Solution

First we perform the differentiation in Equation A.8 to obtain

d2c 2 dc
2
+ − Φ2 c = 0
dr r dr

We define the following Octave function to evaluate this equation at the interior
collocation points. At the two collocation endpoints, we satisfy the boundary
conditions dc/dr = 0 at r = 0 and c = 1 at r = 3. The function is therefore

function retval = pellet(c)

global Phi A Aint Bint Rint n
npts = length(c);
cint = c(2:npts-1);
retval(1) = A(1,:)*c;
retval(2:npts-1) = Bint*c .+ 2*Aint*c./Rint .- \
Phi^2*cint;

57
 retval(npts) = 1 - c(npts);
endfunction

Figure A.5 shows the concentration profiles for different numbers of collocation
points. We require about nc = 30 to obtain a converged concentration profile.
Figure A.6 shows the relative error in the effectiveness factor versus number of
collocation points. If one is only interested in the pellet reaction rate, about 16
collocation points are required to achieve a relative error of less than 10−6.

58
 1

10−2

10−4

10−6
η error 10−8

10−10

10−12

10−14

10−16
5 10 15 20 25 30
nc

Figure 4: Relative error in the effectiveness factor versus number of collocation

points.

59
Example 0.4: A simple fixed-bed reactor problem

Here we set up the collocation approach of this section for the single, second-
order reaction

2
A -→ B r = kcA (14)

Solve the model exactly and compare the numerical solution to the solution ob-
tained with the approximate Thiele modulus and effectiveness factor approaches
shown in Figure 7.26.

60
Solution

The model for the fluid and particle phases can be written from Equations 7.67–
7.77
dNA Sp dceA
= RA = −(1 − B ) DA
dV Vp dr r =R
d2ceA 2 dceA ReA
+ + =0
dr 2 r dr DA

ceA = cA r =R
dceA
=0 r =0
dr

The collocation method produces algebraic equations for the pellet profile as for
the previous example. We use a DAE solver to integrate the differential equation
for NA and these coupled algebraic equations.

61
 10

NA (mol/s)
7

6
1 1 1
 
5 η= −
Φ tanh 3Φ 3Φ
1
4 η=
Φ
3

2
0 50 100 150 200 250 300 350 400
VR (L)

Figure 5: Molar flow of A versus reactor volume for second-order, isothermal

reaction in a fixed-bed reactor; two approximations and exact solution.
3

2.8
NA (mol/s)

2.6 62

2.4
 The results for NA versus reactor length using 25 collocation points for the
pellet are shown in Figure A.7. Also shown are the simplified effectiveness factor
calculations for this problem from Example 7.5. A magnified view is shown in
Figure A.8. Notice the effectiveness factor approach gives a good approximation
for the bed performance. It is not exact because the reaction is second order.


63
 References

[1] K. E. Brenan, S. L. Campbell, and L. R. Petzold. Numerical Solution of Initial-

Value Problems in Differential-Algebraic Equations. Elsevier Science Publish-
ers, New York, 1989.

[2] M. Caracotsios and W. E. Stewart. Sensitivity analysis of initial value problems

with mixed ODEs and algebraic equations. Comput. Chem. Eng., 9(4):359–
365, 1985.

[3] M. Caracotsios and W. E. Stewart. Sensitivity analysis of initial-boundary-

value problems with mixed PDEs and algebraic equations. Comput. Chem.
Eng., 19(9):1019–1030, 1995.

[4] J. W. Eaton. Octave: Past, present and future. In K. Hornik and

F. Leisch, editors, Proceedings of the 2nd International Workshop on Dis-
tributed Statistical Computing, March 15-17, 2001, Technische Universität

64
 Wien, Vienna, Austria, 2001. [Link]
2001/Proceedings/, ISSN 1609-395X.

[5] J. W. Eaton and J. B. Rawlings. Octave—a high level interactive language for
numerical computations. CACHE News, (40):11–18, Spring 1995.

[6] P. E. Gill, W. Murray, M. A. Saunders, and M. H. Wright. User’s guide for

SOL/NPSOL (Version 4.0): A Fortran package for nonlinear programming,
technical report SOL 86-2. Technical report, Systems Optimization Labora-
tory, Department of Operations Research, Stanford University, 1986.

[7] A. C. Hindmarsh. ODEPACK, a systematized collection of ODE solvers. In

R. S. Stepleman, editor, Scientific Computing, pages 55–64, Amsterdam,
1983. North-Holland.

[8] J. Villadsen and M. L. Michelsen. Solution of Differential Equation Models

65
by Polynomial Approximation. Prentice-Hall, Englewood Cliffs, New Jersey,
1978.

66