The statistics of frictional families

The statistics of frictional families

Tianqi Shen    Stefanos Papanikolaou    Corey S. O’Hern    Mark D. Shattuck Department of Physics, Yale University, New Haven, Connecticut 06520-8120, USA Department of Mechanical Engineering & Materials Science, Yale University, New Haven, Connecticut 06520-8260, USA Department of Applied Physics, Yale University, New Haven, Connecticut 06520-8120, USA Benjamin Levich Institute and Physics Department, The City College of the City University of New York, New York, New York 10031, USA

We develop a theoretical description for mechanically stable frictional packings in terms of the difference between the total number of contacts required for isostatic packings of frictionless disks and the number of contacts in frictional packings, . The saddle order represents the number of unconstrained degrees of freedom that a static packing would possess if friction were removed. Using a novel numerical method that allows us to enumerate disk packings for each , we show that the probability to obtain a packing with saddle order at a given static friction coefficient , , can be expressed as a power-series in . Using this form for , we quantitatively describe the dependence of the average contact number on friction coefficient for static disk packings obtained from direct simulations of the Cundall-Strack model for all and .


Granular media are fascinating, complex materials that display gas-, liquid-, and solid-like behavior depending on the boundary and driving conditions. Frictional forces are crucial for determining the structural and mechanical properties of granular media in the solid-like state makse_compress (). For example, friction plays an important role in setting the angle of repose angle (), determining the width of shear bands in response to applied stress shear (); xu (), and enabling arches to form and jam hopper flows pak ().

Figure 1: (Color online) (a) The second invariant of the distance matrix versus packing fraction at jamming onset for static packings of bidisperse disks with generated using the Cundall-Strack model for friction with (circles), (triangles), (squares), and (crosses). Only packings with saddle order and are shown. For , there are packings with  gao (). (b) Close-up of boxed region in (a) with packings (lines) generated using the ‘spring network’ method that originate from the circled packing. (c) All packings (i.e. branch shown as a gray mesh) that are generated from the two highlighted families of first-order saddle packings (dark lines) using the spring network method. packings with the same contact network as that for branch generated using the Cundall-Strack method are also shown as triangles.

For static packings of frictionless spherical particles, it is well known that the minimum contact number required for mechanical stability gao () is , where is number of contacts among particles in the force-bearing backbone of the system and is the spatial dimension. However, at nonzero static friction coefficient , fewer contacts are required for mechanical stability with and in the large- and limits.

Several computational studies have measured the contact number as a function of for packings of frictional disks and spheres using ‘fast’ compression algorithms that generate amorphous configurations silbert1 (); makse (); makse1 (). In particular, these studies find and in the and limits, respectively, for bidisperse disks silbert (); papanikolaou (). For intermediate values of , smoothly varies between and . However, it is not currently known what determines the contact number distribution for each and form of for a given packing preparation protocol. The ability to predict the functional form of the contact number with is important because controls the mechanical somfai () and vibrational bertrand () properties of granular packings.

In this Letter, we develop a theoretical description for static frictional packings at jamming onset in terms of their ‘saddle order,’ or the number of contacts that are missing relative to the isostatic value in the zero-friction limit, . In contrast, previous studies used packings as the reference henkes (). Using a novel numerical procedure (the ‘spring network’ method) that allows us to enumerate packings for each and molecular dynamics (MD) simulations of the Cundall-Strack model cundall () for frictional disks, we show that characterizes the dimension of configuration space that static packings occupy. Frictional packings with contacts form one-dimensional lines in configuration space, packings with populate two-dimensional areas in configuration space, and packings with larger form correspondingly higher-dimensional structures in configuration space. We assume that the probability for obtaining a static packing with saddle order at a given , , is proportional to the volume occupied by force- and torque-balanced th order saddle packings in configuration space. We find that can be written as a power-series in , , where are the normalized coefficients of the power series. Using this form, we are able to quantitatively describe the dependence of the average contact number on the friction coefficient for static disk packings obtained from MD simulations of the Cundall-Strack model over a wide range of and in the large system limit.

We generated static packings of bidisperse (- mixtures of particles with equal mass and diameter ratio ) frictional disks in square cells with periodic boundary conditions using two methods. First, we implemented a packing-generation algorithm in which the system is isotropically compressed or decompressed (followed by energy minimization) to jamming onset gao () at packing fraction that depends on the saddle order. Pairs of overlapping disks and interact via repulsive linear spring forces in the direction of the center-to-center separation vector . We implemented the Cundall-Strack model for the frictional interactions. When disks and come into contact, a tangential spring is initiated with a force that is proportional to the tangential (perpendicular to ) displacement between the disks. The tangential displacement is truncated so that the Coulomb threshold, , is always satisfied. When the disk pairs come out of contact, we set to zero.

Figure 2: (Color online) (a) Zeroth-order saddle () packing of bidisperse frictionless disks with interparticle contacts. To enumerate packings with , the contact between disks and is broken and constrained to have separation (b) , while the other contacts are maintained (visualized as springs) at . The successive compression and decompression packing-generation process gao () is performed with these constraints to create an packing at with only contacts. This process is then repeated as a function of and for the other contacting particle pairs. (c) Contour plot of the minimum friction coefficient required to achieve force and torque balance for a branch of packings in configuration space spanned by the central position of the spring constraining the first broken contact and -component of the spring constraining the second broken contact . The branch of packings emanates from the circled packing. The color scale for increases from dark to light.

We characterize each disk packing in configuration space by plotting the second invariant of the distance matrix, versus , where and are the - and -coordinates of particles and . (Note that is invariant to uniform translations and rotations, as well as particle-label permutations, of the system.) The plot of versus in Fig. 1 (a) for packings with and illustrates several important features. First, in the limit, packings occur as distinct points in configuration space (or versus gao (). Second, as increases, packings form one-dimensional lines in configuration space that emanate from packings. The packings that are stabilized at low are displaced in configuration space from the packings as highlighted in Fig. 1 (b). In contrast, the packings that occur at large approach the packings. Thus, we find that the lengths of the lines increase with . (Fig. 1 (c)) and higher-order saddle packings populate areas and higher-order volumes in configuration space.

We also developed a numerical technique (‘spring network’ method) to enumerate packings at each . The method is best explained using an example. In Fig. 2 (a), we show an packing of frictionless disks with contacts, which corresponds to the packing circled in the lower right corner of the plane in Fig. 1 (a) and (b). To systematically generate packings with contacts, we break one of the contacts in this packing (e.g. the contact between disks and in Fig. 2 (a)) and constrain its separation to be , while the other contacts are constrained to be . With these constraints and as a function of , we implement the successive compression and decompression packing-generation algorithm gao () to find packings at jamming onset, . This procedure is repeated for each of the other contacts in the packing in Fig. 2 (a) to yield the , branches in Fig. 1 (b), and then for each of the packings. As shown in Fig. 1 (b), we find overlap between the branches from the spring network method and the packings generated from simulations of the Cundall-Strack model. and higher-order saddle packings can be obtained using a similar procedure, except multiple contacts are broken, as shown in Fig. 1 (c). Thus, a family emanates from each packing with branches with dimension in configuration space for each saddle order.

Figure 3: (Color online) th root of the volume in configuration space of the collection of th order saddle packings (, squares; , circles) that are stabilized by a friction coefficient . The dashed lines have slope .

We are interested in determining the probability to obtain an th order saddle packing at a given averaged over families. We will assume that is proportional to the volume in configuration space of the th order saddles that can be force- and torque-balanced with friction coefficient or less:


where is a lengthscale required to make Eq. 1 dimensionally correct. To measure , we first employ the spring network method to generate a grid of points for each branch of saddles of order . At each grid point characterized by , we determine the minimum friction coefficient required to achieve mechanical equilibrium for that configuration, using Monte-Carlo moves to search the null-space of the force- and torque-balance matrix shaebani (). The allowed configuration volume is determined by integrating over the contour (as shown in Fig. 2 (c)) such that for a given th order branch. The total admissible volume in configuration space is obtained by summing the volumes over all th order branches. We show in Fig. 3 that scales linearly with for and .

Figure 4: (Color online) (a) The probability to obtain a packing with , , , (peak moving from left to right) as a function of friction coefficient predicted by Eq. 2 with and . (b) for (circles), (squares), (exes), , (diamonds) moving from left to right from the Cundall-Strack model (symbols) for and fits to Eq. 2 (lines) with given by Eq. 3. (c) for the Cundall-Strack model for (squares), (diamonds), and (triangles) and accompanying fits to Eqs. 2 and 3 (lines).

Thus, from the results in Fig. 3, we postulate the following form for the normalized probability for an th order saddle:


where and the highest order saddle is in 2D. The coefficients , where is the number of packings for a given and the normalized coefficients . A reasonable estimate for the number of branches stemming from each packing is the permutations of , . We show in Fig. 4 (a) that Eq. 2 with and yields qualitatively correct results for the measured probabilities to obtain a given th order saddle from MD simulations of the Cundall-Strack model for (Fig. 4 (b)). For example, zeroth order saddles are most highly probable for small , and the highest order saddles are most probable for . However, as shown in Fig. 4 (b), we obtain a much better fit to the data from the Cundall-Strack model using


which indicates an excess of th order saddles for small that likely originates from rattler particles that join the force-bearing network during compression.

The average contact number can be obtained by calculating . Our strategy is to use the results for small systems (i.e. ) to predict and versus for large . We show in Fig. 4 (c) that the predictions from Eqs. 2 and 3 agree quantitatively with with no additional parameters for and obtained from the Cundall-Strack model. Thus, we have developed a method to calculate for large systems by enumerating frictional families in small systems.

The calculations presented in this Letter provide a framework for addressing several important open questions related to frictional packings. For example, why does the crossover from the low- to high-friction limits in the average contact number and packing fraction occur near for disks papanikolaou (); silbert () compared to for spheres silbert (); makse () for fast compression algorithms. In addition, using the methods described above, we will be able to determine how the crossover from low- to high-friction behavior depends on the compression rate and degree of thermalization. Such calculations are crucial for developing the ability to design granular assemblies with prescribed structural and mechanical properties.

We acknowledge support from NSF Grant No. CBET-0968013 (MS) and DTRA Grant No. 1-10-1-0021 (CO and SP). This work also benefited from the facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center and NSF Grant No. CNS-0821132 that partially funded acquisition of the computational facilities.


  • (1) H. A. Makse, D. L. Johnson, and L. M. Schwartz, Phys. Rev. Lett. 84 (2000) 4160.
  • (2) L. Vanel, D. Howell, D. Clark, R. P. Behringer, and E. Clément, Phys. Rev. E 60 (1999) 5040(R).
  • (3) D. Fenistein, J. W. van de Meent, and M. van Hecke, Phys. Rev. Lett. 92 (2004) 094301.
  • (4) N. Xu, C. S. O’Hern, and L. Kondic, Phys. Rev. E 72 (2005) 041504.
  • (5) K. To, P.-Y. Lai, and H. K. Pak, Phys. Rev. Lett. 86 (2001) 71.
  • (6) G.-J. Gao, J. Blawzdziewicz, C. S. O’Hern, and M. D. Shattuck, Phys. Rev. E 80 (2009) 061304.
  • (7) L. E. Silbert, D. Ertas, G. S. Grest, T. C. Halsey, and D. Levine, Phys. Rev. E 65 (2002) 031304.
  • (8) P. Wang, C. Song, C. Briscoe, K. Wang, and H. A. Makse, Physica A 389 (2010) 3972.
  • (9) C. Song, P. Wang, and H. A. Makse, Nature 453 (2008) 629.
  • (10) L. E. Silbert, Soft Matter 6 (2010) 2918.
  • (11) S. Papanikolaou, C. S. O’Hern, and M. D. Shattuck, Phys. Rev. Lett. 110 (2013) 198002.
  • (12) E. Somfai, M. van Hecke, W. G. Ellenbroek, K. Shundyak, and W. van Saarloos, Phys. Rev. E 75 (2007) 020301(R).
  • (13) T. Bertrand, C. F. Schreck, C. S. O’Hern, and M. D. Shattuck, “Vibrations in jammed solids: Beyond linear response,” (unpublished) (2013);
  • (14) S. Henkes, M. van Hecke, and W. van Saarloos, Europhys. Lett. 90 (2010) 14003.
  • (15) P. A. Cundall and O. Strack, Geotechnique 29 (1979) 47.
  • (16) M. R. Shaebani, T. Unger, and J. Kertész, Phys. Rev. E 79 (2009) 052302.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description