Numerical simulations of flow structure and heat transfer in a central jet bubbling fluidized bed

Powder Technology 269 (2015) 139–152

Contents lists available at ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

Numerical simulations of flow structure and heat transfer in a central jet bubbling fluidized bed Musango Lungu, Jingdai Wang ⁎, Yongrong Yang State Key Laboratory of Chemical Engineering, Department of Chemical and Biological Engineering, Zhejiang University, Hangzhou 310027, China

a r t i c l e

i n f o

Article history: Received 11 February 2014 Received in revised form 9 July 2014 Accepted 29 August 2014 Available online 7 September 2014 Keywords: Fluidization Simulation Granular temperature Energy spectrum Heat transfer

a b s t r a c t An attempt has been made to model the flow structure and to predict heat transfer coefficients in a gas– solid bubbling fluidized bed operated with a central jet using a two fluid model with closures from the kinetic theory of granular flow. Quantities such as the fluid-to-particle heat transfer coefficient are not easily obtained practically from experiments with a high degree of accuracy thereby making computational methods attractive. The CFD model has been verified using experimental bubble properties obtained from the literature. Axial and lateral normal Reynolds stresses, energy spectra and granular temperatures have been computed. The simulations show that the maximum local instantaneous fluid-to-particle heat transfer coefficient occurs in the wake of the bubble. Heat transfer coefficients at the center of the bed exhibit high oscillations compared to the near wall region. However time averaged values in the near wall region are larger compared to the bed center. The average heat transfer coefficient exhibits a maximum value with the variation of the jet velocity. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Compared to other modes of contacting, fluidized bed reactors offer high rates of heat and mass transfer between the fluid and solid phases. It is not surprising therefore that fluidized bed technology is widely employed in chemical, petrochemical, metallurgical and energy industries [1]. The excellent features of these reactors are largely attributed to the presence of bubbles [2] and in industrial scale fluidized bed reactor bubbles emerge from discrete holes or other types of orifices in the gas distributor plate. Kuipers et al. [3] note that bubble characteristics are dependent on the initial bubble generated from the orifice. They further state that mass and heat transfer processes are seriously affected by the mechanism of bubble formation. Numerous studies (experimental and theoretical) have been conducted on the formation of bubbles from a single orifice. Harrison and Leung [4] are perhaps one of the first authors to develop a model for the formation of gas bubbles at a single orifice in a fluidized bed. However their model assumed the non-existence of gas leakage at the bubble surface into the emulsion phase an assumption which was later proved otherwise by Nguyen and Leung [5]. Other models describing the formation of gas bubbles at a single orifice are due to Yang et al. [6] and Caram and Hsu [7]. The

⁎ Corresponding author. Tel.:+86 571 87951227; fax: +86 571 87951227. E-mail address: [email protected] (J. Wang).

http://dx.doi.org/10.1016/j.powtec.2014.08.067 0032-5910/© 2014 Elsevier B.V. All rights reserved.

model of Yang and co-workers postulates that the superficial velocity of gas leakage at the bubble-emulsion boundary is equal to the minimum fluidization velocity while that of Caram and Hsu postulates that the superficial exchange velocity is dependent of the pressure gradient. Increased computational capabilities coupled with reduction in hardware costs have given rise to the use of computational fluid dynamics (CFD) in the modeling and simulation of multiphase systems. Therefore for the last two decades or so studies on the formation of bubbles at a single orifice have been executed using CFD. Two basic approaches namely Eulerian–Eulerian (granular flow models) and Eulerian– Lagrangian (discrete particle models) are commonly used to model gas–solid flows. Kuipers and co-workers [3] used a two fluid model (TFM) to study the theoretical bubble formation at a single orifice in a two dimensional bed. The model predicted bubble sizes, formation times and shapes satisfactorily and the results were in good agreement with experimental results and approximate models. Later Nieuwland et al. [2] extended the work of Kuipers and co-workers to investigate the effect of particle properties on the bubble growth process for Geldart B particles. They came to a conclusion that the hydrodynamic model gives a satisfactory good description of the bubble growth process for several particle types. Bokkers et al. [8] used a discrete particle model (DPM) to study the extent of mixing and segregation induced by a single bubble in a monodisperse and bidisperse fluidized bed at incipient conditions and in freely bubbling fluidized beds in addition to experimental studies. They observed that the bubble size predicted by

140

M. Lungu et al. / Powder Technology 269 (2015) 139–152

the DPM and the extent of segregation induced by a single bubble in a bidisperse fluidized bed strongly depends on the drag model used. Patil et al. [9] compared the constant viscosity model (CVM) with the kinetic theory of granular flow (KTGF) in the study of bubble formation at a single orifice with experimental data. Their study revealed that bubble growth is mainly determined by the drag experienced by the gas percolating through the bubble-emulsion boundary. Recently Kumar et al. [10] used the Eulerian–Eulerian modeling approach to study the effect of the wall on the shape of the bubble, bubble wake during formation and the bubble rise for both symmetric and asymmetric injections. In addition to the modeling and simulation of bubble behavior from a single orifice, several other modeling and simulation studies [11–13] have been reported focusing the wall-to-bed heat transfer coefficients in bubbling fluidized beds in which the orifice is located near a heated wall and the bulk of the bed is kept at a constant temperature. Results from such simulations agreed reasonably well with experimental results from carefully controlled experiments. Despite the great advantages offered by CFD very few studies have been reported on fluid-to-particle heat transfer modeling. Probably this is due to the fact that simulation results from such studies are difficult to validate against experimental data. Inaccuracies in the measurement of particle and gas temperatures have made the determination of fluid-to-particle heat transfer coefficients difficult to measure. Furthermore due to the complex flow pattern in fluidized beds, the fluid-toparticle heat coefficients reported in literature vary a great deal depending on the flow assumption used [14]. In addition there are currently no experimental techniques for obtaining energy profiles in fluidized beds without perturbing the system thereby rendering computational methods attractive [15]. The objective of the present work is to analyze the flow structure including the turbulent properties such as the Reynolds stresses, granular temperatures and energy spectra and how they affect the predicted fluid-to-particle heat transfer coefficients at different positions in the bed. Such an approach would give valuable information concerning the coupling between the bed hydrodynamics and the heat transfer coefficients. CFD as a complementary tool to experimentation and reactor models [16] can greatly help in the understanding of transport mechanisms in the reactor and obtain useful information. The flow structure in the fluidized beds can be described in detail using properties such as the solid volume fraction distribution, void fraction distribution and laminar granular temperature as well as turbulent properties such as the turbulent granular temperature and the energy spectra. Examination of the flow structure can give an idea about regions of high and low heat transfer which is necessary for the optimal design and scale up of fluidized bed reactors. The granular temperature gives information about the fluctuation or oscillation of individual particles, bubbles and clusters [17]. Moreover the temporal variation of the heat transfer coefficients at specific positions in the bed can be analyzed and related to the oscillatory behavior of the individual particles, bubbles and/or clusters. Such an approach enables a better and clear understanding of the coupling between the bed hydrodynamics and the predicted heat transfer coefficients. The paper is organized as follows. Section 2 gives the continuum governing equations and constitutive relations. Section 3 describes the simulation set-up together with the procedure for processing the data from the simulations. Section 4 is the results and discussion section. Firstly the model is validated using experimental data and the initialstartup bubble properties are considered for this purpose. The simulated results are also compared with the theoretical model according to Caram and Hsu [18]. A study on the choice of the drag model used has also been done. Next the flow structure of the bed is investigated including turbulent properties such as the turbulent granular temperature and energy spectrums. Finally predicted local and time averaged heat transfer coefficients are discussed with relation to the flow structure. In Section 5 the summary and conclusions of the present work are given.

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 uk Þ ¼ 0 ∂t k k

ð1Þ

in each computational cell the total volume fraction of all phases must be equal to unity, i.e. n X

α k ¼ 1:

ð2Þ

k¼1

Momentum conservation for solid phase: ∂ ðα ρ u Þ þ ∇  ðα s ρs us us Þ ¼ −α s ∇p−∇ps þ ∇  τs ∂t s s s  

ð3Þ

þα s ρs g þ βgs ug −us

Momentum equation for gas phase:    ∂ α g ρg ug þ ∇  α g ρg ug ug ¼ −α g ∇p þ ∇  τ g ∂t  

ð4Þ

þα g ρg g−βgs ug −us

Thermal energy conservation for solid phase:   ∂ ðα s ρs H s Þ þ ∇  ðα s ρs us H s Þ ¼ −∇  α s κ s ∇T s þ hsg T g −T s ∂t

ð5Þ

Thermal energy conservation for gas phase:      ∂ α g ρg H g þ ∇  α g ρg ug H g ¼ −∇  α g κ g ∇Τ s −hsg T g −T s ∂t

ð6Þ

2.2. Granular temperature equation The algebraic form of the granular temperature transport equation was adopted after an initial study to save computational cost and maintain accuracy of results. It is applicable for dense fluidized beds where the convection and the diffusion term can be neglected under the premise that production and dissipation of granular energy are in equilibrium.   0 ¼ −ps I þ τs : ∇us −γ s −3βgs Θs

ð7Þ

2.3. Constitutive equations for interphase momentum transfer Syamlal O'Brien [19]:

βgs ¼

   3 α s α g ρg Res   C ug −us  4 v2rs ds D vrs

12     ρg ug −us ds 4:8 B C ffiA ; Res ¼ C D ¼ @0:63 þ qffiffiffiffiffiffiffiffiffiffiffiffi  Res μg

ð8Þ

0

vrs

ð9Þ

M. Lungu et al. / Powder Technology 269 (2015) 139–152

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

vrs ¼ 0:5 A−0:06Res þ ð0:06Res Þ2 þ 0:12Res ð2B−AÞ þ A2

ð10Þ

where 4:14

A ¼ αg  B¼

Pα 1:28 αg ≤ g α Qg

0:85

α g N 0:85

Q ¼ 1:28 þ

logðP Þ : logð0:85Þ

ps ¼ α s ρs ½1 þ 2ð1 þ eÞα s g 0 Θs :

ð12Þ

2.4.2. Solid shear stresses The solid stress tensor contains shear and bulk viscosities arising from particle momentum exchange due to translation and collision. A frictional component of viscosity to account for the viscous-plastic transition that occurs when particles of a solid phase reach the maximum solid volume fraction has also been included. The bulk solid viscosity which accounts for the resistance of the granular particles to compression and expansion is given by Lun et al. [22] as:

The adjusted Syamlal O'Brien drag function is obtained by an iterative procedure which determines the new model coefficients B, Q and P. The iterative procedure seeks to minimize the difference between the theoretical and experimental minimum fluidization velocities, exp th exp mathematically |uth mf − umf | → 0, where umf and umf are the theoretical and experimental minimum fluidization velocities respectively. The default values for P and Q in the Syamlal O'Brien drag model are 0.8 and 2.65 respectively. The new coefficients after carrying out the iterative procedure are P = 0.633 and Q = 4.08967. Gidaspow [20]:

βWen‐Yu

βErgun ¼ 150

CD ¼

 24

  α s 1−α g μ g α g d2s

Res ½1þ0:15ðRes Þ 0:44;Res N1000

0:687

þ 1:75

    α s ρg ug −us  ds

ð14Þ

for α g ≤0:8

;Res b1000

arctan½262:5ðα s −0:2Þ þ 0:5 π

ð15Þ

ð16Þ

ð17Þ

Thus the momentum transfer coefficient can now be expressed as:   βgs ¼ φgs βErgun þ 1−φgs βWen‐Yu :

rffiffiffiffiffiffi 4 2 Θs λs ¼ α s ρs ds g 0 ð1 þ eÞ : 3 π

ð18Þ

μ s;kin

pffiffiffiffiffiffiffiffi

2 10ρs ds πΘs 4 ¼ 1 þ ð1 þ eÞg 0 α s α s 5 96α s ð1 þ eÞg 0

μ s:col ¼

  h i 2 T τs ¼ α s μ s ∇us þ ð∇us Þ þ α s λs − μ s ð∇  us ÞI 3

rffiffiffiffiffiffi 4 Θs α s ρs ds g 0 ð1 þ eÞ α: 5 π s

ð23Þ

ð24Þ

For the frictional component the expression by Schaeffer [23] is used and it is given as p sinϕ μ s;fr ¼ spffiffiffiffiffiffiffi 2 I2D

ð25Þ

where ϕ is the internal angle of friction and I2D is the second invariant of the deviatoric stress tensor. 2.4.3. Radial distribution function The radial distribution function a correction factor that modifies the probability of collisions between the particles when the solid granular phase becomes dense has been modeled according to Ma and Ahmadi [24] and the expression is 2

The kinetic theory of granular flow analogous to the kinetic theory of gases has used to close the solid phase rheology and the solid phase properties such as the viscosity and solid pressure are expressed in terms of a property called the granular temperature. The phase stress tensors for the respective solid and gas phase are expressed in Newtonian form and given as:

g0 ¼

3

1 þ 2:5α s þ 4:59α s þ 4:52α s :  3 0:678 1− α αs

ð26Þ

s; max

ð19Þ

2.4.4. Collisional dissipation of energy The collisional dissipation of energy represents the rate of energy dissipation within the solid phase due to collisions between particles and the expression due to Lun et al. [22] has been adopted in the present work. It is given as:

ð20Þ

  12 1−e2 g 0 2 3=2 pffiffiffi γs ¼ ρs α s Θs ds π

and  T i 2   − α g μ g ∇  ug I: τg ¼ α g μ g ∇ug þ ∇ug 3

ð22Þ

The shear viscosity of solids is a sum of the kinetic component, collisional component and frictional component i.e. μs = μs,kin + μs,col + μs,fr. The kinetic and collisional viscosity are modeled according to Gidaspow [20] and given as

2.4. Kinetic theory of granular flow (KTGF)

h

ð21Þ

and

To avoid discontinuity in the Gidaspow model, a switch function that gives a rapid but smooth transition from one regime to another has been adopted [21]. φgs ¼

2.4.1. Solid pressure The solid pressure which represents the normal solid phase forces due to particle–particle interactions is modeled according to Lun et al. [22] and is given as:

ð11Þ

ð13Þ

    α s α g ρg ug −us  −2:65 3 ¼ CD αg for α g N0:8 4 ds

141

ð27Þ

142

M. Lungu et al. / Powder Technology 269 (2015) 139–152

recent studies by Deen et al. [30] and Patil et al. [15] have validated the Gunn correlation using direct numerical simulations (DNS) and discrete element method (DEM) respectively and the agreement between the simulations and the empirical equation was excellent.

Table 1 Nusselt number correlations. Ranz and Marshall [26] Fixed bed

Nu ¼ 2 þ 1:8Res 1=2 Pr1=3

BFB

Nu ¼ 2 þ 0:6Res 1=2 Pr1=3 ( )

Nelson and Galloway [29] Fixed bed, BFB

2ζ þ

Nu ¼

2ζ 2 ð1−α g Þ

1=3 2

½1−ð1−αg Þ1=3 

−2

ζ − tanh 1=3 1−ð1−α g Þ

2.5.1. Phase thermal conductivities An expression for the effective thermal bed conductivity which consists of the gas phase conductivity and solid phase conductivity was derived by Kuipers et al. [11] and has been implemented in our model. The gas phase and solid phase conductivities are given by:

tanh ζ ζ

where

1=2 1 ε Pr1=3 ζ¼ 1=3 −1 2 Res ð1−αg Þ Nu = 0.07Res

Cybulski et al. [28] Fixed bed Gunn Fixed bed, BFB Wakao et al. [27] Fixed bed

+ 0.7Re0.2 Pr1/3) s 1.2α2g)Re0.7 Pr1/3 s

5α2g)(1

Nu = (7 − 10αg + + (1.33 − 2.4αg +

κg ¼

 qffiffiffiffiffiffiffiffiffiffiffiffiffi 1− 1−α g κ g;o

ð28Þ

αg

Nu ¼ 2 þ 1:1 Pr1=3 Res 0:6

κs ¼ 2.5. Constitutive equations for the thermal energy balance

ð29Þ

where 2

3

6 A−1 B A B−1 7 −0:5ðB þ 1Þ5 K¼ 4

In − B 1− B 2 A B B A 1− 1− A A 2

κ s;o kg;o

ð31Þ

0.57m

0.5m

A¼

ð30Þ

Freeboard

0.5m

In the energy balance equation the volumetric interphase heat transfer coefficient, hsg, is product of the specific interfacial area and the fluid-to-particle heat transfer coefficient, h. Assuming the particle is spherical which a reasonable assumption for glass beads then the following conversion can be made hsg ¼ ð6α s =ds Þh. The volumetric interphase heat transfer coefficient can be conveniently estimated from empirical correlations. There are numerous empirical correlations reported in literature for the estimation of both packed bed and fluidized bed fluid-to-particle heat transfer coefficients. These correlations relate the Nusselt number to the Reynolds number, Prandtl number and in some cases the voidage. Table 1 below lists some of the correlations. The correlations have been plotted in Fig. 1 using data for the minimum fluidization condition and from the figure two distinct groups can be seen. For the first group which includes correlations due to Gunn [25], Ranz and Marshall [26] and Wakao et al. [27], the Nusselt number approaches a constant value as the Reynolds number approaches zero. For the second group consisting of correlations due to Cybulski et al. [28] and Nelson and Galloway [29], the Nusselt number approaches zero as Reynolds number approaches zero and this is attributed to the negligence of axial dispersion [25]. In the present work the correlation due to Gunn has been adopted as it can be applied over a wide range of porosity i.e. 0.35–1 and for Res of up to 105. The Gunn correlation also satisfies several other asymptotic conditions and furthermore

ðβA þ ð1−βÞK Þκ g;o pffiffiffiffiffiffiffiffiffiffiffiffiffi 1−α s

Bed at minimum fluidization

100

Pr = 0.7155 αg = 0.402

Fluidized bed Fixed bed

ε = 0.6

10

Nu

single sphere

α mf = 0.402

1

Ranz and Marshall [26] Nelson and Galloway [29] Cybulski et al [28] Gunn [25] Wakao et al [27]

0.1

0.01 0.1

1

10

Res Fig. 1. Nusselt number correlations

100

umf

d0

u0 Fig. 2. Geometry of fluidized bed with central jet.

umf

M. Lungu et al. / Powder Technology 269 (2015) 139–152

B ¼ 1:25

1−α g αg

!10=9

0.25

Experimental data [9] Adjusted Syamlal O'Brien Gidaspow Theoretical model [7] Syamlal O'Brien

ð32Þ 0.20

−3

ð33Þ

0.15

De (m)

β ¼ 7:26  10

143

K is the effective thermal conductivity of a cylinder. This cylinder consists of one particle and the fluid phase. β is the ratio of the particle contact area to the particle surface area. B is the factor of deformation.

0.10

0.05

3. Simulation setup

Initially the bed is operated at minimum fluidization condition and air at a velocity many times the minimum fluidization velocity is introduced through the orifice resulting in the formation of bubbles. At the top of the reactor, a pressure outlet boundary condition is assigned. The left and right walls are designated as no slip for air and the partial slip wall boundary condition of Johnson and Jackson has been adopted for the particulate phase with a specularity coefficient of 0.5 which is commonly used in the simulation of bubbling fluidized beds and a particle-wall restitution coefficient of 0.9. In addition the thermal energy boundary condition at the walls is set to adiabatic. Initial values of the pressure profile, solid volume fraction and gas velocity are patched in the flow field.

0.00

0.05

0.10

0.15

0.20

0.25

t (s) Fig. 3. Bubble growth with time.

computed area for which αg N 0.85. A 0.25 second animation video was then created and exported to image processing software, image pro plus 6.0 which was used to calculate the equivalent bubble diameter and other bubble properties. 4. Results and discussion The results reported can be classified into two broad sections, the first section deals with the bubble and particle motion and the second section deals with the interphase fluid-to-particle heat transfer coefficients. Results have mainly been reported for the center and near wall regions at three bed heights. 4.1. Initial startup bubble Fig. 3 shows the bubble diameter as a function of time for the initial startup bubble from numerical simulations using different drag models at a jet velocity of 10 m/s. Also shown is the experimental data extracted from the paper by Patil et al. [9] and the result from the Caram and Hsu approximate model [7]. All the drag models capture the bubble growth with time reasonably well compared with the experimental results. The Caram and Hsu approximate model result is in excellent agreement with the simulation results which is consistent with the findings of Patil [9] and Olaofe [31]. This is expected since the assumptions on

100000

3

3.1. Initial and boundary conditions

0.00

drag function (kg/m s)

The simulated bed geometry is shown in Fig. 2 and is the same as that used by Patil et al. [9]. It is 2D planar with a width of 0.57 m and a height of 1 m with a jet orifice of 0.015 m width installed at the center of the bottom wall. A uniform grid of 0.01 m in the vertical direction and 0.0075 m in the horizontal direction resulted in a grid independent solution. Glass beads with a mean particle diameter of 500 μm and density of 2660 kg/m3 have been used as the granular media and air at ambient temperature is the fluidizing gas. The mesh was created in a CAD program GAMBIT 2.3.16 and exported into FLUENT 6.3.26, a CFD commercial software package used to execute the numerical simulations. A constant time step of 1 × 10−4 seconds has been used throughout. Initially the bed was kept at the minimum fluidization condition and then after 1 second, gas with a velocity higher than the minimum fluidization velocity was injected through the orifice in the center of the bed. Such a procedure minimizes the shock effect of the initial startup bubble and enhances convergence. The start of the bubble injection was considered to be t = 0 second ignoring the first second during which the bed was simulated at incipient conditions and subsequently the simulations were run for 15 seconds of real flow time and quantities were time averaged for the last 10 seconds as it took 5 seconds for the startup transients to subside. It took nearly 10 days to simulate 15 seconds of real flow time on a Sugon I620r-G server employing 16 cores. Simulations with the first order discretization scheme for convective terms resulted in physically unrealistic pointed bubble shapes owing to numerical diffusion; therefore the second order upwind scheme was used to discretize momentum, granular temperature and energy while the QUICK scheme was applied for the spatial discretization of volume fraction, the second order implicit scheme for transient formulation and Phase Coupled SIMPLE for pressure–velocity coupling.

10000

Gidaspow Syamlal O'Brien Adjusted Syamlal O'Brien 1000

3.2. Data processing of simulation results The simulation case and data files were automatically saved every 0.01 seconds and then imported into Tecplot, a CFD visualization software. The minimum voidage contour level was adjusted to 0.85 and the number of levels changed to 1.This converts the contour plot to a binary image. In our study the equivalent bubble diameter is defined as the diameter of a circle possessing the same area as the numerically

100 0.3

0.4

0.5

0.6

0.7

0.8

0.9

voidage Fig. 4. Comparison of the drag functions for different drag models.

1.0

144

M. Lungu et al. / Powder Technology 269 (2015) 139–152

Table 2 Bubble detachment time tb and corresponding equivalent bubble diameter De.

Adjusted Syamlal O'Brien Gidaspow Syamlal O'Brien Caram and Hsu approximate model Patil et al. [9] Experimental data [9]

tb (s)

De (m)

0.20 0.20 0.20 0.21 0.188 0.17

0.173 0.171 0.157 0.188 0.172 0.145

which the Caram and Hsu are based are taken care of in the Euler–Euler approach through the interphase momentum exchange coefficient term. The adjusted Syamlal O'Brien consistently predicts the largest bubble size throughout followed by the Gidaspow model and lastly the Syamlal O'Brien model. In order to understand the observed differences predicted by the drag models we analyze the drag function predicted the drag models. Fig. 4 gives a comparison of the drag coefficient as a function of porosity for the different models at an arbitrary slip velocity of 10 m/s. The adjusted Syamlal O'Brien model predicts the largest drag function throughout the whole voidage range considered. This explains why the adjusted Syamlal O'Brien predicts the largest bubble sizes compared to the other drag models. The Gidaspow and Syamlal O'Brien models intersect each other at void fractions of 0.5 and 0.8. Within this void range the drag function predicted by the Gidaspow model is higher than that due to the Syamlal O'Brien model but still less than that due to the adjusted Syamlal O'Brien model. For bubbling fluidized beds the mean voidage range is reported to be 0.45–0.65 [32]; therefore the porosity range of 0.5–0.8 represents a typical situation prevailing in bubbling fluidized beds. A pattern between the predicted bubble size and the predicted drag function is evident. The higher the predicted drag function, the larger the size of the bubble. It can also be observed from the graph that all drag models approach the single particle drag value with increasing voidage. Table 2 is a summary of the bubble detachment time and the corresponding equivalent bubble diameter at a jet velocity of 10 m/s. From the table it is observed that all three drag models predict the same bubble detachment time but different bubble sizes with the adjusted Syamlal O'Brien predicting the largest bubble size and the Syamlal O'Brien the least. However the bubble detachment time for the drag models as well as the Caram and Hsu model is slightly higher than the experimental value. Bubble detachment occurs when the buoyancy force is balanced by the weight of the bubble. In the simulations the drag force is significantly higher compared to that in the experimental system and thus the higher drag forces keep the bubble attached consequently causing a delay in the bubble detachment from the orifice. The Syamlal O'Brien model gives the closest result to the experimental value for the bubble size among the three models. Finally we compare the bubble shapes predicted by the three drag models to a photograph captured from experimental studies adapted from Patil et al. [9]. Fig. 5 shows the experimental and simulated bubble shapes for different drag models at a jet velocity of 10 m/s and time of 0.2 seconds. The main observable difference is at the bottom of the bubble. The adjusted Syamlal O'Brien shows a pointed bubble base as the bubble is just been detached from the orifice. The base of the bubble near the orifice predicted by the Gidaspow model is rather flat. A close inspection of the

Experiment [9]

Adjusted Syamlal O’Brien

figure reveals that the Syamlal O'Brien compares very well to the experimental photograph even capturing the slit at the bottom of the bubble. The slit is caused by gas leakage as the bubble rises. Based on this study the Syamlal O'Brien model was adopted and used in the remainder of the simulations. 4.2. Bubble frequency Fig. 6 shows simulated as well as experimental porosity fluctuations adapted from Kuipers et al. [33]. The void fraction measurements from the simulations were used to calculate the simulated bubble frequency. A bubble threshold of 0.85 has been used to define a bubble as aforementioned. The bubble frequency is calculated from the equation f b ¼ Nb =t s , where fb is the bubble frequency, Nb is the total number of bubbles monitored at the “numerical probe” and ts is the total sampling time. In both the simulations and the experiment a sampling time of 2 seconds was used. Neglecting the first initial startup bubble, the bubble frequency predicted by the simulations is nearly 10 Hz while the experiment gives a bubble frequency of 6 Hz. Clearly the simulations over predict the bubbling frequency and this can be attributed to the simplification of the simulation domain. The experiments were conducted in a pseudo 2D bed while the simulations were executed in a 2D Cartesian coordinate system. Neglecting front and back wall effects in 2D simulations results in small and fast moving bubbles [34,35]. Since the bubble rise velocities of the bubbles in 2D simulations are higher than those in a pseudo 2D or 3D system, more bubbles are observed at the probe location in the simulation. This explains the observed higher bubbling frequency from the simulation. However if the bubble frequency is considered to be the reciprocal of the bubble detachment time, then the bubble frequency is 5 Hz for a jet velocity of 10 m/s which is pretty close to the observed experimental value. Theoretically from the modified Davidson and Schuler model h 2  i1=3 16C 0 u0 d0 for 2D geometry [33] given by the expression t b ¼ , π g2 where C0 is taken to be one gives a bubble frequency of 5.01 Hz which matches quite well with simulated bubble frequency derived from the reciprocal of the bubble detachment time. The bubble frequency increases with decreasing jet velocity as shown in Fig. 7. This is a result of the reduction of the bubble detachment time with decreasing jet velocity. As aforementioned bubble detachment occurs when the buoyancy force is balanced by the weight of the bubble. However at lower velocities the drag force is less and thus the bubble takes a shorter time to detach and hence the higher frequency. At a jet velocity of 10 m/s the agreement between the simulated and theoretical values is excellent but with decreasing velocity, the simulated values are larger except at a jet velocity of 2 m/s where the theoretical value is larger than the simulated one. 4.3. Effect of simulation time on flow variables An investigation into the effect of different time averaging periods was carried out to identify the quasi-steady state. Therefore the first simulation was run for a total of 20 seconds and numerical results of the solids volume fraction and axial particle velocity were averaged

Gidaspow

Syamlal O’Brien

Fig. 5. Experimental and simulated bubble shapes for different drag models (U0 = 10 m/s, t = 0.2 seconds).

M. Lungu et al. / Powder Technology 269 (2015) 139–152

(a)

145

(b)

1.0 0.9

αg

0.8 0.7 0.6 0.5 0.4 0.3 0.0

0.4

0.8

1.2

1.6

2.0

t (s) Fig. 6. (a) Simulated and (b) experimental porosity fluctuations at the bed centerline 3.13 cm above the jet mouth.

for three time intervals namely 5–10 seconds, 5–15 seconds and 5–20 seconds. Fig. 8 shows the time averaged solid volume fraction at three averaging periods at a jet velocity of 10 m/s. Qualitatively, the flow structure in the bed does not change appreciably though for time averaging period of 5–20 seconds the jet slightly tilts to the left. Three solid flow regions can be observed in the bed namely the jetting region in the vicinity of the central jet, the wall region of the lower half of the bed and bottom corners and the freeboard region. The jet emerging periodically from the orifice accelerates the solid particles in its vicinity and then splits up into bubbles meanwhile some of the particles move downwards in the annular region and are recirculated by the high velocity gas emerging from the orifice. Other particles are carried to the surface of the bed by the rising bubbles and upon bursting of the bubbles, the particles slide down the walls forming a dense phase of solids at the walls and at the bottom corner. Very few particles are entrained from the reactor and hence the freeboard region has little solids concentration. Such a flow structure has been observed experimentally by Yang et al. [36] who used a force probe to study the solids circulation pattern and particles mixing in a large jetting fluidized bed. Fig. 9 shows the simulated mean axial velocity for a jet velocity of 10 m/s at the three time averaging periods. Similarly very slight differences in the simulated axial mean particle velocities are observed. For the time average of 5–10 seconds the maximum positive axial velocity is 1.2 ms while for the averaging periods of 5–15 seconds and 5–20 seconds it is 1.4 m/s. The figure also reinforces the idea of the

10

Simulated value Theoretical value 9

fb (Hz)

8

7

6

5 2

4

6

8

U0 (m/s) Fig. 7. Simulated and theoretical bubble frequency versus jet velocity

10

observed flow pattern in the bed. In the jet region the mean axial velocity is positive due to the upward moving particles and in the wall region it is negative due to the particles sliding down the walls. In parts of the annular region, the axial velocity is 0 m/s indicating recirculation of particles. 4.4. Reynolds stresses The principal characteristic of turbulent flow is the production of additional stresses due to random velocity fluctuations, called Reynolds stresses [37]. In the present work we have adopted the method of Jiradilok et al. [37] to compute the normal Reynolds stresses in both the axial and lateral direction. From the CFD simulations, hydrodynamic velocities u (x, t) were monitored at nine lateral positions for the three bed heights. The normal Reynolds stress is given as: u0i u0i ¼

m 1X ðu ðx; t Þ−ui ðxÞÞðuik ðx; t Þ−ui ðxÞÞ m k¼1 ik

ð34Þ

where ui ðxÞ is the mean hydrodynamic velocity of the particle and is given by: u0i ¼

m 1X u ðx; t Þ: m k¼1 ik

ð35Þ

The hydrodynamic velocity is the ensemble average of the instantaneous particle velocity within a finite volume i.e. computational cell and time period. 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 15 seconds as aforementioned thereby giving a data set of 100,000 data points. Fig. 10 shows the computed axial and lateral normal Reynolds stresses lateral profile at three different heights. The axial normal Reynolds stresses are somewhat a reflection of the axial solids velocity, in the region near the central jet a characteristic peak can be observed at all three heights of 0.1 m, 0.3 m and 0.5 m as a result of maximum axial velocity in this region (see Fig. 9). Furthermore at a height of 0.1 m, two other peaks can be observed near the walls and these can be explained by appreciating the fact that the axial velocity near the walls is negative and Reynolds stresses are computed from the square of fluctuating velocities which ultimately give positive values resulting in the two peaks. At a height of 0.5 m the peak near the center region is slightly tilted to the left reflecting the axial particle velocity (region in red in Fig. 9c) in the bed center which also tilts to the left slightly. Meanwhile the lateral Reynolds stresses also show an increasing tendency with bed height; however no characteristic peaks are readily observable instead the lateral profiles are nearly more less flat in comparison to the axial normal Reynolds

146

M. Lungu et al. / Powder Technology 269 (2015) 139–152

Fig. 8. Time averaged solid volume fraction for time period (a) 5–10 seconds (b) 5–15 seconds (c) 5–20 seconds.

stresses. A comparison of the axial and lateral normal Reynolds stresses reveals the anisotropic characteristics of particle fluctuation. In addition the axial normal Reynolds stress in the axial direction (direction of flow) is larger than that in the lateral direction as expected. A similar trend has been observed in previous works [37,38].

clusters or bubbles. The laminar granular temperature is related to the particle normal stresses and is given as:

4.5. Granular temperature

where i represents the x, y and z directions respectively. For a 2D simulation the fluctuations in the lateral and tangential directions (non-flow directions) are assumed to be equal and thus the laminar granular temperature can be given as:

The concept of granular temperature forms the core of the kinetic theory of granular flow and is analogous to the thermodynamic temperature in the kinetic theory of gases and thus properties such as particulate pressure and viscosity are dependent on the granular temperature. Furthermore the distribution of granular temperature in a fluidized bed is an indicator of the particle mixing and inter-particle heat transfer [39] . Two kinds of granular temperature can be identified in bubbling fluidized beds; the first type is described as a laminar type due to random oscillations of individual particles measured by the classical granular temperature and the other is a turbulent (“bubble-like”) type caused by the motion of clusters of particles or bubbles, measured by the average normal Reynolds stress [40]. The two granular temperatures represent mixing at different scales; the laminar granular temperature represents mixing at the scale of individual particles whereas the turbulent granular temperature represents mixing at the scale of particle

θlaminar ¼

θlaminar ¼

1 ½hC C i 3 i i

E i 1 hD C y C y þ 2hC x C x i : 3

ð36Þ

ð37Þ

The particle normal stresses are calculated from the fluctuating instantaneous velocity in the axial and lateral directions. The fluctuating instantaneous velocity is defined as instantaneous velocity minus the hydrodynamic velocity. The laminar particle stress is defined as: p 1 X ðc ðx; t Þ−ui ðx; t ÞÞðcik ðx; t Þ−ui ðx; t ÞÞ np k¼1 ik n

hC i C i i ¼

Fig. 9. Simulated time averaged axial particle velocity for time period (a) 5–10 seconds (b) 5–15 seconds (c) 5–20 seconds

ð38Þ

M. Lungu et al. / Powder Technology 269 (2015) 139–152

(a)

(b) Lateral normal Reynolds stress (m /s )

2.0

2

1.6

0.8

0.4

-0.5

0.0

0.5

y=0.1m y=0.3m y=0.5m

2

1.2

0.0 -1.0

1.0

2

y=0.1m y=0.3m y=0.5m

2

Axial normal Reynolds stress (m /s )

147

1.0

0.8

0.6

0.4

0.2

0.0 -1.0

-0.5

0.0

0.5

1.0

x/W (-)

x/W (-)

Fig. 10. Lateral distributions of computed normal (a) axial and (b) lateral Reynolds stresses at three different heights

The laminar granular temperature is already programmed into the CFD code and calculated directly by solving the classical partial differential granular temperature equation. However the CFD simulations do not resolve the laminar granular temperature into different components unlike the PIV technique used to calculate the experimental values, instead only the overall laminar granular temperature value is computed [17]. The sum of the two granular temperatures gives the total granular temperature which describes the overall system oscillation [17]:

where p 1 X c ðx; t Þ np k¼1 ik n

ui ðx; t Þ ¼

ð39Þ

where np is the number of particles paer unit volume, k is the particle number in any given volume, x is the position, i refers to the x, y or z coordinate, t is the time interval, cik is the instantaneous velocity of particle k in a given volume, ui is the average velocity per unit volume (hydrodynamic velocity). Meanwhile turbulent granular temperature is due to the oscillation of bubbles and is calculated from the time averaged normal Reynolds stresses and given as follows: θturbulent ¼

1 0 0 uu 3 i i

θtotal ¼ θlaminar þ θturbulent :

Fig. 11 shows the spatially and time averaged laminar, turbulent and total granular temperatures along the height of the reactor. All three granular temperatures exhibit a similar profile. At the bottom of the bed the granular temperatures are low and moving up the bed the granular temperatures show a remarkable increase at a bed height of about 0.4 m and around a bed height of 0.8 m the laminar granular temperature drops to nearly zero due the presence of very little or no particles; meanwhile the turbulent and total granular temperatures also show a decrease due to the weak motion of bubbles and/or particle clusters. The profile of the total granular temperature closely resembles that of the turbulent granular temperature due to the fact that the turbulent granular temperature is much larger than the laminar granular temperature. As with previous studies [17,41], the granular temperature due to bubble motion and/or cluster oscillation is larger than that due to individual particle oscillation an indication that energy required to sustain the bubble motion or particle clusters mainly comes from the relative motion of the gas and solid phase driven by the gas flow [42]. The granular temperature profiles reflect the distribution of the solid volume

ð40Þ

where i represents the x, y or z direction. Similarly the velocity fluctuations in the non-flow directions (x and z directions) are assumed to be equal giving the turbulrnt granular temperature to be: θturbulent ¼

1 0 0 2 0 0 u u þ u u : 3 y y 3 x x

ð41Þ

The turbulent granular temperature as aforementioned is calculated indirectly from the normal Reynolds stresses which are simply the square of time averaged fluctuating velocity components in the different directions. Thus, whereas the laminar granular temperature is a function of the particle normal stresses, the turbulent granular temperature is a function of the time averaged normal Reynolds stresses.

(a)

(b)

1.0

0.4

0.2

0.8

Bed height (m)

Bed height (m)

Bed height (m)

1.0

0.8

0.6

0.0 0.000

(c)

1.0

0.8

0.6

0.4

0.2

0.002

0.004

0.006

0.008 2

2

Granular temperature (m /s )

0.010

0.0 0.0

ð42Þ

0.6

0.4

0.2

0.2

0.4

0.6

0.8 2

2

Granular temperature (m /s )

1.0

0.0 0.0

0.2

0.4

0.6

0.8 2

2

Granular temperature (m /s )

Fig. 11. Axial profiles of spatially and time averaged (a) laminar (b) turbulent and (c) total granular temperatures along the height of the reactor.

1.0

148

M. Lungu et al. / Powder Technology 269 (2015) 139–152

mean free path of the particles as a result of particles being close together and this region is called the collisional regime.

Chalermsinsuwan et al [43] Gidaspow and Mostofi [44] Miller and Gidaspow [45] Polashenski and Chen [46] Present work

1

4.6. Energy spectrum

2

2

Granular temperature (m /s )

10

0.1

0.01

1E-3 0.0

0.1

0.2

0.3

0.4

0.5

0.6

solid volume fraction Fig. 12. Comparison of simulated total granular temperature with literature results.

fraction distributions as shown in the contour plot in Fig. 8. At the bottom of the reactor there is an accumulation of solid material on both sides of the jet and near the bottom walls so the granular temperatures are low since the movement of particles is somewhat restricted. Moving up the bed the jet region expands and the solid phase is dilute in this vicinity meaning that the solid particles have an increased mean free path and therefore the individual particle oscillations increase and at the same time the bubbles grow and coalesce as they move up the bed and thus the turbulent oscillations due to bubble motion increase as well. Fig. 12 shows a comparison of the simulated total granular temperature with derived experimental and simulated data from the literature [43–46]. The agreement between the literature results and those from our simulations is fairly good especially at low solid volume fractions. The simulated total granular temperature can be seen to increase with increasing solid volume fraction passing through a maximum and then decreases. This profile can be observed from the other works presented in the diagram particularly from the works of Chalermsinsuwan et al. [43] and Gidaspow and Mostofi [44]. This can be explained using the analogy of the kinetic theory of gases. In regions of low solid volume fraction i.e. region before the maximum, the granular temperature which is essentially equal to two thirds of the kinetic turbulent energy increases with increasing volume manner in a similar manner to compression of an ideal gas. Thus this region is rightly called the kinetic regime. Meanwhile in regions of high volume fraction after the peak, a reduction in granular temperature is due to the reduced

(a)

∞ 0

Ei ðnÞdn ¼ ui 0 ui 0 :

ð43Þ

If the turbulence contains only large eddies, the distribution function, Ei (n), will exist mainly in the region of low frequencies; alternatively if there are only small eddies, Ei (n) will exist mainly in the region of high frequencies. Figs. 13 and 14 show the computed axial and lateral energy spectra for the bed center and near wall regions at three different heights respectively. From Fig. 13 we observe that the axial energy spectrum at the center of the bed is larger than the lateral energy spectrum at the same position which is consistent with the calculated Reynolds stresses and entails that most of the energy is in the vertical direction i.e. the flow direction. It is also observed from Fig. 14 that the energy spectra

x/W = 0.0

100 10 1 0.1 0.01 0.01

1000

2

1000

y=0.1m y=0.3m y=0.5m

x/W = 0.0

y=0.1m y=0.3m y=0.5m

Radial energy spectrum (m /s)

2

Z

(b)10000

10000

Axial energy spectrum (m /s)

The turbulent eddy cascade theory states that energy is transferred from large unstable eddy to smaller eddies. In turn the smaller eddies break up as well transferring energy to yet even smaller eddies thereby creating an energy cascade. The large eddies contain most of the turbulent energy as a result of external forces acting on the fluid and in addition they possess a higher Reynolds number which entails that inertial effects dominate remembering that the Reynolds number is a ratio of the inertia effects to the viscous effects. This process goes on until such that the viscous forces start to dominate and dissipation of energy occurs. Energy is lost at each hierarchy due to breakage of eddies and by dissipation into heat by molecular viscosity of the fluid which occurs inside the turbulent eddies as they break up. The aim of spectral analysis is to characterize the turbulent energy with various wave numbers which represents different scales of turbulence. In so doing the transfer of turbulent energy between different scales of turbulence as well the dissipation of turbulent energy by viscosity can be physically interpreted. Spectral analysis of turbulent flow has been previously used in our previous work [47] for the characterization of hydrodynamics in the bubbling regime for Geldart B particles as well as other workers [37,38,48] for other regimes. In the present work the energy spectra have been analyzed at the bed center i.e. x/W = 0.0 and near the wall at x/ W = 0.95 at bed heights of 0.1, 0.3 and 0.5 m respectively in order to get information about the dynamic behavior of the particles. The energy spectra can be conveniently obtained from the Fourier transform of u′u i ′i using the fast Fourier transform (FFT) technique. The summation of the energy spectrum function at all frequencies is simply equal to the normal Reynolds stress i.e.

0.1

1

Frequency (Hz)

10

100

100 10 1 0.1 0.01 0.01

0.1

1

Frequency (Hz)

Fig. 13. Simulated (a) axial and (b) lateral energy spectrum at bed center at three different heights

10

100

M. Lungu et al. / Powder Technology 269 (2015) 139–152

(a)

(b) 1000

y=0.1m y=0.3m y=0.5m

2

100 10 1 0.1 0.01 1E-3 1E-4 1E-5 0.01

0.1

1

10

x/W = 0.95

y=0.1m y=0.3m y=0.5m

100

2

x/W=0.95

Radial energy spectrum (m /s)

1000

Axial energy spectrum (m /s)

149

100

10 1 0.1 0.01 1E-3 1E-4 1E-5 0.01

Frequency (Hz)

0.1

1

10

100

Frequency (Hz)

Fig. 14. Simulated (a) axial and (b) lateral energy spectrum at near wall region for three different heights

at the near wall region are lower than that at the bed center due to wall friction which restricts intense particle fluctuation. Further both axial and lateral energy spectra in the near wall region show an increasing tendency with bed height due to lower solid volume fraction and increased oscillation, as we shall later on demonstrate that regions of high energy spectra particularly at the bed center exhibit high oscillations of the heat transfer coefficient as opposed to the region near the wall which exhibits low heat transfer coefficient fluctuations. 4.7. Local heat transfer coefficients The simulations reveal that the maximum local instantaneous fluidto-particle heat transfer coefficient occurs in the wake of the bubble. Fig. 15 shows the local heat transfer associated with the initial startup bubble at different times at a jet velocity of 10 m/s. Fig. 15 clearly shows that the maximum local heat transfer coefficient is in the wake of the bubble close to the orifice. This can be attributed to the turbulence in the localized eddies of the bubble wake [14,49]. Kuipers et al. [11] similarly observed that the maximum instantaneous wall-to-bed heat transfer coefficient was in the wake of the bubble. However the role of the bubble motion in their work is on a macro scale replenishing solid material close to the heated wall while in this work the action of the turbulence in the localized eddies enhances the heat transfer between the two phases. It can also be seen from the figure that the heat transfer coefficient in the freeboard is the least as there are few or no particles in this region to exchange heat with. In the bulk of the bed the heat

t=0.2s

t=0.25s

transfer coefficient is higher in comparison to the freeboard and the region inside the bubble due to the presence of particles. A close inspection of the bubbles in the figure shows that there are green pockets which have a slightly higher heat transfer coefficient compared to the rest of the bubble represented by the bluish color. This is as a result of the background gas which keeps the bed fluidized mostly flowing through the center of the bubble and creating small circulation zones at the sides of the bubble. A recent study of bubble heat transfer by Patil et al. [15] using DEM showed that the hot gas in these circulation zones is stagnant giving rise to localized high temperatures compared to the rest of the bubble and that the heat can only escape by conduction. Therefore the heat transfer coefficient in these circulation zones is lower compared to the rest of the bubble as revealed by the present work. Fig. 16 shows the local instantaneous heat transfer coefficients observed at different bed locations at a jet velocity of 10 m/s. From the figure, the local instantaneous heat transfer coefficient at the bed center decreases with increasing bed height. This can be explained by acknowledging the fact that the bed voidage increases with height and with an increase in voidage, a reduction in the Nusselt number is observed as will be shown later on. As revealed by the normal Reynolds stresses, the particle oscillation increases in the dilute region of the bed due to decreased particle concentration. Thus with the reduction of particle concentration, the oscillation of the heat transfer coefficient at higher height reduces. The figure shows the oscillations of the local heat transfer coefficient at a height of 0.3 m at the bed center and near wall region. Clearly the oscillations in the bed center are larger than

t=0.28s

t=0.33s

Fig. 15. Local instantaneous heat transfer coefficient in the bubble wake at different times

h(W/m2K)

150

M. Lungu et al. / Powder Technology 269 (2015) 139–152

(a)

(b)

1000

y =0.1m y =0.3m y =0.5m

x/W = 0.0 U0= 10m/s

1000

y = 0.3m U0=10m/s

x/W = 0.0 x/W = 0.95

800

h (W/m K)

2

h (W/m K)

800 600

2

600 400

400

200

200 0

0 0

2

4

6

8

10

12

14

0

2

4

t (s)

6

8

10

12

14

t (s)

Fig. 16. Local instantaneous heat transfer coefficients at (a) bed center for different bed heights and (b) a height of 0.3 m for bed center and near wall regions

those in the near wall region. The wall region has a damping effect on the particle velocity fluctuations resulting in reduced fluctuations as shown by both the axial and lateral energy spectra. The turbulent kinetic energy in the bed center is greater than that in the near wall region resulting in greater fluctuations.

of 286.70 W/m2 K at the reactor outlet. This profile is closely related to the solid volume fraction not shown here for the sake of brevity but supplied together with paper in the supplementary data.

4.8. Averaged heat transfer coefficient

Fig. 18 shows the variation of the time averaged heat transfer coefficient with jet velocity. The simulations reveal that the heat transfer coefficient increases with increasing jet velocity reaching a maximum value of 417.06 W/m2 K at a jet velocity of 8 m/s after which it starts to decrease. This plot also shows that increasing the velocity from 2 m/s to 8 m/s results in an increase of only 31 W/m2 K. However owing to the large interfacial areas in fluidized bed typically in the range of 3000–45,000 m2/m3 [50], the heat transfer per unit volume is large despite the small observed change in the fluid-to-particle heat transfer coefficient. In our previous work [51] we have studied and discussed the effects of particle size as well as initial bed height on the predicted volumetric interphase heat transfer coefficients. To get an insight into the reason for the small change in the fluid-particle heat transfer despite a large increase in the jet velocity, we plotted averaged values of the simulated Nusselt number versus the average bed voidage for the five velocities used in the simulation as shown in Fig. 19. From the scatter plot it is observed that the Nu number versus voidage profile is similar for the five velocities without any major appreciable differences. At packed bed condition, the Nu number is high and with increasing voidage the Nu decreases exponentially approaching the minimum theoretical value of 2 as the voidage approaches 1. The Nu number does

Fig. 17 shows the time averaged lateral and axial distribution of the heat transfer coefficient at a jet velocity of 10 m/s. From Fig. 17a the time averaged heat transfer coefficient increases with decreasing height owing to the increase of particle concentration as aforementioned. In addition the lateral profile reveals that the average heat transfer coefficient is highest near the wall region despite the observed high fluctuations in the bed center especially at heights of 0.3 m and 0.5 m. The solid volume fraction is high in the near wall region as shown in Fig. 8 and this explains the high heat transfer coefficients observed in this region. Also the low energy spectrum in the near wall region suggests an accumulation of particles leading to the observed high heat transfer coefficients. Therefore high heat transfer coefficient observations should be an indicator of hot spot formation as in this case. Fig. 17b shows the profile of the time and spatially averaged heat transfer coefficient with height. The profile is not smooth but consists of peaks and toughs with the maximum peak occurring at a height of 0.12 m giving an averaged heat transfer coefficient value of 542.95 W/m2 K. At a bed height of about 0.46 m the heat transfer coefficient exhibits a sharp fall and then finally giving a minimum heat transfer coefficient value

(a)

(b) 1.0

600 y=0.1m y=0.3m y=0.5m

580 560 540

0.8

Bed height (m)

2

Heat transfer coefficient (W/m K)

4.9. Effect of jet velocity on predicted heat transfer coefficient

520 500 480 460 440

0.6

0.4

420 0.2

400 380 360 -1.00

0.0 -0.75

-0.50

-0.25

0.00

x/W (-)

0.25

0.50

0.75

1.00

200

300

400

500

600

2

Heat transfer coefficient (W/m K)

Fig. 17. Computed (a) lateral average heat transfer coefficient profile at three different heights and (b) time and spatially averaged heat transfer coefficient profile along the height of the reactor

M. Lungu et al. / Powder Technology 269 (2015) 139–152 430

instantaneous fluid-to-particle heat transfer coefficient occurs in the wake of the bubble. Local instantaneous heat transfer coefficients in the center of the bed display greater oscillations compared to those in the near wall region. The large fluid-to-particle heat transfer coefficient oscillations in the center of the bed are due to the bubbles and individual particle oscillations which also exhibit high fluctuations at the bed center region as revealed by the granular temperature profiles and energy spectra. In addition heat transfer coefficients have been observed to decrease with increasing height as a result of the decreasing particle concentration. Maximum averaged heat transfer coefficients have been observed to occur near the wall due to the presence of an accumulation of the solid phase. The average heat transfer coefficient increases with increasing jet velocity reaching a maximum value at a jet velocity of 8 m/s after which it starts to decrease. An increase in the jet velocity from 2 m/s to 8 m/s results in an increase of the fluid-to-particle heat transfer coefficient by only 31 W/m2 K. The plot of the variation of the Nusselt number with voidage shows similar profiles at different jet velocities which explain the small change in the heat transfer coefficient.

420

2

h (W/m K)

410

400

390

380

370 2

4

6

8

10

U0 (m/s) Fig. 18. Effect of jet velocity on the average heat transfer coefficient.

not vary much for the five velocities and this explains the small change of the heat transfer coefficient across the five velocities. 5. Conclusions An Eulerian–Eulerian CFD model has been applied to the study of the flow structure and heat transfer in a bubbling fluidized bed operated with a central jet. The model has been validated using experimental bubble properties obtained from literature. The Syamlal O'Brien drag model has been found to be the most appropriate in the prediction of bubble properties. Time averaged solid volume fraction contour plots show the existence of three solid flow regions namely the jetting region in the vicinity of the central jet, the wall region of the lower half of the bed and bottom corners and the freeboard region. These solid flow regions are also reflected by the granular temperature profiles. The granular temperature profile shows an increase in dilute regions particularly in the jetting regions to the increased particle and bubble oscillations and in the dense regions and the freeboard region the granular temperature is low. Meanwhile the normal Reynolds stresses exhibit anisotropy with the axial normal Reynolds stresses being greater than the lateral normal Reynolds stresses, an indication that the flow is predominantly in the vertical direction. The anisotropy in the flow is also revealed by the energy spectra which show that the axial energy spectra are larger than the lateral energy spectra. Simulations also show that the maximum

4.0

Nomenclature CD particle drag force coefficient [−] cp specific heat capacity at constant pressure [J/kg K] ds mean particle diameter [m] g acceleration due to gravity [m/s2] h fluid-particle heat transfer coefficient [W/m2 K] H enthalpy [J/kg] Nu Nusselt number [−] p pressure [pa] Pr Prandtl number Res particle Reynolds number [−] t time [s] U0 jet velocity [m/s] y vertical coordinate [m]

Greek letters α volume fraction [−] βgs gas/solid momentum exchange coefficient [kg/m2s] Θs granular temperature [m2/s2] θturbulent turbulent granular temperature [m2/s2] ρ density [kg/m3] τk stress tensor of phase k [pa] κ thermal conductivity [W/m K]

Subscripts and superscripts g gas k phase k, solid or gas s solid phase

10 m/s 8 m/s 6 m/s 4 m/s 2 m/s

3.5

151

Nu

3.0

Acknowledgements The authors wish to express their sincere gratitude for the financial support provided by The Project of National Natural Science Foundation of China (21176207, 21236007) and the doctoral program of Higher Education Research Fund (20130101110063).

2.5

2.0 0.4

0.5

0.6

0.7

αg

0.8

0.9

Fig. 19. Dependence of the Nusselt number on voidage.

1.0

Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.powtec.2014.08.067.

152

M. Lungu et al. / Powder Technology 269 (2015) 139–152

References [1] V.V. Ranade, Computational Flow Modeling for Chemical Reactor Engineering, Academic Pr, 2002. [2] J. Nieuwland, M. Veenendaal, J. Kuipers, W. Van Swaaij, Bubble formation at a single orifice in gas-fluidised beds, Chem. Eng. Sci. 51 (1996) 4087–4102. [3] J. Kuipers, W. Prins, W. Van Swaaij, Theoretical and experimental bubble formation at a single orifice in a two-dimensional gas-fluidized bed, Chem. Eng. Sci. 46 (1991) 2881–2894. [4] D. Harrison, L. Leung, Bubble Formation at an Orifice in a Fluidized Bed, 1961. [5] X. Nguyen, L. Leung, A note on bubble formation at an orifice in a fluidised bed, Chem. Eng. Sci. 27 (1972) 1748–1750. [6] W.-C. Yang, D. Revay, R. Anderson, E. Chelen, D. Keairns, D. Cicero, Fluidization phenomena in a large-scale, cold-flow model, Fourth International Conference on Fluidization (Engineering Foundation, 1983), 1984. [7] H.S. Caram, K.-K. Hsu, Bubble formation and gas leakage in fluidized beds, Chem. Eng. Sci. 41 (1986) 1445–1453. [8] G. Bokkers, M. van Sint Annaland, J. Kuipers, Mixing and segregation in a bidisperse gas–solid fluidised bed: a numerical and experimental study, Powder Technol. 140 (2004) 176–186. [9] D. Patil, M. van Sint Annaland, J. Kuipers, Critical comparison of hydrodynamic models for gas–solid fluidized beds—part I: bubbling gas–solid fluidized beds operated with a jet, Chem. Eng. Sci. 60 (2005) 57–72. [10] A. Kumar, S. Das, D. Fabijanic, W. Gao, P. Hodgson, Bubble–wall interaction for asymmetric injection of jets in solid–gas fluidized bed, Chem. Eng. Sci. 101 (2013) 56–68. [11] J. Kuipers, W. Prins, W. Swaaij, Numerical calculation of wall-to-bed heat-transfer coefficients in gas-fluidized beds, AIChE J 38 (1992) 1079–1091. [12] D.J. Patil, J. Smit, M. van Sint Annaland, J.A.M. Kuipers, Wall-to-bed heat transfer in gas–solid bubbling fluidized beds, AIChE J 52 (2006) 58–74. [13] R. Yusuf, B. Halvorsen, M.C. Melaaen, An experimental and computational study of wall to bed heat transfer in a bubbling gas–solid fluidized bed, Int. J. Multiphase Flow 42 (2012) 9–23. [14] D. Kunii, O. Levenspiel, Fluidization Engineering, Butterworth-Heinemann, Boston, 1991. [15] A.V. Patil, E.A.J.F. Peters, T. Kolkman, J.A.M. Kuipers, Modeling bubble heat transfer in gas–solid fluidized beds using DEM, Chem. Eng. Sci. 105 (2014) 121–131. [16] J.R. Grace, T. Li, Complementarity of CFD, experimentation and reactor models for solving challenging fluidization problems, Particuology 8 (2010) 498–500. [17] M. Kashyap, D. Gidaspow, W.J. Koves, Circulation of Geldart D type particles: part I — high solids fluxes. Measurements and computation under solids slugging conditions, Chem. Eng. Sci. 66 (2011) 183–206. [18] H.S. Caram, K.-K. Hsu, Bubble formation and gas leakage in fluidized beds.pdf, Chem. Eng. Sci. 41 (1986) 1445–1453. [19] M. Syamlal, T.J. O'Brien, Computer simulation of fluidized beds, AIChE Symp. Ser. 85 (1989) 22–31. [20] D. Gidaspow, Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions, Academic Press, San Diego, 1994. [21] L. Huilin, D. Gidaspow, Hydrodynamics of binary fluidization in a riser: CFD simulation using two granular temperatures, Chem. Eng. Sci. 58 (2003) 3777–3792. [22] C. Lun, S. Savage, D. Jeffrey, N. Chepurniy, Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in a general flowfield, J. Fluid Mech. 140 (1984) 223–256. [23] D.G. Schaeffer, Instability in the evolution equations describing incompressible granular flow, J. Differ. Equ. 66 (1987) 19–50. [24] D. Ma, G. Ahmadi, An equation of state for dense rigid sphere gases, J. Chem. Phys. 84 (1986) 3449–3450. [25] D. Gunn, Transfer of heat or mass to particles in fixed and fluidised beds, Int. J. Heat Mass Transf. 21 (1978) 467–476. [26] W. Ranz, W. Marshall, Evaporation from drops, Chem. Eng. Prog. 48 (1952) 141–146. [27] N. Wakao, S. Kaguei, T. Funazkri, Effect of fluid dispersion coefficients on particle-tofluid heat transfer coefficients in packed beds: correlation of Nusselt numbers, Chem. Eng. Sci. 34 (1979) 325–336.

[28] A. Cybulski, M.J. Van Dalen, J.W. Verkerk, P.J. Van Den Berg, Gas-particle heat transfer coefficients in packed beds at low Reynolds numbers, Chem. Eng. Sci. 30 (1975) 1015–1018. [29] P.A. Nelson, T.R. Galloway, Particle-to-fluid heat and mass transfer in dense systems of fine particles, Chem. Eng. Sci. 30 (1975) 1–6. [30] N.G. Deen, S.H.L. Kriebitzsch, M.A. Van der Hoe, J.A.M. Kuipers, Direct numerical simulation of flow and heat transfer in dense fluid particle systems, Chem. Eng. Sci. 81 (2012) 329–344. [31] O. Olaofe, M. Van der Hoef, J. Kuipers, Bubble formation at a single orifice in a 2D gas-fluidized bed, Chem. Eng. Sci. 66 (2011) 2764–2773. [32] J.R. Grace, Reflections on turbulent fluidization and dense suspension upflow, Powder Technol. 113 (2000) 242–248. [33] J.A.M. Kuipers, H. Tammes, W. Prins, W.P.M. van Swaaij, Experimental and theoretical porosity profiles in a two dimensional gas fluidized bed with a central jet, Powder Technol. 71 (1992) 87–99. [34] N. Xie, F. Battaglia, S. Pannala, Effects of using two- versus three-dimensional computational modeling of fluidized beds, Powder Technol. 182 (2008) 1–13. [35] T.W. Asegehegn, M. Schreiber, H.J. Krautz, Influence of two-and three-dimensional simulations on bubble behavior in gas–solid fluidized beds with and without immersed horizontal tubes, Powder Technol. (2011). [36] W.-C. Yang, B. Ettehadieh, G.B. Haldipur, Solids circulation pattern and particles mixing in large jetting fluidized bed, AIChE J 32 (1986) 1994–2000. [37] V. Jiradilok, D. Gidaspow, S. Damronglerd, W.J. Koves, R. Mostofi, Kinetic theory based CFD simulation of turbulent fluidization of FCC particles in a riser, Chem. Eng. Sci. 61 (2006) 5544–5559. [38] B. Chalermsinsuwan, P. Piumsomboon, D. Gidaspow, Kinetic theory based computation of PSRI riser: part I—estimate of mass transfer coefficient, Chem. Eng. Sci. 64 (2009) 1195–1211. [39] S. Deb, D. Tafti, Investigation of flat bottomed spouted bed with multiple jets using DEM–CFD framework, Powder Technol. 254 (2014) 387–402. [40] J. Jung, D. Gidaspow, I.K. Gamwo, Measurement of two kinds of granular temperatures, stresses, and dispersion in bubbling beds, Ind. Eng. Chem. Res. 44 (2005) 1329–1341. [41] J. Jung, D. Gidaspow, I.K. Gamwo, Bubble computation, granular temperatures, and Reynolds stresses, Chem. Eng. Commun. 193 (2006) 946–975. [42] W. Shuai, Z. Guangbo, L. Guodong, L. Huilin, Z. Feixiang, Z. Tianyu, Hydrodynamics of gas–solid risers using cluster structure-dependent drag model, Powder Technol. 254 (2014) 214–227. [43] B. Chalermsinsuwan, D. Gidaspow, P. Piumsomboon, 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. 171 (2011) 301–313. [44] D. Gidaspow, R. Mostofi, Maximum carrying capacity and granular temperature of A, B and C particles, AIChE J. 49 (2003) 831–843. [45] A. Miller, D. Gidaspow, Dense, vertical gas–solid flow in a pipe, AIChE J. 38 (1992) 1801–1815. [46] W. Polashenski, J.C. Chen, Measurement of particle phase stresses in fast fluidized beds, Ind. Eng. Chem. Res. 38 (1999) 705–713. [47] J. Sun, Y. Zhou, C. Ren, J. Wang, Y. Yang, CFD simulation and experiments of dynamic parameters in gas–solid fluidized bed, Chem. Eng. Sci. 66 (2011) 4972–4982. [48] B. Chalermsinsuwan, D. Gidaspow, P. Piumsomboon, In-depth system parameters of transition flow pattern between turbulent and fast fluidization regimes in high solid particle density circulating fluidized bed reactor, Powder Technol. (2013). [49] W.-C. Yang, Handbook of Fluidization and Fluid-Particle Systems, Taylor and Francis Group LLC, United States of America, 2003. [50] J.S.M. Botterill, Fluid-Bed Heat Transfer: Gas-Fluidized Bed Behaviour and Its Influence on Bed Thermal Properties, Academic Press, London, 1975. [51] M. Lungu, J. Sun, J. Wang, Z. Zhu, Y. Yang, CFD simulations of interphase heat transfer in a bubbling fluidized bed, Korean J. Chem. Eng. 31 (2014) 1148–1161.