Models of particle size distribution

Charles Voivret, Jean-Yves Delenne and Farhang Radjai

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


Figure 1: The size polydispersity of a granular material in which the grains have diameters ranging from 0.001 mm to 0.01 mm: (a) The probability density function and volume fraction functions; (b) The cumulative functions representing the partition function and the grading curve.


Figure 2: Apollonian packing

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 = ∂F-
∂X. In our case, the random variable X is the particle size d in the range [dmin,dmax]. The main steps of the method are the following:

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

        -dmaxdmin-- 1
fuv(d) = dmax - dmin d

           d        d- d
Fuv(d) =----max---- ----min-
        dmax - dmin    d

The inverse of the partition function is given by:

               dmax dmin
F -u1v (y) = d--+-y(d------d---)
          max     min    max

It is then sufficient to generate the random numbers yi between 0 and 1 and compute their images by Fuv-1(yi) 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:

h(d) = -dmin-V(x)f(x)dx-,

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 Np∕c is above a given threshold Np∕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 Np∕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(dmin) = 0 and h(dmax) = 1 which reflect the fact that the whole volume of the distribution is composed of particles of size between dmin and dmax. The first step is to discretize the interval [dmin,dmax] in Nc size classes composed of the sub-intervals [dmini,dmaxi] of amplitudes Δdi dmaxi -dmini with i [1,Nc]. Then, we uncumulate h(d) on each sub-interval in order to obtain the volume fraction fvi of each size class:

  i     i        i
fv = h(dmax)- h(dmin).

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 Fuvi(d). This corresponds to approximating the function h(d) by Nc segments passing through the points (dmini,h(dmini)) and (dmaxi,h(dmaxi)). Considering the mean diameter in each class dmi = Fuv-1 i(0.5), we can evaluate a non-integer amount ni of particles corresponding to a unit total volume:

 i   --fiv--
n  = V(dim ).

We then multiply the quantity ni by the factor Np∕c min∕min(ni) where min(ni) represents the value of ni 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 Ni:

      (   ′ min   )
N′=  E  -Np∕c---ni
 i      M in(ni)

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

      (         )
        N mpin  ′
Ni = E  -N-′- Ni

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

P(d) = Ni.
 i     Np

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

Besides the number of size classes Nc, 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 [dmin,dmax]. 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 Nc. In analogy to the successive sieve meshes used in laboratory, the successive boundaries of size classes may follow a geometric progression of first term dmin and of common ratio r given by:

       {                 }
r = exp         nc         .

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

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:

            1   ∫ x
β(x,a,b) = ------   ta-1(1 - t)b-1 dt,
          B(a,b) 0

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

B(a,b) = Γ (a)Γ (b)∕Γ (a+ b),

where Γ is the gamma function defined by:

       ∫ ∞ x-1 -t
Γ (x) = 0 t   e  dt.

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 [dmin,dmax] with CB, we replace the argument x by the reduced diameter dr defined by:

d (d) = -d---dmin--
 r     dmax - dmin

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:

h(d,a,b) = β(dr(d);a,b).

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 = dmax∕dmin. But as for dr, we can also define the span by the parameter s:

    d   - d
s = -max---min-.
    dmax + dmin

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

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


Figure 3: Examples of the cumulative volume distributions given by the analytic CB model.

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

      d        d2
Cu =  -60Cc = ---30--
      d10     d10d60

where dx 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 [dmin,dmax]. In a real granular material, these extreme size are not well defined and it is easier to use the Cu and Cc coefficients to characterize the polydispersity. Nevertheless, these coefficients can be expressed as a function of the CB model parameters. Since dx is directly given by the inverse function of the grading curve, the Cu and Cc coefficients are given by:

C   =   h-1(0,6)-                         (19)
 u      h-1(0,1)
          {h- 1(0,3)}2
Cc  =   h-1(0,6)-h-1(0,1)-                 (20)

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 Cu is a monotonic function of s. In other words, Cu does represent the ratio of the characteristic sizes of the largest and smallest size classes. On the other hand, the curvature coefficient Cc 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 Cu et Cc.


Figure 4: Evolution of the granulometric coefficients Cu and Cc with size span s and for different values of shape parameters a and b.

2.2 Representativity criteria

The size generation method requires three parameters, namely the number of classes Nc, the minimum number of particles per class Np∕c min and the minimum total number of particles Npmin. 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 yexp(d) with the analytic function h(d) with the help of the regression coefficient Cr defined over k points by:

     ┌│ ∑k-------------------
C =  │∘ --i=1(yexp(di)--h(di))2.
 r         ∑ki=1yexp(di)2

The k sizes di used in the comparison follow a regular geometric discretization between dmin and dmax. A value of Cr 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 Nc = 10, Np∕c min = 10 and Npmin = 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 Nc is sufficiently large, the number of generated particles strongly increases due to the condition of a minimum number of particles per class.


Figure 5: Influence of the number of size classes Nc on Cr (a) and on the total number of generated particles Np (b).

Figure 2.2 shows the evolution of Cr with Np∕c min. The impact of Np∕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 Np∕c min simultaneously satisfies both conditions in extreme classes. In fact, it can be said that Np∕c min acts on the extreme size classes of the size distribution. If Ni is a decreasing (resp. increasing) function of the particle size di, the class of the largest (resp. smallest) particles contains Np∕c min particles. It means that, if the extreme sizes are well represented, then the whole size distribution is also well represented.


Figure 6: Regression coefficient Cr as a function of the minimum number of particles per class Np∕c min (a), semi-log plot of the total number of particles Np as a function of Np∕c min (b).

While Np∕c min sets the local representativity, Npmin 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 Npmin has a weak impact on the quality of the generated distribution. Therefore, we must rather consider Npmin 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 Np is independent of Npmin and a linear evolution where Np = Npmin. The value of Np associated with each plateau corresponds to the number of particles Ni needed to satisfy the relation Ni > Npmin. 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 (Np = Npmin).


Figure 7: Semi-log plot of the regression coefficient Cr as a function of the minimum total number of particles Npmin (a); logarithmic plot of the total number of particles Np as a function of Npmin (b).

2.3 Numerically accessible distributions

Since the total minimum number Npmin is only a lower bound, the total number of generated particles Np may reach very large values depending on the others generation parameters and the shape parameters. It is thus necessary to introduce an “accessibility” limit Npmax 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 Npmin Np Npmax, the distribution is both representative and accessible. For constant shape parameters, we seek a good compromise between Nc and Np∕cmin allowing us to generate a distribution with a total number of particles between Npmin and Npmax.

For constant generation parameters (Nc, Np∕cmin,Npmin et Npmax), 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 Npmax. It is remarkable that the general aspect of the accessible domain is weakly dependent on the parameters and it expands slowly with Np. For Np < 106, 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.


Figure 8: Domains of statistically accessible shape parameters a and b for two different values of size span s = 0.1 (top) and s = 0.97 (bottom) and for several values of the maximum number of particles Npmax. (c) Accessible values of the uniformity coefficient for s = 0.97 and Np = 106.

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 Nc = 10, Np∕cmin = 10, Npmin = 2.5 104 and Npmax = 105. 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 Npmax < 105) and thus these plots could not be traced for the whole range of s [0.02,0.97].


Figure 9: Packing fraction ρ as a function of size span s for several values of shape parameters a and b.

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.


Figure 10: (a) Snapshot of a nearly monodisperse packing s = 0.02. (b) Regions of similar particle arrangements: A) triangular, B) rectangular and C) random.


Figure 11: Snapshot of two packing prepared by the potential-based method with the CB model : (a) s = 0.42 et (b) s = 0.97.

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 ni`t  h size generation in a square lattice, requires a size ratio above Rc = (2,4)n. For insertion in a triangular lattice, the limit size ratio is Rt = 6,46n. 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.