5270
Ind. Eng. Chem. Res. 2010, 49, 5270–5278
Verification of Continuum Models for Solids Momentum Transfer by Means of Discrete Element Method Payman Jalali* and Timo Hyppa¨nen Faculty of Technology, Department of Energy Technology, Lappeenranta UniVersity of Technology, 53851, Lappeenranta, Finland
This paper focuses on the study of momentum transfer rate between granular phases in particulate flows that interact in low and moderate volume fractions. Two systems have been created and simulated. The first system included two granular phases initially mixed homogeneously in a vessel with rectangular cross section. The second system included two granular phases as two concentric cylinders initially separated, but one phase was stationary and the other phase penetrated across a circular opening to the other phase. In one hand, discrete element method (DEM) was employed as the particle trajectory following method, and on the other hand, an Eulerian approach was used as the continuum model to study the momentum transfer between granular phases. The results are also assessed based on the formula derived by Syamlal [Syamlal, M. The Particle-Particle Drag Term in a Multiparticle Model of Fluidization. Topical Report, DOE/MC/213532373, NTIS/DE87006500, National Technical Information Service, Springfield, VA. 1987] for the solid-solid momentum transfer rate. This formula yields solid-solid drag forces larger than those obtained from the DEM simulations especially as volume fraction grows. Our results revealed that in the Eulerian model (in FLUENT) the relative velocity of granular phases during the interacting period may not be influenced only by the solid-solid drag force. Though the two approaches calculate qualitatively similar results in time, there are quantitative differences in calculation of phase volume fractions and the momentum of phases. The interacting phases are highly scattered in the Eulerian model predictions, which may be due to larger drag forces and a lower energy dissipation rate existing in the model. 1. Introduction Numerical simulations of fluidized beds have been increasingly demanded by industrial needs to improve the efficiency of solids mixing, segregation, and transport. Multiscale modeling methods have been introduced by Tsuji1 in simulations of gas-solid flows. Some of the simulation methods combined of the Lagrangian approach for solid phase (microsclae) and the Eulerian approach for gas phase (mesoscale) have been successfully implemented in recent years. For instance, Tsuji et al.2 used nearly 4.5 million particles and over half a million grids for the gas phase in the Lagrangian-Eulerian (LE) simulations of fluidized beds. In the LE simulations of gas-solid flows, the Lagrangian approach is commonly called the discrete element method (DEM), which is used for tracking the path of individual solid particles dictated by contact forces from other solid particles as well as the forces exerted on the particle by the gas phase. The mechanics of contacts between particles is described using the Hertzian contact theory which relates the stiffness of contacting particles to the mechanical properties such as the Young’s modulus and Poisson’s ratio.3,4 In addition to stiffness, the damping factor and the Coulomb friction are taken into account in the formulations of the contact mechanics of solid particles. In fact, the DEM employs very sophisticated theories to describe the instantaneous situation of the contact between solid particles. Though the computational cost of such calculations is incredibly high there are other benefits which cannot be ignored. For instance, the credibility of Eulerian simulations for granular phase under various circumstances is associated with the applicability of kinetic theory formulations used in simulations.5 A main assumption in kinetic theory is the assumption of Maxwell-Boltzmann velocity distribution. * To whom correspondence should be addressed. Tel.: +358 5 621 2791. E-mail:
[email protected].
If such assumption fails for any reason, such as having nonequilibrium state, the consequent calculations of collision frequency and transport properties will not be accurate. However, the DEM simulations do not face such simplifying assumptions, but they can provide necessary data to modify the assumptions of kinetic theory. Since the simulations of largescale fluidized beds are computationally feasible only by Eulerian approach, it is very important to assess the implication of this approach in the computational fluid dynamics (CFD) packages. In this study, the Eulerian approach has been implemented besides the DEM simulations to assess their differences and suggest possible ways to improve the models for better predictions of the momentum transfer rate between granular phases. In this context, the formula derived by Syamlal6 for the solid-solid momentum transfer rate is further utilized to assess the outcome of the simulation methods. In section 2, the methods are described with the governing equations to clarify the type of simulation models used in this study. The systems under study are discussed in details in section 3, which is followed by the results and discussions for each system in section 4. Finally, section 5 discusses the final conclusions of the study. 2. Methods and Governing Equations Discrete Element Method (DEM). In the simulations of granular phases based on the DEM, individual particles are tracked after calculation of acceleration from the net forces exerted on each particle instantaneously. This method is fundamentally a type of molecular dynamics (MD) simulations in which particles are soft spheres interacting only when they contact.7 The forces exchanged between two contacting spheres and the relation between forces and deformation of spheres are
10.1021/ie901538w 2010 American Chemical Society Published on Web 12/16/2009
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010 4
determined based on the Hertzian contact theory, or some modified versions of it.8 In this study, we have developed a computer code for the DEM simulations with the generalization of the model in Tsuji et al.3 to different size and material particles, as described below. Two particles i and j are pressing each other, and the resulting force with corresponding displacements for particle j are shown in Figure 1. The forces in normal and tangential directions can be calculated based on a linear spring-dashpot assembly by FCnij ) (-Knδ3/2 nij - ηnvrij · nij)nij
(1)
FCtij ) -ktδtij - ηtvsij
(2)
Here, the direction of the normal force FCnij is always along the normal unit vector nij. The normal deformation δnij is the overlapping distance between the particles. The tangential displacement vector is denoted by δtij. The effects of spring and dashpot in the model are characterized by the stiffness in normal and tangential directions (Kn and kt) for the spring, and the damping factor in normal and tangential directions (ηn and ηt) for the dashpot. The direction of the tangential (shear) force FCtij is dependent on the relative slip velocity vector vsij. The relative velocity vector vrij, and vsij are defined by vrij ) vi - vj
(3)
vsij ) vrij - (vrij · nij)nij + 0.5(diωi + djωj) × nij
(4)
The angular velocity vector of particle i with diameter di is denoted by ωi. The tangential force is calculated using eq 2 as long as the slip condition is not satisfied. The slip condition is actually the transition from the static friction regime to the dynamic friction regime of the contact between particles. In other words, if the tangential and normal forces calculated from eqs 2 and 3 satisfy the condition of |FCtij| > µf|FCnij|, then the tangential force will be changed to FCtij ) -µf|FCnij| tij, where tij ) vsij/|vsij| and µf represents the dynamic friction coefficient. The tangential displacement is the tangential force divided by the tangential stiffness, δtij ) FCtij/kt. The stiffness is calculated based on the Hertzian contact theory. For two spheres of radii Ri and Rj, the normal stiffness is determined by Kn )
(
4√RiRj
1 - σj2 1 - σi 3 + Ei Ej 2
)√
(5)
Ri + Rj
where σ and E are Poisson’s ratio and Young’s modulus with corresponding indices for each particle to differentiate their material properties.9 This equation can be also used for the contact between a particle and a wall in which the radius of the wall tends to infinity. The tangential stiffness is defined as
kt )
(
8√RiRj 2 - σj 2 - σi + Gi Gj
δnij1/2
)
5271
(6)
√Ri + Rj
where G stands for the shear modulus related to E and σ via G ) E/2(1 + σ). The damping factor has two normal and tangential components that may be assumed identical as j Kn)1/2δnij1/4 ηn ) ηt ) R(m
(7) -1
Note that the reduced mass is m j ) (1/mi + 1/mj) , where mi and mj are the masses of contacting particles. The right form of the coefficient R is related to the coefficient of restitution e as follows:10 R)
- √5 ln(e)
√π2 + ln2 e
(8)
The equation of motion for any individual particle is used to determine the trajectory of particle within consecutive time steps. The contact forces are determined between particles using eqs 1 and 2, which leads to the calculation of net contact force (FCi) and contact torque (TCi) on particle i, as FCi )
∑ (F
Cnij
+ FCtij)
(9)
× FCtij)
(10)
j
TCi )
1 2
∑ (d n
i ij
j
The summation is taken over all particles j in contact to particle i. Using the net forces and torques calculated instantaneously, the translational (ai) and rotational (ri) accelerations are determined by ai )
FCi +g mi
(11)
TCi Ii
(12)
ri )
The vector of gravitational acceleration is g, and the moment of inertia for particle i is denoted by Ii. Eulerian Multiphase Approach. The continuum description of solid phase plays a critical role in simulations of gas-solid 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, ∂ (εF)k + ∇ · (εFu)k ) 0 ∂t
(13)
Figure 1. (a) Tangential and normal contact forces between two particles i,j. Normal and tangential unit vectors are also shown. (b) Normal and tangential deformation of contacting particles.
5272
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
∂ (εFu)k + ∇ · (εFuu)k ) ∇ · (εσ c)k - εk∇p + Fk + εkFkg ∂t (14) where the index k distributes over all the variables inside the parentheses. The first term in the right-hand side of eq 14 represents the stress term within the phase k. Also, the second and the fourth terms correspond to the phase k, which represent the hydrostatic and the gravitational forces acting within the phase, respectively. The third term Fk is the only term that contains various forces exchanged between phases including the drag, lift, and collisional forces. Equations 13 and 14 are solved by means of the computer package FLUENT, v. 6.3.26. In FLUENT, a primary phase should be defined which is assigned as the background phase. In this study, it is defined as air. Additional phases are defined as secondary phase, which in this study are two particulate phases. To suppress the effect of gaseous phase on particulate phases, the mutual drag coefficients are deactivated except between the particulate phases. Moreover, gas pressure is initialized with very low values and there are no boundary conditions that could increase the gas pressure in the system. As a result, the values of maximum gas pressure in simulations remained well below 10-3 Pa through the entire simulation time steps. The model settings are displayed in Table 1. The details about each model can be found in the corresponding references given in the FLUENT user’s manual. The laminar viscous model has been chosen with the transient, pressure-based solver. The second order implicit method has been used for the unsteady formulation. The spatial discretization of momentum and volume fraction is given as the first order upwind. The pressure-velocity coupling is the phase coupled SIMPLE. 3. Studied Systems of Solid Mixtures Homogeneously Mixed Particulate Phases. A system consisting of two identical, homogeneously mixed solid granular phases has been considered in a straight channel with rectangular cross section as shown in Figure 2a. The entire container with dimensions A × B × H is initially filled with two well-mixed particulate phases of volume fractions ε1 and ε2. Velocity of particles in one phase is set to zero and in the other phase is unidirectional along the z-axis with the magnitude of V0. Gravity is not included in simulations because our main purpose is to investigate particle-particle forces exchanged between the particulate phases. Note that the sidewalls of the container are assumed as elastic rigid walls while the upper and lower faces are open boundaries. The presence of sidewalls is essential since they prevent the scattering of solid particles out of the container. However, the analysis of results has been performed within a control volume (CV) located in the core of the container to keep a safe distance from the sidewalls, as shown in Figure 2a. The precise dimensions of the container and the CV are sketched in Figure 2b along with the numerical grids used in the Eulerian simulations. Note that in all Eulerian simulations, the dimensions of the container are set as A ) B ) 0.20 m and H ) 2.0 m. Table 1. The Setting of Eulerian Model in FLUENT modeled quantity
model name
drag coefficient granular viscosity granular bulk viscosity frictional viscosity frictional pressure granular temperature solids pressure radial distribution function
syamlal-obrian-symmetric syamlal-obrian lun-et-al schaeffer syamlal-et-al algebraic syamlal-obrian syamlal-obrian
Figure 2. (a) The setup for the homogeneously mixed solid phases. The large rectangular container with dimensions A × B × H contains particulate phases. The small rectangle (dotted) is an arbitrary open control volume (CV) with dimensions a × b × h for which the average particle-particle drag force is calculated. (b) The outline of the geometry created for Eulerian simulations with the computational grids at a cross section of the container.
The CV is located symmetrically so that a ) b ) 0.10 m, h ) 1.0 m, and h0 ) 0.5 m. In the Eulerian modeling of this system, the interaction term (Fk in eq 14) only consists of the force due to collisions between the particles of granular phases. Syamlal6 presented the following expression for calculating the momentum transfer between two solid phases m and l with instantaneous average velocities of ul and um, as Fml ) -Sml(ul - um)
(15)
Sml ) (3(1 + e)(π/2 + µfπ2 /8)εmεlFmFl(dm + dl)2g0 |ul um |)/(2π(Fmdm3 + Fldl3))(16) g0 )
3dmdl 1 + 2 εg εg (dm + dl)
M
εi
∑d i)1
(17)
i
Here, Sml is the drag coefficient expressed as a function of the material properties of solid phases including restitution coefficient e, friction coefficient µf, solids material densities Fm and Fl, as well as the particles sizes dm and dl, solids volume fractions εm and εl, radial distribution function at contact for the entire mixture g0, and the absolute instantaneous relative velocity of the phases |ul - um|. Note that Fml represents the momentum transfer rate per unit volume. Therefore, its value will be independent of the size of the CV within which Fml is calculated, if all quantities are distributed uniformly inside the CV. In DEM simulations, Fml can be directly calculated within a small averaging time δt by Nml
∑ m (v′ - v ) i
Fml )
i)1
i
i
(18) habδt where h, a, b are the dimensions of the CV in Figure 2a, where there are Nml particles inside the CV used for calculation of Fml within time δt. It is worth mentioning that Nml does not generally
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
5273
Figure 4. Variation of relative velocity between two granular phases (|ul - um|) with time. Various symbols and lines are introduced in the inset. All simulation parameters are presented in Table 2.
Figure 3. (a) The initial setup for the system consisted of a cylindrical granular stream (marked as 1) of uniaxial velocity in the direction of cylinder axis that penetrates into a cylindrically distributed granular medium (marked as 2) having zero initial momentum. (b) Demonstration of a circular element with given dimensions for averaging different quantities. (c) The front view and (d) the isometric view of the structured grid for Eulerian simulations. Region 2 is specified inside the dashed circle.
remain constant in time unless δt is taken small enough so that the variation of Nml almost maintains unchanged. Interpenetration of Two Granular Media. Unlike the above-mentioned studied case, the second case consists of two granular phases which are not initially mixed. Figure 3a depicts the initial configuration of the studied system. The granular stream 1 consists of spherical particles arranged initially in random positions within a cylindrical region of diameter ds and length l. The velocity of particles in this region is downward that is in the negative direction of z-axis. The granular medium 2 consists of spherical particles of the same size and type of the particles in region 1, which are initially located randomly in a cylindrical zone of diameter D and length L with zero velocity. It is worth mentioning that the creation of particles within region 1 is necessary in DEM simulations, while the corresponding phase enters into the granular medium 2 via the velocity inlet boundary condition in Eulerian simulations. The number of particles inside regions 1 and 2 are initially as N1 ) 3ε1ds2l/(2d3) and N2 ) 3ε2D2L/(2d3), respectively. Figure 3b displays a circular element used for averaging purpose in both DEM and Eulerian simulations. The mean value of any quantity in the element can be projected on a vertical plane such as the yz plane. The radial position of the element (r) is not limited to the radius of the granular medium (D/2), but it is also taken in larger radii to study the lateral transfer of momentum due to the penetration of two phases. Therefore, the domain of solution in Eulerian simulations is taken larger than the granular medium 2, as shown in Figure 3c,d. This allows us to see the lateral scattering of granular phases upon penetration at any time. The diameter of the domain of solution is taken three times larger than the initial diameter of the granular medium (D), as indicated in Figure 3c.
The boundary conditions over various surfaces are as follows. Two interior boundaries are defined as the lateral faces of two concentric cylinders with the diameters of ds and D. Except a circular opening with diameter ds at the top of the domain of solution (the cylindrical zone with diameter 3D and length L) the rest of its outer faces are subjected to the pressure-outlet boundary condition. The circular opening is the inlet of the granular phase 2 (in region 1) in Eulerian simulations where the inlet velocity, granular temperature, and phase volume fraction are set to the given values over the face. Obviously, these parameters are set to zero over the inlet for the granular phase 1 (in region 2). The phase interaction between the granular phases is defined via the drag coefficient of the “syamlal-obriensymmetric” scheme in FLUENT and collisions are specified by e. Other model settings are also set as mentioned in Table 1. 4. Results and Discussion Solid-Solid Momentum Transfer Rate. The above-mentioned case consisting of two initially well-mixed granular phases is simulated using the DEM and Eulerian methods to study the solid-solid momentum transfer rate between the phases in time. Though the volume of system in the Eulerian model is kept constant through different simulations, as shown in Figure 2b, the volume of system depends on the phase volume fraction in the DEM model because the number of particles are kept constant in all simulations. This assumption is justified because the momentum is basically transferred by the contacts between particles. Hence the number of particles will be a major factor in determining the rate of momentum transfer between phases. By keeping the particle number constant, the statistics of calculations will be kept in same level for all DEM simulations. The first comparison between the Eulerian and DEM simulations has been made for the decay of the relative velocity of phases |ul - um| in time, which is shown in Figure 4. The results represent three phase volume fractions (εl + εm), namely 0.01, 0.05, and 0.10, where εl ) εm. The initial value of relative velocity is 1 m/s in all simulations for which the particles density is 2500 kg/m3 with the diameter of 5 mm. The total number of particles is 104, which shares 5000 particles for each phase. The restitution coefficient is assumed as 0.80, which along with the density signifies the properties of glass. For both Eulerian and DEM simulations, the relative velocity of phases starts to decay sharply in the beginning of interaction between phases.
5274
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010 Table 2. The Parameters Corresponding to Simulations in Figure 4
Figure 5. (a) Variation of the number of particles in the CV (NCV) scaled by the initial value of particles (N0) with time. Solid, dashed, and dotted lines represent the entire volume fractions (εl + εm) of 0.01 (N0 ) 1691), 0.05 (N0 ) 883), and 0.10 (N0 ) 1564), respectively. (b) Variation of phase volume fraction with time. Open circles, filled circles, and diamonds are symbols representing Eulerian results by FLUENT for εl ) εm ) 0.050, 0.025, and 0.005, respectively. Solid and dashed lines represent the phase volume fractions (εl and εm) from DEM simulations in each set. Other simulation parameters are as described in Figure 4.
The predictions of the two approaches for the decay of relative velocity coincide very well up to about 0.8 m/s, which is 80% of the initial value. Then the deviation of results starts and grows as time moves forward. The predictions of the Eulerian method always remain above the predictions of the DEM model. This may be attributed to the structural changes that occur in particle arrangements with time in particulate models including DEM, which is contrary to the assumption of uniform position distribution for particles in Eulerian models. Such structural changes can increase the momentum transfer rate in the DEM simulations with respect to that estimated in the Eulerian model. The results shown in Figure 4 indicate that the Eulerian and DEM models are in good agreement in short-term. In longterm, the uniform conditions within the CV may alter, which creates considerable differences between the models. Figure 5a displays the changes associated with the number of particles inside the CV in time. The particle number (NCV) is normalized by the initial number of particles (N0) inside the CV that makes the normalized quantity NCV/N0 start from 1 in all cases. Figure 5a shows that the number of particles inside the CV varies below 10% of the initial value in higher volume fractions (0.05 and
parameter
value
restitution coefficient, e particle diameter, d (mm) material density, F (kg/m3) geometric sizes, a ) b (m) geometric size, h (m) geometric size, h (m) Young’s modulus, E (GPa) Poisson’s ratio, σ dynamic friction coefficient, µf
0.80 5 2500 (glass) 0.075 0.5 (εl + εm ) 0.010) 0.05 (εl + εm ) 0.050, 0.100) 70 0.2 0.3
0.10), where it behaves smoothly in the low volume fraction (0.01). Here, the second comparison between the Eulerian and DEM simulations have been made in Figure 5b, which shows the variation of phase volume fraction with time in the three above-mentioned volume fractions. The phase volume fraction from the Eulerian simulations remains strictly constant, while it varies in DEM simulations especially for the phase volume fractions of εl ) εm ) 0.05. For εl ) εm ) 0.005, the DEM shows very little variations within the time t < 0.1 s, though beyond this time significant variations of ε are likely in all cases. In summary, the results of Figure 5b can confirm that the phase volume fractions remain fairly constant in all simulations within the time 0.1s. Note that a characteristic time for the CV can be defined as h/V0, which is 0.5 s for εl ) εm ) 0.005 and 0.05 s for εl ) εm ) 0.025 and 0.050. It represents the time that a particle moves from one end to the other end of the CV (h) with the initial velocity V0. Obviously, the real time for an average particle to displace by h is longer than h/V0. Thus the time t ) 0.1 s is practically shorter than or as same as the time for a particle to displace as much as the length of the CV. Finally, the rate of momentum transfer is calculated from simulations which are compared in Figure 6. The net momentum transfer rate between two granular phases (Fml) has been calculated in three ways as follows. First, it is calculated directly from the DEM simulations using eq 18 (set A), which is symbolized by circles. Each circle stands for the instantaneous value of Fml averaged over 30 cases of DEM simulations started with different initial configurations of randomly distributed particles. The other ways of calculating Fml rely on eq 15 in which the phase volume fractions and the relative velocity between phases are taken either from Eulerian simulations by FLUENT (set B), or from the DEM simulations (set C). Figure 6a shows the results for εl ) εm ) 0.005 in which all three calculations are very close within t < 0.1 s. At the end of this period, the value of Fml is reduced to about one-third of the initial value. The fluctuation of data corresponding to set A covers the discrepancy of predictions by eq 15 between sets B and C. The results of set B are slightly above those of set C. Figure 6b represents the results for εl ) εm ) 0.025 that indicate the set A predictions of Fml are slightly below those calculated for sets B and C. However, all three sets are still in good agreement. These results are presented within a period of time (t < 0.02 s) that the initial value of Fml is decreased by two-thirds. The results corresponding to εl ) εm ) 0.050 are displayed in Figure 6c, which reveals some clear deviation of predictions between different sets. The whole time period is 0.01 s within which Fml has lost about two-third of its initial value. In the first half of this period, all sets decay closely, whereas in the second half of the period sets A and B show a good agreement with each other. However, set C predicts Fml about 25% lower than sets A and B. It should be noted that the fluctuation associated with the direct calculation of Fml from the DEM simulations is high especially as the volume fraction
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
5275
through the instantaneous distribution of momentum for the phases individually or collectively. In this context, the instantaneous directional particulate momentum per volume of the cell containing particles is defined as Nrz
∑mV
i i d
Pd(r, z) )
i)1
2πr∆r∆z
(19)
in the DEM simulations, where it gives the calculation result for a circular element of radius r at the level z, with the thickness ∆r and the height ∆z, sketched in Figure 3b. Note that Nrz represents the number of particles in the element that can contain particles from either all phases or each phase; Vid is the velocity component of particle i in direction d. For the Eulerian simulations, it can be calculated as Nphase
Pd(r, z) ) F(
∑εu) k k d
(20)
k)1
Figure 6. Variation of the magnitude of net momentum transfer rate between granular phases with time. (a) εl ) εm ) 0.005 (δt ) 0.002 s), (b) εl ) εm ) 0.025 (δt ) 0.0004 s), and (c) εl ) εm ) 0.050 (δt ) 0.0004 s). Circles represent the results obtained from eq 18 which are averaged over 30 different cases of DEM simulations. Solid and dash lines are predictions from eq 15 using the instantaneous values of εl, εm, and |ul - um| from DEM simulations, and Eulerian simulations in FLUENT, respectively. Other parameters are as given in Figure 4.
rises, which can be also a result of smaller δt used in higher volume fractions. Lateral Momentum Transfer between Solids. The simulation set for investigation of the interpenetration of two granular media was described in section 3 for the setup demonstrated in Figure 3. The lateral transfer of momentum can be studied
where Nphase indicates the number of phases to be involved in the calculation of Pd. Here, ukd stands for the velocity component of phase k within the given element in direction d. In a small control volume, the calculated values of ε1 + ε2, Pz, and Pr from the DEM simulations may fluctuate considerably from a control volume to another because there may not be sufficient number of particles in each control volume. Therefore, a proper control volume of larger size can provide statistically reliable values. Because of the cylindrical shape of the system, the circular element (ring) shown in Figure 3b can be the proper control volume for which eqs 19 and 20 are indeed written. In the Eulerian simulations, the solution is independent of the polar coordinate of the cylinder and the results can be only shown in any plane which includes z axis. Therefore, the results of DEM and Eulerian simulations for ε1 + ε2, Pz, and Pr can be visualized in rz plane, as shown in Figure 7. In this figure, the parts a, c, and e demonstrate the Eulerian simulation results while the parts b, d, and f represent the DEM simulation results. The granular phase 2 (zone 1) is penetrating into the phase 1 (zone 2) with the velocity of -0.5 m/s with a component only along z direction at the entrance. The volume fractions of both phases are initially given uniformly with the values of ε1 ) ε2 ) 0.05. These values may justify the use of laminar model in Eulerian simulations, though the changes associated with the turbulent viscous models will be examined in the future. Results are obtained by solving eqs 13 and 14 in time. Both the results of DEM and Eulerian simulations in Figure 7 are presented at t ) 0.1 s. Figure 7 panels a and b reveal the differences between DEM and Eulerian simulations for the distribution of total volume fraction. The DEM results show that the penetrating frontier is denser and more concentrated than that from the Eulerian results. The volume fraction distribution reveals that the resulting flow from the Eulerian simulations is largely scattered along the main flow direction (z axis). This is further confirmed by the distributions of Pz and Pr through Figure 7c-f. By comparing Figure 7 panels c and e to Figure 7 panels d and f, it is seen that the momentum in both r and z directions propagates widely in the Eulerian calculations while the DEM calculations display local escalation of momentum with limited scattering in space. The above-mentioned results can be also obtained separately for each phase. Figure 8 demonstrates the profiles of ε1 + ε2, Pz, and Pr at t ) 0.1 s calculated for the two circular elements corresponding to the inlet level (z ≈ 0), and the downstream level z ) -0.020 m (|z|/d ≈ 20), which is considered as the core of the penetration frontier. Figure 8 panels a and b represent
5276
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
tively. For the inlet level shown in Figure 8a, the phase volume fractions are predicted from the two models with some contradictive trends for phase 1. These differences in solid volume fraction of phases are further manifested in Figure 8c,e via the profiles of Pr and Pz. It can be seen in Figure 8c that phase 1 has gained some momentum in radial direction which has propagated farther in the DEM calculation than in the Eulerian simulation. This may explain why the trend of the volume fraction of phase 1 from DEM results contradicts with that from the Eulerian results in Figure 8a. A similar argument can be made on the results shown in Figure 8e in which the momentum gain of phase 1 in the axial direction is propagated farther in the DEM calculation than in the Eulerian simulation. Note that Figure 8e does not include the profiles of Pz for phase 2 as their magnitudes are significantly larger than those for phase 1. In the meantime, the results for calculation of Pz for phase 2 from the DEM and Eulerian models do not differ notably since the interaction between the two phases does not happen at the inlet level by the time of 0.1 s. For the downstream core level shown in Figure 8b, the DEM simulations predict considerably higher volume fractions for both phases than those predicted by the Eulerian simulations. Despite the significant quantitative differences, the profiles predicted by the two models follow similar trends. The volume fraction of phase 1 is about the bulk value (0.05) at the center (r ≈ 0), which increases to a maximum around r ≈ 0.01 m (r/d ≈ 10) in both models. Then it decays gradually to the bulk value at r ≈ 0.02 m (r/d ≈ 20). The volume fraction for phase 2 decays monotonically from the center to outer radii, which depicts the same width in the two models though the calculated values from DEM simulations are twice larger than the Eulerian results. The profiles of radial and axial momenta per volume (Pr and Pz) are shown for the downstream core level in Figure 8d,f. Pr from the DEM calculations can be up to about three times larger than that calculated from the Eulerian model for both phases while they follow same trends. In other words, the width of the radial momentum profile is the same for both the models but the DEM yields larger values. This observation indicates that the Eulerian model employed in this study underestimates the momentum transfer rate between the phases in transverse (radial) direction as compared to the DEM estimations. Moreover, the profiles of axial momentum per volume, as shown in Figure 8f, display considerably lower estimations of Pz in the Eulerian model. As the considered level represents the location of maximum transfer of momentum between the phases, such differences could indicate that the mechanism of momentum transfer modeled in the Eulerian formulations used here is propagating momentum far more than transferring it between the phases. Figure 7. Visualization of various quantities averaged within the circular rings at 0.1 s after starting the interpenetration of granular phases using the DEM and the Eulerian simulations. In both simulations, we set d ) 1 mm and D ) 0.04 m (D/d ) 40), and ds ) 0.01 m (ds/d ) 10). The entering velocity of phase 2 to the phase 1 across the inlet of zone 2 is set to 0.5 m/s. Also, the restitution coefficient is set to 0.9. Initial volume fraction of each phase was set to 0.05. Other material properties are as explained in Figure 4. The Eulerian model settings are given in FLUENT according to Table 1. Distribution of total solid volume fraction (ε1 + ε2) using (a) Eulerian, and (b) DEM simulations with the corresponding color bar. Distribution of Pz using (c) Eulerian, and (d) DEM simulations with the corresponding color bar. Distribution of Pr using (e) Eulerian, and (f) DEM simulations with the corresponding color bar. The coordinates r and z are demonstrated in part a. The extent of the domain in all parts of the figure is from 0 to 0.025 m for r direction (0 < r/d < 25), and from -0.05 m to 0 (-50 < z/d < 0) for z direction.
the volume fraction of the two phases calculated from both methods in the inlet and the downstream core levels, respec-
5. Conclusions The momentum transfer between solids is investigated by the Eulerian approach as the continuum description of granular matter, which is verified by the DEM as a particulate model. The physical systems under study were set in a way that the first system provided a spatial average of the instantaneous momentum transfer rate between two granular phases and the second system displayed spatial distribution of instantaneous momentum for each granular phase. In the first system, the momentum transfer rate is directly calculated from the DEM simulations and is compared to the predictions of the momentum transfer rate using the Syamlal’s formula.6 The results of the DEM and the Eulerian models for instantaneous relative phase velocity and phase volume fractions are employed in the Syamlal’s formula. The results revealed that both methods
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
5277
Figure 8. Instantaneous profiles of solid volume fractions of phases in (a) the inlet level (z ≈ 0), and (b) the downstream level (z ≈ -0.02 m) corresponding to the core of penetration. Instantaneous profiles of Pr of phases in (c) the inlet level (z ≈ 0), and (d) the downstream level (z ≈ -0.02 m). Instantaneous profiles of Pz of phases in (e) the inlet level (z ≈ 0), and (f) the downstream level (z ≈ -0.02 m). All other settings of simulations are as described in Figure 7.
produce close enough estimations in very low volume fraction (0.01 tested). As the volume fraction rises to 0.05 and 0.10, the deviation increases between the direct calculation of momentum transfer rate in the DEM simulations and the predictions of the Syamlal’s formula. However, the deviation can be either positive or negative, which means the Syamlal’s formula predicts above and below the DEM direct results in the volume fractions of 0.05 and 0.10, respectively. It is also worth mentioning that the relative velocity of phases calculated from the DEM and Eulerian methods starts to deviate when it is dropped to about 80% of the initial value regardless of the phase volume fraction. Once their deviation grows, the calculated values from the Eulerian model are always higher than those from the DEM model. In other words, the two phases come to the same velocity slower in the Eulerian model than the DEM model. Thus one may conclude the existence of larger drag force between the phases in the DEM model. However, the Syamlal’s formula shows the opposite according to Figure 6, that is, the drag force is lower in the DEM model. This leads us to conclude that in the Eulerian models the relative velocity
of phases may be boosted more by other quantities than what is decreased by the net momentum transfer rate between granular phases. Further investigation is required to determine such quantities and other factors which influence the relative phase velocity in both the DEM and the Eulerian models. In the second system, the penetration of a moving granular phase into a stationary granular phase was simulated using both the DEM and the Eulerian models. The purpose of making this simulation was to compare the spatial distributions of instantaneous momenta between the two models. Our results revealed that there are certain differences in spatial distributions of volume fraction and momenta so that the DEM produces less scattered configurations while the Eulerian model predicts the granular phases dispersed in a wider region in space as a result of interpenetration. As a result, particles momenta distributions, and thus phase volume fractions, from the DEM simulations are more intense in the penetrating front while they are dilated elsewhere. However, momenta and volume fractions are widely spread within the interacting region of the phases. This leads us to the same conclusion as in the first case that the Eulerian
5278
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
model used in this study underestimates the rate of momentum transfer between granular phases, which can lead to lower energy dissipation and higher scattering of phases. The findings in this paper require further investigation to the roles of different factors influencing the rate of momentum transfer between granular phases in Eulerian models. These factors include the viscous model type, combination of various settings given in Table 1, and the grid size. For the DEM simulations, the future challenges will be basically associated with the higher number of particles through numerous assemblies, as well as applying more advanced models for particle interactions. Acknowledgment Authors would like to acknowledge the financial support of the Academy of Finland under Grant No. 123938 and 124368. In addition, authors gratefully acknowledge Prof. D. Gidaspow’s valuable speech, comments and insights during the summer school in heat and mass transfer 2008 held at Lappeenranta University of Technology in Finland that inspired us to carry out this work. Supporting Information Available: The distributions of phase volume fraction in time from the DEM and Eulerian simulations are presented in video animations. This information is available free of charge via the Internet at http:// pubs.acs.org. Nomenclature a ) acceleration vector, m/s2 A ) width of bed, m B ) width of bed, m d ) particle diameter, m ds ) diameter of granular stream, m D ) diameter of granular medium, m e ) coefficient of restitution E ) Young’s modulus, Pa FCi ) net contact force on particle i, N FCnij ) normal contact force between particles i and j, N FCtij ) tangential contact force between particles i and j, N Fk ) net interaction force on phase k, N Fml ) net momentum transfer rate between phases m and l, N/m3 g ) gravitational acceleration vector, m/s2 g0 ) radial distribution function at contact G ) shear modulus, Pa H ) height of bed, m I ) moment of inertia, kgm2 Kn ) normal stiffness, Pa m1/2 kt ) tangential stiffness, Pa m l ) length of granular stream in DEM simulations, m L ) length of granular medium, m m ) particle mass, kg m j ) reduced mass, kg M ) number of solid phases nij ) unit vector along the centerline of particles i and j Nml ) number of particles in CV Nrz ) number of particles in a circular element Nphase ) number of phases in an Eulerian grid N0 ) initial number of particles in CV p ) hydrostatic pressure, Pa Pd ) particulate momentum per unit volume, kg m-2 s-1 r ) radial position of an element, m R ) particle radius, m Sml ) solid-solid drag coefficient, kg m-3 s-1
tij ) tangential unit vector at contact point between particles i and j t ) time, s TCi ) net torque on particle i, N m u ) velocity vector of continuous phase k, m/s v ) velocity vector of particle, m/s V0 ) magnitude of initial velocity of moving solid phase, m/s vrij ) relative velocity vector between particles i and j, m/s vsij ) slip velocity vector between particles i and j, m/s Greek Letters R ) coefficient relating η to e r ) rotational acceleration, rad/s2 δtij ) tangential deformation vector, m δnij ) normal deformation magnitude, m δt ) averaging time for momentum transfer between phases, s ∆r ) radial thickness of an element, m ∆z ) longitudinal thickness of an element, m ε ) phase volume fraction ηn ) normal damping factor, kg/s ηt ) tangential damping factor, kg/s µf ) dynamic friction coefficient F ) density, kg/m3 σ ) Poisson’s ratio σ c ) phase stress tensor, Pa ω ) angular velocity, rad/s Subscripts d ) spatial direction index (x, y, z) i ) particle index j ) particle index k ) phase index l ) solid phase index m ) solid phase index Superscripts ′ ) postcollisional sign i ) particle index k ) phase index
Literature Cited (1) Tsuji, Y. Multiscale modeling of dense phase gas-particle flow. Chem. Eng. Sci. 2007, 62, 3410–3418. (2) Tsuji, T.; Yabumoto, K.; Tanaka, T. Spontaneous structures in threedimensional bubbling gas-fluidized bed by parallel DEM-CFD coupling simulation. Powder Technol. 2008, 184, 132–140. (3) 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. (4) Cundall, P. A.; Strack, O. D. L. A discrete numerical model for granular assemblies. Geotechnique 1979, 29, 47–65. (5) Gidaspow, D.; Jung, J.; Singh, R. K. Hydrodynamics of fluidization using kinetic theory: An emerging paradigm 2002 Flour-Daniel lecture. Powder Technol. 2004, 148, 123–141. (6) Syamlal, M. The Particle-Particle Drag Term in a Multiparticle Model of Fluidization; Topical Report, DOE/MC/21353-2373, NTIS/DE87006500; National Technical Information Service: Springfield, VA, 1987. (7) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon Press: Oxford, 1987. (8) Zhu, H. P.; Zhou, Z. Y.; Yang, R. Y.; Yu, A. B. Discrete particle simulation of particulate systems: Theoretical developments. Chem. Eng. Sci. 2007, 62, 3378–3396. (9) Po¨schel, T.; Schwager, T. Computational Granular Dynamics; Springer: Berlin, Heidelberg, 2005. (10) Mindlin, R. D. Compliance of elastic bodies in contact. J Appl. Mech. 1949, 16, 259–268.
ReceiVed for reView October 2, 2009 ReVised manuscript receiVed November 24, 2009 Accepted November 30, 2009 IE901538W
