Article pubs.acs.org/IECR
Particle-Cloud Drag Force in Dilute Particle Systems: Discrete Element Method versus Eulerian Simulations Payman Jalali,* Markku Nikku, and Timo Hyppan̈ en Department of Energy, Lappeenranta University of Technology, 53851, Lappeenranta, Finland ABSTRACT: The discrete element method (DEM) is utilized to calculate the drag force on big spherical objects exerted by the cloud of small particles as they are subjected to an initial relative velocity. The drag force is calculated as the average over single collisions within a certain time interval. Then it is comparable to the calculated values obtained from Eulerian simulations, which are performed on the same particulate system by keeping the big spherical objects in discrete form and turning the cloud of small particles into a continuum phase. The DEM simulations indicate that the average drag force on a big particle in a cluster is lower than that on a single big particle. The decrease in drag force is even larger when the big particles in the cluster are arranged in an ordered configuration such as a simple cubic array. In Eulerian simulations, it is found that the model for kinetic viscosity plays an influential role in the calculated drag force. The two models known as Syamlal−O’Brian and Gidaspow models for kinetic viscosity are employed. Considering the DEM simulation results as the reference cases, it is shown that the Gidaspow model is suitable for the clusters of big particles, while the Syamlal−O’Brian model can reproduce the single big particle motion in a cloud precisely. It is discussed how the distribution of shear stress over the surface of particles in Eulerian simulations may explain the observed differences.
1. INTRODUCTION Current computer technology does not allow us to use Lagrangian methods for simulations of large scale fluidized beds. This is mainly due to the computational complexities imposed by the presence of a huge number of particles (several million to several hundreds of million) in real fluidized beds. Thus Eulerian methods1 have provided more practical ways to simulate the behavior of fluidized beds within an acceptable range of accuracy. However, an Eulerian model needs closure models whose accuracy may vary in different conditions. The closure models for a particulate phase which are alternatively called constitutive laws can be based on either empirical models or kinetic theory of granular flow. A constitutive model is needed to describe the drag force imposed on any particulate phase by other particulate phases. This drag force is a result of collisions and contacts between particles of different shapes, sizes, and properties. A model was presented by Syamlal2 for calculating particle−particle drag force which was derived based on the assumption of binary collisions. This formulation has been also implemented in the Eulerian models of computational fluid dynamics (CFD) packages such as FLUENT. In a recent study,3 the Eulerian simulations in FLUENT based on Syamlal’s particle−particle drag term solely showed that the results differ from the calculations by the discrete element method (DEM) in monodisperse systems. This difference grows as the packing density rises.3 Also, Syamlal’s formula yields the closest result to DEM if the spatial distribution of particles remains random and homogeneous. However, it does not remain so either in a typical DEM simulation or in an experiment because particles tend to form clusters of various sizes,4 which makes the mixture become more inhomogeneous. Eulerian simulations are extremely important for simulations of biomass material processed in fluidized beds.5 However, one of the important characteristics of biomass materials is the large diversity of the size and shape of particles, which are not taken © 2013 American Chemical Society
into account in ref 5. The size of particles in a typical biomass material may vary from a few micrometers to several millimeters depending on the type of biomass. Most of the models based on Eulerian or DEM approaches including the current study consider the particles as perfect spheres. This is rarely the case with biomass particles of irregular shapes, which will be the subject of study for upcoming papers. Such a huge variation of particle size and shape may not only affect the drag between particles and gas phase, but it can also alter the drag forces between solid particles of different sizes in the system. Thus one may apply the Eulerian models cautiously in the simulations of biomass materials. Recently biomass combustion and gasification have been widely studied in literature, and one application of interest has been fluidized beds. In order to correctly predict the fuel flow inside a fluidized bed furnace or gasifier, knowledge of the phenomena affecting the fuel particle is needed. While fluid−solid momentum exchange has been widely studied, the literature offers only a few methods to estimate and model solid−solid momentum exchange. In this paper, the importance of particle−particle drag force in spherical biomass material is studied by means of DEM, which is in parallel simulated by Eulerian simulations in FLUENT. The Eulerian method implemented in the present study offers a way directly applicable to particles of irregular shapes, which will be a future attempt to study the effects of particle shapes. Section 2 is related to a short description of the computational approach, which includes both the DEM and the Eulerian methods. In addition, the description of simulated Received: Revised: Accepted: Published: 4342
October 4, 2012 December 17, 2012 February 21, 2013 February 21, 2013 dx.doi.org/10.1021/ie302704j | Ind. Eng. Chem. Res. 2013, 52, 4342−4350
Industrial & Engineering Chemistry Research
Article
Note that δt is the displacement vector in the tangential direction while δn is the magnitude of the displacement in normal direction specified by the unit vector n. The material stiffness is characterized by Kn and kt, where the damping factors are ηn and ηt in normal and tangential directions, respectively. The relative velocity vr and the slip velocity vs for the pair i and j with translational velocities vi and vj and rotational velocities of ωi and ωj are defined as vr = vi − vj (3)
system is in the last part of this section. In section 3, the results are presented with discussions to show first how the collisional forces are averaged to be comparable with the results of Eulerian models. It is followed by the presentation of results of the drag force between the big particle and the cloud of small particles obtained from the DEM, and the Eulerian methods for single big particle and the clusters made of eight big particles. Finally, some conclusions are made in section 4 to explain the differences observed between different models. 2. Computational Approach. Two computational approaches are utilized in this study as described in the following sections. 2.1. Discrete Element Method (DEM). The main computational approach in this study is based on the DEM, which is a Lagrangian method calculating particles trajectory through successive time steps. It has been frequently used in literature for simulations of granular materials3,6,7 and gas−solid flows.5 The DEM code was developed for the present study according to the theory described in ref 3. In this code, a list of potential neighbors is created for each particle which is updated after 100Δt, where Δt is the simulation time step explained more in the rest of this section. The neighborhood of any particle is a spherical shell of the thickness taken as twice as the diameter of cloud particles. The inner surface of the shell is the particle surface. The above-mentioned values of updating frequency and the neighborhood size are given based on experience but they are not optimized. A small neighborhood with low frequency of updating may lead to the crash of the code during execution due to missing some new contacts, while a large neighborhood with high frequency of updating will significantly extend the computational time. For the current application of biomass materials, the list of potential neighbors for particles with large diameters has to be created differently compared to ordinary particles because the typical number of neighbors for such particles (of very large diameter) could be notably higher than that in monodisperse system. This requires an additional usage of memory to store the list of the neighbors of the large particles. Any particle can have multiple contacts with several particles simultaneously. For any pair of contacting particles of radii Ri and Rj such as that shown in Figure 1, the
vs = vr − (vr ·n)n + (R iωi + R jωj ) × n
(4)
The transition from the static friction regime to the dynamic friction regime of the contact between particles is the slip condition. In other words, if the tangential and normal forces calculated from eqs 1 and 2 satisfy the condition of |Ft| > μf|Fn|, then the tangential force will be turned to Ft = −μf|Fn|t, where t = vs/|vs| and μf represents the dynamic friction coefficient. The tangential displacement is the tangential force divided by the tangential stiffness, δt = Ft/kt. The two stiffness parameters Kn and kt are calculated from the formulas resulted from the Hertzian contact theory. These parameters along with the damping factors are tabulated in Table 1. The trajectory of any Table 1. Stiffness Parameters and Damping Factors in Normal and Tangential Directions parameter
formula
normal stiffness
Kn =
4 R iR j ⎛1 − σ 2 3⎜ E i + ⎝ i
tangential stiffness
kt =
shear modulus
1 − σj 2 ⎞ Ej
⎟ R i + Rj ⎠
8 R iR j ⎛ 2 − σi ⎜ + ⎝ Gi
2 − σj ⎞ Gj
⎟
⎠
δn1/2 R i + Rj
G = E/2(1 + σ )
normal and tangential damping factors
ηn = ηt =
⎞1/2 − 5 ln(e) ⎛⎜ 1 ⎟ δn1/4 K n⎟ ⎜ π 2 + ln 2 e ⎝ 1/mi + 1/mj ⎠
particle is obtained by updating its position after any tiny time step. This requires the knowledge of velocity and acceleration vectors. Both linear (ai) and angular (αi) acceleration vectors are linked to the net force and torque on the particle i via Newton’s second law, F ai = i + g mi (5) αi =
Fi =
Ft = −k tδt − ηt vs
(2)
∑ (Fijn + Fijt ) (7)
j
deformations in normal and tangential directions are denoted by δn and δt. The corresponding forces in normal and tangential directions are calculated by (1)
(6)
where the force and torque are calculated as the net value of all contacting particles as follows
Figure 1. Schematic representation of the contact between two spherical particles of different radii Ri and Rj.
Fn = ( −K nδn 3/2 − ηn vr ·n)n
Ti Ii
Ti =
∑ (R i nij × Fijt ) j
(8)
Here, the contact force from neighboring particle j on particle i is denoted by the superscript ij, whereas it represents the direction from i to j as it is applied to the normal vector n. The 4343
dx.doi.org/10.1021/ie302704j | Ind. Eng. Chem. Res. 2013, 52, 4342−4350
Industrial & Engineering Chemistry Research
Article
moment of inertia for the particle is Ii, and g stands for the gravitational acceleration vector. After calculation of acceleration vectors from eqs 5 and 6, the linear (vi) and angular (ωi) velocity vectors as well as the position (ri) vectors are updated after a tiny simulation time step Δt by the following linearized equations, vi = v i0 + a iΔt
with the DEM simulations via the values of the solid volume fraction and the velocity of the cloud phase given at the inlet (Figure 2a) as the boundary condition of eqs 12 and 13. The
(9)
ωi = ωi0 + αiΔt
(10)
ri = r i0 + viΔt
(11)
The superscript 0 refers to the current time step in which all quantities are known. It is taken as the initial time for the linearized equations of 9−11. 2.2. Eulerian Perspective of Solid−Solid Drag. The continuum description of solid phase is crucial for the practical applications in modeling of gas−solid systems such as fluidized beds of industrial scales. The Eulerian multiphase approach describes all the phases as continua governed by the following conservation equations of mass and momentum for the phase k, ∂ (ερ)k + ∇·(ερ u)k = 0 ∂t
(12)
∂ (ερ u)k + ∇·(ερ uu)k = ∇·(εσ ̿ )k − εk∇p + Fk + εkρk g ∂t (13)
where the index k is applied to all the variables inside the parentheses. The first term in the right-hand side of eq 13 represents the stress term within the phase k. The second and the fourth terms also correspond to the phase k, which stand for the hydrostatic and the gravitational forces, respectively. The third term Fk is the only term that contains various forces exchanged between phases including the drag, lift, and collisional forces. The governing eqs 12 and 13 are solved using the computer package ANSYS-FLUENT, v. 13.0. A primary phase is defined as the background fluid, which is air. Additional phases are defined as secondary phases. In this study, one secondary phase is defined as the cloud particles and the big particles are present as excluded spherical volumes whose surfaces are defined with the wall boundary condition. The grid used in simulations is unstructured polyhedral network generated in GAMBIT, v. 2.4.6. The surface of spheres is covered with triangular mesh elements (as will be shown in Figure 7). The number of elements on each sphere varies between 300 and 400, which explains the size of fairly uniform meshes used in the domain. Smaller meshes yield nonsignificant changes in the results. The gas pressure is initialized with very low values about 1 Pa to minimize its effect on the secondary phase though it cannot be totally suppressed due to the problems associated with the solution convergence. On the other hand, the gas phase velocity at the inlet is set to zero. These settings guarantee very low pressure drop across the bed and locally, which keeps the corresponding drag force on particles about 2 orders of magnitude lower than the drag force between solid particles. The effect of kinetic viscosity model for the particulate phase is explored in section 3.2. The laminar viscous model has been chosen with the steady state, pressure-based solver. The spatial discretization of momentum and volume fraction is given as the first-order upwind. The pressure−velocity coupling is the phase coupled SIMPLE. The physical conditions in the Eulerian simulations are consistent
Figure 2. Geometries used in simulations. (a) Eulerian simulations with eight big particles positioned randomly in the box that is not changed in time. The volume fraction of these particles is calculated as 0.025 within the zone they are distributed. The content of the box is set to the particulate phase plus the background phase as air. (b) DEM simulations with eight big particles positioned initially in the same locations as those in part a. During simulations, the positions of these particles changes due to collisions, but it was observed that the displacements are small. The cross-section of the box is similar to the cross-section of the box in part a. The number of particles in the particulate phase is 40 000, where only half of them in the x direction are shown here. The diameters of big and small particles are 3.6 and 0.3 mm, respectively.
boundary condition at the inlet is velocity inlet with constant volume fraction. The outlet is subject to the pressure outlet boundary condition. The four sides of the system are subject to the wall boundary condition. The big particles are located far enough from the side walls (about 2 particle diameters) so that they are not affected by the walls. The drag force between the particulate phases is calculated as the force exerted from the secondary phase (cloud of small particles) over the surface of big particles shown in Figure 2a. The rest of the simulation box is filled with the secondary phase and the primary phase. In DEM simulations, as shown in Figure 2b, the secondary phase is represented by the small particles. Only half of 40 000 particles in x direction are shown in order to clarify the big particles inside the box. The cross-section of the box in both systems for DEM and Eulerian simulations has the dimensions of 2 cm × 2 cm. The drag force between the particulate phases is calculated by averaging the contact forces between big and small particles, as described in section 3. The primary (gas) phase is absent in DEM simulations, which is also ignorable by some settings in FLUENT calculations. 4344
dx.doi.org/10.1021/ie302704j | Ind. Eng. Chem. Res. 2013, 52, 4342−4350
Industrial & Engineering Chemistry Research
Article
into account through “friction viscosity” and “frictional pressure” variables and the “angle of internal friction” parameter. The latter parameter is given in degrees whose tangent results the friction coefficient used in DEM simulations (Table 2).
2.3. Particulate Systems Used in Simulations. Along with the geometry shown in Figure 2, another system was also generated that included only a single big particle of diameter D = 3.6 mm inserted in a cloud of small particles of diameter d = 0.3 mm. Thus the drag force on a particle can be studied without being affected by other particles in a cluster such as that shown in Figure 2. The particle cloud was generated for DEM simulations by inserting the small particles randomly in the box to yield the bulk volume fraction of ϕs. Obviously, it is measured in a volume excluding the big particles. Table 2
3. RESULTS AND DISCUSSION The results of current simulations are organized in different parts to show how the individual collisions can lead us to average forces on a single big particle. Then it will be shown what the effect of clustering of particles is. The results of DEM and Eulerian simulations will be compared to conclude what types of granular models are more appropriate in Eulerian simulations. It should be mentioned that the relative velocity between the two solid phases is created by giving the initial velocity of big particles in the negative z-direction for DEM simulations, while for Eulerian simulations the big particles are fixed and the initial velocity of the phase corresponding to the small particles cloud is given in the positive z-direction. The latter with fixed big particles is a feasible assumption in Eulerian simulations as it is also verified from DEM simulations. In other words, the relative positions of big particles remain practically unchanged within the simulated time. 3.1. Collisional Behavior and Averaging on Collisions for a Single Particle. The motion of a big particle in a cloud of small particles results in the attenuation of its initial velocity in time due to the collisions with small particles. This is similar to the effect of viscosity in an ordinary fluid flow. In Figure 3a, the thick solid line displays the variation of the velocity of single big particle of diameter 3.6 mm as it moves in a cloud of particles of diameter 0.3 mm using the DEM simulations. The dash line is the quartic fit obtained in the form of Vz = 6434.5t4 − 1709.0t3 + 184.5t2 − 11.1t + 0.480. Note that the fitted curve is obtained with no constraints especially at the origin which displays a small discrepancy with the value of V0. This is to ensure that the fitted curve follows the trend of original one over the entire time span. As seen in Figure 3a, the trend of velocity attenuation is such that its rate decreases in time. The big particle velocity would reach zero if its motion continued for long enough time as gravity is absent. An equivalent turbulent flow of ordinary fluids has been considered for comparison with the flow of small particles cloud around the big sphere. The big particle moves in a turbulent flow of an ordinary fluid with the density of ϕsρs. The drag force8 can be expressed as
Table 2. Numerical Values of Parameters Used in DEM Simulations Nb Nb + Ns ϕs D (mm) d (mm) A (m2) E (GPa) σ μf ρb (kg/m3) ρs (kg/m3) e Δt (μs)
1, 8 40 000 0.05 3.6 0.3 4 × 10−4 11 0.2 0.3 500 2500 0.7 0.5
summarizes the important parameters used to build and simulate the cases. In DEM simulations, the initial relative velocity (V0) between big particles and the cloud of small particles is given by assigning it as the initial velocity of big particles in the negative direction of z, while the small particles in the cloud are initially still. The value of V0 is varied from 0.1 to 2 m/s. Moreover, in DEM simulations, the sides of the container are introduced as elastic nondissipative walls. They were also removed to create open sides, which practically resulted in similar values for the drag force. Note that the DEM simulations revealed that the big particles do not displace notably during their motions through the cloud of small particles within the time simulated. This led us to create the big particles as fixed spherical objects for the Eulerian simulations with the same spatial configurations as in the DEM simulations. The initial relative velocity (V0) between the cluster of big particles and the cloud of small particles in Eulerian simulations is created by assigning it as the initial value of the velocity of the cloud phase. For Eulerian simulations, the setting of different models is as shown in Table 3. Note that the dependence of results on the model types for the kinetic viscosity was found significant as discussed in next sections. It is worth mentioning that the effect of rotational motion of particles in Eulerian simulations is taken
Fd =
kinetic viscosity bulk viscosity frictional viscosity frictional pressure solids pressure radial distribution at contact
(14)
2
where the drag coefficient cD is 0.47 for fully turbulent flow, Vz is the instantaneous velocity of the big particle, and A is the projected area that is πD2/4 for sphere. Having an initial velocity of V0 vertically, the instantaneous velocity of big particle can be integrated from the Newton’s second law written for the particle as
Table 3. Models used for Eulerian Simulations in FLUENT physical quantity
c Dϕρ AVz 2 s s
model name in FLUENT
Vz =
(a) Gidaspow (b) Syamlal−O’Brian Lun et al. Schaeffer Syamlal et al. Syamlal−O’Brian Syamlal−O’Brian
1 1 V0
−
3c Dρϕ s s 4ρ b D
t
(15)
The solid line with normal thickness in Figure 3a represents eq 15 in which V0 = 0.5 m/s and the rest of parameters are introduced in Table 2. This indicates that the particulate flow slows down the velocity of big particle at a higher rate than the 4345
dx.doi.org/10.1021/ie302704j | Ind. Eng. Chem. Res. 2013, 52, 4342−4350
Industrial & Engineering Chemistry Research
Article
calculated mean force 2 orders of magnitude lower than the peak values of collisional forces shown in Figure 3b. Obviously, a smaller sampling time will impose more fluctuations in the averaged force, as shown in Figure 4. The thin solid line
Figure 4. Variation of mean force (Fb, in z-direction) with time obtained by averaging Fz. The thick dash line and the thin solid line represent the mean forces with sampling times of 5000Δt (2.5 ms) and 2000Δt (1 ms), respectively. The thick solid line is obtained from eq 16. The thin dash line is calculated from eq 14 for the particle motion in an equivalent turbulent flow. (inset) Thin neighborhood of the big particle in two time instants of 0 and 0.02 s referred by the arrows. Other parameters of the simulation are set as described in Figures 2 and 3. Figure 3. (a) Variation of single big particle velocity in z-direction with time from DEM simulations (thick solid line) and its quartic fit (dash line). The thin solid line represents the particle velocity attenuation in time in a turbulent flow if the density of the fluid is ϕsρs. (b) Variation of the instantaneous force (Fz) with time exerted on the particle by the particles of the cloud. (inset) Zoomed region indicating the existence of individual collisions separated by different quiescent times. The simulation parameters are as indicated in Table 2 with Nb = 1.
represents the mean force with the sampling time of 2000Δt, which clearly shows more fluctuations than the thick dash line corresponding to the sampling time of 5000Δt. On the other hand, the force may be calculated using the kinematics of the big particle as stated by the Newton’s second law, ρb πD3 dVz · 6 dt ρ b πD 3 = (25738t 3 − 5127t 2 + 369t − 11.1) (16) 6 In Figure 4, the thick solid line stands for the variation of force calculated by eq 16 that yields very smooth curve because of the smooth kinematics. However, due to the dependence of this curve on the fitted function (a polynomial) for the particle velocity, it may not represent the force variation accurately in some regions of the time axis. Moreover, it will not be always possible to fit such polynomials to the velocity−time function if other big particles also exist and collide with each other. Hence, eq 16 is practically applicable for single particle. In Figure 4, the thin dash line displays the variation of force with time in a turbulent flow with equivalent density, as described in eq 14. The particulate flow exerts a force on the big particle that is about 2.5 times larger than that in the turbulent flow in the beginning. The force by the particulate flow sharply decays until it reaches to the turbulent flow curve at t = 0.025 s. Then the turbulent flow force remains above the particulate flow force. 3.2. Drag Force between Cloud Particles and Large Particles. Figure 5 depicts the variation of the averaged force exerted by the small particles cloud on a single big particle with instantaneous velocity from DEM simulations (hollow squares), Eulerian simulations in FLUENT with the Syamlal− O’Brian model (filled squares with blue solid line)9 and the Gidaspow’s model (green circles with green solid line)10 of Fb =
equivalent turbulent flow because of the larger drag force exerted by the small particles of the cloud. It is worth mentioning that one can find out the value of cD for the particle cloud (thick solid line in Figure 3a) corresponds to 1 if eq 15 is used. The force fluctuations experienced by the single big particle are shown in Figure 3b, which is transferred via consecutive contacts with the small particles in the cloud. The inset displays that these contacts can be mainly considered as isolated binary collisions between the big particle and the particles of the cloud. The time between two successive collisions is called the quiescent time within which no contact with the particle happens. It should be noted that the isolated collisions do not represent the average force sensed by the particle such as the force described by eq 14. Therefore, an averaging of collisional forces over a sampling time δt is required. The sampling time must include sufficient number of binary collisions in order to result in statistically meaningful averages. If the number of binary collisions within δt is too few, the averaged force will be highly fluctuating. So, δt must be reasonably large in order to capture the trend of force as it is predicted by Eulerian methods or any relevant experimental measurements. However, if it is too large, the function of force with time will lose its trend. In the current study, an appropriate sampling time δt is obtained to be as 5000Δt, which corresponds to 2.5 ms. The big particle spends a great deal of the time quiescently resulting the 4346
dx.doi.org/10.1021/ie302704j | Ind. Eng. Chem. Res. 2013, 52, 4342−4350
Industrial & Engineering Chemistry Research
Article
DEM and Eulerian simulation results. The red dashed line representing eq 17 is obtained using εb = 0.001. It represents very dilute state where big particles can be assumed isolated. The green line stands for the Eulerian simulations with the Gidaspow’s model of kinetic viscosity, which shows significant deviation from the above-mentioned results. On the other hand, the Eulerian simulations were performed along with the DEM simulations for the clusters of big particles shown in Figure 2. The comparison of various simulations is demonstrated in Figure 6. Two different arrangements of
Figure 5. Mean force (Fb, in z-direction) versus the velocity of single big particle. Hollow squares represent the DEM results. Filled squares with blue solid line stand for the results of FLUENT calculations using Syamlal−O’Brian model (model b for kinetic viscosity in Table 3), while the green solid line with circles are the results of calculations using Gidaspow’s model (model a in Table 3). Finally, the red dash line represents the predictions by Syamlal’s formula (eq 17) using εb = 0.001 and εb = ϕs = 0.05. Other parameters are as introduced in Table 2.
kinetic viscosity, as well as the formula presented by Syamlal2 seen as red dash line. The formula is written as follows Fbs = −S bs(us − ub)
(17)
Figure 6. Average force (Fb, in z-direction) on a big particle versus the mean velocity of particles in a cluster of eight particles. The thick solid line represents the curve obtained for single big particle (Figure 5). Red squares and blue circles stand for the force variation in the two random clusters drawn in the background with the colors corresponding to the squares and circles. Black diamonds show the third cluster of regularly arranged particles shown as background in black color. The results of averaged force on a big particle using FLUENT for the cluster of eight particles (shown in red) are plotted as the dashed line (model a for kinetic viscosity) and a dotted−dashed line (model b for kinetic viscosity). The inset shows a focused area on the lower limits of velocity (from −0.4 to −0.1 m/s) from the main graph.
Here, Fbs is the volumetric momentum transfer between two solid phases s and b, having the instantaneous velocities of us and ub, respectively. The drag coefficient Sbs is defined as S bs =
3(1 + e)(π /2 + μf π 2/8)εbεsρb ρs (db + ds)2 g0|u s − ub| 2π(ρb db3 + ρs ds 3) (18)
where g0 is the radial distribution function at contact. For a mixture containing two solid phases b (for big particles) and s (for small particles) of spherical particles and a gaseous phase g, it is calculated as g0 =
⎛ εb ε⎞ 3d d 1 + 2 b s ⎜ + s⎟ εg ds ⎠ εg (db + ds) ⎝ db
randomly distributed big particles (with the same volume fraction of about 0.05) are used in DEM simulations. The mean force on big particles of each cluster is plotted versus velocity as blue circles and red squares. The representation of clusters is shown as blue and red spherical particles in the background. The two data sets basically follow the same curve of force− velocity in average. The thick solid line represents the force on a big single particle. This line lies slightly above the majority of data points corresponding to the clusters. So, any fitted curve of big particle in a cluster would be only slightly below the one of single big particle (thick solid line). Moreover, the difference is more evident in the lower range of absolute velocity than the upper range. Figure 7 shows that the big particles in a cluster experience different volume fractions (Figure 7a) and velocities (Figure 7b), as seen from the Eulerian simulations. For some particles, it may result in lower forces and for some others in higher forces. Nonetheless, the force is decreased for an average big particle in a cluster. In order to show the effect of the shadow of big particles on the forces exerted on the big particles behind, a certain cluster is built with simple cubic (SC)
(19)
In addition to g0 and the magnitude of relative velocity of phases |us − ub|, the drag coefficient contains the collisional properties of solid particles (e and μf), the diameter of particles (db and ds), the material density (ρb and ρs), and the volume fraction of phases (εb and εs). According to the results of Figure 5, the Eulerian simulations (filled squares with blue line) with the Syamlal−O’Brian model of kinetic viscosity follow the DEM results (hollow squares) closely. In addition, the curve resulted from eq 17 is plotted in red dashed line for the given range of velocities and the volume fraction of the cloud equal to 0.05. Note that the value obtained from eq 17 is multiplied by the volume of the container to yield the force in N, which is plotted as the red dash line in Figure 5. Since there is only one big particle in the DEM simulations, the volume fraction of the corresponding phase cannot be calculated precisely. Therefore, the suitable value of εb in eqs 18 and 19 is obtained by testing different numerical values to reach the best agreement with the 4347
dx.doi.org/10.1021/ie302704j | Ind. Eng. Chem. Res. 2013, 52, 4342−4350
Industrial & Engineering Chemistry Research
Article
Figure 7. Distribution of (a) volume fraction and (b) velocity magnitude of the cloud within the central XZ plane and two parallel XY planes. Particles are shown in black. Colorbars corresponding to each part are given with the unit of velocity in meters per second in part b.
Figure 8. Distribution of wall shear stress (Pa) from the particle cloud phase on the surface of big particles obtained from FLUENT simulations using the Gidaspow and Syamlal−O’Brian models for kinetic viscosity. (a) Single particle using the Gidaspow model. (b) Single particle using the Syamlal−O’Brian model. (c) Cluster of eight particles using the Gidaspow model. (d) Cluster of eight particles using the Syamlal−O’Brian model. Colorbars are plotted separately for each part. The unit is pascal for all parts.
arrangement of eight particles, as shown in the background of Figure 6 with black particles. The symbols with diamond shape
in Figure 6 represent the calculated values for force versus velocity on the average big particle in the SC cluster. The 4348
dx.doi.org/10.1021/ie302704j | Ind. Eng. Chem. Res. 2013, 52, 4342−4350
Industrial & Engineering Chemistry Research
Article
DEM simulations reveal that the fuel particle in a cluster experiences, on average, a lower force than the single fuel particle case as a consequence of upstream particles shadow. For the single fuel particle, Syamlal−O’Brian’s model for kinetic viscosity in Eulerian simulations yields the force in the best agreement with the DEM simulations. Gidaspow’s model produces higher values of force for the single fuel particle. Surprisingly, the force calculated from the Eulerian method for a fuel particle in a cluster using Gidaspow’s model exhibits a good agreement with the DEM results, whereas Syamlal− O’Brian’s model considerably underestimates the force on the fuel particle in the cluster. This study attempts to clarify the role of solid−solid drag force in the dynamics of fluidized beds especially when there are two solid phases with a size ratio over 10. DEM simulation has been used as a direct way of measuring the force which can verify the available models in Eulerian simulations, as explained above. The methodology used here can be extended to any other studies of this kind and also deeper insights into the differences between continuum models. For instance, it is shown in this paper that the simulations using the two models of kinetic viscosity lead to different distributions of shear stress over the surface of single fuel particle. The magnitude is also found dissimilar. However, the distribution of wall shear stress is found similar for fuel particles in a cluster while its magnitude differs in the same fashion as for single particle. It is an interesting result that the valid model for kinetic viscosity switches from Syamlal−O’Brian’s model in a single fuel particle to Gidaspow’s model in the cluster of fuel particles. The reason behind this observation, which does not change with the configuration of the cluster, may be associated with the development of different flow fields that can regulate the forces exerted on the surface of particles. On the other hand, further Eulerian simulations can be attempted in which the fuel particle phase is also considered as a continuum medium. In addition, precise experimental methods are highly appreciated in future works for further verification of the models presented here.
average force is considerably decreased due to the exposure of back row particles to lower velocities and volume fractions. The results of the Eulerian simulations by FLUENT are also shown in Figure 6 by the dashed and dashed−dotted lines corresponding to Gidaspow and Syamlal−O’Brian models for kinetic viscosity, respectively. Interestingly, these results reveal that the Gidaspow model produces the forces in good agreement with the DEM results in a cluster, whereas the Syamlal−O’Brian model significantly underestimates the values of average forces in a cluster. In order to elucidate this difference, the distribution of wall shear stress from FLUENT simulations is shown in Figure 8 for the single big particle and the cluster of eight big particles using both models of kinetic viscosity. For the single particle, Gidaspow’s model (Figure 8a) yields a sticky cloud as the wall shear stress covers considerably large area on the particle surface. In contrast, the Syamlal− O’Brian model (Figure 8b) shows a separation of the cloud flow from the surface of the big particle earlier than the former model in which the value of shear stress is also lower. The separation point observed by the Syamlal−O’Brian’s model is qualitatively in agreement with the one observed in 2D hard sphere simulations11 of the encounter of an object with a granular stream. The Syamlal−O’Brian’s model regenerates the values of force closely to the DEM results for the single particle. When the single big particle is replaced by the cluster of big particles, the separation point is predicted identical by both models (Figures 8c and d), and the discrepancy is found only in the magnitude of wall shear stress. This difference is evidently resulting from the higher estimation of kinetic viscosity by Gidaspow’s model than Syamlal−O’Brian’s model especially in the lower range of volume fractions below 0.1. Obviously, a deeper insight into the root of the differences between the models remains a task for future studies. The current paper only tries to compare the results of Eulerian models for the kinetic viscosity with the DEM simulations to validate them.
4. CONCLUSION Using the DEM simulations, the solid−solid drag force is calculated which is exerted by the bed of spherical particles of 300 μm on either single spherical particle of 3.6 mm, or a cluster of eight such particles. The big particles represent biomass fuel particles that are injected into a bed of smaller particles such as sand fluidized by a gas such as air. The gas phase is absent in current simulations to exclude its effects on the results. In parallel, Eulerian simulations are conducted alternative to the system simulated by DEM, assuming the phase of smaller particles as a continuous solid phase while big spherical particles were considered as spherical surfaces arranged in the same configuration as in DEM subjected to the wall boundary conditions. As the volume fraction of the particle cloud is low (about 0.05), the interaction of its particles with the fuel particle (big particle) is seen to be mainly via individual collisions, as DEM simulations reveal. For proper comparisons of the force from DEM to that from the Eulerian simulations, the forces calculated from individual collisions are averaged over windows of time with certain width. The average force is comparable to the values resulting from the continuous phase formulation. The force is also compared to an equivalent turbulent flow if the physical properties of the fluid are assumed similar to the solid phase. The comparison shows faster decay of the force for the granular cloud than the turbulent flow. This paper shows that the granular kinetic viscosity is the key parameter in the simulations using the Eulerian method. The
■
AUTHOR INFORMATION
Corresponding Author
*E-mail: pjalali@lut.fi. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The authors would like to acknowledge the financial support provided by the Academy of Finland under Grant Nos. 123938 and 124368. The authors had very valuable discussions with Dr. Kari Myöhänen from LUT that are kindly acknowledged.
■
4349
NOMENCLATURE a = acceleration vector, m/s2 A = particle projected area, m2 cD = drag coefficient d = particle diameter, m D = diameter of fuel particle, m e = coefficient of restitution E = Young’s modulus, Pa Fi = net contact force on particle i, N Fn = normal contact force between particles i and j, N Ft = tangential contact force between particles i and j, N Fk = net interaction force on phase k, N dx.doi.org/10.1021/ie302704j | Ind. Eng. Chem. Res. 2013, 52, 4342−4350
Industrial & Engineering Chemistry Research
Article
(3) Jalali, P.; Hyppänen, T. Verification of continuum models for solids momentum transfer by means of discrete element method. Ind. Eng. Chem. Res. 2010, 49, 5270−5278. (4) Jalali, P.; Li, M.; Ritvanen, J.; Sarkomaa, P. Intermittency of energy in rapid granular shear flows. Chaos 2003, 13, 434−443. (5) Deza, M.; Franka, N. P.; Heindel, T. J.; Battaglia, F. CFD modeling and X-Ray imaging of biomass in a fluidized bed. J. Fluid Eng. 2009, 131, 111303. (6) Tsuji, Y.; Tanaka, T.; Ishida, T. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technol. 1992, 71, 239−250. (7) Cundall, P. A.; Strack, O. D. L. A discrete numerical model for granular assemblies. Geotechnique 1979, 29, 47−65. (8) Schlichting, H. Boundary Layer Theory, 7th ed.; McGraw-Hill: New York, 1979. (9) Syamlal, M.; Rogers, W.; O’Brien, T. J. MFIX Documentation: Vol. 1, Theory Guide; National Technical Information Service, Springfield, VA, 1993; DOE/METC-9411004, NTIS/DE9400087. (10) Gidaspow, D.; Bezburuah, R.; Ding, J. Hydrodynamics of circulating fluidized beds, kinetic theory approach. In Fluidization VII, Proceedings of the 7th Engineering Foundation Conference on Fluidization, Brisbane, Australia, May 3−8, 1992; pp 75−82. (11) Jalali, P. Flow characteristics and stresses on cylindrical objects immersed in a flow of inelastic hard disks. Powder Technol. 2012, 219, 217−227.
Fbs = net momentum transfer rate between phases b and s, N/m3 Fd = force on fuel particle in an equivalent turbulent flow, N Fb = force on fuel particle exerted by the cloud, N g = gravitational acceleration vector, m/s2 g0 = radial distribution function at contact G = shear modulus, Pa I = moment of inertia, kg m2 Kn = normal stiffness, Pa m1/2 kt = tangential stiffness, Pa m m = particle mass, kg n = unit vector along the centerline of particles i and j N = number of particles p = hydrostatic pressure, Pa r = positional vector of a particle, m R = particle radius, m Sbs = solid−solid drag coefficient, kg/(m3 s) t = tangential unit vector at contact point between particles i and j t = time, s Ti = net torque on particle i, N m u = velocity vector of continuous phase k, m/s v = velocity vector of particle, m/s Vz = magnitude of instantaneous velocity of particle, m/s V0 = initial relative velocity between big particles and the cloud, m/s vr = relative velocity vector between particles i and j, m/s vs = slip velocity vector between particles i and j, m/s Greek Letters
α = rotational acceleration, rad/s2 δt = tangential deformation vector, m δn = normal deformation magnitude, m Δt = time step size, s δt = averaging time span, s ε = phase volume fraction ϕs = bed volume fraction ηn = normal damping factor, kg/s ηt = tangential damping factor, kg/s μf = dynamic friction coefficient ρ = density, kg/m3 σ = Poisson’s ratio σ̿ = phase stress tensor, Pa ω = angular velocity, rad/s Subscripts
b = bed g = gas phase i = particle index j = particle index k = phase index s = fuel Superscripts
0 = corresponding to current time step ij = exchanged between particles i, j k = phase index
■
REFERENCES
(1) Enwald, H.; Peirano, E.; Almstedt, A. E. Eulerian two-phase flow theory applied to fluidization. Int. J. Multiph. Flow 1996, 22, 21−66. (2) Syamlal, M. The Particle-Particle Drag Term in a Multiparticle Model of Fluidization; Topical Report, National Technical Information Service: Springfield, VA, 1987; DOE/MC/21353-2373, NTIS/ DE87006500. 4350
dx.doi.org/10.1021/ie302704j | Ind. Eng. Chem. Res. 2013, 52, 4342−4350
