Suboptimal Design of a n SO?-Oxidation Catalytic Reactor J. David Paynter', Joshua S. Dranoff, and S. George Bankoff2 Chemical Engineering Department, Northwestern University, Euanston, Ill. 60201
A fixed-volume tubular reactor for oxidation of SO2 over commercial VZOScatalyst i s optimized with respect to reactor radius and surface heat flux as a function of axial distance. The suboptimal method employs ordinary functional-differential adjoint equations, rather than partial differential equations. It differs from weighted residual approximation methods in that the dimension of the state and adjoint equations are not the same. The suboptimal design shows a nearly constant radius in the absence of a temperature constraint.
D e s p i t e a fairly extensive literature, dating back more than 15 years, on the optimal design of tubular reactors by variational methods, a recent informal survey showed that there has been almost negligible impact to date upon the actual design procedures used in specifying commercial reactor equipment. I n part, this may be attributed to the usual information gap between researchers and practitioners; but, in large measure, it is also clear that the problems which have been solved to date are vastly oversimplified. Attempts t o solve more realistic, but correspondingly more complex, models have frequently foundered, owing to inefficient computational methods, or have been extremely costly in terms of computer time. Furthermore, the solution surface in the neighborhood ' Present address, Esso Research Laboratories, Baton Rouge, La. 70821 ' To whom correspondence should be addressed.
244
Ind. Eng. Chem. Process Des. Develop., Vol. 10, No. 2, 1971
of the optimum appears often to be fairly flat, so that one can obtain close approximations to the optimal objective functions (usually taken as the yield of a desired product) with quite different control policies. Thus, it is important to develop a range of suboptimal computational techniques which allows the user to choose between economy of effort and proximity to the true optimum, depending upon the nature of the problem which confronts him. A start in this direction has been made in a previous paper, in which a partial averaging technique was used to compute suboptimal heat removal rate policies as a function of axial distance for a two-dimensional catalytic reactor for oxidation of SOZ. This procedure, while specialized, reduces the computational effort by a factor of approximately two compared to the full two-dimensional calculation, and was shown, for the fictitious problem posed by Denn et al. (1966), to give equivalent estimates of the yield.
I n fact, it is rather difficult to calculate the “true” optimal yield. Denn, in a later publication (Denn, 1969), pointed out that slightly higher yields can be obtained by a simple search on the best isothermal jacket temperature for a reactor configuration consisting of an adiabatic section, followed by a jacketed section with constant-temperature cooling water! I t is not uncommon, as in this case, for the two-dimensional gradient calculation to fail to make progress after only a few iterations. This behavior is also shown in the suboptimal calculations (Paynter et al., 1970). T o make further progress toward the optimum, it is necessary to employ ad hoc acceleration methods or restarting methods (Bandy and Bankoff, 1969a; Paynter, 1968), or else one must have recourse to second-order methods (not applicable in this particular case, because of the linearity of the Hamiltonian function with respect t o control) (Bryson and Ho, 1969; Padmanabhan and Bankoff, 1969). I n any case, the uncertainty of the two-dimensional gradient approach t o optimality makes suboptimal methods, such as the partial averaging method employed herein, considerably more attractive. On the other hand, one can oversimplify the problemi.e., a one-dimensional calculation (which neglects all radial gradients) gives a heat flux removal policy which results in substantially lower yields, although predicting higher yields than the two-dimensional result. The partial averaging technique is a hybrid approach, in that the full-state equations are integrated forward, but the adjoint equations are ordinary functional-differential equations, rather than partial differential equations. The computational time is thereby approximately halved, while retaining information concerning the radial gradients throughout the calculation. I t differs thereby from weighted residual approximation methods, such as Galerkin or collocation techniques, in that the dimensions of the state and adjoint equations are not the same. One of the advantages, of course, of the variational method is that one can explore nonconventional designs, such as reactor tubes whose shape is to be chosen optimally, to avoid “hot-spot” problems. With natural-circulation boiling loops, the optimal heater shape gives u p to 5 0 5 greater power as compared to the conventional constantradius design (Bandy and Bankoff, 1969a,b). In the present instance, however, the optimal shape is fairly close to a constant radius, so that the constant-radius design appears to be preferable, from a practical point of view. This information, however, cannot be obtained by recourse to engineering “common sense.” I n this work we examine the variable-radius problem in some detail. Problem Formulation
For a tubular reactor for the oxidation of SO? over a vanadium pentoxide catalyst, the governing dimensionless equation for the extent of reaction is: doc 1 d ax=@’-R dR
doc
RdR)
Here
(3) and it is assumed t h a t radial variations in gas specific volume and mass velocity can be neglected, to simplify the calculations. Similarly, a heat balance gives
(4) The boundary conditions for Equations 1 and 4 are:
oc(R,O)= 0 T(R,O) = 1 doc
d7
doc
dR
dR
dR
- (0,X) = - ( 0 , X ) = - ( 1 , X ) =
0
where q (X), the dimensionless wall heat flux, satisfies upper and lower inequality constraints, corresponding to maximum and zero cooling: Q, 5
q(X) 5 qx
The rate expression, due to Calderbank (1952, 1953), is given in the Appendix. We wish to maximize the radial-average extent of reaction, ;(1), where overbar denotes the averaging operation:
Similarly, one can average Equations 1 and 4 to obtain
dZjdX = 2 d?/dX
=2
p q ( X )+ 2
l1RRdR;
-
a(0) = 0
RRTdR; ?(O) = 1 i1
(7)
(8)
The ordinary functional-differential equations for the adjoint variables are:
+a
where with boundary conditions:
h(1) = 1; Xz(1) = 0 Ind. Eng. Chem. Process Des. Develop., Vol. 10, No. 2, 1971
(11) 245
I.o
I
I
I
I
I
I
I
The following set of equations represents the model of this variable-radius system:
I
dff _ -
0.9
z 0 I-
ax
p"
o
z
r
a7
_ -
ax -
0.8
\
( i - 1) Ax 5 x,< i A X 1 = 1, . . . , N'
w
E ;
where
0.7
0.6 0
2
4
10 12 ITERATION NUMBER 6
8
14
16
18
Figure 1. Profit function vs. iteration number for variableradius reactor
The suboptimal necessary conditions (Paynter, 1968) are
(RT)
=
(-aHd
p C (1 -
4 Lr,
If the number of axial increments is increased, the model can be represented by the following pair of parabolic equations:
except if A2 = 0 over some finite interval. From Equation 10 it follows that a necessary condition for a suboptimal singular arc is that
1'R $ dR
=0
on the arc. Note that a similarly simple switching criterion is not available when the full set of partial differential equations is employed. Variable-Radius Reactor
T o enlarge the space of control functions, and hopefully thereby to improve the average exit yield, we pose the problem of selecting both the optimal shape (reactor radius as a function of axial distance) and the optimal heat removal function. One can approximate the variable-radius reactor as a series of constant-radius reactors, each of different radius. The dimensionless radius is now defined with respect to the outside radius of the reactor which is a function of X. This introduces a theoretical difficulty because N and T are assumed continuous with respect to the dimensionless radius and length but are therefore discontinuous with respect to the real radius. The dimensionless radius could be defined with respect to the radius at X = 0 but this would transfer the discontinuities to the dimensionless variables and would be difficult to handle numerically, owing t o the irregular boundary of the region of interest. 246
Ind. Eng. Chem. Process Des. Develop., Vol. 10, No. 2, 1971
In this case the maximum exit extent of reaction is sought by properly choosing, not only the heat flux profile, but also the reactor radius profile. The heat flux is constrained to lie between 0 and 3, and the radius is constrained to be positive. No constraint is placed on the temperature. T o make valid comparisons, the volume of catalyst must be equal to that of the constant-radius case. This imposes a volume constraint on the system, which is accounted for by introducing a penalty function in the objective function, as follows:
-?[I1{
P=ff(1)
RL,(X)- R:\dX]?
(18)
where R1 is the equivalent constant radius. The state of the system is given by Equations 16 and 17, and the boundary conditions (Equation 5 ) . From these, the following averaged state equations are obtained:
For convenience, a n additional state variable is introduced, satisfying the following equation:
d ~ , / d X= R6 ( X ) - R: ; ~ ~ (=00 )
(21)
This enables the objective function t o be written as a function of the final state only: -
P = a ( 1 ) - ?lx: (1)
(22)
To find the necessary conditions for optimality, the maximum principle is again used. The Hamiltonian function in this case is
H = -2 qxO(l) R i ( X ) - R i ] +
i-,
I
5.01
Y)
0
i-
1
4.0
2.0
3z 0.0 kz -1.0 8
[
w -2.0 J
-3.0
1
1
-4.0
8
4
0
where
12"
I8
ITERATION NUMBER Figure 2. Volume constraint violation vs. iteration number (note that penalty i s negligible by 12th iteration)
The algorithms used for iterating on the radius profile and heat flux were:
(X)= q 1 ( X ) +
q'-l
Ri,'
I
( X ) = R ( X ) + € 2 [-4
QX,(~)
R,(X) +
violation of the volume constraint, as a function of the number of iterations. The coefficient, q , was adjusted from a value of lo6 t o lo8 at the fourth iteration, giving rise to the large decrease in the profit function a t that point. This adjustment was made to reduce the volume constraint violation, with the effect shown in Figure 2. About 39.43 min of CDC-3400 computer time were required to perform the 18 iterations, resulting in a profit function of 0.9410. The evolution of the decision variables is shown in Figures 3 and 4. The suboptimal reactor radius is, in fact, very nearly a constant, the maximum variation being about 2%. Similarly, the profit function is approximately the same as that found earlier for the constant-radius case
3
0.7
k
0.6
W
I
0.4 The initial approach to this problem was to vary both heat flux and reactor radius profiles a t once. Figure 1 summarizes this gradient climb. By the sixth iteration it was apparent that changing both variables a t once is very costly in terms of computer time. I n the next two iterations the heat flux profile was fixed, while the radius profile was allowed to vary. Then the radius function was fixed, and an attempt was made to vary the heat flux profile, but no successful move could be made. From the ninth through the eighteenth iteration only the radius profile was varied in searching for the optimal design. Figure 2 shows ~ ( l which ), is a measure of the
W
P E 0.2 O.I
t
0.0I 0.0
I
II
I
I
0.2
0.4
0.6
0.8
I
I.o
DIMENSIONLESS REACTOR LENGTH Figure 3. Heat flux profile evolution Ind. Eng. Chem. Process Des. Develop., Vol. 10, No. 2, 1971
247
N
P x
3'0F--l
2 2.4
4.30 4'35
(0
4.25 4.20
i
+
$4.I5 P
'4.10 a 0 4.05 d
w a4.00
t
1
1
i
1
3.90 3.95 0.0
0.2 0.4 0.6 0.8 1.0 DIMENSIONLESS REACTOR LENGTH
Figure 4. Reactor shape by gradient method in hybrid space
4 8 12 16 20 ITERATION NUMBER
0
Figure 6. Volume constraint violation vs. iteration number
1 " "
4.20 N
t
INITIAL
0
4.15 c +
vi
1
n
d
18th
4.10
a
P0 9
!
a 4.05
0
I
I
I
I
I
0.0
0.2
0.4
0.6
0.6
4.00
0.90 12 16 20 ITERAT10N NUMB E R
4
8
24
I
DIMENSIONLESS REACTOR L E NGTH
Figure 5. Objective function vs. iteration number, a d hoc gradient procedure
Figure 7. Suboptimal shape; constant-radius initial profile
(Paynter et al., 1970). The reason is that the exit yield is very nearly the equilibrium yield a t the exit temperature, so that the optimum is fairly insensitive to details of the control functions prior to the exit (although quite sensitive to the integrals of these functions). T o examine further the neighborhood of these suboptimal policies, some additional computations were performed. I n Figure 5 the one-dimensional heat flux subspace was searched while keeping the reactor radius profile fixed. After the fifth iteration, no further improvements could be obtained by varying the heat flux profile. From this point to the eighteenth iteration only the reactor radius profile was changed. Thereafter, an attempt was made to vary both che heat flux and reactor radius profiles, but no improve-
ment in the profit function could be found. The difference between the profit function a t the fifth iteration, 0.9439, and the eighteenth iteration, 0.9473, is the improvement in the profit from allowing the reactor radius t o vary as a function of length. The corresponding history of the volume constraint violation is shown in Figure 6. Three additional searches were made, each starting from a different initially assumed profile of reactor radius. For all three searches the heat flux profile assumed was the one found after five iterations in Figure 5 . In each search the reactor radius was allowed to vary with the heat flux profile fixed. Periodically, attempts were made to vary the heat flux profile, but in no case could a successful move be accomplished. The three initially assumed radius
248
Ind. Eng. Chem. Process Des. Develop., Vol. 10, No. 2, 1971
profiles along with the profiles obtained at the end of each search are shown in Figures 7 to 9. The resulting profit functions are summarized in Table I. Because different reactor radius profiles were obtained from the gradient search, depending on the profile initially assumed, four admissible suboptimal designs were found. Because of the singular problem associated with the heat flux profile and the flatness of the optimum, the gradient method cannot discriminate between alternate suboptimal designs. Appendix 0
4a
4.00-
-
0.2 0.4 0.6 0.8 1.0 DIMENSIONLESS REACTOR LENGTH
0.0
Figure 8. Suboptimal shape; step-function initial profile
SOs Oxidation Kinetics. The Calderbank rate expression is used in this study for the kinetics describing the oxidation of SOe.This rate expression is used for the following three reasons: This kinetic model was obtained from data taken using commercial catalyst. Other models reported in the literature are based on data taken using specially prepared catalyst . There is evidence that this model best represents the behavior of many commercial operations (Calderbank, 1953). This model has been used in other optimization studies so that direct comparison of results can be made (King, 1965; Lee and Aris, 1963). The kinetic expression is as follows:
5.0
where
rso?= rate of oxidation of SO, p l = partial pressure of ith species k, = exp [-E,/1.98 T + F , ] i = 1,2
N
0
4.5
X
c +
vi
(1)
In this study i t is more convenient to work in terms of extent of reaction than partial pressure. The partial pressures can be rewritten in terms of extent of reaction as follows:
2U a 4.0 a
E 2a
3.5
where
3.0 0.4 0.6 0.8 1.0 DIMENSIONLESS REACTOR LENGTH
0.0
Figure
0.2
9. Suboptimal shape; tapered initial profile
~
~
~~
Computer
P
Constant radius ( I ) Step profile (2) Tapered profile (3)
0.9473 0.9494 0.9150
X"(1)
-4.13 x lo-* 8.14 x lo-' 2.64 x
= total pressure
Y ?= initial composition of ith species and
~
Table I. Profit Functions for Various Initial Profiles initial radius profile
n
time required
36.30 min 18.10 min 28.20 min
Nomenclature
C,
= specific heat, cal/gmol K d, = catalyst particle diameter, ft D , D T = radial diffusivity for mass and heat, ftL/hr E, = activation energy, cal/gmol K F , = preexponential constant G = mass flux, lb/ft' hr O
Ind. Eng. Chern. Process Des. Develop., Vol. 10, No. 2, 1971 249
Go = mass flux a t bed entrance -AHf$ = heat of reaction, Btu/lb mole H = Hamiltonian function k , = rate constant, lb mole/(lb catalyst) sec k = thermal conductivity, B t u / h r f t K L = reactor length, ft h? = number equal to molecular weight entering reactor M , = molecular weight of ith species N’ = number of radial increments = l / A R N = ipo/p)idDUr,/D) NT = (po/p)idpUo/D~) P, = partial pressure of ith species, atm P = profit or objective function q % ,q” = lower and upper heat flux constraint q = dimensionless heat flux, QR,/kT, Q = heat flux, Btu/hr ft2 r = radial coordinate, ft R , , ( X ) = reactor radius, ft r, = rate of reaction of ith species, lb mole/hr lb catalyst R_ = dimensionless radius, r / R , ( X ) R = M L (1 - e ) prrso.iYSO, G,, l? r = (- H R L) (1 - e ) pcrso,/iT,G,c,) T = temperature, K T , = entering temperature, U,, = average gas velocity a t bed entrance, ftisec X = dimensionless length, z / L Y , = concentration, lb mole i / M lb gas Y , = initial concentration z = axial coordinate, ft OK
= adjoint variable II = total pressure, atm p = density of gas, lb/ft3 p. = density of gas a t bed entrance p. = catalyst density, lb/ft3 T = dimensionless temperature, T / Tu Literature Cited
Bandy, D. B., Bankoff, S. G., “Progress in Heat and Mass Transfer,” Vol. I , U. Grigull and E. Hahne, Eds., Pergamon Press, Elmsford, N. Y., pp 425-71, 1969a. Bandy, D. B., Bankoff, S. G., Nucl. Eng. Des., 11, 77 (196913). Bryson, A. E., Jr., Ho, Y. C., “Applied Optimal Control,” Blaisdell Publishing Co., New York, N. Y., 1969. Calderbank, P. H., Chem. Eng. Progr., 49, 585 (1953). Calderbank, P. H., J . A p p l . Chem. (London), 2,482 (1952). Denn, M. M., “Optimization by Variational Methods,” McGraw-Hill, New York, N. Y., 1969. Denn, M. M., Gray, R. D., Ferron, J. R., Ind. Eng. Chem. Fundam., 5 , 59 (1966). King, R. P., Chem. Eng. Sei., 20, 537 (1965). Lee, K., Aris, R., Ind. Eng. Chem. Process Des. Develop., 2, 300 (1963). Padmanabhan, L., Bankoff, S. G., Proceedings of the Joint Automatic Control Conference, Boulder, Colorado, 1969, pp 34-43. Paynter, J. D., Ph.D. Thesis, Northwestern University, Evanston, Ill., 1968. Paynter, J. D., Dranoff, J. S., Bankoff, S. G., Ind. Eng. Chem. Process Des. Develop., 9, 303 (1970).
GREEKLETTERS N
= extent of reaction, = d,LiNRd ( 0 )
@, 8T
RECEIVED for review March 11, 1970
(Yso - Y s o ) /Ys0
= d,L/NTR;‘(O)
ACCEPTED November 25, 1970 Appreciation is extended to the National Science Foundation (Grant GP-2175) for financial support, and to the Vogelback Computing Center at Northwestern University for a generous donation of computer time.
= step size in the gradient method t H = void fraction of bed p = coefficient of penalty function t
Development of a Model for Kinetics of Olefin Codimeriza tion J. David Paynter’ and William 1. Schwette Esso Research Laboratories, Humble Oil and Refining Co., Baton Rouge, La. 70821 Frequently the process development engineer must analyze and extrapolate data on complex kinetic systems. Often three or more highly coupled reactions are involved. Uncoupling these reactions can be a formidable task. In addition. the data available are generally less than ideal.
’ To whom correspondence 250
should be addressed.
Ind. Eng. Chern. Process Des. Develop., Vol. 10, No. 2, 1971
Often experimentation is designed to develop and demonstrate a process rather than to develop a kinetic model. The press of time and money frequently does not permit model discrimination experiments as suggested by Box and Hill (1967). This paper illustrates techniques for kinetic model building which are useful in the situations mentioned
