Ind. Eng. Chem. Fundam. 1986, 25, 313-321
f = dimensionless rate function, u(l
+ u)2/(1 + U U ) ~
F = function defined by eq 25 F = ignition value of Da, F(zi*;zio) G = function defined by eq 16 hi(y) = static bifurcation values, defined by eq 23 H = function defined by eq 18 Hi-=hi(1) k, k = local and volume average reaction rate constant K = adsorption equilibrium constant L = reactor length n = integer characteristic of pellet shape; = 0 for infinite slab, = 1 for infinite cylinder, = 2 for sphere R = specific exchange rate between the external fluid and solid phases, defined by eq 4 R, = characteristic pellet dimension; half-thickness (n = 0), radius (n = 1, 2) s = x/R, = dimensionless location of the active catalyst within the pellet t = time
u = c/cp a = value of u at s = 8 a,, a* = roots of G'(zi) = 0 ii+, zi- = roots of H'(a) = 0, defined by eq 21 d* = parameter defined by eq 28 u = velocity of external fluid phase V = Liapunov's functional V , = pellet volume x = distance from pellet center X = conversion at reactor outlet, 1 - y ( 1 ) y = z/L z = axial coordinate along the reactor
Greek Letters
B = parameter defined by eq 14 B+, p- = parameters defined by eq 31a
313
Y = Cf/CfO
6 = Dirac delta function 6 = bed void fraction q = effectiveness factor I(. = perturbation variable with respect to steady-state dimensionless concentration, U ( S , T ) - u&) u = dimensionless adsorption constant, KCfo T = dimensionless time, tR 2 / D , 42 = Thiele modulus, hR,1/[(1 + U ) ~ D , ] #,, = function defined by eq 11 Subscripts 0 = pellets at the inlet section 1 = pellets at the outlet section s = steady state Superscript O = feed value
Literature Cited Ark, R. The Methemtkal Theory of Diffusion and Reaction In Permeable Catalysfs;Clarendon: Oxford, 1975: Vol. 2. Berger. A. J.; LapMus, L. A I C M J . 1968, 14, 558. Boyce, W. E.; DiPrlma, R. C. Elementary DifferenHel Equatbns and Bwndaty Value Problems, 3rd ed.; Wiley: New York, 1977. Hegedus, L. L.; Oh, S. H.; Baron, K. AIChE J . 1977, 2 3 , 632. Iooss. 0.;Joseph, D. D. €lemntary Stability and Bifurcation Theory; Springer-Verlag: New York, 1980. Luss, D. Chem. Eng. Scl. 1971. 26. 1713. MorbMelli, M.; Servide, A.; Carra, S.; Varma, A. Ind. Eng. Chem. Fundam. 1985, 2 4 , 116. MorbMeili. M.; Servida, A.; Varma, A. Ind. Eng. Chem. Fundam. 1982. 21, 278. MorbMelli. M.;Varma. A. Ind. Eng. Chem. Fundam. 1982, 21, 284. Varma, A. I d . Eng. Chem. Fundam. 1980, 19. 316.
Received for review October 31, 1983 Revised manuscript received September 19, 1984 Accepted August 23,1985
Optimal Catalyst Activity Profiles in Pellets. 5. Optimization of the Isothermal Fixed-Bed Reactor Masslmo Morbidelli, Albert0 SewIda,+ Serglo
and Arvlnd Varma
Department of Chemical Engineering, Universi3/ of Notre Dame, Notre Dame, Indiana 46556
Nonuniform activity distribution within catalyst pellets can significantly improve performance of an isothermal fixed-bed reactor. For a bimolecular Langmuir-Hinshelwood reaction, the rigorous optimal distribution is found to be a Dirac delta function centered at a location insMe the pellet which depends on position along the reactor axis. The more practical situation where all particles have the 581118 Dirac delta distribution, the location of which is optimized with respect to the outlet conversion, is investigated. This is found to lead to better performance than the usual surface step or uniform activity distributions. However, particularly for high conversions, the rigorous optimal distribution still exhibits better performance, thus indicating the possible existence of other more convenient suboptimal activity distributions.
1. Introduction This paper deals with optimization of the outlet conversion of a fixed-bed catalytic reactor through a proper distribution of an active catalyst within the pellets. This problem is addressed to the case of an isothermal plug-flow reactor where a bimolecular Langmuir-Hinshelwood reaction occurs and interphase transport resistance is neg-
Present address: Dipartimento di Chimica Fisica Applicata, Politecnico di Milano, 20133 Milano, Italy.
ligible. Similar to the case of a single catalyst pellet examined recently (Morbidelli et al., 1982; Morbidelli and Varma, 1982), we shall see that activity distributions of the Dirac delta type actually constitute the rigorous solution of this optimization problem. Once this is shown, the full characterization of reactor behavior reported in part 4 (Morbidelli et al., preceding paper) is applied extensively, following the same nomenclature. It is worth pointing out that most of the published papers in this area have dealt with effects of specific nonuniform activity distributions on performance of single
0196-4313/86/1025-0313~01.50/0 0 1986 American Chemical
Society
314
Ind. Eng. Chem. Fundam., Vol. 25, No. 3, 1986
catalyst pellets (to the list reported by Morbidelli et al. (1982), the following contributions may be added: Oh et ai., 1980; Komiyama et al., 1980; Shyr and Ernst, 1980; Melo et al., 1980; Juang et al., 1981; Johnson and Verykios, 1983,1984). Very recently, Cukierman et al. (1983) have considered the catalyst design problem for an isothermal plug-flow reactor with a van de Vusse reaction network, again limiting the analysis to some specific activity distributions. The rigorous solution to the reactor optimization problem, with internal diffusion-reaction within the pellets, even in the case of a single isothermal reaction, has never been reported before. 2. Rigorous Optimization of Activity Distribution The reactor model considered has been described in detail in part 4, as constituted by the external fluid-phase mass balance and the pellet mass balance (given by eq 5-7 of part 4). The goal is to maximize the reactor outlet conversion, defined as X = 1- y ( l ) , which can be evaluated by integrating once the external fluid-phase mass balance 1
X = D a L Ny) dy
might be adopted in practice can be compared. Solving the pellet mass balance with the activity distribution (2) and condition (3), we can derive the following explicit expressions for the optimal location S = l - A , forn=O = exp(-A), = 1/(1
for n = 1
+ A),
for n = 2
where A =
4(n + l ) ( u y - 1) (1
+ u)2$2
M ( y ) = f ( U M ) = &f
represents the contribution to the reactant conversion at location y along the reactor axis. The optimization variable here is the pellet activity distribution function, a(s,y) (which, in general, depends also on the position along the reactor axis), which must satisfy the usual integral condition (see eq 8 of part 4) deriving from its definition, as well as the constraint imposed by the pellet mass balance (eq 7 of part 4) for the intrapellet concentration profile, 4 s ,Y) The exact analytical expression of the optimal activity distribution, a(s,y), can be readily derived noting that, since My) > 0, the maximum of X is obtained by locally maximizing My) at each y E [0,1]. Thus, the problem is reduced to optimizing, at each y, the single-particle performance, whose solution has been reported previously (Morbidelli et al., 1982). The difference between the two cases is in the boundary condition for u, which now varies with position along the reactor, i.e., u(1,y) = 701)5 1. By use of the same arguments, it can be shown that My) attains its absolute maximum value, M(y) = f(u,) (where u, is the concentration value which maximizes the reaction rate f ( u )in the interval u E [O,y]), if the following activity distribution is adopted
.
(5)
These relationships do not hold in the cases u < 1 or u > 1 and y C uM l/u; in these cases the reaction rate f ( u ) attains its maximum value in the interval u E [O,y] at the point u, = y, so that from eq 3 it follows that S = 1 (i.e., the active catalyst is concentrated at the external pellet surface), and the obtained value of My) is given by
MY) = f ( r ) (la)
(44
so that, since y = yo),s is also a function of position y. The resulting value for My) is given by
where My) = ( n + l)xlsna(s,y)f(u) ds
(4)
(6)
Note also that, even for u > 1and y > uM,eq 4 for the slab geometry ( n = 0) can lead to physically meaningless values of s; i.e., S 0 for $2
ni,*$
= 4(uy - 1)/(1+ u)2
(7)
As discussed in detail for the case of a single pellet before (Morbidelli et al., 1982),the optimal activity distribution is then given by eq 2 with s = 0; in this case My) does not attain its optimum value, i.e., x(y) < &f. The reactant concentration profile along the reactor can be calculated from the external fluid-phase mass balance by substituting the optimal activity distribution (2) dy/dy = -Daf(u(S))
(8)
The integration of this equation depends on the value of S,which is given by eq 4 as a function of y, and on the value of the parameter u, through eq 3. Therefore, we shall analyze separately the two cases a < 1 and u > 1. 2.1. u C 1. In this case S = 1 and u, = y, so that the heterogeneous model under examination reduces to a pseudohomogeneous model, i.e., f(u(S))= f ( y ) . Equation 8 can then be integrated analytically leading to the following implicit expression for y(y) Day = F ( y ) / ( l
+ a)2
(9)
where E ( y ) = -In y
where 6(s - g ( y ) ) is a Dirac delta function centered at a location S(y) within the pellet such that u(s) = u, (3) In other words, the optimal reactor outlet conversion is obtained by concentrating all the active catalyst in a specific location within the pellet support, which changes with pellet position along the reactor axis. (In the sequel, we shall refer to this as the rigorous optimal or type 1 catalyst activity distribution.) Obviously, such a situation is not very practical. Nevertheless, it is worthwhile to analyze this solution further, because it gives the upper bound of the reactor outlet conversion, with which the conversion value from any suboptimal distribution that
+ 2u(1-
y)
+ u2(1 - y2)/2
(9a)
2.2. u > 1; n = 1, 2. A t the reactor inlet, y = 1 > uM so that eq 4 and 5 hold. Since f(u(S))= &f is independent of the reactor axial coordinate, eq 8 leads to y = 1 - D&y
(10)
This solution holds from y = 0 to the value y = YM where y(y) = U M l / u . Beyond this point 701)C u M , SO that = 1and from eq 6 f(u(S))= f(y). Thus, the same procedure used above for the case u < 1may again be adopted. Using the initial condition y = uMat y = yM,we obtain D a b - YM) = [F(r)- F(UM)]/(~ + 0)'
(11)
which holds in the final portion of the reactor, y E [yM,l];
Ind. Eng. Chem. Fundam., Vol. 25, No.
3, 1986 315
Table I. Explicit Expression for the Bulk Reactant Concentration (y) and the Optimal Catalyst Location Profiles (1) along the Reactor Axis for the Slab Pellet (n = 0; Type 1 Distribution)" Q 2 > 4(u - 1)/(1 + U Y 62 < 4(u - 1)/(1 (TI2
+
Y E [YMI,~] y from eq 11 s=1
s o -
0 0
Yy from E [YM,1] eq 11
05
REACTOR AXIS, Y
S = l
"The indicated intervals of y actually exist only if the left extremum is smaller than the right one. bThe solution of this case is discussed in detail in part 4. C?
= 1. The final solution for y and S is reported in Table
I. 0
Da ~ 0 . 2Q; = 2 0
-
0
c.5
I
REACTOR AXIS, y
Figure 1. Rigorous optimal solution: profiles of the active catalyst location B and the external reactant concentration y along the reactor axis: (a) case of YM > 1; (b) case of y~ < 1.
note that this portion exists, obviously, only if Y M = 4(u - l)/Da(u 1)2< 1. With this, we can explicitly evaluate the optimal profile of the catalyst location within the pellet dong the reactor axis; viz., for y E [o,yM],S is given by eq 4 with y given by eq 10, while for y E [yM,l],s = 1. For illustrative purposes the calculated profiles of y and s along the reactor axis, for the case of Y M > 1and Y M < 1,are shown in Figure 1,a and b, respectively. Note that only the optimal location S depends on the Thiele modulus value, while the y profile remains unchanged-a feature already observed for the effectiveness factor in the case of a single pellet (Morbidelli et al., 1982). Moreover, it can also be seen that the optimal location s moves toward the pellet surface while moving downstream along the reactor axis. This is due to the diminished fluid-phase reactant concentration, y, which decreases the reactant concentration within the pellet and, thus, ita inhibition effect on the reaction rate. 2.3. u > I; n = 0. The case of the slab geometry (n = 0) is complicated by the fact that if condition 7 is satisfied, then the optimal location is fixed as B = 0. Noting that dt&,,2/dy < 0, it is seen that if condition 7 is violated at the reactor inlet (y = 0), then it will be violated in the entire range y E [0,1]. On the other hand, if it is satisfied at y = 0, then there might exist a position y = yminwhere 4= so that for y > yminwe have S > 0. The value of ymincan be obtained by using the solution of the reactor model for S = 0, reported in part 4, as the location along the reactor axis where the external fluid-phase concentration y reaches the value ymindefined by eq 7: ymin= [l+ @(1+ u)2/4]/u. In both cmes, when condition 7 is violated (i.e., 8 > 0), the same arguments adopted previously for the cylinder ( n = 1)and the sphere (n = 2) are applied to see if and when the optimal location reaches the pellet surface, i.e.,
+
3. A Suboptimal Catalyst Activity Distribution (Type 2) and Its Optimization As indicated above, the rigorous solution of the optimization problem is cumbersome to apply in practice; it is, therefore, desirable to approximate it with an active catalyst distribution which could be more easily implemented. It appears reasonable, at least as a first trial, to consider the case where all pellets in the reactor have the same activity distribution given by eq 2. (For ease of notation, this will be termed type 2 activity distribution in the sequel.) Also, this is indeed the rigorous optimal solution in some specific situations, particularly for low u values, and moreover it can be approximated reasonably well by a step distribution centered about the optimal location, S, and of a width not larger than ca. 5% of the pellet characteristic dimension (Morbidelli et al., 1982). Such activity distributions can be obtained for noble metals through coimpregnation with dibasic acids (cf. Michalko, 1966). Since the behavior of a reactor packed with nonuniform pellets of this type has been described fully in part 4, we shall focus on just the optimization problem, Le., to define the optimal location, S, which maximizes the reactor outlet conversion X . A comparison of conversion values thus obtained with those for the rigorous optimal distribution (type 1)described above will indicate the extent of need for further improvements of reactor performance. As described in detail in the Appendix, in order to maintain a one-to-one correspondence between B (or equivalently, j3, since the catalyst location S appears only in the parameter j3) and the outlet conversion, X , it is convenient to solve the optimization problem separately for three situations, depending on whether the inlet pellet operates on the low- or high-conversion branches or whether it lies in the region of unique solution for all j3 values, i.e., u C 8 (see Figure 3 of part 4). As shown in part 4,the reactor behavior when ignition does not occur is described by the following set of equations G(ii0) E iio + Bf(ii0) = 1 (12) Da = F(iil;ao) (13) X = 1 - ~ ( 1=) 1 - G(ii1) (14) On the other hand, if ignition does occur along the reactor length, then eq 13 is replaced by Da = F(ii*;iio)+ F(ii,;S*) F + F(iil;G*) (15)
316
Ind. Eng. Chem. Fundam., Vol. 25, No.
3, 1986
Table 11. Classification of Reactor Behavior for Given Values of u and Da (Tme 2 Distribution) u > 8, Inlet Particle on the Low-Conversion Branch (0 IPiopt5 H+) Da > Da*l Plopt= 0 Da*(PiOpt)= Da Dal < Da < Da*l Da, < Da < El Da(Biapt)= Da 0 < Da < Ea Piopt = H+ u
> 8, Inlet Particle on the High-Conversion Branch (Phapt2 H-) Da > Da*3 pho = H0 < Da < Da*3 Da*(PROpt) = Da u < 8 @opt 2 0) Da > Da*l or u < 1 Popt = 0 0 < Da < Da*l Da*(POpt)= Da
where recall that d* and d* are defined as follows (see Figure 2e of part 4). @(a*) = 0 (154
0.6 t
Do I
0.5
o
n
a’ W
m
5
0.4
z a w
-I
I
s
0.3
50 t 0.2 Do2
-Pop,=H+
0.1 0
0. I
0.2
Ht
--P Figure 2. Optimization of the parameter P as a function of Da for the low-conversion regime ( u = 20).
n
u = 20
Let us now first optimize X with respect to (3 and then trace back the optimal location, S. The optimal value, Po t, can be obtained by directly solving the equation dX/c!j3 = 0. However, the possible occurrence of ignition along the bed and the rigorous proof of existence and uniqueness of Bo@for a given set of parameter values make the analysis quite tedious. Therefore, we report here only the final results with brief comments about their practical application and relegate details of the analysis to the Appendix. I t is first convenient to introduce two functions of j3, Da*(@)and S(B), whose origin is clarified in the Appendix, with Da*(P) = P(fio;iio) (16) is that where do = l/do 2 and do is given by eq 12. value of Damkohler number, Da, such that In V ( i i o ) / f ( d l )+ ] T - In T - 1 = 0 (17)
s(@)
where T = aii*/[2(aii* - l ) ] and a*, do,and iil are given by eq 15a, 12, and 13 (or eq 15 in the presence of ignition). In particular, we shall refer to the following quantities Da*l = Da*(O); Da*2 = Da*(&); Da*3 = Da*(H-) (184 Dal = E(@-) = Da*2; E2= Da(H+) (lab) which are all functions of only the adsorption parameter, a.
The solution of the optimization problem is separately reported in Table I1 for the two cases a < 8 and a > 8. For the latter, the expression of the optimal B value is reported separately for the cases where the inlet particle operates on the low- or high-conversion branches, as and Oho t, respectively. If one compares the conversion values ogtained with B = $opt.and phOpt, it is readily seen which one of the two /?values is superior. In the sequel we refer to this value of @ as pope The use of Table I1 is quite straightforward; given a set of values for a and Da, we first evaluate the relevant Da*i and values from eq 18 and then directly select the correct expression for Popt. This is reported in either explicit or implicit form; in the latter case, a numerical solution is straightforward since uniqueness of Bopt is guaranteed. For the specific case of a = 20, optimum /3 values for the inlet particle operating on the low- and high-conversion branches are shown in Figures 2 and 3, respectively. It may be noted that the optimized condition for the reactor, in the case where the inlet particle operates on
si
High-Conversion Regime
2*10-~
P E( H - ,
m)
Figure 3. Optimization of the parameter 6 as a function of Da for the high-conversion regime ( u = 20).
the low-conversion branch, always requires ignition somewhere along the reactor when it is possible, i.e., for p E [P-,H+]. It is shown in the Appendix (section A1.2) that for B E [ / 3 3 + ] if no ignition phenomena occur along the bed, then the outlet conversion increases with 8, Le., dX/dp > 0; thus, in order to maximize X, /3 should be increased. Since for 8 = H+ ignition occurs for all Da values, it follows that under optimal conditions and E [p-JY+] ignition must occur along the reactor. 4. Optimal Catalyst Location for Type 2 Activity Distribution As mentioned above, the optimal catalyst location, S,can be readily calculated once the value of popthas been found by using the definition of p as follows:
s = 1 - Popt/r#?,
n =0
s = e~p[-28,,~/4~],n = 1
(194 (19b)
s = @/(@ + 3POpt), n = 2 (194 I t can be seen that, for n = 0 and @2 < pPpt,eq 19a leads to S < 0. In such cases the optimal location is given by 6 = 0. This is so because it can be shown that dX/d@> 0 for p < Popt (see Appendix), so that the maximum, X,is given by the maximum admissible value of Since dpop,/ds < 0, such a value is obtained for s = 0. In order to determine the optimal catalyst location, it is necessary to compare the reactor outlet conversion corresponding individually to the optimal /3 vr lues for the low- and the high-conversion regimes and then choose the particular p value which provides the larger conversion.
Ind. Eng. Chem. Fundam., Vol. 25, No.
3, 1986 317
x,q t
H.
H-
H,
-B
H,
-4 I Popt
=0
'
rDo*
Pop,= H-
LL W
im
1
0
n
m
I
t
!
I !
aopt
-a X
t -a I
-
IO 100 ADSORPTION P A R A M E T E R ,
1000
u
Figure 4. The master plot for evaluating the optimal fl value for type 2 distribution. Each number corresponds to a side figure which illustrates the X-j3 relationship for the low- and the high-conversion branches.
The final answer is summarized in the master plot shown in Figure 4. Here the u-Da plane is divided into several regions, each characterized by a particular expression for the optimal, Popt. In the smaller figures surrounding the main figure, qualitative behavior of the outlet reactor conversion as a function of P is shown for both the lowand the high-conversion regimes. From these it is clear which of the two regimes gives the optimal performance. Two additional characteristic functions of u have been added to those previously defined by eq 18: Da*4 = Da*(H+) and E3 = =(H-) (20) These are introduced in order to locate Poptwith respect to H+ and H-, respectively, thus leading to a splitting of regions - IV and VI. The argument is simple: since Da*(P) [Da(P)] is a decreasing function of P and POptis defined as Da = Da*(Popt)[Da = &(Popt)] (see Appendix), it follows that for Da > Da*4 (Da > Ea), POpt< H+ (Popt< HJ,and vice versa. This can be observed in the sketches shown around the master plot in Figure 4 for regions IVa,b and VIa,b. For example, in the case of region IVb this observation allows a conclusion that the optimal performance is obtained when the inlet particle is operating on the high-conversionbranch with @ = H-. This is because when, for a given P, both the ignited and the extinguished regimes are available, it is intuitively apparent that the former leads to a larger conversion. In conclusion, the master plot in Figure 4 provides a priori the optimal p value for any pair of u and Da values. The only exception to this is region IVa, where an a posteriori comparison of conversion from the individual optima of the low (i.e., = @lOpt) and high (Le., @ = H-) conversion regimes is required in order to establish the optimal distribution. The fact that in some cases better performance is obtained with the inlet particle on the low-conversion branch should not be surprising. For example, at large Da values
0
02
-
I 04
06
OB
I
OAMKOHLER NUMBER. D o
Figure 5. Optimal catalyst location as a function of the Damkohler number for a spherical pellet (a = 20; type 2 activity distribution).
(i.e., region 111),the external reactant concentration, y, drops quite rapidly after the reactor inlet, so that the rigorous optimal location (type 1distribution; section 2) approaches rapidly the external surface. Thus, most of the rigorously optimized reactor would indeed have S = 1, which is what the optimum for the low-conversion branch provides for the entire reactor. In Figure 5 the optimal catalyst location, B, is shown as a function of the Damkohler number, Da, for a spherical pellet (n = 2) and u = 20. These values are calculated according to Figure 4 by using the relationships reported in Table I1 and eq 19c. I t is clear that for increasing Da values the optimal location moves toward the pellet surface. Increasing Da leads to lower concentration in the fluid phase, thus diminishing the reactant inhibition effect on the reaction rate in a way similar to a direct decrease of the inhibition parameter, u. For a single pellet, we have previously noted that for decreasing u values the optimal catalyst location moves toward the pellet surface (Morbidelli et al., 1982). Note also that for Da = 0 the value of s obtained for the reactor is identical with that for a single pellet.
318
Ind. Eng. Chem. Fundam., Vol. 25, No. 3, 1986
1
$ 1 X
i
06L I
I
04
>
z 0
u
!
t 0.2i
,
a -
0
'
Figure 6. Sensitivity analysis of the optimized reactor performance to the value of 6 (Le., to the catalyst location S) (type 2 activity distribution).
In Figure 5, the discontinuity in the value of g is due to the transition from a situation where the inlet particle operates on the high-conversion branch to one where it operates on the low-conversion branch, i.e., from Phoptto P',, which occurs within region IVa (see Figure 4). Finally, it should be mentioned that, as for the rigorous optimal distribution (section 2), the optimal reactor performance (i.e.,the outlet conversion) depends solely on Da and u, while the Thiele modulus and the pellet geometry affect only the optimal location, S. 5. Discussion and Concluding Remarks
The sensitivity of the optimal outlet conversion to the /3 value has been investigated through the ratio $, defined
as
rc/
= L(P)/UP,,J
(21)
where L is the length of the reactor which, for a given P value, leads to the same conversion as the optimal one, L(@opt).Since all the particles contain the same total amount of active catalyst, the value of $ directly gives the ratio of the total amount of active catalyst to be used in the two reactors. From Figure 6 it is clear that $ = 1at P = Popt,while $ increases as the difference between and Poptincreases. From a practical point of view it is interesting to note that for the reactor small errors in the location of the active catalyst away from the optimal position do not have dramatic consequences on the outlet reactor conversion, as was indeed the case for the single pellet effectiveness (Morbidelli et al., 1982). In fact, from Figure 6 it appears that as Da decreases the reactor conversion becomes more and more sensitive to the difference between ,8 and Popt; specifically, for very small Da values, i.e., short reactors, (cf. Da = 0.1 in Figure 6), the same type of sensitivity as for the single particle is found. This is reasonable because in shorter reactors the averaging effect that a large number of particles provide is not available. Note that the discontinuities in the conversion curve are due to change of the reactor regime, which occurs at P = H-. It is now important to compare reactor performance between the optimized type 2 (section 4) and the rigorous optimal (type 1;section 2) distributions. In Figure 7 the optimal conversion values obtained from the two catalyst distributions are shown as functions of Da for spherical ( n = 2) or cylindrical ( n = 1)pellets. It is clear that as Da 0 or m , the two distributions lead to identical optimal outlet conversions, a conclusion which can be justified on physical grounds. As Da 0, the external fluid concentration, y , remains virtually constant along the reactor and
-
-
1
1
0
0 001
I
I
0.01
0.1
10
D A M K O H L E R NUMBER, D a
Figure 7. Comparison of reactor performance for the rigorous optimum (type 1) and optimized type 2 catalyst distributions as a function of the Damkohler number.
so does the rigorous location, 8, as given by eq 4. Thus, in this case, the rigorous optimal (type 1) distribution reduces to the optimized type 2 distribution. On the other hand, as Da m, y drops rapidly to zero after the reactor 0, and in most of the reactor ( Y M Iy inlet so that Y M 5 1)the rigorous optimal solution is S = 1. From Figure 4, the optimized type 2 distribution for Da m is also 8 = 1. For intermediate values of Da, the rigorously optimized distribution always leads to better performance. For example, suppose we want an outlet conversion, X = 0.90, then the rigorous and optimized type 2 distributions require Da = 0.165 and 0.41, respectively. This means a ratio of reactor lengths (and, therefore, of total catalyst) of about 2.5, thus indicating that at least in this intermediate range of Da values there may be other suboptimal catalyst activity distributions which might be better. Let us now establish the value of adopting the specific suboptimal distribution considered here (Le., type 2) by comparing it with usual catalyst distributions encountered in practice. To this aim, let us consider a reactor packed with pellets whose activity distribution is a step function (termed type 3 distribution, for notation) defined as follows
--
-
E [sJ] = 0, for s E [O,s,]
a(s) = a , for s
(22)
where s, = 1--2A,A is half-width of the activity step, so = 1 - A is the center of the step, and a = 1/(1- s;+l) in order to satisfy the usual activity integral condition (eq 8 of part 4).
The solution of the reactor model with this type of nonuniform pellet requires numerical calculation. The orthogonal collocation method (Villadsen and Michelsen, 1978) was used to discretize the intrapellet mass balance equations along the pellet coordinate, while a fourth-order Runge-Kutta method was used to integrate the resulting system of ODES along the reactor axis. The values of the outlet conversion, X, thus obtained are shown in Figure 8 as a function of the step distribution half-width, A. Note that for A = 0.5 and 0, the particularly common situations of uniform and externally coated pellets are obtained. In the latter case the reactor model reduces to a pseudohomogeneous one, and the value of X becomes independent of 4. It can be seen that for any given Thiele modulus value, 4, there exists a A = Aopt at which the reactor exhibits its optimal performance as long as step distributions of this type are considered. Moreover, unlike
Ind. Eng. Chem. Fundam., Vol. 25, No. 3, 1986 319
n=2;u=20 Do = 0.3
Rigorous Optimum ( Type I )
e
Optimized Type 2
I-I t
02
I
-
0
I
0 25
’,,’, ./’
\/;;’A
0.3
a
u (TypeI) ptimum W
Optimized Type 2
0.01
-
1 \ I
’ ; /
3
I.o
0.1 DAMKOHLER NUMBER,
Do
0
>
-
z
0.01
0 15
0 20
0 25
--P Figure 11. Regions for the occurrence of ignition along the reactor. For fi < fi- or @ > H+ ignition does not occur.
n:2;cr=20
Rigorous Optimum
w
I
\
0 IO
-
I
Figure 10. Sketch of the dimensionless reaction rate, f(u), as a function of concentration u ( u = 10).
/
n=2;u=20 A = 0.5
0.5 -U
Figure 8. Outlet conversion for a reactor packed with pellets having step activity distribution (22), as a function of the activity step half-width. IO
0
05
HALF-THICKNESS, A
0.I
IO
Do Figure 9. Comparison of reactor outlet conversion for rigorous optimal (type l), optimized type 2, and step activity (type 3) catalyst distributions; (a) A = 0.5; (b) A = 0.25. DAMKOHLER NUMBER,
the cases of rigorous optimal (type 1)and optimized type 2 catalyst distributions, the reactor conversion corresponding to A = Aopt.depends on 4. The outlet conversion from such a reactor is compared with that from an equivalent reactor with the optimized type 2 catalyst distribution in Figure 9a,b. Note that the curves for the former may exhibit a discontinuity in the slope due to occurrence of ignition phenomena along the reador. Since the behavior of the step distribution reactor is strongly affected by A and 4, it is difficult to draw general conclusions. In general, in order to obtain the m e conversion, the reactor with the optimized type 2 catalyst distribution requires a smaller total amount of catalyst by a factor between 1.2 and 2.5. This conclusion, however, does not hold for very large conversion (20.801,where all catalyst distributions become substantially equivalent for any Thiele modulus value. In this case, due to the significant variation of y along the reactor, the rigorous op-
timal catalyst location also changes significantly and ita approximation with a single fixed value can become poor. Note that in some particular situations the step distribution (type 3) is even better than the optimized Dirac delta distribution (type 2). In these cases, the rigorous optimal location (type 1 distribution) moves along the reactor in such a way that it happens to be approximated better by a step distribution that by a Dirac delta fixed at some location. In this context, note from Figure 7 that, particularly at high conversions, the gap between the rigorous and the optimized type 2 distributions is large, thus suggesting that better distributions may be found for higher conversions. In particular, a two-zone reactor with two different types of Dirac delta distributions or two step distributions with different centers and widths may well lead to significant improvements.
Acknowledgment We are grateful to the National Science Foundation (Grant No. CPE8308326) and to the Consiglio Nazionale delle Ricerche (Progetto Finalizzato Chimica Fine e Secondaria) for financial support of this work. Appendix The expressions reported in Table I1 for evaluating the optimal j3 value for type 2 activity distributions are developed here. Since for a given u value we optimize the outlet reactor conversion X with respect to the parameter 8, it is convenient to divide the optimization problem into
320
Ind. Eng. Chem. Fundam., Vol. 25, No.
3, 1986
subproblems, where the value of the outlet conversion is uniquely defined. For this, we consider separately the cases u < 8 and u > 8. In the former case (see Figure 3 of part 4), the outlet conversion value is unique for all p values. For u > 8, in order to have a one-to-one correspondence between P and X, we need to separately investigate two subcases. In the first, we limit the search for the optimal 6 value in the range E [O,H+],and for cases where the inlet particle exhibits multiplicity (i.e., E [H-,H+]),we take the low-conversion solution. In the second subcase, we limit the search to the range P E [H-,m], and if the inlet particle exhibits multiplicity, we now take the high-conversion solution. Therefore, for u > 8 we derive two optimal solutions for P (PIoptand Phopt), which we refer to as cases where “the inlet particle is operating on the low-conversion branch (or regime)” and on “the high-conversion branch”, respectively. The actual optimum, reported in Figure 4, for the case where a single Dirac delta distribution is used throughout the entire reactor is then determined by a direct comparison of the outlet conversion values given by these two p values. A l . u > 8. Inlet Particle on the Low-Conversion Branch. In part 4 it was shown that the inlet particle operates on the low-conversion branch (i-e.,large iiovalues) for 0 < p < H+ (in particular, multiplicity occurs when H< P C He) and that ignition along the reactor axis occurs when Da is sufficiently large and P- < p < H+. Let us first examine the region with no ignition phenomena. Al.l. P E [O,B-]. For a given set of values of u, p, and Da the outlet reactor conversion X is given by the system of eq 12-14. In order to optimize the outlet conversion with respect to 0,after a few manipulations we obtain
which vanishes at a value ii, = uosuch that f ( d o ) = f(iio), 1.e.
8, = l/u2iio
(24)
where the locations of iio and do are indicated in Figure 10. Since it can be shown that in the case under examination iio > u M , it can be concluded that 0 < do < U M < ii@ Moreover, since iil < iio,from Figure 10 it follows that
< do for do < iil < iio
f(iil)< f(iio), for iil
> f(iio),
(254
(25b) Introducing now, for any given P, the Damkohler value Da = Da*(/3) such that iil = do, i.e. from eq 13 f(ii,)
Da*(P) = F(do;iio) it can be seen from eq 23 that dX/d@ < 0 as Da > Da*
(16)
dF*/dp = In [ f ( i i O ) / f ( i i * ) ] - 1
(28)
where, since we are considering inlet particles operating in the low-conversion regime, we have ii* < iio. On the other hand since ii* corresponds to the minimum of G(ii), it is larger than iif = 2/u, which corresponds to the inflection point of G(ii),and it is thus also larger than uM = l / u . Thus, U M < uf< ii* < iio which from Figure 10 implies f(iio)< f(ii*),so that from eq 28 it can be concluded that F* decreases for increasing 0;thus F*(P-) > F*(H+)
(29)
where F ( H + )= 0, since iio = ii* at p = H+. The existence of F* implies ignition along the reactor. Now, from definition of F*, if Da < F*, ignition does not occur Da
> F*, ignition occurs
(30a) (30b)
Given values of and u, we can calculate F* as described in detail in part 4. For u = 20, the curve of F* vs. P is shown in Figure 11, where all the qualitative features discussed above can be seen. Now, if ignition does not occur, then the outlet conversion is again given by eq 12-14. However, since in this case U M < ii* f(ao),and thus from eq 23, dX/dp > 0. If ignition occurs, then from eq 12,14,15, and 28 it can be shown that dX/dp = -f(iil){ln [ f ( i i o ) / f ( i i l+) ]7 - In
7
- 1)
(31)
where
where the inequality comes recalling that a* > af = 2/u. Following similar arguments as in the previous section, for each P and u value, we can calculate through eq 12, 15, and 31 the Damkohler number value G(P)such that dX/dP = 0, i.e. In [ f ( i i O ) / f ( i i l ) ] + 7 - In
7
- 1= 0
(17)
It can be shown that Da decreases as P increases, so that
-
asDa=Da* (26) > O asDa Da*2. The curve of Da*(P) is shown in Figure 2 for a value of u = 20. Only the p E [O,&] region has been analyzed so far. For a given Da, it is evident from eq 26 that dX/dp > 0 ( %(H+) = &,
(32)
where note that El = Da*2. The optimal value of P in the range P E [P-,H+]is then given as a function of the Damkohler number by the inverse of the function &(@). By combining the results obtained in the last two sections, we can derive expressions for the optimal value PIopt in the range P E [O,H+J, for the case where the inlet pellets operate on the low-conversion branch; this result is summarized in Table 11. A2. u > 8. Inlet Particle in the High-Conversion Regime. The inlet particle operates on the high-conversion branch for p > H-, and it admits multiplicity for H< p < H+. In this case no ignition can occur, so that the optimization procedure is again based on eq 12-14 and 23.
321
Ind. Eng. Chem. Fundam. 1986, 25, 321-325
Let us introduce the new parameter
6 = 4(a - 1)/(1+ a)2
(33)
such that from eq 12 iio = uM. Thus, noticing that diio/dP < 0 and using eq 23, it appears that, for P > iil iio u M and dX/d@< 0, while, for < P, iio > U M and the same situation encountered in section A l . l is then obtained. Thus, introducizg Da*3 = Da*(H-) and Da*4 = Da*(O) = 0 (since at P = p, iio = uM, do = iio,and from eq 16 Da*4 = 0), the expressions reported in Table I1 for evaluating can be obtained. Note that, for Da = 0, Pwt = ,8 which is identical with the value obtained in part 1 (Morbidelli et al., 1982) for the case of a single catalytic pellet. The value of Phoptas a function of Da is shown in Figure 3 for a = 20. A3. a < 8. In this case all catalyst pellets along the reactor exhibit uniqueness, since y a 8, so that ignition cannot occur. Since the procedure is again the same as adopted in section Al.1, the derivation of the relationships reported in Table 11is not pursued any further. Note that, for a < 1, < 0 and Da*4is not defined any more. In this case, however, for any P > 0, iil < iio U M ,f ( i i 0 ) > f(ii,), and then from eq 23 dX/dP < 0. As reported in Table 11, it can then be concluded that, for a < 1, Popt = 0.
6,
PhF
a
Nomenclature The same symbols as in part 4 have been adopted in addition to the following: Da* = explicit function of 0, defined by eq 16 Da = implicit function of 8, defined by eq 17
M ( y ) = maximum value of f(u) M = f ( u M ) = (1 0)'/4a uf= inflection point of G(u),2/a
+
= maximum point of f(u),1/a a,, iil = value of u at s = s for the pellet at the reactor inlet and outlet, respectively bo = concentration value such that f ( i i o ) = f(b,J, l/iioaz X = outlet reactor conversion, 1 - y(1) PpPt= optimal 8 value j3opt, phopt= optimal 8 value in the case where the inlet particle is operating in the low- and high-conversion regime, respectively A = half-width of the step activity distribution defined by eq UM
22 92,$02 = Thiele modulus, c $ ~ = 902/(1+ aI2,+02 = kR,2/D, Literature Cited Cukierman, A. L.; Laborde, M. A.; Lemcoff, N. 0. Chem. Eng. Scl. 1983, 38. 1977. Johnson, D. L.; Verykios, X. E. J . a t e / . 1983, 79, 156. Johnson, D. L.; Verykios, X. E. A I C M J . 1984, 30, 44. Juang, H. D.;Weng, H. S.; Wang, C. C. I n New Horzons h Catalyds;Seiyama, T.. Tanabe. K., Eds.; Elsevier: Tokyo, 1981; p 866. Komiyama, M.; Merrlll. R. P.; Harnsberger, H. F. J . Catel. 1980, 63, 35. Melo, F.; Cervello, J.; Hermana, E. Chem. Eng. Sci. 1980. 35,2175. Michalko, E. U S . Patent 3 259 589, 1966. Morbldelil, M.; Servlda, A.; Varma, A. Ind. Eng. Chem. Fundam. 1982, 21, 276. Morbldelli, M.; Varma, A. Ind. Eng. Chem. Fundam. 1982. 21, 284. Morbidelli, M.; Servida, A.; Varma, A. Ind. Eng. Chem. Fundam., preceding paper In this issue. Oh, S. H.; Cavendish, J. C.; Hegedus, L. L. AIChE J . 1980, 26, 935. Shyr, Y. S.; Ernst, W. R. J . Catal. 1980, 63, 425. Villadsen, J.; Mlchelsen, M. L. Solution of Dlffmntlal EquatEOn Models by Polynomial Approximation ; Prentictt-!Haii: Englewood Cliffs. NJ, 1976.
Received for review September 19, 1984 Accepted August 23, 1985
Analysis of a Batch Adsorber with Rectangular Adsorption Isotherms Duong D. Do Department of Chemlcal Englneering, University of Queensland, St. Lucia, Queensland 4067, Australla
This paper deals with an analysis of a batch adsorber with a rectangular adsorption isotherm ( i a , a system operating at the flat part of a Langmuir isotherm). This situation is common in many systems of adsorption of organic solutes from aqueous solution onto activated carbon. The theoretical solutions reveal that the transport parameters, the effective diffusivity and the extemal mass transfer coefficient, can be easily extracted from two simple linear plots. Two systems of adsorption of phenol onto activated carbon and Victoria Blue dye onto activated carbon are employed to demonstrate the versatility of the theoretical solutions.
Introduction Activated carbon is one of the most commonly used adsorbents in the petrochemical industry (Barneby, 1971) and in the wastewater treatment area Recently, researches have been focused on removing organic solutes from aqueous streams with activated carbon as the adsorbent medium (Fritz and Schlunder, 1981a; Fritz et al., 1981b,c; Sheindorf et al., 1982; Neretnieks, 1976;Peel and Benedek, 1981; Van Vliet et al., 1980). Activated carbons are also being used by some researchers to remove organics from oil shale retort water (Pedram et al., 1982). Batch adsorbers are very often used to obtain the transport parameters for such systems, namely, the effective diffusivity and the mass transfer coefficient. Traditionally, an adsorption isotherm curve is obtained 0196-4313/86/1025-0327~Q1.SO10
from an equilibrium study and the concentration of the adsorbate is chosen such that the system is operating in the linear part of the isotherm. In this way linear theory can be applied and solutions obtained by Crank (1975) can be used to extract the effective diffusivity. Unfortunately, in the adsorption of organic solutes by activated carbon, there is only a very small Concentration range within which the isotherm is linear. Therefore, the experimental error will be high if a very low concentration is selected for use in the experiment. Moreover, the concentrations of the organic solutes to be treated in practical systems can be well outside the linear range. In fact, they lie in the range in which the solid-phase concentration is almost constant. It is the objective of this paper to present an analysis of a batch adsorber having a rectangular adsorption iso0 1986 American Chemical Society
