Rapid Generation of Particle Packs at High Packing Ratios for DEM Simulations of Granular Compacts 1

In this work we present a simple, fast technique for generating particle packs at high packing ratios aiming at the simulation of granular compacts via the discrete element method (DEM). We start from a random sequence addition particle generation algorithm to generate a “layer” of non-overlapping spherical particles that are let to evolve dynamically in time under the action of “compacting” or “jamming” pseudo forces. A “layer-by-layer” approach is then followed to generate multiple layers on top of each other. In the end, very dense packs with pre-defined bulk shapes and sizes (e.g. rectangles in two dimensions and prisms in three dimensions) are achieved. The influence of rolling motion (with particle rotation and spin) along with inter-particle friction on the density and ordering of the generated packs is assessed. Both congruent and inhomogeneous packs (with respect to particle sizes) are created and their packing properties evaluated. We believe that simple techniques for fast generation of particle packs at high packing ratios are essential tools for the DEM simulation of granular compacts.

Latin American Journal of Solids and Structures 13 (2016) 23-50 within a given bulk volume wherein no two particles overlap.The packing ratio or packing density φ is the part of the bulk space that is occupied by the particles relatively to the total bulk space.For two-dimensional bulks, it is known that the highest possible packing ratio for packs of congruent (same-sized) spheres is 0.91 φ ≈ . This is achieved for an ordered, hexagonal arrangement of spheres (here we say "spheres" in a general sense, since in a 2D setting they are actually discs).For three-dimensional bulks, in turn, the densest packing of congruent spheres is the one which exhibits an ordered, either face-centered cubic or hexagonal close-packed lattice (the so-called "fcc" and "hcp" lattices, respectively), for which 18 0.74 φ π = ≈ (Torquato et al., 2000).This has been known for centuries as the Kepler´s conjecture but was proved only very recently (Hales, 2005).For irregular (disordered) packs of same-sized spheres, however, the packing density limit is known to be about 0.64 φ ≈ (Song and Wang, 2008), corresponding to the usually called random close-packed state.For dimensions greater than three, the study of dense packs is a difficult task and has been a matter of intense research for many years now.In the authors´ opinion, the contributions from Kansal and Torquato (2002); Scardicchio, Stillinger and Torquato (2008); Batten, Stillinger and Torquato (2011); Zhang and Torquato (2013), along with Conway and Sloane (1998), appear as the most notorious ones on this regard, especially for setting theoretical upper bounds on the densities and also for bringing out that disordered packs may indeed be the densest packs in higher dimensions.
Many different methods of particle packing have been developed and proposed in the literature.One of the earliest and most notorious contributions is the so-called Lubachevsky-Stillinger generation algorithm (Lubachevsky and Stillinger, 1990), which inaugurated the approach of running a pseudo dynamics simulation of disperse (randomly generated) particles.Various algorithms followed thereafter.The methods of Zinchenko (1994) and Speedy (1998), e.g., were very popular in the 1990s (and indeed are still being used by many researchers nowadays).In Ferrez (2001), the idea of vibrating the particles to increase the density of the arrangement was explored.In 2003, Owen and coworkers proposed an advancing front strategy, wherein an initial assembly of three discs was incrementally increased by adding new discs tangentially to the preexisting ones, until the whole domain was filled.In the same year, Cui and O'Sullivan (2003) proposed a triangulation-based algorithm in which a system of points inside the domain was triangulated (forming a mesh of triangles) and particles generated as the incircles of the triangles (extension to three-dimensions was possible by using a mesh of tetrahedra and inspheres).Numerous variations of these ideas are seen, e.g.Han, Feng, Owen (2005); Stroeven P and Stroeven M (1999); Jiang, Konrad and Leroueil (2003); Siiria and Yliruusi (2007) (the three latter being particularly focused on packs of granular compacts), to cite just a few.Recently, algorithms for more complex shapes and applications have been proposed, like the one by Latham et al. (2008) for arbitrary angular particles (the authors used 3D imaging to capture realistic rock aggregate geometries and create a corresponding mesh of particles) and the one by Thornton, Gong and Chan (2012) for samples aimed at numerical triaxial tests.
In this context, this work presents a simple technique for generating dense hard particle packs of two-dimensional and three-dimensional bulk shapes that can be used in the simulation of granular compacts via the discrete element method (DEM).Our approach starts from a random sequence addition (RSA) particle generation algorithm to generate a layer of non-overlapping spherical particles within the bulk space, which are then let to evolve dynamically in time under the action of "compacting" or "jamming" pseudo forces.A DEM model is used to simulate the pseudo dynamics of Latin American Journal of Solids and Structures 13 (2016) 23-50 this compaction.After the layer has settled down, a "layer-by-layer" approach is followed to generate multiple "jammed" layers on top of each other.In the end, very dense packs with pre-defined bulk shapes and sizes (e.g.rectangles in two dimensions or prisms in three dimensions) are achieved.Following this technique, we construct several sets of packs of a given particle size or size distribution and perform statistical assessment of the attained densities by computing their mean value, standard deviation and coefficient of variation.Reliability of the technique for generating more or less repeatable (yet randomly created) dense packs is thereby evaluated.Also, the influence of interparticle rolling motion and friction during the compaction stage is assessed to verify whether these physics improve or spoil the quality of the results w.r.t.density and ordering.We emphasize that in this work we are concerned primarily with the density of the packs, and thus other pack properties such as coordination number, fabric tensor and contact tensor are not assessed (though they could have been easily computed).
Our technique has connections with the well-known, long-established Lubachevsky-Stillinger algorithm of Lubachevsky and Stillinger (1999), in the sense that a pseudo dynamics simulation of randomly generated particles is performed to squeeze the particles, but differs from it in that (i) we do not employ particle growth (thus congruent packs are made possible in a more natural, easier way), (ii) we consider inter-particle friction with rolling motion in the dynamics (this is possible since we adopt a DEM description instead of a molecular dynamics one), (iii) our pseudo jamming forces are of a different nature and (iv) we follow a layer-wise procedure.In this latter aspect, our technique may be viewed as a variation of the procedure proposed by Jiang and coworkers in 2003.It is important to remark that our goal here is not to devise a technique that generates the densest packs possible for a given bulk shape and particle size, but instead one that creates packs that are dense enough and randomly (dis)ordered enough as to serve as samples of granular compacts that can be utilized for computational simulations.The possibility to generate large numbers of different samples with such characteristics is of primary interest, since it allows us e.g. to run multiple simulations for statistical assessment of mechanical responses.Simplicity and rapidity are thereby crucial matters for this technique.It is worth mentioning also that, though our original motivation is the simulation of granular compacts, the technique proposed herein can be used to study many other different systems and industrial processes such as adsorption of colloidal particles onto surfaces from dense suspensions (thin-film deposition processes), formation of self-assembling biological membranes, deposition of sintering powder for laser sintering and selective laser sintering manufacturing processes, and many others.For an early history of the particle packing subfield, we refer the reader to Torquato, Truskett and Debenedetti (2000); Kansal, Torquato and Stillinger (2002); Speedy (1998); Han, Feng and Owen (2005) and the comprehensive work of Donev (2006), and references therein.For reviews on the discrete element method and on the physics of granular media, see e.g.Cundall and Strack (1979); Bicanic (2004); Zhu et al. (2008); O´Sullivan (2011); Zohdi (2012);and Duran (1997); Pöschel and Schwager (2004) and Thornton (1999), respectively.
The paper is organized as follows.In Section 2 we present the technique including an algorithmic overview.In Section 3 we provide a brief description of the DEM model that is adopted to simulate the particles´ dynamics within the compaction stages, showing in Section 4 its related time integration scheme (also with an algorithmic overview).In Section 5 we use the technique to generate multiple packs of both congruent and inhomogeneous particles, assessing the quality of the results with respect to attained densities and orderings.In Section 6 we derive our conclusions.Throughout the

THE PROPOSED TECHNIQUE
Our technique generates packs in a layer-by-layer fashion beginning from the bottom of the bulk domain until a desired total height is achieved.The bulk domain (shape and sizes) has to be specified a priori, as e.g. a rectangle for 2D packs or a prism for 3D ones.For each new layer, the positions of the particles are treated as random variables whose values are sampled from underlying probability distribution functions (pdfs), following a random sequence addition (RSA) approach.More precisely, the new layer is created by (1) defining a region of thickness (or height) 0 h within the bulk domain on top of the previously generated layers and then (2) filling it with particles in an RSA fashion until a pre-specified initial value of packing density 0 φ is achieved.The layer, at this stage, consists of relatively disperse particles.A jamming force is then activated and the motion acquired by the particles is resolved through a DEM simulation.This jamming causes the new layer to settle down and interact with the previous layers via mechanical contact, leading the whole bulk to a denser state that possess significantly higher packing density φ than the initial value 0 φ of the top layer.The procedure is repeated for new layers until a desired total height d H for the pack is achieved.Figure 1 provides a schematic illustration.In each new layer, the radii of the generated particles can either be the same (a priori specified) for all particles or be sampled from ad-hoc pdfs, e.g. from Gaussian distributions wherein a corresponding mean value and standard deviation must be provided (in such a case, a "truncated" distribution is employed, with given minimum and maximum values for the radii).This latter option is useful for generating packs of granular materials according to their granulometric curves (which, of course, must be known in advance), such as soils, sintering powders, etc.We remark that we do not adopt a particle growth scheme in the jamming stage of our technique, as opposed to the Lubachevsky-Stillinger algorithm and its many variations.This allows us to generate congruent packs in a more natural way, and also inhomogeneous packs that entirely preserve the original granulometric distribution given as input data.
Four aspects must be controlled for the technique to attain density effectiveness and computational efficiency.They are: (1) the use of appropriate material parameters for the particles within the DEM simulations, so that the contacts and collisions due to the jamming forces result inelastic enough (this enforces the system to settle down very rapidly for each new layer, preventing undesired bouncing of particles as well as high frequency motions, which would require small time steps to be resolved adequately); (2) the use of an explicit instead of implicit time integration scheme to perform the DEM simulations (this boosts up computational efficiency); (3) the start of each new layer at moderate initial densities (like 0 0.4 φ = in 2D or 0 0.25 φ = in 3D) rather than at lowdensity states like the ideal gas state as done in other generating techniques (this feature is inherited from the RSA algorithm); and (4) consideration of rotational motion with particle spin and inter-particle friction in the DEM simulations (this, though requiring more computations and memory space in the computer, leads to slightly denser states, as will be shown in Section 5).It is important to remark that one does not need to be concerned with accuracy issues on the DEM simulations, since these are only pseudo-dynamics simulations of the system.Thereby, moderate to relatively large time steps for the explicit integration may be adopted, the only concern on this regard being the avoidance of numerical instabilities (i.e.system blow-up).
The technique can be implemented into a script code to generate multiple packs of particles, or even sets of multiple packs, which can be used for multiple DEM simulations allowing for statistical assessment of system responses.It is summarized in Algorithm 1 below.

SYSTEM PSEUDO DYNAMICS
Following a DEM approach, we treat the generated particles within the bulk domain as a discrete dynamical system in which each particle interacts with the others and the surrounding media via a combination of gravity forces, drag forces, near-field (attractive and repulsive) forces, and contact and friction forces due to touching and collisions.The time evolution of the system is described by the equations of motion from classical dynamics, which are in turn solved via a numerical (timestepping) integration scheme.The particles are allowed to have both translational and rotational motions (in this sense, the DEM model shown in this Section is a generalization of the models presented by Campello and Zohdi (2014) in two articles; Campello (2014), wherein rotations and spins were not considered).).We denote the position vector of a particle by i x , the velocity vector by i v and the Algorithm 1. Layer-by-layer generation of packs of particles.Script for sets of multiple packs.
1. Get bulk domain data (bulk shape and side dimensions) 2. Get layers' properties (initial height 0 h , initial density 0 φ , particle size or size distribution and particles' material properties) spin vector by i ω , as depicted in Fig. 2. The rotation vector relative to the beginning of the motion is denoted by i α , whereas the incremental rotation vector (rotation vector relative to two consecutive configurations at discrete time instants) by i ∆ α (we remark that, for spherical particles, though the particle´s spin is relevant to the particle´s motion, its rotation may be totally irrelevant, even when i ≠ 0 ω ; exception holds e.g. when stick-slip friction with other particles is to be expected, since the rotation of the contact point between the contacting particles needs to be mapped.This type of friction will not be considered in this work).
We denote the total force vector acting on particle i by tot

, ,
where the superposed dot denotes differentiation with respect to time and i j is the particle´s rotational inertia relative to its center of mass, i.e., . The total force vector reads as in which g is the gravity acceleration vector, drag i f is the drag force vector (standing for viscous effects of the surrounding media on the motion of the particle), nf i f are the forces due to near-field interactions with other particles, con i f the forces due to mechanical contacts (or collisions) with other particles and/or obstacles, and fric i f the forces due to friction that arise from these contacts or collisions.The total moment vector, on its turn, has contributions only from the friction forces, since all other forces are assumed to be central forces (i.e. they act with no eccentricity relatively to the center of the particle), such that , where fric i m is shown later in equation ( 15).We adopt standard expressions for the force contributions in equation ( 2).The drag force, for example, is given by ( ) , where fluid c is a damping parameter depending on the viscosity of the surrounding fluid and fluid v is the (local) velocity of the fluid.The forces due to near-field interactions with other particles, on their turn, are given by ( ) where nf ij f is the near-field force between particle i and particle j , in which the α ´s and β ´s are scalar parameters dictating the intensity of the force for the pair { , i j } ( 1 α and 1 β are related to the attractive part whereas 2 α and 2 β to the repulsive part) and ij n is the unit vector that points from the center of particle i to the center of particle j , i.e., which is referred to as the pair´s central direction.The forces due to contact/collisions with other particles and/or rigid walls are described following Hertz´s elastic contact theory (see e.g.Johnson, Kendall and Roberts, 1971), according to which where is the contact force between particle i and particle (or wall) j , C N is the number of particles (and/or walls) that are in contact with particle i , 2 2 (1 ) (1 ) are the effective radius and the effective elasticity modulus of the contacting pair { , i j } (in which ,  ( ) is the geometric overlap (or penetration) between the pair in the pair´s central direction, ij δ ɺ is the rate of this penetration, and Latin American Journal of Solids and Structures 13 (2016) 23-50 is a damping constant that is introduced to allow for energy dissipation in the pair´s central direction.This constant is taken here following the ideas of Wellmann and Wriggers (2012), wherein ξ is the damping rate of the collision (which must be specified) and m * is the effective mass of the contacting pair, i.e., .
Fig. 3 (top part) provides a schematic illustration of the contact/collision for a contacting pair.The forces due to friction are given by assuming that the friction coefficients of all contacting pairs are small enough so that a continuous sliding (with an opposing dynamic friction force and a corresponding moment) is to be expected during the entire duration of a contact/collision (see Fig. 3, bottom part).By continuous sliding, we mean that there is to be no sticking between a contacting pair.Although a stick-slip model could be considered (e.g.following the scheme recently proposed by Campello in 2015), we find it unnecessary for the dynamics we are concerned with in this work (stick-slip models are much more expensive when compared to continuous sliding ones, since they need to keep track of the contact history for every contacting pair and check the Coulomb´s stickslip condition at all iterations; considering that our goal here is merely to generate dense random packs of particles, irrespective of the trajectory that the particles undergo throughout the generation process, a continuous sliding model suffices).Thereby, here we write is the friction force between contacting particles i and j , d µ is the coefficient of dynamic friction for the contacting pair and , , , , is the tangential direction of the contact (or sliding direction), which is the direction of the tangential relative velocity of particles i and j , wherein one has ) (subscript " τ " above stands for tangential direction).The moment generated by the friction forces on particle i (relatively to the center of the particle), on its turn, is given by where n is the vector that connects the center of particle i to the contacting point j C with particle j , as indicated in Fig. 3 (bottom part).
Remark 1.In the context of particle packing, one could argue that the solution of the system dynamics could be speed up by adopting an impulse-based scheme for solving the contacts/collisions, as proposed in Campello and Zohdi (2014), or even by an event-driven scheme.In principle, this would allow the use of much larger time steps within the "WHILE" loop of Algorithm 1 (Section 2).This is, however, only partially true.As the packing becomes more and more compact due to the jamming forces, the collision rates of particles are increased and a state of "enduring contact" for all or nearly all particles is approached.This requires small times steps irrespective of the contact scheme that is being adopted.Indeed, by following an impulse-based or an event-driven scheme, the so-called "stuck-in-time failure" or "inelastic collapse" could eventually show up and prevent the algorithm to progress toward a jammed state (McNamara and Young, 1994).Thereby, in our technique we opt for an overlap-based scheme.This will also be the scheme that we will use in the future for simulating the real physics of the packs.

TIME INTEGRATION SCHEME FOR SOLUTION OF THE SYSTEM´S PSEUDO DYNAMICS
To solve the system´s pseudo dynamics, our scheme starts by performing time integration of equation (1) between time instants t and t t + ∆ , which furnishes The integrals on the right-hand side of ( 16) are then approximated by means of a generalized trapezoidal rule: in which 0 1 θ ≤ ≤ .When 0 θ = , the integration amounts to an (explicit) forward Euler scheme, which is the one we use in Algorithm 1 of Section 2; when , to an (implicit) classical trapezoidal rule.By inserting ( 17) into ( 16) we have On the other hand, by time integration of the velocity and incremental rotation vectors between t and t t + ∆ we have Algorithm 2. Explicit time integration scheme for solution of the system´s pseudo dynamics 1. Initialize solution with known (given) quantities: ii. Loop over particles: FOR 1,..., P i N = DO Update velocity and spin vectors: Update position and incremental rotation vectors: The generalized trapezoidal rule is then invoked again to approximate the integrals on the righthand side of ( 19), rendering By introducing ( 20) into ( 19) we arrive at Expressions ( 18) and ( 21) constitute a set of equations for 1,..., P i N = particles, with which the velocity, spin, position and incremental rotation vectors of each particle at t t + ∆ may be computed once ( ) and ( ) i t x are known.When 0 θ ≠ , this computation cannot be performed directly, since (18) requires the evaluation of ( ) , which in turn are functions of all unknown position, velocity, spin and incremental rotation vectors at t t + ∆ 1 .On the contrary, when 0 θ = , all equations are uncoupled (this can be seen by doing 0 θ = in ( 18) and ( 21)) and the solution turns to be very efficient.Accordingly, the scheme is as summarized in Algorithm 2 above.

EXAMPLES
In this section, we employ our technique to generate multiple packs of particles of both 2D and 3D bulk shapes.A gravity force of magnitude g is adopted as the jamming force, acting in the direction from top to bottom of the bulk domain.A slightly viscous fluid is considered within the bulk space to enforce a small drag on the particles, allowing for some energy dissipation.For each generated pack, we compute the attained density and assess (qualitatively) the ordering of the particles' arrangement.The influence of rotational motion on the results is also investigated.
We remark that for computation of the densities we do 1 This means that, in the implicit versions of the scheme, all equations are strongly coupled and a recursive solution strategy is necessary.Here, since we are interested in an explicit integration (for reasons pointed out in Section 2), we will omit details on such case and refer the interested reader to Campello (2015).
where P V is the sum of the volumes (or areas) of the particles inside the bulk domain and V is the volume (or area) of the bulk domain, computed with the lateral dimensions of the bulk and its desired height d H .Here it is important to mention that, since some particles may end up slightly above d H at the end of the pseudo dynamics, they do not contribute to ( 22), whereas those that end up part above and part below d H contribute only partially.The following data are common for all examples: • Mass density of the particles = 1000 kg/m 3 ; • Elasticity modulus and Poisson coefficient of the particles (needed to resolve contacts/collisions): • Near field forces: • Initial velocity and spin of each newly generated particle: ( • Time step size for the DEM simulations: • Duration of the jamming stage for each newly generated layer (i.e., DEM simulation times): • Typical total number of particles in a pack: We recall that selection of the time step size does not need to be based on accuracy arguments, since the DEM simulations here are only pseudo-dynamics simulations of the system.Thereby, t ∆ may be chosen solely as to avoid numerical instabilities (i.e.system blow-up).We consider the contact force of equation ( 7), which is a highly nonlinear force w.r.t.time, and ensure stability by guaranteeing that it be sufficiently well integrated.In other words, time-steps that are moderately large may be selected, but not too large such that both ij δ and ij δ ɺ (see ( 7)) be satisfactorily computed.Here we use 0.0001 t ∆ = s, chosen so as to ensure at least five time steps per contact/collision.This value is arrived at by using the Hertz´s formula for the duration of elastic collisions in Johnson (1985), and then dividing the result by five, i.e., Latin American Journal of Solids and Structures 13 (2016) 23-50 where rel v is the relative velocity of a typical contacting pair in the pair´s central direction at the beginning of the contact/collision.

2D Rectangular Packs of Congruent Particles
Rectangular packs of length .With these data, ten packs are generated by means of Algorithm 1.The results are shown in Fig. 4, wherein snapshots of the final configuration of each pack are depicted.The final densities attained for the packs are presented in Table 1, middle column.We repeat the procedure and generate another ten packs but now turning off the rotational degrees of freedom (spins and rotations) during the DEM simulations of the compaction stages.The results are shown in Fig. 5, and the corresponding densities are depicted in Table 1, right column.
It is interesting to notice the strong differences observed in the ordering of the particles between the packs of Figs. 4 and 5.When the rotational motion is considered, the packs tend towards a regular, ordered arrangement (apart from wall effects near the boundaries, as is expected) that resembles the hexagonal lattice of the densest congruent pack in two dimensions.Moreover, the attained densities are closer to the maximum possible density for these types of packs (which, as mentioned in Section 1, is ).On the other hand, by inhibiting the rotations, the packs attain disordered arrangements with corresponding smaller densities.We found this to be an interesting result.Our explanation is that, within the dynamics of the compaction stage, the particles will look for states that are more stable and thereby of lower energy, and this is more easily achieved through combinations of translational and rotational motions -these last being triggered by interparticle friction.If rotations are not allowed, the particles will arrange themselves in equilibrium states of slightly higher energy levels.Also, it is interesting to notice the very small value of the standard deviation of the densities in both cases.This indicates that our technique, though based on random generation, is able to generate packs that are satisfactorily repeatable w.r.t.density and order.The small value of the coefficient of variation (the ratio of standard deviation to mean, which provides a non-dimensional measure of dispersion of the obtained φ ´s, or in other words, a zerothorder estimate of the extent to which variabilities in the input parameters propagate to the output results) also supports this conclusion.

2D Rectangular Packs of Inhomogeneous Particles
Here we repeat the input data of the previous example but now generate packs of inhomogeneous particles.Accordingly, the particles´ radii are enforced to follow a Gaussian distribution of mean m, which encompass more than 99% of the distribution.The results obtained with consideration of rotational motion are depicted in Fig. 6, whereas the ones for inhibited rotations are shown in Fig. 7. Table 2 summarizes the attained densities for both cases.
It is noteworthy to observe that the densities are slightly smaller than those of the previous example.This is somewhat counterintuitive.The explanation is that, though the mean value of the radii is here equal to the radii of the previous example, the fluctuations on this geometric parameter spoil the ability of the system to arrange according to the ordered, hexagonal pattern of the maximum possible density state.The packs are all disordered and thereby of smaller densities, even when there is rotational motion (though the presence of rotational motion slightly increases the density).Albeit this aspect, the standard deviation and the coefficient of variation of the attained densities are again very small, indicating the good repeatability of the technique.

3D Brick Packs of Congruent Particles
This is the three-dimensional version of example 5.1, with some slightly modified data.Accordingly, the packs have a rectangular brick shape with horizontal base lengths .The results obtained with consideration of rotational motion are depicted in Fig. 8, whereas the ones for inhibited rotations are shown in Fig. 9. Table 3 summarizes the attained densities for both cases.
One can observe the same general behavior as in the previous examples: the presence of rotational motion helps to increase the density of the packs.One aspect, however, must be realized: the arrangement of the particles within the packs is not completely ordered.It does not fully resemble the fcc (face-centered cubic) or the hpc (hexagonal close-packed) lattices of the densest congruent pack in three dimensions, as it would be expected.This is explained by the perturbations induced by the walls, which have a more pronounced extent here since the size of the particles is not negligible when compared to the side lengths of the packs.(Indeed, in a separate simulation, we have generated one pack with particles of smaller size -same radius as in example 5.1 -and verified that the fcc lattice is obtained in the interior of the pack; for the case of inhibited rotations, regions of fcc lattice are also observed, however not prevalent throughout the whole domain of the pack).As a consequence, the attained densities in Table 3 are far from the maximum possible density for congruent 3D packs (which, as mentioned in Section 1, is 0.74 φ ≈ ) and instead are closer to the packing density limit for irregular packs (which is 0.64 φ ≈ ).All the same, the standard deviation and the coefficient of variation obtained are again very small, indicating the good repeatability of the technique also for three-dimensional bulks.

3D Brick Packs of Inhomogeneous Particles
Here we consider the same data of the previous example, except that the packs are now generated with inhomogeneous particles.The particles´ radii follow a Gaussian distribution with mean  4 summarizes the attained densities for both cases.Once again, the presence of rotational motion is seen to increase the density of the packs.However, contrary to what is expected, the attained values are (slightly) higher than those of the congruent packs (i.e.those of the previous example) -though only slightly.This is explained not by an increase in the ordering of the packs (as it can be seen, the system still faces difficulty in arranging according to the fcc or hpc lattices), but instead by the void-filling effect that builds up due to the different sizes of the particles, which compensates part of the wall perturbations and marginally increases the densities.Still, the standard deviation and the coefficient of variation of the densities are once more very small, showing the good repeatability of the technique.

CONCLUSIONS
A good algorithm for random packing of particles can considerably reduce the cost of creation of compact initial configurations frequently required in discrete element simulations.Motivated by this idea, the main purpose of this work was to present a simple, fast technique for generation of packs of spherical particles at high packing densities with which DEM simulations of granular compacts may be conducted.Our approach relies on random generation of particles in a layer-wise fashion coupled with a compaction scheme to jam the particles to denser states, resolving the compaction with a robust DEM model.It can be implemented with small effort by researchers interested in the field.The technique proved to generate packs at sufficiently high densities for granular compacts, and also to be very reliable w.r.t.repeatability of the packing properties, i.e., it is able to generate packs with densities that are very much repeatable though being based on random procedures.One interesting conclusion is that the consideration of rotational motion during the compaction stages improves the density and ordering of the packs.Indeed, attaining the ordered lattice of the densest possible pack for a given problem dimension seems to be possible only through its consideration.We believe that simple techniques for fast generation of particle packs at high packing ratios are essential tools for the DEM simulation of granular compacts.

Figure 1 :
Figure 1: Schematic illustration of the packing procedure (case of rectangular 2D packs).
E.M.B. Campello and K.R. Cassares / Rapid Generation of Particle Packs at High Packing Ratios for DEM Simulations of Granular Compacts Latin American Journal of Solids and Structures 13

Figure 2 :
Figure 2: Description of a single particle.
the elasticity modulus and the Poisson coefficient of i and j , respectively),
encompasses more than 99% of the distribution).The results obtained are depicted in Fig.10(with consideration of rotational motion) and Fig.11(for inhibited rotations).Table