1 Size polydispersity

1.1 Polydispersity and packing

The size polydispersity is characterized by the proportions of particles of different sizes. Depending on the experimental method, the polydispersity is measured either by the number f(d) or the volume V (d) of particles as a function of their size d. The polydispersity can be represented by a continuous function, but in practice the fraction of particles belonging to a finite size interval [d - δd,d + δd] is determined and referred to as a size class. Integrating the size distribution f(d) we get the partition function F(d) whereas the integration of V (d) yields the cumulative volume fraction or the classical grading curve h(d).

Two approaches should be distinguished for the management of particle size distribution. The first approach consists of adapting the particle size to the assembling method. An example of this approach is the construction of apollonian packings. This geometric construction is based on the fact that three tangent and non-intersecting circles have two other tangent circles. Starting with three identical and tangent circles and iterating indefinitely this construction rule, one obtains a packing presented in figure 2 which tends to fill the space with increasingly smaller particles inserted. The resulting size distribution is a power law and the structure is fractal.

The second approach is to adapt the assembling method to a given size distribution. This is the natural way to construct random packings from a given size distribution. For example, the grains are poured into a box by pluviation to obtain a loose packing and the box is shaken to get a denser packing. Such a method may lead to size segregation with inhomogeneous distribution of different particle sizes inside the packing. In the same way, the geometrical assembling methods may lead to size segregation. For example, the application of the potential minimization method or ballistic deposition to a 3D system leads to the accumulation of small particles at the bottom of the packing as the small particles can flow in the pore space between large particles. Moreover, the computational cost increases with size span. Since segregation is a real physical effect, the MC method and dynamic simulations lead to similar segregation effects. A 3D packing devoid of segregation can be obtained only with cohesion or by means of the CA method.

1.2 Generation of a discrete ensemble of particle sizes

During a granulometric analysis or a laboratory mechanical test, the sample must be representative of particle properties including the particle size distribution both in volume (mainly for the smallest particles) and in number (mainly for the largest particles). This representativity of PSD is actually one of the conditions for a sample to be a representative volume element (RVE). However, in discrete numerical simulations of polydisperse granular samples this condition can not be always satisfied because of the computational cost imposing a limited an upper limit on the number of particles (presently 10⁵ particles for month-long simulations). For this reason, one needs appropriate methodology for optimal control of the statistical representativity of size classes. The issue is to extract a discrete ensemble of particle sizes from a measurement or an analytical function with a given number of particles and such that the distribution is optimally represented.

A classical approach is the inverse transformation method which can be used to generate a sample of n events from a continuous random variable X with a given cumulative distribution function F and a probability density function f = . In our case, the random variable X is the particle size d in the range [d_min,d_max]. The main steps of the method are the following:

• Let Y be a random variable whose distribution can be described by the cumulative distribution function Y = F(X) . This is a random variable with a uniform distribution in [0,1]. Its partition function G is defined by G(y) = P[Y < y] = y.
• Let y₁,...,y_n a sample of size n of a random variable uniformly distributed in [0,1]. The values y_i can be considered as a realization of the random variable Y .
• To compute the realizations x_i, it is sufficient to compute the value corresponding to the image of y_i by its partition function:

  (1)

For example, let us apply this method to a size distribution uniform by volume fraction. This distribution f_uv(d) is defined by f_uv(d′) d′² = constant, i.e. all size classes have the same volume of particles. Thus, the grading curve is a straight line passing through the points (d_min,0) and (d_max,1). Moreover, f_uv is a probability density function so that its integral is 1. We can easily deduce the expressions (2) of f_uv(d) and (3) from its partition function F_uv(d):

The inverse of the partition function is given by:

It is then sufficient to generate the random numbers y_i between 0 and 1 and compute their images by F_uv^-1(y_i) to obtain n particle sizes satisfying a distribution uniform by volume fraction.

This method is easy to implement since it only requires a generator of random numbers. Other more efficient methods exist for the generation of nonuniform distributions ?. But, the simple method presented above is sufficient since the generation is not a critical stage in terms of numerical efficiency. The point is that the application of this method requires the analytical expression of the size partition function F(d) and this function should be invertible. Moreover, the statistical nature of this method ensures only the statistical representativity in the number of particles without care for particle volumes. The number of n generated particles must be large enough to ensure also a representative volume for each size class.

To generate a discrete set of particle sizes from the cumulative volume distribution h(d) (grading curve), we first need to evaluate its partition function. Except for simple distributions such as the uniform distribution by particle sizes, this stage can be rather complicated for some distributions or when no analytical expression is available. By definition, we have:

where V (x) is the volume of a particle of size x. We clearly see that relating F(d) and h(d) explicitly is not straightforward in most cases. To tackle this delicate step, we present below a generation method based on the inverse sampling method.

The main goals of this method are 1) to evaluate implicitly the partition function F(d) of a given cumulative volume distribution and 2) to generate a set of particle sizes satisfying criteria of representativity in volume and in number, both at the level of the distribution (global scale) and at the level of each size class (local scale). We consider two explicit representativity criteria:

1. The number of particles in each class N_p∕c^′ is above a given threshold N_p∕c^′ min;
2. In each class, the volume of one particle is a small fraction of that of the whole class.

For a nearly monodisperse size distribution, these two conditions are equivalent. But, in a polydisperse case these conditions are critical for the less populated classes which are often the largest size classes. A high value of N_p∕c^′ min ensures that the rare particles are statistically represented and the volume of a single particle is small compared to that of the whole class.

We assume that the cumulative volume distribution h(d) is an increasing function satisfying h(d_min) = 0 and h(d_max) = 1 which reflect the fact that the whole volume of the distribution is composed of particles of size between d_min and d_max. The first step is to discretize the interval [d_min,d_max] in N_c size classes composed of the sub-intervals [d_minⁱ,d_maxⁱ] of amplitudes Δd_i ≡ d_maxⁱ -d_minⁱ with i ∈ [1,N_c]. Then, we uncumulate h(d) on each sub-interval in order to obtain the volume fraction f_vⁱ of each size class:

In the case of a discontinuous volume distribution h(d), the discretization step can be avoided if the discontinuities can be considered as size classes. For example, the grading curve obtained by sieving is, by essence, discontinuous and represented by a stair function. The consecutive sizes of the sieve meshes appears then as a natural discretization of the particle size span.

To evaluate the number of particles in each class, we assume that the size distribution in each size class follows a distribution uniform by volume fraction normalized on each sub-interval F_uvⁱ(d). This corresponds to approximating the function h(d) by N_c segments passing through the points (d_minⁱ,h(d_minⁱ)) and (d_maxⁱ,h(d_maxⁱ)). Considering the mean diameter in each class d_mⁱ = F_uv^-1 i(0.5), we can evaluate a non-integer amount nⁱ of particles corresponding to a unit total volume:

We then multiply the quantity n_i by the factor N_p∕c^′ min∕min(n_i) where min(n_i) represents the value of n_i in the less populated class. We extract the integer part in order to obtain a first estimation of the number of particles in each size class N_i^′:

This step ensures that the less populated size class contains at least N_p∕c^′ min particles. Therefore, the estimated total number of particles N_p^′ = ∑ _i=1^N_c N_i^′ depends both on the number of classes N_c and on the minimum number of particles per class N_p∕c^′ min. For a couple of chosen values of N_c and N_p∕c^′ min, N_p^′ can be inadequate to ensure a given level of global representativity. To control this point, we introduce the minimal total number of particles N_p^min. Hence, if N_p^′ < N_p^min, we multiply the number N_i^′ by the factor N_p^min∕N_p^′ to obtain the number of particles N_i in each size class i:

The N_i values finally represent an estimation of the number of particles needed to satisfy an imposed level of representativity. Then, we deduce the probability density P_i(d) of the sizes d for each class i:

The corresponding partition function F(d) is obtained by summing up P_i(d) for all size classes below d. In order to generate the final discrete size set, we perform in each class N_i inverse sampling of a uniform distribution by volume fraction normalized on each size class F_uvⁱ(d).

Besides the number of size classes N_c, the discretization properties used for the uncumulation can impact the quality of the generation procedure. In fact, one type of discretization can help to satisfy the representativity conditions while another type may lead on the contrary to its violation. For example, we consider a regular discretization with n classes of equal widths over [d_min,d_max]. A large number of size classes will provide a good approximation of the curve but the mass fraction represented by each class will decrease as n increases. This point is crucial for the class of largest particles because the whole mass fraction can be represented by a unique particle and thus contradict the second representativity criteria. To overcome this difficulty, we have introduced the minimal number of particles per class. One way to improve this point is to use an irregular discretization in which the class width increases with particle size. For the same distribution and number of classes, with an irregular discretization the number of largest particles increases as the mass fraction of largest classes increases.

A genuine irregular discretization is one which provides, after uncumulation, an equal mass fraction for each size class. We call this discretization isochoric. Taking care of the choose of a reasonable mass fraction, the most populated classes will be well represented. On the contrary, for the size distributions where the class of the largest particles represents a small mass fraction, this kind of discretization is problematic. In this case, the largest particles will be represented by a unique class of large width in which the statistical representativity will be unsatisfactory. Another inconvenient aspect of this discretization is its dependence with respect to the shape of the size distribution. An alternative method consists of using of an irregular discretization depending only on the number of size classes N_c. In analogy to the successive sieve meshes used in laboratory, the successive boundaries of size classes may follow a geometric progression of first term d_min and of common ratio r given by:

2 Polydisperse packings

2.1 Analytical model of polydispersity

For the numerical simulations, an analytic model of the grading curve is useful. Such model must allow one to vary at least two basic properties of polydispersity: the size span and the shape. The size span corresponds to the ratio of the largest and smallest sizes. The shape describes the relative weights of these two classes. We will specify in more detail these features below.

For practical reasons, a model of grading curve should

• be simple, i.e. involving a small number of parameters to change the size span and the shape;
• be defined over a finite interval;
• contain classical distributions such as power laws;
• be capable of representing double-curved distributions as observed in applications, for example the double curved shape frequently observed in geomaterials.

Voivret et al. have shown that the β distribution in its cumulative form β distribution is a good model meeting all the above requirements ?. It is defined by:

where a > 0 and b > 0 are the two parameters of the distribution. The B(a,b) function is given by:

where Γ is the gamma function defined by:

The cumulative β distribution (CB) is defined and normalized on the interval [0,1] with β(0) = 0 and β(1) = 1. In order to represent the grading curve of particle size in the range [d_min,d_max] with CB, we replace the argument x by the reduced diameter d_r defined by:

which varies in the range [0,1]. With this definition, the model grading curve h(d) is given by the CB distribution of the reduced diameters:

The parameters a and b of the CB model allow us to easily vary the shape of the distribution h(d). On the other hand, the size span can be characterized by the ratio R = d_max∕d_min. But as for d_r, we can also define the span by the parameter s:

The case s = 0 corresponds to a strictly monodisperse distribution while the case s = 1 corresponds to an infinitely polydisperse size distribution.

The CB model is simple since it only requires three parameters to vary the shape and span of the size distribution. Moreover, only two a and b are used to produce a rich set of various distribution shapes. Depending on the values of a and b, we have a grading curve which is

• a < 1, b < 1 U-shaped
• a < 1, b ≥ 1 ou a = 1, b > 1 strictly decreasing
  • a = 1, b > 2 strictly convex
  • a = 1, b = 2 a straight line
  • a = 1, 1 < b < 2 strictly concave
• a = 1, b = 1 uniform ove [0;1]
• a = 1, b < 1 ou a > 1, b ≤ 1 strictly increasing
  • a > 2, b = 1 strictly convex
  • a = 2, b = 1 a straight line
  • 1 < a < 2, b = 1 strictly concave
• a > 1, b > 1 unimodal
• a = b symmetric with respect to 1∕2

Moreover, the CB model contains classical size distributions some of which are displayed in figure 3 ):

• a = b = 1 uniform by volume fraction
• a = 3, b = 1 uniform by particle sizes
• a = 1, b > 1 strictly convex power law
• a > 1, b = 1 strictly concave power law
• a > 1, b > 1 double curved
• a < 1, b < 1 tends to the bi-disperse case as a and b decrease, the particles sizes being concentrated around d_r = 0 et d_r = 1.

In soil mechanics, two descriptors are frequently used to characterize the soils: the uniformity coefficient C_u and the curvature coefficient C_c also called the Hazen coefficients. They are defined by:

where d_x designs the diameter for which the volume fraction of smaller particles is equal to x%. These coefficients are the lowest-order descriptors of the grading curve such like they can be easily extracted from the relations (18). It is clear that the descriptors of the CB model are appropriate for a discrete size set defined over an interval [d_min,d_max]. In a real granular material, these extreme size are not well defined and it is easier to use the C_u and C_c coefficients to characterize the polydispersity. Nevertheless, these coefficients can be expressed as a function of the CB model parameters. Since d_x is directly given by the inverse function of the grading curve, the C_u and C_c coefficients are given by:

Figure 5 shows these coefficients as a function of s for different values of shape parameters. We remark that for all values of shape parameters, the uniformity coefficient C_u is a monotonic function of s. In other words, C_u does represent the ratio of the characteristic sizes of the largest and smallest size classes. On the other hand, the curvature coefficient C_c is a complex function of s and shape parameters a and b. Nevertheless, this relation between the granulometric coefficients and the CB model parameters allows one to determine the values of the shape and span parameters from the measured values of C_u et C_c.

2.2 Representativity criteria

The size generation method requires three parameters, namely the number of classes N_c, the minimum number of particles per class N_p∕c^′ min and the minimum total number of particles N_p^min. As our main goal is to generate representative sets of particle sizes with a minimum number of particles, we now focus on the influence of these parameters on the quality of the generated set. To evaluate the quality, we will compare the resulting grading curve y^exp(d) with the analytic function h(d) with the help of the regression coefficient C_r defined over k points by:

The k sizes d_i used in the comparison follow a regular geometric discretization between d_min and d_max. A value of C_r close to 0 indicates a good agreement between the analytic and generated distributions. To investigate the influence of the generation parameters, we focus on three different size distributions (a = 1;b = 3; a = b = 3;a = 3;b = 1) with the same size span s = 0.98, which correspond to R = 100. We fix two parameters while the third varies. The default values of the fixed parameters are N_c = 10, N_p∕c^′ min = 10 and N_p^min = 2000.

The number of size classes represents the number of segments with which we will approximate the cumulative volume fraction h(d). For a given discretization, it is clear that a large number of classes will provide a good approximation of the continuous curve. Figure 2.2 shows that ten size classes provide a good agreement between the discrete set and h(d). Figure 2.2 shows the total number of particles as a function of the number of size classes. For each distribution, we observe the effect of the minimum total number of particles for a small number of classes. When N_c is sufficiently large, the number of generated particles strongly increases due to the condition of a minimum number of particles per class.

Figure 2.2 shows the evolution of C_r with N_p∕c^′ min. The impact of N_p∕c^′ min mainly shows up on the distribution a = 1;b = 3. Since the volume fraction of the largest particles is small, a large value of N_p∕c^′ min simultaneously satisfies both conditions in extreme classes. In fact, it can be said that N_p∕c^′ min acts on the extreme size classes of the size distribution. If N_i is a decreasing (resp. increasing) function of the particle size d_i, the class of the largest (resp. smallest) particles contains N_p∕c^′ min particles. It means that, if the extreme sizes are well represented, then the whole size distribution is also well represented.

While N_p∕c^′ min sets the local representativity, N_p^min fixes the total number of particles for which we consider that a distribution is globally representative. We see in figure 2.2 that the value of N_p^min has a weak impact on the quality of the generated distribution. Therefore, we must rather consider N_p^min as the lower bound depending on the usage of the distribution. In figure 2.2 we see two variation steps for each distribution: a plateau where N_p is independent of N_p^min and a linear evolution where N_p = N_p^min. The value of N_p associated with each plateau corresponds to the number of particles N_i^′ needed to satisfy the relation N_i^′ > N_p^min. The linear evolution can be explained by the fact that when this relation is no more fufilled, the total number of particles is equal to its minimum value (N_p = N_p^min).

2.3 Numerically accessible distributions

Since the total minimum number N_p^min is only a lower bound, the total number of generated particles N_p may reach very large values depending on the others generation parameters and the shape parameters. It is thus necessary to introduce an “accessibility” limit N_p^max defined as the maximum number of particles acceptable for a given usage of the distribution. Presently, ten thousands seems to be an acceptable number of particles for performing dynamic simulations on a weekly basis with a moderate-power computer. If N_p^min ≤ N_p ≤ N_p^max, the distribution is both representative and accessible. For constant shape parameters, we seek a good compromise between N_c and N_p∕c^min allowing us to generate a distribution with a total number of particles between N_p^min and N_p^max.

For constant generation parameters (N_c, N_p∕c^min,N_p^min et N_p^max), all values of the distribution shapes a and b are not accessible. The accessibility limit defines thus the space of the accessible values of the shape parameters. Figure 8 shows the space of accessible values for two values of size span and for different values of the accessibility limit N_p^max. It is remarkable that the general aspect of the accessible domain is weakly dependent on the parameters and it expands slowly with N_p. For N_p < 10⁶, the values larger than a = b = 10 can not be realized. On the other hand, for strictly convex distribution (a = 1), the limit value of b is about 4 whereas for strictly concave distribution (b = 1) the limit value of a is about 7. The shape of the accessible domain is not symmetric about the line a = b. The accessible values of a and b are all the more reduced as the small particle classes are more populated. This restriction is more important for the points located above the line a = b, i.e. for a < b. We note that uniform distributions in diameter (a = 3 and b = 1) and by volume fractions (a = b = 1) are inside the accessible domain.

2.4 Size polydispersity and packing fraction

The texture of a packing assembled by a geometric method depends on the polydispersity. We focus here on the influence of size polydispersity on the space-filling properties in 2D and optimal values of the packing fraction ??. The generation parameters are set to N_c = 10, N_p∕c^min = 10, N_p^min = 2.5 10⁴ and N_p^max = 10⁵. The size span varies from 0.02 to 0.97. We use the potential minimization method to assemble the particles.

We consider six size distributions with the following shape parameters : (a = b = 1); (a = 1;b = 3); (a = 2;b = 4); (a = 3;b = 1); (a = 4;b = 2) and (a = b = 4) (see figure 3). The accessible space parameter is similar as that shown in figure 8. Each plot shown below represents about 50 data points, each point corresponding to the average of two independent depositions. For some plots (9), the data points for large values of s are not statistically accessible (for N_p^max < 10⁵) and thus these plots could not be traced for the whole range of s ∈ [0.02,0.97].

The packing fraction ρ is shown in Fig. 9 as a function of size span s for different values of shape parameters. We observe similar trends for all values of shape parameters a and b. The packing fraction declines in the range s < 0.1. Then, it increases slowly with s in the range 0.1 < s < 0.4 and finally increases strongly beyond s ≃ 0.4. Figures 2.4, 2.4 and2.4 show snapshots of the samples for three different values of s generated for a = b = 1. At very low values of s (s < 0.1), the arrangement is quasi-monodisperse involving long-range order and dislocations. The coordination number is about 4. The larger values of ρ (ρ > 0.82) in this range are due to local triangular structures inside a globally square lattice of particles. The decrease of rho with s in this range corresponds to the well-known order-disorder transition observed in 2D monodisperse granular media due to the perturbation of an initially triangular lattice in the presence of a small amount of polydispersity ?.

The initial decrease of ρ can be explained by the fact we leave a monodisperse size distribution. In 2D a strictly monodisperse granular material spontaneously arranges itself in an ordered and highly compact (ρ ≃ 0.906) structure with triangular cells. A small dispersion in particle size is sufficient to break this crystalline order in favor of a random close packing (order-disorder transition) ?? with a packing fraction close to 0.81. The values of the packing fraction between these two values are due to local triangular structures inside a globally disordered microstructure. Figure 2.4 displays a snapshot of a nearly monodisperse packing s = 0.02. Figure 2.4 shows the same map in which the zones with similar arrangements are marked. We observe large regions of rectangular arrangement (B) with long range correlations. The compact triangular regions (A) are less extended.

Apart from the uniform distribution by volume fractions, we see that the shape parameter has a weak influence on the packing fraction for the size span ranging from s = 0.1 to s = 0.4. Figure 2.4 shows a weakly polydisperse packing. The previous crystallized regions have been totally replaced by a dense arrangement in which only a short-range order persists. The increase of packing fraction is more due to the increasing disorder than to an efficient pore filling. Indeed, in this range of size span, the smallest particles are not small enough to fit into the pores between larger particles. This condition occurs precisely at s = 0.41 for a square lattice. This condition can be considered as a geometric transition in the texture of the packing. Of course, the largest particles do not form a square lattice but we see that this phenomenon limits the increase of the packing fraction below s ≃ 0.4. With the increase of the size span, this insertion mechanism occurs with generations of increasingly smaller particles. Inserting the n^ih size generation in a square lattice, requires a size ratio above R_c = (2,4)ⁿ. For insertion in a triangular lattice, the limit size ratio is R_t = 6,46ⁿ. However, due to the disorder only the effect related to the first size generation for a square lattice is clearly observed.

Beyond this value of size span, we distinguish two groups of size distributions for which the packing fraction does not increase in the same way. The distributions favoring the number of small particles (a = 1;b = 3, a = 2;b = 4) ensure a faster increase of packing fraction with s than those favoring the large particles (a = 3;b = 1, a = 4;b = 2). When the smallest particles prevail in number, they fill efficiently the voids. On the contrary, when they are in minority, the pores are closed before they are filled by small particles during the assembling process. In figure 2.4 we see that the voids between the largest particles are well filled and that the void sizes are much smaller than the surrounding particles. The distribution uniform by volume fractions is a particular case. For all values of size span s, the packing fraction is larger than that obtained by the others size distributions. This means that, this distribution ensures the best compromise between the volume of pores generated by the largest particles and the volume of the smallest particles to fill these pores.