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

Abstract

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.

Keywords:
Particles; particle packing; granular materials; discrete element method

1 INTRODUCTION

The study of particle (especially hard spheres) packing and the development of particle packing algorithms have been topics of great interest among the computational physics and computational mechanics communities. They are important in the simulation of a myriad of systems such as (but not limited to) colloids (Bolintineanu et al., 2014), crystals (Tan et al., 2009), polymers (Kroupa et al., 2012), powders (David et al., 2007), particulate flows (Kamrin and Koval, 2014) and granular materials (Xu and Chen, 2013). Packs of "hard" particles are defined as arrangements of particles 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 ϕ = π/ ≈ 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 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 inter-particle 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 text, plain italic letters (a, b, ..., α, β, ..., A, B, ...) denote scalar quantities, whereas boldface italic ones (a, b, ..., α, β, ..., A, B, ...) denote vectors in a three-dimensional Euclidean space. The (standard) inner product of two vectors is denoted by u v , and the norm of a vector by ||u|| = .

2 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) h0 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 Hd for the pack is achieved. Figure 1 provides a schematic illustration.

Figure 1
Schematic illustration of the packing procedure (case of rectangular 2D packs).

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 low-density 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.

Algorithm 1
Layer-by-layer generation of packs of particles. Script for sets of multiple packs.

3 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 (time-stepping) 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).

Let the system be comprised of Np particles, each one with mass mi and radius ri (i = 1,...,Np). We denote the position vector of a particle by xi, the velocity vector by vi and the 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 (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).

Figure 2
Description of a single particle.

We denote the total force vector acting on particle i by and the total moment (with respect to the particle´s center) by . From the Euler´s laws, the following equations must hold for each particle at every time instant t:

where the superposed dot denotes differentiation with respect to time and ji 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, is the drag force vector (standing for viscous effects of the surrounding media on the motion of the particle), are the forces due to near-field interactions with other particles, the forces due to mechanical contacts (or collisions) with other particles and/or obstacles, and 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 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 cfluid is a damping parameter depending on the viscosity of the surrounding fluid and vfluid is the (local) velocity of the fluid. The forces due to near-field interactions with other particles, on their turn, are given by

where 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 nij 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, NC is the number of particles (and/or walls) that are in contact with particle i,

are the effective radius and the effective elasticity modulus of the contacting pair {i, j} (in which Ei, Ej and vi, vj are the elasticity modulus and the Poisson coefficient of i and j, respectively),

is the geometric overlap (or penetration) between the pair in the pair´s central direction, is the rate of this penetration, and

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, 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 stick-slip 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

Figure 3
Contact/collision between two particles.

where 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 "r" 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 is the vector that connects the center of particle i to the contacting point Cj 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.

4 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 θ =1, to an (implicit) backward Euler one; and when θ = 05, 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

The generalized trapezoidal rule is then invoked again to approximate the integrals on the right-hand side of (19), rendering

By introducing (20) into (19) we arrive at

Expressions (18) and (21) constitute a set of equations for i = 1, ..., Np particles, with which the velocity, spin, position and incremental rotation vectors of each particle at t + Δt may be computed once vi(t), ωi(t) and xi(t) are known. When θ ≠ 0, this computation cannot be performed directly, since (18) requires the evaluation of and , which in turn are functions of all unknown position, velocity, spin and incremental rotation vectors at t + Δt . 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.

5 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

where Vp 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 Hd. Here it is important to mention that, since some particles may end up slightly above Hd at the end of the pseudo dynamics, they do not contribute to (22), whereas those that end up part above and part below Hd contribute only partially. The following data are common for all examples:

• Mass density of the particles = 1000 kg/m3;

• Elasticity modulus and Poisson coefficient of the particles (needed to resolve contacts/collisions): E = 108 N/m2 and v = 0.25;

• Damping rate for contacts/collisions: ξ = 0.1;

• Coefficient of dynamic friction for contacts/collisions: μd = 0.2;

• Magnitude of gravity force (jamming pseudo-force): g = 2 m/s2;

• Drag force parameters: cfluid = 0.005 N∙s/m and vfluid = o;

• Near field forces: = o (α1 = α2 = β1 = β2 = 0);

• Initial velocity and spin of each newly generated particle: vi(0) = ωi(0) =o

• Time step size for the DEM simulations: Δt = 0.0001 s;

• Duration of the jamming stage for each newly generated layer (i.e., DEM simulation times): t = 1.0 s.

• Typical total number of particles in a pack: Np = 1300 (2D packs) and Np = 7500 ~ 8000 (3D packs).

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 (see (7) be satisfactorily computed. Here we use Δt=0.0001 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.,

where vrel is the relative velocity of a typical contacting pair in the pair´s central direction at the beginning of the contact/collision.

5.1 2D Rectangular Packs of Congruent Particles

Rectangular packs of length L = 1.0 m, desired height Hd = 0.5 m and congruent particles of radius ri = 0.001 m (i = 1, ..., Np) are considered in this example. The initial thickness of each layer is set to h0 = 0.2Hd, whereas the initial packing density is ϕ0 = 0.4. 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.

Figure 4
2D packs of congruent particles (rotational motion is ON).

Figure 5
2D packs of congruent particles (rotational motion is OFF).

Table 1
Packing densities (ϕ) obtained for example 5.1.

Pack	rotations ON	rotations OFF
1	0.854	0.813
2	0.854	0.838
3	0.854	0.841
4	0.863	0.830
5	0.854	0.838
6	0.864	0.837
7	0.853	0.840
8	0.863	0.839
9	0.892	0.846
10	0.862	0.846
mean	0.861	0.837
st. deviation	0.012	0.010
coef. variation	0.014	0.012

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 ϕ ≈ 0.91). 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 inter-particle 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 zeroth-order estimate of the extent to which variabilities in the input parameters propagate to the output results) also supports this conclusion.

5.2 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 rm = 0.001 m and standard deviation σr = 0.1rm m. The truncation limits of the distribution are set to rmin = rm - 3σr = 0.007 m and rmax = rm + 3σr = 0.013 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.

Figure 6
2D packs of inhomogeneous particles (rotational motion is ON).

Figure 7
2D packs of inhomogeneous particles (rotational motion is OFF).

Table 2
Packing densities (ϕ) obtained for example 5.2.

Pack	Rotations ON	Rotations OFF
1	0.826	0.803
2	0.830	0.827
3	0.825	0.814
4	0.819	0.815
5	0.822	0.814
6	0.825	0.819
7	0.826	0.813
8	0.822	0.810
9	0.827	0.806
10	0.826	0.806
mean	0.825	0.813
st. deviation	0.003	0.007
coef. variation	0.004	0.009

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.

5.3 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 B = L = 1.0 m and desired height Hd = 0.5 m. Congruent particles of radius ri = 0.02 m (i = 1,...,Np) are considered to generate the packs. The initial thickness of each layer is h0 =0.2Hd and the initial packing density is set to ϕ0 = 0.25. 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.

Figure 8
3D packs of congruent particles (rotational motion is ON).

Figure 9
3D packs of congruent particles (rotational motion is OFF).

Table 3
Packing densities (ϕ) obtained for example 5.3

Pack	Rotations ON	Rotations OFF
1	0.560	0.540
2	0.576	0.538
3	0.572	0.539
4	0.608	0.542
5	0.564	0.534
6	0.559	0.541
7	0.563	0.539
8	0.560	0.540
9	0.564	0.539
10	0.561	0.536
mean	0.569	0.539
st. deviation	0.015	0.002
coef. variation	0.026	0.004

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.

5.4 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 rm =0.02 m, standard deviation σr = 0.1rm m and truncation limits rmin = rm - 3σr = 0.014 m and rmax = rm + 3σr = 0.026 m (this 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 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.

Figure 10
3D packs of inhomogeneous particles (rotational motion is ON).

Figure 11
3D packs of inhomogeneous particles (rotational motion is OFF).

Table 4
Packing densities (ϕ) obtained for example 5.4.

Pack	Rotations ON	Rotations OFF
1	0.578	0.548
2	0.563	0.543
3	0.566	0.552
4	0.571	0.547
5	0.567	0.542
6	0.574	0.552
7	0.578	0.553
8	0.575	0.542
9	0.578	0.546
10	0.578	0.538
mean	0.573	0.546
st. deviation	0.006	0.005
coef. variation	0.010	0.009

6 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.

Acknowledgements

This work was supported by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) and FAPESP (São Paulo State Science Foundation), under the grants 303793/2012-0 and 2014/19821-7, respectively. Material support and the stimulating discussions in CMRL (Computational Materials Research Laboratory, University of California at Berkeley, USA) are also gratefully acknowledged.

References

• Batten, R.D., Stillinger, F.H., Torquato, S. (2011). Inherent structures for soft long-range interactions in two-dimensional many-particle systems, J Chem Phys 135:054104.
• Bicanic, N. (2004). Discrete Element Methods, Encyclopedia of Computational Mechanics, Volume 1: Fundamentals. Stein E, de Borst R, Hughes TJR (eds), John Wiley & Sons (Chichester).
• Bolintineanu, D.S., Grest, G.S., Lechman, J.B., Pierce, F., Plimpton, S.J., Schunk, P.R. (2014). Particle dynamics modeling methods for colloid suspensions, Comp Part Mech; 1(3):321-356.
• Campello, E.M.B., Zohdi, T.I. (2014). A computational framework for simulation of the delivery of substances into cells, Int J Numer Meth Biomed Engng 30:1132-1152.
• Campello, E.M.B., Zohdi, T.I. (2014). Design evaluation of a particle bombardment system used to deliver substances into cells, Comp Model Engng Sci 98(2):221-245.
• Campello, E.M.B. (2014). Computational modeling and simulation of rupture of membranes and thin films, J Braz Soc Mech Sci Engrg, DOI 10.1007/s40430-014-0273-5.
  » https://doi.org/10.1007/s40430-014-0273-5
• Campello, E.M.B. (2015). A description of rotations for DEM models of particle systems, Comp Part Mech 2:109-125.
• Conway, J., Sloane, N. (1998). Sphere Packings, Lattices and Groups, Springer (Berlin).
• Cui, L., O'Sullivan, C. (2003). Analysis of a triangulation based approach for specimen generation for discrete element simulations, Granular Matter 5:135-145.
• Cundall, P.A., Strack, O.D.L. (1979). A discrete numerical model for granular assemblies, Géotech 29:47-65.
• David, C.T., Garcia-Rojo, R., Herrmann, H.J., Luding, S. (2007). Powder flow testing with 2D and 3D biaxial and triaxial simulations, Part and Part Sys Char 2007; 24:29-33.
• Donev, A. (2006). Jammed packings of hard particles, PhD Dissertation, Princeton University.
• Duran, J. (1997). Sands, Powders and Grains: An introduction to the physics of granular matter, Springer (New York).
• Feng, Y.T., Han, K., Owen, D.R.J. (2003). Filling domains with disks: an advancing front approach, Int J Numer Mth Engng 56:699-713.
• Ferrez, J.A. (2001). Dynamic Triangulations for Efficient 3-D simulation of Granular Materials, Ph. D. thesis, Ecole Polytechnique Federal de Lausanne, Lausanne.
• Hales, T.C. (2005). A proof of the Kepler conjecture, Ann Math 162:1065-1185.
• Han, K., Feng, Y.T., Owen, D.R.J. (2005) Sphere packing with geometric based compression algorithm, Powder Technology 155:33-41.
• Jiang, M.J., Konrad, J.M., Leroueil, S. (2003). An efficient technique for generating homogenous specimens for DEM studies, Comput Geotech 30:579-597.
• Johnson, K.L., Kendall, K., Roberts, A.D. (1971). Surface energy and the contact of elastic solids, Proc R Soc Lond A 324:301-313.
• Johnson, K.L. (1985). Contact Mechanics, Cambridge University Press (Cambridge).
• Kamrin, K., Koval, G. (2014). Effect of particle surface friction on nonlocal constitutive behavior of flowing granular media, Comp Part Mech 1:169-176.
• Kansal, A.R., Torquato, S., Stillinger, S.H. (2002). Diversity of order and densities in jammed hard-particle packings, Phys Rev E 66:041109.
• Kroupa, M., Klejch, M., Vonka, M., Kosek, J. (2012). Discrete element modeling (DEM) of agglomeration of polymer particles, Proc Eng 42:58-69.
• Latham, J-P., Munjiza, A., Garcia, X., Xiang, J., Guises, R. (2008). Three-dimensional particle shape acquisition and use of shape library for DEM and FEM/DEM simulation, Minerals Engineering 21:797:805.
• Lubachevsky, B.D., Stillinger, F.H. (1990). Geometric properties of random disk packings, J Stat Phys 60:561-583.
• O'Sullivan, C. (2011) Particle-based discrete element modeling: Geomechanics perspective, Int J Geomech 11:449-464.
• McNamara, S., Young, W.R. (1994) Inelastic collapse in two dimensions, Phys Rev E 50:R28-R31.
• Pöschel, T., Schwager, T. (2004). Computational Granular Dynamics, Springer (Berlin).
• Scardicchio, A.,Stillinger, F.H., Torquato, S. (2008). Estimates of the optimal density of sphere packings in high dimensions, J Math Phys 49:043301.
• Siiria, S., Yliruusi, J. (2007). Particle packing simulations based on Newtonian mechanics, Powder Tech 174:82-92.
• Speedy, R.J. (1998). Random jammed packings of hard discs and spheres, J Phys: Condens Matter 10:4185.
• Song, C., Wang, P., Makse, H.A. (2008) A phase diagram for jammed matter, Nature 453(7195): 629-632.
• Stroeven, P., Stroeven, M. (1999). Assessment of packing characteristics by computer simulation, Cem Concr Research 29:1201-1206.
• Tan, Y., Yang, D., Sheng, Y. (2009). Discrete element method (DEM) modeling of fracture and damage in the machining process of polycrystalline SiC, J Euro Ceram Soc 29(6):1029-1037.
• Thornton, C., Gong, G., Chan, A.H.C. (2012). DEM Simulations of Undrained Triaxial Behavior of Granular Material, J Eng Mech 138:560-566.
• Torquato, S., Truskett, T.M., Debenedetti, P.G. (2000). Is random close packing of spheres well defined?, Phys Rev Lett 84(10):2064-2067.
• Thornton, C. (1999). Three-dimensional behavior of granular materials, In: M. Oda and K. Iwashita (Eds.), Introduction to mechanics of granular materials, Chapter 3.3.5, p. 187-196, Balkema (Roterdam).
• Xu, W.X., Chen, H.S. (2013). Numerical investigation of effect of particle shape and particle size distribution on fresh cement paste microstructure via random sequential packing of dodecahedral cement particles, Comp Struct 114-115:35-45.
• Wellmann, C., Wriggers, P. (2012). A two-scale model of granular materials, Comput Methods Appl Mech Engrg 205-208:46-58.
• Zhang, G., Torquato, S. (2013) Precise algorithm to generate random sequential addition of hard hyperspheres at saturation, Phys Rev E 88:053312.
• Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, A.B. (2008). Discrete particle simulation of particulate systems: a review of major applications and findings, Chem Eng Sci 63:5728-5770.
• Zinchenko, A. (1994).Algorithm for random close packing of spheres with periodic boundary conditions, J Comput Phys 114:298-315.
• Zohdi, T.I. (2012). Dynamics of charged particulate systems: Modeling, theory and computation, Springer (New York).

Publication Dates

• Publication in this collection
  Jan 2016

History

• Received
  13 Nov 2014
• Reviewed
  21 Sept 2015
• Accepted
  21 Oct 2015