Dynamic Simulation of the Packing of Ellipsoidal Particles - American

ARTICLE pubs.acs.org/IECR

Dynamic Simulation of the Packing of Ellipsoidal Particles Zong-Yan Zhou,† Rui-Ping Zou,† David Pinson,‡ and Ai-Bing Yu†,* †

Laboratory for Simulation and Modelling of Particulate Systems, School of Materials Science and Engineering, The University of New South Wales, Sydney, NSW 2052, Australia ‡ BlueScope Steel Research, P.O. Box 202, Port Kembla, NSW 2505, Australia ABSTRACT: This paper presents a numerical study of the packing of nonspherical particles by the use of the discrete element method. The shapes considered are oblate and prolate spheroids, with the aspect ratio varying from 0.1 to 7.0. It is shown that the predicted relationship between packing fraction and aspect ratio is consistent with those reported in the literature. Ellipsoids can pack more densely than spheres. The maximum packing fraction occurs at an aspect ratio of 0.6 for oblate spheroids, and 1.80 for prolate spheroids. The packing characteristics with aspect ratio are further analyzed in terms of structural parameters such as coordination number and radial distribution function. It is shown that ellipsoids with small or large aspect ratios tend to give a locally ordered structure. The results demonstrate that DEM provides a useful method to investigate the packing dynamics of ellipsoidal particles.

1. INTRODUCTION Particle packing is important and widely found in nature and industry. The packing of uniform spheres has been widely used to model the structures of liquids, colloids, amorphous, and crystals materials.1 Proper description of particle packing is on the other hand fundamental to many industrial processes ranging from raw material preparation to advanced material manufacturing in many industries.2 It is important to understand and quantify the structure of particle packing in order to generate information for general application. In the past, extensive research has been carried out in this direction.3,4 Early studies are mainly based on laboratory experiments and focused on the packing of uniform spheres, producing some structural properties for general application.36 However, experimental extraction of structural information is laborious and involves relatively large errors and uncertainties. To overcome this difficulty, various numerical algorithms have been developed to simulate particle packing. Most of them are Monte Carlo algorithms and based on either sequential addition or collective rearrangement. As these algorithms usually involve various assumptions about particle motion and/or stability criteria which may not represent the reality truly, the resulting structural information is often not so comparable with that measured.7,8 For example, such an algorithm often can only generate a final and static structure which may satisfactorily represent some but not all structural properties. Forming a packing is a dynamic process which involves various forces, and an ideal simulation method should take into account all dynamic factors related to not only geometry but also force. From this perspective, the use of the so-called discrete element method (DEM) is inherently superior to the early Monte Carlo algorithms. In the past decade or so, DEM has been increasingly used to study the packing of particles under various conditions as summarized by Zhu et al.9 The packing of particles is controlled by many variables among which particle shape is important. Particle shape affects packing structures and hence transport properties such as permeability related to pore connection and thermal conductivity related to r 2011 American Chemical Society

particle connection in a packing. The packing of nonspherical particles has been studied probably since 1930s or even earlier.10,11 One important topic in this area is about the relationship between packing fraction and particle shape. Brown et al.12 produced the first general diagram showing that the deviation from spheres would increase porosity (= 1 packing fraction). They used sphericity, defined as the ratio of surface area between sphere and nonspherical particle with the same volume, to measure particle shape. Their relationship, although largely qualitative, is widely accepted.2 However, Yu and colleagues showed that this relationship is not fully correct.1315 In particular, realizing that for convex nonspherical particles, small sphericity can only be obtained with platy or elongated particles, they mainly used cylinders and disks of different aspect ratios in their experiments under the random close and loose packing conditions as used by Scott.16 Their results indicate as sphericity decreases, porosity decreases to a minimum and then increases, and the relationship between porosity and sphericity varies with packing method and particle shape. To be more quantitative, recent studies are focused on particles of well-defined shapes.1729 In particular, ellipsoids have attracted a lot of attention in recent years. Ellipsoids can represent a wide range of shapes from very platy to elongated. Donev et al.17 used the so-called LubachevskyStillinger (LS) simulation scheme30 and established a relationship between packing fraction and aspect ratio, and demonstrated that the packing fraction of ellipsoids shows a cusplike increase as the aspect ratio deviates from 1 (for spheres), followed by a maximum (∼0.71 at aspect ratio of 0.67 for oblate spheroids and 1.5 for prolate spheroids) and then a strong density decrease at a higher aspect ratio. Similar relationships have also been reported by other investigators using different simulation techniques.28 Received: April 22, 2011 Accepted: July 5, 2011 Revised: July 5, 2011 Published: July 05, 2011 9787

dx.doi.org/10.1021/ie200862n | Ind. Eng. Chem. Res. 2011, 50, 9787–9798

Industrial & Engineering Chemistry Research

ARTICLE

Figure 1. Two-dimensional illustration of forces acting on particle i in contact with particle j.

Efforts have also been made to understand the structure of ellipsoidal particles in terms of coordination number (CN)17,19 and other properties.24,29 However, to date, it is not clear if the packing dynamics of ellipsoidal particles can be properly reproduced by DEM, although such studies become more and more popular for spherical particles.9 In fact, the reported attempts were made in two dimensions for ellipses.22 In this work, a three-dimensional DEM model is developed and validated to study the packing of ellipsoids. It will be shown that the method can indeed simulate the dynamics of a packing process. The relationships between packing fraction and sphericity or aspect ratio are then re-examined, by comparing the previous experimental and simulated ones in the literature. The differences in results are identified and discussed. The packing structure of ellipsoids is further analyzed in terms of CN and radial distribution function (RDF). Such information is useful for structural understanding and further development in this area.

2. SIMULATION METHOD 2.1. Governing Equations. Techniques for DEM modeling of nonspherical particles have been reported in the literature, for example, as reviewed by Dziugys and Peters.31 Various investigators have contributed to the modeling of ellipsoidal particles.3237 For completeness, a brief description of the DEM method for the present work is given below. According to the DEM originally proposed by Cundall and Strack,38 a particle can have two types of motion: translational and rotational. The governing equations for the translational and rotational motion of particle i with mass mi and moment of inertia Ii can be written as

mi

dv i ¼ dt

ki

∑ ðf c, ij þ f d, ij Þ þ mi g j¼1

ð1Þ

and Ii

dω i ¼ dt

ki

∑ ðMt, ij þ Mn, ij þ Mr, ij Þ j¼1

ð2Þ

where vi and ωi are the translational and angular velocities of the particle, respectively, and kc is the number of particles in interaction with the particle. As shown in Figure 1, the forces involved are the gravitational force mig and interparticle forces between particles which include elastic force fc,ij and viscous damping force fd,ij. The torques acting on particle i by particle

j include: Mt,ij generated by the tangential force, Mr,ij commonly known as the rolling friction torque, and also the torque Mn,ij generated by the normal force when the normal force does not pass through the particle center. Equations used to calculate the interaction forces and torques between two spheres have been well-established in the literature.39 The ellipsoidal particles used in our work have smooth/continuous surfaces like a sphere. Hence, the Coulomb condition or sliding/rolling friction models are the same for these spheres and ellipsoids. In this work, we simply extend the nonlinear model to ellipsoids, and those equations are listed in Table 1. This approach was also used by other investigators.4043 It should be noted that an additional torque is introduced, which is caused by the normal force. This is because the normal direction of contact force does not necessarily pass through the center of an ellipsoid, which, together with the tangential forces and rolling torque, governs the rotational motion of the particle. Another parameter is the calculation of the so-called reduced radius R* in the calculation of the contact forces between particles i and j. For spheres, R* is the particle radius. But for ellipsoidal particles, R* = 1/2(A0 B0 )1/2 where A0 and B0 are related to the radii of the particle shape curvature at a contact point.31 2.2. Solution Technique. The explicit time integration method is widely used to solve the translational and rotational motions of a system of discrete particles in the DEM simulations.38 Although it is established for spheres, such a method can be extended to ellipsoids. The difficulties in association with the extension mainly lie in two aspects: particleparticle detection and particle orientation. The methods to overcome the difficulties are described as follows. The detection of particle contacts for nonspherical particles is much more complicated than spherical particles. Various analytical methods have been proposed to detect the contacts between ellipses or ellipsoids, including the so-called intersection algorithm,32 geometric potential algorithm,33,35 and common normal algorithm.35 As noted by Dziugys and Peters,31 the geometric potential algorithm is more reliable. Therefore, it is used in the present work. It should be pointed out that the computational time for contact detection is huge. This is largely because the algorithm used to determine one contact point involves the numerical solution of a sixth-order polynomial equation.35 For a system with a large number of particles, millions of contact points are required to determine at each time step, which significantly extends the simulation time. Particle orientation is another parameter that must be considered for nonspherical particles. It is generally described by three Euler angles (ϕ,θ,ψ).44 Briefly, at each time step, for the convenience to determine the inertia tensor Ii of an ellipsoid, the rotational equation expressed by eq 2 in the space-fixed coordinate system (x,y,z) is converted to the body-fixed coordinate system (x0 ,y0 ,z0 ) which is a moving Cartesian coordinate system fixed with the particle and whose axes are superposed by the principal axes of inertia. Thus, in this converted coordinate system, the angular velocities ω0 i of particles can be calculated as used for spheres, they are then used to determine the new three Euler angles on the basis of the so-called quaternion method. More details about the method can be found elsewhere.31,44 2.3. Simulation Conditions. A packing can be formed under different conditions. One of them is the so-called poured packing.2,45 This rain-dropping procedure has been used by various investigators in their numerical simulations of the packing 9788

dx.doi.org/10.1021/ie200862n |Ind. Eng. Chem. Res. 2011, 50, 9787–9798

Industrial & Engineering Chemistry Research Table 1. Components of Forces and Torques Acting on Particle ia forces

equations

normal

/3E*(R*)

elastic force (fcn,ij) normal

cn(8mijE*(R*δn)1/2)1/2vn,ij

4

1/2

δn3/2n

damping force (fdn,ij) tangential

μs|fcn,ij|(1 (1 δt/δt,max)3/2)δ^t

elastic force (fct,ij) tangential

ct(6 μsmij|fcn,ij|(1 δt/δt,max)1/2/δt,max)1/2vt,ij

damping force (fdt,ij) coulumb

μs|fcn,ij|δ^t

friction force (ft,ij) torque by

Rc,ij (fct,ij + fdt,ij)

tangential forces (Mt,ij) torque by

Rc,ij (fcn,ij + fdn,ij)

normal force (Mr,ij) rolling friction

μr,ij|fn,ij|ω^nij

torque (Mr,ij)

Note: 1/mij = 1/mi + 1/mj, E* = E/2(1 v2), ω^nij = ωnij/|ωnij|, δ^t = δt/ |δt|, δt,max = μs(2 v)/2(1 v) 3 δn, vij = vj vi + ωj Rc,ij ωi Rc,ij, vn,ij = (vij 3 n) 3 n, vt,ij = (vij n) n. Note that tangential forces (fct,ij + fdt,ij) should be replaced by ft, ij when δt g δt,max. a

of uniform spheres.46 It is used in the present work, in order to produce results comparable to those measured under similar conditions. Thus, a packing is formed by pouring particles at a certain height into a cylindrical container, as shown in Figure 2a for example. A small velocity with random direction and random orientation are assigned to each ellipsoid to be dropped. During the packing process, particles may collide with neighboring particles and bounce back and forth. This dynamic process proceeds until all particles reach their stable positions with an essentially zero velocity as a result of the damping effect for energy dissipation, as illustrated in Figure 2b. Figure 2c further shows the variation of the energy of particles with time during the packing process, where the kinetic energy of particles is defined as ∑(mivi2/2 + Iiωi2/2) and the potential energy as ∑(migHi) where H is the height of mass center of particles from the container bottom. It can be seen that the kinetic energy approaches to zero and the total energy reaches a constant after around 1.1 s in the simulation, indicating that the packing structure has been stable, or in a state of mechanical equilibrium. Spheroids with aspect ratios varying from 0.1 to 7.0 are used in this work to produce shapes from oblate to prolate spheroids. Figure 3 shows the packings for three different aspect ratios of spheroids, and their force networks where each stick represents a normal contact force between particles. Detailed analysis of these forces will be studied in our future work. Here we are mainly concered with the common packing properties. The definitions of particle size, aspect ratio, and sphericity are shown in Table 2. In the present work, b = c for spheroids. dv is the equivalent diameter of a sphere with the same volume as an ellipsoid, and set to 10 mm in this work. All the particles with different aspect ratios have the same volume. The table also shows other physical properties of particles and parameters in the DEM simulation. These properties may affect the packing behavior of ellipsoids, as the case for spheres.4650 This will be demonstrated in our comparison with the experimental measurements. However, a systematic

ARTICLE

study is beyond the scope of this work which is focused on the effect of aspect ratio. The exact simulation conditions for the generation of a packing are as follows. First, all particles are dropped from the height of 400 mm layer-by-layer, and particles in each layer are uniformly distributed across the section. Then an initial velocity of 1.25 m/s is added to each particle, randomly distributed among the three directions. The three Euler angles (ϕ,θ,ψ), randomly generated in the range of 090°, are also assigned to each ellipsoid as its initial angles. During the simulation, the variation of the resulting vector OA shown in Figure 4 will be traced all the time, and then the particle orientation can be known.

3. RESULTS AND DISCUSSION 3.1. Macroscopic Packing Properties. To validate the DEM model, comparisons are made between the measured and simulated packing results. The particles used are spherical glass beads (diameter = 8 mm; density = 2500 kg/m3) and regular M&M’s chocolate candies (2a = 6.93 mm, 2b = 2c = 13.4 mm; aspect ratio = 0.517; density = 1349 kg/m3). Following the work of Zou and Yu,15 the measurements are made under the loose and dense random packing conditions. The latter involves vertical tapping. To compare this packing process, one-dimensional vibration in the vertical direction is also implemented in the present work, as used by An et al.51 To be consistent, in the DEM simulation, the bottom of the container is vibrated according to Z(t) = 0.2dv sin[200(t t0)] where Z(t) is the vertical displacement at time t, and t0 is the starting time of vibration. The vibration amplitude and frequency are respectively set to 0.2dv and 200, as used by An et al.51 to generate the random close packing for spheres. After 20 s, the vibration stops, and the bed is allowed to settle for a few seconds to reach its static state. Containers of different diameter D are used. This way, the effect of side wall can be eliminated by extrapolation.16,52 Figure 5 shows the comparison between the measured and simulated results. The agreement is satisfactory considering that variables related to the packing method and material properties are not adjusted in the simulation. It can be seen that with the increase of container diameter, the packing fraction increases as a result of the decreased effect of wall. This can be described reasonably by a linear relationship between packing fractions and dv/D ratio. It should be noted that with the further decrease of container diameter, such a linear relationship will disappear. The presence of a wall can lead to different packing structures. This topic has been extensively discussed.52,53 For example, for uniform spheres, three regions have been identified according to the mechanisms of the wall effect or the roles played by the ordered and disordered packings, that is, dv/D < 0.25, 0.25 < dv/D < 0.53, and dv/D > 0.53. It is not clear if such an approach can be applied to ellipsoids, which should probably be investigated in the future. Figure 5 also shows that packing fraction varies with packing method. Tapping or vibration can transform a poured packing to a denser packing. This is the case for both spheres and ellipsoids. The results also confirm that M&M ellipsoids can give a packing fraction higher than spheres.17 The packing fractions after eliminating the wall effect are consistent with those reported in the literature. For example, the packing fraction is around 0.627 under poured packing condition and 0.651 under tapped or vibrated conditions for spheres.2,51,52 For M&M particles, it is around 0.682 under the current tapped or vibrated conditions, which is very close to 0.685 measured by Donev et al.17 or 0.69 by 9789



ARTICLE

Figure 3. Final packings of ellipsoids with aspect ratios of (a) 0.2; (b) 2.0; and (c) 4.0; and their corresponding force structures (d, e, and f, respectively). Note that the force network only shows a slice of packing structure at |y| e 20 mm, and a thicker/darker stick represents a larger normal contact force between two particles.

Table 2. Particle Parameters Used in the DEM Simulationa parameters container diameter, D dropping height

values 300 mm 400 mm

particle number

25 000

particle size (2a, 2b = 2c, dv)

dv = 2(abc)1/3 = 10 mm

particle aspect ratio (η = a/b)

0.17.0

particle density

1439 kg/m3

sliding friction coefficient, μs

0.3

rolling friction coefficient, μr

0.1 mm

normal damping coefficient, cn tangential damping coefficient, ct

0.3 0.3

Young’s modulus, E

1 107 kg/(m 3 s2)

Poisson ratio, ν

0.3

time step, Δt

3.2 105 s

a

Note: (i) a is the principal radius in the polar direction, b and c are principal radii in the equatorial plane; (ii) for a prolate spheroid, aspect ratio η > 1; for a sphere, η = 1; and for an oblate spheroid, η < 1; (iii) wall is assumed to have the same physical properties as particles.

Figure 2. Formation of a packing of oblate spheroids (aspect ratio η = 0.6) with different microscopic parameters (for visualization of inside the packing, only halves are shown at y g 0): (a) coordination number; (b) particle vertical velocity; and (c) variation of kinetic and total energy of particles with time.

Man et al.18 Plastic beads with prolate shape (2a = 4.86 mm, 2b = 2c = 3.07 mm, aspect ratio = 1.58, density = 910 kg/m3) are also used to make a comparison, and the measured and the simulated packing fractions are, respectively, at 0.674 and 0.672, when dv/D ratio = 0.02. Clearly, the proposed DEM can produce packing results comparable to those measured.

As mentioned earlier, material properties will affect the packing results, as demonstrated in the study of sphere packing.4650 Their effects can also be observed in the packing of ellipsoids. As an example, Figure 6 shows the results for regular M&M ellipsoids, with reference to three material properties: Young’s modulus, sliding, and rolling friction coefficients. It can be observed that in general, increasing these properties can decrease packing fraction, qualitatively consistent with that for spheres.4650 Note that the effect of Young’s modulus becomes insignificant when it is larger than 10 MPa. Compared to particle shape, rolling friction coefficient just gives secondary, if not negligible effect. However, the effect of sliding friction coefficient is always there and comparable to that for spheres. For ellipsoids of different aspect 9790



ARTICLE

Figure 4. Illustration of (a) a prolate spheroid with aspect ratio η = 2; and (b) an oblate spheroid with aspect ratio η = 0.5. The vector OA is used to represent the particle’s orientation (O, particle mass center; A, polar apex).

Figure 6. Effects of particle properties on packing fraction for the regular M&M ellipsoids: (a) Young’s modulus and (b) sliding friction and rolling friction coefficients. Figure 5. Comparison between the simulated and measured packing fractions (solid symbols are experimental data and empty symbols are DEM data); in the inset table, A is the slope of a fitted line, and F0 is the extrapolated packing fraction for simulated data.

ratios, the effects of three material properties may be different. Their quantification needs a systematic study which is beyond the current effort. Nonetheless, the results in Figure 6 confirm that material properties represent one factor responsible for the difference between the measured and simulated results observed in Figure 5. In the following, the DEM results will be discussed against various macroscopic findings in the literature and, at the same time, to develop microscopic understanding. To eliminate the wall effect, the data used for discussion will be only for the inner packing samples at least 23dv away from external surfaces. Unless otherwise specified, all packings are obtained under poured packing conditions. The effect of particle shape on packing fraction has been studied in the past. In particular, Zou and Yu15 established a relationship between packing fraction (1-porosity) and particle sphericity for disks/cylinder-type of particles. Figure 7 shows their results, together with the present DEM ones. Clearly, both sets of results show that packing fraction increases to a maximum and then decreases with the decrease of sphericity. Elongated and platy particles follow different trends. Quantitatively, the two sets

Figure 7. Relationships between packing fraction and sphericity under the poured packing used in the DEM and dense packing used in the experiment.15

of results differ in various ways. First, the two curves for ellipsoids have an intersection at sphericity = 0.83 (corresponding to aspect ratio of 0.36 for oblate spheroids and 3.3 for prolate spheroids). 9791



Figure 8. Variation of packing fraction with aspect ratio under poured conditions.

For sphericity larger 0.83, prolate spheroids have a higher packing fraction than oblate spheroids. For sphericity smaller than 0.83, oblate spheroids give a higher packing fraction than prolate spheroids. Such an intersection is not reported in the work of Zou and Yu,15 where the packing fraction of disks is always higher than that of cylinders. Second, the highest packing fraction occurs at different sphericities. The maximum packing fraction for ellipsoids is at sphericity = 0.95 (corresponding to aspect ratio of 0.6 for oblate spheroids and 1.8 for prolate spheroids). But for cylinders and disks, the maximum packing fraction occurs at sphericity = 0.84. Consequently, the relationship between packing fraction and sphericity differs significantly. The differences are understandable because the packing method and particle shape used are different. However, they indicate that although elongated and platy particles are the right choices for creating particles of sphericity in a wide range, different shaped particles may give different relationships between packing fraction and sphericity. Cylinders/disks must have sphericity less than 0.874, where the height of cylinders/disks is equal to their diameter. They are not as good as ellipsoids with spheres as a special case. The data used in the analysis of Zou and Yu are mainly for sphericities less than 0.88. There is probably a need to experimentally investigate the packing in the vicinity of sphericity = 0.95 where packing fraction is very sensitive to sphericity, as shown for ellipsoids. An alternative way to examine the macroscopic results is the plot of packing fraction against aspect ratio. This can be done for elongated and platy particles of different cross-sectional shape. For ellipsoids, such a relationship has been established by different researchers.1719,22 It would be of interest to see if there is any difference for results generated by different simulation methods. Here we focused on that obtained by Donev et al.17 by means of the LS scheme. Figure 8 shows the DEM results in comparison with those of Donev et al. It can be observed that they are qualitatively comparable. The relationship between packing fraction and aspect ratio gives an M-shaped curve. Ellipsoidal particles can pack more densely than spheres in a certain range of aspect ratios, for example, from ∼0.25 to ∼4.0. If particles are too flat or too elongated, the packing becomes looser than spheres. Quantitatively, there are a few significant differences between the present and previous simulated results. First, the present

ARTICLE

DEM simulation gives packing fractions lower than those by Donev et al.,17 particularly in the vicinity of aspect ratio = 0.6 and 1.8. Also, the aspect ratios at the maximum packing fraction are slightly different. The values reported by Donev et al.17 are 0.67 for oblate spheroids and 1.5 for prolate spheroids. In the present DEM simulation, the two aspect ratios are 0.6 and 1.8, respectively. The accuracy of a simulation technique can be examined against measurements. It appears that the simulation of Donev et al.17 consistently produces higher packing fractions. For example, for regular M&M particles (aspect ratio = 0.517), the packing fraction is 0.68517 or 0.69.18 The simulated result by Donev et al.17 is 0.70. The present DEM gives 0.669 for poured packing and 0.682 for vibrated packing, as shown in Figure 5. Different simulation methods are actually linked to different packing mechanisms. In DEM, particles are permitted to suffer minute deformations, and these deformations are used to calculate the contact forces and torques between colliding particles. The force models employed are based on the classical contact mechanics. The motion of a particle under these forces/ torques is described according to Newton’s laws of motion. Thus, DEM simulation can provide dynamic information, such as the trajectories of and transient forces acting on individual particles. On the other hand, the LS approach, here as an example for comparison because of its wide application, is a hard-particle molecular dynamics algorithm mainly for producing dense packings. In this approach, as described by Donev et al.,17 initially small ellipsoids are randomly distributed and oriented in a box with periodic boundary conditions and without any overlap. The ellipsoids are then given velocities and their motion followed as they collide elastically and also expand uniformly. After some time a jammed state with a diverging collision rate is reached and the density reaches a maximal value. An event-driven MD algorithm is used to implement this process efficiently. This approach can be extended to soft particles, as done by Schreck and O0 Hern.29 Clearly, although both DEM and LS are dynamic, they are quite different, as reflected by the following aspects: (i). Contact Forces between Particles. The DEM method considers the physical contacts between particles, with all the known forces and torques, including the gravity, taken into account explicitly. LS only considers the elastic collisions, assuming that particles are frictionless and the driving force in forming a packing is not necessarily the gravitational force. (ii). Formation of a Packing. DEM treats a packing as a dynamic process and aims to reproduce a real, physical process. Like many previous simulation methods, LS introduces some treatments in its event-driven algorithm to generate a stable packing, and (some of) these treatments may be devised based on assumptions. (iii). Energy Dissipation. DEM dissipates the energy by damping the energy arising from the contact interaction between particles, giving the collisional energy dissipation, frictional energy dissipation, and rolling energy dissipation. LS does not explicitly consider these contact mechanisms but assumes certain rules for particles to rearrange themselves to achieve their jammed state which is considered to correspond to their minimum energy and mechanical stability. Because of these differences, it is not surprising that DEM and LS may generate different results. However, it should be noted that under some conditions, the two methods can still provide consistent results. For the packing of ellipsoidal particles, they can produce qualitatively consistent results as shown in Figure 8. As recently pointed out by Jia et al.,54 this is because for such 9792


Industrial & Engineering Chemistry Research particles, geometrical constraints play their dominant role. Therefore, it is expected that the difference between DEM and LS schemes can be observed for packings formed under strong interparticle forces. Such packing systems are widely found for granular materials, for example, the mechanical forces induced by vibration for the transition from a loose to dense packing,51,55 and the van der Waals or capillary force associated with the packing of fine8 or wet56 particles. DEM has demonstrated its effectiveness for such systems, but it is yet not clear how to take into account the effect of these forces in the algorithms such as LS or geometrical packing. On the other hand, it should be noted that DEM and LS are developed for different research purposes. DEM is for studying granular materials and can handle complicated particle systems as demonstrated by Zhu et al.9 In this connection, probably particle packing is one of the simplest applications. To date, LS is devised for studying packing systems which may not be necessarily relevant to granular materials, although granular spheres are often the starting point of such studies. For example, it has been extended to multidimensional space filling.57 It is probably for this reason that LS and some simulation schemes for packing studies do not consider the frictional and/or rotational motion of particles. As shown in Figure 6, material properties such as Young’s modulus, sliding, and rolling friction coefficients do affect packing behavior. As noted by Liu et al.,7 for granular particles, the previous algorithms, either sequential addition or collective rearrangement, experience a difficulty to produce structural results fully comparable with those measured, even for uniform spheres. Although the LS algorithm is a dynamic simulation, it has some features similar to a collective method without considering the forces and the resulting dynamics, and hence may have the same problem when applied to ellipsoids. As demonstrated in the recent study of sphere packing,7,8,4651,55,56,5860 such difficulty can be readily overcome by DEM-based dynamic simulation where interparticle forces of various types, including dynamic variables related to these forces and packing conditions, can all be explicitly taken into account. In fact, the technique has been widely used to study the dynamics of fluidized ellipsoids.61 In the following, some structural properties of ellipsoids will be examined on the basis of the present DEM simulated results. 3.2. Microscopic Packing Properties. CN, defined as the number of particles in contact with a considered particle, is one of the commonly used parameters describing the structure of a packing. It has been used for many years to distinguish ordered and disordered packings for spheres. For example, for sphere, CN has a single value for an ordered packing but shows a distribution for a disordered packing.62 The typical contacts between spheres in a packing change from 1D chains (CN = 2 when F < 0.3) to 2D triangles (CN = 3 when F ∈ (0.3, 0.40)), and then to 3D tetrahedra (CN = 56 when F ∈ (0.4, 0.66)), and finally to mixed ordered-disordered packing (CN = 612 when F ∈ (0.66, 0.74)).60 It is therefore of interest to examine the CN for ellipsoids. It should be noted that CN is sensitive to the critical separation distance less than which two particles are in contact. In an experiment, a cutoff value has to be assumed due to relatively large measurement errors.5,16,62 Numerical experiments can be more precise. Nonetheless, a critical distance is needed to reflect numerical errors. This distance is normally set to 0.01d (d is particle diameter) in the study of sphere packing.8,59 In this work, unless otherwise specified, this value is also adopted, so that two

ARTICLE

Figure 9. Effect of aspect ratio on the mean CN for different critical distances (real contacts are counted when the distance is 0.0dv).

Figure 10. (a) Effect of critical separation distance on the mean CN under poured and vibrated packing conditions and (b) frequency density distribution for different critical separation distances for the vibrated packing.

particles are considered in contact if their closest surface distance is less than 0.01dv. Different critical distance may give different CN results but as long as it is within a small range, say, less than 9793



ARTICLE

Figure 11. Effects of Young’s modulus and sliding friction coefficient on the overall mean CN for regular M&M ellipsoids.

0.005d, they are qualitatively similar;8 that is, the selection of the critical distance would not affect the dependence of CN on aspect ratio, the major concern in the present study. To test this consideration, Figure 9 shows the variation of the mean CN with aspect ratio for different critical distances varying in a large range from 0.0 to 0.05dv. It can be seen that the critical distance within this range indeed does not affect the trend, but the magnitude. Similar to packing fraction, spheres have the least CN for cases considered in this work. The maximum CN for oblate spheroids occurs when the aspect ratio = 0.45. However, this is not the case for prolate spheroids. As aspect ratio increases, CN increases to the maximum with the corresponding aspect ratio at around 3.5, and is then decreasing very slowly. The trends are similar to those reported by means of the LS algorithm.17,19 However, quantitatively the DEM and LS algorithms produce different CN (the comparison should be made when the critical distance is zero). In general, the CN from DEM simulations is much lower, particularly for small critical distances. Donev et al.17 described a measured CN distribution in a dense packing for regular M&M particles. DEM simulations have been conducted under poured and vibrated packing conditions for such particles as described earlier. Figure 10a shows the variation of mean CN with different critical distance under poured and dense packing conditions. It can be seen that if the critical separation distance is zero which represents the true contacts between particles, the mean CN is 8.19 for dense packing and 8.02 for poured packing. With the increase of the critical distance, CN increases significantly. The reported mean CN in the work of Donev et al.17 is 9.82, which corresponds to the critical separation distance of 0.35 mm in the present work. Figure 10b shows the frequency distribution of CN. It can be seen that the CN frequency distribution depends on the critical separation distance. The agreement between the simulated and measured frequency distributions is satisfactory when the critical separation distance is 0.350.4 mm. As mentioned above, material properties such as Young’s modulus, sliding, and rolling friction coefficients affect packing fraction. They will also affect structural properties such as CN, as shown in Figure 11 (note that the effect of rolling friction is insignificant and hence not shown in the figure). However, for simplicity, the comparison in

Figure 12. Frequency distribution of CN for different aspect ratios when the critical separation distance is 0.01dv: (a) oblate spheroids and (b) prolate spheroids.

Figure 13. Relationship between mean CN (determined when the critical separation distance is 0.01dv) and packing fraction.

Figure 10 does not take into account these effects, because this work is focused on the effect of aspect ratio. 9794



ARTICLE

Figure 14. Radial distribution function for oblate spheroids with different aspect ratios.

Figure 12 shows the CN frequency density distributions for different aspect ratios when the critical separation distance is 0.01dv. It can be seen that as aspect ratio decreases or increases from 1.0, the distribution first moves to the right and then to the left, corresponding to the increased and then decreased mean CN as shown in Figure 9. A decrease in the peak (most probable) value indicates a wider CN frequency distribution. It is of interest to note that for ellipsoids, a high packing fraction does not necessarily correspond to a high mean CN. As shown in Figure 13, for a given packing fraction, the packing structures may be very different. For example, a packing fraction at 0.66 corresponds to four different CNs which correspond to four structures: two of them for oblate spheroids when the aspect ratio is 0.4 or 0.8; and two of them for prolate spheroids when the aspect ratio is 1.2 or 3.0. Such a feature is quite different from that for spheres which show a one-to-one correlation between packing fraction and mean CN.8 RDF is another parameter commonly used to describe the structure of a packing. It is defined as the probability of finding one particle center at a given distance from the center of a given particle, given by g(r) = n(r)/(4πr2ΔrF), where n(r) is the number of particle centers situated at a distance between r and r + Δr from the center of a given particle. F is the number of particles per unit volume in the packing. For uniform spheres, RDF has been well established. Here we try to do something similar for ellipsoids. In the calculation, Δr is set to 0.02dv, and the so-called radial distance is made with reference to the minimum of a, b, or c. Figure 14 shows the results for different aspect ratios. The RDF for uniform spheres has been well established. For the random close packing, there is a split second peak in the RDF, with its first component at 31/2d and the second component at 2d.6 This feature has been reproduced in the previous DEM simulations.7,8,46 As expected, it is also reproduced in the present work as shown in Figure 14a. The peaks in the RDF correspond

Figure 15. Ordered packing structures for (a) spheres; (b) oblate spheroids when aspect ratio = 0.4; and (c) prolate spheroids when aspect ratio = 3.0.

to certain structures. As shown in Figure 15a, the first peak corresponds to particles in contact, and the two components in the second peak respectively correspond to the so-called edgesharing in-plane equilateral triangles and three spheres along a line.63 Such RDF features will change gradually when particles become more ellipsoidal as discussed below. Figures 14 panels b, c, and d show the RDFs for oblate spheroids. It can be observed that with the decrease of aspect ratio, the split second peak observed for spheres vanishes, and the first peak decreases significantly and shifts away from 1.0. Strictly speaking, the first peak is not a peak and featured with a span which increases with aspect ratio. For example, as shown in Figure 14c when the aspect ratio = 0.6, the RDF increases first at the distance of 1.0 (corresponding to the equatorial diameter 2a), then reaches the maximum at d1 = 1.33 and keeps almost a constant, and finally has a sharp decrease at d2 = 1.67. A similar feature can also be seen for the case when the aspect ratio = 0.4 shown in Figure 14d. The RDFs indicate that the packing structures of oblate spheroids are more complicated than spheres. Figure 15b gives the possible structures of oblate spheroids corresponding to the peaks observed in the RDF for the case of aspect ratio = 0.4. Figure 16 shows the RDF for prolate spheroids, which is quite different from oblate spheroids but to some degree similar to spheres. The first peak always occurs at the distance of 1.0, indicating particles prefer to contact at the equator plane as shown by d1 in Figure 15c. The second peak always occurs at the distance of 2.0, reflecting the structures of three particles contacting in a column as shown by d2 in Figure 15c. Some of the key 9795



ARTICLE

’ AUTHOR INFORMATION Corresponding Author

*Tel: 61 2 9385 4429. Fax: 61 2 9385 5956. E-mail address: a.yu@ unsw.edu.au.

’ ACKNOWLEDGMENT The authors are grateful to Australian Research Council (ARC) and BlueScope Steel Research for the financial support of this work, and the NCI National Facility for the support in computation.

Figure 16. Radial distribution function for prolate spheroids with different aspect ratios.

features in sphere packing can also be observed in prolate spheroids packing. The relatively slow variation with aspect ratio may simply reflect the fact that for prolate spheroids, the change of aspect ratio gives a slow deviation for spheres as compared to oblate spheroids. The RDFs in Figures 14 and 16 suggest that ellipsoids have local ordered structures to some degree, which should be investigated further in the future.

4. CONCLUSIONS A DEM model has been developed and validated to study the packing of ellipsoidal particles. On this basis, the effect of aspect ratio has been examined at both macroscopic and microscopic levels. The results from the present study are summarized as follows: (1) The relationship between packing fraction and sphericity varies with packing method and particle shape. Ellipsoids give quantitatively different relationships from those of cylinders/disks. The relationship between packing fraction and aspect ratio shows an M-shaped curve as observed by others. However, the maximum packing fraction occurs at an aspect ratio of 0.6 for oblate spheroids, and 1.80 for prolate spheroids. (2) The variation of CN with aspect ratio is different for oblate and prolate spheroids. For oblate spheroids, CN has a similar feature to packing fraction. But for prolate spheroids, CN increases first, and then tends to be saturated when the aspect ratio is larger than about 3.0. For ellipsoids, a high packing fraction does not necessarily have a high CN. (3) RDF distribution features of uniform spheres vanish gradually with aspect ratio away from 1.0. The first peak in the distribution for oblate spheroids has a span which increases with decreasing aspect ratio. For prolate spheroids, RDF is similar to the spheres to some degree. Ellipsoids tend to have local ordered structures.

’ NOMENCLATURE a, b, c = principal semidiameters of an ellipsoid, m A0 , B0 = defined in eqs 3 and 4, m cn = normal damping coefficient, dimensionless ct = tangential damping coefficient, dimensionless D = container diameter, m dij = distance between the centers of particle i and j, m dp = particle diameter, m dv = equivalent volume diameter, m E = Young’s modulus, Pa fc,ij = particleparticle contact force, N fd,ij = particleparticle damping force, N g = gravitational acceleration, m/s2 H = height of mass center of particles from the container bottom, m I = moment of inertia of particle, kg 3 m2 ki = number of particles in interaction with particle i, dimensionless Mn = torque generated by the normal force, N 3 m Mr = rolling torque, N 3 m Mt = tangential torque, N 3 m m = mass of particles, kg R = vector from the mass center of the particle to the contact point, m R* = reduced particle radius, m t = time, s v = particle translational velocity, m/s x = X direction in a coordinate system, m y = Y direction in a coordinate system, m z = Z direction in a coordinate system, m Greek

δt,max = maximum δt when sliding starts, m δn = relative normal displacement at contact, m δt = relative tangential displacement at contact, m δ^t = unit vector of δt, dimensionless εf = porosity, dimensionless ϕ,θ,ψ = Euler angles, degree η = aspect ratio, dimensionless μr = rolling friction coefficient, m μs = sliding friction coefficient, dimensionless ν = Poisson ratio, dimensionless F = density, kg/m3 ω = angular velocity, 1/s ω^ = unit vector of ω, dimensionless Subscripts

c = contact d = damping i = particle i 9796


Industrial & Engineering Chemistry Research ij = between particles i and j j = particle j n = normal component r = rolling t = tangential component

’ REFERENCES (1) Stillinger, F. H.; Weber, T. A. Packing structures and transitions in liquids and solids. Science 1984, 225, 983–989. (2) German, R. M. Particle Packing Characteristics; Metal Powder Industries Federation: Princeton, NJ, 1989. (3) Bideau, D.; Hansen, A. Disorder and granular media. Elsevier Science Publishers, North-Holland: Amsterdam, 1993. (4) Mehta, A. Granular Matter, an Interdisciplinary Approach; Springer-Verlag: New York, 1994. (5) Scott, G. D.; Knight, K. R.; Bernal, J. D.; Mason, J. Radial distribution of random close packing of equal spheres. Nature 1962, 194, 956–958. (6) Finney, J. L. Random packings and structure of simple liquids. 2. Molecular geometry of simple liquids. Proc. R. Soc. London, Ser. A 1970, 319, 495–507. (7) Liu, L. F.; Zhang, Z. P.; Yu, A. B. Dynamic simulation of the centripetal packing of mono-sized spheres. Phys. A 1999, 268, 433–453. (8) Yang, R. Y.; Zou, R. P.; Yu, A. B. Computer simulation of the packing of fine particles. Phys. Rev. E 2000, 62, 3900–3908. (9) Zhu, H. P.; Zhou, Z. Y.; Yang, R. Y.; Yu, A. B. Discrete particle simulation of particulate systems: A review of major applications and findings. Chem. Eng. Sci. 2008, 63, 5728–5770. (10) Fraser, H. J. Experimental study of the porosity and permeability of clastic sediments. J. Geol. 1935, 43, 910–1010. (11) Brown, R. L.; Hawksley, P. G. W. Packing of regular (spherical) and irregular particles. Nature 1945, 156, 421–422. (12) Brown, G. G. Unit Operations; Wiley: New York, 1950; pp 210. (13) Yu, A. B.; Zou, R. P.; Standish, N. Packing of ternary mixtures of nonspherical particles. J. Am. Ceram. Soc. 1992, 75, 2765–2772. (14) Yu, A. B.; Standish, N. Characterization of nonspherical particles from their packing behavior. Powder Technol. 1993, 74, 205–213. (15) Zou, R. P.; Yu, A. B. Evaluation of the packing characteristics of mono-sized non-spherical particles. Powder Technol. 1996, 88, 71–79. (16) Scott, G. D. Packing of equal spheres. Nature 1960, 188, 908–909. (17) Donev, A.; Cisse, I.; Sachs, D.; Variano, E.; Stillinger, F. H.; Connelly, R.; Torquato, S.; Chaikin, P. M. Improving the density of jammed disordered packings using ellipsoids. Science 2004, 303, 990–993. (18) Man, W. N.; Donev, A.; Stillinger, F. H.; Sullivan, M. T.; Russel, W. B.; Heeger, D.; Inati, S.; Torquato, S.; Chaikin, P. M. Experiments on random packings of ellipsoids. Phys. Rev. Lett. 2005, 94, 198001. (19) Chaikin, P. M.; Donev, A.; Man, W. N.; Stillinger, F. H.; Torquato, S. Some observations on the random packing of hard ellipsoids. Ind. Eng. Chem. Res. 2006, 45, 6960–6965. (20) Wouterse, A.; Williams, S. R.; Philipse, A. P. Effect of particle shape on the density and microstructure of random packings. J. Phys.: Condens. Matter 2007, 19, 376108. (21) Garcia, X.; Akanji, L. T.; Blunt, M. J.; Matthai, S. K.; Latham, J. P. Numerical study of the effects of particle shape and polydispersity on permeability. Phys. Rev. E 2009, 80, 021304. (22) Guises, R.; Xiang, J. S.; Latham, J. P.; Munjiza, A. Granular packing: Numerical simulation and the characterisation of the effect of particle shape. Granul. Matter 2009, 11, 281–292. (23) Williams, R. A.; Jia, X. D. A new method for prediction of bulk particle packing behavior for arbitrary-shaped particles in containers of any shape. Part. Sci. Technol. 2003, 21, 195–205. (24) Zeravcic, Z.; Xu, N.; Liu, A. J.; Nagel, S. R.; van Saarloos, W. Excitations of ellipsoid packings near jamming. Europhys. Lett. 2009, 87, 26001. (25) Sacanna, S.; Rossi, L.; Wouterse, A.; Philipse, A. P. Observation of a shape-dependent density maximum in random packings and glasses of colloidal silica ellipsoids. J. Phys.: Condens. Matter 2007, 19, 376108.

ARTICLE

(26) Wills, J. M. An ellipsoid packing in E3 of unexpected high density. Mathematika 1991, 38, 318–320. (27) Abreu, C. R. A.; Tavares, F. W.; Castier, M. Influence of particle shape on the packing and on the segregation of spherocylinders via Monte Carlo simulations. Powder Technol. 2003, 134, 167–180. (28) Li, S. X.; Zhao, J.; Lu, P.; Xie, Y. Maximum packing densities of basic 3D objects. Chin. Sci. Bull. 2010, 55, 114–119. (29) Schreck, C. F.; O’Hern, C. S. Computational Methods to Study Jammed Systems, in Experimental and Computational Techniques in Soft Condensed Matter Physics; Olafsen, J. S., Ed.; Cambridge University Press: New York, 2010; pp 2561. (30) Lubachevsky, B. D.; Stillinger, F. H.; Pinson, E. N. Disks vs spheres—Contrasting properties of random packings. J. Stat. Phys. 1991, 64, 501–524. (31) Dziugys, A.; Peters, B. An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers. Granul. Matter 2001, 3, 231–265. (32) Rothenburg, L.; Bathurst, R. J. Numerical simulation of idealized granular assemblies with plane elliptical particles. Comput. Geotech. 1991, 11, 315–329. (33) Ting, J. M. A robust algorithm for ellipse-based discrete element modelling of granular material. Comput. Geotech. 1992, 13, 175–186. (34) Ting, J. M.; Khwaja, M.; Meachum, L. R.; Rowell, J. D. An ellipse-based discrete element model for granular-materials. Int. J. Numer. Anal. Met. 1993, 17, 603–623. (35) Lin, X.; Ng, T. T. Contact detection algorithms for threedimensional ellipsoids in discrete element modelling. Int. J. Numer. Anal. Met. 1995, 19, 653–659. (36) Lin, X.; Ng, T. T. A three-dimensional discrete element model using arrays of ellipsoids. Geotechnique 1997, 47, 319–329. (37) Vu-Quoc, L.; Zhang, X.; Walton, O. R. A 3-D discrete-element method for dry granular flows of ellipsoidal particles. Comput. Methods Appl. Sci. 2000, 187, 483–528. (38) Cundall, P. A.; Strack, O. D. L. Discrete numerical-model for granular assemblies. Geotechnique 1979, 29, 47–65. (39) 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. (40) Langston, P. A.; Tuzun, U. Continuous potential discrete particle simulations of stress and velocityfields in hoppers-transition from fluid to granular flow. Chem. Eng. Sci. 1994, 49, 1259–1275. (41) Langston, P. A.; Tuzun, U.; Heyes, D. M. Discrete element simulation of granular flow in 2d and 3d hoppers—Dependence of discharge rate and wall stress on particle interactions. Chem. Eng. Sci. 1995, 50, 967–987. (42) Zhou, Y. C.; Wright, B. D.; Yang, R. Y.; Xu, B. H.; Yu, A. B. Rolling friction in the dynamic simulation of sandpile formation. Phys. A 1999, 269, 536–553. (43) Zhu, H. P.; Yu, A. B. Averaging method of granular materials. Phys. Rev. E 2002, 66, 021302. (44) Goldstein, H. Classical Mechanics; Addison-Wesley Publishing Company: Boston, MA, 1980. (45) Macrae, J. C.; Gray, W. A. Significance of the properties of materials in the packing of real spherical particles. Brit. J. Appl. Phys. 1961, 12, 164–172. (46) Zhang, Z. P.; Liu, L. F.; Yuan, Y. D.; Yu, A. B. A simulation study of the effects of dynamic variables on the packing of spheres. Powder Technol. 2001, 116, 23–32. (47) Yen, K. Z. Y.; Chaki, T. K. A dynamic simulation of particle rearrangement in powder packings with realistic interactions. J. Appl. Phys. 1992, 71, 3164–3173. (48) Silbert, L. E.; Ertas, D.; Grest, G. S.; Halsey, T. C.; Levine, D. Geometry of frictionless and frictional sphere packings. Phys. Rev. E 2002, 65, 031304. (49) Yang, R. Y.; Zou, R. P.; Yu, A. B. Effect of material properties on the packing of fine particles. J. Appl. Phys. 2003, 94, 3025–3034. (50) An, X. Z.; Yang, R. Y.; Zou, R. P.; Yu, A. B. Effect of vibration condition and interparticle frictions on the packing of uniform spheres. Powder Technol. 2008, 188, 102–109. 9797



ARTICLE

(51) An, X. Z.; Yang, R. Y.; Dong, K. J.; Zou, R. P.; Yu, A. B. Micromechanical simulation and analysis of one-dimensional vibratory sphere packing. Phys. Rev. Lett. 2005, 95, 205502. (52) Zou, R. P.; Yu, A. B. The packing of spheres in a cylindrical container—The thickness effect. Chem. Eng. Sci. 1995, 50, 1504–1507. (53) Zou, R. P.; Yu, A. B. Wall effect on the packing of cylindrical particles. Chem. Eng. Sci. 1996, 51, 1177–1180. (54) Jia, X.; Caulkin, R.; Williams, R. A.; Zhou, Z. Y.; Yu, A. B. The role of geometric constraints in random packing of nonspherical particles. Europhys. Lett. 2010, 92, 68005. (55) Yu, A. B.; An, X. Z.; Zou, R. P.; Yang, R. Y.; Kendall, K. Selfassembly of particles for densest packing by mechanical vibration. Phys. Rev. Lett. 2006, 97, 265501. (56) Yang, R. Y.; Zou, R. P.; Yu, A. B. Numerical study of the packing of wet coarse uniform spheres. AIChE J. 2003, 49, 1656–1666. (57) Skoge, M.; Donev, A.; Stillinger, F. H.; Torquato, S. Packing hyperspheres in high-dimensional Euclidean spaces. Phys. Rev. E 2006, 74, 041127. (58) Yu, A. B.; Liu, L. F.; Zhang, Z. P.; Yang, R. Y.; Zou, R. P. Computer simulation of the packing of particles. Int. J. Mater. Prod. Technol. 2003, 19, 324–336. (59) Dong, K. J.; Yang, R. Y.; Zou, R. P.; Yu, A. B. Role of interparticle forces in the formation of random loose packing. Phys. Rev. Lett. 2006, 96, 145505. (60) Dong, K. J.; Yang, R. Y.; Zou, R. P.; An, X. Z.; Yu, A. B. Critical states and phase diagram in the packing of uniform spheres. Europhys. Lett. 2009, 86, 46003. (61) Zhou, Z. Y.; Pinson, D.; Zou, R. P.; Yu, A. B. CFDDEM simulation of gas fluidization of ellipsoidal particles. Proc. 7th International Conference on CFD in the Minerals and Process Industries, Melbourne, Australia, 911 December, 2009, Paper No. 201 Zho (CD Rom ed.). (62) Smith, W. O.; Foote, P. D.; Busang, P. F. Packing of homogeneous spheres. Phys. Rev. 1929, 34, 1271–1274. (63) Clarke, A. S.; Jonsson, H. Structural-changes accompanying densification of random hard-sphere packings. Phys. Rev. E 1993, 47, 3975–3984.

9798


Dynamic Simulation of the Packing of Ellipsoidal Particles - American

Recommend Documents