This chapter describes numerical and simulation techniques that are used in computational studies of dense granular materials. In particular,we will focus on computational methods to generate static packings of model granular materials and to measure their response to compression and simple shear deformation modes.
This chapter describes numerical and simulation techniques that are used in computational studies of dense granular materials. In particular,we will focus on computational methods to generate static packings of model granular materials and to measure their response to compression and simple shear deformation modes.
There are two main classes of computational methods for generating static packings of granular materials: hard- and soft-particle methods. For the first, particles interact via a hard-core potential for which the pair interaction energy is infinite if the particles overlap, whereas it is zero otherwise. The packing-generation protocol can involve a number of different types of boundary conditions and driving mechanisms for obtaining a broad ensemble of static packings including isotropic compression, simple, and pure shear. The hard-particle methods include various geometrical [7,15,40] and Monte Carlo  techniques and molecular dynamics (MD) of hard particles  coupled with a series of compressions of the system (or increases in particle sizes), for example, the Lubachevsky– Stillinger algorithm . Soft-particle methods employ energy minimization techniques or dissipative MD simulations  to relax the forces between overlapping particles after each successive compression of the system for dilute initial configurations (or each successive decompression of the system for dense, overcompressed initial configurations). Many of the geometrical and Monte Carlo packing-generation methods employ “single-particle” moves, displacing one particle at a time independently, while energy minimization and MD techniques employ “collective” moves involving the cooperative motion of many particles. As particular examples, we describe packing-generation algorithms for (1) fric-tionless disks and spheres, (2) frictional spherical particles, (3) frictionless non-spherical composite particles such as dimers, and elongated particles including (4) ellipsoids, (5) spherocylinders, and (6) superquadrics.
The basic Monte Carlo method for generating hard-particle packings combines single particle movements with system compressions, as illustrated in Figure 3.1 for a bidisperse collection of frictionless disks. In Figure 3.1 half of the disks are large and half are small, with diameter ratio σ l /σ s = 1.4(σ l,s are the diameters of the large/small disks). The system is initialized with N randomly placed points and the disks grow uniformly until the closest pair is at contact (panel a). At this point, each disk is given a random displacement in the x-and y-directions (with the root-mean-square displacement for a given particle set by the average distance between its three closest neighbors). The move is accepted if it does not give rise to particle overlaps as shown in panel (b). After N attempts, the particle sizes are again increased uniformly until the nearest pair is at contact (panel c). This process of moving and growing the particles is repeated for N I iterations to yield a single static packing. Further, N C configurations are generated by seeding the Monte Carlo method with N C independent sets of N initial points.
Figure 3.1 Illustration of the Monte Carlo packing-generation method for a system of N = 6 bidisperse frictionless disks (half small and half large with diameter ratio 1.4) in 2D. (a) The x-and y-coordinates for N = 6 random points are first generated in a square cell with periodic boundary conditions and the particles are grown uniformly until the closest pair of disks is in contact. (b) An attempt is made to move each particle randomly from the original (dashed outline) to the new position (shaded disk), and the move is accepted if it does not give rise to particle overlap. (c) The disks are expanded uniformly from the original (dashed outline) to the new size (shaded disk) until the two closest disks touch. This process is repeated for N I iterations to obtain a single static packing.
The results for N = 6 bidisperse hard frictionless disks in a unit square with periodic boundary conditions are shown in Figure 3.2. For this system, the packing fraction converges to its final value within N I ≈ 103 iterations. Figure 3.2b compares the distribution of final packing fractions from the Monte Carlo packing-generation method (for N C = 105 static packings) to that for the soft-particle MD method. (See Section 3.2) The Monte Carlo method generates each of the mechanically stable (MS) packings obtained from the soft-particle MD method, but with different probabilities, as well as many additional MS configurations (Figure 3.2c) that are not obtained from the soft-particle MD method.
Figure 3.2 (a) The deviation in the packing fraction ф from its final value ф f as a function of the number of iterations N I of the Monte Carlo packing-generation method for three independent trials (different symbols) using bidisperse frictionless disks as in Figure 3.1. The dashed line indicates an average over N C =103 independent trials. (b) The probability distribution of final packing fractions after N I = 105 iterations obtained from N C = 104 (filled circles) and 105 (rectangles) independent trials. (c) The probability distribution of final packing fractions (rectangles) after N I = 105 iterations and N C = 105 trials. The open circles give the probabilities for obtaining the 20 distinct mechanically stable packings of frictionless bidisperse disks for N = 6 using soft-particle molecular dynamics with the packing fractions indicated by the solid vertical lines.
To characterize the structural properties of the packings generated via the Monte Carlo method, we calculated the number of interparticle contacts (defined as T ij ≤ δ c ,where δ c is a small cutoff, typically such that T ij /σ ij −1 ≤ 10−3,where σ ij = (σ i + σ j )/2 and σ i is the diameter of disk i), the N × N distance matrix for each configuration
whereis the location of particle i, and its second invariant
The cumulative distribution of interparticle separationsscaled by N for packings of bidisperse frictionless hard disks generated using the “basic” Monte Carlo method is shown in Figure 3.3a and b for several values of the total number of iterations N I and two system sizes N = 6 and 12. As N I increases, a plateau develops in the cumulative distribution, which indicates a bimodal distribution of separations r ij . The probability distribution of contact number differences over all packings obtained from the Monte Carlo method is shown in Figure 3.3c. is the isostatic number of contacts for frictionless disks, where the number of contacts matches the number of degrees of freedom and N r is the number of rattler particles with two or fewer interparticle contacts. For N = 6, the Monte Carlo method generates similar probabilities for MS packings and first-order saddle points , whereas higher-order saddle points have frequencies of 1% or less. As N increases, the probability for obtaining first- and second-order saddles increases above that for MS packings. Thus, at large N, the basic Monte Carlo method becomes inefficient at generating MS packings.
A scatter plot of the second invariant of the distance matrix q 2 for the N c = 105 packings obtained from the “basic” Monte Carlo packing-generation method (after N I =105 iterations) versus the packing fraction is shown in Figure 3.4a. The first-order saddle packings form 1D lines (or geometrical families [12,13]) in the space of q 2 versus ф, whereas MS packings occur as discrete points. Two first-order saddle packings that belong to the same geometrical family with the same network of interparticle contacts (but different q2) are shown in Figure 3.4b and c. As q 2 varies from −3.1to −3 along this geometrical family (cf. Figure 3.4a), particle 1 moves toward 2. This geometrical family terminates in a MS packing when particle 1 comes into contact with 2. The packings that are missing two contacts (second-order saddle packings) form more diffuse lines in the q 2−ф plane. Improved sampling of the higher-order saddle packings is necessary to determine to what extent they form lines or fill in 2D areas.
Figure 3.3 Cumulative distribution C ¯ ( r i j / σ i j − 1 ) of interparticle separations r ij /σ ij − 1 normalized by N for packings generated using the “basic” Monte Carlo method for N I = 104 (triangles), 105 (diamonds), 106 (squares), 107 (circles), and 108 (stars) and two system sizes (a) N = 6 and (b) 12. The dashed vertical line indicates the interparticle separation δ c below which we define two particles to be in contact. (c) The distribution of contact numbers N c iso − N c using the criterion that contacting particles satisfy r ij /σ ij −1 ≤ δ c = 10−3 for N = 6 (circles) and 12 (squares).
Over the past 10 years, the O’Hern group has developed a suite of discrete element simulations to generate static granular packings using several packing-generation protocols including slow and fast isotropic compression [10,11,33], gravitational deposition , simple shear , and vibration . O’Hern and coworkers have studied packings of frictionless [10–13] and frictional  disks, particles with anisotropic shapes including dimers , ellipses [20,31,35], ellipsoids , and bumpy particles , and particles with purely repulsive and adhesive contact interactions .
Figure 3.4 (a) The second invariant q 2 = ((TrD)2 − (TrD)2)/2 of the distance matrix D plotted versus the packing fraction φ for N c = 105 bidisperse frictionless hard disk packings with N = 6 generated via the “basic” Monte Carlo method (after N I = 105 iterations) with one (exes) or two (triangles) missing contacts N c iso − N c = 1 or 2 relative to the isostatic value for frictionless disks. The MS packings with N c = N c iso are indicated by the filled circles. Two configurations shown in panels (b) and (c) with N c iso − N c = 1 that belong to a single geometrical family are labeled (b) and (c) in panel (a).
We describe two numerical procedures for creating isotropically compressed MS packings of frictionless particles at jamming onset: (1) a Lubachevsky–Stillinger (LS) [18,19] method for hard particles and (2) a dissipative MD method for soft particles that interact via purely repulsive contact forces .
For the LS method, N particles with mass m are placed with random positions in a square box with unit length (L = 1) in 2D (or cube with unit length in 3D) with no particle overlaps and a Gaussian velocity distribution at temperature T . The particle positions are updated using energy- and momentum-conserving event-driven dynamics for hard spherical particles for a time period t, after which all particle sizes are increased by, where is the smallest separation between all pairs of particles. The particles are moved forward in time at constant velocity until the next interparticle collision. The velocities of the two colliding particles are updated and all of the particles are again moved forward in time at constant velocity. The process is repeated for a given number of collisions [18,19].
Figure 3.5a shows the collision frequency f (in units ofversus packing fraction φ for the LS method for a 2D bidisperse system with N = 6 disks. The collision frequency diverges as f ~ (ф J − ф)−1 near a MS packing with ф J ; the compressions are stopped when the collision frequency exceeds a large value, f =106
For the MD soft-particle packing-generation method, N particles are again given random, nonoverlapping positions in the simulation cell, this time with zero velocities. Collisions are now not assumed instantaneous, but rather the particles interact via the purely repulsive linear spring potential (with ∝ =2):
Particle sizes are again increased uniformly by Δф, and the total potential energy V = ∑ ij V(r ij ) is minimized using overdamped MD simulations. If the system is below jamming onset, the potential energy can be minimized to zero. In this case, we terminate the energy minimization when V/Nϵ <Vtol = 10−16.Incon-trast, when the system is above jamming onset, the minimized potential energy is nonzero, and we terminate the minimization when (V t − V t+1)/(V t + V t+1) <V tol = 10−16,whereV t is the total potential energy at time t, as shown in Figure 3.5b. Each time the system switches from below to above jamming onset or vice versa, the packing fraction increment is decreased by a factor of 2. We terminate the MD packing-generation procedure when the minimized total potential energy obeys V tol <V/Nϵ < 2V tol. The MS packings obtained from this method are accurate to 10−8 in packing fraction and particle positions.
Using the LS and MD packing-generation methods, we performed more than 106 independent trials for each system size to enumerate nearly all MS static packings in small 2D bidisperse and 3D monodisperse systems. Figure 3.6 shows a 2D projection of the paths in configuration space for the LS and MD packing-generation methods with two different initial conditions for 2D bidisperse systems with N = 6. We plot, where particle 1 is fixed at the origin. For the calculation of , the configurations from the LS and MD methods were aligned so that the particle positions and labels were identical in the final MS packings. For the MD and LS methods, the packing-generation process was initialized with the configurations in shapshots (a) and (∝), respectively. For the MD method, the system populates configurations (a), (b), and (c), and terminates at a particular MS packing (center of Figure 3.6). In contrast, for the LS method, the system visits different configurations (∝), (β), and (γ), but ends at the same MS packing obtained from the MD method.
Figure 3.5 (a) Collision frequency f for the LS packing-generation procedure as a function of packing fraction ф for two compression rates (triangles: higher, diamonds: lower) for a 2D bidisperse system with N = 6. The inset shows f versus ф J − ф on a log10 − log10 scale for the data in the main panel. The dashed line has slope −1. (b) The total potential energy V/Nϵ as a function of time during energy minimization for the MD packing-generation procedure at three packing fractions, two above jamming onset (squares and triangles) and one below (circles) for 2D bidisperse systems with N =6.
Figure 3.6 Two-dimensional projection of the paths in configuration space for the LS (circles) and MD (squares) packing-generation methods with two different initial conditions for 2D bidisperse systems with N = 6. For these initial conditions, both methods generate the same MS packing. For the calculation of r → p = 1 5 ∑ i = 2 6 [ ( x 1 − x 1 ) x ^ + ( y i − y 1 ) y ^ ] , where particle 1 was fixed at the origin (in the lower left corner), the configurations from the LS and MD methods were aligned so that the particle positions and labels were identical in the final MS packing.
Even though the LS and MD packing-generation methods yield the same ensemble of MS packings at jamming onset, the probabilities with which the MS packings occur can vary strongly with the method. Figure 3.7 shows the number of times that each of the 20 distinct MS packings for N = 6 bidisperse disks was obtained over 106 trials. The MS packing probabilities increase with packing fraction, although there are fluctuations. The LS and MD methods give similar results for the MS packing probabilities for highly probable packings, but not the low-probability packings. Thermal fluctuations cause the low-probability packings to become more rare for the LS method compared to the MD method with overdamped dynamics. The MS packing probabilities for the LS and MD methods will approach each other in the fast compression rate limit for LS. In contrast, only the densest packing will be obtained using the LS method in the slow compression rate limit.
Figure 3.7 Number of times out of N t = 106 trials that each of the 20 distinct MS packings at ф J is obtained using the LS (squares) and MD (circles) packing-generation methods for 2D bidisperse systems with N =6.
To provide a specific example of soft-particle packing-generation methods, we describe discrete element simulations of bidisperse mixtures of frictional disks subjected to isotropic compression in 2D, implementing the Cundall–Strack model for static friction . The position of the center of mass of each disk i (with mass m) obeys the equation of motion
which is the sum of an elastic contact force, a viscous damping (relative to a stationary background) force, and a contact friction force. (is the separation vector between the centers of disks i and j, and the sum over j indicates a sum over particles in contact with disk i). The elastic force ,and V is the purely repulsive elastic energy (Equation 3.3, here with ∝ = 2) between overlapping disks i and j. The viscous damping force is , with a damping coefficient chosen to be in the overdamped limit. The force on disk i arising from static friction with contacting disks j is with pair tangential forces at each contact, given by
The relative tangential displacementbetween disks i and j is obtained by solving
with, and . During the lifetime of each contact, the Coulomb sliding condition is enforced, where μ is the static friction coefficient. When a contact between disks i and j breaks, is reset to zero. This model for static friction has been implemented in a number of discrete element simulations of static and flowing granular media [8,38].
The rotational degree of freedom for each disk θ i obeys
where the torque on disk i
the moment of inertia, and the rotational damping coefficient is in the overdamped limit.
The soft-particle MD packing-generation algorithm for frictional particles proceeds as follows. First, disks are placed at random positions in a square cell of length L at low packing fraction ф = 0.01 with zero translational and rotational velocities. The disks are enlarged by an initial packing fraction increment Δф =10−3 and, after each compression, Equations 3.4 and 3.7 are integrated until the maximum particle kinetic energy, net force, and net torque are sufficiently small (here K max = max i K i /ϵ < 10−15,, and τmax =max i |τ i |/ϵ < 10−11,where ). These conditions ensure that all particles in the system are in force and torque balance after each compression.
A single MD relaxation during the packing-generation procedure is shown in Figure 3.8 for N = 6 bidisperse systems of frictional disks with μ =0.05 at packing fractions (a) below and (b) above the onset of jamming. For both (a) and (b), the maximum kinetic energy and torque on a single disk reach their terminal values before the maximum force on a disk does. In panel (a), the final total potential energy per particle, whereas in panel (b). To generate “just-touching” jammed packings, we successively compress relaxed “unjammed” packings, or decompress relaxed “jammed” packings, with successively smaller packing fraction increments Δф each time the process switches from compression to decompression or vice versa. The packing-generation algorithm is stopped when the relaxed total potential energy per particle satisfies .
Figure 3.8 Maximum kinetic energy K max (circles), force magnitude F max (squares), and torque τmax (diamonds) on a given disk, and total potential energy per particle V/ϵN (triangles) during a single MD relaxation during the packing-generation procedure for a packing fraction (a) below and (b) above the onset of jamming for N = 6 bidisperse frictional disks with μ =0.05. The dashed (dotted) horizontal line indicates F max, τmax = 10−11 (K max, V ¯ = 10 − 16 ).
All of the packing-generation methods discussed to this point in the chapter employed periodic boundary conditions in cubic cells with side length L. An illustration of periodic boundary conditions for a 2D MS packing of bidisperse disks with N = 6 is provided in Figure 3.9. Eight image cells surround the central cell to specify the interparticle interactions near the boundaries. (In 3D, 26 image cells surround the central cell.) We see that disk 5 interacts with disks 3, 4, and 6 in the central cell. In contrast, disk 6 interacts with disk 5 in the central cell, 2′ is the lower image, 4″ in the right image, andin the lower right image. (Note that if the diameter of one of the particles is larger than half of the box size, σ >L/2, it can interact with more than one image of the same particle.) To maintain continuous evolution of the particle trajectories, the particles may move out of the central cell during the packing-generation process, but this does not affect the interparticle separations and forces.
Figure 3.9 Illustration of periodic boundary conditions for a mechanically stable packing of N = 6 bidisperse disks in a square cell. The disks in the central cell are labeled 1 to 6. The disks are replicated in the eight adjacent image cells. The disks in the bottom image cell are labeled 1′ to 6′. Note that disk 5 in the central cell is in contact with disks 3, 4, and 6 in the same cell, while disk 6 in the central cell is in contact with disk 5 in the central cell, 2′ in the bottom image, 4″ in the right image cell, and 1 ‴ in the bottom right image cell.
These systems can also be used to study the linear (i.e., the shear modulus) and nonlinear response (e.g., shear-induced jamming) of static packings to applied simple shear strain. To implement simple shear, an affine shear is applied to all particles i in the central cell,
Figure 3.10 Illustration of the implementation of simple shear in Lees-Edwards boundary conditions for a packing of bidisperse disks. The particles in the main cell are given the affine deformation in Equation 3.9. The image cells above and below in the central cell in the shear gradient direction are also shifted by ±LΔγ in the shear flow direction.
Shear-periodic boundary conditions are illustrated in Figure 3.10 for 2D systems. The top/bottom row of image cells is shifted to the right/left by ΔγL. As in the quiescent case, we identify the particles in the main or image cells that contact each of the particles in the main cell and, following the application of shear strain, relax the system using energy minimization at fixed locations of the image cells. This is followed by the application of successive shear strain increments and energy minimization to a prescribed total shear strain γ.
Simple shear deformation can also be implemented by shifting top and bottom rigid boundaries by ±ΔγL in the shear flow direction and applying periodic boundary conditions in the directions perpendicular to the shear gradient direction with or without affine distortions in the bulk, as shown in Figure 3.11. After each movement of the boundaries, the system can be relaxed using energy minimization. Boundary-driven shear allows one to model the interactions between the sheared material and the container and set up strongly nonlinear velocity profiles .
In this section, we focus on nonspherical composite particles in 2D formed from collections of disks that are rigidly connected together, for example, dimer- and trimer-shaped particles [34,35]. Packings of these anisotropic particles display interesting structural and mechanical properties, for example, shear strengthening and large yield strength . Composite particles with circular asperities have also been employed to model friction in granular packings  and flows .
Figure 3.11 Illustration of the implementation of simple shear in 2D systems with fixed walls (dark shading) in the shear gradient direction and periodic boundary conditions in the shear flow direction. In this example, only the boundaries were moved with no affine bulk (light shading) deformation, followed by energy minimization at fixed locations of the boundaries. The particle positions before and after a shear strain step (center) and before and after the energy minimization (right) are shown.
Figure 3.12a through c shows static packings composed of bidisperse (50–50 by number) dimer- and trimer-shaped particles in 2D. The schematic in Figure 3.12d defines the geometries. For dimers, the geometry is defined by the aspect ratio ∝ D = L D /σ. Trimers are specified by three shape parameters, sin(θ/2) = βT/∝ T , β T = h T /σ − 1, and ∝ T = (L T /σ − 1)/2. To frustrate positional and orientational ordering, we consider bidisperse mixtures of composite particles that satisfy(for dimers) and (for trimers). Figure 3.13 shows the distribution of packing fractions at jamming onset for systems composed of N = 256 bidisperse, frictionless disks as well as elongated shapes, including dimers, trimers, and ellipses at roughly the same aspect ratio. The location of the most probable packing fraction at jamming onset increases in the following order for this particular set of shapes and aspect ratios: disks, linear trimers, nonlinear trimers, dimers, and ellipses. The widths of the packing fraction distributions are ~1% − 2% for all shapes at N = 256 and decrease as a power law with increasing system size.
Figure 3.12 Mechanically stable static packings composed of 50–50 bidisperse mixtures (by number) of N = 256 frictionless (a) dimers with aspect ratio ∝ D = L D /σ =1.5and size ratio α D l / α D s = 1.4 , trimers with (b) θ = 180°, ∝ T = 0.5, and β T = 0.0, and (c) θ = 90°, ∝ T = 0.5, and β T = 0.5, both with size ratio α T l / α T s = β T l / β T s = 1.4 . (d) Schematic that defines the shape parameters for rigid dimer- and trimer-shaped particles.
Figure 3.13 Distribution of the packing fractions P(ф J ) at jamming onset for 50–50 bidisperse mixtures by number of N = 256 frictionless disks with diameter ratio σ l /σ s = 1.4 (circles), trimers with θ = 180°, ∝ T = 0.5, and β T = 0.0 (squares) and θ = 90°, ∝ T = 0.5, and β T = 0.5 (downward triangles) both with size ratio α T l / α T s = β T l / β T s = 1.4 , dimers with ∝ D = L D /σ =1.5and α D l / α D s = 1.4 (diamonds), and ellipses with ratio of the major to minor axes ∝ =1.5and α l /∝ s = 1.4 (upward triangles).
Identifying physically reasonable interaction laws between dimer-shaped particles is not trivial. When dimer packings are generated at jamming onset with only infinitesimal overlaps (c.f. Figure 3.14a), overlaps between pairs of lobes on different dimers are distinct. (See Figure 3.15b). Computational studies have shown that, at the onset of jamming, dimer packings are typically isostatic with z = z iso =2d f contacts per particle, where d f = 3 in 2D . However, when dimer packings are sufficiently overcompressed (Figure 3.14b), the area of overlap between one pair of lobes merges with that between another pair of lobes (Figure 3.15c and d). Thus, two different contact numbers may be defined: a “lobe-lobe” number z ll that treats all lobe-lobe interactions as distinct, and a single contact number z s in which all merged overlaps are considered as a single contact. These need not be identical, and the difference between the two reveals information about the packing. For the dimer packing Figure 3.14b, the lobe-lobe contact is z ll = 13.5, much greater than the isostatic number z iso = 6, whereas Zs = 5.67 <Ziso.
Figure 3.14 Static packing of N = 48 bidisperse dimers with aspect ratio ∝ D =1.1 and compression (a) Δф = 10−4 and (b) 0.2. The contact number is z = Z iso = 2(3N − 1)/N ≈ 5.96 in (a) and z ll = 13.58 >Ziso (counting all lobe-lobe interactions) or z s = 5.67 <Z iso (merged overlap areas counted as a single contact) in (b). (c) Contact number z ll (dashed line) and z s (solid line) versus overcompression Δф averaged over ≈50 configurations. The dotted line indicates z = z iso =6.
Figure 3.14c shows the two definitions of the contact number z ll and z s as a function of compression Δф for N = 48 bidisperse dimer packings with aspect ratio ∝D =1.1. For low compressions (Δ∝ < 10−3) z ll = z s = z iso.ForΔф > 10−3, z ll begins to increase, reaching ≈ 8 near Δф = 3 × 10−2 when all lobe-lobe contacts have formed between dimers initially in contact at Δф =0. z ll continues to increase to 24 when all lobe-lobe contacts have formed between Voronoi neighbors (defined at Δф = 0). In contrast, z s decreases below Z iso for Δф > 10−3 as some of the lobe-lobe interactions merge into extended contacts. z s returns to the isostatic value 6 at large Δφ when all Voronoi neighbors (defined at Δф = 0) are in contact.
Figure 3.15 (a) Schematic of the possible repulsive interactions between dimer lobes. F → i k i , j k j is the pair force on lobe k i on dimer i arising from lobe k j on dimer j.(b) Two overlapping dimers (i and j) with aspect ratio ∝ D =1.5. The area of overlap between lobes i1 and j1 (dark gray) is separate from the area of overlap between lobes i2 and j2 (light gray). (c) and (d) Two overlapping dimers at ∝ D = 1.5 and 1.1, respectively, for which the areas of overlap between the lobes have merged together (dark gray).
The multiple counting of lobe interactions for compressed dimer-shaped particles gives rise to anomalous behavior for the structural and mechanical properties, which do not approach those for disks in the ∝ D → 1 limit. We now describe a micromechanical model based on the area of overlap between compressed dimer-shaped particles that result in structural and vibrational properties that converge smoothly to those in the ∝ D → 1 limit at all values of the compression Δф.
As shown in Figure 3.16, the total overlap area between two dimers with the same sized constituent disks can be expressed in terms of the area of overlap between two (A 2a , A 2b , A 2c, A 2d ), three (A 3a , A 3b , A 3c , A 3d ), and four lobes (A4)
Figure 3.16 The total ove lap a ea A T fo the pai of dime s in Figu e 3.15d can be decomposed into the areas of overlap between two (A 2a , A 2b , A 2c ,and A 2d ), three (A 3a , A 3b , ,and A 3d ), and four lobes (A4).
Each of these areas (formed from two, three, and four lobe intersections) can be further decomposed as the sum of n circular segments (1 per pair of overlapping lobes) with areas
where, ,and are the locations of the points of intersections of the lobes. (See Figure 3.17.)
For disks that are not the same size (diameters σ i and σ j and center-to-center separation r ij ), the area of overlap A ij (r ij ) is
Figure 3.17 The area of intersection of three dimer lobes, 1, 2, and 3, can be decomposed into a triangle (dark gray) and three circular segments (light gray). (The area of intersection between two and four dimer lobes can be decomposed into a triangle and two and four circular segments, respectively.) The locations of the lobe intersection points are r → 12 , r → 23 and r → 31 , and the lengths of the lines connecting these points are a 1, a 2,and a 3.
and σ ij = (σ i + σ j )/2. To lowest order in r ij /σ ij , this is.
Given a particular amount of overlap A between a given number of dimer lobes, for example, A2a between two dimer lobes, we calculate their effective separation(assuming the lobes are individual disks) by inverting Equation 3.13, that is, . We assume that this effective separation between overlapping lobes gives rise to a purely repulsive linear spring potential, (Equation 3.3), and do this for each of the nine area terms in Equation 3.10 to obtain the total potential energy between dimers i and j:
Keeping only leading order terms in the area of overlap reduces this to
The potential energy in Equation 3.15 converges to the repulsive linear spring interaction potential used in Refs. [34,35] to generate static dimer packings in the Δф → 0 limit for all ∝ D . Equation 3.15 allows the numerical calculation of the total potential energy V = ∑ i>j V ij , force on dimer i,, torque on dimer i, T i = −dV/dθ i ,where θ i is the angle that dimer i makes with the horizontal axis, and the dynamical matrix elements ,where ξ i = x i , y i ,and σ i θ i give the center of mass location and orientation of dimer i.
Figure 3.18e shows the effective force law (force versus displacement) between two parallel dimers with aspect ratio ∝ D =1.3 undergoing compression for the micromechanical model in which (1) lobe interactions are multiply counted or (2) the interaction potential is given by Equation 3.15. The two force laws are the same as long as overlaps between lobes have not merged, or δm/σ < 0.021 for the configuration in Figure 3.18a. Beyond δm, the two force laws differ. The force law based on the total area of overlap converges to linear behavior F ~ δ more quickly than the one that multiply counts lobe interactions, for example, it is not sensitive to the formation of the fourth lobe contact at δ/σ = δ4/σ =0.075. In future studies, these results can be compared to finite element analyses of linear elastic particles with complex shapes.
Figure 3.19a shows the two definitions of the contact number, z ll and z s ,for static packings of N = 48 bidisperse dimers versus Δф over a wide range of aspect ratios. Similar to Figure 3.14c, z ll increases with Δф and z s initially decreases. For Δф ≲ 0.03, changes in contact number are mainly due to the merging of areas of overlap, while for Δф > 0.03, changes in contact number are due to additional contacts from neighbors that were not in contact at Δф = 0. The packings may be divided into three regimes, defined by the degree of multiple counting of lobe contacts: z ll = z s = z iso (no multiple counting), z s <z ll < 2z s (some lobes multiply counted, and z ll >z s (all lobes multiply counted). In Figure 3.19b, we show the different regimes for z ll and z s as a function of Δф and ∝ D . The regime z ll = z s is bounded by (∝D − 1) > Δф1/2. In contrast, the z ll > 2z s regime occurs when (∝D − 1) ≲ Δф. The z s ≥ z ll < 2z s regime occurs for Δф < (∝ D − 1) < Δф1/2.
Figure 3.20 shows the vibrational density of states D(ω) (in the harmonic approximation) for static packings of dimer-shaped particles with and without multiple counting of the lobe interactions as a function of Δф and ∝ D . With multiple counting, the maximum frequency in D(ω) can increase significantly above the disk value to ωmax ≈ 4 (even for nearly circular particles with ∝ D ≈ 1) as lobe-lobe interactions are added. Without multiple counting, the maximum frequency for small aspect ratio dimer packings remains near the disk value ωmax ≈ 2 with increasing Δф until next-nearest neighbors are made, after which it increases to ωmax ≈ 2.5.
Figure 3.18 Snapshots during the compression of two parallel dimers i and j with aspect ratio ∝ D = 1.3 (a) at δ = 0.009σ, two distinct contacts between lobe pairs i 1 and j1 and i 2 and j 2 exist and at (b) δ = 0.017σ, three distinct contacts between lobes i 1 and j 1, i 2 and j 2,andi 1 and j 2 exist. (c) A single extended area of overlap is found at δ =0.045σ. (d) For δ =0.098σ, the extended overlap area includes interactions between two, three, and four lobes. (e) Magnitude of the inter-dimer normal force F versus displacement δ/σ when compressing dimers i and j. The vertical solid lines correspond to the configurations in (a–d). The dashed vertical lines (from small to large δ/σ) correspond to the displacement at which the third interaction between dimer lobes is formed, the contact areas between different lobe pairs merge, and the fourth contact between lobes is formed.
Figure 3.21 shows the vibrational density of states D(ω) for static packings of bidisperse dimers as a function of compression Δφ on a logarithmic scale to focus on the low-frequency behavior. D(ω) contains a strong peak of mainly rotational modes at low frequency. With multiple counting of the lobe interactions in the total potential energy, the location of the low-frequency peak ωmin does not depend on Δф. However, with the interaction potential in Equation 3.15, for the regime in Δф where the overlaps between lobes have merged, the location of the low-frequency peak scales as ωmin ~ (Δф)1/4. Figure 3.22 shows that, with multiple-counting of the lobe interactions, the low-frequency peak ωmin scales linearly with aspect ratio ∝ D − 1. In contrast, implementing the interaction potential in Equation 3.15 yields a different scaling behavior, ωmin ~ (∝ D − 1)1/2, after the overlap areas merge.
Figure 3.19 (a) Contact numbers z ll (gray) and z s (black) versus Δф for static packings of bidisperse dimers with aspect ratio ∝ D =1.0 (solid line), 1.01 (dashed line), 1.03 (dash-dash-dotted line), 1.1 (dash-dotted line), 1.3 (dot-dot-dashed line), and 1.5 (dotted). (b) The Δф and ∝ D − 1 parameter space decomposed into regions where z ii = z s (circles), z s <z ll ≤ 2z s (squares), and z ll > 2z s (triangles). The solid (dashed) line has slope 1/2 (1).
Figure 3.20 The vibrational density of states D(ω) (in the harmonic approximation) for static packings of N = 48 bidisperse dimers at aspect ratio ∝ D = 1 (solid line), 1.01 (dashed line), 1.1 (dot-dashed line), and 1.3 (dot-dot-dashed line) at compression Δф =0.03 (a) with and (b) without multiple counting of the interactions between lobes. The vibrational density of states D(ω) (in the harmonic approximation) for static packings of N = 48 bidisperse dimers at aspect ratio ∝ D = 1 (solid line), 1.01 (dashed line), 1.1 (dot-dashed line), and 1.3 (dot-dot-dashed line) at compression Δф = 0.03 (c) Largest frequency ωmax at which D(ωmax) = 0.1versus Δф with (gray) and without (black) multiple counting of the interactions between lobes for aspect ratio ∝ D = 1.0 (solid line), 1.01 (dashed line), 1.03 (dash-dash-dotted line), 1.1 (dash-dotted line), 1.3 (dot-dot-dashed line), and 1.5 (dotted line).
The structural and mechanical properties of static packings of elongated particles differ strongly from those for static packings of spherical particles. For example, a number of studies have shown that the packing fraction of static packings first increases from random close packing for spherical particles  with increasing aspect ratio ∝ and then decreases as 1/∝ in the large aspect ratio limit . Also, there are data that suggest that static packings of elongated particles are hypostatic, with fewer contacts than necessary to stabilize all degrees of freedom [5,9]. However, additional studies are required to identify rattler particles, particles with multiple contacts, and contacts between particles that are nearly parallel.
In what follows, we describe methods to calculate contact distances and forces between elongated particles including ellipsoids, spherocylinders, and superquadrics.
For completeness, we summarize the techniques for generating static packings of ellipsoidal particles and studying their mechanical properties here. More detailed expositions may be found in Refs. [20,31,32].
It is assumed that, in both 2D and 3D, particles interact via the pairwise, purely repulsive linear spring potential (Equation 3.3 with ∝ = 2). This equation is generalized so that r ij is the center-to-center separation between nonspherical particles i and j and σ ij is the orientation-dependent center-to-center separation at which particles i and j come into contact as shown in Figure 3.23.
Perram and Wertheim developed an efficient method for calculating the exact contact distance between ellipsoidal particles with any aspect ratio and size distribution in 2D and 3D [3,14,27]. In their formulation, the contact distance between particles i and j is obtained from
Figure 3.21 (a) The vibrational density of states D(ω) for static packings of N = 48 bidisperse dimers with aspect ratio ∝ D = 1.01 and compression Δф = 10−4 (solid line), 10−3 (dashed line), 10−2 (dotted line), and 10−1 (dot-dashed line) with (a) and without (b) multiple counting of the lobe interactions. (c) Smallest frequency ωmin at which D(ωmin = 0.1) versus Δф with (gray) and without (black) multiple counting of the lobe interactions for aspect ratio ∝ D = 1.0003 (solid line), 1.01 (dashed line), 1.1 (dash-dash-dotted line), and 1.5 (dash-dotted line). The dotted line has slope 1/4.
where a and b are the major and minor axes of the ellipsoids, respectively. The expression for σ ij (λ) is minimized over the parameter λ to find the minimum contact distance σ ij = σ ij (λmin) between ellipsoidal particles i and j.
Since σ ij depends explicitly on particle positionsand orientations, we can calculate the forces and torque ,where is the point of contact between ellipsoidal particles i and j. The forces and torques are then used in the soft-particle MD packing-generation method described in Section 3.2.2 to obtain static packings of ellipsoidal particles.
To investigate the mechanical properties of static packings of ellipsoidal particles, one can calculate the eigenvalues of the dynamical matrix and the resulting density of vibrational modes in the harmonic approximation. The dynamical matrix is defined as
Figure 3.22 The vibrational density of states D(ω) for static packings of N = 48 bidisperse dimers at compression Δф =0.03 with aspect ratio ∝ D =1.0003 (solid line), 1.003 (dashed line), 1.03 (dotted line), 1.1 (dot-dashed line), and 1.2 (dash-dash-dotted line) with (a) and without (b) multiple counting of the lobe interactions. The vibrational density of states D(ω) for static packings of N = 48 bidisperse dimers at compression Δф = 0.03 with aspect ratio ∝ D = 1.0003 (solid line), 1.003 (dashed line), 1.03 (dotted line), 1.1 (dot-dashed line), and 1.2 (dash-dash-dotted line) with (c) The location of the low-frequency peak in D(ω), ωmin,versus Δф with (gray) and without (black) multiple counting of the lobe interactions for compressions Δф =0.0003 (solid line), 0.003 (dashed line), 0.03 (dash-dash-dotted line), 0.1 (dash-dotted line), 0.3 (dot-dot-dashed line), and 0.5 (dotted line). The solid (dashed) line has slope 1 (1/2).
Figure 3.23 Definition of the contact distance σ¿j for ellipsoidal particles i and j with unit vectors μ ^ i and μ ^ j that characterize the orientations of their major axes. σ ij is the center-to-center separation at which ellipsoidal particles first touch when they are brought together along r ^ i j at fixed orientation.
In 2D, d f = 3 with
and in 3D for prolate ellipsoids d f = 5 with
See Ref.  for complete expressions for the forces, torques, and dynamical matrix entries in 2D.
A spherocylinder is a cylinder (length l diameter σ) with hemispherical endcaps (see Figure 3.24a). Because of this construction, all points on the surface lie a distance σ/2 from the line segment along the cylinder axis. The question of whether two spherocylinders overlap, and therefore interact through contact forces, is thus reduced to finding the shortest distance between two line segments. If this distance is less than σ (for equal-sized spherocylinders), the spherocylinders overlap and experience equal and opposite repulsive forces. Pournin  has shown that the distance between two points on different line segments is given by
The shortest distance is found by minimizing Equation 3.21 with respect to λ i and λ j . This minimization can be determined analytically and is thus computationally quite efficient. Purely repulsive linear spring forces between spherocylinders can be implemented using Equation 3.3 with d ij replacing r ij . Packings of frictionless spherocylinders can be generated using the successive compression plus energy minimization technique, similar to the algorithms developed in Refs. [24,25] and discussed earlier in this chapter. Sample static packings of frictionless spherocylinders with aspect ratio α = l/σ = 5 and 50 are shown in Figures 3.24b.
Figure 3.24 (a) Two spherocylinders with lengths I i and l j and diameters σ i and σ j .The shortest distance between the central line segments d ij is indicated. When d ij < (σ i + σ j )/2, the two spherocylinders interact. (b and c) Static packings of spherocylinders with aspect ratio ∝ = l/σ = 5 and 50.
Superquadrics are a family of shapes whose surfaces are defined by the equation
In particular, if n x = n y = n z = 2 the particle is elliptical, and as n x ,n y,and n z →∞ the particle takes the shape of a rectangular box with sharp corners. Superquadrics can model particles with concave or convex shapes as well as blunt or sharp edges. Several examples of superquadrics are shown in Figure 3.25.
Figure 3.25 Three examples of superquadric shapes. (a) n x = n y = 2, n z = 50 gives a cylinder; (b) n x = n y = n z = 50 gives a rectangular box; (c) n x = n z = 3, n y = 2 gives a briquette shape.
Superquadrics allow for a far more flexible range of particle shape at the cost of computational efficiency, as the detection of contacts cannot be calculated analytically. Studies of 2D superquadrics [22,23] have established an efficient numerical method that involves minimization in one fewer dimension than the system. Essentially, the equation that defines the surface of the i th particle is rewritten as f i (x, y, z) is less than 1 for all points in the particle, equal to 1 for all points on the surface, and greater than 1 for all points outside the particle. This function is then minimized on the surface of the j th particle, that is, on a 2D surface. The two particles are in contact if. The forces and torques can then be determined with the discrete element modeling (DEM) technique discussed earlier, as with spheres and spherocylinders.
This chapter described computational methods to generate static packings of fric-tionless spheres as well as nonspherical composite particles such as rigid dimers and elongated particles including ellipsoids, spherocylinders, and superquadrics. These methods generate force- and torque-balanced, MS particle packings at jamming onset with only infinitesimal interparticle overlaps. Both hard- (Monte Carlo and Lubachevsky-Stillinger) and soft-particle (dissipative MD) methods were reviewed, and we showed that hard- and soft-particle methods generate the same MS packings of spherical particles at jamming onset although with different probabilities of occurrence. We also reviewed methods to prepare static packings of frictional spherical particles using the Cundall-Strack model. Finally, we presented novel results that identify the regime of aspect ratio and compression where the composite particle model for nonspherical particle shapes yields consistent results for the contact number and density of vibrational modes (in the harmonic approximation).