Assembling methods

Charles Voivret, Farhang Radjai and Jean-Yves Delenne

1 Introduction

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.

2 Ballistic deposition

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

2.1 Random deposition with relaxation

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,b1960], 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,b1996], 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.


Figure 1: (a) Ballistic deposition: first contact followed by steepest descent. (b) Small-scale periodic sample.

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.

2.2 Angle of capture

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 ⃗nij 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.


Figure 2: Geometrical configuration of a falling particle i in contact with an already deposited particle j.

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 Mr 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 3: (a) Packing composed of 500 monodisperse particles prepared by the capture angle method; (b) The same packing after dynamic relaxation by the contact dynamics method.

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.


Figure 4: (a) Evolution of the packing fraction ρ as a function of the capture angle for the CA and CD methods; (b) Evolution of the coordination number Z for the same methods.

3 Assembling method with configurational potential

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.

3.1 Potential minimization 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 Ψ({⃗r(i)}), a scalar function of particle positions ⃗r(i). By definition, different equilibrium states of the system correspond to different positions of the particles such that the bulk force f⃗(i) = -Ψ∕∂⃗r(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.


Figure 5: (a) Potential minimization method: placing directly a particle at the lowest possible position given its size, (b) Small periodic sample prepared by minimizing a potential.

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:

Ψ =  g    m (i)z (i)

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 6: Effect of the size of the deposited particle on the “precision” of the function Ψ : (a) for a large particle; (b) pour a small particle; the black particles represent the two local equilibrium positions.

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


Figure 7: (a) Packing fraction ρ as a function of size span s for a uniform size distribution by volume fraction and for the two deposition methods. (b) Packing fraction ρ as a function of the period L rescaled by the maximum diameter dmax for two different particle size distributions but with the same size span.

3.2 Monte-Carlo Method

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 Np non-overlapping particles inside a given domain. To this system we attribute a scalar potential Ψ({⃗r(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:

⃗r = ⃗r + A  (⃗r - ⃗r0)

where A is a diagonal matrix with elements (ax,ay,az) which are the three dilatation coefficients in the three space directions, and ⃗r0 is the origin of dilatation. The origin may coincide with the position of one particle in the system. For an isotropic transformation ax = ay = az = 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, ax = ay = 0 and az = 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 ⃗r”(i) given by:

 ′′      ′
⃗r (i) = ⃗r (i) + δ(1 - 2ξ)⃗n

where ⃗n 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 ⃗r”(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(⃗r”) - W(⃗r) < 0, or accepted with a probability PW) = βe-βΔW if ΔW > 0. The parameter 1∕β represents an “activation temperature”. In practice, when the temperature is small, the probability PW) 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 NMC 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 k13 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 PW) = 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.

4 Conclusion

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

   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

   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

   R. Jullien, P. Meakin, and A. Pavlovitch. Three-dimensional model for particle-size segregation by shaking. Phys. Rev. Lett., 69(4):640–, 1992a. URL

   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

   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

   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

   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.