Discrete Element Study of Solid Mixing Behavior with Temperature

Article pubs.acs.org/IECR

Discrete Element Study of Solid Mixing Behavior with Temperature Difference in Three-Dimensional Bubbling Fluidized Bed Shiliang Yang, Kun Luo,* Mingming Fang, Jianren Fan, and Kefa Cen State Key Laboratory of Clean Energy Utilization, Zhejiang University, Hangzhou 310027, People’s Republic of China ABSTRACT: This work numerically evaluates the mixing behavior of a binary mixture with the same physical properties but different temperatures in a three-dimensional bubbling fluidized bed based on the coupling of computational fluid dynamics and discrete element method. General mixing procedure and mixing index are adopted to investigate the macroscopic mixing behavior of the solid phase. Moreover, the effects of superficial velocity and the initial temperature configuration of particle groups on the temperature evolution of the system are discussed. The results show that three mixing mechanisms induced by the bubbles can be identified. The time required for a system to reach the temperature dynamic equilibrium is different from that required for the spatial mixing. With an increase in the superficial velocity, the mixing procedure is enhanced, and the temperature of each particle group distributes more uniformly. In addition, the effect of temperature configurations of different particle groups appears at the initial stage of solid mixing. This initial effect weakens slowly and even disappears with time evolution.

1. INTRODUCTION Solid mixing widely exists in the applications of technology and nature, especially commonly encountered in many chemical processes due to its easy operation. In the fluidized bed, particles interact vigorously with the fluidizing medium due to the momentum, energy, and even composition exchanges. If the mixing procedure is nonuniform or not well controlled in some harsh environments, the agglomerates and catalyst deactivation are formed in the chemical reactors, which strongly influence the operation performance of the apparatus. Thus, the prediction of solid mixing performances under different operating conditions is extremely necessary for the design and operation of a dense fluidized bed. Plenty of experiments have been carried out to investigate the mixing behavior of solid phase in the gas−solid fluidized bed.1−3 It is well-known that difficulties are involved in experimentally determining the mixing behavior of particles, especially quantitatively evaluating the mixing intensity of the solid phase in the dense fluidized bed. With the development of computation capacity, numerically investigating the solid mixing has been carried out for the purpose of obtaining high-quality products of chemical processes in past years. As we have known, there are mainly two approaches used for modeling the dense gas−solid flows, namely, the Eulerian−-Eulerian method and the Eulerian− Lagrangian method. In the framework of the Eulerian−Eulerian approach, both the fluid and solid phases are treated as continuous mediums. Cooper and Coronella4 conducted 2-D numerical simulation based on a two-fluid model (TFM) to study the particle mixing in a binary of a fluidized bed. Their results demonstrated that the solid wake trailing and mixing mechanisms can be identified from the vertical solid flux. Chang et al.5 recently studied the particle−particle heat transfer in the binary mixture of a fluidized bed using the TFM approach. They showed that the heat transfer coefficient between different particle groups increases with enlarging of the larger particle size or superficial velocity. Liu et al.6 investigated the © 2014 American Chemical Society

gas and solid mixing behavior in the catalytic cracking strippers based on the TFM. On the other hand, the Eulerian−Lagrangian method, namely, the computational fluid dynamics combined with discrete element (CFD-DEM) coupling approach, has been frequently adopted to study gas−solid hydrodynamics due to its advantage of tracking the trajectory of each particle.7−10 Compared with TMF, dynamical information on particles such as velocities, positions, and collision information can be easily obtained in the tracking procedure. The CFD-DEM method is attractive in exploring the particle mixing behavior especially when the particle−particle or particle−wall conductive heat transfer is taken into account. By means of this approach, plenty of works on solid mixing have been carried out. Feng et al.11 investigated the segregation and mixing of a binary mixture in a fluidized bed with discrete particle simulation. Deen et al.12 discussed various approaches that are available to give quantitative information on the solid mixing state in the granular systems based on the discrete element method in a pressurized fluidized bed. The results showed that the mixing quality improves with the operating pressure due to the increased porosity in the dense phase flow. Gui and Fan13 conducted a numerical study on the particle mixing based on the fractal and entropy analysis in the bubbling fluidized bed. They proposed a dimensional function to study the microscopic characteristic of the mixing interface. Nevertheless, almost all of these reports on the solid mixing are carried out in the cold model. In the particle handling procedure, one of the important operation parameters influencing the product quality of the chemical process is the temperature. In some cases, a high-quality product requires that attention be Received: Revised: Accepted: Published: 7043

November 13, 2013 February 14, 2014 April 2, 2014 April 2, 2014 dx.doi.org/10.1021/ie4038399 | Ind. Eng. Chem. Res. 2014, 53, 7043−7055

Industrial & Engineering Chemistry Research

Article

εf is gas voidage calculated as

paid not only on the mixing time but also on the temperature evolution tendency of the whole system. If the mixing procedure is not controlled properly, some harsh environment appears such as hot formation with agglomeration and cohesion, which are not beneficial to the outcome percentage of the desired production. However, limited literature focusing on the solid mixing behavior with temperature difference can be found. Thus, the current work analyzes the spatial mixing behavior of solid phase in a 3-D bubbling fluidized bed and gains insight into the heat transfer and mixing behaviors of particle groups with different temperatures using the sophisticated CFD-DEM coupling approach. The gas motion is resolved in the framework of CFD, while solid motion is tracked with a discrete element method. The hydrodynamics of gas and solid phases of this apparatus operated in the cold model have been investigated in our previous work.14 In addition, the dispersion and circulation patterns of solid phase are explored.15 This work focuses on the mixing behavior of the solid phase with the temperature considered in the system. The CFD-DEM coupling approach is extended with consideration of the heat transfer between gas and solid phases to investigate the mixing properties of particle groups with different temperatures. The influence of bubble movement on the particle mixing is studied, followed by the quantitative and qualitative descriptions of both the spatial and temperature mixing procedures of the solid phase. Furthermore, the effects of superficial velocity on solid mixing and the temperature distribution of the solid phase are discussed. Finally, the influence of initial temperature configurations of particle groups on the heat transfer behavior of the system is analyzed.

n

εf = 1 −

∂t

+

∂(ρf εf u f,̃ i) ∂t

∂(ρf εf u f,̃ i) ∂xi

+

⎛ ∂u ∂u f, j ⎞ ∂u f, k f, i ⎟⎟ − 2 μf τf,̃ ij = μf ⎜⎜ + δij ∂xi ⎠ 3 ∂xk ⎝ ∂xj

⎛ ⎞ 1 τf,̃ ij SGS = ρf ⎜2νf, t Sf,̃ ij + τf,̃ kkδij⎟ ⎝ ⎠ 3

∂t

+

∂xi

+

(5)

(6)

where vf,t = (CsΔ)2 (2S̃f,ijS̃f,ij)1/2 and S̃ij = (1/2)((∂ũi/∂xj) + (∂ũj/∂xi)) are the deformation tensor of the filtered field and the eddy viscosity coefficient, respectively. Cs is the Smagorinsky constant coefficient, which is given by Lilly17 as Cs =

1⎛ 2 ⎞ ⎜ ⎟ π ⎝ 3C k ⎠

3/4

(7)

where Ck = 1.4 is the Kolmogorov constant. Δ (=(ΔxΔyΔz)1/3) is the characteristic length scale used for filtering variables, which is equivalent to the mesh spacing. In order to capture the computational details near the wall region, the wall function approach is used to calculate the instantaneous velocity at the first node near the wall. The wall function of a three-layer log-law with the wall-shear velocity proposed by Breuer and Rodi18 is adopted in the current study. Sp and Q in eqs 2 and 3 are the volumetric momentum source from solid motion and the volumetric heat flux source for the current computational cell, respectively. They are calculated as

(1)

n

∂xj ∂pf̃

(4)

where μf is the fluid viscosity coefficient and δij is the Kronecker function. The subgrid stress tensor τ̃f,ijSGS (=ρf(ũiũj − uiu͠ j )) in eq 2 is modeled based on the eddy viscosity concept, for which the well-known Smagorinsky model16 is adopted in the calculation. The equation for the subgrid stress tensor is formulated as

∂(ρf εf u f,̃ iu f,̃ j)

= ρf εf g − εf ∂(ρf εf Cp ,f Tf )

=0

ΔV

where n, ΔV, and Vp,i are the total number of particles locating in the current computational cell, the volume of the current cell, and the volume of particle i in the current cell, respectively. τ̃f,ij in the momentum equation is the filtered stress tensor which is calculated as

2. GOVERNING EQUATIONS 2.1. Governing Equations for Gas Motion. Gas phase in the bubbling fluidized bed is treated as the continuous medium and governed by the Navier−Stokes equations considering the interphase exchanging terms with solid phase. The conservation equations of the mass, momentum, and energy can be formulated as ∂(ρf εf )

∑i = 1 Vp, i

∂(εf (τf,̃ ij + τf,̃ ij SGS)) ∂xj

Sp = −∑ Fd, i /ΔV i=1

+ Sp

n

(2)

Q = (∑ Q f, i + Q f,w )/ΔV i=1

∂(ρf εf u f,̃ jCp ,f Tf )

(9)

where Fd,i is the drag force exerted on particle i by the gas phase. Qf,i and Qf,w are the heat transfer from fluid to particle and fluid to wall, respectively. The calculation details of Fd,i, Qf,i, and Qf,w will be discussed later. 2.2. Governing Equations for Solid Motion. In the modeling process, the discrete element method taking the temperature into account is adopted to investigate the solid motion. Newton’s second law is used for tracking the translational and rotational motions of each individual particle, and the energy balance equation is adopted to predict the particle’s temperature.

∂xj

⎡⎛ Cp ,f νf, t ⎞ ∂Tf ⎞⎤ ∂ ⎢⎜ ⎛ ⎟⎥ − Q ⎜ ⎟⎟ = εf ⎜k f + Prf, t ⎠ ∂xj ⎟⎠⎥⎦ ∂xj ⎢⎣⎜⎝ ⎝

(8)

(3)

where ũf,i, ρf, Tf, p̃f, and Cp,f are the velocity component, the density, the temperature, the pressure ,and the specific heat capacity of the fluid phase, respectively. kf and Prf,t are the thermal diffusivity and turbulent Prandtl number, respectively. g is the gravitational acceleration. 7044

dx.doi.org/10.1021/ie4038399 | Ind. Eng. Chem. Res. 2014, 53, 7043−7055


Article

Particles exchange heat flux with the surrounding environment due to the temperature difference between them. The existing heat transfer patterns for a particle in the dense fluidized bed can be categorized as three mechanisms, namely, the conductive heat transfer between the colliding particle− particle/particle−wall pair, the convective heat transfer between the particle and fluid phase, and the radioactive heat transfer between the particle and its surroundings. Radioactive heat transfer is ignored in the current research since its effect is relatively explicit only when the temperature is higher than 600 K, whereas the maximum temperature involved in the current work is only 400 K. The governing equations for particle i are formulated as mi

Ii

dvi = Fd, i + Fp, i − dt

dωi = dt

miCp , i

F3(εp) = 0.0673 + 0.212εp +

Rep =

(17)

Here, uf and dp are the velocity vector of gas phase in the current cell and the particle diameter, respectively. 2.3.2. Collision Dynamic Model. When a particle collides with its neighboring particles or the bed walls, the contact force exists due to the deformation produced in the contact procedure. Two calculation models are available to evaluate the

k

∑ Fc,ij + mi g (10)

j=1

(11)

j=1 k

= Q i ,f +

∑ Q i ,j + Q i ,wall (12)

j=1

(u f − vp)

Figure 1. Geometrical configuration of a 3-D bubbling fluidized bed.

Table 1. Physical and Numerical Parameters for Simulation gas phase inlet temperature, K viscosity, Pa·s superficial gas velocity, m/s pressure, atm density, kg/m3 thermal conductivity, W/(m·k) specific heat capacity, J/(kg·k) particles number diameter, m density, kg/m3 Young’s modulus, Pa Poisson’s ratio restitution coefficient friction coefficient thermal conductivity, W/(m·k) specific heat capacity, J/(kg·k) initial temperature, K geometry width, mm height, mm thickness, mm cells in x-direction cells in y-direction cells in z-direction

(13)

where εp (=1 − εf) is the volume fraction of the solid phase in the current cell. β is the interphase momentum exchanging coefficient, for which the correlation proposed by Koch and Hill19 is adopted in this work. It is formulated as β=

μf

k

dt

εp

εf ρf |u f − vp|d p

∑ Mij

dTp, i

Vpβ

(16)

where Rep is the Reynolds number of the particle and can be calculated as

where mi, vi, ωi, Fd,i, Fp,i, Fc,ij, Ii, Mij, Cp,i, Tp,i, and k represent the particle mass, the translational velocity, the rotational velocity, the drag force exerted by the fluid phase, the pressure gradient force, the collision force between the collision pair of particle i and particle j or the particle with the wall, the moment inertia of the particle, the torque, the specific heat capacity, the temperature of particle i, and the contact number of the neighboring particles with the current one, respectively. Qi,f, Qi,j, and Qi,wall are the convective heat flux between the particle and fluid phase, the conductive heat flux through the contact area of colliding particles, and the conductive heat flux between the particle and the wall when a particle collides with the wall, respectively. 2.3. Submodels for Two-Phase Coupling. 2.3.1. Drag Force Model. The drag force is one of the most important interaction forces between a particle and fluid phase due to the velocity difference between them. It is explicitly affected by the local concentration of the solid phase and can be formulated as Fd =

0.0232 εf 5

18μf εf 2εp ⎛ ⎞ 1 ⎜ F (ε ) + F3(εp)Rep⎟ 0 p ⎠ 2 d p2 ⎝

⎧ εp 135 ⎪ 1 + 3 2 + 64 εp ln(εp) + 16.14εp , ⎪ ⎪ 1 + 0.681εp − 8.48εp2 + 8.16εp3 F0(εp) = ⎨ ⎪ 10ε p ⎪ ⎪ ε 3 , εp ≥ 0.4 ⎩ f

(14)

εp < 0.4

(15) 7045

298 1.511 × 10−6Tf3/2(Tf + 120) 0.6−1.2 1 PM/RTf 7.76 × 10−5Tf + 2.873 × 10−3 1.232 × 10−2Tf + 1002.74 42 970 8 × 10−4 2490 5.5 × 108 0.25 0.97 0.1 1.1 840 298, 400 60 600 6.4 20 8 200



Article

phase, which is calculated from the correlation proposed by Li and Mason.27 The equation is formulated as

interaction force between the colliding pair, namely, the softsphere contact model proposed by Cundall and Strack20 and the hard-sphere contact model proposed by Hoomans et al.21 The main advantage of the soft-sphere contact model as compared with the other one lies in the fact that it can be used to describe the multiple collisions of a particle with its neighboring particles at the same time. Due to the contact-dominated flow characteristics of the solid phase in the fluidized bed, the softsphere contact model is adopted in the current research to capture the collision between particles. The interaction force between particle i and particle j can be evaluated as Fc, ij = (kn , ijδn , ij n − γn , ij vn , ij) + (kt , ijδt , ij t − γt , ij vt , ij)

hi ,f = Nu p, iλf /di

Nu p, i

(18)

⎧2 ⎪ ⎪ ⎪2 =⎨ ⎪ ⎪ ⎪2 ⎩

+ 0.6εf nRep, i1/2Pr1/3 ,

(22)

Rep, i < 200

+ 0.5εf nRep, i1/2Pr1/3 + 0.02εf nRep, i 0.8Pr1/3 , 200 < Rep, i < 1500 + 0.000045εf nRep, i1.8 ,

Rep, i > 1500 (23)

where δn,ij and δt,ij are the relative displacements of colliding particles in the normal and tangential directions, respectively. The subscripts n and t stand for the variable components in the normal and tangential directions, respectively. n and t represent the unit vectors in normal and tangential directions of a colliding pair, respectively. k and γ are the stiffness coefficient and the damping coefficient, respectively. Detailed calculation equations for k and γ can be found in the previous work.8,22 2.3.3. Conductive Heat Transfer Model. Conductive heat transfer occurs through the contact area of colliding particles. Two mechanisms evaluating the conductive heat flux are considered in the present work. One is the static contact between the colliding pair with zero relative velocity, in which the conductive heat flux is obtained using the equation proposed by Batchelor and O’Brien23 and then modified by Cheng et al.24 The other conductive mechanism for the colliding particle pairs with a nonzero relative velocity. The conductive heat flux is determined by Sun and Chen25 and then modified by Zhou et al.26 For the first way of heat transfer, the conductive heat flux between a colliding pair with zero relative velocity can be obtained from the equations formulated as Q i,j =

4k p, ik p, j k p, i + k p, j

(Acontact, i − j )1/2 (Tp, j − Tp, i)

Figure 2. Snapshot of solid distribution of the packed bed formed in the investigated bubbling fluidized bed used for model validation. (19)

For the second way of heat transfer, the conductive heat flux between the particle pair with nonzero relative velocity can be calculated as Q i,j =

cAcontact, i − j (Tp, j − Tp, i)tc

Table 2. Physical and Numerical Parameters Adopted To Validate the Proposed Model simulation domain width (W) × depth (D) × height (H), m3 grid: W × D × H grid number solid property particle number static bed height, m particle density, kg/m3 particle diameter, mm Young’s modulus, MPa Poisson ratio restitution coefficient friction coefficient time step of solid motion, s gas property gas density, kg/m3 gas viscosity, kg/(ms) gas velocity, m/s

−1/2

(ρp, i Cp , ik p, i)−1/2 + (ρp, j Cp , jk p, j)−1/2

(20)

where kp and ρp are the thermal conductivity and density of the particle, respectively. Acontact,i−j is the contact area between the particle−particle or particle−wall pair, which can be obtained from the deformation information in the collision pair. tc is the particle−particle or particle−wall collision time. c is a correction coefficient which was described in Zhou et al.26 2.3.4. Convective Heat Transfer Model. Convective heat transfer appears due to the temperature difference between a particle and the continuous phase. Mathematically, the quantity of convective heat flux between a particle and the gas phase can be formulated as Q i ,f = hi ,f Ai (Tf, i − Tp, i)

(21)

where Ai and Tp,i are the overall surface area and the temperature of the particle, respectively. Tf,i is the local temperature of current cell that the particle locates in. hi,f is the convective heat transfer coefficient between a particle and the continuous

minimum fluidizing velocity of experiment data,28 m/s time step of gas motion, s 7046

0.044 × 0.01 × 0.2 15 × 3 × 68 3060 9240 0.03 1000 1.2 0.12 0.33 0.97 0.1 1 × 10−5 1.225 1.8 × 10−5 0−0.4 m/s (for testing minimum fluidizing velocity); 0.9 m/s (for investigating the time-averaged property) 0.3 1 × 10−4



Article

Figure 3. Comparison of the simulated results and the experimental data in the bubbling fluidized bed: (a) pressure drop profiles in the systems with increasing or decreasing superficial velocity; (b) transverse distribution of solid vertical velocity at the height of Z = 0.01 m, Ug = 0.9 m/s; (c) transverse distribution of time-averaged voidage at the the height of Z = 0.0164 m, Ug = 0.9 m/s; (d) transverse distribution of the solid vertical velocity at the height of Z = 0.02 m, Ug = 0.9 m/s.

a diameter of 0.8 mm and density of 2490 kg/m3. Details of the physical parameters of gas−solid phases are listed in Table 1. A uniform structured grid is used for the discretization of the calculation domain with its width, depth, and height divided into 20 × 8 × 200 elements, respectively. The grid independent study and the influence of bed thickness on the simulation results have been performed in our previous work.14 For the discretization of the governing equations of fluid phase, the Crank−Nicolson scheme is adopted for the time discretization of the unsteady items; the central difference scheme is applied for the diffusion terms. Convective items of the equations are discretized with the first-order up-wind scheme. To track the position and temperature of each individual particle, Newton’s governing equations are integrated explicitly with fourth-order Runge−Kutta algorithm. 3.3. Boundary and Initial Conditions. In order to reduce the computational task, the cyclical boundary condition is assigned to the front and back walls of the bed for all of the parameters involved in both the CFD and DEM calculations. For the boundary condition of velocity, uniform distribution of the gas phase is chosen for the bed inlet, while the left and right walls of the bed are assigned with a nonslip boundary condition. For the boundary condition of pressure, the fixed value of ambient atmosphere is assigned to the bed outlet. The zerogradient boundary condition is used for the inlet and the left and right walls. The initial temperature of gas introduced into the bed is fixed with 298 K, while other boundaries are assigned with zero-gradient boundary condition. For the solid phase, a total number of 42970 particles are generated randomly in the bed and then fall to accumulate on the bed bottom under the effect of gravity. A packed bed with a height of 0.05 m and voidage of 0.42 is formed until the

where n is the model parameter, which is chosen to be 3.5 in the current work. Nup, λf, and Pr are the Nusselt number, the thermal conductivity of the fluid, and the Prandtl number, respectively. As can be captured from the relation of hi,f, the local voidage of the fluid phase is introduced into the correlation to take account of the influence of other particles in dense gas−solid flow.

3. NUMERICAL DETAILS 3.1. Simulation Description. The current work extends the CFD-DEM coupling with the temperature being taken into account, and the whole calculation process can be described generally. After the start of the calculation, the cell of the Eulerian grid that a particle locates in is obtained with its position. Then, the interphase exchanging terms including the drag force and the heat flux exerted on the particle can be obtained. Subsequently, collision events for all of the particles of the whole system are detected, followed by quantifying the collision forces and the conductive heat fluxes between the colliding pair. Afterward, explicit integration of the governing equations of solid motion is applied in the DEM framework, followed by the information on the position, velocity, and temperature of each particle updated. Meanwhile, the CFD solver fetches the DEM data to calculate the newest information on the voidage, the momentum, and energy sources for the governing equations of the fluid phase. After solving the partial differential equations describing the gas motion, the velocity and the temperature of gas phase are fetched again by the DEM solver to update the solid motion for the calculation of the next time step. 3.2. Computation Conditions. As shown in Figure 1, the bubbling fluidized bed has a cross-section of 60 mm × 6.4 mm (width × depth) and a height of 1000 mm. The solid phase has 7047



Article

Figure 4. Time evolution snapshots of solid mixing in the bubbling fluidized bed, Ug = 1.2 m/s.

given with 298 K. In the calculation, the time steps for gas and solid motions are 1 × 10−5 and 1 × 10−7s, respectively. Each simulation case is performed with a physical time of 15 s.

4. RESULTS AND DISCUSSION 4.1. Model Validation. The model validation is carried out based on the experimental work of Müller et al.28 The investigated bubbling fluidized bed has a geometrical configuration of 0.044 m × 0.01 m × 0.2 m along the width, depth, and height directions, respectively. As demonstrated in Figure 2, a static bed with a height of 0.03 m is formed before the calculation. The system consists of 9240 particles with diameter of and density of 1.2 mm and 1000 kg/m3, respectively. The physical and numerical parameters adopted in the simulation are listed in Table 2. In order to fluidize the particles in the bubbling fluidized bed, the gas vertically introduced into the system must exceed a critical value, namely, the minimum fluidizing velocity. To obtain this critical value, the gas velocity increases initially from 0 to 0.4 m/s and then decreases to zero. The cross-section averaged value of the pressure drop between the inlet and outlet of the system is recorded. As can be observed from Figure 3a, a sudden change of the profile occurs when increasing the superficial velocity to 0.3 m/s, which corresponds to the breakup of the packed bed. Thus, this predicted value is

Figure 5. Time evolutionary profiles of the mixing index in the systems under different superficial velocities.

kinematic energy of the system decays to zero. In order to identify the mixing process of the binary mixture, different labels are assigned to represent the particle groups before the calculation. Particles are divided into two groups with different temperatures and labels. The temperature of a moving particle is changing due to the heat exchange with its surroundings. However, the label is constant in the tracking procedure, which can be used to identify the spatial mixing level of the solid phase. The particles locating in the upper half of the bed are assigned with 400 K, while the particles in the lower group are 7048



Article

Figure 6. Time evolutionary snapshots of the gas−solid temperature distribution in the slice Y = 0.0032 m of the system with Ug = 1.2 m/s.

identified as the minimum fluidizing velocity, which is the same as the value obtained in the experiment.28 Moreover, the timeaveraged property of gas−solid phases in the bubbling fluidized bed operating in the superficial velocity of 0.9 m/s (Ug = 0.3Umf) are obtained and compared with the experimental data. As shown in Figure 3b−d, the gas phase mainly goes from the central region of the bed due to the nonslip wall effect. The upward motion of the solid phase mainly concentrates in the central region of the bed while the vigorous downward flow appears in the near wall region. Moreover, the predicted quantities of voidage and vertical solid velocity agree well with the experimental data both qualitatively and quantitatively, indicating the internal dynamical property of the gas−solid phase can be captured by the proposed model. In addition, in order to validate the correctness of the coupling model on predicting the mixing behavior, the mixing of the solid phase in the spout-fluid bed experimentally investigated by Jin et al.1 has been studied and the comparison details can be found in our previous work.8 The mixing index of the system reaching dynamical equilibrium agrees well with the experimental data, indicating the key feature of solid mixing behavior can be obtained with the coupling approach. Based on the proposed model, the mixing behavior of the solid phase in a baffle-type internally circulating fluidized bed has been explored.8 Thus, it is confident that the proposed coupling model can capture the important properties of gas−solid hydrodynamics in the bubbling fluidized bed. 4.2. Influence of Bubbles on Solid Mixing. The turbulent and chaotic motion of the solid phase in the bubbling fluidized bed is mainly induced by the rising bubbles. Figure 4

shows the mixing process of two particle groups with different labels under Ug = 1.2 m/s. Before the calculation, the particles in the bed are vertically divided into two groups with the upper one labeled with red color and the lower one labeled with blue, respectively. The effect of bubble motion on the solid mixing can be captured from the rising process of an isolated bubble formed at the initial stage of calculation, as shown in Figure 4b. When the bubble is generated at the bed bottom, the particles with blue color are entrained into its wake and transported upward by the rising bubble, which leads to the vertical mixing between the particles in the lower part of the bed and those in the upper region. In the rising procedure, some of the blue particles drifted outside the bubble wake while others migrate continuously with a significant vertical distance, resulting in the redistribution of solid phase in the bubble wake and a tail with blue particles following the trajectory of the rising bubble. The appearance of the blue tail is due to the solid dispersion and exchange between the blue particles and the red particles, which originally results from the entrainment and drift effects of bubbles. The local exchange of particles from these two groups indicates the strong effect of rising bubbles on the local mixing of the solid phase. After the rising of the bubble, the interface of these two particle groups is distorted with the return flow of particles to compensate for the void produced by the particles leaving with the bubble, giving rise to the vertical mixing of the solid phase. In the following procedure, particles in the nose of the rising bubble are injected into the free domain due to the bubble eruption, leading to the intensive mixing of particles in the horizontal direction. With this process repeated, the 7049



Article

Figure 7. Temperature histogram of the solid phase at several time instants, Ug = 1.2 m/s: (a)t = 0.2 s; (b) t = 0.6 s; (c) t = 1.0 s; (d) t = 2.4 s; (e) t = 4.8 s; (f) t = 6.4 s.

particles of the two groups mix with each other and the boundary between them diminishes gradually. A mixing equilibrium can be established with an operation time of nearly 2.4 s, followed by a stable dynamic mixing state, as shown in Figure 4g,h. Thus, it can be concluded that three dominant mixing mechanisms caused by the rising bubbles can be captured from the snapshots of the particle distribution, namely, the gross mixing process due to the solid transportation with the bubble wake, the local mixing between particles moving outside of the bubbles and the local emulsion phase due to the drift effect of the rising bubbles, and the return flow of particles generated to compensate the void due to the leave of the bubbles. This is in accordance with the conclusions discussed in the literature.29 4.3. Quantitative Description of the Mixing Process. Quantitative investigation of solid mixing behavior is useful for controlling the mixing process in industry. In order to describe the mixing degree, the Lacey mixing index30 is applied to represent the mixing level of the solid phase. In order to calculate the Lacey index of the system, the bed is divided into many small sampling cells whose sizes are larger than the

computational grid but smaller than the bubble size in the bed, then the concentration of a specific particle group in each individual sampling cell is obtained on the basis of the calculation results.31 In the current work, the fluidized bed is divided into 150 cubic sampling cells, with the lengths in the x, y, and z directions divided into 6 × 1 × 25 elements, respectively. Figure 5 shows the time evolutionary profiles of mixing degree under different superficial velocities. In general, the profile slope of the mixing index is found to be the largest at the initial stage, reflecting the mixing of the solid phase in this stage is the most intensive. This largest mixing intensity of the solid phase is mainly due to the complete segregation of the two particle groups. With the continuous mixing of the solid phase, the slope decreases gradually due to the migration of particles from one group to the other one. Finally, the mixing profile fluctuates around a constant value, which corresponds to the micromixing stage of the solid phase. A similar ultimate value of mixing index is reached for the systems of all the superficial velocities. However, these cases reach a desired mixing degree with different times. In general, the enlargement of the inlet velocity reduces the time needed to 7050



Article

over, the conductive heat transfer occurs through the contact area of colliding particles. This can be clearly observed from the temperature distribution of particles locating in the boundary of two particle groups at the initial stage of calculation. With the advancement of the mixing procedure, temperatures of these two particle groups change gradually toward each other, and the temperature variance of the whole system becomes small. Even so, the particles in the system are cooled by the introduced gas of low temperature, which can be inferred from the heated gas flow demonstrated in Figure 6h. As the mixing time is adopted to quantify the time required to achieve a proper spatial mixing status, similarly the time required for the temperature of the system to be equilibrated can be used as an important parameter describing the solid mixing behavior with temperature effect. As vividly shown in Figure 6, the temperature of the system becomes relatively uniform at time t = 6 s, which is larger than the mixing time obtained for the spatial mixing. This raises a caution in the design of the industrial mixing apparatus that the attention should not only be paid on the spatial mixing time but also the time needed to reach a relatively uniform status of the system. Furthermore, it can be captured from Figure 6c,d that the cluster consisting of high temperature particles is formed in some places of the bed, reducing the contact efficiency between two groups with temperature difference. The formation of the cluster can be dispersed with the proper design of the system and the adjustment of important operating parameters, such as the superficial velocity and the operating pressure. Quantitative description focusing on the dynamic evolution of the system temperature is shown in Figure 7. At the initial stage of the calculation, two different parts appear with the temperature of one group mainly concentrating in the nearby of 300 K and the other in 400 K. Temperatures of different groups show different shapes in the procedure toward each other. At time t = 1.0 s, the temperature of the lower part distributes relatively uniformly as compared with that of high temperature. With the heat exchange between the particles, the temperature profiles of these two parts interact with each other. From the time instant of t = 4.8 s to t = 6.4 s, the temperatures of the two groups are more and more close to a central value. As expected, dynamic balance of particle temperature is established after a long time, but it can be inferred that the mean temperature of the whole system will be in a decreasing tendency due to the heat removed from the particles. 4.5. Effect of Superficial Velocity on Temperature Evolution. Superficial velocity has a great influence on not only the spatial mixing but also the heat transfer behavior. The larger the superficial velocity is, the less time is needed for the solid mixing. Moreover, enlarging the superficial velocity increases the convective heat transfer with the fluid medium but decreases the conductive heat flux due to the larger gap between the particles. Figure 8 illustrates the time evolutionary profiles of the average temperature and its standard variance in the systems under different superficial velocities. It can be observed that the standard variance of system temperature decreases more sharply at the initial mixing stage under larger inlet velocity, leading to a smaller time required for the standard variance of the system reaching zero. Small variation of system temperature in the mixing procedure indicates that the temperature of each particle group distributes relatively uniformly at each time instant. Furthermore, it can be found in Figure 8b that the average temperatures of the investigated systems decrease gradually due to the heat removed from the system. Increasing the

Figure 8. Time evolutions of the variance and average of system temperature under different superficial velocities: (a) average temperature; (b) temperature variance.

reach the dynamic mixing equilibrium. The time required for the system reaching the dynamic equilibrium is small at the superficial velocity of 1.2 m/s and perhaps is three times of that under Ug = 0.6 m/s. The mixing index of system under Ug = 0.6 m/s is the lowest in the time interval of t = 0−4 s. This is mainly due to the fact that a large cluster consisting of particles from the same group is formed under a low inlet velocity since few bubbles are formed and the solid motion is not vigorous. Consequently, more time is needed for the particles in the cluster to disperse into the other one. In addition, the mixing index fluctuates around a constant value even after the dynamic equilibrium is reached, indicating that there are many places composed of particles from the same group. Thus, the mixing index is only a systematic parameter adopted to represent the mixing level in a macroscopic view, and the local mixing characteristics of the solid phase cannot be described from the time evolutionary profile. 4.4. Heat Transfer of the Solid Phase. Figure 6 illustrates the snapshots of the temperature distribution of particles in the system with Ug = 1.2 m/s. The internal circulation of the solid phase is triggered by the rising bubbles, resulting in the vigorous mixing of particles. Spatial mixing of the solid phase results in the heat transfer occuring between particles and its surrounding environment. Due to the temperature difference between the gas phase and particles of high temperature, the convective heat flux is transferred from high-temperature particles to the gas phase, which can be captured from the contour plot of the gas temperature in the regions with a cluster of hightemperature particles. This is extremely obvious in Figure 6a−c. After passing the particles of high temperature, the gas temperature is removed due to the heat transferred to the particles of low temperature. This can be observed in Figure 6c. More7051



Article

Figure 9. Time evolutionary profiles of the system temperature histogram, Ug = 1.2 m/s: (a) t = 0.2 s; (b) t = 0.6 s; (c) t = 1.0 s; (d) t = 2.4 s; (e) t = 4.8 s; (f) t = 6.4 s.

taking the position into account will be the same since the label used to identify the different particle group has not been changed. The time evolution of the system temperature is demonstrated in Figure 9. The effect of the initial configuration of different particle groups can be captured from the comparison between Figure 7 and Figure 9. Particles with low temperature are heated easily, and particles with high temperature are cooled correspondingly. The reason is that the cold gas is heated after being introduced into the bed, and convective heat transfer occurs because of the temperature difference between the continuous phase and the high temperature particles, which results in the heat removed from the particles and the rise of gas temperature. With the upward motion of the gas phase, particles with low temperature are encountered, and heat transfer appears due to the temperature difference. This is not the case in the initial configuration with the high-temperature group in the upper region since the heat flux is lost after being removed from the particle group in the upper part of the bed.

superficial velocity lowers the average temperature of solid phase due to more heat flux transferred from particles to the gas phase. Hence, it is advantageous for the system to operate with a high superficial velocity because of the uniform distribution of particle temperature in the mixing procedure, especially for those conditions in which the particle temperature in a specific particle group should behave in a uniform pattern. 4.6. Effect of Initial Configuration on Heat Transfer Behavior. Different configurations of particle groups result in different mixing behaviors and outcome percentages of desired production. The initial configurations of particle groups with different temperatures are important for the mixing procedure when the temperature effect is considered. In order to check the influence of initial temperature configurations of particle groups, the temperatures of two groups are exchanged with each other, namely, the temperature of the upper group is changed to be 298 K and the temperature of the lower one is 400 K. It can be expected that the spatial mixing behavior only 7052



Article

due to the drift effect of bubble phase, and the return flow generated to fulfill the gap left by the rising bubble. (2) The mixing behavior is the most vigorous at the initial stage because of the complete segregation of the system. Then, the mixing intensity decays until a dynamic equilibrium is reached. With increasing superficial velocity, the mixing time needed to reach the dynamic equilibrium decreases, and the procedure is more excellent due to the little formation of particle cluster with particles from the same group. (3) The temperature variance of the system decreases due to the proceeding of the mixing behavior, and average temperature of the system decreases due to the heat transfer to the fluid phase. The time required for the system reaching the temperature dynamic balance is different from the mixing time obtained from the spatial mixing. This is useful in the design of apparatus in industry. Increasing the superficial velocity results in a relatively uniform temperature distribution with the time evolution. (4) The initial temperature configurations of different particle groups have little influence on the mixing result except for the initial mixing stage. This initial effect disappears after enough time. The average temperature of the system is lower if the high-temperature particle group lies in the lower part of the bed. This is helpful if the initial stage is important for the outcome of the product, such as in the quick chemical reaction process.

■

Figure 10. Time evolutionary profiles of variance and average temperature under different initial temperature configurations of particle groups, Ug = 1.2 m/s.

AUTHOR INFORMATION

Corresponding Author

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

The authors declare no competing financial interest.

■

Qualifying this effect of temperature configuration on the mixing behavior can be captured from the time evolutionary profiles of the average temperature and its standard variance. As shown in Figure 10a, time evolutionary profiles of the standard variance of system temperature in the beds with two different temperature configurations nearly coincidence with each other and only minor differences appear at the initial calculation stage. Figure 10b compares the time evolutionary profile of the average temperature of the system. The difference is the largest at the initial stage and then decreases with the time evolution. These properties of the profiles tell the truth that the effect of different initial conditions on the temperature evolution of the system mainly lies in the initial mixing stage. With the process of the mixing behavior, this initial effect is weakened and even disappears after an enough time. This is very useful for the design and optimization of the mixing process which has a strict requirement on both the mixing time and the operating temperature, such as quick chemical reaction encountered in many industries.

ACKNOWLEDGMENTS Financial support from the Major Program of the National Natural Science Foundations of China (Grant Nos. 51390491 and 51390493) and the National Program for Special Support of Top-Notch Young Professionals are sincerely acknowledged.

■

5. CONCLUSION The CFD-DEM coupling approach combined with heat transfer model is used to study the particle mixing behavior in a 3-D bubbling fluidized bed, in which the gas motion is solved in the framework of large eddy simulation while the solid motion is modeled with the discrete element method. Based on the simulation results, several conclusions can be drawn: (1) Three mechanisms can be captured for the effect of bubbles on solid mixing, namely, the gross mixing induced by the particle transported with the bubble wake, the local mixing between moving particles in the wake and the emulsion phase 7053

NOMENCLATURE a, b, c = constant for the calculation of convective heat flux A = surface area of particle, m2 Acontact,i−j = contact area between collision particles i and j, m2 Ck = Kolmogorov constant Cp,f = fluid specific heat capacity, J/(kg·K) Cp = particle heat capacity, J/(kg·K) Cs = Smagorinsky constant coefficient dp = particle diameter, m Fc,ij = contact force between particle i and particle j, N Fd,i = drag force exerted on particle i, N Fp,i = far field pressure force exerted on particle i, N g = gravitational acceleration, m/s2 hi,f = heat transfer coefficient between the particle and fluid, W/(m2·k) I = particle moment of inertia, kg·m2 k = total number of particles in contact with the current one kf = fluid effective thermal diffusivity, W/(m·k) kn,ij = normal stiffness of the solid phase, N/m kp = particle thermal conductivity, W/(m·k) kt,ij = tangential stiffness of the solid phase, N/m M = air molar mass (0.29), kg/mol M = torque on particle, N·m dx.doi.org/10.1021/ie4038399 | Ind. Eng. Chem. Res. 2014, 53, 7043−7055


Article

n = normal component of a variable p = particle phase t = tangential component of a variable

m = particle mass, kg Nu = Nusselt number of particle n = normal unit vector between colliding particles n = number of particles locating in the current cell and the model parameter pf = pressure, Pa Pr = Prandtl number Prf,t = turbulent Prandtal number Q = volume heat flux rate, J/(s·m3) Qi,f = convective heat flux between particle i and fluid, W Qi,j = conductive heat flux between particle i and j, W Qi,wall = conductive heat flux between particle i and the wall, W R = gas constant (8.314), J/(mol·K) Rep = particle Reynolds number Sp = momentum sinking item, kg/(m·s2) Sij = deformation tensor of the filtered field or resolved strain rate, 1/s t = time instant, s tc = particle−particle or particle-wall collision time, s Tf = temperature of fluid phase, K t = tangential vector between two colliding particles Tp = temperature of particle, K Δt = time step, s t = tangential unit vector between colliding particles tc = time of particle−particle collision, s Ug = superficial gas velocity, m/s uf = fluid velocity in the computational cell, m/s Umf = minimum fluidization velocity, m/s uf,i, uf,j = components of gas velocity, m/s vp = particle velocity, m/s vn,ij = normal component of relative velocity between colliding pair, m/s vt,ij = tangential component of relative velocity between a colliding pair, m/s Vp = particle volume, m3 ΔV = volume of the current cell, m3 x, y, z = coordinate axes

Abbreviations

■

CFD = computational fluid dynamics DEM = discrete element method TFM = two-fluid model 2-D = two-dimensional 3-D = three-dimensional

REFERENCES

(1) Jin, B.; Zhang, Y.; Zhong, W.; Xiao, R. Experimental Study of the Effect of Particle Density on Mixing Behavior in a Spout-Fluid Bed. Ind. Eng. Chem. Res. 2009, 48, 10055−10064. (2) Hadi, B.; van Ommen, J. R.; Coppens, M. Enhanced Particle Mixing in Pulsed Fluidized Beds and the Effect of Internals. Ind. Eng. Chem. Res. 2011, 51, 1713−1720. (3) Abanades, J. C.; Grasa, G. S. Modeling the Axial and Lateral Mixing of Solids in Fluidized Beds. Ind. Eng. Chem. Res. 2001, 40, 5656−5665. (4) Cooper, S.; Coronella, C. J. CFD simulations of particle mixing in a binary fluidized bed. Powder Technol. 2005, 151, 27−36. (5) Chang, J.; Wang, G.; Gao, J.; Zhang, K.; Chen, H.; Yang, Y. CFD modeling of particle−particle heat transfer in dense gas-solid fluidized beds of binary mixture. Powder Technol. 2012, 217, 50−60. (6) Liu, Y.; Lan, X.; Xu, C.; Wang, G.; Gao, J. CFD simulation of gas and solids mixing in FCC strippers. AIChE J. 2012, 58, 1119−1132. (7) Yang, S.; Luo, K.; Fang, M.; Fan, J. Discrete element simulation of the hydrodynamics in a 3D spouted bed: Influence of tube configuration. Powder Technol. 2013, 243, 85−95. (8) Fang, M.; Luo, K.; Yang, S.; Zhang, K.; Fan, J. Computational Fluid Dynamics-Discrete Element Method Investigation of Solid Mixing Characteristics in an Internally Circulating Fluidized Bed. Ind. Eng. Chem. Res. 2013, 52, 7556−7568. (9) Luo, K.; Yang, S.; Zhang, K.; Fang, M.; Fan, J. Particle Dispersion and Circulation Patterns in a 3D Spouted Bed with or without Draft Tube. Ind. Eng. Chem. Res. 2013, 52, 9620−9631. (10) He, Y.; Deen, N. G.; Annaland, M. V. S.; Kuipers, J. A. M. Gas− Solid Turbulent Flow in a Circulating Fluidized Bed Riser: Numerical Study of Binary Particle Systems. Ind. Eng. Chem. Res. 2009, 48, 8098− 8108. (11) Feng, Y. Q.; Xu, B. H.; Zhang, S. J.; Yu, A. B.; Zulli, P. Discrete particle simulation of gas fluidization of particle mixtures. AIChE J. 2004, 50, 1713−1728. (12) Deen, N. G.; Willem, G.; Sander, G.; Kuipers, J. A. M. Numerical Analysis of Solids Mixing in Pressurized Fluidized Beds. Ind. Eng. Chem. Res. 2010, 49, 5246−5253. (13) Gui, N.; Fan, J. Numerical study of particle mixing in bubbling fluidized beds based on fractal and entropy analysis. Chem. Eng. Sci. 2011, 66, 2788−2797. (14) Luo, K.; Yang, S.; Fang, M.; Fan, J.; Cen, K. LES−DEM investigation of the solid transportation mechanism in a 3-D bubbling fluidized bed. Part I: hydrodynamics. Powder Technol. 2013, DOI: http://dx.doi.org/10.1016/j.powtec.2013.11.039. (15) Yang, S.; Luo, K.; Fang, M.; Fan, J. LES−DEM investigation of the solid transportation mechanism in a 3-D bubbling fluidized bed. Part II: Solid dispersion and circulation properties. Powder Technol. 2013, DOI: http://dx.doi.org/10.1016/j.powtec.2013.12.049. (16) Smagorinsky, J. General circulation experiments with the primitive equations. I. The basic experiment. Mon. Weather Rev. 1963, 91, 99−164. (17) Lilly, D. K. A proposed modification of the Germano subgridscale closure method. Phys. Fluids 1992, 4, 633−635. (18) Breuer, M.; Rodi, W. Large eddy simulation of complex turbulent flows of practical interest. In Flow Simulation with HighPerformance Computers II, Notes on Numerical Fluid Mechanics, 2nd ed.; Hirschel, E. H., Ed.; Vieweg Verlag: Braunschweig, Germany, 1996; Vol. 52, pp 258−274.

Greek Symbols

β = interphase momentum transfer coefficient, kg/(m3·s) γn = damping coefficient in normal direction, kg/s γt = damping coefficient in tangential direction, kg/s δij = Kronecker delta δn,ij = normal displacements between particle i and particle j, m δt,ij = tangential displacements between particle i and particle j, m Δ = subgrid characteristic length scale, m Δx, Δy, Δz = mesh size in the x, y, and z dimensions, m εf = voidage εp = solid concentration μf = gas dynamic viscosity, kg/(m·s) νt = eddy viscosity coefficient, kg/(m·s) ρf = gas density, kg/m3 ρp = particle density, kg/m3 τf,ij = viscous stress tensor, Pa τf,ijSGS = subgrid stress tensor, Pa ω = particle angular velocity, 1/s Subscripts

c = contact force d = drag force f = fluid phase i = particle i j = particle j 7054



Article

(19) Koch, D. L.; Hill, R. J. Inertial effects in suspension and porousmedia flows. Annu. Rev. Fluid Mech. 2001, 33, 619−647. (20) Cundall, P. A.; Strack, O. D. L. A discrete numerical model for granular assemblies. Geotechnique 1979, 9, 47−65. (21) Hoomans, B. P. B.; Kuipers, J. A. M.; Briels, W. J.; van Swaaij, W. P. M. Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidised bed: A hard-sphere approach. Chem. Eng. Sci. 1996, 51, 99−118. (22) Luo, K.; Fang, M.; Yang, S.; Zhang, K.; Fan, J. LES-DEM investigation of an internally circulating fluidized bed: Effects of gas and solid properties. Chem. Eng. J. 2013, 228, 583−595. (23) Batchelor, G. K.; O’Brien, R. W. Thermal or Electrical Conduction Through a Granular Material. Proc. R. Soc. London, Ser. A 1977, 355, 313−333. (24) Cheng, G. J.; Yu, A. B.; Zulli, P. Evaluation of effective thermal conductivity from the structure of a packed bed. Chem. Eng. Sci. 1999, 54, 4199−4209. (25) Sun, J.; Chen, M. M. A theoretical analysis of heat transfer due to particle impact. Int. J. Heat Mass Transfer 1988, 31, 969−975. (26) Zhou, J. H.; Yu, A. B.; Horio, M. Finite element modeling of the transient heat conduction between colliding particles. Chem. Eng. J. 2008, 139, 510−516. (27) Li, J.; Mason, D. J. A computational investigation of transient heat transfer in pneumatic transport of granular particles. Powder Technol. 2000, 112, 273−282. (28) Müller, C. R.; Holland, D. J.; Sederman, A. J.; Scott, S. A.; Dennis, J. S.; Gladden, L. F. Granular Temperature: Comparison of Magnetic Resonance Measurements with Discrete Element Model Simulations. Powder Technol. 2008, 184, 241−253. (29) Eames, I.; Gilbertson, M. A. Mixing and drift in gas-fluidised beds. Powder Technol. 2005, 154, 185−193. (30) Lacey, P. M. C. Developments in the theory of particle mixing. J. Appl. Chem. 1954, 4, 257−268. (31) Ren, B.; Shao, Y.; Zhong, W.; Jin, B.; Yuan, Z.; Lu, Y. Investigation of mixing behaviors in a spouted bed with different density particles using discrete element method. Powder Technol. 2012, 222, 85−94.

7055


Discrete Element Study of Solid Mixing Behavior with Temperature

Recommend Documents