Numerical Prediction of the Bulk Density of Granular Particles of Different Geometries

The Discrete Element Method (DEM) is increasingly used to simulate the motion of granular material in engineering devices. The most widely used shape to model granular particles is the sphere due to its simplicity which facilitates computations. However, with increasing advancement and popularity of the DEM approach, a more accurate representation of real granular shapes is becoming intensely important. This paper is an attempt to predict the bulk density of granular particles of cylindrical shape by DEM. The method used to model tube-like particles is the composite method where the tube-like particle is composed of a number of spherical sub-shapes. The results show that the DEM prediction of the bulk density of tube-like particles were accurate provided that a sufficient number of spherical sub-shapes was used to construct the particle. In addition, the bulk density of granular particles in the shape of 32-face polyhedrons was also investigated. It was shown that 32-face polyhedrons can be approximated as spheres provided that the rolling resistance of the shape is taken into consideration.


Introduction
Granular materials are complex systems of a large number of particles of various sizes, shapes and materials. Under different circumstances they reveal different behavior such as solid-, liquid-or gas-state behavior. In fact, the science of granular flow is not yet as well understood and well developed as other classes of materials (Pöschel and Schwager, 2005;Campbell, 2006). In order to predict and optimize the behavior and motion of granular matter in engineering devices, numerical simulation tools are increasingly employed (Cleary, 2004). The continuous increase in computing power is now enabling researchers to implement numerical methods that deduce their global characteristics from analyzing the individual behavior of each grain. This forms the foundation of the Discrete Element Method (DEM) which is, to date, the leading approach to simulate the dynamics of granular media.
However, predicting the behavior of such a complex physical system by computer models has its own challenges. One of the main challenges of DEM is in modeling the complex shapes which the granular particles often exhibit in reality. The most widely used particle shape in DEM is the sphere due to its simplicity. With increasing advancement and popularity of the DEM approach, the expectations of its capabilities rise, too. More realistic modeling of granular particles of various geometries is becoming increasingly necessary.
Different methods already exist to model complex shapes in DEM which will be briefly noted here. More complete reviews of different methods in DEM to model non-spherical shapes are reported by Hogue (1998), Latham andMunjiza (2004), Haff (1993) and Ristow (1994).
One approach is to model the complex shape as a cluster of spherical sub-shapes. Different variations of the so-called composite method can be found in the literature. The sub-shapes can be connected to one another by damped springs (Buchholtz and Pöschel, 1994;Pöschel and Buchholtz, 1993;Maeno, 1996;Lu and McDowell, 2007), or else the aggregate shape can be considered as rigid (Gallas and Sokolowski, 1993;Nolan and Kavanagh, 1995;Mutze, 2006). The sub-shapes used to construct the complex shape might also overlap each other. The subshapes do not necessarily need to be spheres; for example, triangular sub-shapes have been developed which are useful to model shapes with sharp edges (Pöschel and Buchholtz, 1995), and arc sub-shapes can be used to construct oval shapes (Wang et al., 1999). Another approach to model non-spherical shapes is to address the shape as it is and develop analytical methods for the detection of contacts. Such analytical methods are already developed for some † shapes such as ellipsoids (Lin and Ng, 1995;Lin and Ng, 1997;Allen et al., 1989;Perram et al., 1984;Perram and Wertheim, 1985;Ouadfel and Rothenburg, 1999;Johnson et al, 2004) and polyhedrons (Cundall, 1988;Zhao et al., 2006). Superquadrics are also proposed (Barr, 1981) and used in DEM studies, e.g. superellipsoids (Wellmann et al., 2008;Cleary et al., 1997). A different method is the virtual space method where the space is divided into discrete voxels (Munjiza, 2004;Jia and Williams, 2001). Each particle takes up some space in the digitized space for which the corresponding voxels would be identified. Two particles are in contact if they both occupy the same voxel. Each of these methods has its own advantages and disadvantages but a common challenge remains the computational costs and complexity of modeling complex shapes. The most widely accepted method, at least in commercial software, is the composite approach. However, according to Kruggel-Emden et al. (2008), the validity of the method has not been sufficiently tested and verified. This paper, therefore, is an attempt to further study the validity of the composite approach and to investigate its capabilities and limitations.
The application chosen to verify the method is the poured bulk density of granular particles. The dependency of the bulk density on the shape of the individual particles has long been recognized (Riley and Mann, 1972;Fonner et al., 1966;Ridgway and Rupp, 1969), which makes the case suitable for validating the composite method.

DEM study
Specific boundary conditions for the DEM study are presented for each case in Table 1. The friction and restitution coefficients are based on simple lab tests conducted on sample particles, and are verified by the available data in the literature (Engineers handbook; Engineering toolbox; Bennett and Meepagala, 2006;Beardmore, 2012). The time step was fixed at 10 -5 and the end time of the simulation was set to 8 seconds when all particles come to rest and the kinetic energy tends to zero. Since the size of the container in all simulations and experiments was chosen to be identical, it suffices to compare the height of the packed bed. The height of the packed bed is calculated as the distance of the highest point of the packed bed after complete settlement to the bottom of the container. It is worth noting that the height of the packed bed in each experiment reported here is not comparable because in each case, a different number of particles and sometimes different size of container is used. The objective is therefore to compare the numerical results with the experiments in three independent cases presented in Table 1.
The rolling friction is included in order to account for the resistance of the spheres against rolling due to small deviations from the spherical shape. An assumption taken regarding the shape of the tube-like particles is that the hollowness of the tubes was neglected in the DEM study, i.e. the particles were assumed to be full. The explanation is that the hollowness of the tubes does not affect the packing of the particles for small hole diameters because other particles cannot penetrate through the inner hole or the penetration is negligible. Another simplifying assumption is made regarding the Young modulus of the moving particles, whereas the Young modulus is reduced one order of magnitude from the original value. Conducting a sensitivity analysis confirmed the findings of several other studies (Kawaguchi et al., 1997;Ketterhagen et al., 2007;Ketterhagen et al., 2005;Tuley et al., 2010;Hogue and Newland, 1994;Pérez, 2008;Silbert et al., 2001;Kruggel-Emden et al., 2007) that lowering the stiffness within a certain range does not affect the results of the DEM simulation. This is due to the fact that although changing the stiffness constant does affect the duration of a single collision, the velocity and trajectory of the particle after the collision is only dependent on the coefficient of restitution and the coefficients of friction. However, the particle's hardness should be large enough so that the maximum overlap does not exceed the particle's radius, in which case the particles pass through each other. Without this assumption, the time step should have been reduced to accurately resolve the impact of small sub-shapes used to construct the tubular particles. This assumption is therefore necessary to reduce the high computational costs of the simulation of tube-like particles. The tube-like particles are modeled by the composite method as an aggregate of spherical sub-shapes. In this case, each tube is represented by 20 spheres which cluster and form the cylindrical shape. The same coefficient of rolling friction as spheres is assigned to the cylinders when they roll around the principal axis of rotation. The 32-face polyhedrons used in the experimental study (Gan et al., 2004;Jia et al., 2007) are quite similar in shape to spheres and are therefore represented by spheres in this work. The coefficient of rolling friction is raised slightly for the polyhedral particles because their divergence from perfect spheres is an intrinsic characteristic of the shape. The diameter of the representative sphere can be taken as either the outer or inner diameter of the polyhedron. Both assumptions are considered in this study. The outer diameter of the polyhedron is defined as the diameter of the smallest sphere which can be drawn that contains the whole polyhedron. Alternatively, the inner diameter of the polyhedron is defined as the diameter of the biggest sphere which can be drawn that will be fully contained inside the polyhedron. The outer diameter of the 32-face polyhedron used in this study is 8.2 mm and the inner diameter is 7.5 mm. Similar to the tube-like particles, the hollowness of the polyhedral particles is neglected and the spheres are assumed to be full.

Results
The results of the DEM study are presented in this section along with the corresponding experimental results to validate the numerical approach. Also for comparison, the results of the digital approach applied by Gan et al. (2004) and Jia et al. (2007) based on the method of virtual space are presented. In their method called DigiPac, a particle is just a coherent collection of digitized pixels which are stored as integers or bits. Although DigiPac is able to model particles of arbitrary shapes, consideration of particle interactions is limited to their geometric constraints.
As an initial validation case, the packing behavior of 660 spheres was predicted by DEM. Validation of the simple case of spheres is necessary before advancing to more complex particle geometries. The results, demonstrated in Fig. 1, show that DEM prediction of the bulk density of spherical particles is in agreement with the experimental results.
The results of the prediction of the bulk density of tubular particles are shown in Fig. 2. Despite the slight over-prediction of the height of the packed bed, the DEM prediction of the bulk behavior of the tube-like particles is quite acceptable. The small error can be explained by the approximation of the tube-like shape as a cluster of spheres in the DEM study. Fig. 3 shows the prediction of the bulk behavior of polyhedral particles. It is evident that the inner diameter approximation is a good approximation and provides results in quite good agreement with the experimental results. On the other hand, the outer diameter approximation does not yield a good approximation of the polyhedral shape.

Conclusions
DEM modeling of granular grains of complex shapes is still a young research topic. This study is an attempt to verify the most commonly used method to model non-spherical shapes, i.e. the composite method. The composite method is used to predict the bulk density behavior of granular particles of cylindrical geometry. The numerical results are compared with a digital  method called DigiPac applied by Gan et al. (2004) and Jia et al. (2007). The results of both methods are compared with the experimental results provided by the same authors. The comparison verifies the competence of the composite method in predicting the bulk density of cylindrical particles. The capabilities of the DigiPac method has already been discussed in the author's work and is not repeated here. In addition to cylinders, the bulk density of 32-face polyhedrons is also investigated and it is verified that a 32-face polyhedron can be approximated by a sphere with a diameter being the inner diameter of the polyhedrons. Future work should attempt to predict the bulk behavior of other geometries such as conical, cubic and ellipsoidal particles.

Kasra Samiei
Kasra Samiei graduated from Amirkabir University of Technology, Tehran, Iran as a mechanical engineer specialized in solids design. He later completed his master's degree at the University of Boras, Sweden in "waste management and resource recovery technology". He holds a Ph.D. in engineering science from the University of Luxembourg. His main research focus is on computational granular dynamics in conjunction with computational fluid dynamics. He is currently a post-doctoral researcher at the University of Leoben, Austria.

Girma Berhe
Girma Berhe received his bachelor's degree in computer science (1996), and master's degree in information science (2001)

Bernhard Peters
A graduate in mechanical engineering (Dipl.-Ing.) and with a Ph.D. from the Technical University of Aachen, Bernhard Peters is currently head of the thermo-/fluid dynamics section at the University of Luxembourg and an academic visitor to the Lithuanian Energy Institute. After being a post-doctoral research associate at the Imperial College of London, he established a research team dedicated to the thermal conversion of solid fuels at the Karlsruhe Institute of Technology (KIT), and worked thereafter in industry at AVL List GmbH, Austria. His research activities include thermo-/fluid dynamics, in particular multi-phase flow, reaction engineering, numerical modeling, high-performance computing and motion and conversion of particulate materials.