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.

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^{5} 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
_{1},...,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′^{2} = 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):

| (2) |

| (3) |

The inverse of the partition function is given by:

| (4) |

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:

| (5) |

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:

- The number of particles in each class N
_{p∕c}^{′}is above a given threshold N_{p∕c}^{′ min}; - 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}^{i},d_{max}^{i}] of amplitudes Δd_{i} ≡ d_{max}^{i} -d_{min}^{i} with i ∈ [1,N_{c}]. Then,
we uncumulate h(d) on each sub-interval in order to obtain the volume fraction f_{v}^{i} of
each size class:

| (6) |

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}^{i}(d). This corresponds to approximating the
function h(d) by N_{c} segments passing through the points (d_{min}^{i},h(d_{min}^{i})) and
(d_{max}^{i},h(d_{max}^{i})). Considering the mean diameter in each class d_{m}^{i} = F_{uv}^{-1 i}(0.5),
we can evaluate a non-integer amount n^{i} of particles corresponding to a unit total
volume:

| (7) |

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}^{′}:

| (8) |

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}^{Nc} 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:

| (9) |

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:

| (10) |

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}^{i}(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:

| (11) |

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:

| (12) |

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

| (13) |

where Γ is the gamma function defined by:

| (14) |

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:

| (15) |

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:

| (16) |

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:

| (17) |

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:

| (18) |

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

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:

| (21) |

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}).

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^{6}, 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.

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^{4} and
N_{p}^{max} = 10^{5}. 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^{5}) 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)^{n}. For insertion in a triangular
lattice, the limit size ratio is R_{t} = 6,46^{n}. 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.