Force networks and jamming in shear-deformed sphere packings

The emergence of rigidity upon changes of temperature, density, or applied stresses in disordered assemblies of particles is of interest in a wide range of soft matter, from glass formers, gels, foams, and granular matter. Shear jamming of frictional grains presents an interesting special case wherein the application of shear stress leads to rigidity rather than its loss. The formation of self-organized structures that resist shear deformation offers an appealing geometric picture of shear jamming, which nevertheless is incompletely developed, and not well integrated with ideas concerning rigidity in frictionless systems. Exploiting the observation that athermally sheared sphere assemblies develop structural features necessary for shear jamming even in the absence of friction [H. A. Vinutha and S. Sastry, Nature Physics 12, 578 (2016)], we analyze conditions for jamming in such assemblies computationally. Solving force and torque balance conditions for their contact geometry, we show, and validate with frictional simulations, that the mean contact number $Z$ equals $D+1$ (for spatial dimension $D=2,3$) at jamming for both finite and infinite friction, above the “random loose packing” limit density, at variance with previous analyses of frictional jamming. We show that the shear jamming threshold satisfies the marginal stability condition recently proposed for jamming in frictionless systems. Along lines explored in studying covalent glasses, we perform rigidity percolation analysis for $D=2$ and find that rigidity percolation precedes shear jamming, which, however, coincides with the percolation of over-constrained regions, leading to the identification of a regime analogous to the intermediate phase observed in covalent glasses. Together, these results provide a geometric description of shear jamming that relate closely with analyses of jamming, rigidity, and the glass transition in frictionless systems, and thus help develop a unified description of jamming phenomenology in diverse disordered matter.

Figure 1

(a) The average contact number $Z$, (b) stress ${\sigma }_{xz}$, and (c) the average contact force $\langle f_c\rangle$ ($f_c$ includes both normal and tangential components), as a function of strain, shown for different densities. $Z$ values are obtained from sheared configurations, and ${\sigma }_{xz}$ and $f_c$ from force balance solutions. The vertical lines show strain values corresponding to the shear jamming (SJ) transition, where ${\sigma }_{xz}$ and $\langle f_c\rangle$ show discontinuous jumps, which also correspond to $Z=4(=D+1)$. The maroon horizontal lines indicates $Z=4$ and the cutoff used to identify the shear jamming strain in ${\sigma }_{xz}$ and $\langle f_c\rangle$ plots. (d) Threshold shear strains as a function of density, which shows that the shear jamming (SJ) transition, in the limit of infinite friction and finite friction, occurs at $Z=D+1$.

Figure 2

The average contact number $Z$ and the average contact force $\langle f_c\rangle$, for (a) BCM and (b) DEM, $\mu =1$. $Z$ values for BCM ($\mu =1$) are obtained after removing rattler contacts, i.e., contacts with $f_c\le {10}^{-10}$. For the DEM case, only contacts that remain after the simulations are considered. The vertical (dashed) lines shows strain values corresponding to the shear jamming (SJ) transition, where $\langle f_c\rangle$ show discontinuous jumps, and the bold vertical lines show strain values corresponding to $Z=4(=D+1)$. The maroon horizontal lines indicates $Z=4$ and the cutoff used to identify SJ strain in $\langle f_c\rangle$ plots.

Figure 3

(a) Mean angle (${\alpha }_{\text{mean}}$) averaged over all pairs of independent solutions as a function of density. For comparison, we show the mean angle ${\alpha }_{\text{mean}}=41.4^{\circ }$ (maroon horizontal line) between two vectors chosen randomly with all the components positive (first quadrant), which is greater than the mean angle obtained from force solutions at $\varphi =0.627$. As the density is decreased along with the decrease in the null space dimension, see Fig. 8, ${\alpha }_{\text{mean}}$ mean also decreases and approaches zero at the isostatic limit. (b) Stress anisotropy as a function of strain for different densities, indicating saturation across the strain value where $Z$ equals $D+1$ (indicated by dashed vertical lines). Stresses are obtained from the BCM method in the limit of infinite friction.

Figure 4

(a) Contact force distributions obtained from BCM for the steady-state configurations at different densities. Solutions are obtained starting with random initial guesses for the contact forces. Above $\varphi =0.55$, the force distributions acquire jammed-like character i.e., display finite force peaks. (b) Comparison of the jamming strain and $Z=D+1$ with the strain values where the system becomes marginal [shown by a black arrow in (f)]. Different thresholds are consistent with each other. (c) Small-force distributions ($f_c<\langle f_c\rangle$) for different densities in the steady state. (d) Exponent value $\theta$ of small-force distributions shown as a function of density and compared with $\frac{1}{{\gamma }_g}-2$, which needs to be smaller than $\theta$ for stability of the packings. (e) Small-force distribution as a function of strain, shown for $\varphi =0.61$, for different windows $Z$ (or $\gamma$) values. SS in the legend indicates steady-state $Z$ values. (f) Exponent value $\theta$ of the small-force distributions shown as a function of strain and compared with $\frac{1}{{\gamma }_g}-2$, which needs to be smaller than $\theta$ for stability of the packings. The strain values used to plot small-force exponents are the strain value at the lower end of the $Z$ window over which they are averaged. ${\gamma }_g$ is the power-law exponent of $g\left(r\right)$, computed from configurations in a window of $Z$ (or $\gamma$) values (see Fig. 9). Supporting data for densities $\varphi =0.57$ and $\varphi =0.58$ are shown in Fig. 10 and 11.

Figure 5

Shear jamming behavior for the two-dimensional soft disk system. (a) Threshold shear strains as a function of density. Data points marked “$Z=3$” corresponds to strain values where $Z$ reaches $D+1=3$ for each density, and “BCM” corresponds to the shear jamming strain (see Fig. 13). Data points “rigid” and “rigid + stress” correspond to strain values where the percolation probability, for percolation along both $x$ and $y$, reaches the value 0.5 (see Fig. 15). (b) Percolation probability as a function of $Z$, showing the presence of an intermediate phase between rigidity percolation and rigid+stress percolation. The blue shaded region indicates the floppy (F) phase, the green region indicates the intermediate (R) phase, and the red region indicates the stressed or over-constrained phase (R+S). (c) An overlay of the rigid, rigid+stress percolating clusters and the strong force network ($f_c>\langle f_c\rangle$; white bonds), at $\varphi =0.82,\gamma =0.082,Z=2.928$. Orange discs belong to the floppy regions (or small rigid clusters of size smaller than 5), green discs belong to the percolating rigid cluster, and blue discs belong to over-constrained regions. The graphic illustrates a strong correlation between the strong force network and the over-constrained region. (In Fig. 19, we show force network and the over-constrained regions for a series of strain values, illustrating this further.) (d) Distribution of rigid+stress clusters, shown for $\varphi =0.8$, for different values of $dZ=Z_c-Z$ (with which the different curves are labeled). $Z_c$ is computed for each initial configuration corresponding to the percolation of the rigid+stress cluster. The maroon curve represents a slope of $-1.75$, shown for reference. For comparison, for a self-organized rigidity percolation model $\alpha =-1.94$ [27].

Figure 6

(a) Illustration of two particles in contact and the forces exerted by particle $j$ on particle $i$. (b) Average contact and total forces obtained using the null space method and direct minimization for different densities. The plot shows that the force solutions obtained using the null space method are better force balanced.

Figure 7

Comparison of forces from DEM simulations and BCM, shown for different densities. $f_n$ and $f_t$ are the normal and the tangential contact forces.

Figure 8

The null space dimension $D_{\text{ns}}$, scaled with $6N$ (maximum dimension of the matrix $M$), as a function of density for two different system sizes. The green line is a fit to $N=2000$ data (black squares). The plot shows that the fit curve and the $Z$ data of SS configurations identifies $\varphi =0.55$ as the lower density limit, as $Z$ approaches the isostatic limit at $\varphi =0.55$, $Z_c=4(=D+1)$.

Figure 9

Near contact power-law form of $g\left(r\right)$ for different windows of $Z$ (or $\gamma$) values, (a) $Z=[3–4]$ and $\gamma =[0.055–0.088]$, (b) $Z=[4–4.3]$ and $\gamma =[0.088–0.101]$, (c) $Z=[4.3–5]$ and $\gamma =[0.101–0.175]$, (d) $Z=[5–5.175(SS\left)\right]$ and $\gamma >0.175$, at $\varphi =0.61$. The strain values, for data $1/{\gamma }_g-2$, in Fig. 4 of the paper is the strain value at the lower end of the $Z$ window. Data is shown for two values of strain steps ($d\gamma$) used during AQS method. The power-law form extends to at least three decades in all the regions. The power-law exponent does not vary much with densities for the same window of $Z$ values.

Figure 10

(a) Small-force distribution as a function of strain, shown for $\varphi =0.57$, for different windows of $Z$ (or $\gamma$) values. SS implies steady-state $Z$ value. (b) Exponent values of small-force distributions should be than greater than $\frac{1}{{\gamma }_g}-2$, for stability of the packings. The strain values for small-force exponents are the strain value at the lower end of the $Z$ window, obtained from Fig. 2 of the paper. The strain value indicated by the black arrow is where the system becomes marginal.

Figure 11

(a) Small-force distribution as a function of strain, shown for $\varphi =0.58$, for different windows of $Z$ (or $\gamma$) values. SS implies steady-state $Z$ value. (b) Exponent values of small-force distributions should be than greater than $\frac{1}{{\gamma }_g}-2$, for stability of the packings. The strain values for small-force exponents are the strain value at the lower end of the $Z$ window, obtained from Fig. 2 of the paper. The strain value indicated by the black arrow is where the system becomes marginal.

Figure 12

(a) Stress-stress spatial correlation as a function of density. (b) Spatial contact stress correlation at $\varphi =0.61$ along different planes, $xz$ being the shear plane, in the steady state. Shear jammed packings have nonzero stress correlations along the shear plane and zero in other planes. The data are averaged over 500 independent solutions for a given contact network and density.

Figure 13

Mean contact number $Z$, contact force $\langle f_c\rangle$, and stress as a function of strain for the infinite friction case for the two-dimensional packings. $\langle f_c\rangle$ shows a jump at the shear jamming transition. The contact forces are computed from the force solutions obtained using the BCM method. The maroon horizontal line in the $Z$ vs. $\gamma$ plot marks $Z=3$, and the horizontal lines in the stress and $\langle f_c\rangle$ plots marks the cutoff used to identify the shear jamming strain. The vertical lines marks shear jamming (SJ) strain values corresponding to sudden increase in $\langle f_c\rangle$ value, shown in Fig. 5 of the paper. The shear jamming strain values identified using $\langle f_c\rangle$ also marks the strain at which the stress shows a sudden jump.

Figure 14

(a) Jamming phase diagram for the two-dimensional system at finite friction $\mu =1$, which indicates that $\varphi =0.81$ is the lower density limit for the physical value of friction coefficient of $\mu =1$. $Z$ and $\langle f_c\rangle$ as a function of strain, shown for different densities, for (b) BCM ($\mu =1$) and (c) DEM ($\mu =1$) protocols. The shear jamming strain values obtained from BCM and DEM methods are close to each other and closely corresponds to $Z=3(=D+1)$. The vertical lines (dashed) shows strain values corresponding to the shear jamming (SJ) transition, $\langle f_c\rangle$ show discontinuous jumps, and the bold vertical lines correspond to $Z=4(=D+1)$.

Figure 15

(a) Percolation probability of rigid clusters as a function of strain $\gamma$ for different densities for the two-dimensional system for $N=2000$. (b) Percolation probability of over-constrained regions (rigid+stress) as a function of strain $\gamma$ for different densities. The label one refers to the case where the percolation probability is computing considering spanning cluster along one direction, whereas all refers to the case where the percolation probability is computed requiring that a spanning cluster is system spanning in both the $x$ and $y$ directions. The rigid+stress percolation strain values in Fig. 5 of the paper correspond to a percolation probability of 0.5 (maroon horizontal line) of the all case.

Figure 16

Percolation probability as a function of $Z$, shown for (a) $\varphi =0.79$, (b) $\varphi =0.81$, and (c) $\varphi =0.82$ for the two-dimensional system, $N=2000$. The data shown is for the all case, where percolation probability is computed requiring a spanning cluster along both the directions. These plots show the presence of an intermediate phase for all densities $\varphi \ge {\varphi }_{\text{RLP}}=0.79$.

Figure 17

Rigid percolation and rigid +stress percolation data along one and all ($x$ and $y$) directions occurs at different $Z$ values. The percolation of rigid+stress clusters that percolate along one direction can support stress transmission along the direction of percolation, which indicates the presence of fragile force networks in the region between rigid percolation and rigid+stress percolation or shear jamming transition.

Figure 18

(a) Stress anisotropy (SA) and (b) fabric anisotropy (FA) as a function of strain, shown for different densities. The force balance solutions are obtained from BCM method in the limit of infinite friction. The stress anisotropy plot shows that in the strain window of rigid percolation and $Z=3(D+1)$ the stress anisotropy begins to saturate, which suggests the presence of fragile force networks in this region or strain window.

Figure 19

The left panels show the evolution of rigid and over constrained (rigid+stress) clusters as a function of strain (red network—biggest rigid cluster, other rigid clusters are not shown; blue filled circles—over-constrained discs; and gray network—floppy regions and other smaller rigid clusters). The right panels show the evolution of contact forces $\langle f_c\rangle$ as a function of strain ($f_c>\langle f_c\rangle$—red bonds; $f_c<\langle f_c\rangle$—green bonds; and $f_c<{10}^{-14}$—gray bonds). Data shown are for $\varphi =0.82$. The force balance solutions (null space method) and rigidity analysis is shown for the same configurations. Contact forces $f_c>\langle f_c\rangle$ are mostly concentrated on the particles that belong to the over-constrained regions and there is a strong spatial correlation between the force network and the rigid network.