Ind. Eng. Chem. Res. 2008, 47, 481-485
481
RESEARCH NOTES Master Curve for the Discrete-Element Method Eldin Wee Chuan Lim* DSO National Laboratories, 20 Science Park DriVe, Singapore 118230
The discrete-element method (DEM) is a numerical model consisting of the fundamental laws of motion and a suitable force-displacement model, such as the linear spring-and-dashpot model, as closure. The damping coefficient found in the latter is related to the phenomenological concept of coefficient of restitution of solid particles. However, the relationship between these two coefficients is strongly dependent on various material properties of the solid particles to be simulated, and so its determination often involves a numerical trialand-error procedure. This paper addresses this problem by presenting a master calibration curve and equation that can be used to calculate the value of the damping coefficient directly from values of the coefficient of restitution and other material properties. The Buckingham Pi theorem was applied to derive a functional relationship between the two coefficients and other material properties, which allowed a large number of calibration curves to be collapsed onto a single master calibration curve. Introduction The discrete-element method (DEM) has, in recent years, become a viable complement and even alternative to experimental research in such areas as granular physics, particulate science, and powder engineering. In particular, with the advent of powerful computing facilities in many parts of the world, it has now become possible to perform accurate computer simulations of various types of granular flow processes or solid-fluid multiphase systems in order to gain a better understanding of the mechanics involved at either the process or phenomenon level.1-3 In recent years, the technique of combining DEM with computational fluid dynamics (CFD) for numerical studies of multiphase flow systems has also gained much popularity.4-14 DEM was originally developed by Cundall and Strack15 for describing the mechanical behavior of assemblies of discs and spheres. The translational and rotational motions of individual solid particles are governed by Newton’s laws of motion,
dvi mi
Following the originators of DEM, a linear spring-and-dashpot model may be applied as a closure to the governing equations. With such a closure, interparticle collisions are modeled as compressions of a perfectly elastic spring while the inelasticities associated with such collisions are modeled by the damping of energy in the dashpot component of the model. Collisions between particles and a wall may be handled in a similar manner but with the latter not incurring any change in its momentum. In other words, a wall at the point of contact with a particle may be treated as another particle but with an infinite amount of inertia. Contact and damping forces are, thus, calculated according to the following equations,
fcn,ij ) -κn,iδn,ij
(3)
fct,ij ) -κt,iδt,ij
(4)
fdn,ij ) -ηn,i(Vr‚ni)ni
(5)
fdt,ij ) -ηt,i{(Vr‚ti)ti + (ωi × Ri - ωj × Rj)}
(6)
N
)
dt
(fc,ij + fd,ij) + mig ∑ j)1 dωi
Ii
dt
(1)
N
)
∑ j)1
Tij
(2)
where mi and Vi are the mass and velocity of particle i, N is the number of particles in contact with this particle, fc,ij and fd,ij are the contact and viscous contact damping forces, respectively, Ii is the moment of inertia of particle i, ωi is its angular velocity, and Tij is the torque arising from contact forces, which will cause the particle to rotate. The relationships between contact and damping forces with relative displacements and velocities of particles have to be determined using an appropriate force-displacement model. * To whom correspondence should be addressed. Tel.: (65) 6871 2967. Fax: (65) 6873 0742. E-mail:
[email protected].
where fcn,ij, fdn,ij and fct,ij, fdt,ij are the normal and tangential components of the contact and damping forces, respectively; κn,i, δn,ij, ni, ηn,i and κt,i, δt,ij, ti, ηt,i are the spring constants, displacements between particles, unit vectors, and damping coefficients in the normal and tangential directions, respectively; Vr is the relative velocity between particles; and Ri and Rj are the radius vector (from particle center to a contact point) for particles i and j, respectively. If |fct,ij| > |fcn,ij| tan φ, then “slippage” between the two contacting surfaces is simulated by a Coulomb-type friction law, |fct,ij| ) |fcn,ij| tan φ, where tan φ is analogous to the coefficient of friction. Other than the linear spring-and-dashpot model applied in the present study, other nonlinear force-displacement models include that derived based on a Hertzian theory and that used by Zhou et al.1 in their study of particle heap formation and the relationships between angle of repose, rolling friction coefficient, and particle size.
10.1021/ie0709032 CCC: $40.75 © 2008 American Chemical Society Published on Web 12/15/2007
482
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008
One important feature of DEM in comparison with so-called hard-sphere models is that damping forces are calculated rigorously with the spring-and-dashpot model, giving rise ab initio to the phenomenological concept of the coefficient of restitution of particles, instead of it being imposed arbitrarily as a parameter. Hard-sphere models are essentially event-driven algorithms that permit only binary, instantaneous collisions between particles. The advantage of such algorithms is the possibility of larger temporal and spatial scale simulations, which are necessary for the study of flow instabilities.16 In DEM, the damping coefficient determines the amount of energy lost during a collision in analogy with the way energy is damped in a dashpot and so affects the elasticity of particle-particle or particle-wall collisions. The coefficients of restitution for these two types of collisions are material properties of both particles and wall and may be obtained experimentally. However, values of coefficients of restitution have to be “translated” to equivalent values of damping coefficients before DEM can be applied. In other words, the damping coefficient becomes a parameter for which an appropriate numerical value has to be assigned prior to a DEM simulation. Cundall and Strack15 proposed expressions relating the damping coefficient to the value of the spring constant. These were derived based on the condition of critical damping of a system consisting of a mass, m, spring, and dashpot with a single degree of freedom:
ηn ) 2 xmκn
(7a)
ηt ) 2 xmκt
(7b)
Tsuji et al.4 were the first to recognize the correspondence between the damping coefficient in the force-displacement model and the coefficient of restitution of real particles. They derived an expression through a heuristic approach, which relates the damping coefficient to the particle mass, spring constant, displacement, and coefficient of restitution,
ηn ) R xmκn δn,ij1/4
(8)
where R is an empirical parameter correlated with the coefficient of restitution through a graphical relationship by the authors. Tsuji et al.5 also derived an expression relating the two coefficients by solving the equation of motion for an oscillation system consisting of a mass, spring, and dashpot,
η ) 2γ xmκ
(9)
where γ can be determined from the coefficient of restitution, , as γ ) R/(x1+R2) and R ) -ln /π. Another common approach used for the calculation of the damping coefficient is to conduct what is referred to as a numerical experiment to obtain a value for the damping coefficient, which corresponds to a desired value of the coefficient of restitution for the specific type of particles to be simulated.6,7,10-14 In the present study, this procedure will be referred to as calibrating the damping coefficient in the springand-dashpot model of DEM to the phenomenological coefficient of restitution of physical, solid particles. This involves a trialand-error procedure whereby single particles with assigned values of damping coefficient are simulated to fall from rest and rebound from a solid surface to obtain the corresponding value of coefficient of restitution. The advantage of this approach lies in its capability to calculate the value of the coefficient of restitution from that of the damping coefficient through a direct
numerical simulation of the actual physical phenomenon of particle-particle or particle-wall collision, which provides the fundamental basis for the definition of the former. The damping coefficient so obtained is then a direct and first principle correspondence of the coefficient of restitution, the accuracy of which is only limited by that of the direct numerical simulation performed. However, the disadvantage of using this approach for calculating the value of the damping coefficient is the iterative nature of the numerical procedure involved in its implementation. It is the primary objective of this study to construct a universal calibration curve, which may be used by all future DEM users for carrying out such calibrations without having to go through this trial-and-error numerical procedure. Numerical Method The most fundamental method of determining a numerical value for the damping coefficient that corresponds to a desired value for the coefficient of restitution involves simulating a single particle to fall from rest under the effects of gravity and the rebound from a flat surface. The coefficient of restitution is then given by xh/ho, where ho and h are the maximum heights of the particle above the surface before and after the collision, respectively. This type of simulation is repeated with different values of the damping coefficient until the required value of coefficient of restitution is obtained from the calculation to a desired level of accuracy. Generally, the higher the damping coefficient, the lower is the coefficient of restitution. However, as will be shown in the following section of this Research Note, there exists a nonlinear relationship between the two coefficients, which is strongly dependent on other parameters used in the simulation. Specifically, these include material-related properties of the solid particles simulated such as the value of the spring constant, which determines the stiffness of the spring used in modeling particle-particle or particle-wall collisions and so is an indication of the hardness of the colliding entities and the solid density and diameter of the particles. As such, and has been implied in the Introduction, it becomes necessary for each DEM user to carry out his or her own calibration for the particular type of solid particles that are to be simulated. To the knowledge of the author, there has been no attempt to address this issue in the scientific literature to date. In the present study, the above calibration procedure was repeated for values of the spring constant, solid density, and particle diameter in the range 103-105 N m-1, 500-5000 kg m-3, and 10-410-2 m, respectively, and the relationship between the damping coefficient and the coefficient of restitution for each case was recorded. The ranges of values considered for the material properties were deemed sufficient to include most types of noncohesive particulate materials used in experimental or numerical research. The initial height, ho, from which a particle was simulated to undergo freefall determines its impact velocity upon reaching the flat surface. However, this was found to have negligible effect on the relationship between the two coefficients for ho in the range 0.1-10 m. For example, for a particle with κ ) 105 N m-1, Fp ) 1000 kg m-3, and D ) 1 mm, the coefficient of restitution varies in the range 0.940.96 for impact velocities in the range 1.4-14.0 m s-1. Without loss of generality, the value of ho was fixed at a constant value of 1 m for all simulations conducted in this study. The time step used in all cases was 10-8 s. This has also been verified during the preliminary phase of this study to be sufficiently small such as to cause negligible numerical errors even with the stiffest particles simulated.
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008 483
Figure 2. Master calibration curve for relating the coefficient of restitution, , to the dimensionless group, Π ) η/xFpκD3. A total of 19 145 data points obtained from conducting the numerical calibration experiment in the present study are plotted.
using DEM. This explains why each DEM simulation has to be preceded by the numerical experiment to calibrate the damping coefficient to the coefficient of restitution described previously. It would be desirable to construct a calibration curve that can be used by all DEM users to derive the value of the damping coefficient without having to go through such a numerical procedure, regardless of the specific values of the material properties of their particles. This implies that a universal way of representing all data points on the various calibration curves as one single curve or mathematical expression is required. To this end, a dimensional analysis of the various parameters using the Buckingham Pi theorem was applied. The three material-related parameters, namely, D, κ, and Fp, were selected as repeating parameters to construct the dimensionless groups. It would then be straightforward to derive the following functional relationship,
(x )
)f
η
(10)
FpκD3
where is the coefficient of restitution and η is the damping coefficient. Figure 2 shows a plot of coefficient of restitution against the
Figure 1. Relationship between the damping coefficient and the coefficient of restitution for (a) κ ) 5 000 N m-1, Fp ) 1 000 kg m-3, and D ) 2, 4, 6, 8, and 10 mm; (b) D ) 1 mm, Fp ) 2 500 kg m-3, and κ ) 2 000, 4 000, 6 000, 8 000, and 10 000 N m-1; and (c) D ) 3 mm, κ ) 1 000 N m-1, and Fp ) 1 000, 2 000, 3 000, 4 000, and 5 000 kg m-3.
Results and Discussion Parts a, b, and c of Figure 1 show the relationships between the damping coefficient and the coefficient of restitution for different values of the particle diameter, D, spring constant, κ, and solid density, Fp, respectively, keeping the values of the other two parameters fixed in each case. It may be seen that there exists a nonlinear, nonincreasing relationship between the two coefficients, thus verifying the point made in the previous section. It is strongly dependent on materialrelated properties of the solid particles to be simulated
dimensionless group, Π ) η/xFpκD3, seen in eq 10 for all data points collected using the numerical calibration procedure in the present study. A total of 19 145 data points are presented. Remarkably, all points seem to fall on a single smooth curve. This would represent a master calibration curve, which can be used to derive the value of the damping coefficient given the material properties of the solid particles to be simulated and the desired value for their coefficient of restitution. With such a curve, it is no longer necessary to carry out a numerical calibration each time one needs to do a DEM simulation. To facilitate calculation or computer programming, it would be more convenient to have such a calibration curve in the form of a mathematical expression. Thus, by simple curve fitting of the data points in Figure 2 to an analytical function, the universal relation between the two dimensionless groups may be expressed as the following master calibration equation:
) -0.8
(x ) (x ) (x ) η
FpκD
3
3
+ 1.9
η
FpκD
2
3
- 2.1
η
FpκD3
+1 (11)
484
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008
It can be easily shown that the expression derived by Tsuji et al.5 (eq 9) for relating the damping coefficient and the coefficient of restitution may be rearranged into the following form:
(
) exp -
)
2.2η
xFpκD3 - 0.5η2
(14)
Parts a and b of Figure 3 show a quantitative comparison of the accuracy of both equations in predicting values of coefficients of restitution generated from numerical experiments. It may be seen that, in the range of 0.0-0.2, both eq 11 and eq 14 tend to underpredict the values of coefficient of restitution. On the other hand, within the range of values of 0.2-1.0, both expressions show excellent agreement with the numerically generated data. As the expression due to Tsuji et al.5 shown in eq 14 was derived by analytical solution of the equation of motion for a linear oscillation system, this provides strong validation of the adequacy of the numerical approach presented in this study for deriving such relationships between the damping coefficient and the coefficient of restitution. In cases of nonlinear systems, such as those described by the Hertzian theory mentioned earlier where no such analytical solution is possible, the present numerical method would provide a reliable and, in fact, only means of deriving such relationships. The exploration of this will be the subject of a future study. Conclusions
Figure 3. Comparison of values of coefficient of restitution generated from numerical experiments with values predicted from eq 11 (present study) and eq 14 (Tsuji et al.)5 in the range (a) 0.0-0.2 and (b) 0.2-1.0.
A further analysis of the dimensionless group, Π, in the calibration equation will now be made to gain some possible further insights to its physical significance:
Π2 )
η2 FpκD3
(12)
By applying the force-displacement model given in eq 5, it may be shown that
Π2 )
fd,ijD FpVrD 2 κD µ
(13)
where µ is a characteristic viscosity of the dashpot in the springand-dashpot model defined here as the damping coefficient, η, divided by the particle diameter, D, which is taken to be a characteristic length scale. Here, a characteristic viscous dissipation, ED, is defined as the product of the characteristic damping force, fd,ij, and the characteristic length scale, D. Similarly, a characteristic elastic potential energy, EE, is defined as the product of the spring constant, κ, and D2. The dimensionless quantity, FpVrD/µ, in eq 13 resembles a Reynolds number and will be referred to as a DEM-Reynolds number, ReDEM. Finally, it may then be observed that the dimensionless group, Π, in the calibration equation may be interpreted as a function of a ratio between two energy-related quantities and a modified Reynolds number, Π ) xE/ReDEM, where E ) ED/EE.
The numerical value of the damping coefficient of the linear spring-and-dashpot model often applied in DEM has to be supplied a priori to the actual DEM simulation. It may be associated with the more common phenomenological concept of coefficient of restitution of solid particles. However, because of the fact that the relationship between the two coefficients is strongly dependent on various material properties of the particles to be simulated, a direct association between them would, on first sight, be a difficult multidimensional, multivariable problem. This study has attempted to address this problem by presenting a master calibration curve and equation that can be used to calculate directly the damping coefficient from values of the coefficient of restitution and other material properties without going through such an iterative process. The Buckingham Pi theorem was applied to derive a functional relationship between the coefficient of restitution and a dimensionless group consisting of the damping coefficient as well as material properties of the solid particles. A large number of calibration curves depicting the various relationships between the two coefficients for a practical range of values for all material properties were collapsed onto a single master calibration curve, which could also be described by a simple analytical function. These would be useful references for future research workers doing DEM simulations. It would also be an interesting extension to the present study to develop similar calibration curves for other nonlinear force-displacement models of DEM. Nomenclature D ) particle diameter E ) energy ratio ED ) characteristic viscous dissipation EE ) characteristic elastic potential energy fc,ij ) contact force fcn,ij ) normal component of contact force fct,ij ) tangential component of contact force
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008 485
fd,ij ) viscous contact damping force fdn,ij ) normal component of viscous contact damping force fdt,ij ) tangential component of viscous contact damping force ho ) maximum height of particle above surface before collision h ) maximum height of particle above surface after collision Ii ) moment of inertia mi ) particle mass ni ) unit vector in normal direction N ) number of particles Ri ) radius vector from particle center to a contact point ReDEM ) DEM-Reynolds number ti ) unit vector in tangential direction Tij ) torque Vi ) particle velocity Greek Symbols R ) empirical parameter δn,ij ) displacements between particles in normal direction δt,ij ) displacements between particles in tangential direction ) coefficient of restitution φ ) angle of friction κn,i ) spring constant for normal collisions κt,i ) spring constant for tangential collisions ηn,i ) damping coefficient in normal direction ηt,i ) damping coefficient in tangential direction µ ) characteristic viscosity Π ) dimensionless group Fp ) particle density ωi ) angular velocity Literature Cited (1) Zhou, Y. C.; Wright, B. D.; Yang, R. Y.; Xu, B. H.; Yu, A. B. Rolling friction in the dynamic simulation of sandpile formation. Physica A 1999, 269, 536. (2) Ning, Z.; Ghadiri, M. Distinct element analysis of attrition of granular solids under shear deformation. Chem. Eng. Sci. 2006, 61, 5991.
(3) Jayasundara, C. T.; Yang, R. Y.; Yu, A. B.; Curry, D. Discrete particle simulation of particle flow in the IsaMill process. Ind. Eng. Chem. Res. 2006, 45, 6349. (4) Tsuji, Y.; Tanaka, T.; Ishida, T. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technol. 1992, 71, 239. (5) Tsuji, Y.; Kawaguchi, T.; Tanaka, T. Discrete particle simulation of two-dimensional fluidized bed. Powder Technol. 1993, 77, 79. (6) Xu, B. H.; Yu, A. B. Numerical simulation of the gas-solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics. Chem. Eng. Sci. 1997, 52, 2785. (7) Xu, B. H.; Yu, A. B.; Chew, S. J.; Zulli, P. Numerical simulation of the gas-solid flow in a bed with lateral gas blasting. Powder Technol. 2000, 109, 13. (8) Li, J.; Webb, C.; Pandiella, S. S.; Campbell, G. M.; Dyakowski, T.; Cowell, A.; McGlinchey, D. Solids deposition in low-velocity slug flow pneumatic conveying. Chem. Eng. Process. 2005, 44, 167. (9) Pandit, J. K.; Wang, X. S.; Rhodes, M. J. On geldart group A behaviour in fluidized beds with and without cohesive interparticle forces: A DEM study. Powder Technol. 2006, 164, 130. (10) Lim, E. W. C.; Wang, C. H.; Yu, A. B. Discrete element simulation for pneumatic conveying of granular material. AIChE J. 2006, 52, 496. (11) Lim, E. W. C.; Zhang, Y.; Wang, C. H. Effects of an electrostatic field in pneumatic conveying of granular materials through inclined and vertical pipes. Chem. Eng. Sci. 2006, 61, 7889. (12) Lim, E. W. C.; Wang, C. H. Diffusion modeling of bulk granular attrition. Ind. Eng. Chem. Res. 2006, 45, 2077. (13) Lim, E. W. C.; Wong, Y. S.; Wang, C. H. Particle image velocimetry experiment and discrete-element simulation of voidage wave instability in a vibrated liquid-fluidized bed. Ind. Eng. Chem. Res. 2007, 46, 1375. (14) Lim, E. W. C. Voidage waves in hydraulic conveying through narrow pipes. Chem. Eng. Sci. 2007, 62, 4529. (15) Cundall, P. A.; Strack, O. D. L. A discrete numerical model for granular assemblies. Geotechnique 1979, 29, 47. (16) Conway, S. L.; Liu, X.; Glasser, B. J. Instability-induced clustering and segregation in high-shear Couette flows of model granular materials. Chem. Eng. Sci. 2006, 61, 6404.
ReceiVed for reView July 2, 2007 ReVised manuscript receiVed August 26, 2007 Accepted November 29, 2007 IE0709032
