Ind. Eng. Chem. Res. 2010, 49, 10641–10650
10641
Catalyst Deactivation in 3D CFD Resolved Particle Simulations of Propane Dehydrogenation Mohsen Behnam,† 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, PO Box 1, Belasis AVenue, Billingham, CleVeland TS23 1LB, U.K.
Catalyst deactivation by carbon deposition has been investigated for the dehydrogenation of propane to propene on a Cr2O3/Al2O3 catalyst. Computational fluid dynamics was used to couple the 3D transport and reaction processes occurring inside the cylindrical pellet to the gas flow around the pellet. The pellet scale reaction and carbon laydown are shown to be strongly affected by the bed scale tube wall heat flux supplied for the endothermic reactions, and the species distributions on the pellet surface are also affected by the ease of reactant access to the particle. The development of particle internal gradients and carbon accumulation are illustrated for the early stages of deactivation. Carbon deposition is initially strongest in the high temperature regions close to the tube wall. As time progresses, the increased deactivation caused by the carbon acts to reduce all rates of reaction, and propene production and coke formation shift to other regions of the pellet. 1. Introduction The deactivation of catalysts is an important problem in many processes, such as steam reforming and the catalytic dehydrogenation of alkanes, which feature highly endothermic reactions. The local carbon laydown can cause particle breakage or strongly decrease reaction rates, both of which remove the heat sink and lead to local hot spots and reactor tube failure. This is especially a problem for reactors with a low tube-to-particle diameter ratio (N), in the range 3 e N e 10, where a large fraction of the catalyst particles are located next to the heated tube wall. Conventional packed-bed reactor modeling involves several simplifying assumptions, such as the disregard of actual bed structure in pseudo continuum approaches, the utilization of correlation-based effective transport parameters which lump disparate transport phenomena, and reduction of the true threedimensional velocity field to a one-dimensional (sometimes position-dependent) velocity component. When the rate of reaction is inhibited by diffusion limitations inside the catalyst pellet, the source terms in the tube-scale balances use effectiveness factors based on single particle models with uniform catalyst pellet surroundings. These simplifications have been motivated by a desire to reduce the dimensionality and complexity of the governing equations, and by the difficulty of addressing the (semi-) random complex structure of the packed tubes. The problem with such simplifying assumptions is that local phenomena, resulting from subtle changes such as catalyst pellet design, are averaged out. The challenge we address here is to develop a better understanding of the local interactions among flow patterns, species pellet diffusion, thermal gradients, and catalytic reaction, with special attention to catalyst deactivation by coke deposition. The analysis of fluid flow, heat and mass transfer processes, and coupled chemical conversion in fixed bed reactors is an area that has received a great deal of attention from reaction engineers. For many decades, engineers have striven to obtain * To whom correspondence should be addressed. E-mail: agdixon@ wpi.edu. † Worcester Polytechnic Institute. ‡ Johnson Matthey.
better insight into reactor systems. As one approach to that goal, computational fluid dynamics (CFD) has been adopted as a suitable tool by many researchers,1,2 and the application of CFD to packed tube modeling has recently been reviewed.3 Although the use of CFD to simulate geometrically complex flows in randomly packed tubes 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 reactor tubes. CFD can be used to obtain detailed three-dimensional velocity, species, and temperature fields that are needed to improve engineering approaches. Frequently, CFD can be used to investigate situations that are not amenable to experimental measurement. In the present paper, the focus is on the development and use of CFD to obtain better understanding of the effects of local flow, heat, and mass transfer on rates of catalyst deactivation in cylindrical catalyst pellets. The inclusion of heterogeneous chemical reaction inside catalyst particles in CFD is a capability that is still an active area of research. Our group and others began by attempting to model the catalyst particles as porous regions with the velocity components constrained to zero.4,5 In subsequent work,6 we demonstrated that the porous particle approach had a serious deficiency, a convective flux across the particle-fluid interface which was an artifact of the numerical method, and which overestimated interphase transport. To address this problem, we developed an improved method that employed user-defined scalars to mimic the species mass fractions inside the catalyst pellets. The pellets could then be defined as true solid regions with “no-slip” boundary conditions. This method did not suffer from the porous method convective flow deficiency. Illustrative computations of methane steam reforming in spherical catalyst pellets were reported. We then presented7 extended applications of our solid particle method to investigate local interactions of particle processes and external flow for cylindrical particles. We illustrated the complexity of local distributions of species, temperature, and reaction rates for methane steam reforming (MSR) and propane dehydrogenation (PDH). It was shown that the local variations in species and reaction rate depended on the tube-wall-induced gradients in temperature in the near-wall particles, on the flow
10.1021/ie100456k 2010 American Chemical Society Published on Web 08/19/2010
10642
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
patterns around the particle exterior, and on the accessibility of fresh reactants to the particle surface. These results also changed according to total mass flow rate, temperature level, and diffusivity, among other variables. In none of the cases studied did we consider catalyst deactivation. The direct dehydrogenation of alkanes to alkenes has been developed and studied since the early 1930s.8 In particular, the dehydrogenation of propane to propene is becoming a commercial reaction of growing importance, due to the world demand for propene monomer for polypropene manufacture. Traditionally, propene has been mainly produced as an ethylene byproduct in steam cracking (∼ 75% in 1998), fluid catalytic cracking units (∼ 24%), and only 1% by on-purpose PDH manufacture. This picture has changed over the last ten years due to the more rapidly developing market for polypropene, compared to the market for ethylene derivatives, with stricter environmental regulations. This has increased interest in catalytic processes for propane dehydrogenation to propene, as an economically attractive alternative.9 A major problem is the fast deactivation of the catalyst by coke formation, spurring interest in different concepts to enable periodic regeneration of the catalyst by oxidation. Current processes employ either alumina-supported promoted chromia (Lummus/Houdry CATOFIN process,10 adiabatic fixed beds; Snamprogetti/Yarsintez process,11 interconnected fluidized beds) or alumina-supported platinum catalysts (UOP Oleflex process,12 adiabatic moving beds; BASF/Linde process,13 parallel tubular fixed beds); the Pt catalysts may be Sn-promoted to reduce coking (Phillips STAR process, parallel tubular fixed beds).14 A recent comparison of the various processes is available.15 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 pressure close to atmospheric. Some novel concepts to manage both the deactivation and the thermal problems include the use of a two-zone fluidized bed,16 membrane reactors,17 a rotating monolith reactor,18 and reverse flow reactors.19,20 To date, the economic and practical problems raised by these more sophisticated designs have prevented their commercial implementation. 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. The objective of this work is to use CFD to examine the interaction between flow and reaction/deactivation for the case of propane dehydrogenation (PDH) to propene, for which suitable reaction and deactivation kinetics are available. Detailed, fundamental studies of butane dehydrogenation deactivation kinetics over alumina-supported Cr2O3 are available,21,22 and similar mechanistic models have been developed for propane dehydrogenation over Pt-Al2O3.23 For propane dehydrogenation over aluminasupported Cr2O3 we have chosen to use the simplified kinetic and deactivation model of Jackson and Stitt24,25 to illustrate our methodology for incorporating unsteady catalyst deactivation into CFD simulations. This model has the virtues of simplicity and agreement with available experimental data. It is however noted that, while Jackson and Stitt do present mechanistic studies to support their kinetic model, it may be considered deficient26 in that it considers neither propene as a potential coke precursor nor the suppression effect of hydrogen on coke formation and, further, it neglects time effects on the structure of the coke and the effect this has on the deactivation impact of the surface deposits. It may also be noted, maybe more pertinently, that while Jackson and Stitt use the model to interpret reactor axial
Figure 1. (a) The 120° wall segment model with cylindrical packing and test particle, showing contours of propane dehydrogenation rate on fresh catalyst (kmol/m3(solid) · s); (b) details of mesh with thin prism cells on tube wall and both inside (thicker) and outside (thinner) of the pellet.
and temporal profiles, this model is in fact derived from global data obtained from integral reactor studies and its use to interpret local information should be treated with caution.27 2. Simulation Model Development CFD simulations of coupled flow, heat transfer, diffusion, and reaction were carried out in a 3D 120° periodic wall segment of an N ) 4 tube, as shown in Figure 1a. Focus was on a central full cylinder test particle at a 45° angle to flow, in realistic surroundings of particle segments. The simulations were not intended to represent the overall behavior of a reactor tube, for which a larger assembly of particles would be necessary. The cylindrical particles were placed by hand to reproduce the most common arrangement seen in a series of experimental packings in clear columns. 2.1. Geometry. A 120° wall segment (WS) modeling approach was used3,4 to investigate changes in performance for a single particle at the tube wall. The wall segment is illustrated in Figure 1a for a packing of equilateral full cylinder pellets. The bottom and top surfaces were identical, so that the model represented a geometry that varied in a repeating manner along the direction of flow. This gives rise to a situation for which nonreacting streamwise periodic flow conditions exist, in which the flow pattern repeats over a length L, with a constant pressure drop across each repeating section along the streamwise direction. This allowed the computation of a streamwise periodic flow from which a developed velocity profile could be obtained at the inlet of the wall segment, which was subsequently used for the nonperiodic reacting case simulations instead of taking a uniform inlet velocity, which would have been less realistic. The side walls of the WS model had symmetry conditions imposed on them. For cylindrical particles, this symmetry does not hold, and results near the sides of the model will not be realistic. An earlier study was carried out28 to address this point, which compared a WS model with cylinders to a full bed (CW) model built to extend the 120° segment. It was shown that flow and heat transfer for the middle 60° of the segment, which contained the test particle, were unaffected by the artificial symmetry conditions imposed at the side walls. An extension of this comparison was made for the present work, for the reacting case with flow, heat transfer, and diffusion. Several quantities were compared for the test particle, as listed in Table 1: temperature, mass fractions, heat uptake, reaction rates, and surface mass flow. The values for the CW model and the WS model were in very good agreement, confirming that the simulations for the test particle were not adversely affected by the presence of the symmetry planes in
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
the WS model. The WS model was therefore used for the main study, as the focus was on the behavior of the test particle at the wall. The near-wall center test particle 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, and were not intended to represent the performance of a realistic bed of cylinders. The WS model height was 0.0508 m, the tube radius was 0.0508 m, and the equilateral cylinder particles were 0.0254 m. 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. 2.2. Governing Equations. The conventional governing equations are the equations of conservation of total mass, momentum, and energy. These are well-described in any standard reference for CFD, and have been given in our previous work.6,7 We shall describe here only the solid particle method that we introduced recently, as this connects with the newer developments here on deactivation. In our solid particle method we define Nsp - 1 user-defined scalars φk, in both the fluid and solid phases, 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)
Nsp
1 ) Dek
∑ ∆1
kl
l)1
(
yl - y k Nsp
1 - yk
∑ l)1
N _l N _k
10643
)
(5)
N _l N _k
where 1 1 ) e + ∆kl Dkl
1 DeKkDeKl
(6)
ym
∑D
e Km
m
This approximate method assumes pressure variation inside the catalyst pellet is small compared to the fluid pressure, without making the stricter assumption of constant particle pressure. The diffusion coefficients were assumed to be unaffected by the deposition of coke. Equations 1, 3, and 4 are solved for the user-defined scalars in the fluid phase, then we set Yk ) φk 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. A user-defined code was created to describe the sinks/source terms in the catalyst particles. In the code, the species source/ sinks terms were defined as:
for k ) 1, 2, ... Nsp - 1, assuming isotropic diffusion. In the solid phase the kth scalar satisfies
NR
Si ) Fs
∑ (R r )M ij j
(7)
i
j)1
-∇ · (Γ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. In our work it represents production of species k due to reaction. In both phases the equations above are supplemented by
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
Q ) Fs
Nsp-1
φNsp ) 1 -
∑φ
{
j
(8)
j
j)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 )
∑ r (-∆H )
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 derived by Hite and Jackson29 from the dusty gas model, and is given by
where ∆Hj is the heat of reaction of reaction j, in J/kmol. The fluxes of 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. Care was taken to properly specify boundary conditions there to ensure equality of fluxes. 2.3. Meshing. To develop a suitable mesh for cylindrical particles, a study was made7 of the number and thickness of boundary prism layers on a single full cylinder in a box, at an angle of 45° to flow. The study was carried out for heat transfer, with heat sinks enabled in the test particle. The results showed that the most important factor in computing particle heat transfer
Table 1. Model Verification and Mesh Refinement quantity
WS mesh 030 (coarse)
WS mesh 025 (medium)
WS mesh 020 (fine)
CW
cell size (m) cell count (M) -∆P/L (Pa/m) Tav (K) Ypropene Yhydrogen Q (W) r1 (kmol/s) r2 (kmol/s) r3 (kmol/s) npropene (kmol/s)
0.0762 2.26 7.91 855.3 0.4138 0.01341 6.5441 4.8770 × 10-8 0.4281 × 10-8 0.1188 × 10-8 4.8864 × 10-8
0.0635 3.78 7.93 855.6 0.4142 0.01347 6.5258 4.8966 × 10-8 0.4298 × 10-8 0.1194 × 10-8 4.9065 × 10-8
0.0508 6.78 7.56 855.3 0.4134 0.01339 6.5091 4.9366 × 10-8 0.4285 × 10-8 0.1188 × 10-8 4.9512 × 10-8
0.0762 6.71 6.70 854.6 0.4121 0.01325 6.4848 4.8925 × 10-8 0.4240 × 10-8 0.1174 × 10-8 4.9194 × 10-8
10644
Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010
C3H8 f 3C + 4H2 -∆H2 ) -103.9 kJ/mol
Table 2. Reactor Conditions and Properties for PDH phase fluid solid
Tin [K] 873.15 873.15
qwall P F cp k [kW/m2] [kPa] [kg/m3] [J/kg · K] [W/m · K] 6.0
101
0.5847 1947
3700.9 1000
0.1210 1.0000
µ [Pa · s] 1.97 × 10
C3H8 f C + 2CH4 -∆H3 ) 45.9 kJ/mol -5
was a sufficiently thin first layer, and that the size of the UNS mesh was of secondary importance. The mesh consisted of boundary layer prism cells at both the inside and outside particle surfaces and at the tube wall, with tetrahedral cells in the main fluid volume. Three prism layers were implemented on the particle outer surfaces and four on the tube wall, of first layer thickness 2.54 × 10-5 m (0.001 in.) and increase ratio 1.2, to create a fine enough grid structure in the wall-fluid contact regions, to allow the laminar sublayer to be resolved. Additional prism layers were introduced on the interior particle surfaces to obtain a finer grid structure to capture any steeper temperature and species gradients near the particle surface. The mesh structure for a typical particle near a tube wall is shown in Figure 1b. Four prism layers were applied inside the particles with the first layer height of 7.6 × 10-5 m (0.003 in.) and an increase ratio of 1.2. The regions outside the prism layers were meshed by unstructured tetrahedral cells of 7.6 × 10-4 m (0.03 in.) for both fluid and particle volumes. Mesh-independence of the results was examined for the present work in the WS model using the same boundary-layer structure and three sizes of tetrahedral cell in the unstructured mesh. Results are shown in Table 1 for simulations at base conditions, i.e., three reactions but with no deactivation, which corresponds to the highest reaction rates and is the most stringent test. The test simulations therefore included simultaneous flow, diffusion, heat transfer, and reaction. Three sizes of mesh were consideredscoarse, medium, and finesand the corresponding cell sizes and cell counts are given in Table 1. Since the focus of the study was particle 2, the metrics chosen for comparison were all for that particle. The surface species flow and volumetric production by reaction of the propene were in excellent agreement in each case, showing that steady-state had been reached. Differences between the results produced by the three meshes were small (
