Jamming of Bidisperse Frictional Spheres

By generalizing a geometric argument for frictionless spheres, a model is proposed for the jamming density $\phi_J$ of mechanically stable packings of bidisperse, frictional spheres. The monodisperse, $\mu_s$-dependent jamming density $\phi_J^{\mathrm{mono}}(\mu_s)$ is the only input required in the model, where $\mu_s$ is the coefficient of friction. The predictions of the model are validated by robust estimates of $\phi_J$ obtained from computer simulations of up to $10^7$ particles for a wide range of $\mu_s$, and size ratios up to 40:1. Although $\phi_J$ varies nonmonotonically with the volume fraction of small spheres $f^s$ for all $\mu_s$, its maximum value $\phi_{J,\mathrm{max}}$ at an optimal $f^{s}_{\mathrm{max}}$ are both $\mu_s$-dependent. The optimal $f^{s}_{\mathrm{max}}$ is characterized by a sharp transition in the fraction of small rattler particles.

≃ 0.55 [14,16,17]. Conversely, bidispersity in particle sizes can increase φ J to a theoretical maximum φ J ≃ 0.87 for frictionless spheres [4-6, 9, 11, 12, 18]. However, only limited studies [19] have explored the combined role of friction and dispersity on the jamming of spheres. A key impediment is the difficulty of simulating large-scale, mechanically stable, jammed packings of frictional particles in the limit of marginal rigidity [20] using previously established methods that render such packings prone to instabilities at low confining pressures near φ J [13]. Additional challenges include inefficient neighbor finding algorithms for large particle size ratio systems [21,22].
Dispersity in particle sizes results in a rich structural diversity in particulate materials such as granular materials [11,[23][24][25], colloids [26,27], emulsions [28] and geophysical materials [29]. Maximizing φ J in such materials through size dispersity is important in various applications such as battery electrodes [30], cement [31] and chocolate [32]. The mechanics of such materials strongly correlates with φ J [12,33,34], and size dispersity provides a powerful knob to optimize their mechanical properties. Furthermore, φ J also critically governs the equilibrium [24,35,36] and nonequilibrium [37] rheology of dense, bidisperse particulate materials. Therefore, beyond their fundamental jamming characteristics, a robust estimation of φ J for bidisperse particulate materials with realistic interparticle interactions is important towards advancing our knowledge of their mechanics and * isriva@lbl.gov rheology.
Here, we use pressure-controlled simulations to simulate mechanically stable jammed packings of frictional, bidisperse spheres for a wide range of particle size ratios 2 ≤ α ≤ 40 and volume fraction of small particles , where d s and d l are the diameters, and N s and N l are the numbers of small and large spheres, respectively. The effect of µ s on φ J is demonstrated through µ s -dependent state diagrams on (φ J , f s ) axes by varying the coefficient of friction over several orders of magnitude 0 ≤ µ s ≤ 0.5, where µ s = 0 corresponds to frictionless particles. Similar to previous findings [9,11,12], we find that φ J varies nonmonotonically with f s and attains a maximum φ J,max at f s max . At large α, this variation exhibits a sharp transition at f s max . Although φ J varies similarly with f s for all µ s , the optimal φ J,max and f s max depend significantly on µ s . The optimal φ J,max decreases systematically with µ s , while requiring a larger f s max to attain its optimal value. Through these findings, we propose a generalized Furnas model [38]-originally proposed for notionally placing small particles in the available volume left by large particles without reference to mechanical constraints or particle properties other than infinite size ratio-that accurately predicts µ s -dependent φ J for large α, although small deviations from the model are observed for frictional packings at f s < f s max and very large α. Remarkably, the generalized model only requires µ s -dependent monodisperse jamming density φ mono J (µ s ) as an input to accurately predict φ J ; thus φ mono J (µ s ) appears to encode the mechanical stability constraint of the packing in such a way that the available volume for small frictional particles within the space left over by large frictional particles is correspondingly adjusted. Lastly, we find that for f s < f s max all the small particles rattle in the void space of the jammed network of large particles, and the fraction of such rattlers drops sharply at f s max . Unlike frictionless particles, the rattler fraction does not drop to zero for frictional particles.
Three key features distinguish the present simula-tions from similar previous works: (i) we use pressurecontrolled simulations that guarantee mechanical stability, unlike commonly used volume controlled methods [2,39]; (ii) we create packings up to an unprecedentedly high α = 40 in order to probe the large size ratio limit; as a result, we have created the densest known bidisperse jammed packings for frictionless and frictional spheres; (iii) we provide statistically robust estimates of φ J by fixing a large value of N l = 1000 in all simulations. This constraint resulted in a large number of small particles up to N s ≃ 3×10 7 , which were simulated through large-scale discrete element method using the software LAMMPS [40]. A typical jamming simulation starts with a dilute nonoverlapping collection of spheres of material density ρ = 1 at φ = 0.05 containing N l large particles and N s small particles determined from the chosen values of α and f s . The initial simulation cell is periodic and cubic; however, its triclinic nature allows shear distortions in addition to volume deformation. At t = 0 a pressure p a is applied such that the total stress on the system is σ a = p a I, where I is the identity matrix. Particles interact with damped Hookean and the tangential frictional forces. The tangential spring stiffness k n is set equal to the normal spring stiffness k t , and a Coulomb coefficient of interparticle friction µ s sets the sliding frictional force [41,42]. The damping at contacts is set by normal and tangential velocity damping coefficients γ n,t = 0.5. Time is normalized by the characteristic timescale of collision t c = π(2k n /m s − γ 2 n /4) −1/2 , where m s is the mass of the small particle. The simulation time step is set to 0.02t c . The stress in this setup is scaled by k n /d s and energy is scaled by k n d 2 s . Under the action of p a , the system steadily compacts and its internal pressure p steadily increases towards p a , as shown in Fig. 1(b). The applied pressure is fixed at a low value p a = 10 −4 in all simulations to model jamming in the asymptotic hard-particle regime [43]. The inertial equations of motion for the particles and the simulation cell are described in detail in Ref. [41]. After transient evolution, the system eventually achieves a jamming density φ J when p balances p a , as shown in Fig. 1(b). We set a criterion to terminate a simulation when the average kinetic energy per particle falls below 10 −11 or if the simulation has proceeded for at least 5 × 10 6 and up to 10 7 steps, which were found sufficient to achieve a stable jammed packing within a reasonable computational time. Beyond pressure, the pressure-controlled simulation method also ensures that the deviatoric internal stress is equal to its externally applied value upon jamming, i.e., zero [2,14]. This is achieved by allowing shear distortions of the simulation cell [14,41]. As a result, mechanically stable jammed packings are created at a specified pressure p a , similar to soft particle jamming using the variable-cell structural optimization method [2,44], and strictly-jammed hard particle packings using the adaptive shrinking cell method [6].
Traditional neighbor finding and interprocessor com- munication algorithms in particle-based simulations are inefficient for modeling systems with large size dispersity [21,22]. Here, we adapt and implement an algorithm originally developed for colloidal mixtures [21] to simulate bidisperse granular systems for large size ratios up to α = 20. We also include results from jamming simulations at very large α = 40 by adapting a recently developed neighbor finding algorithm that efficiently simulates granular systems of large size ratios [45]. An example of a α = 40 packing that would have been computationally prohibitive to simulate using traditional neighbor finding methods is shown in Fig. 1(a). About a century ago, Furnas predicted a nonmonotonic relationship φ J (f s ) in the limit α → ∞ based on a simple geometric model [38]. At low f s up to a critical The predictions of this model have been tested in simulations of bidisperse, frictionless spheres up to α = 12 [9,11]. Using large scale simulations, we created jammed, bidisperse, frictionless packings up to α = 40. Figure 2(a) shows the variation of φ J with f s for various α along with the predictions from the Furnas model. At low α, φ J varies continuously and nonmonotonically with f s , and this variation exhibits a sharp transition at high α, similar to previous observations [9,11]. For α = 20, φ J follows the predictions of the Furnas model and closely approaches φ J,max = 0.869 at f s max = 0.265, which were estimated by the Furnas model using φ mono J ≃ 0.638 obtained at f s = 0. We also simulated a packing for α = 40 containing N s ≃ 2.3×10 7 particles at the purported op- timal f s = 0.265, as shown in Fig. 2(a). The value φ J = 0.857 obtained in this large simulation is, to our knowledge, the densest bidisperse jammed packing created using computer simulations.
A similar variation of φ J with f s is observed for frictional particles, as shown in Fig. 2(b) for µ s = 0.3. Similar to frictionless particles, φ J varies continuously at low α, whereas the variation exhibits a sharp transition at higher α. However, φ J at all α and f s is substantially lower than φ J for frictionless particles. For monodisperse spheres (f s = 0), we obtained φ mono J ≃ 0.588. In addition to lower φ J compared to frictionless particles, the maxima in φ J occurs at a higher f s . Following the data for α = 20 in Fig. 2(b), we find that the maximum φ J = 0.823 occurs at f s = 0.29, which is higher than f s required to maximize φ J for frictionless particles.
Such differences in φ J (f s ) based on µ s are explained by generalizing the Furnas model to include µ s -dependence of φ mono J . At low f s , large particles form a jammed network, but at a reduced φ mono J (µ s ) [13]. Upon increasing f s , small particles fill the voids of the jammed network, but only up to the same reduced φ mono J (µ s ), at which point φ J,max is obtained at f s max . Using φ mono J = 0.588 for µ s = 0.3 in the generalized Furnas model, f s max = 0.292 and φ J,max = 0.83 are predicted, which are remarkably close to the values obtained from simulations shown in Fig. 2(b) for α = 20. A large jamming simulation for α = 40 containing N s ≃ 2.6 × 10 7 particles for f s = 0.292 results in φ J = 0.83, as shown in Fig. 2(b), which confirms these predictions.
For the largest simulated frictional packings at α = 40, φ J is slightly greater than the model prediction for f s < f s max , and the deviation from the model increases with f s , as shown in Fig. 2(b). A similar but reduced deviation from the model is also observed for α = 20 packings. This is caused by the large particles packing at a density greater than φ mono J (µ s ), but still bounded by the frictionless φ mono J . This effect vanishes for f s > f s max when small particles also participate in jamming, and φ J agrees well with the model. As monodisperse particles can jam at a range of densities depending on the protocol [17,46,47], we expect that in the region near f s max , the same would be true for bidisperse packings.
To test the predictions of the generalized Furnas model for φ J,max as a function of µ s , we first obtain monodisperse φ mono J (µ s ) shown in Fig. 3(a). These values are consistent with previous studies [13,14,42]. Upon substituting φ mono J in the generalized Furnas model, theoretical estimates for φ J,max are obtained, as depicted in Fig. 3(a). The theoretical predictions are tested against φ J,max obtained from simulations for various α, as shown in Fig. 3(a). At lowest α = 2, the calculated φ J,max are much lower than the theoretical maxima, but still substantially higher than φ mono J . As α increases, the calculated φ J,max steadily increases towards its theoretical maxima and is within 2% of its maximum value for α = 20. For high friction, α = 20 appears to be sufficient to nearly attain the theoretical maxima φ J,max . Jamming at low friction requires much longer simulations to completely remove the kinetic energy of the particles. The presence of friction between particles provides additional constraints to the particle motion [14], thus enabling a quicker approach to jamming. Similar to φ J,max , the theoretical predictions for µ sdependent f s max compare well with simulations for α = 20, as shown in Fig. 3(b). As µ s increases, a larger fraction of small particles is required to achieve the highest packing density, with the highest friction case requiring 3% higher volume of small particles than the frictionless case.
Previous studies on jammed bidisperse frictionless particles have demonstrated a sharp transition in the structural properties across f s max [9,11]. Upon increasing f s at high α, the fraction of small particles that are rattlers drops rapidly at f s max [9]. We analyze the fraction of small rattler particles to highlight a similar structural transition for frictional bidisperse packings. A particle is a rattler if it has too few contacts to contribute to the mechanical stability of the packing, i.e., the particle is lo- cally unjammed [48]. We follow a recursive method [48] to identify if a particle i is a rattler based on the criteria: Z i < 3 for frictionless particles and Z i < 2 for frictional particles, which emerges from constraint counting arguments [14]. In Fig. 4, the fraction of small rattler particles N r s /N s is plotted as a function of f s for four values of µ s , where N r s is the total number of small rattler particles. Below a critical µ s -dependent f s that corresponds well with f s max shown in Fig. 3(b), almost all of the small particles are rattlers within a jammed network of large particles. As f s is increased beyond f s max , the fraction of small rattlers rapidly drops, and the majority of small particles participate in the mechanical backbone of the packing. Nearly all small particles are locally jammed for frictionless packings, whereas up to 22% small particles are rattlers at high friction. These findings are consistent with previous studies on rattlers in monodisperse frictional packings near the jamming point [14,49]. The fraction of rattlers at high friction is similar to its value for monodisperse packings [14], thus indicating that small particles form a percolating jammed network with suspended (and locally jammed) large particles.
We have demonstrated that φ J for bidisperse frictional particles is well-predicted by a simple µ s -dependent model for large size ratios. The model compares well with simulations that were used to create large mechanically stable bidisperse frictional packings at unprecedent- edly high φ J . For frictional packings at f s < f s max and α 20, φ J is slightly larger than the model prediction, possibly due to the lubricating effect of the small particles. Friction is found to reduce φ J for all f s . At large α, φ J varies nonmonotonically with f s and exhibits a sharp transition at f s max , where its value is maximized. The f s max that maximizes φ J increases with µ s , with frictional contacts requiring up to 3% more small particles to achieve φ J,max . The sharp transition in φ J is also observed structurally, where the fraction of small rattler particles rapidly drops across f s max . The results presented here open avenues to extend our understanding of the physics of jamming-particularly the structure and mechanical stability of frictional packings [50,51]-to bidisperse frictional packings at large size ratios. Our ongoing work on including additional modes of friction, such as rolling and twisting, support an extended generalization of the Furnas model that takes into account φ mono J in the presence of these additional frictional modes [14]. This suggests that a century-old model may be a valuable predictor of φ J in mechanicallystable bidisperse particulate packings spanning a wide variety of interparticle interactions. Lastly, a recent study has discovered an additional jamming transition at low f s for large α in bidisperse frictionless particles [11]. Although such a transition was not observed in our pressure-controlled jamming simulations, we can not preclude its presence without a careful analysis of the path dependence of jamming in our systems.
Helpful discussions with Andrew Santos and Dan Bolintineanu are acknowledged. I. S. acknowledges support from the U.S. Department of Energy (DOE), Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program under contract No. DE-AC02-05CH11231. This work was performed at the Center for Integrated Nanotechnologies, a U.S. DOE and Office of Basic Energy Sciences user facility. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. DOE's Na-tional Nuclear Security Administration under Contract No. DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. DOE or the United States Government.