Numerical Simulations of a Bubbling Fluidized Bed Reactor with an

Article pubs.acs.org/IECR

Numerical Simulations of a Bubbling Fluidized Bed Reactor with an Energy Minimization Multiscale Bubble Based Model: Effect of the Mesoscale Musango Lungu, Jingdai Wang,* and Yongrong Yang State Key Laboratory of Chemical Engineering, Department of Chemical and Biological Engineering, Zhejiang University, Hangzhou, Zhejiang 310027, P.R.China S Supporting Information *

ABSTRACT: The present work explores the influence of the mesoscale structures on simulated hydrodynamics including turbulent properties in a bubbling fluidized bed. A bubble based energy minimization multiscale (EMMS) approach has been used to account for heterogeneous structures. Simulations have been performed using the Gidaspow drag model and an EMMS bubble based drag model which accounts for the mesoscale structures via the heterogeneous index. A comparative study between the models in terms of the predicted fluidization behavior, solids volume fraction distribution, bed pressure drop, intermittency index, granular temperatures, and turbulent properties has been done in order to understand the effect of the mesoscale on the predicted quantities. Our study reveals that a multiscale modeling approach is necessary to predict correct bubbling fluidized bed hydrodynamics and transport properties. Simulated particle velocities and granular temperatures agree fairly well with experimental data and empirical correlations from the open literature.

1. INTRODUCTION Gas-particle flows in fluidized beds are chaotic and inherently unstable exhibiting fluctuations spanning a wide range of spatial and temporal scales. This instability (typically associated with the relative motion between the gas and particles) is responsible for the formation of mesoscale structures, i.e., particle clusters, streamers, voids, and bubbles. Further, it has been reported that bubbles in bubbling fluidized beds and clusters in circulating fluidized bed reactors (CFBs) are as a consequence of the same instability.1 With this thought in mind, it seems natural that the key to successful modeling and capturing of these structures is to provide an appropriate drag model which is able to account for heterogeneous structures. McKeen and Pugsley2 demonstrated that the Gibilaro drag model3 in its original form could not capture the bed expansion in a freely bubbling fluidized bed containing FCC particles (Geldart A type). Instead, they introduced a scaling factor using a trial and error method until the scaling factor which gave the best prediction in comparison to the experimental results was obtained; other researchers4,5 have also applied such an ad-hoc procedure. It is not surprising therefore that such homogeneous drag correlations fail to capture typical flow structures in gas− solid fluidized beds in their original form because they were originally developed from experimental data of homogeneous systems like pressure drop data for fixed beds, thereby failing to capture the heterogeneous structures which exist at the scale of the computational cell. Wang et al.6 showed that even in the absence of a scaling factor simulations using traditional drag models in their original form can predict correct hydrodynamics of bubbling fluidized beds employing Geldart A particles so long as sufficiently fine grids and a small time step are used. Even with the current advancements in computational hardware, simulations involving fine grids are costly in computational time and not ideal for industrial scale fluidized beds; © 2014 American Chemical Society

the same is true for direct numerical simulations (DNS) and Lattice-Boltzmann simulations which instead are useful for obtaining accurate correlations used in coarse grid simulations. DNS simulations solve unsteady Navier−Stokes equations on spatial grids that are sufficiently fine to resolve the Kolmogorov Length scales at which energy dissipation takes place and with time steps sufficiently small to revolve the period of the fastest fluctuation. In the light of limited computational resources, practical design and scale-up of industrial scale fluidized beds will still rely on the continuum approach for some time to come. The continuum approach treats the gas phase and solid phase as interpenetrating continua and hence the name twofluid model. Since the solid phase is said to behave as a fluid, expressions are required to describe the solid pressure, solid viscosity, and other properties analogous to gas pressure, gas viscosity, and so on in gases and to this effect the kinetic theory of granular flow (KTGF)7 which is widely used and has been adopted in the present work. Industrial and large scale modeling of fluidized bed reactors using the continuum approach has been done by some research workers, and the results have been quite promising. In these studies, however, coarse grids were employed and to that effect modified drag correlations were used so as to preserve the heterogeneous structures which would otherwise be lost as a result of using coarse grids. These mesoscale structures have been said to affect the overall performance of fluidized beds vis-à-vis heat transfer, mass transfer, and reactions (ref 8 and references cited therein) and in bubbling fluidized beds they may relate to bubbles. Ge and co-workers9 give a comprehensive list of Received: Revised: Accepted: Published: 16204

April 10, 2014 August 25, 2014 September 29, 2014 September 29, 2014 dx.doi.org/10.1021/ie501492y | Ind. Eng. Chem. Res. 2014, 53, 16204−16221

Industrial & Engineering Chemistry Research

Article

bubbling to fast fluidization; moreover, bubbles and/or voids are easier to describe in comparison to clusters. A bubbleemulsion hydrodynamic model for modeling bubbling fluidized bed reactors recently advanced by Gao and co-workers27 also reinforces the idea of using bubble based models for the successful modeling of bubbling fluidized bed reactors. Coarse grid simulations performed by these authors agreed quite well with experimental results from the literature. For the sake of a clear and concise presentation, the bubble based EMMS drag model is abbreviated as B while the Gidaspow drag model is abbreviated as G. Thus, hereafter, abbreviations G and B will be used. The aim or goal of the work is to get insight on the role of the mesoscales on the predicted hydrodynamics in the fluidized bed including granular temperatures and turbulent properties. The origin, evolution, and influence of the mesoscales are far from being well understood,28 and it is hoped that our work will advance the understanding of mesoscale modeling. The remainder of the paper is organized as follows: in Section 2, the CFD model used is presented. Section 3 describes the simulation setup. Section 4 combines the results and discussion. Finally, the conclusions drawn from this work are given in Section 5.

different multiphase industrial scale reactors modeled using such an approach. Different approaches have been reported in open literature on the modeling mesoscale structures in gas−solid systems. One such approach is the use of the so-called filtered models developed to account for the subgrid closures of the fluidparticle drag force.10−13 In the filtering approach, results obtained from highly resolved simulations of model flow problems, each representing possible flow conditions occurring in the industrial scale reactor, are filtered (averaged) to construct subgrid corrections.14 However, another approach is the energy minimization multiscale (EMMS) method originally developed to model particle clustering in high velocity CFBs. The EMMS model has been successfully used to model turbulent fluidized beds15 and CFB risers.16,17 In the recent past, the concept of the EMMS has been applied in the development of bubble based drag models applicable to the bubbling regime of fluidization.18−20 From the literature surveyed, it is clear that most of the modeling effort using the EMMS based drag models has been devoted to circulating fluidized bed risers and downers. Scale up and design of bubbling fluidized beds still possess a major challenge partly due to the lack of proper understanding of the mesoscale, and in view of this challenge, EMMS based modeling seems to be a promising tool in understanding and tackling this problem.9 Schneiderbauer et al.21 have compared different subgrid drag modifications in bubbling fluidized beds and have come to a conclusion that, though these modifications differ in terms of their dependencies on void fraction and slip velocity, they are all able to resolve the main features of the gas and particle flow. However, the focus of this work is the EMMS and its application to bubbling fluidized beds due to its simplicity compared to other subgrid modeling approaches. Loha et al.22 investigated the effect of an EMMS cluster based model among other drag models on the hydrodynamics of a 2D bubbling fluidized bed reactor. In comparison to other models that were tested, they observed that the EMMS based model could not capture the core-annulus structure of the bed but predicted the lateral distribution of the granular temperature reasonably well. Similarly, Wang and Liu23 applied an EMMS based model modified by the introduction of an implicit cluster diameter description in the modeling and simulation of a bubbling fluidized bed of FCC particles. The work demonstrated that the model was able to predict the experimentally observed axial and radial solid concentration profiles as well as the radial particle velocities quite well. In the current work, a two-fluid model (TFM) with closures from the kinetic theory of granular flows (KTGF) has been used to simulate a bubbling fluidized bed reactor. Our work presents an effort in the understanding of the role of mesoscale structures on hydrodynamics in bubbling fluidized bed reactors including turbulent properties such as Reynolds stresses, granular temperatures (laminar and turbulent), and energy spectra. A bubble based EMMS drag model and the Gidaspow drag model have been used to close the interphase momentum exchange coefficient. The simulated results are compared with experimental results from the literature24,25 and empirical correlations. The bubble based EMMS model18 is a derivative of the original EMMS model26 used for mesoscale modeling of clusters in high velocity circulating fluidized beds. In the bubble based EMMS model, the bubbles replace the clusters as the characteristic structures thereby allowing the EMMS model to be applied to a wider range of flow regimes ranging from

2. CFD MODEL The continuum governing equations, constitutive relations, and boundary conditions for a monodisperse system are presented below: 2.1. Governing Equations. Mass conservation for phase k (k = g for gas and s for solid phase) ∂ (αkρk ) + ∇ × (αkρk uk) = 0 ∂t

(1)

where n

∑ αk = 1 k=1

(2)

Momentum conservation for the solid phase: ∂ (αsρs us) + ∇ × (αsρs usus) = ∇ × τs − αs∇p − ∇ps ∂t + βgs(u g − us) + αsρs g

(3)

Momentum equation for the gas phase: ∂ (αgρg u g) + ∇ × (αgρg u gu g) = ∇ × τg − αg∇p ∂t − βgs(u g − us) + αgρg g

(4)

Granular temperature equation: ⎤ 3⎡ ∂ ⎢⎣ (αsρs Θs) + ∇ × (αsρs usΘs)⎥⎦ = ( −ps I + τs) 2 ∂t : ∇us + ∇ × (κs∇Θs) − γs − 3βgsΘs

(5)

2.2. Constitutive Equations for Interphase Momentum Transfer. Drag Model G.7 ⎧ α (1 − αg)μg αsρg |u g − u s| ⎪150 s , αg ≤ 0.8 + 1.75 2 ⎪ ds αgds ⎪ βgs = ⎨ ⎪ 3 αsαgρg |u g − u s| ⎪ CD αg −2.65 , αg > 0.8 ⎪4 ds ⎩ 16205

(6)

dx.doi.org/10.1021/ie501492y | Ind. Eng. Chem. Res. 2014, 53, 16204−16221

Industrial & Engineering Chemistry Research ⎧ 24 0.687 ⎪ [1 + 0.15(Res) ], Res < 1000 C D = ⎨ Res ⎪ ⎩ 0.44, Res > 1000 Res =

Article

Supporting Information and the bubble diameter, db, from equation TS1-6 in the Supporting Information. 7. Compute the mass specific energy consumption for suspending the particles, Ns, from equation TS1-10 in the Supporting Information and store other parameters if Ns is the minimum. 8. Compute βEMMS from eq 11 and Hd from eq 10. The solution of the bubble based EMMS model yields a two-dimensional matrix of the voidage and the heterogeneous index among other parameters. The resulting two-dimensional matrix is regressed using an exponential fitting function of the form HD = a + beαg/c where a, b, and c are constants in the Origin-Lab graphing software with a correlation coefficient of over 0.99 achieved for all the cases. Several such fitting functions were obtained at various superficial gas velocities based on the literature data24 and have been presented in Table TS-2 in the Supporting Information. Figure 1 is the plot of the fitted

(7)

αgρg |u g − us|ds μg

(8)

18

Drag Model B.

⎧ α (1 − αg)μg αsρg |ug − us| ⎪150 s + 1.75 , αg ≤ αmf 2 ⎪ ds αgds ⎪ ⎪ αsαgρg |ug − us| βgs = ⎨ 3 C D αg −2.65Hd , αmf ≤ αg < αd ⎪4 ds ⎪ ⎪ 3 αsαgρg |ug − us| −2.65 αg , αg ≥ αd ⎪ CD ds ⎩4

(9)

Hd is the heterogeneous index which is a dimensionless drag coefficient scaled with the Wen and Yu drag coefficient and represented as Hd = βEMMS /βWenandYu

(10)

where βEMMS =

αg 2 Uslip

[(1 − δ b)(1 − αe)(ρs − ρg )(g + ae)

+ δ b(ρe − ρg )(g + ab)] βWenandYu =

3 ρg αgαs|u g − us| −2.65 CD αg 4 ds

(11)

(12)

The Wen and Yu drag model can be viewed as a constant factor without considering the heterogeneous structure. Though the Wen and Yu drag model has been used and accepted as the traditional model to which EMMS based models are scaled, it has now been recently shown29 that the EMMS based models are insensitive to the scaling traditional drag model used. However, in the present work, the Wen and Yu drag model has been maintained. Equation 11 is the structure dependent drag coefficient, and the structure parameters in the equation are obtained by solving a set of nonlinear equations and an objective function given in Table S1 in the Supporting Information using a matlab code. A global search scheme is used to calculate the parameters and is as follows: 1. Input physical parameters ρg, ρs, μg, and superficial gas velocity, ug, and traverse over the trial value of αe within the range of [αmf, αg]. 2. Compute volume fraction of bubbles, δb, from equation TS1-1 in the Supporting Information. 3. Compute superficial slip velocity in the emulsion phase, Use, from equation TS1-5 in the Supporting Information. 4. Compute the inertial term in the bubble phase, ab, from equation TS1-7 in the Supporting Information. 5. Compute superficial particle velocity in the emulsion phase, Upe, and gas superficial gas velocity in the emulsion phase, Uge, from equations TS1-3 and TS1-4, respectively, in the Supporting Information. 6. Compute the superficial bubble velocity relative to bubble phase, Ub, from equation TS1-2 in the

Figure 1. Plot of voidage versus fitted heterogeneous index functions at different superficial gas velocities.

heterogeneous function at different superficial gas velocities against voidage. The following can be deduced from the figure: at low fluidization velocities, for instance at a superficial gas velocity of 0.2 m/s, the Hd is close to unity within the entire range indicating that the flow structure is more or less homogeneous, but with increasing superficial gas velocity, it reduces further and further as the flow becomes increasingly inhomogeneous. This explains why the bubble based EMMS model predicts reduced drag coefficients in comparison to the standard model. Thus, HD is a correction factor applied to the Wen and Yu drag correlation within a certain voidage range that accounts for the heterogeneity in the bed. Constitutive Equations Based on the Kinetic Theory of Granular Flow. Solid phase stress tensor: ⎛ 2 ⎞ τs = αsμs [∇us + (∇us)T ] + αs⎜λs − μs ⎟∇ × usI ⎝ 3 ⎠

(13)

Gas phase stress tensor: τg = αgμg [∇u g + (∇u g)T ] − 16206

2 αgμ (∇ × u g)I 3 g

(14)



Article

Solids pressure: ps = αsρs [1 + 2(1 + e)αsg0]Θs

Table 1. Summary of Numerical Parameters Used in the Simulations

(15)

λs =

Θs 4 2 αs ρs dsg0(1 + e) 3 π

air density air viscosity solids density particle diameter void fraction at minimum fluidization minimum fluidization velocity superficial gas velocity terminal velocity30 wall boundary conditions

(16)

Shear viscosity of solids: μs = μs,kin + μs,col + μs,fr

(17)

Kinetic viscosity, μs,kin: μs,kin =

10ρs ds π Θs ⎡ ⎤2 4 1 (1 e) g + + α α ⎢ 0 s⎥ ⎦ s 5 96αs(1 + e)g0 ⎣

(18)

particle−particle restitution coefficient particle−wall restitution coefficient packing limit

Collisional viscosity, μs,col: μs,col =

Θs 4 αsρ dsg (1 + e) αs 5 s 0 π

(20)

where ϕ is the internal angle of friction and I2D is the second invariant of the deviatoric stress tensor. Diffusion of granular temperature: κs =

150ρs ds Θsπ ⎡ ⎤2 Θs 6 2 ⎢1 + αsg0(1 + e)⎦⎥ + 2ρs dsαs g0(1 + e) ⎣ 384(1 + e)g0 5 π

(21)

Collisional dissipation of solid particle fluctuating energy: γs =

12(1 − e 2)g0 ds π

2 3/2 ρα Θ s s s

(22)

Radial distribution function: g0 =

1 + 2.5αs + 4.59αs 2 + 4.52αs 3 ⎡ ⎢1 − ⎣

3 ⎤0.678

( ) ⎥⎦ αs

αs,max

Jung et al.25

1.21 kg/m 1.81 × 10−5 kg/m s 2500 kg/m3 310 μm 0.46

1.206 kg/m3 1.82 × 10−5 kg/m s 2500 kg/m3 530 μm 0.4

0.15 m/s 0.4 m/s 2.31 m/s no slip for both gas and solid phases 0.9

0.229 m/s 0.587 m/s 3.99 m/s no slip for air and specularity coefficient of 0.6 for solid phase 0.99

0.9

0.95

0.55

0.63

Detailed numerical parameters for both cases have been given in Table 1. The CFD commercial software package Fluent 6.3.26 was used as the platform to execute the simulations. For the first case, an initial grid sensitivity study showed that a grid size of 32 by 320 was sufficient to give a grid independent numerical solution. A detailed grid study is presented later on in the paper. The QUICK scheme is selected to discretize the convective terms while the second order implicit scheme is adopted for the transient formulation. For the pressure−velocity coupling, the Phase Couple SIMPLE algorithm is used. Under relaxation factors of 0.5 for pressure, 0.4 for momentum, 0.2 for volume fraction, and 0.2 for granular temperature were used. The simulations are run for 20 s in an unsteady mode with a constant time step of 1 × 10−4 s, and the numerical solution is obtained within a residual tolerance of 1 × 10−4 for all equations. For the purpose of obtaining time averaged results, quantitates are averaged for the last 15 s as it takes nearly 5 s for the initial startup effects to subside. For the second case, a grid size of 5 mm in both directions was used following that in ref 31 in which the same case was simulated. The simulations were run for 40 s, and time averaged quantities were obtained between 15 and 40 s. Similar settings for the discretization scheme, velocity coupling, under relaxation factors, time step, and residual tolerance were adopted for the second case.

ps sin ϕ 2 I2D

3

(19)

Frictional viscosity, μs,fr: μs,fr =

Yuu et al.24

description

Bulk solid viscosity:

(23)

3. SIMULATION SETUP The simulated system is based on the experimental setup of Yuu et al.24 The rectangular fluidized bed has a height of 0.806 m and a width of 0.0806 m, and the initial static bed height is 0.104 m. Geldart B particles were employed with a Sauter mean diameter of 310 μm and density of 2500 kg/m3. Additionally, a second system based on the experimental work of Jung et al.25 was simulated to validate the present model and demonstrate that it can predict the granular temperature with reasonable accuracy and consequently the turbulent properties. Their experimental setup consisted of a two-dimensional rectangular bed with a height of 0.58 m and width of 0.154 m. Similarly, glass particles (Geldart B) were employed with air as the fluidizing medium. The initial static bed height was 0.14 m. Granular temperatures of glass bead particles were measured indirectly from instantaneous particle velocities using a CCD camera technique. Nonetheless, most of the results reported in the paper are based on the first case.

4. RESULTS AND DISCUSSION The results in this work have been mostly obtained at two dimensionless heights, i.e., y/D = 0.375 near the distributor plate and y/D = 1.25 corresponding to the heights used in the experimental setup. Axial profile results along the height of the reactor were spatially (lateral direction) and time averaged (5−20 s). In addition, transient local axial and particle velocities as well as the particle volume fraction were also monitored at 10 lateral positions at these two heights also corresponding to the positions used in the experimental setup. In their work, Yuu et al.24 obtained instantaneous particle velocities by recording colored particle trajectories using a video camera allowing them to obtain lateral profiles of time averaged axial and lateral particle velocities. Their results have been compared with our simulation results. For the second case of Jung et al.,25 only the laminar granular temperature results have been reported for the 16207



Article

pressure drop; then, the time averaged pressure drop over a given time period can be expressed mathematically as

sole purpose of validating that the CFD model is capable of predicting granular temperature and therefore turbulent properties. 4.1. Grid Study. The sensitivity of the numerical solution to the grid size is one of the first steps in carrying out a successful numerical simulation as it quantifies the errors due to spatial discretization. Ideally, the spatial discretization error should be kept to a minimum. In the present work, three grid resolutions were tested for drag model G, 16 × 160 (coarse grid), 32 × 320 (medium grid), and 64 × 520 (fine grid) corresponding to a refinement factor of 2 in both the lateral and vertical directions; meanwhile, for drag model B, only the coarse and medium grids were used due to the fact that drag models employing subcorrections such as EMMS based models have demonstrated their capability to capture the detailed flow structure with relatively coarse grids. The results from the grid refinement study are reported in terms of average macroscopic fluidization properties in the bed, namely, the height, voidage, pressure drop, and pressure drop gradient and are presented in Table 2.

S ̅ (t , Δt ) =

H (m)

αg,ave

ΔP (Pa)

ΔP/H (kPa/m)

coarse grid (G) medium grid (G) fine grid (G) coarse grid (B) medium grid (B)

0.1467 0.1413 0.1447 0.1326 0.1273

0.6176 0.6102 0.6382 0.5777 0.5591

1372.80 1360.31 1328.98 1376.17 1352.63

9.36 9.63 9.18 10.38 10.63

t +Δt

S dt /Δt

t≥5

(24)

where t is the initial point for averaging and Δt is the time interval for averaging. The mean and standard deviation of the time average pressure drop can be calculated by varying the starting point t over the entire data range and calculated using the following expressions, respectively: ⟨S ̅ (Δt )⟩ =

1 Nt

Nt

∑ S ̅(ti , Δt )

(25)

i=1

and ⟨S ̅(Δt )⟩std

⎡ 1 Nt ⎤0.5 2⎥ ⎢ = ∑ (S ̅(ti , Δt ) − ⟨S ̅(Δt )⟩) ⎢⎣ Nt i = 1 ⎥⎦

(26)

where Nt is the number of data samples. The ratio ⟨S̅(Δt)⟩std/ ⟨S(̅ Δt)⟩ gives the time averaging error as a function of Δt. The results of this analysis (see Table 3) reveal that the difference in

Table 2. Averaged Macroscopic Fluidization Properties for Drag Models G and B at Different Grid Resolutions grid size

∫t

Table 3. Time Average Pressure Drop Statistics at Different Averaging Time Intervals averaging time interval, Δt (s) 5 time-average pressure drop

⟨S̅(Δt)⟩ (Pa) ⟨S̅(Δt)⟩std (Pa) relative error (%)

It can be observed that generally model B predicts relatively lower average bed height and voidage for both the grid sizes in comparison to drag model G as expected as it accounts for the bed heterogeneity. The average pressure drop decreases with mesh refinement for both drag models probably due to better resolution of the near wall region. The average pressure drop gradient in the bed is a good indicator of the average bed solid volume fraction. An increase in the average pressure drop gradient signifies an increase in the solids volume concentration. However, the average pressure drop gradient predicted by drag model B for both grid sizes is larger in comparison to drag model G for all grid sizes demonstrating the effect of heterogeneous structures. Interestingly, the bed height and average bed voidage increases when a fine grid is used for drag model G in comparison to the medium grid. Cloete and coworkers32 have previously observed an increase in the expanded bed height with increasing mesh refinement for Geldart B particles, and the reason for this behavior is not particularly well-known. For the respective drag models, the averaged macroscopic fluidization properties do not show large variations over the different mesh resolutions used and, thus, the medium grid was adopted in the simulations as a compromise between computational time and accuracy. As aforementioned in Section 3, the simulations were run for a total of 20 s and by 5 s the system had reached a steady state by monitoring the transient pressure drop at the reactor inlet (see Figure 6a). Following Gel et al.,33 the uncertainty caused by time averaging has been quantified. The pressure drop is averaged at different time intervals, and the results are compared taking into account statistical analyses, namely, the standard deviation and the relative error. The transient pressure drop predicted by drag model B has been used arbitrarily. Let S be the quantity of interest (Qoi) which in this case is the

10

15

1351.9238 1351.8686 1352.6265 2.1274 0.1574

0.8073 0.0597

0 0

the time averaged pressure drop values at different averaging time intervals is negligible thereby showing that the system has reached a pseudosteady state by 5 s, and furthermore, the relative error for the three averaging time intervals does not vary appreciably and is less than 1 % for all the cases. For the averaging period of 15 s, only 1 data set is available thus giving a standard deviation and relative error of 0 however it can be observed that the relative error decreases with increasing averaging period. 4.2. Fluidization Behavior. Figure 2 displays the instantaneous and time average gas volume fraction contour plots for models G and B at a superficial gas velocity of 0.4 m/s. Both models capture bubble growth with height with small bubbles forming near the distributor plate and growing as they rise due to coalescence. However, the relative large size of bubbles and large average voidage predicted by drag model G in comparison to drag model B is easily observable which is consistent with the findings from the grid size study. Large bubble size predicted as a consequence of using drag model G has been previously reported by Bokkers and co-workers34 who observed drag model G over predicted experimental values. Better agreement with experimental results was obtained when these authors used a lattice-Boltzmann derived drag model which accounted for the bed heterogeneity. As was aforementioned in the introductory section, the mesoscale structure in bubbling fluidized beds may relate to bubbles and by controlling or manipulating this property the performance of the fluidized bed can be altered and all together be enhanced or optimized. As we shall demonstrate in the remainder of this paper, the mesoscale structures significantly alter the hydrodynamics in the fluidized bed. Second, the simulated bed height predicted using model G is greater probably due to the presence of large voids in the bed. 16208



Article

significantly lower than that predicted using model G. This result confirms our earlier suspicion from the observations of the gas volume fraction contour plots. The bed expansion frequency is calculated using the FFT technique with a sampling frequency of 100 Hz for data between 5 and 20 s to reduce the influence of startup transients. Interestingly, model G predicts a higher average expanded bed height and bed expansion ratio relative to model B, and yet, the expansion frequency is lower. This can be explained by appreciating the fact that the fluctuation period is related to the ratio of the expanded bed height to the rise velocity.36 From the numerical simulations, it is clear that the choice of using model B results in a compact bed as reflected by the bed expansion ratio. A compact bed presents a smaller voidage through which the gas can pass thereby increasing the velocity of the gas. An increased velocity entails a shorter residence time and consequently a higher bed expansion frequency. Therefore, it is reasonable that model B gives a higher bed expansion frequency. The bed height fluctuations are closely associated with large scale oscillations due to bubbles or clusters and can be interpreted as gravity waves.37 Gidaspow and co-workers37 obtained an analytical solution of the dominant frequency from the basic equations of motion. The basic frequency is the gravity divided by the initial bed height, hence the name gravity wave. It is corrected for by approximately the square root of the solids volume fraction and is given as

Figure 2. (a) Instantaneous (20 s) and (b) time averaged (5−20 s) void fraction contour plots for drag models G and B, respectively.

The dynamic bed height is given by 2⟨h⟩ where ⟨h⟩, the average bed height, is given by Goldschmidt and co-workers35 as ⟨h⟩ =

∑# cells αs,cellhcell ∑# cells αs,cell

(27)

1/2 1/2 1 ⎛ g ⎞ ⎡ (3αs/αg + 2)αs ⎤ ⎥ f= ⎜ ⎟ ⎢ ⎥⎦ αs,max 2π ⎝ Ho ⎠ ⎢⎣

where αs,cell is the solid volume fraction in each computational cell and hcell is the height of each computational cell, respectively. Figure 3 displays the simulated transient expanded bed

(28)

where g is the gravitational acceleration and H0 is the initial bed height. Plugging in the relevant values in the formula gives a fundamental frequency of 2.58 Hz using model G and 2.89 Hz using model B meanwhile the expansion frequencies using FFT technique are 2.25 and 2.4 Hz for drag models G and B, respectively. It can be seen that the bed expansion frequencies obtained from the simulations using the FFT technique and the fundamental frequencies from the analytical expression are in close agreement though the values obtained using the analytical expressions are higher for both models. This result therefore supports the idea of large scale oscillations of bubbles or particle clusters being responsible for the observed bed height fluctuations. It can also be seen that accounting for mesoscale structures increases the dominant frequency in the system due to the reduced voidage. These frequencies calculated are typical25 of the bubbling regime. A summary of the bed dynamics results is given in TableS3 in the Supporting Information. 4.3. Solid Volume Fraction Distribution. Figure 4 shows the lateral profiles of the time averaged solids volume fraction predicted by the two models at dimensionless heights of y/D = 0.375 and y/D = 1.25. The effect of model B on the predicted solids volume fraction is evident. The model predicts a higher time averaged solids volume fraction especially at the lower height. At the higher height, the difference is not that much, and in fact between x/D values of 0.13 and 0.34, model G gives greater values. Such a solids volume fraction profile resulting from using an EMMS based model has been observed elsewhere in a recent study.38 For spherical particles, the terminal velocity, ut, can be expressed as

Figure 3. Simulated expanded bed height for models G and B.

height obtained using the two drag models. Qualitatively, both models exhibit a similar behavior with the bed height rising sharply at the start of the fluidization process due to the initial startup rising bubbles carrying along with them the bed material, and with the bubbles reaching the bed surface and bursting, the solid particles are released falling back into the bulk of the bed under the action of gravity only to be brought up to the surface again by further rising bubbles. This motion brings about the observed bed height fluctuations. Quantitatively, the dynamic bed height predicted by model B is

⎛ 4 ρs g ⎞1/2 ut = ⎜ ⎟ ⎝ 3 CD ⎠ 16209

(29)



Article

Figure 4. Lateral profiles of time averaged solids volume fraction for models G and B at (a) y/D = 0.375 and (b) y/D = 1.25.

4.4. Pressure Drop. Figure 6a displays the transient pressure drop for drag models G and B. Both models correctly predict the rapid drop in pressure at the start of the simulation due to the unlocking of particles in the packed bed, and as the time elapses, both models exhibit the fluctuations due to the passage of bubbles about the classical equation, i.e., eq 30; by 5 s, the fluidized bed reaches a statistical steady state for both models as the pressure drop does not significantly vary with time, although a closer inspection of the plot reveals that the fluctuations from drag model G are somewhat larger reflecting the passage of the larger bubbles. The predicted spatially and time averaged bed pressure drop profile along the dimensionless height of the reactor for the two models is displayed in Figure 6b. The pressure drop profile reflects both the macro- and mesoscale features of how the solids are distributed along the axial direction.39 At the reactor inlet, model G gives a pressure drop of only 1233.14 Pa while model B predicts a pressure drop of 1417.47 Pa. The classical expression for predicting the pressure drop across the bed is given as

Equation 29 shows that the terminal velocity is inversely proportional to the square root of the drag coefficient, CD; therefore, lower CD values predicted by EMMS based drag models will give higher terminal velocities for a constant particle diameter meaning that particles are therefore not easily carried into the dilute phase or entrained all together with, and this brings about the observed high solids concentration in the dense phase of the reactor predicted by model B. A plot of a spatially (along bed width) and time averaged (5−20 s) profile of the solid volume fraction with dimensionless bed height is shown in Figure 5. The higher expanded bed height predicted

ΔP = (1 − αmf )(ρs − ρg )g Hmf

(30)

It can be seen that the pressure drop profile is a function of the solid volume fraction when other quantities are constant. The pressure drop profile is a reflection of the solid volume fraction as aforementioned with model B predicting a higher pressure drop between the y/D range of 0 and 0.76, but beyond y/D = 0.76, the two pressure profiles now cross each other with model G, predicting higher pressure drop values. Yuu et al.24 did not obtain any experimental data for the pressure profile in the bed, but they did however obtain plots of pressure drop versus superficial gas velocity from both experiments and DEM simulations. At a velocity of 0.4 m/s, the pressure drop predicted from the DEM simulations is 1433 Pa which is quite close to the value predicted by model B. The findings from this study are in agreement with those of Benyahia40 who observed that the ultimate effect of subgrid drag models in the dense phase is an increase in the solids hold up and consequently the pressure drop.

Figure 5. Axial profile of a spatially (along the bed width) and time averaged (5−20 s) solids volume fraction along the height of the reactor.

by model G is once again obvious. In the dense part of the bed, specifically in the solid volume fraction range of 0.3−0.55, the concentration of particles predicted by model B is higher which results in the compact dense bed height. The solid concentration profiles also have an effect on the bed pressure drop profile as we shall later see. 16210



Article

Figure 6. (a) Variation of pressure drop with time and (b) spatially and time averaged pressure drop profile along the height of the reactor for drag models G and B.

Figure 7. Simulated (a) axial and (b) lateral mean velocities at y/D = 0.375.

4.5. Axial and Lateral Velocity Distributions. Figures 7 and 8 display the simulated axial and lateral distributions of the mean air and particle velocities predicted by models G and B together with the experimental particle velocity data of Yuu et al.24 Both drag models under predict the axial and lateral mean particle velocities, although qualitatively model B predictions give better trends with the experimental results for the most part. At the lower height, i.e., y/D = 0.375, the under prediction is quite significant and this is probably due to the exclusion of the distributor plate in the simulations. In addition, the kinetic theory of granular flow in its original form also has limitations as was demonstrated by Lijie et al.41 Hydrodynamics near the reactor inlet are seriously affected by the distributor plate. However, at the higher height, the predictions from both models are significantly improved especially the prediction of the particle axial velocity and this is due to the fact that at this height the flow is fully developed. As expected, the axial gas velocities predicted by both the models are larger

than the particle velocities, and the difference between the gas and particle velocity is termed as the slip velocity and is the driving force for momentum exchange, as in the same way concentration difference drives mass transport and temperature difference drives energy transport. Model B predicts larger slip axial velocities at both heights as opposed to model G as a result of the reduced bed voidage which increases the gas velocity and consequently the slip velocity. Larger slip velocities using an EMMS based model have been previously reported for a circulating fluidized bed reactor.40 CFBs usually employ Geldart A particles, and therefore, the higher slip velocities observed using drag models incorporating subgrid drag closures are due to gas bypassing of particle clusters, thus increasing the slip between the two phases. In bubbling fluidized beds employing Geldart B particles like in the present work, it is unlikely that particle clustering is significant, and therefore, the observed higher slip velocities using the EMMS based model is due to the reduced voidage as aforementioned. A summary of 16211



Article

Figure 8. Simulated (a) axial and (b) lateral mean velocities at y/D = 1.25.

Figure 9. Axial particle velocity versus solids volume fraction for (a) model G and (b) model B.

distinguish the flow structure between a circulating-turbulent fluidized bed (C-TFB), turbulent fluidized bed (TFB), circulating fluidized bed (CFB), and high density circulating fluidized bed (HDCFB) for Geldart A particles by simultaneously measuring the instantaneous particle velocities and the particle solid volume fraction at different spatial positions. These authors showed that collisions dominated the motion of particles in the C-TFB and TFB as opposed to the CFB and the HDCFB in which the particle motion is controlled by drag forces. This method has for the first time been used in the present work to distinguish the flow structure predicted by the two models in the bubbling regime using CFD. Figures 9 and 10

the axial slip velocity can be viewed from Table S4 in the Supporting Information. In our study, the lateral slip velocity is not as pronounced as the axial slip velocity probably due to the fact that the flow is predominantly in the axial direction. 4.6. Coupling of Instantaneous Solids Concentration and Particle Velocity. Plots of local instantaneous particle velocities against local instantaneous solids volume fraction are very useful in providing an insight into the dynamic flow structure of the gas−solid system under consideration as well as identifying the dominant particle interaction mechanism taking place, i.e., particle−particle interactions (collisions) or gas− particle interactions (drag). Qi et al.42 used this methodology to 16212



Article

Figure 10. Lateral particle velocity versus solids volume fraction for (a) model G and (b) model B.

consistent time averaged solids concentration at the same position given as

show plots of local instantaneous axial and lateral particle velocities versus local instantaneous solid volume fraction at three different positions along the width of the reactor at a dimensionless height of 1.25 for models G and B, respectively. For both models, it is observed that there is no strong correlation between the particle velocities and the solids volume fraction, an indication that gas−particle interactions have a weak effect on the particle motion which is typical and expected in dense gas solid systems. In fact, it is now well established in the modeling of dense gas solid systems that gas turbulence has little or no effect on the modeling results and can be neglected while particle−particle interactions (adjusted via the coefficient of restitution in CFD simulations) have a profound effect on the modeling results. Both models G and B demonstrate that collisions are the dominant mechanism of solids motion as manifested by the wide range distributions of solids concentration and particle velocities implying a random flow structure. However, a closer inspection of the plots reveals some interesting differences. First, the distribution of the axial particle velocities for both models G and B show a wider variation in comparison to the lateral particle velocities showing that the motion in the bed is predominantly in the vertical direction. Second, the particle velocities for model G exhibit a wider scatter in comparison to model B, an indicator that the particle motion predicted by model G is more energetic and a fact we will later demonstrate when we analyze the granular temperature and turbulent properties in the subsequent sections. 4.7. Intermittency Index. The intermittency index43 is defined as γ=

σ σs

σs =

αs(αs,mf − αs)

(32)

It is a parameter that describes the segregations of the gas phase from the solids phase. This parameter assumes values between 0 and 1 with 0 assigned to ideal core-annulus flow and 1 to ideal cluster flow, and in practical cases, all flows vary between these two limits. This criterion is somewhat similar to the compromise between the solid and gas phase in the EMMS formation. According to the EMMS concept at fixed bed operation, the particles dominate the fluid in the sense that the fluid cannot suspend the particles which seek a condition of minimum potential energy; meanwhile, for the dilute transport operation, the fluid completely dominates the particles inhibiting the particle movement tendency. However, fluidized beds are typically operated between these two cases implying that neither condition dominates the flow and instead the actual flow is between these two extremes; that is, the particles and the fluid compromise with each other. Thus, the compromise leads to the formation of heterogeneous structures. Du et al.43 obtained experimental intermittency indices for various fluidization regimes, and they observed that the bubbling regime gave the highest intermittency index at the center of the bed and that it decreased monotonically along the radial distance moving from the bed center to the wall. In the current study, we have compared the radial profile of the intermittency index at two different dimensionless bed heights in a bubbling fluidized bed predicted by models G and B using CFD. Figure 11 displays the intermittency index for drag models G and B at the two heights. The intermittency index for both models is between 0.6 and 0.1 showing that the flow in the bubbling fluidized reactor is neither completely core-annular nor completely clustering but in between. It is also quite evident from Figure 11 that model G predicts a higher intermittent index

(31)

where σ is the standard deviation of the solids concentration fluctuations at a given point in the bed and σs is the standard deviation of the solids concentration fluctuations with a 16213



Article

Figure 11. Simulated lateral profiles of the intermittency index for models G and B at (a) y/D = 0.375 and (b) y/D = 1.25.

also adopted in the present work. The normal Reynolds stress is mathematically defined as

at both heights in comparison to model B. However, the difference in the intermittency index at a dimensionless height of y/D = 1.25 between the two models is less compared with that at a height of y/D = 0.375. This result further validates our earlier postulate that model G predicts larger bubble sizes compared to model B. A higher intermittency index predicted by model G is a reflection of the higher bubble induced particle motion by that model. The high bubble induced motion also explains the previously observed wide scatter of particle velocities using model G in Figures 9 and 10. Another point worth mentioning is that the intermittency index at the lower bed height, i.e., y/D = 0.375, is not uniform as the flow at this height is not fully developed. However, at a height of y/D = 1.25, the dome shaped lateral profile of the intermittency index in most of the core of the bed can be readily observed. At this dimensionless bed height, model G gives a maximum intermittency index of 0.80 at x/D = 0.36 while for model B it is 0.78 and occurs at x/D = 0.54. Further, for model B, the intermittency index is generally reduced with lateral distance moving from the bed center toward the wall except for a peak which appears at x/D = 0.17. Thus, model B predicts an intermittency index lateral profile similar to the one observed by Du et al.43 for the bubbling regime. For model G, the intermittency index slightly increases near the wall which perhaps reveals that the model predicts a greater accumulation of particles near the wall region. 4.8. Reynolds Stresses. Reynolds stresses are the extra stress terms that arise from the derivation of Reynolds decomposed and time averaged forms of the Navier-stokes equations and are characteristic of turbulent motion. The fluctuations of bubbles or clusters give rise to Reynolds stresses. CFD codes such as the IIT code, Fluent, and MFIX have demonstrated their capabilities in predicting Reynolds stresses (both normal and shear) with reasonable accuracy. Turbulent properties such as turbulent granular temperature and turbulent dispersion coefficients can in turn be computed from the normal Reynolds stresses. Jiradilok and co-workers15 have presented a procedure for computing the Reynolds stresses from CFD simulations. Their method has been used previously by various workers44−46 for different fluidization regimes and is

ui′ui′ =

1 m

m

∑ (uik(x , t ) − ui(x))(uik(x , t ) − ui(x)) k=1

(33)

where ui(x) is the time averaged velocity of the particle and is expressed as ui′ =

1 m

m

∑ uik(x , t ) k=1

(34)

The subscript i represents the x or y directions and m is the total number of data over a given time period. In the present work, we considered the time period from 5 to 20 s which is sufficient for the calculation and estimation of turbulent properties. There is no direct way of computing the Reynolds stresses from simulations; instead, axial and radial particle velocities are monitored at 10 different lateral positions at dimensionless bed heights of 0.375 and 1.25, thereby giving a time series of the so-called hydrodynamic velocities at each position. From the time series of the hydrodynamic velocity, other quantities, i.e., the mean particle velocity, fluctuating velocity, and normal Reynolds stresses, are computed. Figures 12 and 13 show axial and lateral normal Reynolds stresses profiles along the width of the reactor at the two heights predicted by the two drag models. From both figures, we yet again observe that model B predicts relatively lower axial and lateral normal Reynolds stresses at all the positions and heights in comparison with model G. A look back at the plots of the instantaneous particle velocities versus instantaneous solid volume fraction, i.e., Figures 10 and 11, provides a clue to this observation. As was earlier demonstrated, the particle velocities predicted by model G exhibit a wider scatter as opposed to model B; naturally, the normal Reynolds stresses predicted by model G, which are primarily computed from the particle velocities, will be larger than those of model B. The axial and lateral normal Reynolds stresses profile along the width of the reactor from the present study closely resemble those computed in the turbulent FCC riser.15 The axial profiles of the normal Reynolds stresses are generally flat in the dense 16214



Article

Figure 12. Axial Reynolds normal stress at (a) y/D = 0.375 and (b) y/D = 1.25.

Figure 13. Lateral Reynolds normal stress at (a) y/D = 0.375 and (b) y/D = 1.25.

characterizes the fluctuation of individual particles, particle clusters, or bubbles in a gas solid system. Research conducted at the Illinois Institute of Technology (IIT) has shown the existence of two types of granular temperatures:47 the first kind is the laminar granular temperature due to individual particle motion computed directly by the CFD code by solving the fluctuating energy equation. The second type is the turbulent (bubble-like) granular temperature computed indirectly from the hydrodynamic velocities and is a result of particle clusters or bubble motion. Thus, this type of granular temperature actually represents the pseudoturbulent fluctuation characteristic of the system under investigation. The sum of the two granular temperatures gives the total granular temperature of the system. Determination of the distribution of granular temperature in the fluidized bed gives valuable information related to particle mixing and transport phenomena. The two types of granular temperatures give rise to two types of mixing: the first one is on the particle level manifested as the laminar granular temperature and the other is on the level of particle clusters and/or

region and increase with height. Meanwhile, the lateral profiles of the normal Reynolds stresses exhibit a dome shaped profile with high stresses predicted in the core of the bed and low stresses near the reactor walls. The anisotropy of the particle fluctuations in the present work is not readily observable and is significantly smaller than that reported in previous works.15,44 This might be due to the lower operating velocities in bubbling fluidized beds as opposed to the much higher velocities used in turbulent and circulating fluidized beds. In the present work, the superficial gas velocity used is only 0.4 m/s as opposed to the range of 3.25−6.1m/s in the turbulent fluidized bed15 and 5.2 m/s in the PSRI riser.44 These high operating superficial gas velocities undoubtedly bring about the observed large anisotropy of the particle fluctuations in these contactors. 4.9. Granular Temperature and Dispersion Coefficients. The concept of granular temperature forms the very core of the Kinetic Theory of Granular Flow as quantities such as particulate pressure, particulate viscosity, and diffusivity are functions of the granular temperature. Granular temperature 16215



Article

bubbles manifested as the turbulent granular temperature. Though both types of granular temperatures have been computed routinely using CFD simulations,15,44 our study is one of the few that tries to explore the influence of the mesoscale on the predicted distributions of the granular temperature and other turbulent parameters as well. The laminar granular temperatures predicted by models G and B are compared with experimental results for the case of Jung and co-workers25 and displayed in Figure 14. The figure

Θturbulent =

1 2 u′yu′y + ux′ux′ 3 3

(36)

Meanwhile, the turbulent kinetic energy in the system is given as Eturbulent =

1 1 1 ux′ux′ + u′yu′y + uz′uz′ 2 2 2

(37)

The turbulent granular temperature is related to the turbulent kinetic energy through a simple relationship written as Θturbulent =

2 Eturbulent 3

(38)

Figure 16 displays the simulated turbulent granular temperature lateral profiles for the two models at dimensionless heights of 0.375 and 1.25 for the case of Yuu et al.24 The profiles bear a similar resemblance to the radial normal stresses. Compared to the laminar granular temperatures, the turbulent granular temperatures are higher which is consistent with previous findings.5,15,25,44 This is as a result of the bubble motion which acts on a global scale unlike individual particle motion which is localized, and since the bubbles may relate to the mesoscale, altering the size of the bubble produces different results. In particular, in this study, the bubble size predicted by model B is smaller which entails that the turbulent granular temperature, which is due to the bubble or particle cluster motion, will also be smaller. Thus, this result indirectly confirms our earlier visual observations of the smaller bubble size predicted by model B. Cody and co-workers48 were perhaps one of the first groups to experimentally measure granular temperature of monodispersed glass spheres using the acoustic shot noise technique. For Geldart B glass spheres fluidized at superficial gas velocities greater than twice the minimum fluidization velocity, these researchers derived an empirical relationship for the granular temperature given as

Figure 14. Comparison of simulated time averaged laminar granular temperature for models G and B with experimental data25 at a bed height of 0.14 m.

clearly demonstrates the need for including a correction factor to the Gidaspow drag model even in the case of Geldart B particles. Model B predicts lower granular temperature values in comparison to model G and gives a better agreement with the experimental data of Jung and co-workers.25 A similar trend is observed for the case of Yuu and co-workers.24 Granular temperature is a function of the fluctuating velocity, and as observed from Figures 9 and 10, drag model G exhibits a wider scatter in the particle velocity. The normal Reynolds stresses, which are the square of the fluctuating velocity, also prove that the fluctuating velocities predicted by drag model G are larger which explains the observed higher granular temperature predicted by the model. Qualitatively, both models exhibit a similar profile predicting nearly constant values in the core of the bed and a sharp increase near the walls. High solids concentration translates to a reduction in the granular temperature since the distance the particles have to travel before they collide (mean free path) is now reduced. Granular temperature increases with a reduction in particle concentration due to an increase in collisions between particles. The turbulent granular temperature is defined as the average of the three squares of the velocity components in the three directions using the following definition: 1 1 1 Θturbulent = ux′ux′ + u′yu′y + uz′uz′ (35) 3 3 3

granular temperature ≈ (110 μm/ds)2 ug 2

(39)

In the present study, glass beads were used as the granular material thereby making the use of this empirical correlation justifiable. From Figure 15, it can be observed that the simulated turbulent or “bubble-like” granular temperature using model B agrees fairly well with the empirical correlation of Cody and co-workers48 for the physical parameters and operating conditions of the present study at both heights. Figure 16 shows the spatially and time averaged profiles of the laminar, turbulent, and total granular temperatures predicted by the two models along the height of the reactor for the present work. From Figure 16a, the laminar granular temperature is low at the bottom of the reactor, i.e., between y/D = 0 and y/D = 2 due to the high solids volume fraction for both models; meanwhile, within this range of the reactor height, the turbulent granular temperature is much higher than the laminar granular temperature (Figure 16b). This is an indication that in the dense phase of the fluidized bed the granular temperature due to the oscillation of bubbles or particle clusters is much more significant compared to that of individual particle oscillation. At higher reactor heights, specifically between y/D = 2 and y/D = 3.5, it can be observed that the laminar granular temperature is now much higher than the turbulent granular temperature for both models due to the fact that this region is dilute, and the so-called mean free path of the particles is increased; therefore, the particle oscillations increase. At even higher heights (i.e., the freeboard region),

The velocity fluctuations in the flow direction (y direction) dominate, and assuming that the velocity fluctuations in the nonflow directions (x and z directions) are equal, the granular temperature is then given as 16216



Article

Figure 15. Lateral profiles of turbulent granular temperature for models G and B at (a) y/D = 0.375 and (b) y/D = 1.25.

Figure 16. Axial profile of spatially and time averaged (a) laminar, (b) turbulent, and (c) total granular temperatures along the height of the reactor for models G and B.

granular temperature. There is a laminar type of the dispersion coefficient due to random oscillations of individual particles associated with the laminar granular temperature and a turbulent dispersion coefficient associated with bubble motion or particle clusters. Kashyap et al.49 derived a formula to calculate the laminar dispersion coefficient for particles in dense conditions. The expression is given as

both the laminar and turbulent granular temperatures drop to zero. As aforementioned, the sum of the two granular temperatures gives the total granular temperature which is displayed in Figure 16c for both drag models. At the bottom of the reactor, the total granular temperature profile resembles that of the turbulent granular temperature profile while at upper heights it resembles the laminar granular temperature profiles because the respective granular temperatures dominate in the respective regions. For all the three granular temperatures, model G predicts higher values in comparison to model B as expected. Summing up this discussion, we can safely say that the effect of the heterogeneous structures is to reduce the granular temperature of the system. This dispersion coefficient measures the quality of mixing which is related to the transfer of mass. Kashyap and coworkers49 (including references cited therein) note that the dispersion coefficient in the fluidized bed varies by five or more orders of magnitude owing to the dependence of the dispersion coefficient on the hydrodynamics which in turn depend on system geometry and properties. Therefore, it seems reasonable that a structure dependent drag coefficient which takes into account the mesoscale structures is necessary for the correct prediction of the dispersion coefficients and hence the mass transfer. Two types of dispersion coefficients exist just like the

D laminar =

π ds 1 Θs 8 αs g0(1 + e)

(40)

It is worth noting that CFD simulations fail to resolve the laminar granular temperature into radial and axial components, but instead, the overall laminar granular temperature is calculated from the code. Thus, from the simulations, only one value is obtained at a particular location. The laminar dispersion coefficients computed at the two dimensionless heights using the two drag models have been presented in Figure 17. In general, the dispersion coefficient rises or increases with height, and it is high in the core of the bed due to low particle volume fraction concentration and low near the walls due to the increased particle concentration. In a similar fashion to the granular temperature profiles, model G predicts higher dispersion coefficients as opposed to model B. In turbulent15 and circulating fluidized bed49 reactors, the laminar dispersions are much 16217



Article

Figure 17. Laminar dispersion for models G and B at heights of (a) y/D = 0.375 and (b) y/D = 1.25.

higher compared to the bubbling fluidized bed reactor owing to higher operating velocities which allow the particles to form dilute regions along most of the height of the reactor. Turbulent dispersion coefficients have been computed previously for bubbling fluidized beds,5 turbulent fluidized beds,15 and CFBs44,49 using the autocorrelation method. In this method, the turbulent dispersion coefficient is estimated by multiplying the normal Reynolds stress with the Lagrangian integral time scale. In turn, the Lagrangian integral time scale is estimated as the area under the curve of the autocorrelation of the fluctuating velocity versus the time lag. The use of this method to gas−solid fluidized beds is not fully understood,49 and neither is there a reliable predictive theory for estimation of turbulent coefficients15 which explains why values reported in the literature vary by orders of magnitude. Thus, no deliberate attempt was made to calculate the turbulent dispersion coefficient using this method; however, the turbulent dispersion coefficient can be estimated using the expression Dturbulent ≅ ((1 − αs)/f)Θturbulent.25 Substituting the relevant values in the expression yields average turbulent dispersion coefficients of 1.24 × 10−3 and 6.85 × 10−4 m2/s for drag models G and B, respectively. As expected, the turbulent dispersion coefficient values are higher than the laminar dispersion coefficients. Once again, the higher values predicted by drag model G are obvious. Dispersion coefficients reported in the literature as reported by Jung and co-workers25 are in the range of 1.0 × 10−4 to 1 m2/s, and our simulated values lie within this range, although it is worth noting that most dispersion coefficient data reported in the literature are in fact the total dispersion coefficient data consisting of the laminar and turbulent components. 4.10. Energy Spectrum. According to the energy cascade concept, energy is transferred from large unstable eddies to smaller eddies. In turn, the smaller eddies break up and transfer energy to yet even smaller eddies, thereby creating an energy cascade. The large eddies contain most of the turbulent energy due to the fact that they have a higher Reynolds number which entails that inertial effects dominate. This process goes on until the viscous forces start to dominate and dissipation of energy occurs. The energy spectrum function, Ei(n), describes how the turbulent kinetic energy is distributed among these eddies of different sizes. The energy spectrum can be conveniently

computed from the Fourier transform of ui′ui′ using the fast Fourier transform (FFT) technique.15 The summation of the energy spectrum function at all frequencies is simply equal to the normal Reynolds stress, i.e.,

∫0

∞

Ei(n)dn = ui′ui′

(41)

where n is the frequency. The energy spectrum has been computed for fluidized beds operating in different regimes.15,44,49 In the present work, the energy spectra predicted by models G and B have been analyzed at two lateral locations near the bed center, i.e., x/D = 0.54, and near the wall region, x/D = 0.92, and at the two dimensionless bed heights, i.e., y/D = 0.375 and y/D = 1.25. Figures 18 and 19 display the simulated axial and lateral energy spectra, respectively, near the bed center region and near the wall region for the two drag models. From both figures, two distinct segments of the energy spectra are identifiable, i.e., the energy containing segment and the inertial segment also referred to as the Kolmogorov range due to the fact that the −5/3 Kolmogorov power law is observed in this segment. At low frequencies, the turbulent energy possessed by the large unstable eddies is high, and with increasing frequency, this unstable eddy breaks and transfers energy to smaller eddies with high frequencies. It has also been reported that at low frequencies gravity waves and internal solids circulation dominate.15 Gravity waves are due to large scale oscillations from the passage of bubbles or particles clusters. A closer inspection of the plots also reveals that the turbulent kinetic energy predicted by model G is higher at all the heights and positions as expected which is also consistent with other parameters tested so far. The −5/3 Kolmogorov power law is obeyed fairly well in the inertial segment for both the lateral and axial directions. Another interesting observation is that the energy spectra increases with height for both drag models at all positions due to the reduction in the solids concentration which is a consequence of the increased oscillations and turbulent kinetic energy. The fact that model G predicts a higher energy spectrum consolidates the earlier observation of the large bubble sizes predicted by the model. The larger bubble sizes predicted by drag model G are responsible for the observed higher energy spectra and other parameters studied so far. Also, from 16218



Article

Figure 18. Simulated axial energy spectrum using drag models G and B at different bed heights for the (a) near central region and (b) near wall region.

Figure 19. Simulated lateral energy spectrum using drag models G and B at different bed heights for the (a) near central region and (b) near wall region.

solids volume fraction in the dense phase predicted by model B results in an increased pressure drop profile in the dense fluidized bed, and the overall pressure drop measured at the reactor inlet agrees very well with the DEM simulation results of Yuu et al.24 Axial profiles of the time averaged particle and gas velocities reveal that the slip velocity predicted by drag model B is slightly higher than that predicted by drag model G. The particle velocities predicted by model G show a wider scatter when compared to the particle velocities predicted by model B as revealed by the plots of local instantaneous particle velocities against local instantaneous solids volume fraction. This observation was later confirmed by the higher predicted Reynolds stresses, granular temperatures, and other turbulent properties. The simulated classical and turbulent granular temperatures using model B agrees very well with the experimental data of Jung et al.25 and the empirical correlation of Cody et al.,48 respectively. The present work has illuminated

Figure 19, it is observed that the lateral energy spectra do not exhibit much variation in comparison with the axial energy spectra, a testimony that most of the energy is in the direction of flow, i.e., the vertical direction.

5. SUMMARY AND CONCLUSIONS In the present work, there was a comparative study of hydrodynamic parameters using the drag models G and B. The first and foremost important conclusion drawn from this work is that the mesoscale structures (accounted for using the heterogeneous index in drag model B) significantly affect predicted fluidized bed hydrodynamics including turbulent properties. In particular, the bubbles, which relate to the mesoscales structures predicted by model B, are smaller in size compared to the bubble size predicted by model G. The mesoscale structures have an effect of lowering the predicted bed height and increasing the solids volume fraction. In addition, the high 16219



Article

(12) Milioli, C. C.; Milioli, F. E.; Holloway, W.; Agrawal, K.; Sundaresan, S. Filtered two-fluid models of fluidized gas-particle flows: New constitutive relations. AlChE J. 2013, 59, 3265. (13) Schneiderbauer, S.; Pirker, S. Filtered and heterogeneity-based subgrid modifications for gas−solid drag and solid stresses in bubbling fluidized beds. AIChE J. 2014, 60, 839. (14) Miller, D. C.; Syamlal, M.; Mebane, D. S.; Storlie, C.; Bhattacharyya, D.; Sahinidis, N. V.; Agarwal, D.; Tong, C.; Zitney, S. E.; Sarkar, A.; Sun, X.; Sundaresan, S.; Ryan, E.; Engel, D.; Dale, C. Carbon capture simulation initiative: A case study in multiscale modeling and new challenges. Annu. Rev. Chem. Biomol. Eng. 2014, 5, 301. (15) Jiradilok, V.; Gidaspow, D.; Damronglerd, S.; Koves, W. J.; Mostofi, R. Kinetic theory based CFD simulation of turbulent fluidization of FCC particles in a riser. Chem. Eng. Sci. 2006, 61, 5544. (16) Lu, B.; Wang, W.; Li, J. Eulerian simulation of gas−solid flows with particles of Geldart groups A, B and D using EMMS-based mesoscale model. Chem. Eng. Sci. 2011, 66, 4624. (17) Benyahia, S. Analysis of model parameters affecting the pressure profile in a circulating fluidized bed. AlChE J. 2012, 58, 427. (18) Shi, Z.; Wang, W.; Li, J. A bubble-based EMMS model for gas− solid bubbling fluidization. Chem. Eng. Sci. 2011, 66, 5541. (19) Hong, K.; Shi, Z.; Wang, W.; Li, J. A structure-dependent multifluid model (SFM) for heterogeneous gas−solid flow. Chem. Eng. Sci. 2013, 99, 191. (20) Liu, X.; Jiang, Y.; Liu, C.; Wang, W.; Li, J. Hydrodynamic modeling of gas−solid bubbling fluidization based on energyminimization multiscale (EMMS) theory. Ind. Eng. Chem. Res. 2014, 53, 2800. (21) Schneiderbauer, S.; Puttinger, S.; Pirker, S. Comparative analysis of subgrid drag modifications for dense gas-particle flows in bubbling fluidized beds. AIChE J. 2013, 59, 4077. (22) Loha, C.; Chattopadhyay, H.; Chatterjee, P. K. Assessment of drag models in simulating bubbling fluidized bed hydrodynamics. Chem. Eng. Sci. 2012, 75, 400. (23) Wang, J.; Liu, Y. EMMS-based Eulerian simulation on the hydrodynamics of a bubbling fluidized bed with FCC particles. Powder Technol. 2010, 197, 241. (24) Yuu, S.; Umekage, T.; Johno, Y. Numerical simulation of air and particle motions in bubbling fluidized bed of small particles. Powder Technol. 2000, 110, 158. (25) Jung, J.; Gidaspow, D.; Gamwo, I. K. Measurement of two kinds of granular temperatures stresses and dispersion in bubbling beds. Ind. Eng. Chem. Res. 2005, 44, 1329. (26) Li, J.; Kwauk, M. Particle Two Phase Flow: The EnergyMinimization Multi Scale Method; Metallurgical Industry Press: Beijing, 1994. (27) Gao, X.; Wang, L.-J.; Wu, C.; Cheng, Y.-W.; Li, X. Novel Bubble−Emulsion Hydrodynamic Model for Gas−Solid Bubbling Fluidized Beds. Ind. Eng. Chem. Res. 2013, 52, 10835. (28) Yang, N.; Wang, W.; Ge, W.; Li, J. Modeling of meso-scale structures in particle-fluid systems: The EMMS/CFD approach. China Particuol. 2005, 3, 78. (29) Zhou, Q.; Wang, J. Coarse grid simulation of heterogeneous gassolid flow in a CFB riser with EMMS drag model: Effect of inputting drag correlations. Powder Technol. 2014, 253, 486. (30) Khan, A. R.; Richardson, J. F. The resistance to motion of a solid sphere in a fluid. Chem. Eng. Commun. 1987, 62, 135. (31) Syamlal, M.; Pannala, S. Multiphase continuum formulation for gas-solids reacting flows. In Computational Gas-Solids Flows and Reacting Systems - Theory, Methods and Practice; Pannala, S., Syamlal, M., O’Brien, T. J., Eds.; IGI Global: Hershey, PA, 2011; pp 1. (32) Cloete, S.; Zaabout, A.; Johansen, S. T.; van Sint Annaland, M.; Gallucci, F.; Amini, S. The generality of the standard 2D TFM approach in predicting bubbling fluidized bed hydrodynamics. Powder Technol. 2013, 235, 735. (33) Gel, A.; Li, T.; Gopalan, B.; Shahnam, M.; Syamlal, M. Validation and uncertainty quantification of a multiphase computational fluid dynamics model. Ind. Eng. Chem. Res. 2013, 52, 11424.

the concept of energy minimization and has helped to shed more light on mesoscale modeling in bubbling fluidized beds, hopefully contributing to improved design, scale up, and optimization of gas−solid bubbling fluidized beds.

■

ASSOCIATED CONTENT

S Supporting Information *

The unsteady EMMS bubble based model equations, summary of simulation results in tabular form, and nomenclature. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.:+86-0571-87951227. Fax: +86-0571-87951227. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge the support and encouragement of National Natural Science Key Foundation of China 21236007 and National Basic Research Program of China 2012CB720500.

■

REFERENCES

(1) 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. (2) McKeen, T.; Pugsley, T. Simulation and experimental validation of a freely bubbling bed of FCC catalyst. Powder Technol. 2003, 129, 139. (3) Gibilaro, L. G.; Di Felice, R.; Waldram, S. P.; Foscolo, P. U. Generalized Friction Factor and Drag Coefficient for Fluid Particle Interactions. Chem. Eng. Sci. 1985, 40, 1817. (4) Li, T.; Pougatch, K.; Salcudean, M.; Grecov, D. Numerical simulation of horizontal jet penetration in a three-dimensional fluidized bed. Powder Technol. 2008, 184, 89. (5) Chalermsinsuwan, B.; Gidaspow, D.; Piumsomboon, P. Two and three dimensional CFD modeling of Geldart A particles in a thin bubbling fluidized bed comparison of turbulence and dispersion coefficients. Chem. Eng. J. 2011, 171, 301. (6) 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. (7) Gidaspow, D. Multiphase flow and fluidization: Continuum and kinetic theory descriptions; Academic Press: San Diego, 1994. (8) Hong, K.; Wang, W.; Zhou, Q.; Wang, J.; Li, J. An EMMS-based multi-fluid model (EFM) for heterogeneous gas−solid riser flows: Part I. Formulation of structure-dependent conservation equations. Chem. Eng. Sci. 2012, 75, 376. (9) Ge, W.; Wang, W.; Yang, N.; Li, J.; Kwauk, M.; Chen, F.; Chen, J.; Fang, X.; Guo, L.; He, X.; Liu, X.; Liu, Y.; Lu, B.; Wang, J.; Wang, J.; Wang, L.; Wang, X.; Xiong, Q.; Xu, M.; Deng, L.; Han, Y.; Hou, C.; Hua, L.; Huang, W.; Li, B.; Li, C.; Li, F.; Ren, Y.; Xu, J.; Zhang, N.; Zhang, Y.; Zhou, G.; Zhou, G. Meso-scale oriented simulation towards virtual process engineering (VPE)The EMMS Paradigm. Chem. Eng. Sci. 2011, 66, 4426. (10) Igci, Y.; Andrews, A. T.; Sundaresan, S.; Pannala, S.; O’Brien, T. Filtered two-fluid models for fluidized gas-particle suspensions. AlChE J. 2008, 54, 1431. (11) Igci, Y.; Sundaresan, S. Constitutive models for filtered two-fluid models of fluidized gas−particle flows. Ind. Eng. Chem. Res. 2011, 50, 13190. 16220



Article

(34) Bokkers, G.; van Sint Annaland, M.; Kuipers, J. Mixing and segregation in a bidisperse gas−solid fluidised bed: A numerical and experimental study. Powder Technol. 2004, 140, 176. (35) Goldschmidt, M. J. V.; Beetstra, R.; Kuipers, J. A. M. Hydrodynamic modelling of dense gas-fluidised beds: Comparison and validation of 3D discrete particle and continuum models. Powder Technol. 2004, 142, 23. (36) Ku, X.; Li, T.; Løvås, T. Influence of drag force correlations on periodic fluidization behavior in Eulerian−Lagrangian simulation of a bubbling fluidized bed. Chem. Eng. Sci. 2013, 95, 94. (37) Gidaspow, D.; Jung, J.; Singh, R. K. Hydrodynamics of fluidization using kinetic theory: An emerging paradigm. Powder Technol. 2004, 148, 123. (38) Xiong, Q.; Kong, S.-C. Modeling effects of interphase transport coefficients on biomass pyrolysis in fluidized beds. Powder Technol. 2014, 262, 96. (39) Panday, R.; Shadle, L. J.; Shahnam, M.; Cocco, R.; Issangya, A.; Spenik, J. S.; Ludlow, J. C.; Gopalan, B.; Shaffer, F.; Syamlal, M.; Guenther, C.; Karri, S. B. R.; Knowlton, T. Challenge problem: 1. Model validation of circulating fluidized beds. Powder Technol. 2014, 258, 370. (40) Benyahia, S. On the effect of subgrid drag closures. Ind. Eng. Chem. Res. 2009, 49, 5122. (41) Lijie, Y.; Shuyan, W.; Huilin, L.; Shuai, W.; Pengfei, X.; Lixing, W.; Yurong, H. Flow of gas and particles in a bubbling fluidized bed with a filtered two-fluid model. Chem. Eng. Sci. 2010, 65, 2664. (42) Qi, X.; Zhu, H.; Zhu, J. Demarcation of a new circulating turbulent fluidization regime. AlChE J. 2009, 55, 594. (43) Du, B.; Warsito, W.; Fan, L.-S. Behavior of the dense-phase transportation regime in a circulating fluidized bed. Ind. Eng. Chem. Res. 2006, 45, 3741. (44) Chalermsinsuwan, B.; Piumsomboon, P.; Gidaspow, D. Kinetic theory based computation of PSRI riser: Part IEstimate of mass transfer coefficient. Chem. Eng. Sci. 2009, 64, 1195. (45) Shuai, W.; Guangbo, Z.; Guodong, L.; Huilin, L.; Fixiang, Z.; Tianyu, Z. Hydrodynamics of gas-solids risers using cluster structuredependent drag model. Powder Technol. 2014, 254, 214. (46) Sun, J.; Wang, J.; Yang, Y. CFD investigation of particle fluctuation characteristics of bidisperse mixture in a gas−solid fluidized bed. Chem. Eng. Sci. 2012, 82, 285. (47) Tartan, M.; Gidaspow, D. Measurement of granular temperature and stresses in risers. AlChE J. 2004, 50, 1760. (48) Cody, G. D.; Goldfarb, D.; Storch, G. V., Jr.; Norris, A. N. Particle granular temperature in gas fluidized beds. Powder Technol. 1996, 87, 211. (49) Kashyap, M.; Gidaspow, D.; Koves, W. J. Circulation of Geldart D type particles: Part I − High solids fluxes. Measurements and computation under solids slugging conditions. Chem. Eng. Sci. 2011, 66, 183.

16221


Numerical Simulations of a Bubbling Fluidized Bed Reactor with an

Recommend Documents