5122
Ind. Eng. Chem. Res. 2010, 49, 5122–5131
On the Effect of Subgrid Drag Closures Sofiane Benyahia* National Energy Technology Laboratory, Morgantown, West Virginia 26505
The effect of two subgrid drag closures on the flow of air and Geldart group A particles is presented in this study. A subgrid drag model based on fitting simulation data obtained from finely resolved simulations and a drag model based on the energy minimization approach are both used to solve a gas-solids flow in the riser section of a circulating fluidized bed. The numerical results using a coarse computational grid obtained with these subgrid models are compared with those using a standard drag model as well as experimental data obtained in a pilot-scale riser. Numerical predictions using both subgrid models showed higher solids holdup in the riser indicated by the radial solids density and axial pressure drop profiles in the 2D and 3D system geometries considered in this study. These subgrid models are demonstrated to be both needed and useful as large-scale numerical simulations commonly use coarse computational grids that are unable to resolve the smallest heterogeneous structures observed in the fluidization of small particles. Introduction Current continuum gas-solids flow models rely on two major constitutive relations: one for the solids-phase stresses derived from granular kinetic theory1-4 and another for gas-solids momentum exchange. Traditionally, gas-solids friction coefficients have been expressed using semiempirical correlations such as the well-known Wen-Yu and Ergun correlations,2 but recently, similar expressions were fitted to simulation data obtained using physical models based on first principles.5-9 These recently derived drag models were obtained by finely resolving the flow of the fluid around the particles, and the friction or drag can be obtained by integrating the fluid viscous stress acting on the particles. These continuum models are able to predict qualitative gas-solids flow phenomena, such as core-annulus structures in circulating fluidized beds, bubble formation, and bed expansion in bubbling fluidized beds. However, the quantitative numerical predictions of these models are not very accurate, especially for the fluidization of small particles of Geldart10 group A classification of powders, such as fluid catalytic cracking (FCC) catalyst. This has been recognized in the fluidization literature11,12 for FCC particles, and the proposed simple solution to this issue is to use an effective agglomerate size, of the order of a few particle diameters, in the standard drag formula in order to match experimental data for pressure drop or bed height expansion. The reason for these discrepancies between numerical results and experimental data can be traced to the fact that simulation results of large-scale fluidized beds are usually conducted with coarse computational grids13 where the continuum models are discretized and solved. Furthermore, fully resolved numerical simulations of even small-scale fluidized beds are still prohibitively expensive for Geldart A particles.14 By using fine computational grids of the order of a few particle diameters, it is possible to resolve small-scale heterogeneous structures, such as small clusters and bubbles in fluidized beds, that are usually unresolved with coarser computational grids.13,14 The origin of these heterogeneous structures that occur in fluidization is due to local instabilities and was summarized by Agrawal et al.13 as (1) relative motion between the fluid and particles, which can also be traced to the nonlinearity of the gas-solids drag correlation15 and (2) energy dissipation due to inelastic collisions * To whom correspondence should be addressed. E-mail:
[email protected]. 10.1021/ie900658k
between particles and viscous damping. The origin of small agglomerates of particles that occur in fluidization of Geldart A particles may also be due to the fact that these small particles are slightly cohesive, requiring corrections to some constitutive relations.16,17 However, a recent study that included the effect of cohesive forces resulted in an effective agglomerate size a few times greater than the individual particle diameter.18 Although the use of a relatively narrow agglomerate size seems adequate for the dense bubbling regime resulting in a simple correction to drag correlations,12 it has been shown that more complex corrections to the standard drag correlations are required for more dilute flows, such as those in risers.13,19 The purpose of this study is to investigate the effect of subgrid corrections to gas-solids friction coefficients by conducting numerical simulations using coarse grids and comparing results with those obtained using a standard drag correlation as well as with experimental data obtained in the Particulate Solid Research Inc. (PSRI) riser for the flow of air and FCC particles. The subgrid drag models obtained recently by Igci et al.20 and a modified drag model obtained using the energy minimization multiscale (EMMS) technique described by Yang et al.19 were used in this study. Gas-Solids Flow Model The governing equations, constitutive relations, and boundary conditions are summarized in Table 1. The kinetic theory of granular flows1-3 was modified to include the effect of the interstitial fluid in the expression of solids viscosity and conductivity.4,13,21 These modifications describe the correct behavior of the solids viscosity and conductivity for very dilute flow conditions.22 The granular pressure and collisional dissipation are obtained using standard kinetic theories,1-3 and the radial distribution function at contact (eq 13) can be simply derived from the compressibility factor presented by Carnahan and Starling.23 The rate of dissipation of granular energy due to viscous damping was modeled following the derivation of Gidaspow,2 and the production due to gas-solids slip shown in eq 21 was modeled following a simplified expression presented by Agrawal et al.13 using their eq 16. Near the solids packing limit, the kinetic theory of granular flows is not valid because enduring contacts between particles occur. In this dense flow regime, the frictional pressure and viscosity derived by Srivastava and Sundaresan24 were used in this study and were
This article not subject to U.S. Copyright. Published 2010 by the American Chemical Society Published on Web 07/23/2009
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
5123
Table 1. Gas-Solids Hydrodynamic Model Governing Equations Conservation of mass for phase m (s solids and g gas): ∂(Fmεm) + ∇ · (Fmεmvm) ) 0 ∂t Conservation of linear momentum:
[
(1)
]
∂(Fmεmvm) + ∇ · (Fmεmvmvm) ) -εm∇Pg + ∇ · (τm) - βmn(vm - vn) + εmFmg ∂t
Translational granular energy conservation equation: 3 ∂(εsΘs) F + ∇ · (εsΘsvs) ) -∇ · q + τc : ∇vs - Js + ΠΘ 2 s ∂t
[
]
(2)
(3)
Interphase Momentum Exchange Term (βgs ) βsg ) β) Wen-Yu correlation:
3 Fgεgεs |vg - vs | -2.65 β ) CD εg 4 dp CD )
{
Fgεg |vg - vs |dp 24/Re(1 + 0.15Re0.687) Re < 1000 , Re ) µg 0.44 Re g 1000
EMMS drag model:
{
150
β)
where
{
(4)
εs(1 - εg)µg εgdp2
Fg + 1.75 εs |vg - vs | εg < 0.7 dp (6)
3 Fgεgεs |vg - vs | C f(εg) 4 D dp
εg g 0.7
-15.63 + 16.63εg
f(εg) )
(5)
516.1 - 2243εg + 3658εg2 - 2652εg3 + 721.6εg4
εg > 0.94 0.94 e εg > 0.84
17127 - 82203εg + 148026εg2 - 118519εg3 + 35599εg4
0.84 e εg > 0.79
(7)
-39898 + 199506εg - 372449εg2 + 307602εg3 - 94803εg4 0.79 e εg > 0.72 0.72 e εg > 0.7 -276.9 + 764.0εg - 520.1εg2
{ {
Filtered drag model: for dimensionless filter size (or computational grid size): g∆/Vt2 ) 2.056 4.277εg4 - 2.363εg3 - 4.324εg2 - 2.474εg + 4.904 βVt ) Fsεsεgg
εg2 - 5.411εg + 4.427
0.4820 εg3.574
0.04083 < εs e 0.2589
148.0εs3 - 124.8εs2 + 37.50εs - 3.642 εgεs
εs > 0.2589
For dimensionless filter size (or computational grid size): g∆/Vt2 ) 4.112 5.301εg3 - 5.473εg2 - 5.723εg + 5.906 βVt ) Fsεsεgg
εs e 0.04083
εg3 - 3.820εg2 + 1.145εg + 1.686 0.3349 εg3.872
(8)
εs e 0.0533 0.0533 < εs e 0.2922
(9)
-839.0εs + 1574εs - 1012εs + 277.3εs - 27.51 εs > 0.2922 εgεs 4
3
proven reliable by recent studies.25,26 For the solids phase, the frictional-collisional wall boundary conditions derived by Johnson and Jackson27 were used in this study along with a free-slip wall boundary condition for the gas phase. Following previously published studies,2,12-15 the gas-phase stress is modeled using only the laminar viscosity. Furthermore, a recent study by Benyahia et al.22 showed that, under moderately dense flow conditions, the gas-phase turbulent stresses could be neglected without any noticeable effect on the time-averaged flow profiles.
2
The standard Wen-Yu drag correlation is considered in this study because it is commonly used in the fluidization literature either in itself or in combination with the Ergun equation.2 It will also serve as a base case for comparison with subgrid drag models to highlight their effect on the flow behavior. The filtered drag models presented in this study were fitted to data obtained by Igci et al.20 and were obtained from Igci and Sundaresan (private communications, 2008). These drag closures were obtained using finely resolved simulations with the same flow model described in Table 1 along with Wen-Yu drag correla-
5124
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
Table 1. Continued Constitutive Relations for Gas and Solids Phases 1 1 Sm ) (∇vm + (∇vm)T) - ∇ · vmI 2 3
Gas-phase stress tensor: τg ) 2µgSg, where
Solids-phase stress tensor:
τs ) [-(Ps + Pf) + ηµb∇ · vs]I + 2(µs + µf)Ss
(11)
Ps ) εsFsΘs[1 + 4ηεsg0]
(12)
Solids kinetic-collisional pressure: Radial distribution function at contact:
1 - 0.5εs
g0 )
Solids viscosity model:
µs )
(10)
(13)
(1 - εs)3
( 2 +3 R )[ g η(2µ*- η) (1 + 58 ηε g )(1 + 58 η(3η - 2)ε g ) + 53 ηµ ] s 0
s 0
(14)
b
0
[
µ* ) µ 1 +
2βµ (εsFs)2g0Θs
]
-1
Granular energy flux and conductivity: κs )
5 F d πΘ , 96 s p√ s
µ)
( )[( κ* g0
1+
256 2 µε g 5π s 0
(15)
)(
(
)
6βκ 5(εsFs)2g0Θs
Frictional solids pressure and viscosity: Pf ∇ · vs ) 1Pc n√2 sin(δ)√Ss : Ss + Θs /dp2
{
(16)
12 12 64 ηε g 1 + η2(4η - 3)εsg0 + (41 - 33η)η2(εsg0)2 5 s 0 5 25π
[
)
n-1
,
]
-1
µf )
,
κ)
Fr
(εs - εsmin)r (εsmax
- ε s)
s
Collisional dissipation of granular energy and viscous damping terms: Fsεs2g0 3/2 48 η(1 - η) Θs , Js ) dp √π
(17)
{
(18)
() }
Pf Pf sin(δ) n - (n - 1) P 2 c √2 √S : S + Θ /d s s s p
εsmax g εs > εsmin , εs e εsmin
0
]
75Fsdp√πΘs 48η(41 - 33η)
1025(εs - εsmax)10 εs > εsmax
Pc )
µb )
q ) -κs∇Θs
κ* ) κ 1 +
where
,
and
n)
{
√3 sin(δ) ∇ · vs g 0 2 ∇ · vs < 0 1.03
ΠΘ ) -3βΘs + 81
εsµg2 |vg - vs | 2 g0dp3Fs√πΘs
1/n-1
(19)
(20)
(21)
Gas-Solids Wall Boundary Conditions Gas phase: free slip Solids phase: Johnson-Jackson partial slip27
φπFsεsg0√Θs vsl · (τk + τf) · n + vsl + (n · τf · n)tan δw ) 0 |vsl | 2√3εsmax
n·q)
φπ|vsl | 2Fsεsg0√Θs 2√3εsmax
tion. The study of Agrawal et al.13 showed an increase in gas-solids slip velocity as the computational grid was refined and more flow structures were resolved, which prompted these authors to propose filtered drag models to be used with coarsegrid simulations of gas-solids flows. The interested reader is referred to the published literature13,20 for more details. Also, the gas and particles physical properties used in this study were similar to those used by Igci et al.20 with the exception of particle-particle restitution coefficient that was slightly different. As indicated by Igci et al.,20 the most important factor in determining the filtered drag coefficient is the averaged solids volume fraction as the dependence on Reynolds number is weak.
-
√3πFsεsg0(1 - ew2)√Θs 4εsmax
Θs
(22)
(23)
Thus, it is not important in this study to use similar gas and solids flow rates as those used by Igci et al.20 Equation 8 expresses the filtered drag model used in the 2D riser simulation corresponding to a dimensionless grid size g∆/Vt2 ) 2.056 or a grid size ∆ ) 0.014 m, and eq 9 was used with the relatively coarser grid for the 3D riser corresponding to a dimensionless grid size g∆/Vt2 ) 4.112 or a grid size ∆ ) 0.028 m. It should be noted that these filtered drag closures were obtained in fully periodic domains and that corrections due to the effect of wall boundaries are currently unavailable as will be discussed later in this paper. Finally, eq 7 shows the correction to the EMMS drag model that applies to a modified Wen-Yu drag correlation
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
5125
Figure 1. Instantaneous values of friction coefficient obtained throughout the 2D riser for the three drag relations used in this study. The parameters used to make drag correlations dimensionless were also used by Igci et al.20
as shown in eq 6. Similar to the approach of Yang et al.,19 the correction to the drag applies only to the relatively dilute flow conditions as denser flow regions are modeled using the Ergun equation. This approach is a simplification of the original EMMS approach as it does not require a solution of a set of nonlinear equations and the subsequent optimization.19 Thus, it is relatively easy to apply this technique in a computational fluid dynamics (CFD) model with minimum changes to computer code and with low additional computational cost. However, it should be noted that the simplified EMMS model presented by Yang et al.19 can only be used for specific gas and solids flow rates. Therefore, the EMMS model presented in this study by eq 7 is different from that of Yang et al.19 and applies to the PSRI riser under the flow conditions and gas-solids physical properties presented in the next section (The current EMMS model was fitted to data obtained by Yang at the Institute of Process Engineering (www.ipe.ac.cn), Chinese Academy of Sciences, 2008). In order to show quantitative differences between these drag models, Figure 1 provides instantaneous data obtained using the Wen-Yu correlation (eq 4), EMMS model (eqs 6 and 7), and the filtered drag model (eq 8) in a 2D riser. The parameters used to make these drag expressions dimensionless were taken from the study of Igci et al.20 Thus, the scatter in the data obtained using the Wen-Yu and the EMMS approaches is an indication that these models depend not only on the solids volume fraction but also on the slip velocity. Through most of the range of solids volume fraction (between 5% and 25%), the EMMS model shows the lowest values of friction coefficient, indicating that it would probably predict the highest solids holdup and pressure drop in the riser. The filtered drag model is roughly a factor of 2 lower than the Wen-Yu correlation and is, thus, expected to yield higher solids hold-up in the riser. The values of the EMMS drag model are shown to be higher than those of the Wen-Yu correlation for solids fractions higher than 25%, which is anomalous behavior due to the simplifications to the EMMS approach presented by Yang et al.19 Furthermore, it should be noted that the EMMS drag model presented by Yang et al.19 in their Figure 3b does not correspond to their “model B” (see their Table 2) but was rather obtained using the full EMMS approach. Nevertheless, this anomalous behavior occurs only at high solids volume fraction, and thus,
Table 2. Gas and Solids Physical Parameters process temperature process pressure air density air viscosity solids density particle diameter particle-particle restitution coefficient particle-wall restitution coefficient specularity coefficient void fraction at maximum packing
300 K 101325 Pa computed using ideal gas law (about 1.2 kg/m3) 1.8 × 10-5 Pa · s 1712 kg/m3 76 µm 0.95 0.7 0.0001 0.4
the simplified EMMS model provides reasonable results for the riser flow case presented in this study. Simulation Conditions The physical properties of FCC particles (Geldart type A powder) and air are summarized in Table 2. Some physical properties of the FCC particles used in this study were estimated from previously published work. Particle-particle and particlewall restitution coefficients as well as the specularity coefficient used in the wall boundary conditions were taken from the work of Benyahia et al.22 The value of the particle-particle restitution coefficient for FCC particles varies in the open literature (typically from 0.9-0.99), but a value of 0.95 has been used by several researchers12,29-31 in the past. The system geometry considered in this study is a 0.2 m diameter and 14.2 m length vertical riser located at the PSRI experimental research facility in Chicago (www.psrichicago.com). The description of the PSRI riser section of the circulating fluidized bed (CFB) is shown in Figure 2 for the 2D Cartesian and the 3D cylindrical geometries, and other details of the experimental setup are available elsewhere.28,29 Similar to previous studies,29-31 the solids enter the system from both sides for the 2D case although, in this study, the solids inlets (each 0.1 m wide) were extended 0.2 m horizontally, as shown in Figure 2, to improve gas-solids mixing in the inlet region of the 2D riser.32 Gas enters from the bottom of the riser, and gas and solids can leave the system from the top of the riser, where 50% of the available exit area was blocked symmetrically by setting walls in order to avoid excessive gas backmixing. The computational grid consists of 14 cells along the width of a 0.2 m riser and 1017 cells along the height of the 2D riser resulting in a cell size of roughly 1.4
5126
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
Figure 2. Schematic drawing of the 3D cylindrical and 2D Cartesian geometries of the PSRI riser with detailed inlet and outlet designs used in the computer simulations.
cm with a width-to-height ratio equal to 1. A 3D cylinder is also considered in this study to better approximate the real geometry of the experimental riser. Only one solids inlet and one outlet, set on the side of the riser, were considered for this 3D case, and the main gas flow was fed from the bottom of the riser. The solids inlet consists of a patch centered at a height of 0.15 m with a diameter of 0.1 m; it azimuthally extends from π/5 to 4 π/5. The outlet patch was centered at a height of 14 m with a diameter of 0.2 m; it azimuthally extends from π/4 to 3 π/4. The computational grid for the 3D riser consists of 4 cells along the radial direction (half the diameter), 10 cells azimuthally, and 508 uniform cells along the axial direction, which results in a cell size of about 2.5-3 cm (roughly doubled the cell size for the 2D case). Initially, some solids are present in the system (1% per volume) with gas and solids velocity components set to zero and granular temperature set to 0.004 m2/s2. The main gas flow from the bottom of the riser is set at a velocity of 5.2 m/s. At the main solids inlet, the solids volume fraction is set at 40%, and a mass flux of 489 kg/s · m2 is achieved in the riser. The solids inlet velocity is set to 0.714 m/s for each of the two solids inlets in the 2D geometry and is set to 1.42 m/s in the 3D geometry. The flow is isothermal at an ambient temperature of 300 K. The MFIX CFD code (https://mfix.netl.doe.gov) was used to solve the model equations shown in Table 1 following the finite volume numerical technique. The convective terms for all
governing equations were discretized using the superbee second order scheme along with deferred correction.33 The set of discretized and linearized equations are then solved by a linear solver to a low tolerance (10-8). The nonlinear numerical residuals for all equations are set to 10-4, and the gas and solids pressure correction residuals are normalized using the first iteration residual at each time-step. An adaptive time stepping algorithm was used in this study with an average time-step of about 2 × 10-4 s. Results and Discussion Simulation results presented in this section consist of numerical data obtained in 2D followed by results in the 3D riser. Simulations are typically run for 15 s, and then time-averaged data are collected for an additional 30-50 s at a frequency of 50 Hz. The time-averaging period is judged sufficient when horizontal profiles of flow variables are symmetrical and axial profiles of the gas pressure gradient are relatively smooth. The overall gas-solids instantaneous flow patterns in the 2D riser are shown in Figure 3. The solids concentration in the two horizontally extended solids inlets is high due to the fact that no aeration is added in that flow region and that some solids occasionally fall below these inlets and can reach the bottom of the riser where the main gas flow is fed. Large clusters of solids are present throughout the length of the riser and tend to congregate near the walls where the momentum of the carrier
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
Figure 3. Instantaneous (after about 30 s of simulation) contour profile of solids volume fraction in the 2D riser (to scale, showing the bottom half and top sections separately) for the three drag models used in this study. Red and blue indicate dense regions with a maximum solids volume fraction of about 0.5 and dilute regions, respectively.
gas is usually lower. It is apparent in Figure 3 that the predictions using the standard Wen-Yu drag model result in more dilute flow. The solids concentration in the riser increases with the use of the filtered drag model (eq 8), and even higher solids concentrations were obtained using the EMMS drag model. This is clearly shown in Figure 3 as the EMMS drag model predicts larger and denser clusters not only in the lower sections but also near the top regions of the 2D riser. The flow profiles in the riser may be asymmetrical at a given instant as clusters of solids tend to form on one side and then on the other side of the riser as found in previous studies,22 but the time-averaged profiles of flow variables are symmetrical as presented next, even though the experimental measurements were only conducted along one half of the riser diameter28 (experimental data plotted along the full width of the riser assumed symmetry). Figure 4 shows the time-averaged solids apparent density (Fsεs) and axial mass flux profiles along the width of the riser at a height of 3.9 m above the solids side inlets where experimental measurements are available. The choice of using a standard drag or subgrid drag model has a profound effect on the solids hold-up in the riser. Subgrid models that include the effect of unresolved agglomerates show large values of solids hold-up in the riser as previously indicated by Figure 3. The high solids concentration near the walls of the riser is correctly predicted by both the filtered and the EMMS drag models. The low solids concentration at the walls predicted by the standard Wen-Yu drag model was obtained in a previous study31 using a slightly different drag model.
5127
Figure 4. Time-averaged profiles of solids apparent density and axial mass flux at a height of 3.9 m above the solids inlet computed using the three drag models used in this study and compared with experimental measurements obtained at the same height in the riser.
However, both subgrid drag models overpredict, by approximately the same amount, the solids concentration at the center of the riser. Such behavior was computed previously using a similar EMMS drag model for the PSRI riser under the same flow conditions,29 and a grid refinement study showed that similar solids concentrations were obtained at the center of the riser. Computationally, it seems that such high solids concentrations are necessary in order to predict the high solids mass flux at the core of the riser. In fact, the low solids concentration predicted by the standard drag model, albeit higher than the experimental measurements, has a direct effect on the low solids mass flux predicted at the core of the riser. The large downflow of solids near the walls of the riser, shown by the negative solids mass flux in Figure 4, is due to the high solids concentrations in that flow region where the weight of particles locally exceeds the gas pressure drop, which is the main driving force of the gas-solids flow. Figure 5 shows the time-averaged gas and solids axial velocity profiles computed along the width of the riser at a height 3.9 m above solids inlets. The standard Wen-Yu drag model predicts a small slip velocity, which averages to 0.23 m/s, approximately equal to the terminal velocity Vt ) 0.26 m/s calculated using the Khan and Richardson34 correlation. On the other hand, the averaged slip velocity computed using the EMMS drag model was about 1.13 m/s, which is greater than four times the individual particle terminal velocity. This is an indirect proof of the formation of clusters of particles (a direct proof is indicated in Figure 3) that cause gas to bypass them, which increases the slip velocity between the gas and solids phases. The downflow of solids at the walls is indicated by the negative velocity computed near the walls.
5128
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
Figure 5. Time-averaged gas and solids axial velocity profiles along the width of the riser at a height of 3.9 above the solids inlet computed using the standard Wen-Yu drag correlation and the EMMS approach.
The low magnitude of this negative velocity cannot account for the large negative flux shown in Figure 4. This behavior was found in earlier work22,35 and is due to the fact that the time-averaged solids flux is not equal to time-averaged velocity times time-averaged solids density (FsεsVs * Fsεs Vs), so that the solids volume fraction and velocity correlate in the presence of clusters. Figure 6 shows the time-averaged gas pressure gradient profiles along the height of the 2D riser at a location close to the wall. Fluctuations in the gas pressure gradient profiles are also an indication of the presence of large clusters and can be smoothed further by taking longer averaging periods. In fact, the longest averaging period used (50 s) was for the EMMS simulation case that, as shown in Figure 3, produced the largest clusters. The lowest pressure gradients along the riser height were computed using the standard Wen-Yu drag model that clearly underpredicts the experimental measurements. The most accurate predictions were obtained using the filtered drag model for both the magnitude and shape of the pressure gradient’s profile. The relatively high values of pressure gradients measured at the bottom of the riser are an indication of large solids hold-up in that region, and a flat profile at the upper region of the riser is referred to as a “fully developed” flow region. The pressure drop in the upper region is accurately predicted by most drag models used in this study, although the EMMS model predicts slightly higher pressure gradients due to the large clusters computed in that flow region (see Figure 3). The rest of this section is devoted to results obtained in a 3D cylinder with a solids inlet and outlet design similar to those used experimentally but with a coarser mesh than that used in the 2D riser case presented earlier. It should be noted, however, that although the size of the computational grid is double that used in 2D, the number of cells in the 3D simulation is more than 20 000 and takes about 2 months to compute 30 s of simulation time using an Intel Xeon 2.66 GHz CPU. It should also be pointed out that such a coarse mesh may be considered extremely fine for simulations of large-scale industrial risers. For example, a typical industrial riser of 1-2 m diameter and 30 m height would require a prohibitively large mesh (1-6 million cells) even with the coarse mesh used in this 3D case. So, it is essential to develop subgrid models for even coarser mesh sizes in order to be able to conduct simulations of large-scale industrial units in a reasonable time. Figure 7 shows the time-averaged solids apparent density and axial mass flux at a height of 3.9 m above the solids
Figure 6. Time-averaged gas pressure gradient profiles plotted along the height of the 2D riser. Numerical predictions obtained for the three drag models used in this study are compared with experimental data.
inlet and compared with experimental data taken at the same height. The numerical data were plotted across half the diameter of the riser and averaged azimuthally even though the results obtained for this 3D case were nearly symmetrical. The solids density near the walls was computed accurately for both subgrid models, and the standard Wen-Yu drag model underpredicts the solids concentration at the walls similar to the previous 2D findings. The difference between 2D and 3D simulation results can be seen at the core of the riser where the solids density was predicted to be lower for the 3D case and is in better agreement with experimental data. It should be noted that the filtered drag model is expressed by eq 9 and is, thus, optimized for the grid size used in this 3D case, unlike the EMMS drag model that does not depend on the grid size. Overall, the solids mass flux is underpredicted for all three drag models used in this 3D case.
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
5129
Figure 9. Comparison of time-averaged solids axial mass flux obtained in 3D and 2D riser geometries. The same coarse computational grid was used in both simulations.
Figure 7. Time-averaged profiles of solids apparent density and axial mass flux at a height of 3.9 m above the solids inlet computed using the three drag models used in this study and compared with experimental measurements obtained at the same height in the riser. Numerical simulations were conducted in 3D cylindrical coordinates with a coarse mesh.
Figure 8. Time-averaged gas and solids axial velocity profiles plotted along the radius of the riser at a height of 3.9 m above the solids inlet computed using the standard Wen-Yu and filtered drag correlations.
Near the wall, the downflow of solids is not predicted due to the fact that the annulus is not fully resolved with this coarse computational mesh. Due to mass conservation, the lack of downflow of solids in the annulus results in a low solids mass flux computed in the core region of the riser. Further proof of this behavior can be seen in Figure 8 where an upflow of solids and gas is computed in the annular region using both the Wen-Yu and the filtered drag models. The positive solids velocity at the walls resulted in a large solids mass flux at the walls due to the large solids density computed
Figure 10. Time-averaged gas pressure gradient profiles plotted along the height of the 3D riser. Numerical predictions obtained for the three drag models used in this study are compared with experimental data taken near the wall of the riser (results using Wen-Yu and EMMS drag are represented with black and gray continuous lines, respectively).
(and measured) in this flow region. Experimentally, it is possible to obtain an upflow of solids in the annular region for the case of high density and high gas velocity flow regime.36 However, the upflow computed in this 3D case is clearly due to the unresolved flow in the annular region. Figure 9 shows that both 2D and 3D simulations predict similar solids mass flux profiles for the same coarse computational grid. This is a proof that the solids upflow in the annulus is due, not to a 3D effect, but to the unresolved flow of solids in the annular region. As mentioned earlier, the coarse mesh used in this 3D simulation may still be considered as extremely fine for other cases that involve large-scale industrial risers. Therefore, it is necessary to develop near-wall corrections to the subgrid models as well as corrections to the boundary conditions that may need adjustments in the case of a coarse mesh, as recent results obtained by Igci and Sundaresan37 indicate. Another approach is to use a relatively fine radial computational mesh to resolve the annular flow region and a much coarser mesh along the height of the riser. A similar approach using computational cells of height-to-width ratio equal to 4 was successfully used for a large-scale riser with subgrid models developed specifically for such a nonuniform grid.38 Figure 10 shows the time-averaged gas pressure gradient profiles along the height of the 3D riser taken near the wall opposite to the solids inlet and outlet patches. The constant profile computed using the standard Wen-Yu drag model is
5130
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
an indication of the absence of dense clusters throughout the height of the riser. Near the top of the riser, both subgrid drag models predict higher pressure gradients due to solids backmixing that can be tracked to the abrupt blind-tee used computationally instead of the smooth elbow exit design used experimentally.32 Nevertheless, this backmixing of solids affects only the top region of the riser since solids hold-up agrees reasonably well with experimental data at the bottom region of the riser as seen in Figure 7. Summary and Conclusion In this study, we evaluated the effect of two subgrid drag models based on a filtered drag model obtained by finely resolved continuum simulations20 using the same flow model summarized in Table 1 and a simplified version of the EMMS approach19 that applies to FCC particles under specific flow conditions used in the PSRI riser. Results obtained with these subgrid models were compared with numerical data using the standard Wen-Yu drag model as well as experimental data obtained for the flow of air and Geldart type A particles. Both 2D and 3D coarse numerical simulations using the standard Wen-Yu drag model underpredicted the solids hold-up and gas pressure drop in the riser. The main effect of these subgrid drag models is an increase in the solids hold-up, which leads to an increase in gas pressure drop, and a better agreement with available experimental data obtained in the PSRI riser. It is nearly impossible to resolve all heterogeneous flow structures in large-scale industrial risers using a computational grid size of the order of a few particle diameters as proposed by previous modelers.14,22,38 For example, current coal gasification research conducted at NETL using available supercomputers at Oak Ridge National Laboratory (ORNL) with thousands of processors and a computational grid consisting of millions of cells is still not able to resolve all flow structures due to the large scale of the reactors that require billions of computational cells. Further research in this area should address the observed fact that the annular flow region in the riser must be resolved in order to quantitatively predict the large solids mass flux gradients near the wall. This can be achieved by considering numerical cells with large height-to-width aspect ratios as well as additional corrections to the subgrid drag closures near wall boundaries. Acknowledgment The author gratefully acknowledges the help received from Yesim Igci and Sankaran Sundaresan who developed the filtered drag models used in this study. He is also thankful to useful discussions and the development of the modified EMMS model specifically tailored for PSRI riser flow conditions by Ning Yang from the research group of Jinghai Li at the Institute of Process Engineering, Chinese Academy of Sciences. Nomenclature dp ) particle diameter e, ew ) particle-particle and particle-wall restitution coefficients Fr, r, s ) constant in the Srivastava-Sundaresan24 model, equal to 0.5, 2, and 5 dyn/cm2, respectively g ) acceleration of gravity g0 ) radial distribution function at contact I ) identity tensor Js ) granular energy dissipation due to inelastic collisions n ) coefficient in the frictional model n ) unit vector normal to wall surface
P ) pressure q ) flux of granular energy Sm ) strain-rate tensor Vm ) velocity vector of phase m Greek Letters R ) constant in solids viscosity model, equal to 1.6 δ, δw ) angle of internal friction and wall friction of about π/6 and π/10, respectively εs ) solids volume fraction φ ) specularity coefficient used in wall boundary condition η ) constant depending on particle restitution coefficient equal to (1 + e)/2 κ ) solids phase dilute granular conductivity κ* ) granular conductivity with effect of interstitial fluid κs ) conductivity of solids granular energy µ ) solids phase dilute granular viscosity µ* ) granular viscosity with effect of interstitial fluid µb ) bulk viscosity of the solids phase µs ) granular viscosity ΠΘ ) production/dissipation of granular energy due to solids-fluid interaction Θs ) granular temperature Fs ) solids density τ ) stress tensor Indices c ) critical f ) frictional g ) gas phase k ) kinetic-collisional m, n ) indices for the phases max ) maximum packing min ) minimum frictional solids fraction (εsmin ) 0.5) s ) solids phase sl ) particle-wall slip t ) terminal w ) wall
Literature Cited (1) Ding, J.; Gidaspow, D. A bubbling fluidization model using kinetictheory of granular flow. AIChE J. 1990, 36, 523. (2) Gidaspow, D. Multiphase Flow and Fluidization: Continuum and Kinetic Theory Description; Academic Press: San Diego, 1994. (3) Lun, C. K. K.; Savage, S. B.; Jeffrey, D. J.; Chepurniy, N. Kinetic theories for granular flow - Inelastic particles in couette-flow and slightly inelastic particles in a general flowfield. J. Fluid Mech. 1984, 140, 223. (4) Lun, C. K. K.; Savage, S. B. Kinetic theory for inertial flows of dilute turbulent gas-solids mixtures. In Granular Gas Dynamics; Poeschel, T., Brilliantov, N., Eds.; Springer-Verlag: Berlin Heidelberg, 2003; pp 624, 267. (5) Hill, R. J.; Koch, D. L.; Ladd, J. C. The first effects of fluid inertia on flows in ordered and random arrays of spheres. J. Fluid Mech. 2001, 448, 213. (6) Hill, R. J.; Koch, D. L.; Ladd, J. C. Moderate-Reynolds-number flows in ordered and random arrays of spheres. J. Fluid Mech. 2001, 448, 243. (7) Benyahia, S.; Syamlal, M.; O’Brien, T. J. Extension of Hill-KochLadd drag correlation over all ranges of Reynolds number and solids volume fraction. Powder Tech. 2006, 162, 166. (8) Beetstra, R.; van der Hoef, M. A.; Kuipers, J. A. M. Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres. AIChE J. 2007, 53, 489. (9) Yin, X. L.; Sundaresan, S. Drag law for bidisperse gas-solid suspensions containing equally sized spheres. Ind. Eng. Chem. Res. 2009, 48, 227. (10) Geldart, D. Types of gas fluidization. Powder Tech. 1973, 7, 285. (11) Arastoopour, H.; Gidaspow, D. Analysis of IGT pneumatic conveying data and fast fluidization using a thermo-hydrodynamic model. Powder Tech. 1979, 22, 77.
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010 (12) McKeen, T.; Pugsley, T. Simulation and experimental validation of a freely bubbling bed of FCC catalyst. Powder Tech. 2003, 129, 139. (13) Agrawal, K.; Loezos, P. N.; Syamlal, M.; Sundaresan, S. The role of meso-scale structures in rapid gas-solid flows. J. Fluid Mech. 2001, 445, 151. (14) Wang, J.; van der Hoef, M. A.; Kuipers, J. A. M. Why the twofluid model fails to predict the bed expansion characteristics of Geldart A particles in gas-fluidized beds: A tentative answer. Chem. Eng. Sci. 2009, 64, 622. (15) Li, J.; Kuipers, J. A. M. Gas-particle interactions in dense gasfluidized beds. Chem. Eng. Sci. 2003, 58, 711. (16) Gidaspow, D.; Huilin, L. Equation of state and radial distribution functions of FCC particles in a CFB. AIChE J. 1998, 44, 279. (17) Ye, M.; van der Hoef, M. A.; Kuipers, J. A. M. From discrete particle model to a continuous model of Geldart A particles. Chem. Eng. Res. Des. 2005, 83, 833. (18) van Wachem, B.; Sasic, S. Derivation, simulation and validation of a cohesive particle flow CFD model. AIChE J. 2008, 54, 9. (19) Yang, N.; Wang, W.; Ge, W.; Wang, L.; Li, J. Simulation of heterogeneous structure in a circulating fluidized-bed riser by combining the two-fluid model with the EMMS approach. Ind. Eng. Chem. Res. 2004, 43, 5548. (20) Igci, Y.; Andrews, A. T.; Sundaresan, S.; Pannala, S.; O’Brien, T. Filtered two-fluid models for fluidized gas-particle suspensions. AIChE J. 2008, 54, 1431. (21) Balzer, G.; Simonin, O.; Boelle, A.; Lavieville, J. A unifying modeling approach for the numerical prediction of dilute and dense gassolid two phase flow. Preprints of the 5th International Conference on Circulating Fluidized Beds, Beijing, China, May 28-31,1996. (22) Benyahia, S.; Syamlal, M.; O’Brien, T. J. Study of the ability of multiphase continuum models to predict core-annulus flow. AIChE J. 2007, 53, 2549. (23) Carnahan, N. F.; Starling, K. E. Equation of state for nonattracting rigid spheres. J. Chem. Phys. 1969, 51, 635. (24) Srivastava, A.; Sundaresan, S. Analysis of a frictional-kinetic model for gas-particle flow. Powder Tech. 2003, 129, 72. (25) Benyahia, S. Validation study of two continuum granular frictional flow theories. Ind. Eng. Chem. Res. 2008, 47, 8926. (26) Passalacqua, A.; Marmo, L. A critical comparison of frictional stress models applied to the simulation of bubbling fluidized beds. Chem. Eng. Sci. 2009, 160, 2795. (27) Johnson, P. C.; Jackson, R. Frictional-collisional constitutive relations for granular materials, with application to plane shearing. J. Fluid Mech. 1987, 176, 67.
5131
(28) Karri, S. B. R.; Knowlton, T. PSRI Challenge Problem 1, Workshop 3 Modeling Test. Presented at the Eighth International Conference on Fluidization. Tour, France, May 14-19, 1995. (29) Chalermsinsuwan, B.; Piumsomboon, B.; Gidaspow, D. Kinetic theory based computation of PSRI riser: Part I-Estimate of mass transfer coefficient. Chem. Eng. Sci. 2009, 64, 1195. (30) Chalermsinsuwan, B.; Piumsomboon, B.; Gidaspow, D. Kinetic theory based computation of PSRI riser: Part II-Computation of mass transfer coefficient with chemical reaction. Chem. Eng. Sci. 2009, 64, 1212. (31) Benyahia, S.; Arastoopour, H.; Knowlton, T. M.; Massah, H. Simulation of particles and gas flow behavior in the riser section of a circulating fluidized bed using the kinetic theory approach for the particulate phase. Powder Tech. 2000, 112, 24. (32) Benyahia, S.; Arastoopour, H.; Knowlton, T. M. Numerical simulation of a large-scale circulating fluidized bed. In Circulating Fluidized Bed Technology VII; Grace, J. R., Zhu J., de Lasa, H., Eds.; Canadian Society for Chemical Engineering: Ottawa, Canada, 2002; p 451. (33) Guenther, C.; Syamlal, M. The effect of numerical diffusion on simulation of isolated bubbles in a gas-solid fluidized bed. Powder Tech. 2001, 116, 142. (34) Khan, A. R.; Richardson, J. F. The resistance to motion of a solid sphere in a fluid. Chem. Eng. Commun. 1987, 62, 135. (35) Benyahia, S. A time-averaged model for gas-solids flow in a onedimensional vertical channel. Chem. Eng. Sci. 2008, 63, 2536. (36) Karri, S. B. R.; Knowlton, T. M. Wall solids upflow and downflow regimes in risers for group A solids. In Circulating Fluidized Bed Technology VII; Grace, J. R., Zhu J., de Lasa, H., Eds.; Canadian Society for Chemical Engineering: Ottawa, Canada, 2002, 310. (37) Igci, Y.; Sundaresan, S. Coarse-grid simulation of fluidized gasparticle flows. Presented at the 2008 AICHE annual meeting, Philadelphia, PA, Nov 16-21, 2008. (38) Andrews, A. T.; Loezos, P. N.; Sundaresan, S. Coarse-grid simulation of gas-particle flows in vertical risers. Ind. Eng. Chem. Res. 2005, 44, 6022.
ReceiVed for reView April 27, 2009 ReVised manuscript receiVed July 7, 2009 Accepted July 8, 2009 IE900658K
