Flow, Transport, and Reaction Interactions for Cylindrical Particles

9026

Ind. Eng. Chem. Res. 2010, 49, 9026–9037

Flow, Transport, and Reaction Interactions for Cylindrical Particles With Strongly Endothermic Reactions M. Ertan Taskin,†,§ Alexandre Troupel,† Anthony G. Dixon,*,† Michiel Nijemeisland,‡ and E. Hugh Stitt‡ Department of Chemical Engineering, Worcester Polytechnic Institute, 100 Institute Road, Worcester, Massachusetts 01609-2280, and Johnson Matthey, P.O. Box 1, Belasis AVenue, Billingham, CleVeland TS23 1LB, United Kingdom

Interactions between reaction rates, conduction, and diffusion inside catalyst particles can be complex, especially when influenced by nonuniform surface conditions produced by the flow field external to the particle, or by the highly directional temperature field near a heated tube wall. In this work, a three-dimensional, realistic flow field was coupled to species and energy balances in cylindrical catalyst particles using computational fluid dynamics (CFD). Two strongly endothermic reactions were studied: methane steam reforming and propane dehydrogenation. Detailed pellet surface and intraparticle temperature, species, and reaction rate distributions were obtained for a near-wall particle. Nonuniform and nonsymmetric surface and intraparticle variations were observed. These effects are primarily attributed to the steep temperature gradients at the tube wall, as well as depletion of the reactants in regions of low or stagnant flow where particles approach each other closely. 1. Introduction Packed bed multitubular reactors with low tube-to-particle diameter ratios (N) are frequently selected for strongly endothermic reactions such as methane steam reforming (MSR) and propane dehydrogenation (PDH). For low-N tubes, the presence of the tube wall causes significant changes in bed structure, flow patterns, transport rates, and the amount of catalyst per unit volume, compared to high-N tubes. In particular, the particles close to the wall may be expected to behave differently from those inside the bed. Such particles make up a significant fraction of the total, ranging from 69% of all particles at the wall for N ) 4 to 30% for N ) 10. In conventional packed bed reactor modeling approaches, the transport phenomena have been described on the basis of simplifying assumptions, such as pseudohomogeneity, effective transport parameters, and plug flow. For reactions with strong diffusion limitations, reaction source terms use effectiveness factors based on averaged particle reaction rates with uniform surface conditions. Due to these simplifying assumptions, the effects of strong local gradients and catalyst pellet geometry changes in the near-wall environment may be lost. Computational fluid dynamics (CFD) has been considered as a useful tool by many researchers,1,2 and the application of CFD to packed tube modeling is rapidly increasing.3 Although the use of CFD to simulate geometrically complex flows is too expensive and impractical currently for routine design and control of fixed-bed reactors, the real contribution of CFD in this area is to provide more fundamental understanding of the transport and reaction phenomena in such reactors. CFD can supply the detailed three-dimensional velocity, species, and temperature fields that are needed to improve engineering approaches. In the present paper, the focus is on using CFD to * To whom correspondence should be addressed. E-mail: agdixon@ wpi.edu. † Worcester Polytechnic Institute. ‡ Johnson Matthey. § Current address: Artificial Organs Laboratory, Department of Surgery, University of Maryland School of Medicine, MSTF-454 10 S. Pine St., Baltimore, Maryland 21201.

obtain a better understanding of the local interactions between flow patterns, heat transfer, species pellet diffusion, and reaction in cylindrical catalyst pellets in the near wall region, which could help to inform efforts to develop simplified models of such pellets.4 In CFD, the Navier-Stokes momentum equations and, optionally, the energy and species balances are solved numerically. Studies of packed tube performance have been made5 using the lattice Boltzman method (LBM) with reaction on solid surfaces, but these have been restricted to low Re values and isothermal conditions. Simulations using a finite volume CFD code have been done more recently,6 again with surface reactions only. The inclusion of heterogeneous chemical reaction inside catalyst particles has not yet become a routine capability. Our group has recently included the heat effects of an endothermic heterogeneous reaction, methane stream reforming (MSR), in spherical, cylindrical, and shaped catalyst particles.7–10 Our method allowed a more realistic, although still limited, evaluation of different particle shapes and structures for reactor performance. In that approach, species were not available inside the solid catalyst particles, and the very strong diffusion limitations in the catalyst were introduced in the CFD model by an approximation, in which a step change was applied in species partial pressures, so that the species were at the fluid composition in a thin region near the surface of a particle, and zero further inside the particle. The diffusion of the chemical species inside the catalyst pellets and its effect on the chemical reactions, however, could not be coupled to the 3D external flow field, so that this approach was limited to only the thermochemical aspects of the MSR process. To overcome these problems, and to allow for the presence of species gradients within the catalyst pellets, we and others attempted to model the catalyst particles as porous regions with the velocity components constrained to zero.11,12 In ensuing work,13 we have found that the porous particle approach had a serious deficiency, namely, a convective flux across the particlefluid interface which was an artifact of the numerical method, and which led to an overestimation of interphase transport. We put forward an alternative method in which user-defined scalars

10.1021/ie1003619  2010 American Chemical Society Published on Web 08/19/2010

Ind. Eng. Chem. Res., Vol. 49, No. 19, 2010

were used to mimic the species mass fractions inside the catalyst pellet regions of the computational domain, which were defined as true solid regions. This “solid particle” method did not suffer from the convective flow deficiency of the porous particle method. We provided illustrative computations of methane steam reforming in spherical catalyst pellets. In the present work, we have extended applications of our solid particle method to investigate local interactions of particle processes and external flow for more realistic cylindrical particles. We considered two industrially important strongly endothermic reactions in heated multitubular reactors and restricted attention to the detailed study of a single catalyst cylinder in plausible surroundings near the tube wall. The operating conditions chosen were MSR under heat-transferlimited conditions near the entrance of a typical tube, and under equilibrium-limited conditions at midtube, and propane dehydrogenation (PDH), which is less strongly diffusion-limited, to provide a comparison to MSR. The strongly endothermic methane steam reforming reactions are the basis of widely used processes for the manufacture of synthesis gas, which may be used as a source for hydrogen production for the Fischer-Tropsch processes, or for methanol or ammonia synthesis.14 A major design consideration has been the efficient supply of the high amounts of heat required to maintain the reaction.15 The effects of excessive temperatures on the tubes could become a serious problem and lead to early tube failure. A second consideration is the limitations imposed by equilibrium on conversion, and efforts have been made to overcome these by using membrane reactors.16 The modern improved design of catalyst pellets can mitigate both of these problems, by allowing a closer approach to equilibrium, lower methane slip, and improved heat transfer.14,17 The developing market for polypropylene has increased the interest in catalytic processes for propane dehydrogenation to propene. Current processes employ mainly fixed bed operation, over supported chromia or platinum catalysts. A major problem is the fast deactivation by coke formation, spurring interest in different concepts to enable regeneration of the catalyst.18 Further difficulties are the highly endothermic nature of the reaction, and the equilibrium limitations that require operation at high temperatures (800-950 K) and at pressures close to atmospheric. In both of these processes, close attention is paid during reactor design and engineering to the catalyst pellet shape, size, and properties. The development of hot spots and areas of deactivation in the catalysts are examples of local phenomena which require examination at a detailed level, taking into account all the transport and reaction phenomena present in the reactor. Our aim in the present work is to construct a detailed picture for a typical case of a cylindrical particle near the reactor tube wall, using CFD, for several industrially important scenarios. 2. Simulation Model Development The simulations performed in this work focused on the behavior of a typical full cylinder catalyst particle, at the wall of a heated reactor tube. This configuration allowed the detailed study of the transport and reaction in and around the particle. Local variations in temperature, species, and reaction rate could be obtained, at a high degree of resolution. The simulations were not intended to represent the overall behavior of a reactor tube, for which a larger assembly of particles would be necessary. 2.1. Geometry. A 120° wall segment (WS) modeling approach was developed7,8 to investigate changes in particle shape and reaction conditions for either a small cluster of

9027

spherical particles or a single cylindrical particle at the tube wall. The segment is illustrated in Figure 1a for the case of equilateral full cylinders. The bottom and top (not shown) surfaces were identical and had periodic flow conditions imposed on them. The side walls of the WS model are shown in black in Figure 1a and had symmetry conditions imposed on them. For the case of spheres, this reflects the true angular symmetry of the near-wall particles. For cylinders, this symmetry does not hold, and results near the sides of the wall segment will not be realistic. We carried out a study to address this point,19 which compared a WS model with cylinders to a full bed model built to extend the 120° segment. It was shown that flow and heat transfer for the middle 60° of the segment were unaffected by the artificial symmetry conditions imposed at the side walls. The near-wall center test particle was located tangentially at the center of the model and was rotated 45°. This configuration was chosen as the most common that we observed in a visualization study of many cylinder packings.8 It was the only particle completely in the wall segment. The other particles in the geometry were not entirely in the WS model, and they were placed accordingly to provide typical surroundings for the test particle. The WS model height was 0.0508 m, the tube radius was 0.0508 m, and the equilateral cylinder particles were of size 0.0254 m. This gave a nominal tube-to-particle diameter ratio of N ) 4. The particles were placed so as to approach each other and the tube wall closely but did not touch, so no reduction in size was necessary to avoid contact points, as is usual in packings of spheres.7,11,13 2.2. Governing Equations. The governing equations for the fluid phase are the conventional equations of conservation of mass, momentum, and energy. These are well-described in any standard reference for CFD. Here, we provide a brief description of the new solid particle method that we introduced in our previous work, where a more complete description may be found.13 In our solid particle method, we define Nsp - 1 userdefined scalars φk, in both the fluid and solid phases. These scalars are used to represent the species mass fractions. Under steady-state conditions, the kth scalar in the fluid phase satisfies the equation ∇ · (Fu_φk - Γk∇φk) ) 0

(1)

for k ) 1, 2, ..., Nsp - 1, assuming isotropic diffusion. In the solid phase, the kth scalar satisfies -∇ · (Γk∇φk) ) Sφk

(2)

again for k ) 1, 2, ..., Nsp - 1, and where the source term on the right-hand side is a generation term for scalar φk, and in our work represents production of species k due to the reaction. In both phases, the equations above are supplemented by Nsp-1

φNsp ) 1 -

∑φ

k

(3)

k)1

as the mass fractions sum to unity. The fluid-phase and solid-phase diffusivities for each species are given by Γk )

{

FDek (solid) FDk + µt /Sct (fluid)

(4)

The fluid value is the sum of molecular and turbulent diffusivities; the solid value is an effective Fickian diffusivity

9028


derived by Hite and Jackson20 from the dusty gas model and is given by

1 ) Dek

(

N _l 1 y - yk ∆kl l N _k N _ Nsp l 1 - yk∑l)1 N _k

Nsp ∑l)1

)

NR

Si ) Fs

∑ (R r )M ij j

i

(7)

j)1

(5)

where Fs is the pellet density, Ri,j represents the stoichiometric coefficient of component i in the reaction j, and Mi is the molecular weight of species i. The heat generation/consumption by the reactions was calculated as NR

where

Q ) Fs 1 1 ) e + ∆kl Dkl

j

j

(8)

j)1

1 DeKkDeKl∑m

∑ r (-∆H )

ym

(6)

e DKm

This approximate method assumes that pressure variation inside the catalyst pellet is small compared to the fluid pressure, without making the stricter assumption of constant particle pressure. Equations 1, 3, and 4 are solved for the user-defined scalars in the fluid phase, in place of the usual mass balances in terms of the mass fractions Yk. The fluid velocity field, however, is computed on each iteration using the mass fractions Yk to obtain the average molecular weight, and both fluid heat capacity and density depend on Yk. We therefore set Yk ) φk in the fluid phase at the start of each iteration to ensure that the velocity field is modified (via the average molecular weight) as the composition of the gas phase changes due to reaction. Equations 2, 3, and 4 are solved for the user-defined scalars in the solid phase. A user-defined code was created to describe the sink/source terms in the catalyst particles. In the code, the species source/ sinks terms were defined as

where ∆Hj is the heat of reaction of reaction j, in joules per kilomole. For the industrial flow rates of the present study, the flow was usually turbulent, and the Reynolds-averaged Navier-Stokes (RANS) approach was used, with the shear stress transport (SST) k-ω model for closure, which reduces to a traditional twolayer zonal model if the near-wall mesh is fine enough to resolve the laminar sublayer (y+ ≈ 1). The fluxes for the user-defined scalars (species mass fractions) across the particle-fluid conjugate boundary were computed via user-defined code using the standard expressions for diffusive flow across cell interfaces in a nonorthogonal mesh, taking care to properly specify boundary conditions there to ensure equality of fluxes. This approach can be applied for laminar flows or for turbulent flows with fine meshes giving y+ ≈ 1 only, as it does not take into account wall laws that are used for coarse meshes in turbulent flow. 2.3. Meshing. Prism layers were introduced inside each solid particle at the surface to obtain a finer grid structure to capture the steep temperature and species gradients near the particle surface which are known to occur, especially in MSR. The mesh

Figure 1. (a) The wall segment model with cylindrical packing and test particle. (b) Details of mesh with thin prism cells on solid surfaces. (c) Visual planes inside the test particle. (d) Flow patterns around the test particle.


Figure 2. (a) Pressure drop (Pa/m) and (b) y+ as functions of first layer height on cylindrical particles.

structure for a typical particle near a tube wall is shown in Figure 1b. Due to the steep gradients near the particle surface, the model utilized for the MSR simulations had a different prism structure for the test particle than the other particles. Approximately, the outer 10% of the particle diameter was meshed by 15 prism layers with a first layer height of 3.8 × 10-5 m (0.0015 in.) and an increase ratio of 1.2. For the other particles, four prism layers were applied with a first layer height of 7.6 × 10-5 m (0.003 in.) and an increase ratio of 1.2. The remaining regions in the particles were meshed by unstructured tetrahedral cells of 7.6 × 10-4 m (0.03 in.) in size. For the higher-catalyst-activity PDH reaction, the above-described different prism structure for the test particle was not applied, and the intraparticle grid structure was kept the same for all of the particles. In order to achieve an accurate representation of the flow in the near-wall region, prism layers were implemented in the fluid domain on the tube wall and particle surfaces to create a fine enough grid structure to allow the laminar sublayer to be resolved. In a test for mesh-independence, we have previously reported10 a study of the number and thickness of boundary prism layers on a single full cylinder in a box, at an angle of 45° to the flow. The study was carried out for heat transfer, with heat sinks enabled in the test particle, and all other conditions the same as for MSR. The results showed that the most important factor in computing particle heat transfer is a sufficiently thin first layer, to obtain y+ e 1, and that the size of the UNS mesh was of secondary importance. An extension of that mesh-independence study is given here, where we consider the pressure drop in the full wall segment. Periodic flow simulations were carried out under the inlet conditions of steam reforming. The first prism layer height was varied from 10-6 m to 10-4 m. The UNS grid was applied for the rest of the fluid volume with a size of 0.000762 m. The results confirmed those of the single-particle study in general, and sample results for the variation of the first-layer thickness are shown in Figure 2. In this study, although the first layer height was changed by a factor of 100, the change in pressure drop was only around 7% on the basis of the data points lying on the trend line. The pressure drop difference for the models with different numbers of prism layers, shown inside the dashed ovals in Figure 2a, was a maximum of 3%. Note that the y+ values for those models were very much lower than the

9029

recommended value of 1.0. Since the prism layer structures were regular, compared to UNS grids, the y+ values varied linearly with the first prism height, as shown in Figure 2b. From these results, it is clear that for the wall segment cluster of particles, the pressure drop change is sufficiently small if we refine the mesh on the solid surfaces sufficiently to maintain y+ ≈ 1, which was obtained with a first layer prism thickness of 2.54 × 10-5 m (0.001 in.). 2.4. MSR Operating Conditions and Settings. The operating conditions considered for the methane steam reforming simulations were based on the industrial conditions obtained from a Johnson Matthey detailed reformer model of a methanol plant steam reformer, at the tube inlet and midtube locations.7 A Reynolds number of approximately 9500 based on superficial velocity and the particle diameter of a sphere of equivalent volume to the cylindrical particle was utilized. The operating conditions are shown in Table 1, and the initial values of species mass fractions are given in Table 2. The solid particles had the properties of alumina, with density Fs ) 1947 kg/m3, specific heat cps ) 1000 J/kg · K, and effective thermal conductivity keff ) 1.0 W/m · K. Species transport in the pellets was calculated using eqs 5 and 6 from straight-pore Knudsen and molecular diffusion coefficients corrected using pellet porosity and tortuosity values of ε ) 0.44 and τ ) 3.54. Under MSR tube inlet conditions, these gave values in the range 5 × 10-7 m2/s (CO2) to 2.5 × 10-6 m2/s (H2), while under MSR midtube conditions the range increased to 9 × 10-7 m2/s (CO2) to 3.4 × 10-6 m2/s (H2). The source terms in the solid phase were based on the kinetic model of Hou and Hughes21 for methane steam reforming over a Ni/RAl2O catalyst: CH4 + H2O ) CO + 3H2 -∆H1 ) -206.1 kJ/mol

(r1)

CO + H2O ) CO2 + H2 -∆H2 ) 41.15 kJ/mol

(r2)

CH4 + 2H2O ) CO2 + 4H2 -∆H3 ) -165.0 kJ/mol

(r3) The reaction rates and parameter values are given in the original reference and were put on a basis of unit volume of the solid particle and combined with the stoichiometric coefficients to obtain source terms for the user-defined scalar equations, as in eqs 7 and 8. 2.5. PDH Operating Conditions and Settings. The PDH reactor conditions and the fluid properties are also given in Table 1. The inlet mole fractions were 0.90 for C3H8, 0.05 for C3H6, and 0.05 for H2. The pellet properties were taken to be the same as for the MSR particles. The diffusivities were calculated with the same procedure as the MSR calculations and gave higher values than for methane steam reforming due to the lower operating pressure: 2.8 × 10-6 m2/s (C3H6) to 1.4 × 10-5 m2/s (H2). Due to the importance of catalyst deactivation in alkane dehydrogenation, the PDH reaction kinetic expressions usually include coke formation reactions.22 In many experimental studies at relatively lower operating temperatures and shorter contact times, the propene production reaction [C3H8 ) C3H6 + H2] is the only one considered. In this work, for simplicity, we have focused only on that reaction, and we will defer consideration of the coke deposition reactions to a later time. The steadystate flow solution was obtained for 4030 h-1 gas hourly space velocity (GHSV), which corresponds to a Reynolds number of 143 based on superficial velocity and the particle diameter of a sphere of equivalent volume to the cylindrical particle. Kinetics

9030


Table 1. Reactor Conditions and Fluid Properties model

Tin [K]

qwall [kW/m2]

P [kPa]

F [kg/m3]

cp [J/kg.K]

kf [W/m.K]

µ [Pa.s]

MSR - tube inlet MSR - midtube PDH

824.15 1019.15 873.15

113.3 117.3 6.0

2159 2110 101

6.1616 3.8615 0.5847

2395.38 2658.56 3700.9

0.0876 0.1685 0.1210

3.00 × 10-5 3.78 × 10-5 1.97 × 10-5

Table 2. Bulk Fluid Inlet Species Mass Fractions for MSR Reactions location

CH4

H2

CO

CO2

H2O

tube inlet midtube

0.1966 0.0926

0.0005 0.0442

0.0007 0.1181

0.1753 0.2771

0.6269 0.4680

for the reversible propene reaction were obtained from previous published work on chromia-catalyzed PDH by Jackson and Stitt:22

(

r ) k1 CCsHs -

CCsHsCHz Keq

)

(9)

and the rate and equilibrium constant values presented in that reference were used without alteration. The tube wall heat flux was adjusted by trial and error in the absence of any definite information, to give an essentially zero temperature rise in the fluid, as has been suggested for industrial processes such as the BASF/Linde PDH process in a fired multitubular reactor. 2.6. Computational Procedure. The governing equations described above were solved using the finite volume commercial CFD code Fluent, version 6.3. The pressure-based segregated solver was used, with the SIMPLE scheme for pressure-velocity coupling. All test results presented here used first-order upwind interpolation for the convection terms; all diffusion terms used second-order discretization. The second-order upwind interpolation for the convection terms frequently became unstable at the pressure outlet surface, due to the unavoidably large amount of recirculating flow there. Under-relaxation factors were usually left at the Fluent default settings, unless some instability was observed in the iterations, when they were occasionally reduced. The simulations were first run to determine an initial isothermal constant-composition flow solution in the segment with periodic top and bottom conditions. This flow field was used subsequently to provide a nonuniform inlet velocity profile for the reacting case in which the wall segment top and bottom were set as velocity inlet and pressure outlet boundaries, respectively, with the periodic restriction removed. For the reacting case, continuity, momentum, turbulence, energy, and species (user-defined scalar) equations were solved, giving 11 simultaneous PDEs for the MSR runs and nine simultaneous PDEs for the PDH runs. The flow simulations were not decoupled from the energy and species simulations in this work, as we wished to allow the changes in moles and temperature due to reaction to have their proper effect on the flow. The convergence was monitored by the pressure drop value for the flow runs, and by the energy balance and the reaction rates in the test particle for energy and species simulations, in addition to the residuals. Simulations were run on a Sun Microsystems X2200 M2 x64, a 64-bit server with two dualcore processors (4 CPU total) at 2.6 GHz each with 8 GB of RAM. Convergence took from 15 000 to 25 000 iterations, due to the slow diffusion processes occurring within the solid particles. Computation times were on the order of 24 h for 1000 iterations. 3. Results and Discussion Methane steam reforming simulations were run corresponding to two positions in the reactor tube. The role of heat transfer

under tube inlet conditions is crucial, as it must supply the needs of the endothermic reactions and at the same time increase the particle temperature up to the range for which reforming to CO and H2 (reaction r1) dominates over reforming to CO2 and H2 (reaction r3). In addition, there are strong diffusion limitations in the particles. Under midtube conditions, about 200 K higher, intrinsic reaction rates are much higher, but diffusivities have increased to a lesser extent, and the process has even stronger diffusion limitations. Reaction r1 dominates over reaction r3 at the midtube. Both of these conditions give rise to reaction rates that are significant only in the outer shell of the catalyst, so that surface distributions become of primary interest. For the propane dehydrogenation process, reaction rates are also high, but intraparticle diffusivity is much higher, so that reaction occupies a larger fraction of the pellet volume. For this case, distributions of species and energy through the particle volume take on higher importance. Considering the above, the results presented below are based on both intraparticle variations and surface variations. Intraparticle resolved gradients were examined on two internal pellet planes, shown in Figure 1c. Plane 1 cuts the pellet lengthwise and is roughly perpendicular to the tube wall. Therefore, the effect of wall heat flux and near-wall temperature gradients on the pellet performance can be shown using this plane. Plane 2 was perpendicular to plane 1, also cutting the particle lengthwise, but was approximately parallel to the tube wall and was used to show gradients that were much less affected by wall heat flux. Since the test particle side surface was exposed to large temperature gradients, we have focused on the temperature, species, and reaction rate variations on that surface. To examine the surface variations, the lateral surface of the test particle was unwrapped. To do this, we performed geometrical transformations on the surface coordinates which mapped the test particle surface to a standard cylinder centered on the origin and aligned with the z axis. Any point on the surface could then be identified using only two parameters in cylindrical coordinates: the angle θ and the height Z, the radius r remaining constant. Finally, a surface coordinate s was created as the product of the angle and the radius (s ) θr), and the unwrapped surface was plotted in s-Z coordinates. The flow pattern around the test cylinder is shown in Figure 1d and is similar to those presented previously.8,19 Flow accelerates in the region between the pellet and the tube wall. On the upper curved surface, there is a vortex type of flow feature, caused by the blocking effect of the particle located above the test particle. The test particle bottom flat surface causes significant lateral displacement of the flow, which may be seen toward the center of the picture. This displacement helps cause the region of recirculation directly above. These flow features can be related to the transport and reaction plots as discussed below. 3.1. MSR Inlet Conditions. Intraparticle Variations. The temperature, methane and hydrogen mass fraction, and rate of reaction r3 contours obtained are shown in Figure 3, on planes 1 and 2. Most of the particle was at a lower temperature than the surrounding fluid, due to the endothermic reforming reactions. Close to the tube wall region at r/rp ) -1.0, the particle


9031

Figure 4. Temperature and species mass fraction contours of the test particle unrolled side surface for MSR reaction tube inlet conditions.

One of the important reaction engineering parameters is the effectiveness factor, η, to express the averaged intraparticle behavior. It was calculated for the ith reaction by the following expression:

∫

Figure 3. Intraparticle temperature (K), species mass fraction, and rate of reaction r3 (kmol/m3 · s) contours on planes 1 and 2 for MSR reaction tube inlet conditions.

was at a higher temperature due to the wall heat flux. This caused a temperature gradient across the particle of up to 35 K. On plane 1, higher reactant methane depletion was observed corresponding to the higher temperature. Correspondingly, product hydrogen and carbon monoxide contents were high at the same region. Plane 2 was located perpendicular to the wall heat flux, and so there were no significant temperature and species variations across the plane. This demonstrates that for the near wall particle the section closest to the tube wall was highly affected by the wall heat flux. From the species mass fraction contours, a sudden change can be noticed at the edges of both planes due to the relatively fast reaction and slow diffusion, which illustrate the “egg shell” type of steam reforming pellet behavior. From previous work,7 it was established that reaction r3 was the dominant one for the tube inlet operating conditions, being approximately an order of magnitude greater than reaction r1 near the surface. Rate contours for reaction r3 are also provided in Figure 3. As can be seen, the reaction was confined to the outer region of the catalyst particle, approximately 5%, in agreement with the active region used in our heat effects studies.9,10 The thin reaction zone insulated the pellet center from external mass transfer, so that the species mass fractions attained values consistent with reaction equilibrium at the prevailing local temperature. Inside the pellet, there was essentially no methane consumption, and almost zero net reaction rate was seen for the entire inner particle.

1 r dV Vp Vp i ηi ) ri(Ts,aV, Yis,aV) The volume averaged rates for the test particle were calculated by summing volume times the reaction rate in the individual cells of the particle. The reaction rates at surface conditions were calculated using average surface values of temperature and species. The effectiveness factors for reactions r1 and r3 were found to be 0.178 and 0.086, respectively. These values indicate that the MSR reaction is essentially on the particle surface, confirming that catalyst activity is proportional to pellet geometric surface area.14 Surface Variations. Figure 4 presents the temperature and species mass fraction contours on the test particle side surface. As a consequence of the supplied heat from the tube wall, a hot spot was noticed on the particle surface near the closest approach of the pellet to the tube wall (s ) 0.04 m). Due to the convective flow, which accelerates in the region between the pellet and the tube wall, the hot spot is displaced toward the right (s ) 0.05 m). Relatively cooler regions were observed for s > 0.06 m and s < 0.03 m. Furthermore, the lowest temperature was seen at around s ) 0.07 m. The methane and hydrogen mass fraction contours were very similar in shape; low YCH4 and high YH2 roughly corresponded to each other, but not to high T. There was considerable local variation on the surface, even away from the heated wall, up to 24 K for temperature and somewhat less for the mass fraction in methane (0.190-0.195). If these variations can be ascribed to the influence of flow on external transport, then this difference makes sense as the external heat transfer resistance is known to be significant compared to intraparticle conduction resistance,

9032


Figure 6. (a) Detailed contours of the section in between particle 2 (test particle) and particle 8. (b) Special viewing plane cutting through regions of low methane.

Figure 5. Reaction rate (kmol/m3 · s) contours of the test particle unrolled side surface for MSR reaction tube inlet conditions.

unlike external mass transfer, which is usually minor compared to intraparticle diffusion barriers. These observations can be coupled to the reaction rate contours as shown in Figure 5. All three reactions had similar contours and showed maxima and minima in the same locations. Reaction rates on the surface varied by factors of 2 (reaction r3) to 7 (reaction r1), following the greater sensitivity to the temperature of reaction r1 due to its higher activation energy. The rate of the water-gas shift (WGS) reaction r2 was quite low as compared to the other two reactions. Under the tube inlet conditions, the WGS reaction proceeded in the forward direction to produce H2. In general, on the pellet surface, the reaction rates followed the temperature variation. Regions of highest temperature were also where reaction rates were highest. The regions of lowest temperature, however, did not necessarily correspond to the lowest reaction rates. The regions of lowest CH4/highest H2 also did not correspond to regions of highest temperature/reaction rate. These observations suggest that mass transfer/diffusion plays a role in the surface distribution of species and influences the location of lower surface reaction rates. Methane Depletion Regions. As discussed above and shown in Figure 4, surface species mass fractions varied away from the near-wall region. Figure 6a highlights some zones of methane depletion between the test particle (part 2) and one of its neighbors (part 8). As can be seen from the figure, these particles were positioned close to each other. In order to study these depletion regions in more detail, a cross-sectional plane was created which passed through some of these zones, as shown in Figure 6b. The first parameter considered on that cross-sectional plane was the methane mass fraction. In Figure 7, two zones of relatively low methane concentration could be observed, marked by the two dark ellipses. A higher hydrogen mass fraction (not shown) also corresponded to the two zones of methane depletion. A careful look showed that the methane depletion zone appeared in a narrow zone between particles, where the

flow velocity was low. The velocity magnitude is shown in that same cross-sectional plane, marked by two white ellipses, and low velocity corresponded well with the depleted methane. Reaction rate 3 was also low in those locations, indicated by the darker regions in the surface cells. The temperature, however, was less than 5-10 K lower in the two regions than elsewhere, which confirmed that the low methane regions were not due to thermal effects. As we saw previously in Figure 4, the temperature and mass fraction extremes were not located exactly in the same place. At first glance, this seems contradictory, since a high temperature increases the reaction rate, thus more methane should be consumed and more hydrogen produced. So, it is natural to expect the maximum temperature, the minimum methane mass fraction, and the maximum hydrogen mass fraction to be located at the same place. To examine this, a new cross sectional plane was created, passing through the temperature maximum and the methane mass fraction minimum. The plots for temperature, velocity magnitude, methane mass fraction, and reaction rate 3 are given in Figure 8. It may be noticed that the pellet maximum temperature was not achieved at the section of that surface which was the closest to the heated tube wall, but in fact slightly further away. This may be explained by the velocity profile. There was a strong flow in the narrow gap between the pellet and tube wall, cooling down that section of the pellet. On the contrary, the hottest section was near a low flow zone. Consequently, heat exchanged by convection was low, and there was less cooling of the particle surface. However, the region where the velocity was lowest was not near that hot zone but was actually between the two particles; this was where the methane concentration was the lowest. Methane mass transport to the particle surface by turbulent flow was greatly reduced, so the renewal of reactants in the nearby section of the pellet was low. The resulting reaction rate was lower and occupied a thinner region near the particle surface between the two particles. In addition to the effects of velocity, the confinement of fluid between the two particles reduced the availability of fresh reactant. Usually, methane can diffuse to the surface (dark arrows) from five different directions, as illustrated by the cartoon in Figure 9a; for a zone between the two pellets, it could only diffuse from four directions as shown in Figure 9b. Moreover, the pellets were zones of methane consumption, and therefore, the methane was drained from two surfaces (light arrows, case b) as opposed to only one surface for a free particle (light arrows, case a). In other words, in a confined fluid region, less methane diffused in and more diffused out than for a free fluid next to a surface. As a result, the surface methane concentration for a confined region was lower than for the rest of the fluid.


9033

Figure 7. Contours of temperature, velocity, methane mass fraction, and reaction r3 on the plane defined in Figure 6.

Figure 8. Close-up of contours in the vicinity of low methane regionsstemperature, velocity, methane mass fraction, and reaction r3.

Figure 9. Schematic of geometric restrictions on methane supply.

3.2. MSR Midtube Conditions. Intraparticle Variations. The contours of T, YCH4, and YH2 for planes 1 and 2 in Figure 10 under midtube conditions are qualitatively similar to those under

tube inlet conditions, but closer examination reveals some interesting differences. The temperature penetration into the pellet was less intrusive, as shown on plane 1, which could be attributed to the increased thermal demand of the endothermic reforming reactions at higher temperatures. A small amount of heat and a slightly increased temperature may be seen at the top flat surface of the particle, as a result of the unimpeded flow conditions there. A temperature difference of about 28 K existed between the front and back of the particle. The methane mass fraction showed a more developed profile inside the pellet than for inlet conditions, and the gradient from the unheated surfaces toward the minimum YCH4 at the tube wall surface can be clearly seen. YH2 was more uniform within the particle, probably due to its higher diffusivity. Plane 2 showed no

9034


Figure 11. Temperature (K) and species mass fraction contours of test particle unrolled side surface for MSR reaction midtube conditions.

Figure 10. Intraparticle temperature (K), species mass fraction, and rate of reaction r1 (kmol/m3 · s) contours in planes 1 and 2 for MSR reaction midtube conditions.

intraparticle gradients parallel to the tube wall, aside from the usual sharp changes at the particle surfaces. Reforming to CO and H2 is known to be the dominant reaction under these conditions. This time, the reactions were found to be confined to the outer 1-2% of the particle radius, in agreement with literature reaction engineering results under similar midtube conditions.15 At the midtube temperatures, r1 increased 50-fold compared to the colder tube inlet, whereas r3 remained at comparable values, due to decreased partial pressure of reactant CH4 and increased partial pressures of products CO, CO2, and H2, which competitively adsorb on the catalyst. The effective diffusivities for the pellet increased only on the order of 50% from the tube inlet values, so that the reactions were even more strongly diffusion limited. The effectiveness factors of reactions r1 and r3 were calculated to be 0.0092 and 0.0096, respectively, which shows the lower utilization of the test particle as compared to the tube inlet simulation results. These values were in agreement with the industrial observations,14 and with the pseudocontinuum modeling results,15 where the authors focused on the axial middle location of the reactor. In our study, we have utilized the 3D flow field around the explicitly positioned catalyst particles and considered the local nonsymmetric variations. Surface Variations. Figure 11 represents the temperature and species surface contours for MSR midtube operating conditions. A high temperature region (at s ) 0.05 m) was observed, qualitatively similar to the one for the inlet operating conditions. Temperatures were high in the region 0.03 m < s < 0.06 m,

whereas for the other regions temperatures closer to the bulk values were observed. Additionally, a slightly higher temperature region was noticed at s ) 0.07 m and Z ) 0. The species mass fraction contours corresponded to the temperature field so that higher T, lower YCH4, and higher YH2 coincided. In these plots, patches of depleted methane corresponded to areas of increased temperature, as expected. Overall, the distributions of temperature and species over the surfaces appeared to be more uniform, which may have reflected the closer approach to equilibrium of the reaction system at the midtube conditions. The patch of higher T, lower YCH4, and higher YH2 at s ) 0.06-0.07 did reflect the lower external heat and mass transfer rates in the recirculation flow region. The three reaction rates are shown in Figure 12, where r1 and r2 were considerably increased over the values shown earlier in Figure 5. The two reforming reactions followed temperature and were highest in the hottest regions, and all three reaction rates showed similar contour patterns. Additionally, the WGS reaction rate contour showed that, at the hot spot, the reverse reaction proceeded, producing higher levels of carbon monoxide. In general, methane consumption rate and WGS reaction rate contours were opposite each other; for example, it was noted that, at the back of the particle (s ) 0 and s ) 0.08), the rates of reactions r1 and r3 were lower than the front values, whereas the WGS rate was higher at the back. The observed surface reaction rate distributions may be explained by considering that the reaction rates were higher than the diffusion rates for all locations, and therefore temperature and species variations were in quite good agreement with the reaction rate variations for the midtube operating conditions. 3.3. PDH Reaction. Intraparticle Variations. The temperature, species, and reaction rate contours are presented in Figure 13 for planes 1 and 2. For both planes, higher temperature regions were observed for the dimensionless radial region -1 < r/rp < -0.2, corresponding to proximity to the tube wall in the case of plane 1. The propane mass fraction YC3H8 was quite high on the particle surfaces, except on the surface next to the tube wall, where the reactant was depleted. The hydrogen mass fraction YH2 was high near the tube wall, and showed a decrease


9035

Figure 12. Reaction rate (kmol/m3 · s) contours of test particle unrolled side surface for MSR reaction midtube conditions.

Figure 14. Temperature (K), species mass fraction, and reaction rate (kmol/ m3 · s) contours of test particle unrolled side surface for PDH reaction operating conditions.

Figure 13. Intraparticle temperature (K), species mass fractions, and reaction rate (kmol/m3 · s) contours in planes 1 and 2 for PDH reaction operating conditions.

toward the back of the particle on plane 1, and decreased from particle center to edge on plane 2. The species mass fraction

and reaction rate contours showed that the reaction took place further inside the particle and was not as confined to the outer region of the pellet as observed in the MSR reaction simulations. On plane 2, the temperature and propane were higher, and the hydrogen lower, on the left side of the pellet. This is attributable to good convective transport, which increased the heat and mass transfer coefficients there, compared to the other surfaces. This combination resulted in a higher reaction rate at the dimensionless radial location -1 < r/rp < -0.8 of plane 2. The effectiveness factor for the test particle was found to be 0.247 for the PDH reaction operating conditions, which was also a supportive finding for the larger particle activity observation discussed above, and in quite good agreement with the literature.22 Surface Variations. The detailed side surface temperature, species mass fraction, and PDH reaction rate contours are presented in Figure 14. Generally the region of 0.02 m < s < 0.06 m was hotter than the other parts. The highest temperature was at s ) 0.04 m, where the particle approached the tube wall. There was no displacement due to flow, in contrast to the MSR case. Also unlike the previous MSR operating conditions, hotter regions were observed on both top and bottom of the side surface at s ) 0.04 m. Those top and bottom edges were the closest parts of the catalyst to the tube wall. In this region, the lower flow rate provided less convective heat transfer. The wall heat flux to the test particle was not reduced by a high rate of heat transfer to the flow along the pellet side surface, so the hot spots were due to the proximity of the edges to the wall, instead of being found on the midsections of the side surface, as they were for MSR. Similarly to the MSR case, relatively cooler regions were observed for s > 0.06 m and s < 0.03 m, and the lowest temperature was seen at around s ) 0.07 m. The species mass fraction side surface contours showed strong variations with position. The propane (C3H8) and propene (C3H6)

9036


mass fractions were opposite each other; i.e., the propene product mass fractions were low at the locations where the propane reactant mass fractions were high. For this reason, the propene contours are not shown here. At 0.02 m < s < 0.04 m, the propane content was high for the entire height of the pellet, whereas at 0.04 m < s < 0.06 m, there were sections with a low propane mass fraction. At s ) 0.06 m and in the upper half of the particle (Z > 0 m), the reactant propane mass fraction and product hydrogen mass fraction took their lowest and highest values, respectively, but this location did not correspond to the highest temperature. This was similar behavior to that observed for MSR, suggesting that species transport again plays a role. The reaction rate took its highest values at the top and bottom for s ) 0.04 m, corresponding well to the temperature maxima. In this case, the reaction rate minimum was at s ) 0.06, whereas the temperature minimum was closer to s ) 0.07 m. The reaction rate followed the temperature distribution more closely for PDH than for MSR but again appeared to be affected by the reactant distribution. The final distribution of the surface reaction rate for PDH was therefore determined by a combination of temperature and propane mass transfer. 4. Conclusions Three-dimensional flow fields have been coupled to diffusion/ reaction in a nonisothermal full-cylinder catalyst particle for different reactions with different activity levels. Three dimensional temperature, species, and reaction rate fields have been obtained within the particle without assuming symmetry or uniform surface conditions. Contours of these quantities have been presented on internal and external surfaces to illustrate that there is significant local variation present in a reacting catalyst particle near a heated wall. The primary cause of variations within the catalyst particle is, of course, the strong heat flux from the tube wall. This causes strong temperature, species, and reaction rate variations in the direction perpendicular to the wall, so that nonsymmetric scalar fields and nonuniform surface distributions are found, which is quite contrary to the conventional assumption in catalyst pellet models. It was also seen that the flow field has a strong influence on the diffusion/reaction limitations, being associated with regions of reactant depletion on the particle surface. From maps of the test particle surfaces, species, and temperature, extremes do not always have the same location; high reaction rates and temperature maxima matched for MSR and PDH but low reaction rates and temperature minima were displaced. Furthermore, it has been noticed for MSR inlet conditions that the temperature maximum and reactant minimum are always located near a low velocity region. This is especially true for temperature, since the area which gets the hottest is not the one the closest to the wall but, rather, located near the wall and near a low velocity region. Therefore, for MSR, both temperature and reactant concentration minima can be explained by their proximity to a low velocity zone, where convective transport for both energy and species are low. Nevertheless, an extra parameter comes into play; otherwise the reactant minimum concentration would still be located at the same place as the temperature maximum. It appears that the fluid near the reactant minimum is also located in a confined region between two particles. Hence, species diffusion can only come from a limited number of directions. Moreover, since this confined zone is between two particles, the reactant is depleted by diffusion toward two different particle surfaces. Compared to the fluid in more open regions, less can diffuse in, but more has to diffuse out. The reactant concentra-

tion minimum can therefore be explained by its close location to a low flow region and a confined zone. When the diffusion limitations are strong, as under the MSR operating conditions, the temperature and species variations are affected so that the reactions are confined to an outer region of 1% to 5% of particle radius, depending on the operating conditions. For less significant diffusion limitations as in the PDH reaction, the intraparticle variations were shifted toward the inside of the particle as a result of the higher particle activity. The methodology developed in our group has allowed full 3D simulations of conduction, diffusion, and reaction in a catalyst pellet to be coupled to 3D external flow, heat, and mass transfer in the fluid. The variations of species have not been presented before, and the previous approximate analysis of heat effects of reacting particles is considerably extended by the present work. A detailed picture of the local variations of temperature, species, and reaction rates has emerged, as well as their connection to local flow features. The reaction rates for both the MSR and PDH processes have been represented by simple overall kinetics; the methodology presented here could equally well be applied to microkinetics models, at the cost of considerable computing expense. The qualitative 3D behavior of the particle has been illustrated here; future work will assess the quantitative significance of the nonsymmetric and nonuniform temperature, species mass fractions, and consequently reaction rates.

Acknowledgment Acknowledgement is made to the Donors of the American Chemical Society Petroleum Research Fund for support of this research.

Appendix Nomenclature cp ) fluid specific heat, J/kg · K cps ) solid specific heat, J/kg · K dp ) particle diameter, m k ) turbulent kinetic energy, J/kg kf ) fluid thermal conductivity, w/m · K keff ) effective particle thermal conductivity, w/m · K Lp ) particle height, m Mi ) molecular weight of species i N ) tube-to-particle diameter ratio Nsp ) number of species P ) static pressure, kPa Q ) heat generation/consumption term qwall ) wall heat flux, w/m2 r ) particle radial coordinate, m rt ) tube radius, m ri ) reaction rate i, kmol/m3 · s Re ) Reynolds number, Fdpvz/µ s ) arc length, m Si ) species source/sinks term T ) temperature, K Vz ) axial velocity, m/s |V| ) velocity magnitude, m/s y+ ) dimensionless wall coordinate yk ) mole fraction of species k

Ind. Eng. Chem. Res., Vol. 49, No. 19, 2010 Yk ) mass fraction of species k Z ) particle height coordinate, m Greek Letters Ri,j ) stoichiometric coefficient of i in reaction j ∆Hj ) heat of reaction j, kJ/mol ε ) pellet porosity Γ ) mass diffusion coefficient, kg/m · s φk ) kth user-defined scalar µ ) fluid viscosity, kg/m · s F ) fluid density, kg/m3 Fs ) solid density, kg/m3 τ ) pellet tortuosity ω ) specific dissipation rate, s-1

Literature Cited (1) Kuipers, J. A. M.; van Swaaij, W. P. M. Computational Fluid Dynamics Applied to Chemical Reaction Engineering. AdV. Chem. Eng. 1998, 24, 227. (2) Ranade, V. Computational Flow Modeling for Chemical Reactor Engineering; Academic Press: New York, 2002. (3) Dixon, A. G.; Nijemeisland, M.; Stitt, E. H. Packed Tubular Reactor Modelling and Catalyst Design Using Computational Fluid Dynamics. AdV. Chem. Eng. 2006, 31, 307. (4) Mariani, N. J.; Keegan, S. D.; Martinez, O. M.; Barreto, G. F. A One-Dimensional Equivalent Model to Evaluate Overall Reaction Rates in Catalytic Pellets. Trans. Inst. Chem. Eng. Part A 2003, 81, 1033. (5) Freund, H.; Zeiser, T.; Huber, F.; Klemm, E.; Brenner, G.; Durst, F.; Emig, G. Numerical Simulations of Single Phase Reacting Flows in Randomly Packed Fixed-Bed Reactors and Experimental Validation. Chem. Eng. Sci. 2003, 58, 903. (6) Kuroki, M.; Ookawara, S.; Ogawa, K. A High-Fidelity CFD Model of Methane Steam Reforming in a Packed Bed Reactor. J. Chem. Eng. Jpn. 2009, 42 (Supplement 1), s73. (7) Dixon, A. G.; Nijemeisland, M.; Stitt, E. H. CFD Simulation of Reaction and Heat Transfer near the Wall of a Fixed Bed. Int. J. Chem. Reactor Eng. 2003, 1, A22. (8) Nijemeisland, M.; Dixon, A. G.; Stitt, E. H. Catalyst Design by CFD for Heat Transfer and Reaction in Steam Reforming. Chem. Eng. Sci. 2004, 59, 5185. (9) Taskin, M. E.; Dixon, A. G.; Stitt, E. H.; Nijemeisland, M. Approximation of Reaction Heat Effects in Cylindrical Catalyst Particles with Internal Voids Using CFD. Int. J. Chem. Reactor Eng. 2007, 5, A41.

9037

(10) Taskin, M. E.; Dixon, A. G.; Nijemeisland, M.; Stitt, E. H. CFD Study of the Influence of Catalyst Particle Design on Steam Reforming Reaction Heat Effects in Narrow Packed Tubes. Ind. Eng. Chem. Res. 2008, 47, 5966. (11) Dixon, A. G.; Taskin, M. E.; Stitt, E. H.; Nijemeisland, M. 3D CFD Simulations of Steam Reforming with Resolved Intraparticle Reaction and Gradients. Chem. Eng. Sci., 2007, 62, 4963. (12) Kolaczkowski, S. T.; Chao, R.; Awdry, S.; Smith, A. Application of a CFD Code (Fluent) to Formulate Models of Catalytic Gas Phase Reactions in Porous Catalyst Pellets. Chem. Eng. Res. Des. 2007, 85, 1539. (13) Dixon, A. G.; Taskin, M. E.; Nijemeisland, M.; Stitt, E. H. A CFD Method to Couple 3D Transport and Reaction in Fixed Bed Catalyst Pellets to the External Flow Field. Ind. Eng. Chem. Res. Submitted (doi: 10.1021/ ie100298q). (14) Stitt, E. H. Reactor Technology for Syngas and Hydrogen. In Sustainable Strategies for the Upgrading of Natural Gas: Fundamentals, Challenges and Opportunities; Derouane, E., Parmon, V., Lemos, F., Ramoa-Ribiero, F., Eds.; NATO Science Series II, Kluwer Press: Dordrecht, The Netherlands, 2005; Vol 191, p 185. (15) Pedernera, M. N.; Pina, J.; Borio, D. O.; Bucala, V. Use of a Heterogeneous Two-dimensional Model to Improve the Primary Steam Reformer Performance. Chem. Eng. J. 2003, 94, 29. (16) Gallucci, F.; Paturzo, L.; Basile, A. A Simulation Study of the Steam Reforming of Methane in a Dense Tubular Membrane Reactor. Int. J. Hydrogen Energy. 2004, 29, 611. (17) Kagyrmanova, A. P.; Zolotarskii, I. A.; Vernikovskaya, N. V.; Smirnov, E. I.; Kuz’min, V. A.; Chumakova, N. A. Modeling of Steam Reforming of Natural Gas Using Catalysts With Grains of Complex Shapes. Theor. Found. Chem. Eng. 2006, 40, 155. (18) Stitt, E. H.; Jackson, S. D.; Shipley, D. G.; King, F. Modelling Propane Dehydrogenation in a Rotating Monolith Reactor. Catal. Today 2001, 69, 217. (19) Taskin, M. E.; Dixon, A. G.; Stitt, E. H. CFD Study of Fluid Flow and Heat Transfer in a Fixed Bed of Cylinders. Num. Heat Tr. 2007, 52, 203–218. (20) Hite, R.; Jackson, R. Pressure Gradients in Porous Catalyst Pellets in the Intermediate Diffusion Regime. Chem. Eng. Sci. 1977, 32, 703. (21) Hou, K.; Hughes, R. The Kinetics of Methane Steam Reforming Over a Ni/R-Al2O Catalyst. Chem. Eng. J. 2001, 82, 311. (22) Jackson, S. D.; Stitt, E. H. Propane Dehydrogenation over Chromia Catalysts: Microcatalysis, Deactivation, Macrokinetics, and Reactor Modeling. Curr. Top. Catal. 2002, 3, 245.

ReceiVed for reView February 16, 2010 ReVised manuscript receiVed June 30, 2010 Accepted July 19, 2010 IE1003619

Flow, Transport, and Reaction Interactions for Cylindrical Particles

Recommend Documents