9712
Ind. Eng. Chem. Res. 2005, 44, 9712-9717
Self-Oscillation Liquid Front inside a Partially Wetted Catalyst Pellet under r-Methylstyrene Hydrogenation: Experiment and Theory Alexey B. Shigarov,* Alexander V. Kulikov, and Valery A. Kirillov Boreskov Institute of Catalysis, Pr. Acad. Lavrentiev 5, Novosibirsk 630090, Russia
Self-oscillations on a partially wetted cylindrical catalyst pellet 1% Pd/γ-Al2O3 in R-methylstyrene (AMS) hydrogenation accompanied by liquid evaporation are studied. By NMR tomography (magnetic resonance microimaging), images of the liquid-phase distribution within the porous pellet are obtained and the liquid-phase redistribution is monitored immediately during the process. The physical mechanism and simplified one-dimensional model are suggested to explain the observed phenomenon. The model predicts well the experimental period (2-3 min) and the amplitudes of the temperature oscillations at the lower (dry) part of the pellet and also the liquid front oscillations at the upper part. The instability mechanism (negative feedback with time lag) is based on complex interaction among the processes of imbibition, evaporation, and gas-phase hydrogenation. Introduction It is known that interaction of physical processes (imbibition, evaporation, and heat and mass transfer) with exothermic reaction can result in the formation of various types of instability and critical phenomena in multiphase catalytic systems. The investigation of such regimes is interesting from both academic and industrial points of view with regards to hot spots and runaway phenomena in trickle-bed reactors. Liquid hydrocarbon hydrogenation presents an important example on this item. It is worth mentioning here the early research,1 where the crucial roles of gas-phase reaction and moving hot spots were experimentally found in a laboratory trickle-bed reactor under R-methylstyrene (AMS) hydrogenation. Slow progress in the understanding of various intriguing phenomena observed in laboratory trickle-bed reactors is partly explained by complex hydrodynamics of multiphase flow. So, the focus of research on a single particle scale seemed promising. For example, AMS hydrogenation was used to demonstrate the acceleration of drying on a single catalyst pellet.2 The multiplicity of catalyst pellet liquid filling and temperature steady states under cyclohexene hydrogenation in a single-pellet reactor was reported.3 Ignition and extinction phenomena on different types of catalyst pellets under variation of the liquid AMS feed were experimentally studied.4,5 Theoretical considerations are rather scarce. Careful analysis of the model of a half-wetted and partially liquid-filled catalytic porous slab was made6 but without its experimental verification. Physical mechanism and mathematical models were proposed7 that successfully explained the hysteresis behavior observed earlier4,5 by the same authors. A novel promising approach to the experimental study is application of in situ NMR imaging,8 accompanied by accurate intrapellet thermocouple measurements. This method allows one to obtain a picture of the liquid distribution inside the porous pellet and to in situ monitor its redistribution in time. * To whom correspondence should be addressed. E-mail:
[email protected].
The first brief report on this item appeared recently.9 Further development of this approach enabled one to find experimental self-oscillations on a single cylindrical catalytic pellet under AMS hydrogenation.10 The goal of our current work is coupling of experimental and theoretical approaches in order to explain this unusual self-oscillating behavior. Experimental Setup As a model reaction, we chose AMS hydrogenation:
C9H10 + H2 ) C9H12 This reaction is widely used as a model reaction in the experimental investigation of processes in trickle-bed reactors. Figure 1 presents a schematic of the singlepellet reactor unit we used. This setup contained a reactor consisting of three vertical coaxial glass tubes. In the inner tube (9 mm i.d.), a catalyst pellet was placed and blown with a preheated hydrogen flow. The 10-mm-high cylindrical catalyst pellet of 5 mm diameter contained 1% Pd/γ-Al2O3. The pellet was kept by a chromel/copel thermocouple of 0.2 mm diameter. The earlier experiments9 revealed that a thermocouple itself significantly distorts NMR microimages of the liquid distribution inside a porous catalyst. So, we were forced to place the thermocouple at the bottom part of the pellet in order to avoid distortions at least for the rest pellet volume. The thermocouple was introduced into a carefully drilled channel in the lateral surface of the cylinder. The channel was at a distance of ∼6 mm of the pellet height from the upper base of the pellet. The depth of the channel was equal to the cylinder radius. The gap between the inner and middle tubes was used to control the reactor temperature by blowing with a hot air flow. The gap between the middle and outer tubes was evacuated to prevent the NMR tomograph from being heated. Liquid AMS was fed to the upper base of the pellet through a glass capillary of 0.2 mm i.d. Simultaneously, the catalyst pellet was blown with a pure hydrogen flow. The temperature of the gas
10.1021/ie050269x CCC: $30.25 © 2005 American Chemical Society Published on Web 05/26/2005
Ind. Eng. Chem. Res., Vol. 44, No. 25, 2005 9713
Figure 1. Schematic of a single pellet reactor unit placed in the probe of the NMR tomograph.
Figure 3. Experimental catalyst pellet temperature oscillations and NMR microimages of liquid distribution inside the pellet at discrete moments of time (Gliq ) 0.8 mg/s).
Figure 2. Steady-state hysteresis curve for AMS hydrogenation on a single catalyst pellet (from ref 7) with the region of sustained self-oscillations.
(hydrogen) flow was measured with a second thermocouple, which was placed upflow at 8 mm from the pellet top. The experimental setup also included hydrogen and liquid AMS supply systems, temperature control and measurement units, and air heaters for controlling the reactor temperature. All of the experiments were carried out at atmospheric pressure, a hydrogen flow rate of 10.9 cm3/s, and a gas and liquid inlet flow temperature of 87 °C. The rate of feed of liquid AMS onto the pellet was 0.8 mg/s. The original details of the NMR-microimaging technique can be taken from a previous report.10 Results of Experiments At the first time, self-oscillations under AMS hydrogenation in a single-pellet reactor were experimentally found by temperature measuring (without the use of NMR imaging). In those experiments, the cylindrical catalytic pellets were half as short (D ) H ) 5 mm), and thermocouple was placed at the center and top of the particle. A typical steady-state curve7 displayed the hysteresis behavior (Figure 2). Self-oscillations appeared on the higher branch of this curve at sufficiently large liquid flow rates. This phenomenon seems so unbelievable that it was decided to apply NMR imaging in further research and to take cylindrical pellets with larger height (H ) 10 mm). Figure 3 presents the results of the temperature measurement at the catalyst pellet bottom in the
process of self-oscillations and also two-dimensional images, which were obtained by NMR microimaging and characterized the filling of the catalyst pellet with the liquid (dark zones). The numbers at the temperature curve indicate the moments of time at which the microimages were made. It is seen that liquid AMS imbibes into the porous structure from the wetted pellet top and leads to the formation of a wetted zone in the catalyst pellet with clear-cut contours of the liquid distribution within the porous structure. The consecutive magnetic resonance microimages in Figure 3 show that the size and boundary of this zone continuously pulsate. The rest of the pellet down to the liquid front remains dry and filled with the gas phase (white zones). Note a time shift between the maximal temperature and the minimal volume of the wetted zone. When the movie was made from the set of microimages, it became clear that oscillation of the liquid front resembles the motion of a liquid piston. Mathematical Model By now it seems extremely difficult just to formulate a three-dimensional model of the process. Therefore, a simplified one-dimensional model was developed that meanwhile includes (as we believe) an intrinsic mechanism of self-oscillations. The basic simplification is that liquid forms a cylindrical zone of constant radius but variable height inside the pellet (Figure 4). This zone is surrounded by the dry porous layer of constant thickness ∆rdry. The bottom of the liquid-filled zone tends to move down owing to liquid imbibition (capillary forces) from the top. However, this tendency is counteracted by the process of evaporation. A part of the liquid evaporates without hydrogenation from the wetted top of the pellet. Another part evaporates in a narrow front within the pellet at the boundary of the liquid-filled zone, and the forming AMS vapor reacts with hydrogen in the dry porous layer surrounding this zone. The external surface of the pellet exchanges heat with the hydrogen flow. The heat transfer within the pellet occurs by heat conduction in the axial (vertical) direction.
9714
Ind. Eng. Chem. Res., Vol. 44, No. 25, 2005
The evaporation rate at the lateral surface is determined by the internal-diffusion mass transfer from the liquid-phase boundary through the dry layer of thickness ∆rdry. The internal-diffusion mass-transfer coefficient is calculated as
DAH
(τ)∆r
βpel )
, DAH ) 6 × 10-5
dry
Figure 4. Schematic of physical assumptions that are adopted in the mathematical model of the self-oscillation liquid front inside a catalyst pellet.
The model is written as a set of two differential equations:
0 e h e H, t g 0 2 liq π(R2Fcatccat p + Rliq Fliqcp )
Rliq2λliq)
πRliq2HFliq
2
∂T + 2πR[Wint ev (QrXcon - Qev) 2 ∂h Rgas(T - Tgas)] (1)
dφ ) Gliq - mliq(πRliq2Wext ev + dt 2πR
∫0H
liq(φ)
Wint ev dh) (2)
h ) 0 (pellet top)
Sh ) 2 + 0.6xReSc1/3, Nu ) 2 + 0.6xRePr1/3
Rgas )
∂T (3) ∂h
(4)
h ) H (pellet bottom) ∂T )0 ∂h
(5)
Initial conditions: t ) 0, T(h,0) ) T0, φ(0) ) φ0. At the top (h ) 0, 0 < r < Rliq) and lateral (r ) Rliq, 0 < h < Hliq) boundaries of the liquid-filled cylindrical zone, the conditions for equilibrium saturation of hydrogen with AMS vapor at the temperature T(h,t) are met. The AMS saturation vapor pressure and the specific heat of evaporation are found from the empirical correlations11
Psat ) 1.3367 × 1010
( )
Nuλgas Tgas , λgas ) 0.216 deq 373
-Hev , RgasT
1.9
(9)
(10)
The thickness of the lateral dry catalyst layer was considered as a single adjustable parameter of the model, and its value was chosen to fit the experimental data. It is important that this value appears realistic.
∆rdry ) R - Rliq ) 0.00045 m, ∆rdry/R ) 0.18 (11)
(6) Qev ) RgasTcTb ln Pc (Tc - Tb)1.38
(12)
The evaporation rate from the internal liquid front (across the dry porous layer) to the hydrogen flow is controlled by internal and/or external mass transfer. This rate is easily derived if we compare the internal molar vapor flux (from the liquid front to the pellet surface) to the external one (from the pellet surface to the hydrogen flow):
(
)
Ps Psat - Ps ) βpel Wint ev ) βgas RgasT RgasT
(13)
and, thus, the expression for the lateral surface partial pressure of vapor is derived from eq 13:
( )
βpel Ps Psat ) RgasT βpel + βgas RgasT
(14)
When eq 14 is substituted into eq 13, the final expression for the evaporation rate at the lateral surface is obtained:
Wint ev ) (Tc - T)0.38
(8)
0.78
φ ) Hliq/H
where
Rgas(T - Tgas) + λcat
( ) ( )
ShDAH Tgas , DAH ) 6 × 10-5 deq 373
πR2Rgas(T - Tgas) + πRliq2Wext ev Qev +
Psat Wext ev ) βgas RgasT(0,t)
(7)
The degree of the catalyst pellet liquid filling φ can be expressed through the ratio between the height of the cylindrical liquid-phase region and the pellet itself:
The boundary conditions for eq 1 are as follows:
2 2 cliq p Gliq(T - Tliq) ) (πR λcat + πRliq λliq)
1.9
Because the coefficient DAH of molecular diffusion of AMS vapor with hydrogen is close to that for cumene with hydrogen, the effect of multicomponent diffusion can be ignored and the DAH value can be taken to be equal to the coefficient of binary diffusion of vapor (mixture of AMS and cumene) with hydrogen. This allows one to calculate the interphase heat- and masstransfer coefficients using the conventional empirical correlations12
βgas )
∂T ) π(R2λcat + ∂t
T (373 )
(
)( )
βgasβpel Psat , if 0 < h < Hliq (15) βgas + βpel RgasT
For the dry (in our case bottom) part of the pellet, the evaporation rate is zero.
Ind. Eng. Chem. Res., Vol. 44, No. 25, 2005 9715
Wint ev ) 0, if Hliq < h < H
(16)
To estimate the degree of conversion, we assume that the rate of AMS hydrogenation in the liquid phase can be neglected in comparison with the rate of AMS hydrogenation in the gas phase; the AMS hydrogenation in the gas phase takes place in the dry surface layer of the pellet; the change in the hydrogen partial pressure along the thickness of the reaction catalyst layer is insignificant and is equal to the equilibrium pressure:
PH2 ) P - Psat(T)
(17)
and the hydrogenation kinetics is described by the power law1
Wr ) k0PH20.64 exp(-E/RgasT)
(18)
where the order in AMS vapor is zero. Under these assumptions, from known local values of AMS intrapellet evaporation rate (15) and (16) and hydrogenation rate (18), one can find the thickness ∆rreac of the dry layer (surrounding the liquid-filled cylindrical zone inside the catalyst pellet) that is necessary for 100% conversion of AMS vapor:
∆rreac ) Wint ev /Wr
(19)
The degree of AMS conversion is calculated from expressions (17)-(19) as follows:
if ∆rreac < ∆rdry, then Xcon ) 1
(20)
if ∆rreac > ∆rdry, then Xcon ) ∆rdry/∆rreac
(21)
Results of Modeling Numerical analysis of the mathematical model showed that, in the physically realizable parameter range, there is a unique steady-state solution. After a certain boundary in the space of the model parameters is crossed, this state can become unstable simultaneously with the “soft” birth of a self-oscillation mode (so-called Hopf bifurcation). Here we illustrate this by the influence of the AMS liquid flow rate. As has been mentioned already in the Experimental Section, the steady-state stability at low AMS flow rates was found experimentally. The model successfully confirms this observation. The modeling results at Gliq ) 0.7 mg/s are given in Figure 5. The phase trajectories are presented in Figure 5a. Figure 5b presents the dynamics of attainment of a stable steady-state mode of the catalyst pellet. In the initial moment of time (at t ) 0), the pellet is absolutely dry (φ ) 0) and its temperature is equal to the temperature of gas flow of 87 °C. Oscillations emerging after the start-up of the liquid feed are gradually damped, and a stable steady-state mode is established. The calculated temperature of the wetted zone (Figure 5b) is more than 35 °C higher than the temperature of the dry zone and the gas flow. The calculated degree of the catalyst pellet liquid filling φ is 0.28. The resulting stable steady-state temperature profile along the pellet height is shown in Figure 5c. The maximal temperature is achieved at the middle of the liquid-filled zone. AMS conversion for this and all other calculations was complete (Xcon ) 1). Taking into account that the specific heat of reaction exceeds the specific heat of evaporation
Figure 5. Calculated (a) phase trajectory, (b) dynamics of attainment of a stable steady state, and (c) steady-state temperature profile at Gliq ) 0.7 mg/s.
(Qr > Qev), one may conclude that a positive heat source arising at the lateral pellet surface (in eq 1) is distributed along the whole liquid-filled zone (0 < h < Hliq). The conductive heat flux to the left-hand side of the maximum (Figure 5c) is mainly consumed by the evaporation at the wetted top of the pellet (at h ) 0). The conductive heat flux from the maximum to the right-hand side transfers heat to the dry zone (Hliq < h < H), which is cooled by hydrogen flow. When the liquid flow rate has been increased from 0.7 to 0.8 mg/s with all other parameters remaining constant, then the new steady state became unstable. Figure 6a illustrates the development of self-oscillations after the start-up of the liquid feed and attainment of a periodic mode (cycle) at Gliq ) 0.8 mg/s. It is seen that the trajectory in the phase space winds, approaching an attracting orbit. Figure 6b compares the calculated dynamics for the upper (h ) 2 mm) and lower (h ) 6
9716
Ind. Eng. Chem. Res., Vol. 44, No. 25, 2005
feedback. In our case, this feedback acts through a change in the degree of liquid filling φ and the pellet temperature T. For example, if because of a random perturbation the heat release exceeds the heat removal, the catalyst pellet temperature begins to increase. This causes an increase in the partial pressure of vapor within the pellet and, hence, in the total evaporation rate, which leads to a decrease in the liquid filling and the evaporation surface area. A consequence of these processes is the deceleration with a time lag in the increase in the total evaporation rate, which leads to a decrease in the partial pressure of vapor and, consequently, in the heat release rate. If this time lag, which is caused by the delay in change of φ, is sufficiently large, then the system may “jump” over its steady state and, furthermore, start to swing. In our case, the delay in negative feedback is ensured by the fact that the characteristic time it takes for the pellet to be filled with the liquid is long enough in comparison with the thermal relaxation time of the pellet. The increase in the amplitude of such oscillations can be held back only by nonlinear mechanisms or additional restrictions (0 < φ < 1). Conclusions The self-oscillations of the temperature and liquid fraction within the partially wetted catalyst pellet in AMS hydrogenation were experimentally detected and explained using a simplified one-dimensional model. The results of the calculations agreed with the experimental data on the dynamics of the temperature in the bottom part of the cylindrical pellet and the dynamics of the liquid distribution within the pellet. This confirmed the previously proposed mechanism7 of critical phenomena and hysteresis of temperature and liquid filling for a partially wetted pellet during exothermic reactions accompanied by evaporation of the liquid reactant. In the proposed mechanism, a decisive contribution to the total conversion rate is made by the gasphase catalytic reaction, which is controlled by internaldiffusion evaporation. Note also the experimentally determined high nonuniformity of the temperature and liquid-phase distributions over a single catalyst pellet. The challenging future task is to find out how single oscillating catalytic particles may cooperate and whether such self-oscillation phenomena should be avoided in an industrial trickle-bed reactor or it may give rate enhancement. Figure 6. (a) Calculated phase trajectory approaching a closed orbit, (b) the comparison of the calculated and experimentally measured dynamics of the catalyst pellet temperature, and (c) oscillating temperature profiles at Gliq ) 0.8 mg/s.
mm) parts of the pellet with the experimental thermocouple measurement. Note the satisfactory agreement between the calculation results and the experimental data on the period (2-3 min) and the amplitude (85-93 °C) of the temperature oscillations at the bottom part of the pellet. Figure 6b also demonstrates a phase lag between the oscillations of the temperature and degree of the catalyst pellet liquid filling φ. The selfoscillations of a temperature profile along the pellet height at the discrete moments of time are shown in Figure 6c. It is known that self-oscillations are possible in the dynamic system in which there is delayed negative
Acknowledgment The financial support from NWO Grant 047.014.004 is greatly acknowledged. The authors are also thankful to Dr. I. V. Koptyug and Dr. A. Lysova for kindly presenting NMR microimages of the intrapellet liquid distribution. Notations ccat p ) 1070 J/(kg K), specific heat capacity of the material of the pellet (Al2O3) cliq p ) 2000 J/(kg K), specific heat capacity of liquid AMS DAH ) coefficient of binary diffusion of AMS vapor with hydrogen, m2/s deq ) 0.005 m, equivalent pellet diameter E ) 3.78 × 104 J/mol, activation energy of hydrogenation of AMS vapor on a Pd/γ-Al2O3 catalyst
Ind. Eng. Chem. Res., Vol. 44, No. 25, 2005 9717 Gliq ) (0.7-0.95) × 10-6 kg/s, liquid AMS mass flow rate h ) axial coordinate for a cylindrical pellet, m H ) 0.01 m, catalyst pellet height Hliq ) height of the cylindrical liquid-phase region within the pellet, m mliq ) 0.118 kg/mol, molar weight of AMS P ) 105 Pa, total pressure PH2 ) hydrogen partial pressure in the reaction zone, Pa Psat ) AMS saturation vapor pressure, Pa Ps ) vapor partial pressure on the pellet external surface, Pa Pc ) 33.6 bar, AMS critical pressure Qev(T) ) specific heat of AMS evaporation, J/mol Qr ) 1.16 × 105 J/mol, specific heat of reaction (hydrogenation of AMS vapor) at 110 °C R ) 0.0025 m, catalyst pellet radius Rgas ) 8.31 J/(mol K), ideal gas constant Rliq ) 2.05 × 10-3 m, radius of the cylindrical liquid-phase region inside the pellet, m T(h,t) ) temperature of the catalyst pellet, °C or K Tb ) 438.5 K, AMS boiling temperature Tc ) 654 K, AMS critical temperature Tgas ) 87 °C, inlet gas (H2) flow temperature Tliq ) 87 °C, inlet liquid (AMS) flow temperature t ) time, s Wr ) rate of hydrogenation of AMS vapor per unit catalyst pellet volume, mol/(m3 s) int Wev ) specific rate of AMS evaporation from the lateral surface of the cylindrical catalyst pellet, mol/(m2 s) Wext ev ) specific rate of AMS evaporation from the wetted (top) surface of the cylindrical catalyst pellet, mol/(m2 s) Xcon ) degree of conversion of AMS vapor Greek Symbols Rgas ) 132 W/(m2 K), coefficient of heat transfer between the pellet surface and the gas flow βgas ) 0.042 m/s, coefficient of mass transfer between the pellet surface and the gas flow βpel ) internal-diffusion mass-transfer coefficient, m/s ∆rdry ) 0.45 × 10-3 m, thickness of the dry subsurface catalyst pellet layer surrounding the lateral surface of the liquid-filled cylindrical zone ∆rreac ) thickness of the dry subsurface catalyst pellet layer surrounding the liquid-filled cylindrical zone that is necessary for complete conversion of AMS vapor, m φ(t) ) degree of the catalyst pellet liquid filling ) 0.5, coefficient of intrapellet porosity λcat ) 0.2 W/(m K), thermal conductivity coefficient of the material of the pellet (γ-Al2O3) λliq ) 0.2 W/(m K), thermal conductivity of liquid AMS Fcat ) 800 kg/m3, pellet density Fliq ) 900 kg/m3, liquid AMS density τ ) 3, coefficient of pore tortuosity
Dimensionless Groups Nu ) 3.1, Nusselt number Pr ) 0.7, Prandtl number Re ) 4.4, Reynolds number Sc ) 2.8, Schmidt number Sh ) 3.8, Sherwood number Subscripts and Superscripts ev ) evaporation ext ) external int ) internal pel, cat ) catalyst pellet liq ) liquid s ) steady-state value
Literature Cited (1) Germain, A. H.; Levebvre, A. G.; L’Homme, G. A. Experimental Study of Catalytic Trickle Bed Reactor. Adv. Chem. Ser. 1974, 133, 164-178. (2) Watson, P. C.; Harold, M. P. Dynamic Effects of Vaporization with Exothermic Reaction in a Porous Catalytic Pellet. AIChE J. 1993, 39, 989-1006. (3) Watson, P. C.; Harold, M. P. Rate Enhancement and Multiplicity in a Partially Wetted and Filled Pellet: Experimental Study. AIChE J. 1994, 40, 97-111. (4) Kulikov, A. V.; Kuzin, N. A.; Shigarov, A. B.; Kirillov, V. A.; Kronberg, A. E.; Westerterp, K. R. Experimental Study of Vaporization Effect on Steady State and Dynamic Behavior of Catalytic Pellets. Catal. Today 2001, 66, 255-262. (5) Kirillov, V. A.; Kulikov, A. V.; Kuzin, N. A.; Shigarov, A. B. Experimental Study of Heat Regimes on a Dry, Partially or Completely Wetted and Liquid Filled Catalyst Particle. Chem. Eng. J. 2004, 98, 219-235. (6) Harold, M. P. Bimolecular Exothermic Reaction with Vaporization in the Half-wetted Slab Catalyst. Chem. Eng. Sci. 1993, 48, 981-1004. (7) Shigarov, A. B.; Kulikov, A. V.; Kuzin, N. A.; Kirillov, V. A. Modeling of Critical Phenomena for Liquid/Vapor-Gas Exothermic Reaction on a Single Catalyst Pellet. Chem. Eng. J. 2003, 91, 205213. (8) Koptyug, I. V.; Sagdeev, R. Z. Non-Traditional Applications of NMR Tomography. Russ. Chem. Rev. 2003, 72, 165-191. (9) Koptyug, I. V.; Kulikov, A. V.; Lysova, A. A.; Parmon, V. N.; Sagdeev, R. Z. NMR Imaging of the Distribution of the Liquid Phase in a Catalyst Pellet During R-Methylstyrene Evaporation. J. Am. Chem. Soc. 2002, 124, 9684-9685. (10) Koptyug, I. V.; Lysova, A. A.; Kulikov, A. V.; Kirillov, V. A.; Parmon, V. N.; Sagdeev, R. Z. Functional Imaging and NMR Spectroscopy of an Operating Gas-Liquid-Solid Catalytic Reactor. Appl. Catal. A 2004, 267, 143-148. (11) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. Properties of gases and liquids, 4th ed.; McGraw-Hill: New York, 1987. (12) Hugmark, G. A. Mass and Heat Transfer from Rigid Spheres. AIChE J. 1967, 3, 1219-1227.
Received for review February 28, 2005 Revised manuscript received April 14, 2005 Accepted May 2, 2005 IE050269X
