A finite element solution for diffusion and reaction in partially wetted

A finite element solution for diffusion and reaction in partially wetted catalysts with power-law kinetics. Patrick L. Mills, Steven Lai, Milorad P. D...
1 downloads 0 Views 536KB Size
Znd. Eng. Chem. Res. 1988,27, 191-194 Registry No. Fe(OH),, 18624-44-7.

Literature Cited Ames, L. L.; Salter, P. F.; McGarrah, J. E.; Walker, B. A. “Selenate-Selenium Sorption On a Columbia River Basalt“. Chem. Geol. 1979,43, 287. Baldwin, R. J.; Acres, W.; Stauter, J. C.; Terrell, D. L. U S . Patent 4 405 464, 1983. Chau, Y. K; Riley, J. P. Anal. Chim. Acta 1965,33, 36-49. Cooper, W. C.; Westbury, R. A. In Selenium; Zingaro, R. A., Cooper, W. C., Ed.; Van Nostrand Reinhold: New York, 1974; pp 109-111. De Carlo, E. H.; Zeitlin, H.; Fernado, Q. Anal. Chim. 1981, 53, 1104-1107. Doran, J. W.; Alexander, M. “Microbial Transformations of Selenium”. Appl. Enuiron. Microbiol. 1977, 33, 31-37. Doumas, A. C. US.Patent 3387928, 1968. Feigl, F.; Anger, V. Spot Tests in Inorganic Analysis, 6th ed.; Elsevier: New York, 1972; pp 408-409. Gersberg, R. M.; Brenner, R.; Elkins, B. V. “Removal of Selenium Using Bacteria”. California Agricultural Technology Institute: 1986; Pub. No. CATI/860201. Handbook of Chemistry and Physics, 57th ed.; CRC: Boca Raton, FL, 1976-1977; pp D-67-D-77. Hartmanis, M. G.; Stadtman, T. C. “Isolation of a Selenium-Containing Thiolase From Clostridium Kluyuen’: Identification of the Selenium Moiety as Selenomethionine”. Proc. Natl. Acad. Sci. U.S.A. 1982, 79, 4912-4916. Hingston, F. J.; Atkinson, R. J.; Possner, A. M.; Quirk, J. P. “Specific Adsorption of Anions”. Nature (London) 1967,215,1459-1461. Hingston, F. J.; Possner, A. M.; Quirk, J. P. “Adsorption of Selenite by Geothite”. Adu. Chem. Ser. 1986, 70,82-90.

191

Jackson, T., Research Leader, Fish and Wildlife Service, Denver Field Office, Columbia National Fisheries Research Laboratory, unpublished results, 1986. Jennings, P. H.; Yannopoulos, J. C. In Selenium: Zingaro, R. A., Cooper, W. C., Eds.; Van Nostrand Reinhold New York, 1974; pp 48-59. Kauffman, J. W.; Laughlin, W. C.; Baldwin, R. A. U S . Patent 4 519 912, 1985. Lee, G. F.; Stumm, W. “Determination of Ferrous Iron in the Presence of Ferric Iron With Bathophenantholine”. J . Am. Water Works Assoc. 1960,52, 1567-1574. Maneval, J. E.; Klein, G.; Sinkovic, J. “Selenium Removal From Drinking Water by Ion Exchange”. Municipal Environmental Research Laborator.7, Office of Research and Development, EPA-81-02-5401, Cincinnati, OH, 1983. Marchant, W. N. US.Patent 3 933 636, 1976. Oppenheimer, J. A.; Eaton, A. D.; Kreft, P. H. “Speciation of Selenium in Groundwater”. Water Engineering Research Laboratory, EPA-600/S2-84-190, Cincinnati, OH, 1985. Schwertmann, U.; Taylor, R. M. In Minerals in Soil Enuironments; Dixon, J. B., Weed, S. B., Eds.; Soil Science Society of America: Madison, WI, 1977; pp 167-172. Siu, K. W.; Berman, S. S. “Determination of Selenium(1V) in Seawater by Gas Chromatography after Coprecipitation with Hydrous Iron(II1) Oxide”. Anal. Chem. 1984, 56, 1806-1808. Sorg, T. J.; Logsdon, G. S. “Treatment Technology to Meet the Interim Primary Drinking Water Regulations for Inorganics: Part 2”. J. Am. Water Works Assoc. 1978, 70, 379-393. Received for review September 17, 1986 Revised manuscript received August 25, 1987 Accepted October 3, 1987

COMMUNICATIONS A Finite Element Solution for Diffusion and Reaction in Partially Wetted Catalysts with Power-Law Kinetics The Galerkin finite element method (FEM) is used to obtain approximate solutions to the problem of diffusion and reaction in partially wetted catalyst pellets for power-law reaction rate forms. For the case of a first-order reaction, it is shown that effectiveness factors developed from the FEM solution are in good agreement with those obtained from an integral equation and a least-squares solution. The effects of reaction order and other associated parameters on the concentration profiles and catalyst effectiveness factor are briefly illustrated. Fixed beds of solid catalyst with gas-liquid cocurrent flow are commonly used in commercial processes such as hydroprocessingof various petroleum feedstocks, selective partial hydrogenation of chemicals, and oxidation of dilute contaminants that may be present in aqueous process streams (Germain et al., 1979; Ramachandran and Chaudhari, 1984). Under certain process conditions for a particular reacting system, the catalyst pores can be either totally or partially liquid-filled and the catalyst surface can be covered by a combination of actively flowing liquid films and stagnant liquid films or contain dry patches that are in direct contact with the gas or possibly volatile liquid (Mills and DudukoviE, 1980). The ability to evaluate the catalyst effectiveness factor for these situations is an important step that should be performed when designing new reactors or during interpretation of performance for existing reactor systems (cf., Beaudry et al. (1985)). Over the past decade, various models for evaluation of catalyst effectiveness factors for partially wetted catalysts have 0888-5885/88/2627-0191$01.50/0

been developed (DudukoviE, 1977; DudukoviE and Mills, 1978; Goto et al., 1981; Herskowitz et al., 1979; Herskowitz, 1981a,b; Lee and Smith, 1982; Levec et al., 1980; Martinez et al., 1981; Mills and DudukoviE, 1979, 1980; Mills et al., 1981; Ramachandran and Smith, 1979; Sakornwimon and Sylvester, 1982; Tan and Smith, 1980). With only one exception (Goto et al., 1981),these models assumed linear first-order reaction kinetics which limits their application since many systems having practical significance as described by nonlinear kinetics (cf., Ramachandran and Chaudhari, 1984). For the above case where nonlinear kinetics were assumed, a finite-difference scheme with a fixed grid size was used to solve the resulting split boundary-value problem for regular slab geometry. This scheme is not necessarily the best choice for this type of problem since detailed studies on related problems described by Laplace’s equation have shown that the solution accuracy can be quite low in the vicinity of the split boundary discontinuity (cf., Whiteman, 1968). Steep 0 1988 American Chemical Society

192 Ind. Eng. Chem. Res., Vol. 27, No. 1, 1988

concentration profiles that occur for certain ranges of parameters require the use of fine mesh sizes to place a sufficient number of grid points within the profile (Goto et al., 1981; Loffler and Schmidt, 1975). Finally, the diffusional characteristics of catalyst particles of irregular geometry such as trilobes, stars, macaroni, and wagon wheels, which are patented shapes used in residuum hydrotreating (Dautzenberg, 1984), cannot be readily simulated by finite-difference methods. The objective of this paper is to obtain solutions for diffusion and nonlinear reaction in a partially wetted catalyst by the Galerkin finite elements method (FEM) since it is ideally suited for nonlinear problems that are posed on either regular or irregular geometry. Since the FEM has not been applied to the particular nonlinear split boundary value problem treated here, this represents a novel application of the method and an extension of our previous work on linear diffusion and diffusion-reaction systems (Mills et al., 1987).

Model Equations If one assumes the catalyst surface is partially wetted and slab geometry is used, then the governing equations in dimensionless form for the case of power-law reaction kinetics are (s/W)za2u/ap

5 =o, = 1, 7 = 0,

- l/Biw

= 0, - i/Bid

+ a2U/aT2

- 4 2 ~ a=

o

a u / a g = 0;

oI

I1

(2)

a u / a g = 1;

0 1 tl I1

(3)

+ u = 1; 0 I5 e p/w aulaT + u = p/w c 5 I1

= 1,

(1)

au/atl = 0;

o I4 1

1

(4a) (4b) (5)

The value u represents the dimensionless reactant concentration u = C/Cb, while the parameters Bi, and Bid are the Biot numbers for the actively wetted and inactively wetted or dry catalyst surfaces, respectively. The parameter E denotes the ratio of bulk concentrations at the inactively wetted and actively wetted surfaces, i.e., t = C,/Cw. The aspect ratio of the pellet is denoted by 6/w,while the fraction of actively wetted external catalyst surface is given by plw. It is assumed that E = 1and 6/w = 1for simplicity so that the bulk concentration of reactant is uniform and the pellet has a square shape. One quantity of interest is the catalyst effectiveness factor which corresponds to the ratio of the observed rate to the intrinsic or maximum possible rate. It is given by

Mehta and Aris (1971) developed solutions to eq 1-5 for the limiting case of p/w = 1corresponding to the totally wetted or totally dry catalyst pellet. These solutions are useful for obtaining approximate values for the effectiveness factor of partially wetted catalysts by the weighting factor method of Ramachandran and Smith (1979). Solution by Galerkin Finite Elements The use of the FEM for solution of diffusion-reaction equations with discontinuous boundary conditions and linear first-order reaction kinetics has been described in a previous publication (Mills et al., 1987). An extension of this approach to the nonlinear power-law kinetic scheme given in eq 1 is briefly outlined below. The Galerkin finite-element method requires that the residual of eq 1be orthogonal to each trial function N Lfor

the whole area of the element. Application of the Green Gauss theorem to eq 1 gives I1

I

I = i j , k (7 I11 A nonlinearity appears in term I11 of eq 7 which originates from the nonlinear reaction rate form. This term can be linearized, but an iterative procedure has to be adopted to find the finite-element solution. Two methods for linearization of the nonlinear term were examined. Method 1 simply calculates the nonlinear part by using the previously determined values, while method 2 uses the Taylor series expansion (including only first-order terms) about the previous value as method 1: method 2: Ui+l(e)"

+ [Ui+l(e) - Ui(e)~aUi(e)"-'= (1- a ) u p n+ a[Ui(e)a-']Ui+l(e) (9)

where i means the previous value and i + 1 means the current value. The results given here are based on method 2. A more detailed comparison of both methods is given elsewhere (Mills et al., 1987). Substitution of the expression for the dependent variable given by eq 16 in our previous publication (Mills et al., 1987) and its derivatives into eq 7 gives the following set of linear equations after performing a few elementary operations: [K(e)( u ~ ( ~[u) ~) ]+ ~ = ( ~[fie)( )] u~(~))] (10) where the matrix K(e)is given by the sum K(e)(~i(e)) = Klfe) K2(e) K3(e)(~i(e)) (11)

+

+

and the matrix P ( e )is given by the sum fie)(ui(e)) = p2(d- p 3 ( e ) ( U i ( e ) )

(12)

The matrix K3(e)(ui(e))and the matrix P3(e)(~ife)) result from term I11 and are given by

for method 1 and

p 3 ( e ) ( U i ( e ) )=

42

1L(e)(l ui(e)"N dA - a)

(16)

for method 2 For the case of a first-order reaction, P3(e) = 0 and K3(e)= @ 11 N NT dA are obtained as limiting values from both iterative schemes. The various matrices that appear above in eq 11-16 such as Kl(e),K2(e),and Pz(e)are defined by eq 23, 26, and 27, respectively, in our previous publication and will not be repeated here. Equation 10 can be reduced to the following linear system of banded equations:

KNNXNNi= uNNi+' = PNN' (17) where NN is the number of nodes. Once an initial estimate for the concentration profile at the various nodes is assumed, the current values for the concentration ui+' are

Ind. Eng. Chem. Res., Vol. 27, No. 1, 1988 193 Table I. Comparison of Effectiveness Factors Obtained by Three Methods. Parameters: Bi, 0.5, a = 1

6 1

5 20 50 a

exact 6.5470 (-1) 1.1967 (-1) 2.6220 (-2) 1.0207 (-2)

effectiveness factor, TTB least squares, integral eq, finite element, N = 200 N=6 NT = 150 6.5361 (-1) 6.5311 (-1) 6.5705 (-1) 1.1910 (-1) 1.1892 (-1) 1.2081 (-1) 2.6040 (-2) 2.6038 (-2) 2.6510 (-2) 1.0131 (-2) 1.0152 (-2) 1.0349 (-2)

- m,

Bid

0,6/w = 1, p / w =

% relative errorn

least squares, N = 200 0.166 0.476 0.686 0.539

integral eq, N=6 0.243 0.627 0.694 0.745

finite element, NT = 150 0.359 0.953 1.106 1.391

Defined as 100 ITTB(approximate) - vTB(exact)l/qTB(exact).

'0

a

C

'0 '0

b

d

Figure 1.

determined by the Choleski decomposition method. The iteration is continued until the relative Euclidean norm or total residual (Carey and Oden, 1985) is less than a specified error criterion, e.

Results and Discussion Figure 1 gives a perspective view of the dimensionless concentration profiles obtained by the FEM for Bi, m and Bid = 0 at reaction orders of a = 'Iz and a = 2 and Thiele moduli of 4 = 1 (parts a and b of Figure 1)and 4 = 5 (parts c and d of Figure 1). The fraction of externally wetted catalyst surface was plw = 0.5. This choice of Biot numbers corresponds to the case of a nonvolatile liquid limiting reactant as might be encountered, for example, in residuum hydroprocessing. In the viewing angle given here, the center of the pellet is nearest the observer and the reactant concentrations lie on the vertical plane formed at 17 = 1for 0 If I1. The reactant concentrations on the external surface of the pellet are on the vertical plane at 7 = 0 for 0 If < 1 and are farthest from the observer. On the external surface, u = 1 for 0 If < 0.5 which corresponds to the actively wetted fraction where reactant enters the pellet interior. Since the reactant is nonvolatile, &/at = 0 over the inactively wetted fraction of the catalyst surface, corresponding to 0.5 < E I1. No reactant enters the pellet on this surface, and the reactant concentration sharply decreases with increasing distance from the wetted boundary. Additional inspection of Figure 1shows that for a given value of the Thiele modulus, an increase in the reaction to a = 2 produces a lesser demand for order from a = 'Iz the reactant so that the centerline concentration of reac-

-

tant increases. For a given reaction order, increasing values of the Thiele modulus results in a greater demand for the reactant so that the concentration profiles are steeper. These results are intuitively expected, but Figure 1provides insight into the results for a partially wetted catalyst. Particularly noteworthy is that the FEM produces smooth profiles with no spurious behavior near the discontinuous boundary. Results at larger values of the modulus, such as 4 = 50 (not shown for brevity), have a similar behavior. Profiles for the case where both Bi, and Bid are finite corresponding to volatile liquid-limiting or gaseous reactant limiting reactions are given elsewhere (Mills et al., 1987). As shown earlier by eq 6, the volume beneath the uu(f,7) surface for 0 I f I1 and 0 Iq I 1 corresponds to the catalyst effectiveness factor. Table I gives a comparison of catalyst effectiveness factors obtained by three independent methods for the case of first-order kinetics. The techniques used to obtain results by the least-squares and integral equation methods have been previously described (Mills et al., 1987). These are based upon solutions of the dual-series equations derived from the separation of variables solution of eq 1-5. The exact results were obtained by an extrapolation procedure which uses the least-squares results as the basis (Mills et al., 1987). Inspection of the table shows that all three methods produce relative errors that are within 0.1-1% of each other. The errors increase with increasing values of the Thiele modulus, since the concentration gradients are more steep and the degree of approximation for each method is fixed. The magnitude of the errors associated with the finite-element method are within reasonable limits for application purposes. Proper selection of the triangle density for the finite-element method will ensure that the relative errors will fall within acceptable limits. One of the key assumptions used in the above model is that the external surface of the catalyst is contacted by flowing liquid over the surface fraction 0 I$. < p / w ,while the remaining surface fraction plw < E I 1 is covered by a stagnant liquid film or is in direct contact with the gas. Furthermore, the external wetted fraction plw or so-called external contacting efficiency is assumed to be constant and not fluctuate in time. Under appropriate reaction conditions, it is possible that the local liquid flow hydrodynamics on the catalyst particle surfaces change in such a way that the external contacting efficience plw is not a constant, but varies with respect to both time and the position of the catalyst particle in the reactor. In this case, it would appear that a complete description of catalyst wetting dynamics on the particle scale is necessary before the catalyst effectiveness factor can be locally evaluated. Experimental methods that allow local values for the external contacting efficiency to be determined as a function of reactor operating parameters have not yet been developed. In addition, a model which is based upon first principles that permits an a priori prediction of the local catalyst external wetting is also not available. It would be useful to know, however, how the performance of a

194 Ind. Eng. Chem. Res., Vol. 27, No. 1, 1988

trickle-bed reactor is affected by the catalyst effectiveness factor when this latter quantity is calculated using a constant value for the external catalyst wetting efficiency versus one that accounts for a distribution of external contacting efficiencies within the catalyst bed. To obtain an idea on whether or not a complete description of wetting on the particle scale is necessary, consider a simple model for prediction of trickle-bed reactor performance that assumes plug flow of liquid, isothermal operation, uniform global liquid distribution, and first-order reaction with respect to a nonvolatile liquid limiting reactant. If the fraction of catalyst particles in the reactor having an external contacting efficiency between p / w and p / w + d ( p / w ) is represented by a density function f ( p / w ) d ( p / w ) then , the reactant conversion will be XB = 1 - exp

[

-DUO

qTB('#',P

/ w , * * . ) f (/pw )d ( p/ w )

]

(18)

where Duo = p (1- EB)k&R/UL is the Damkoehler number. Since aqTB/&/w) 2 0 and a2qTB/a2(p/w)< 0, i.e., the overall catalyst effectiveness factor for the above case is a monotonic, increasing function of the external wetting efficiency p l w , then

This implies that the catalyst effectiveness factor evaluated __ at an average external wetting efficiency q o ( p / w )is less than the average effectiveness factor. To obtain a lower bound or conservative estimate of reactor performance, only an average value for the external wetting efficiency is necessary from which an effectiveness factor a t this average can be evaluated. The reactant conversion will then be given by XB 5 1 - e-DQO9TB (20) The above development, while not necessarily general for all cases of interest, suggests that use of an average value for the catalyst external contacting efficiency p / w is justified. When reliable values for the local external wetting efficiency become available, a more rigorous approach that accounts for the local variation in the effectiveness factor can be considered.

Summary and Conclusions The Galerkin finite-element method was used to obtain numerical solutions for a novel type of nonlinear mixed or split boundary-value problem that describes diffusion and reaction in partially wetted catalyst pellets with power-law kinetics. For the case of linear first-order kinetics, the FEM results were in good agreement with results obtained by two other independent methods. Results for other reaction orders were shown to give the correct trends. Investigation of novel catalyst shapes having irregular geometry as part of a larger study to design optimal catalysts for trickle-bed reactor applications represents a possible application of this work. Nomenclature A(")= area of a particular finite element first appearing in eq 7

Bid, Bi, = Biot numbers for the dry and wetted catalyst surfaces, respectively e = domain of a particular finite element first appearing in eq 7

K ( e ) ( ~ ; (=e )linearized ) finite-element matrix first appearing

in eq 10 and defined by eq 11

Kl("),K2(e),K3(e)= matrices that occur in the finite-element method first appearing in eq 11 KNNxNNi = diagonal matrix appearing in eq 17 Nl = shape function first appearing in eq 7 p = x coordinate where the boundary condition discontinuity occurs, cm P(e)(ui(e)) = linearized finite-element matrix first appearing in eq 10 and defined by eq 12 P 2 ( eP ) , 3 ( e ) ( ~= i (matrices e)) that occur in the finite-element method first appearing in eq 12 u = dimensionless reactant concentration first appearing in eq 1 de)= dimensionless reactant concentration in a particular finite element first appearing in eq 7 w = pellet half-width, cm Greek Symbols a = reaction order 7 = dimensionless y coordinate, y l w qTB = catalyst effectiveness factor for a partially wetted pellet

defined by eq 6 6 = pellet thickness, cm

5 = dimensionless x coordinate, x/6

4 = Thiele modulus defined as 6(hb/D,)1/2,dimensionless

Literature Cited Beaudry, E.; DudukoviE, M. P.; Mills, P. L. "Diffusional Effects for Gas Limited Reaction in a Trickle-Bed Reactor", 84th National AIChE Meeting, Chicago, Nov 10-15, 1985; paper 131c. Carey, G. F.; Oden, J. T. Finite Elements: A Second Course; Prentice-Hall: Englewood Cliffs, NJ, 1985. Dautzenberg, F. M. "Advances in Hydrotreating Technologies", paper presented at the Catalytica Associates Annual Science and Technology Seminar, Santa Clara, CA, April 15-17, 1984. DudukoviE, M. P. AZChE J . 1977, 23, 940. DudukoviE, M. P.; Mills, P. L. Adv. Chem. Ser. 1978, 65, 387. Germain, A,; L'Homme, G.; Lefebvre, A. In Chemical Engineering of Gas-Liquid-Solid Catalyst Reactions; L'Homme, G. A., Ed.; CEBEDOC: Liege. 1979. Goto, S.; Lakota, A ;: Levec, J. Chem. Eng. Sci. 1981, 36, 157. Herskowitz, M. Chem. Eng. Sci. 1981a, 36, 1665. Herskowitz, M. Chem. Eng. Sci. 1981b, 36, 1099. Herskowitz, M.; Carbonell, R. G.; Smith, J. M. AIChE J. 1979,25, 272. Lee, H. H.; Smith, J. M. Chem. Eng. Sci. 1982, 37, 223. Levec, J.; Pavko, S.; Dobovisek, M. Chem. Eng. Sci. 1980, 35, 180. Loffler, D. G.; Schmidt, L. D. AZChE J . 1975,21, 786. Martinez, 0. M.; Barreto, G. F.; Lemcoff, N. 0. Chem. Eng. Sci. 1981, 36, 901. Mehta, B. N.; Aris, R. J . Math. Anal. Appl. 1971, 36, 611. Mills, P. L.; DudukoviE, M. P. Znd. Eng. Chem. Fundam. 1979, 18, 139. Mills, P. L.; DudukoviC, M. P. Chem. Eng. Sci. 1980, 35, 2267. Mills, P. L.; Erk, H. F.; Evans, J.; DudukoviE, M. P. Chem. Eng. Sci. 1981, 36, 950. Mills, P. L.; Lai, S.; DudukoviE, M. P.; Ramachandran, P. A. Comp. Chem. Eng. 1987, in press. Ramachandran, P. A,; Smith, J. M. AZChE J . 1979,25, 538. Ramachandran, P. A.; Chaudhari, R. V. Three Phase Catalytic Reactors; Gordon and Breach: New York, 1984. Sakornwimon, W.; Sylvester, N. D. Znd. Eng. Chem. Process Des. Deu. 1982, 21, 16. Tan, C. S.; Smith, J. M. Chem. Eng. Sci. 1980, 35, 1601. Whiteman, J. R. Q. J . Mech. Appl. Math. 1968, 21, 41. t

Monsanto Company. Washington University.

Patrick L. Mills,*+ S t e v e n Lai* Milorad P. DudukoviE; Palghat A. R a m a c h a n d r a n f Monsanto Company S t . Louis, Missouri 63167 and Department of Chemical Engineering Washington University S t . Louis, Missouri 63130 Received for review July 17, 1986 Accepted September 9, 1987