The rheology of granular materials is largely dominated by geometrical constraints arising from steric exclusions and involving particle shapes and size distributions. This prevailing role of geometry permits one to simplify sometimes the dynamics in favor of a better description of the geometry or/and higher numerical efficiency. For example, a dense granular packing may be efficiently constructed by replacing the equations of dynamics by simple displacement rules satisfying the geometrical constraints. Such purely geometrical procedures can be much simpler and numerically faster than more realistic procedures employed in dynamic or quasistatic methods.

The issue of the assembling methods is to construct configurations of particles as close as possible to a state of mechanical equilibrium with built-in packing properties. This can be, for example, a target packing fraction for a given particle size distribution. In the same way, the average connectivity of the particles (coordination number) and the anisotropy of the contact network are basic geometrical properties that control the mechanical response of a packing and may be built into a packing by an appropriate method. The homogeneity of the particle assembly in terms of packing fraction and connectivity is another important property which depends on the assembling rules. We introduce here several basic assembling methods by simple geometrical rules.

Contrary to dynamic simulation methods, the geometrical methods allow for quick
assembling of a large number of particles. Such packings may then be used as the
initial state for dynamic simulations. The geometrical methods help in this way to
improve numerical efficiency in the preparation phase. For example, gravitational
deposition of 10^{4} particles located initially on a regular grid can require
hours of computation whereas a nearly similar result may be obtained by
means of a geometrical method in only a few minutes. The drawback
is that the resulting sample will not be in mechanical equilibrium for
sure and no information is available on the contact forces. Nevertheless,
depending on the relaxation rule, the sample may still be sufficiently close
to equilibrium to be considered as a good starting point for mechanical
simulations.

In this section, we present the random deposition method which is the most intuitive method of particle assembling. This method, initiated by the work of Vold Vold [1959a,b, 1960], has inspired several algorithms among which one finds Visscher et al Visscher and Bolstreli [1972] and the duo Jullien-Meakin Jullien and Meakin [1987], Jullien et al. [1992a,b, 1996], Jullien and Meakin [2000] who popularized this method.

The general principle of this method is quite simple. As shown in figure 2.1, it consists of placing the particles consecutively on a substrate or a layer of already deposited particles. Each particle first touches the substrate or a deposited particle, then undergoes a relaxation process (single-particle restructuring) along the steepest descent (steepest descent model) until a more stable position according to a stability criterion is reached. As the construction of the packing proceeds layer by layer from the substrate, this deposition model is also known as bottom-to-top restructuring model.

The first step is to release a particle from a random position above the substrate. Upon contact with the first deposited particle, the particle rolls following the steepest descent until a new contact is formed with a second particle. In 2D, two contacts are sufficient to balance a particle if its center of gravity lies between the two contacts. This corresponds to a position of local stable equilibrium. If this criterion is not met, the particle continues to roll and the procedure is iterated until a local stable position is reached. The wall effects can be eliminated by requiring periodic boundaries in the horizontal direction (perpendicular to that of deposition). Figure 2.1 shows a small sample prepared by this method (the grey particles are periodic images of the black ones). In 3D, the same procedure works but a deposited particle needs three contacts to be stabilized.

In this method, the order of deposited particles is generally random and independent of their sizes. In figure 2.1, we see that in a sample prepared by ballistic deposition large voids are created below larger particles. This feature is partially due to the local character of relaxation. In the section 3.1, we will see that this feature can be reduced by employing global relaxation criteria.

An interesting variant of the deposition method, to which we will refer as
“capture angle” (CA) method, consists of making the local relaxation depend on
the direction of the normal at the contact point. Figure 2.2 shows the geometrical
configuration of a falling particle i and an already deposited one j. The contact
normal makes an angle _{ij} with the vertical direction. The CA method is defined
as follows Watson et al. [1997], Bratberg et al. [2002]. If the angle θ_{ij} is below a
critical angle θ_{c} in the range [0,π∕2], the particle sticks irreversibly and does not
move any more. On the contrary, if θ_{ij} > θ_{c}, particle i rolls down until a new
contact is formed with a particle k of contact norman making an angle θ_{ik} with
the vertical. Then, two cases are possible: 1) if |θ_{ik}| < θ_{c} or the new position
reached by i is stable in the sense defined previously, then the particle i sticks;
2) if neither of these criteria is fulfilled, then the particle continues to
relax until a stable position satisfying one of the two preceding criteria is
reached.

Physically, the irreversible capture of particle i for a contact angle below the
critical angle corresponds to a rolling resistance due to a cohesion. The
equilibrium of particle i subjected to his own weight is ensured by a torque M
exerted at the contact by particle j. The particle i is equilibrated is the torque can
balance the weight. Generally, in analogy with sliding friction, the rolling
resistance is characterized by a threshold M_{r} that cannot be exceeded. If this
threshold is proportional to the normal force, the rolling condition is reduced to a
limit angle of the normal beyond which particle i rolls on particle j. Consequently,
the critical angle introduced in the CA method corresponds to a rolling threshold
Bratberg et al. [2002].

The CA model has two particular limits. If θ_{c} = π∕2, the particles stick at the
first contact and thus no relaxation occurs. This limit corresponds to
the model of random sequential adsorption Vold [1959a], Meakin and
Jullien [1985], Lubachevsky et al. [1993], Tang and Liang [1993]. On the other
hand, we have the limit θ_{c} = 0 which corresponds to the steepest descent model
Visscher and Bolstreli [1972], Jullien and Meakin [1987].

Figure 2.2a displays a packing of 500 monodisperse particles prepared by the
CA method with periodic boundary conditions along the horizontal direction. The
critical angle is θ_{c} = 40^{∘}. We observe several particle chains due to irreversible
sticking and dense zones due to geometrical relaxation. Figure 2.2b shows
the same packing obtained by contact dynamics (CD) simulation under
gravity starting from the packing shown in figure 2.2a with an angle of
rolling friction equal to the critical angle θ_{c}. Although all particles in
packing obtained by the CA method satisfy the local stability criterion,
rearrangements occur during CD simulation, showing that the packing is
mechanically unstable. This difference is not only due to inertial effects,
which have been minimized during the CD simulation, but result rather
from the fact that the local stability criteria in the CA method are not
sufficient to ensure the stability of the same particles supporting consecutive
deposited particles. This shows the limit of geometrical methods, but it is
noteworthy that the rearrangements in the packing constructed by the
CA method and the variation of packing fraction are quite small and,
therefore, the geometrical packing is a good starting point for dynamic
simulations.

Figure 2.2a shows the evolution of the packing fraction ρ as a function of
the capture angle θ_{c} for the packings prepared alternatively by the CA
method and CD method. In both cases, ρ declines from 0.82 for θ = 0
to 0.37 for θ = π∕2.. However, the packing fraction is always higher for
the packing prepared by the dynamic method than by the geometrical
method. This difference is a consequence of collective rearrangements which
are absent from geometrical method. Nevertheless, the values of packing
fraction are comparable in the range 20^{∘} < θ_{
c} < 80^{∘}. The coordination
number Z, shown in figure 2.2b exhibits the same trends. We see that Z
decreases from a value slightly above 4 for the purely ballistic deposition
(θ_{c} = 0) to a value slightly above 2 for deposition with irreversible adhesion
(θ_{c} = 90) in which the particle chains largely dominate the texture of the
packing.

We have seen that the ballistic deposition of highly polydisperse particles leads to significant voids around the largest particles. This effect, which reflects the local nature of the relaxation process in this protocol, reduces the homogeneity and packing fraction of the samples. To solve this problem, we need to replace ballistic deposition by a method allowing for a better filling of the space. In this section, we present two methods of this type: an alternative deposition method and a Monte Carlo method.

The ballistic deposition method has a random feature in the choice of particle positions. An alternative way to the deposit the particles is to place the particles in such a way to minimize directly a potential Ψ({(i)}), a scalar function of particle positions (i). By definition, different equilibrium states of the system correspond to different positions of the particles such that the bulk force (i) = -∂Ψ∕∂(i) is equilibrated by the contact forces acting on particle i. Since we do not care about the computation of forces in this geometrical approach, this potential can be any decreasing function of the packing fraction so that its minimum coincides with the maximum packing fraction for each deposited particle.

Since numerical efficiency is at the heart of the geometric method, the potential must be simple and computationally low-cost. For example, assuming that the particles are subjected to a gravitational field, the potential energy of the system is simply given by:

| (1) |

where m(i) and z(i) are the mass and vertical position (parallel to the local gravity direction), respectively, and g is the acceleration of gravity. Obviously, the minimum of Ψ coincides with the lowest position of a deposited particle. Thus, in order to minimize the potential of a particle of size R, one just has to find the lowest stable position (in sense defined previously) depending on R.

Figure 3.1 shows the function Ψ for the same geometrical configuration as in the example of ballistic deposition. Figure 3.1 shows the same particles as in figure 2.1 but deposited with the minimization method. We see that there is now much less void around the big particles. This is not due to a better placement of the big particles but rather to an efficient placement of the smallest ones. Indeed, the function Ψ is all the more sensitive to the details of the the substrate that the size of the particle to be placed is small. The number of possible positions increases and leads to a higher probability of filling, rather than obstructing the pores. This phenomenon is illustrated in figure 3.1 where the function Ψ is evaluated for a big particle 3.1 and a small one 3.1 for the same substrate. Despite the minimization, we see that a large particle can obstruct a half-filled pore. But, minimizing the potential for the small particles rather than local relaxation (shown in black in the figure 3.1) favors pore filling.

Figure 3.1 shows the packing fraction ρ as a function of size span s for the two methods (local relaxation and potential minimization). The span s represents the extension of the PSD. In this example, the PSD is uniform by volume fraction. As expected, higher levels of packing fraction are obtained with the minimization method for all size spans s. Moreover, the increasing difference with s between the two methods indicates that the minimization method is increasingly more efficient as the size span increases. Another important point is that both methods lead to the same trends of packing fraction with size span: ρ first declines, then increases after passing by a minimum.

The periodic boundary condition in the horizontal direction limits the boundary effects but it does not remove the finite size effects. Indeed, in the minimization method the number of possible positions essentially depends on two parameters: the size of the deposited particle and the width L of the periodic domain. As we have seen above, the small particles see more details of the surface and find an optimal pore-filling position. Increasing L leads to an increase of the number of possible potential sites. Therefore, the probability to place a particle correctly without generating too much void around it, increases with L.

Figure 3.1 shows the influence of the width L on the packing fraction ρ for the
two deposition methods and for the same size span s = 0.97 with 10^{5}
particles. For the potential minimization, we see that the packing fraction
increases logarithmically over one decade with L. Independently of the
PSD, the packing fraction increases about 0.02 per decade. With the
ballistic method, the packing fraction is nearly independent of L. This can
be explained by the fact in this method the position of the deposited
particle is guided by the local relaxation process rather than a global
criterion.

A way of achieving dense packings consists of applying a dilatation-relaxation cycle on the particles Metropolis et al. [1953], Barker and Mehta [1992]. The two steps of dilatation and relaxation can be applied in the form of probabilistic geometrical transformations as in the Monte Carlo (MC) method imitating the physical process of densification of a granular packing under vibrations Knight et al. [1995]. In contrast to sequential methods, the dilatation-relaxation method involves periodic rearrangements of the whole packing. Since the particle positions are not frozen, this protocol can lead to higher levels of packing fraction than in the sequential methods where the relaxation is applied to a single particle at a time. The relaxation-dilatation method can be used to achieve dense packings from an initial configuration prepared by the sequential method.

We start with a collection of N_{p} non-overlapping particles inside a given
domain. To this system we attribute a scalar potential Ψ({(i)}), a decreasing
function of the packing fraction. With respect to this configurational potential,
the relaxation process consists of a global evolution of the system towards lower
potential states until a “ground” state with no evolution is attained. This is the
state of mechanical equilibrium. The most direct and natural potential is the
porosity 1 -ρ. But, the local porosity can be hard and computationally expensive
to evaluate for each MC relaxation step. Thus, as in deposition method by
potential minimization, the gravitational potential (1) can be a reasonable
choice. Due to particle displacements, the variations of Ψ are easy to
evaluate and its decrease leads on the average to the increase of packing
fraction.

During the dilation step, the particles should be rearranged such that the gaps between the particles and thus the pore space increases. For example, for particles confined in a box we may simply move the whole packing a distance δ upward. This creates a gap between the bottom of the packing and the base of the box. The relaxation process will fill this gap with a propagation of particle rearrangements upwards into the packing.

The dilatation can also be applied in a homogenous way in the bulk by an affine transformation:

| (2) |

where A is a diagonal matrix with elements (a_{x},a_{y},a_{z}) which are the
three dilatation coefficients in the three space directions, and _{0} is the
origin of dilatation. The origin may coincide with the position of one
particle in the system. For an isotropic transformation a_{x} = a_{y} = a_{z} = a,
two contiguous particles with a distance d between their centers will be
separated by a gap ad after the transformation. Thus, the coefficient a
corresponds to the dilatation amplitude normalized by the particle size. The
transformation can also be unidirectional such that, for example, a_{x} = a_{y} = 0 and
a_{z} = a.

Unlike dilatation, the relaxation stage is sequential. It consists of the application of a MC step to one randomly chosen particle. This particle is moved to a new position ”(i) given by:

| (3) |

where is a random direction, δ is the maximum displacement and ξ is a
random number drawn from a uniform distribution in the range [0, 1].
The position ”(i) of the particle i is rejected if it creates an overlap
with a neighboring particle. Otherwise, it is accepted if the potential
declines, i.e. ΔW = W(”) - W(′) < 0, or accepted with a probability
P(ΔW) = βe^{-βΔW } if ΔW > 0. The parameter 1∕β represents an “activation
temperature”. In practice, when the temperature is small, the probability
P(ΔW) ≃ 0 so that the displacements are basically guided by potential
decrease.

The relaxation process continues by applying the MC step a large number of
times. As a result, the packing fraction increases and tends asymptotically to a
constant value. The relaxation stage stops either when a maximum number of MC
steps N_{MC} is performed or if a precision criterion for the value of the packing
fraction is satisfied. The precision criterion can be the relative variation of the
packing fraction Δρ∕ρ. One cycle, composed of one dilatation step and a
sequence of MC trials (relaxation), leads to compaction by a factor k. Since
the packing fraction increases during a cycle, the dilatation amplitude
for the next cycle must be reduced at least by a factor k^{1∕3} in 3D. In
this way, a large number of cycles can be performed by decreasing the
amplitude in each cycle. The compaction rate depends on the parameters
δ (initial value), β and ϵ. The two parameters δ and β play a similar
disordering role, allowing the particles to “breath” and thus to explore nearby
configurations. The asymptotic packing fraction increases with β and tends to its
highest accessible value in the limit β →∞, i.e. for P(ΔW) = 0 when
ΔW > 0.

The potential used for the simulation of compaction process is often the gravitational potential Ψ of equation (1). As mentioned above, the local porosity can also be used as potential. The local porosity can be evaluated by means of a Voronoi diagram which permits one to attribute a porosity to each particle and thus to measure its variation in each MC trial. It should be noted that the connectivity of the particles cannot be controlled in the MC method. In fact, the particles never touch during relaxation unless a small overlap is allowed between the particles.

The presented method are the most generic geometric methods capable of handling all kinds of size distributions. There can be found other methods in the literature mostly dedicated to particular systems. For example, there are particle insertion models inspired by the parking lot model Renyi [1958]. Another example is the random sequential addition algorithm Widom [1966], Torquato [2002], Torquato and Stillinger [2006], Zhang et al. [1997] that consists of placing consecutively and randomly the particles in the voids between particles avoiding overlap. It should also be remarked that most methods are not appropriate for nonspherical particles.

The geometrical methods can be used as an independent tool for the investigation of the granular media. But in association with dynamic methods, they can be used as a preparation tool to construct genuine starting configurations. The constructed packings are not mechanically stable and dynamic relaxation is required to obtain suitable texture properties (isotropy, connectivity...).

G. C. Barker and Anita Mehta. Vibrated powders: Structure, correlations, and dynamics. Phys. Rev. A, 45(6):3435–, 1992. URL http://link.aps.org/abstract/PRA/v45/p3435.

I. Bratberg, F. Radjai, and A. Hansen. Dynamic rearrangements and packing regimes in randomly deposited two-dimensional granular beds. Phys. Rev. E, 66(3):031303–, 2002. URL http://link.aps.org/abstract/PRE/v66/e031303.

R. Jullien and P. Meakin. Simple three-dimensional models for ballistic deposition with restructuring. EPL (Europhysics Letters), 4(12): 1385–1390, 1987. ISSN 0295-5075.

R. Jullien and P. Meakin. Computer simulations of steepest descent ballistic deposition. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 165(1-3):405–422, 2000. URL http://www.sciencedirect.com/science/article/B6TFR-3YF9X8C-T/2/a8c9ae4c28941312482fef9f076df0a2.

R. Jullien, P. Meakin, and A. Pavlovitch. Three-dimensional model for particle-size segregation by shaking. Phys. Rev. Lett., 69(4):640–, 1992a. URL http://link.aps.org/abstract/PRL/v69/p640.

R. Jullien, P. Meakin, and A. Pavlovitch. Random packings of spheres built with sequential models. J. of Phys. A, 25:4103, 1992b.

R. Jullien, P. Jund, D. Caprion, and D. Quitmann. Computer investigation of long-range correlations and local order in random packings of spheres. Phys. Rev. E, 54(6):6035–, 1996. URL http://link.aps.org/abstract/PRE/v54/p6035.

James B. Knight, Christopher G. Fandrich, Chun Ning Lau, Heinrich M. Jaeger, and Sidney R. Nagel. Density relaxation in a vibrated granular material. Phys. Rev. E, 51(5):3957–, 1995. URL http://link.aps.org/abstract/PRE/v51/p3957.

B. D. Lubachevsky, V. Privman, and S. C. Roy. Morphology of amorphous layers ballistically deposited on a planar substrate. Phys. Rev. E, 47(1):48–, 1993.

P. Meakin and R. Jullien. Structural readjustment effects in cluster-cluster aggregation. Journal de Physique, 46(9):1543–1552, 1985. doi: 10.1051/jphys:019850046090154300. URL http://dx.doi.org/doi/10.1051/jphys:019850046090154300.

Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller. Equation of state calculations by fast computing machines. J. Chem. Phys., 21(6):1087–1092, 1953.

A. Renyi. On a one-dimensional problem concerning random space-filling. Publications of Mathematical Institute of Hungarian Academy of Sciences, 3:109–127, 1958.

Chao Tang and Shoudan Liang. Patterns and scaling properties in a ballistic deposition model. Phys. Rev. Lett., 71(17):2769–, 1993.

S. Torquato. Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer, 2002.

S. Torquato and F. H. Stillinger. Exactly solvable disordered sphere-packing model in arbitrary-dimensional euclidean spaces. Phys. Rev. E, 73(3):031106–8, 2006.

W. M. Visscher and M. Bolstreli. Random packing of equal and unequal spheres in two and three dimensions. Nature, 239(5374):504–507, 1972.

M. J. Vold. A numerical approach to the problem of sediment volume. Journal of Colloid Science, 14(2):168–174, 1959a.

M. J. Vold. Sediment volume and structure in dispersions of anisometric particles. J. Phys. Chem., 63(10):1608–1612, 1959b.

M. J. Vold. The sediment volume in dilute dispersions of spherical particles. J. Phys. Chem., 64(11):1616–1619, 1960.

P.K. Watson, H. Mizes, A. Castellanos, and A. Pérez. The packing of fine cohesive powders. In R.P. Behringer and J.T. Jenkins, editors, Powders and Grains 1997, pages 109–112, Rotterdam, 1997. Balkema.

Widom. Random sequential addition of hard spheres to a volume. J. Chem. Phys., 44:3888–3894, 1966.

Z. P. Zhang, A. B. Yu, and J. A. Dodds. Analysis of the pore characteristics of mixtures of disks. Journal of Colloid and Interface Science, 195(1):8–18, 1997.