Effect of particle size distribution on 3D packings of spherical particles

. We use molecular dynamics simulations of frictionless spherical particles to investigate a class of polydisperse granular materials in which the particle size distribution is uniform in particle volumes. The particles are assembled in a box by uniaxial compaction under the action of a constant stress. Due to the absence of friction and the nature of size distribution, the generated packings have the highest packing fraction at a given size span, deﬁned as the ratio α of the largest size to the smallest size. We ﬁnd that, up to α = 5, the packing fraction is a nearly linear function of α . While the coordination number is nearly constant due to the isostatic nature of the packings, we show that the connectivity of the particles evolves with α . In particular, the proportion of particles with 4 contacts represents the largest proportion of particles mostly of small size. We argue that this particular class of particles occurs as a result of the high stability of local conﬁgurations in which a small particle is stuck by four larger particles.


Introduction
Granular materials with broad particle size distributions are very common in nature and industrial applications.Broad size distributions often arise as a result of particle breakage [1][2][3][4] and strongly influence the space-filling and strength properties of granular materials [5][6][7][8][9][10][11].Despite their importance, most studies of packing structure are devoted either to weakly polydisperse systems [12,13] or to ideal apollonian or random apollonian constructions characterized by the assumption that there is no lower bound on particle size [7,8,14].Real granular materials may, however, have a variety of size distribution functions and bounds on the particle size.Many geometrical algorithms have been proposed to generate jammed packings of particles [7,9,10,[15][16][17][18].In contrast, dynamical simulations of polydisperse packings are quite rare as the generation of statistically representative samples of different size classes in a packing requires increasingly large numbers of particles as the ratio α between the largest and smallest sizes increases [18].
In this paper, we investigate polydisperse packings of spherical particles by means of molecular dynamics simulations up to a ratio α = 5.Since we are interested in well-defined granular states that do not depend on the construction method, we consider only frictionless particles, that assemble into an isostatic packing with the largest possible packing fraction up to possible finite size effects.The absence of friction leads also to a nearly isotropic packing.We first describe the system characteristics and numerical procedures.Then, we analyze the geometrical properties of the packings such as packing fraction and connectivity of the particles.

Method and system parameters
We consider size distributions with uniform volume fractions in a range [d min , d max ] of particle diameters.This means that the normalized distribution P d (d) of particle diameters varies as d −3 so that the volume V d is the same and equal to the total volume V divided by the total number N p of particles for each size d.Hence, the normalized distribution is Since only size ratios matter, we normalize all particle diameters by d min : where d r is the reduced size and α is the 'size ratio'.The distribution of reduced diameters is given by In practice, for a discrete system with a finite number of particles, this distribution must be discretized by dividing the range d r ∈ [1, α] into N c size classes.Each class i contains N p/c (i) particles.The total number of particles is N p = i N p/c (i).Since the system is polydisperse, for a given number N p of particles, the number N c of size classes and the values of N p/c (i) must be such that in each class both particle volumes and their numbers are correctly represented.This requires a large enough number of largest particles N min p/c and large enough volume of smallest particles.By choosing N min p/c = 50, and given that the volume is the same in all size classes, we get the minimum number of particles N min p required to construct the packing by summing up the values of N min p/c .The values are given in Table 1 where we see that the largest number of particles requires to ensure a minimum of 50 particles in each class is below 11,000 particles for α = 5 and N c = 10.Hence, for a representative size distribution, we used 25,000 particles in our simulations for 5 values of α varying from 1 (monodisperse) to 5.
The particles are smooth spheres (friction coefficient μ = 0).They are initially placed randomly inside a rectangular cell and subjected to uniaxial compaction by applying a constant load σ zz on the top wall, the bottom wall imposed to be immobile, and using periodic boundary conditions along x and y directions.The normal contact force f n between particles is governed by a linear visco-elastic law.Hence, the energy is only dissipated by inelastic collisions between particles and with the top and bottom walls.The compaction is fairly fast at the beginning but considerably slows down in approaching the jamming point.We monitor the total kinetic energy, contact network and packing fraction Φ.When kinetic energy is low compared to σ zz V, which represents the potential energy of the system, and the packing fraction and contact network are stable, we stop the simulation.The simulation parameters are given in Table 2.An important parameter is the maximum contact deflection δ max .We set the applied pressure to σ zz = 2 × 10 5 Pa, so that δ max /d max ∼ 10 −3 .This means that the particles can be considered as nearly undeformable and the calculated values of the packing fraction are not influenced by contact deflections.The 5 samples generated for increasing value of α are denoted as S1, S2, S3, S4 and S5.A snapshot of a portion of the packing S4 is displayed in Fig. 1.The force chains are displayed for a portion of the packings S1 and S4 in Fig. 2.

Packing fraction
The packing fraction Φ is calculated by simply dividing the total volume of the particles by the total volume of the packing.Fig. 3 shows Φ as a function of α.We see that for considered range of size polydispersity, Φ is a nearly linear function of α.Consistently with the measured values of the random close packing fraction, we get Φ 0.643 in the monodisperse case (α = 1) [19][20][21][22].For α = 5, Φ is as large as 0.71.Obviously, the linear relationship can not be extended to much larger values of α as the packing fraction can not exceed 1.It can thus be the beginning of an exponential relation, for example, which tends to 1 as α → ∞.Note also that for a size ratio α ≥ 4, a sphere of size d r = 1 can be fit inside a pore between four spheres of size d r = α.But a close examination of the samples reveals that there is almost no such floating particles in samples S4 and S5.This is mainly because, given the much lower proportion of large particles, the probability of a configuration with four such touching particles is quite low.But, it is expected that for larger values of α, the number of small floating particles increases, as it was also shown in 2D packings of disks [18].

Contact networks
The contact network can be analyzed using various descriptors.
We analyze only partially here the connectivity of the particles described by the proportion P dc (d r , c) of particles of reduced size d r and having exactly c contacts.The connectivity of the particles with no distinction between particle sizes, is P c (c) = d r P dc (d r , c), which is simply the proportion of particle with c contacts.The coordination number is given by Z = c P c (c).In the absence of friction, our packings are isostatic and thus Z is nearly equal to 6 in all our packings.But the connectivity depends on α.Fig. 4 displays the distributions P c (c) for all values of α.In S1, we observe a peak on c = 6, but we have a large number of particles with c = 4, c = 5 and c = 7 contacts.We also observe that P c (3) 0 for all values of α.Indeed, for c = 3, force balance is possible only if the three force vectors are on the same plane, which is quite improbable.In the same way, force balance for c = 2 implies equal forces on diametrically opposite points, that is a configuration of zero chance to occur.
As α increases, we observe an decreasing number of particles with c > 6.Interestingly, the peak value of P c tends to c = 4 and P c (4) increases from α = 3 to α = 5.For α = 5, the proportion of particles with c = 4 contacts is almost half the total number of contacts.This observation indicates that the configurations in which a particle is equilibrated by 4 neighboring contact particles  are stable in the mechanical sense.But since the peak value increases with α, it is expected that all such particles are of small size.This is indeed the case, as we see in Fig. 5 where for α = 5 we have plotted the distribution P dc (d r , c = 4), representing the proportion of particles with 4, 5 and 6 contacts as a function of their reduced size d r .The largest proportion (above 30%) of all particles is for c = 4 and for small particles (d r = 1).
The presence of such a high number of small particles stuck by four larger neighboring particles can be understood by the following argument.When a small particle (the one that is stuck between four particles) is in contact with three larger particles, the orientations of the normals of the three contacts (that are also the orientations of the forces in the absence of friction) define a broader range than when the particle is of the same size as its neighbors; see Fig. 6.For force balance, the force vector of the forth contact should be in this range but in the opposite direction.In the limit where the diameter of the central particle is smaller than 1/4 the large particle diameters, it can perfectly fit into the pore.The evolutions of P c (4), P c (5) and P c (6) are displayed as a function of α in Fig. 7.We see that P c (5) is nearly constant and represents  20% of particles whereas P c (4) increases at the expense of P c (6) and particles of higher number of contacts.

Conclusions
In this paper, we studied polydisperse packings of frictionless spherical particles by means of molecular dynamics simulations.We showed that, for a uniform size distribution in particle volume fractions, the packing fraction is a nearly linear function of size ratio between largest and smallest particles up to a ratio 5.While the coordination number is nearly 6 due to the isostatic nature of the packings, the structure of the contact network evolves with size ratio.In particular, we evidenced the existence of a large proportion of particles with exactly four contacts for size ratios above 3.The force distributions and their correlations with the microstructure are other interesting features of polydisperse packings that we will investigated in a future paper.

Figure 1 .
Figure 1.Snapshot of the S4 packing.Gray level is proportional to particle size.

Figure 2 .
Figure 2. Force chains in packings S1 (a) and S4 (b).Line thickness is proportional to normal force.

Figure 3 .
Figure 3. Packing fraction as a function of size ratio α for the packings.

Figure 4 .
Figure 4.The connectivity P(c) of the contact network for the packings S1-S5.

Figure 5 .
Figure 5.The proportions of particles with c = 4, 5 and 6 contacts as a function of their reduced size d r for α = 5.

Figure 6 .
Figure 6.A small particle stuck by four larger particles.The range of possible positions of the top particle for the equilibrium of the central particles depends on the size ratio.

Figure 7 .
Figure 7.The connectivity of the contact network for the packings S1-S5, illustrating the variation of P c (4), P c (5), and P c(6).

Table 1 .
Dispersity properties of the samples considered in this study.
† For polydisperse samples, uniform distributions in volume is used.‡ Min.number of particles needed for satisfying N min p/c = 50 and N c = 10.

Table 2 .
Parameters used in the DEM simulations.