Ind. Eng. Chem. Res. 2004, 43, 7049-7056
7049
CFD Flow and Heat Transfer in Nonregular Packings for Fixed Bed Equipment Design Alfredo Guardo,† Miguel Coussirat,‡ M. Angels Larrayoz,† Francesc Recasens,† and Eduard Egusquiza*,‡ Chemical Engineering Department and Fluid Mechanics Department, ETS d’Enginyeria Industrial de Barcelona, Universitat Polite` cnica de Catalunya, Av. Diagonal 647, ETSEIB, 08028 Barcelona, Spain
This work aims to test the application of computational fluid dynamics (CFD) modeling to fixed bed equipment design. Studies of CFD with a fixed bed design commonly use a regular packing approach to define bed geometry. However, assuming nonregular packing is a more realistic way to simulate the behavior of a fixed bed and therefore to estimate important design parameters. As a fluid flow simulation tool, CFD allows us to obtain a more accurate view of the fluid flow and heat transfer mechanisms present in fixed bed equipment. Forty-four spheres stacked in a nonregular maximum-space-occupying arrangement in a cylindrical container were used as the geometrical model. Estimates of the pressure drop along the bed, and wall heat transfer parameters were chosen as validation parameters. ∆P, Nuw, and kr/kf are given for different values of Re (transition and turbulent flow), and they are compared to commonly used correlations. Air was chosen as the flowing fluid. Cases of laminar and turbulent flows are presented, and their results are compared. To account for the fluid flow and thermally fluctuating components in the turbulent cases, one- and two-equation turbulence models were used for simulation. Introduction In recent years, there has been considerable interest in optimizing global efficiency in production processes and minimizing waste generation due to the general trend of market globalization, environmental actions, higher client expectations, and increased profit revenue for manufacturing companies. The traditional approach of taking a product from the laboratory to a pilot plant and then to production scale is no longer attractive due to the high costs involved. Product and process development is carried out almost simultaneously, and fast analysis and prototype design capability are required to meet expectations.1 Computational fluid dynamics (CFD) is one of the critical “enabling technologies” for achieving this. It allows process engineers to predict, manipulate, and design the desired fluid dynamics in process equipment. In the modeling and design of fixed bed equipment, CFD can be used to simulate singlephase and multiphase flow through porous media and to perform detailed modeling of a packed bed, and it is used to design equipment with single-phase flow through a porous medium.2 The models assumed for flow patterns in a fixed bed are one of the main limitations when modeling a fixed bed reactor. It is known that the global behavior of a fluid in a transport system depends directly on the local flow structure.3 However, in almost all cases, hydrodynamic models for fixed beds apply plug flow restrictions. In recent developments, the plug flow model has been improved by adding a radial distribution to the axial velocity component.4 The classical approach to modeling * To whom correspondence should be addressed. Tel.: +34934016714. Fax: +34-934015812. E-mail:
[email protected]. † Chemical Engineering Department. ‡ Fluid Mechanics Department.
the rate of heat transport in fixed bed equipment has been to use an effective thermal conductivity in the radialsand, if necessary, axialsdirection together with an apparent heat transfer coefficient for the wall.5 These parameters are obtained from experimental data fitted to data regression analysis of actual models, usually assuming plug flow. However, there are many differences between data published by different authors.6 Thermal energy in a fixed bed is transported by strong convective radial flows while the fluid moves through the packing elements. Experimental studies have been developed for determining flow maldistribution in packed columns, considering factors such as flow redistributors design, bed height, flow rate, and physical properties of the fluid. They have found that the distribution of fluid inside the column is far from uniform, showing channeling of flow at the column wall.7 Dispersion and flow in porous media have been studied using magnetic resonance imaging (MRI), allowing stagnant flow regions, radial flows, and back-flow regions within the fixed bed to be identified.8,9 This confirms the presence of a strong radial flow at the fluid inlet into the bed, developing preferential flow zones where heat transfer could be accentuated, creating nonuniform temperature gradients along the bed. CFD applications to simulate fluid flow in porous media have a broad range of applications from oil basin simulation10 to modeling corn seed drying.11 In all cases, the results depend on an appropriate geometrical model, mesh definition, and the selection of a turbulence model (when considered necessary). There have been few CFD studies in fixed beds due to complex factors such as geometry definition and heat and mass transfer modeling, in addition to limitations in computational power. The first CFD approaches included 2D studies that resolved flow patterns and heat transfer around pro-
10.1021/ie034229+ CCC: $27.50 © 2004 American Chemical Society Published on Web 09/24/2004
7050
Ind. Eng. Chem. Res., Vol. 43, No. 22, 2004
posed ideal geometries.12,13 Other studies have used variational methods for evaluating steady-state velocity profiles for isothermal, incompressible flow in 2D packed beds, showing that there is channeling at the wall when the particle-to-channel diameter ratio is low.14 Although this study was developed for an ideal case, it has been the basis for the development of many other studies, contributing to improved theoretical predictions for the behavior of catalytic reactors.15 The unit-cell approach is another option for CFD modeling of fixed bed reactors. In this case, a periodic and repetitive structure within the packing is identified and subsequently used to define the computational domain limits. This geometrical strategy has been used for flow and heat transfer modeling in catalytic reactors,16,17 where it is used to predict heat transfer parameters for different values of Re. Yet these studies have been limited to regular packing schemes. The approach presented here is the simulation of fluid flow and heat transfer in nonregular packing. Laminar and turbulent flow solutions are obtained and compared. For turbulent solutions, one- and two-equation Reynolds-averaged Navier-Stokes (RANS) models are used. The results obtained are compared to state-of-the-art correlations for estimating pressure drop and heat transfer parameters in a fixed bed. Conduction/convection mechanisms are taken into account. A gas-solid system was chosen for this purpose.
Figure 1. 2D model obtained from a transversal cut through a unit cell.
Geometrical Model Geometrical modeling is one of the most critical stages in CFD simulation; correct definition of the geometry provides a more realistic scenario for the simulation, and the technique used for constructing the geometry will ensure the feasibility of generating a mesh good enough to capture all of the phenomena involved in the problem. The first step was to select a proper arrangement for the fixed bed. In a bed with a mixture of particle sizes if average particle size is used in the calculations, heat and mass dispersion follow the predictions for a bed of monosized particles.18 Therefore, a homogeneous sphere stack was selected for this study. Regular packing (i.e., a simple cubic or hexagonal lattice) does not offer good wall-to-particle contact unless the tube-to-particle diameter ratio is 2. Other diameter ratios lead to large empty spaces near the wall, which generate high-speed flow channels at the wall and, therefore, a bad flow distribution and incorrect calculation of parameters. Previous work using a finite-element (FEM) code has been done on this subject by Guardo et al.19 A maximumspace-occupying arrangement must be used to simulate a more realistic case. So, a 2D approach to the circle packing theory is used to build a unit-cell and then extrapolate it to 3D geometry. Circle packing is an arrangement of circles inside a given boundary such that no two of them overlap and some (or all) of them are mutually tangent. Solutions for the smallest diameter circles into which n unit circles can be packed have been proved optimal for 1 e n e 65.20 Furthermore, for certain values of n, several distinct optimal configurations are possible. In this case, following the instructions of Melissen21 for building the 2D model, an 11-circle arrangement with a diameter ratio equal to 3.923+ was chosen for the 3D extrapolation. Figure 1 shows the 2D model mentioned above.
Figure 2. Generated geometric model (lateral view).
Figure 3. Generated geometric model (isometric view).
An 11-sphere arrangement, with 9 particle-to-wall contact points and 14 particle-to-particle contact points, was built on the basis of the preceding 2D approach. A four-layer arrangement (44 spheres) with a 60° rotation around the reactor axis within each layer was chosen as the geometrical model for CFD simulations. Figures 2 and 3 show a lateral and isometric view of the constructed geometry, respectively. The construction of wall-to-particle and particle-toparticle contact points is also an important subject in model generation. Previous work reports no contact points between surfaces,12,13 or the emulation of contact
Ind. Eng. Chem. Res., Vol. 43, No. 22, 2004 7051 Table 1. Dimensionless Groups’ Orders of Magnitude 101 101 101-103 10-4-10-2 10-1 10-10-10-8 10-10-10-8
Re [dpuF/µ] St [L/tu] Eu [∆P/Fu2] Fr [u2/Lg] Pr [Cpµ/k] Br [µu2/kT] Ec [Br/Pr]
102 101 10-1-101 10-2-100 10-1 10-8-10-6 10-8-10-6
103 101 10-1 100 10-1 10-6 10-6
points (leaving small gaps between surfaces and assuming zero velocity in the gap) to avoid convergence problems.16,17,19 In this study, to include real contact points, the spheres were modeled overlapping by 1% of their diameters with the adjacent surfaces in the geometric model. Convergence problems were not detected during simulation runs. The modeled geometry was constructed following the bottom-up technique (generating surfaces and volumes from nodes and edges) to control the mesh size around critical points (i.e., particle-to-particle and particle-towall contact points). This was necessary to avoid grid element skewness and also to gain computational resources by reducing the number of elements in zones of low interest (i.e., geometrical zones away from contact points or away from the walls). Mesh Design and CFD Modeling To properly design a mesh capable of capturing the transport mechanisms present in the study in detail, a dimensionless analysis of Navier-Stokes equations under simulation conditions was developed. The dimensionless equations corresponding to mass, momentum, and energy balances are as follows:
()
u0 δF0 ∂Fˆ u0 bˆ ) bˆ ‚∇ B ‚u + δF0(u B Fˆ ) ) F0 (-Fˆ ∇ t0 ∂tˆ L L St
()
(1)
bˆ ∂u bˆ ‚∇ bˆ ) + Fˆ (u B )u ∂tˆ 1 ∇ B ‚{µˆ [∇ Bb u + (∇ Bb u)T]} + Re 1 1 (Fˆ ˆf ) (2) ∇ B ‚{µˆ T [∇ Bb u + (∇ Bb u)T]} + ReT Fr b
Eu(-∇ B pˆ ) +
St
()
∂T ˆ bˆ ‚∇ + (u B )T ˆ ) ∂tˆ 1 1 1 kˆ 1 1 kˆ T Ec [vˆ (τc:∇ Bb uˆ )] + ∆T ˆ + ∆T ˆ (3) Re Re Pr Fˆ Re PrT Fˆ
( )
( )
The orders of magnitude of the dimensionless groups were estimated by taking physical-chemical property values for air from experimental data and empirical correlations available in the literature.22,23 Reynolds number was calculated using particle diameter as characteristic length. Reynolds’ analogy was used to estimate values of PrT from ReT (see Table 1). The dimensionless analysis shows that turbulent forces make little contribution to eq 2 at low Re, but their contribution becomes more important as Re increases. Buoyancy forces and pressure drop are the most important terms in eq 2 at low Re, but their contribution decreases as Re increases. In eq 3, the diffusive term becomes important for energy balance. Steady-state analysis was considered for both equations.
Dimensionless analysis allows us to identify the problem as heat transfer in transition or turbulent flow. Accurate modeling of the flow mixing in this case implies using a highly refined homogeneous mesh, but this requires many computational resources. One of the aims of this work is to verify whether accurate modeling is possible if the homogeneous mesh requirement is relaxed, to minimize simulation times. Theoretically, the mesh should be able to appropriately define the boundary layer around the geometry present in the model for the laminar solution. In the case of a turbulent solution, mesh density will depend on the near-wall modeling strategy adopted for resolving the problem and will be determined by the characteristic y+ parameter.24 Transition from laminar to turbulent flow in fixed beds has not been extensively studied with numerical CFD simulation, and there are still doubts about when the turbulent model should be activated, because there are no reliable guidelines to predict the flow transition in complex geometries, such as fixed bed reactors or extraction equipment. Experimental studies have found that a transition from laminar to turbulent flow in a damped bed of spheres occurs over the range from 110 to 150 for the Reynolds number, and that around Re ) 300 the flow pattern is turbulent.25 Other authors have stated that a transition from laminar to turbulent flow occurs at Re ) 100.26 These results should be used as an indicator of when it is necessary to activate a turbulent model. In this work, laminar and turbulent flow solutions are calculated and shown to compare their performances. For the turbulent solutions, the turbulence intensity boundary condition at the bed inlet was set to 4%, as estimated from common pipe flow equations. To verify the choice of the turbulence intensity parameter, a simulation was run applying a lower value of this parameter, with no significant change in the results obtained. Although physically the variation of the turbulence intensity at the flow inlet can generate changes in flow behavior, numerically these changes are not detected by most of the turbulence models when the difference between the imposed value and the real value is reasonable.24 The boundary conditions for the model equations are as follows: (i) constant temperature and velocity of the fluid at the inlet, (ii) constant temperature at the wall, (iii) constant pressure at the fluid outlet, and (iv) nonslip conditions at the wall and particle surfaces. Model Analysis Navier-Stokes equations and energy balance were solved using commercially available finite volume code software Fluent 5.0. The fluid was taken to be incompressible, Newtonian, and in a laminar or turbulent flow regime due to the impossibility of RANS models to capture flow transition. Air was chosen as the simulation fluid, for which the constants were available in the software database. The incompressible ideal gas law (for density) and the power law (for viscosity) were applied to the model to make these variables temperaturedependent. Simulations were run on an HP C3000 workstation, and simulation times ranged from 12 to 48 h depending on the case studied. Second-order upwind schemes were selected to compute the field variables. The pressure-velocity coupling algorithm was the SIMPLE scheme. This scheme derives an equation for the pressure from the discrete continuity equation.27 Numerical convergence of the
7052
Ind. Eng. Chem. Res., Vol. 43, No. 22, 2004
Table 2. Velocity Conditions Applied to the CFD Model and G, µ, and Re Obtained for Each Case laminar u [m/s] 0.075 0.1125 0.150 0.225 0.375 0.525 0.75
F
[kg/m3] 1.0126 1.0666 1.0761 1.0688 1.0881
Spalart-Allmaras
µ [kg/m‚s]
Re
2.12 × 10-5 2.04 × 10-5 2.03 × 10-5 2.04 × 10-5 2.01 × 10-5
85 133 177 265 447
F
[kg/m3] 1.0253 1.0430 1.0551 1.0729 1.0944 1.1075 1.1192
standard k-
µ [kg/m‚s]
Re
2.10 × 10-5 2.08 × 10-5 2.06 × 10-5 2.04 × 10-5 2.00 × 10-5 1.99 × 10-5 1.97 × 10-5
86 130 175 265 448 633 912
F
[kg/m3] 1.0108 1.0276 1.0398 1.0580 1.0814 1.0960 1.1105
µ [kg/m‚s]
Re
2.12 × 10-5 2.10 × 10-5 2.08 × 10-5 2.06 × 10-5 2.02 × 10-5 2.00 × 10-5 1.98 × 10-5
84 127 170 259 438 619 893
model was checked on the basis of the residuals of all computed variables. For more complete convergence verification, the drag force over particle surfaces and the average static temperature at the bed outlet were also chosen as monitors. Results and Discussion The objective of this work is to test CFD as a design tool for fixed bed reactors. Simulations were run for several values of Re (see Table 2) with constant temperatures at the bed inlet and wall. Slight differences in the values of Re, F, and µ can be observed among laminar and turbulent simulations with an identical inlet velocity; the values of Re were first estimated at normal conditions and then, after the simulations, corrected with the obtained numerical average values of F and µ. The real problem is approached by means of nonregular packing modeling. Standard correlations for pressure drop and heat transfer parameters were selected as reference values to be compared against the numerical results generated. The use of “real” contact points was successfully included during model generation. The CFD model was solved for both laminar and turbulent flow situations, and buoyancy terms were activated for low Re simulations. For turbulence modeling, the Spalart-Allmaras and the standard k- turbulence models were used. Details of turbulence models can be found in the literature.24 The results of the simulations are discussed below. Velocity Profiles. To study the velocity distribution along the packed bed, a vertical cut was made along the cylinder’s diameter to generate velocity vector plots. Velocities at the inlet were set between 7.5 × 10-2 and 7.5 × 10-1 m/s. The inlet temperature was taken as 298.15 K, and the wall temperature was taken as 423.15 K. Air was used as the fluid, and density and viscosity were taken as temperature-dependent properties. Re was calculated using mass flow rate and average viscosity given by the software. Velocity profiles were also observed in the near-wall region of the modeled arrangements. As expected, in all of the cases analyzed, flow channeling took place near the wall and inside the bed, due to the presence of constrained flow areas. Strong radial flow from the middle to the wall was also noticeable. Due to the channeling of the flow (strong axial flow and reduced radial flow) at the wall, the local radial heat transfer rate decreases, causing the well-known “temperature jump” near the wall. Figures 4 and 5 illustrate the near-wall flow and the velocity profiles along a transversal cut in the fixed bed. It was also noticed that velocity increased by up to 4 times the inlet velocity in some constrained areas of the fixed bed. Stagnation points and secondary flows were noticed near the contact points (see Figure 6). As Re increases, within the range of velocities studied, eddy
Figure 4. Velocity vectors profile near the wall for Re ) 633. The velocity profile is expressed in m/s.
Figure 5. Velocity vectors profile along a transversal cut (z ) 0.0375 m) for Re ) 448. The velocity profile is expressed in m/s.
flow becomes easier to identify around the spheres, for both simulated turbulence models. All of the features mentioned above are also found by Suekane et al.,9 who used an MRI technique to directly measure the velocity of flow in a pore space that models a simply packed bed. Temperature Contours. Figure 7 shows temperature contour plots for Spalart-Allmaras model simulations to illustrate how the temperature field changed with flow rate. It shows a strong influence of the fluid flow pattern on temperature profiles. The temperature profiles penetrate the bed faster for low Re than for high Re flow patterns, because the mixing in the zones nearest the wall increases as Re increases. Analyzing the kinetic energy (k) profiles obtained with the k- model simulations, the aforementioned idea becomes clear (see Figure 8). A comparison between two selected cases shows increasing kinetic energy inside the packed bed with increasing Re. The temperature profiles show that the idea of thinking of one single value of Nuw for
Ind. Eng. Chem. Res., Vol. 43, No. 22, 2004 7053
Figure 6. Velocity vectors profile in a cross section (y ) 0) for Re ) 633. The velocity profile is expressed in m/s.
the whole fixed bed is not realistic; however, it can be useful for equipment design. Pressure Drop along the Bed. The results obtained using CFD simulation for the pressure drop along the bed were compared to Ergun’s correlation for pressure drop in beds packed with spherical particles.28 (Average density and viscosity of the fluid within the bed were taken for Ergun’s equation calculations.) See Figure 9 for details. All of the models used (laminar, Spalart-Allmaras, and standard k-) show good agreement with Ergun’s equation. This can be explained by the fact that the velocity fields obtained for laminar and turbulent solutions are similar. The pressure drop calculation is intrinsically related to the velocity field and is not
affected by mixing parameters or additional diffusive terms included within the different models’ equations. In the case of turbulence modeling, the results from the Spalart-Allmaras model show better agreement with Ergun’s prediction of the frictional pressure drop in the fixed bed than those from the standard k- model. The better near-wall treatment in the SpalartAllmaras model under the meshing conditions used for simulations favors velocity-pressure coupling and the estimation of drag coefficients over the involved surfaces. The situation is similar when estimating heat transfer parameters. A more complete discussion of the influence of the near-wall treatment on parameter estimates follows. Heat Transfer Parameter Estimates. To calculate the values of Nuw, kr, and Bi, kf was taken as a constant reference value (kf ) 0.0242 W/m‚K). The mass flow rate given by the software and the average values for the viscosity along the bed were used to calculate Re. Figures 10 and 11 show the values of Nuw and kr/kf, respectively, for different Re values, which was calculated according to previously published methods.29 To show the influence of the choice of flow model on the values of Nuw and kr/kf, values for the laminar and both the Spalart-Allmaras and the standard k- turbulent solutions are shown for Re between 84 and 912. The CFD results in Figure 10 for the SpalartAllmaras and standard k- turbulent solutions show good agreement with empirical models such as those proposed by Olbrich and Potter30 or Dixon and Cresswell5 for the Nuw estimate, as opposed to the simple model proposed by Li and Finlayson.31 These complex models take into account important factors such as pressure drop or the geometrical characteristics of the bed that are also included within the calculations of the CFD code for modeling the heat transfer phenomena. It is important to notice that in the transition zone, as Re decreases the values obtained for Nuw tend to be underestimated if compared to the values obtained with the empirical correlations shown in the case of a
Figure 7. Temperature contours in a cross section (y ) 0) for (a) Re ) 86 and (b) Re ) 912 applying the Spalart-Allmaras turbulence model. The temperature profile is expressed in K.
7054
Ind. Eng. Chem. Res., Vol. 43, No. 22, 2004
Figure 8. Kinetic energy contours in a cross section (y ) 0) for (a) Re ) 84 and (b) Re ) 893 applying the standard k- turbulence model. The kinetic energy profile is expressed in m2/s2.
Figure 9. Pressure drop versus Re: (×) laminar; (O) SpalartAllmaras; (+) standard k-; (---) Ergun equation.28
Figure 10. Nusselt number versus Re: (×) laminar; (O) SpalartAllmaras; (+) standard k-; (s - s) Dixon and Creswell (laminar);5 (- s -) Dixon and Creswell (turbulent);5 (s) Li and Finlayson;31 (---) Olbrich and Potter.30
turbulent solution. In the case of the laminar solution, it should be noticed that the first three values of the series (Re < 177) agree excellently with the values predicted by Dixon and Creswell5 for laminar flow. When Re ) 265 and greater, the laminar solution overestimates the value of Nuw.
Figure 11. Effective radial conductivity (kr/kf) versus Re: (×) laminar; (O) Spalart-Allmaras; (+) standard k-; (s) Yagi and Kunii;33 (---) Yagi and Wakao.34
The behavior of the solutions for the laminar and turbulent cases is in total concordance with the expected behavior of the transition from laminar to turbulent flow in a fixed bed. In the transition flow regime (approximately 110 < Re < 300 for fluid flow through a fixed bed),25 turbulent results must be carefully studied due to the inability of RANS models to predict the transition from laminar to turbulent flow.32 Estimates of effective radial conductivity (kr/kf) values from the CFD numerical results also agree well with accepted correlations, such as those proposed by Yagi and Kunii33 or Yagi and Wakao34 (see Figure 11). Effective radial conductivity is strongly influenced by the flow velocity profile. Other factors, like the geometrical characteristics of the bed or wall coupling functions, may influence the radial conductivity estimates, but in this case their influence is not as strong as in the estimation of Nuw. CFD results obtained using the Spalart-Allmaras model seem to agree with the empirical correlations better than the results obtained with the k- model for the turbulent solution. A possible explanation for this must be based on two facts:
Ind. Eng. Chem. Res., Vol. 43, No. 22, 2004 7055
(i) It is well known that the standard k- model does not work in the near-wall region. Normally, for modeling near-wall flow with a k- model, it is necessary to use either a wall function or a two-layer modeling (TLM) scheme. The choice of one or the other option is governed by the analysis of the y+ parameter. The version of the Spalart-Allmaras model used incorporates coupling between wall functions and damping functions. This model automatically discriminates the use of a wall function or a damping function for modeling near-wall flow, which makes it unnecessary to clearly define a priori the mesh density near the wall. Studies in this area indicate that values of 30 < y+ < 60 allow the use of a wall function and values of 1 < y+ < 5 allow the use of a TLM.24 (ii) The standard k- model presents a stagnation point anomaly. This means that in stagnation points kinetic energy is overpredicted, affecting the heat transfer evaluation in these zones. The results obtained with the Spalart-Allmaras model allowed y+ values near the wall to be analyzed before using the k- model. Values of y+ for the mesh obtained were in the range of 4 < y+ < 12, making the mesh inappropriate for the use of a wall function or a TLM for near-wall modeling. Several attempts at varying the mesh density were made, without obtaining any important change in y+ values. The fact that there is no good coupling for the k- model explains the better agreement obtained using the Spalart-Allmaras model. Conclusions As a design tool for fixed bed equipment, CFD proves to be useful in estimating wall to fluid heat transfer parameters, and also for calculating pressure drop along the bed. It was possible to model a realistic case of a fixed bed using a nonregular arrangement of spheres and including contact points on the surfaces involved in the geometry. There is an increased difficulty in modeling a nonregular arrangement of spheres because of the complexity involved in geometry generation (i.e., defining coordinates of center points of spheres). However, working in a randomized fashion allows the possibility of building a nonregular arrangement of spheres by extrapolating from the 2D circle packing theory to 3D geometries. Mathematicians have been researching this subject, and circle packing theory works as a resource for the construction of nonregular packing schemes for CFD simulation of packed beds. The inclusion of real contact points is also possible in geometry construction and mesh design. Controlling node distribution over edges and surfaces allows a correct mesh to be defined, avoiding element skewness and highly variable gradients over calculations. The definition of a good mesh allows fluid dynamics variables to be calculated, such as velocity and pressure, by solving the Navier-Stokes 3D equations. However, in our case, the proposed geometry has a greater influence than mesh density and element size in the near-wall area, and this fact affects the definition of an appropriate y+ parameter when a near-wall treatment function has to be applied. It is not an easy task to define an adequate y+ for correct coupling for the standard k- model under nearwall treatment. The y+ parameter is crucial during the selection of the appropriate turbulence model to be applied in the simulation. Good near-wall modeling is
fundamental to obtain more accurate results in pressure drop and heat transfer calculations, and the selection of the right turbulent model will depend on the geometry proposed and the values of y+ at the wall. The results obtained for all of the cases studied (laminar and turbulent) agree among themselves and with the selected empirical and semiempirical correlations when analyzing pressure drop and effective radial conductivity. This can be explained by the similarity in the velocity field obtained for each simulation. The calculation of these parameters is more closely related to velocity fields than to mixing parameters. The prediction of the mixing rate within the bed along with the near-wall treatment appreciably affects the estimate of Nuw. Flow regime zones can be identified using the heat transfer coefficient estimate. A laminar solution overestimates the value of this coefficient in the turbulent flow zone, and turbulent solutions tend to underestimate the value of the coefficient in the laminar transition zone. Turbulence models used for simulations do not predict the transition regime well, and this can be seen in the discrepancy between numerically obtained results and correlations in the low Re and transition range. Results obtained using the Spalart-Allmaras turbulence model show better agreement than those obtained using the standard k- turbulence model for estimating pressure drop and heat transfer parameters. This could be due to the fact that the Spalart-Allmaras model uses damping functions for near-wall treatment and does not present the stagnation point anomaly. The calculated velocity profiles fit the expected results qualitatively, and the calculated values of pressure drop along the bed fit in quite well with previously published data and accepted correlations. Flow structures within the bed (i.e., wall channeling, stagnant points, eddy flows) are easily identifiable. Temperature profiles obtained from inside the bed allowed wall heat transfer parameters such as Nuw and kr/kf to be estimated. These estimated parameters are similar to those obtained using well-known empirical correlations that take pressure drop and geometrical characteristics into account as important factors in the calculations. To ensure the applicability of CFD to fixed bed reactor design through flow and heat transfer simulation, a parametrical study must be carried out to determine the sensitivity of the response of the numerical method to the variation of important factors such as velocity, temperature gradients, bed porosity, etc. Acknowledgment A fellowship to A. Guardo from the FI program (AGAUR, Generalitat de Catalunya, Spain) is acknowledged. Funding from the Spanish Ministry of Science and Technology (Madrid, Spain, Grant AGL2003-05861) is also appreciated. Nomenclature Bi ) Biot number [hwRt/kr] Br ) Brinkman number [µu2/kfT] Cp ) specific heat, J kg-1 K-1 D ) vessel diameter, m dp ) particle diameter, m Ec ) Eckert number [Br/Pr] Eu ) Euler number [dP/Fu2] ˆfb ) dimensionless body forces
7056
Ind. Eng. Chem. Res., Vol. 43, No. 22, 2004
Fr ) Froude number [u2/Lg] g ) Gravity forces, m s-2 hw ) wall heat transfer coefficient, W m-2 K-1 k ) kinetic energy, m2 s-2 kf ) fluid conductivity, W m-1 K-1 kr ) radial conductivity, W m-1 K-1 kˆ ) dimensionless thermal conductivity kˆ T ) dimensionless turbulent thermal conductivity L ) characteristic length, m Nuw ) wall Nusselt number [hwdp/kf] P ) pressure, Pa pˆ ) dimensionless pressure Pr ) Prandtl number [Cpµ/kf] PrT ) Prandtl turbulent number Rt ) column radius, m Re ) Reynolds number [dpuF/µ] ReT ) Reynolds turbulent number St ) Strouhal number [L/t0u] T ) temperature, K T ˆ ) dimensionless temperature ˆt ) dimensionless time t0 ) characteristic time, s u ) velocity, m s-1 b uˆ ) dimensionless velocity vector u0 ) characteristic velocity, m s-1 y+ ) characteristic dimensionless parameter in boundary layer flow Greek letters µ ) viscosity, kg m-1 s-1 µˆ ) dimensionless viscosity µˆ T ) dimensionless turbulent viscosity vˆ ) dimensionless kinematical viscosity F ) density, kg m-3 Fˆ ) dimensionless density τc ) stress tensor, N m-2
Literature Cited (1) Pordal, H. S.; Matice, C. J.; Fry, T. J. Computational Fluid Dynamics: a Key Analytical Tool. Hydrocarbon Process. 2001, August, 85-91. (2) Joshi, J.; Ranade, V. Computational Fluid Dynamics for Designing Process Equipment: Expectations, Current Status, and Path Forward. Ind. Eng. Chem. Res. 2003, 42, 1115-1128. (3) Leal, L. G. Challenges and Opportunities in Fluid Mechanics and Transport Phenomena. In Perspectives in Chemical Engineering Research and Education; Colton, C., Ed.; Advances in Chemical Engineering 16; Academic Press: San Diego, CA, 1991; pp 61-79. (4) Cybulski, A.; Eigenberger, G.; Stankiewicz, A. Operational and Structural Nonidealities in Modeling and Design of Multitubular Catalytic Reactors. Ind. Eng. Chem. Res. 1997, 36, 31403148. (5) Dixon, A.; Cresswell, D. Theoretical Prediction of Effective Heat Transfer Parameters in Packed Beds. AIChE J. 1979, 25, 663-676. (6) Tobis, J.; Zio´lkowski, D. Modeling of Heat Transfer at the Wall of a Packed-Bed Apparatus. Chem. Eng. Sci. 1988, 43, 30313036. (7) Yin, F.; Wang, Z.; Afacan, A.; Nandakumar, K.; Chuang, K. T. Experimental Studies of Liquid Flow Maldistribution in a Random Packed Column. Can. J. Chem. Eng. 2000, 78, 449457. (8) Manz, B.; Gladden, L. F.; Warren, P. B. Flow and Dispersion in Porous Media: Lattice-Boltzmann and NMR Studies. AIChE J. 1999, 45, 1845-1854. (9) Suekane, T.; Yokouchi, Y.; Hirai, S. Inertial Flow Structures in a Simple-Packed Bed of Spheres. AIChE J. 2003, 49, 10-17.
(10) Allen, M.; Furtado, F. In Fluid Transport in Porous Media; du Plessis, J. P., Ed.; Advances in Fluid Mechanics; Computational Mechanics Publications: Southampton, U.K., 1997; pp 255-302. (11) Prakash, M.; Turan, O ¨ . F.; Li, Y.; Thorpe, G. R. CFD Modeling of Natural Convection Heat and Mass Transfer in Hygroscopic Porous Media. Drying Technol. 2000, 18, 2175-2201. (12) Dalman, M. T.; Merkin, J. H.; McGreavy, C. Fluid Flow and Heat Transfer Past Two Spheres in a Cylindrical Tube. Comput. Fluids 1986, 14, 267-281. (13) Lloyd, B.; Boehm, R. Flow and Heat Transfer around a Linear Array of Spheres. Numer. Heat Transfer, Part A 1994, 26, 237-252. (14) Vortmeyer, D.; Schuster, J. Evaluation of Steady Flow Profiles in Rectangular and Circular Packed Beds by a Variational Method. Chem. Eng. Sci. 1983, 38, 1691-1699. (15) Kuipers, J. A. M.; Van Swaaij, W. P. M. In Computational Fluid Dynamics Applied to Chemical Reaction Engineering; Wei, J., Anderson, J., Bischoff, K., Denn, M., Seinfield, J., Stephanopoulos, G., Eds.; Advances in Chemical Engineering 24; Academic Press: San Diego, CA, 1998; pp 227-328. (16) Logtenberg, S. A.; Nijemeisland, M.; Dixon, A. G. Computational Fluid Dynamics Simulations of Fluid Flow and Heat Transfer at the Wall-Particle Contact Points in a Fixed-Bed Reactor. Chem. Eng. Sci. 1999, 54, 2433-2439. (17) Dixon, A. G.; Nijemeisland, M. CFD as a Design Tool for Fixed-Bed Reactors. Ind. Eng. Chem. Res. 2001, 40, 5246-5254. (18) Guedes de Carvalho, J. R. F.; Delgado, J. M. P. Q. Lateral Dispersion in Liquid Flow through Packed Beds at Pem < 1400. AIChE J. 2000, 46, 1089-1095. (19) Guardo, A.; Larrayoz, M. A.; Recasens, F. In Regular Packing Types for CFD Simulation of SCF Extraction and Reaction Equipment; Brunner, G., Kikic, I., Perrut, M., Eds.; Proceedings of the 6th International Symposium on Supercritical Fluids; Institut National Polytechnique: Lorraine, France, 2003; pp 643648. (20) Graham, R. L.; Lubachevsky, B. D.; Nurmela, K. J.; O ¨ stergard, P. R. J. Dense Packings of Congruent Circles in a Circle. Discrete Math. 1998, 181, 139-154. (21) Melissen, H. Densest Packing of Eleven Congruent Circles in a Circle. Geom. Dedic. 1994, 50, 15-25. (22) Reid, R.; Prausnitz, J.; Poling, B. The Properties of Gases and Liquids; McGraw-Hill: Boston, MA, 1987; pp 95-205. (23) Yaws, C. Chemical Properties Handbook; McGraw-Hill: New York, 1999; pp 1-55. (24) Coussirat, M. Ph.D. Thesis, Universitat Polete`cnica de Catalunya, Barcelona, Spain, 2003. (25) Jolls, K. R.; Hanratty, T. J. Transition to turbulence for flow through a dumped bed of spheres. Chem. Eng. Sci. 1966, 21, 1185-1190. (26) Tobis, J.; Zio´lkowski, D. Modelling of heat transfer at the wall of a packed-bed apparatus. Chem. Eng. Sci. 1988, 43, 30313036. (27) Fluent, Inc. Fluent 5 User’s Guide, 1998. (28) Ergun, S. Fluid flow through packed columns. Chem. Eng. Prog. 1952, 48, 89-94. (29) Wakao, N.; Kaguei, S. Heat and Mass Transfer in Packed Beds; Gordon and Breach: New York, 1982; pp 161-204. (30) Olbrich, W.; Potter, O. Mass transfer from the wall in small diameter packed beds. Chem. Eng. Sci. 1972, 27, 1733-1743. (31) Li, C.; Finlayson, B. A. Heat transfer in packed beds - a reevaluation. Chem. Eng. Sci. 1977, 32, 1055-1066. (32) Durbin, P. A.; Petterson Reif, B. A. Statistical Theory and Modeling for Turbulent Flows; John Wiley & Sons: New York, 2001; pp 138-139. (33) Yagi, S.; Kunii, D. Studies in heat transfer near wall surface in packed beds. AIChE J. 1960, 6, 97-104. (34) Yagi, S.; Wakao, N. Heat and mass transfer from wall to fluid in packed beds. AIChE J. 1959, 5, 79-85.
Received for review October 31, 2003 Revised manuscript received July 19, 2004 Accepted August 13, 2004 IE034229+
