Shear-Jamming in Two-Dimensional Granular Materials with Power-Law Grain-Size Distribution

Abstract

Although substantial progress has been made in recent years in research on sheared granular matter, relatively few studies concentrate on the behavior of materials with very strong polydispersity. In this paper, shear deformation of a two-dimensional granular material composed of frictional disk-shaped grains with power-law size distribution is analyzed numerically with a finite-difference model. The analysis of the results concentrates on those aspects of the behavior of the modeled system that are related to its polydispersity. It is demonstrated that many important global material properties are dependent on the behavior of the largest grains from the tail of the size distribution. In particular, they are responsible for global correlation of velocity anomalies emerging at the jamming transition. They also build a skeleton of the global contact and force networks in shear-jammed systems, leading to the very open, “sparse” structure of those networks, consisting of only $\sim 35$% of all grains. The details of the model are formulated so that it represents fragmented sea ice moving on a two-dimensional sea surface; however, the results are relevant for other types of strongly polydisperse granular materials, as well.

1. Introduction

Granular materials are an example of systems in which relatively simple interactions between similar discrete objects (grains, or particles) produce very complex emergent behavior. Extensive experimental and numerical research on granular materials in recent years produced many important insights into the dynamics of those systems. One group of studies has concentrated on the jamming phase transition, revealing new details of the (relatively well understood) isotropic jamming (e.g., [1,2]), as well as the existence of previously unexplored jammed states in systems subject to shear strain [3,4,5,6,7,8,9,10]. However, the behavior of very strongly polydisperse materials in those settings remains very poorly understood. Most works, including those cited above, concentrate on materials with narrow grain-size distributions (GSD). How polydispersity influences the system dynamics close to and at the jamming phase transition remains an open question.

An example of a granular material with a very wide GSD is sea ice, especially close to the ice edge (the so-called marginal ice zone) or, more generally, in regions where, due to the action of wind, ocean surface waves and currents, the ice cover is fragmented into separate floes. A typical example of this ice cover type is shown in Figure 1. Because the vertical dimension (thickness) of the floes is much smaller than their horizontal dimension (diameter), sea ice can be regarded as two-dimensional (2D). The shape of the ice floes may vary from very irregular through polygonal to nearly circular, depending on the external forcing (especially waves) and the ice age and thickness. However, in most situations, the geometrical properties of the floes, like, e.g., the aspect ratio, remain within a relatively narrow range independently of the area and conditions of observation [11,12,13]. More importantly, the observed floe-size distributions (FSDs) are very wide and have power-law tails with an exponent $\alpha <2$ [11,12,13,14,15,16]. Although it is generally acknowledged that the granular nature of fragmented sea ice influences its dynamics (see, e.g., [17]), most large-scale sea ice models treat ice as a viscous-plastic continuum; our knowledge of how and when the processes taking place at a floe level influence the large-scale behavior of sea ice is very limited.

Figure 1. Fragment of a satellite image of fragmented sea ice in the marginal ice zone off the Antarctic Peninsula (source: Landsat [18]).

Figure 1. Fragment of a satellite image of fragmented sea ice in the marginal ice zone off the Antarctic Peninsula (source: Landsat [18]).

This work is a continuation of previous numerical studies on sea ice composed of disk-shaped floes with power-law size distribution [19,20,21]. It examines the behavior of a 2D polydisperse granular material composed of frictional grains under pure-shear deformation (constant packing fraction, or, in the sea-ice nomenclature, ice concentration A). The grains are placed on a frictional substrate (representing the ocean) and interact with each other by means of Hertzian contact forces. Although the details of the model are formulated so that it can represent sea ice moving on the sea surface, the results are relevant in a more general context of sheared, strongly polydisperse granular materials. Therefore, the specific sea-ice terminology is generally avoided in the rest of this paper, with an exception yo some parts of the discussion in the last section.

The paper is structured as follows: the next section contains the description of the model—its assumptions, governing equations and numerical formulation. The results are presented and discussed in Section 3, with an emphasis on those aspects of the model behavior that are related to the polydispersity of the material. In particular, it is demonstrated that grains from the tail of the GSD play a crucial role in the development of the force and contact networks during the jamming phase transition and are responsible for the emergence of domain-wide correlations between velocity anomalies of individual grains. Finally, conclusions are formulated in Section 4.

2. Model Description

2.1. Model Equations

The modeled system consists of $i=1,\dots ,N$ disk-shaped grains with radii $r_i$, thickness $h_i$ and density ${\rho }_i$, occupying a certain two-dimensional region, $\mathcal{S}$. Let us denote the surface area, volume and mass of the i-th grain with $S_i$, $V_i$ and $m_i$, respectively. Obviously, $m_i={\rho }_i{V}_i=\pi {\rho }_i{h}_i{r}_i^2$. The grains move within $\mathcal{S}$, due to both external forcing (e.g., friction against the underlying material) and interactions with neighboring grains. The external forcing acting on the individual grains can be expressed in terms of the density of the surface and body forces, denoted with ${\check{\mathbf{f}}}_{s,i}$ and ${\check{\mathbf{f}}}_{b,i}$, respectively. The net interaction force acting on grain i at a given time instance, t, is a sum of all pairwise interaction forces with grains that are in contact with i at time t. The set of those grains will be denoted with ${\mathcal{C}}_i\left(t\right)$. For $j\in {\mathcal{C}}_i$, ${\widehat{\mathbf{F}}}_{ij,n}$ and ${\widehat{\mathbf{F}}}_{ij,t}$, denote the normal and tangential components, respectively, of the grain-grain interaction force. Under these assumptions, the general form of the equations for the linear and angular momentum of the i-th grain is:

$$\begin{array}{c}\hfill m_i\frac{d{\mathbf{u}}_i}{dt}={\int }_{S_i}{\check{\mathbf{f}}}_{s,i}\mathrm{d}s+{\int }_{V_i}{\check{\mathbf{f}}}_{b,i}\mathrm{d}v+\sum _{j\in {\mathcal{C}}_i\left(t\right)}{\widehat{\mathbf{F}}}_{ij,n}\end{array}$$

(1)

and:

$$\begin{array}{c}\hfill m_i\frac{r_i^2}{2}\frac{d{\omega }_i}{dt}=\mathbf{k}·\left[{\int }_{S_i}\mathbf{r}\times {\check{\mathbf{f}}}_{s,i}\mathrm{d}s+{\int }_{V_i}\mathbf{r}\times {\check{\mathbf{f}}}_{b,i}\mathrm{d}v+\sum _{j\in {\mathcal{C}}_i\left(t\right)}{\mathbf{r}}_{ij}\times {\widehat{\mathbf{F}}}_{ij,t}\right]\end{array}$$

(2)

where ${\mathbf{u}}_i$ denotes the velocity of the grain’s mass center, ${\omega }_i$ (its angular velocity), $\mathbf{k}$ (a unit vector pointing vertically upward), $\mathbf{r}$ (the horizontal distance from the grain’s center) and ${\mathbf{r}}_{ij}$ (a vector pointing from the center of grain i to the contact point with grain j). The interaction forces are calculated based on the Hertzian contact model. The normal force, ${\widehat{\mathbf{F}}}_{ij,n}$, has two components, a contact force (dependent on the overlap between grains) and a damping force (dependent on the relative normal velocity between grains). The tangential force, ${\widehat{\mathbf{F}}}_{ij,t}$, has two component,s as well, namely the shear force (the so-called “history effect” that accounts for the tangential displacement of the interacting grains during contact) and the damping force (dependent on the relative tangential velocity between grains). All four of those forces depend on the effective radius of grains i and j, $r_{ij}=r_i{r}_j/(r_i+r_j)$, and on their material properties, i.e., the elastic modulus, E, and Poisson’s ratio, ν, assumed constant for all grains. Details concerning the formulation of ${\widehat{\mathbf{F}}}_{ij,n}$ and ${\widehat{\mathbf{F}}}_{ij,t}$ can be found, e.g., in [22,23] and in the documentation of the numerical model (see below).

Further details concerning the model formulation are given in [21]. The previous works [19,20,21] stressed the importance of the size-dependent response of individual grains (ice floes) to the forcing acting on them, relevant at low and medium packing fractions ($A\ll 1$). However, the focus of this paper is on a slow deformation of a compact material in or close to the jammed state. Therefore, the simulations described further were performed with a simplified set of equations, without the form-drag terms responsible for the size-dependent response (see [21] for comparison). Furthermore, it is assumed that the frictional substrate is at rest (i.e., both the wind speed and the current speed are zero); the inertial effects (Coriolis term) are omitted, as well. Linearized formulae are used for the grain-substrate (ice-ocean) friction term. Thus, Equations (1) and (2) simplify to:

$$\begin{array}{c}\hfill m_i\frac{d{\mathbf{u}}_i}{dt}=-\pi r_i^2{C}_f{\mathbf{u}}_i+\sum _{j\in {\mathcal{C}}_i\left(t\right)}{\widehat{\mathbf{F}}}_{ij,n}\end{array}$$

(3)

and:

$$\begin{array}{c}\hfill m_i\frac{r_i^2}{2}\frac{d{\omega }_i}{dt}=-\pi \frac{r_i^4}{2}{C}_f{\omega }_i+\mathbf{k}·\sum _{j\in {\mathcal{C}}_i\left(t\right)}{\mathbf{r}}_{ij}\times {\widehat{\mathbf{F}}}_{ij,t}\end{array}$$

(4)

where $C_f$ denotes the friction coefficient (in the context of sea ice, $C_f={\rho }_w{C}_{hw}$, where ${\rho }_w$ is the water density and $C_{hw}$, the water-ice drag coefficient).

As described in [21], the model is based on the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) library [24,25], designed for simulating large systems of interacting objects (particles, molecules, etc.). For the purpose of sea ice modeling, LAMMPS has been extended to disk-shaped particles moving within two-dimensional domains. The Hertzian contact model, available in the official version of LAMMPS, but only for spherical particles, has been modified in order to account for non-spherical grain shape. The modification concerns the relationship between the overlap between the grains and the shape and size of the contact area between them, which, in turn, determines the resulting interaction force. In the case of spherical particles, the contact area is circular, in the case of cylindrical particles, rectangular.

For all N grains, Newton equations of motion (1) and (2) (or, in this particular case, (3) and (4)) are solved by means of the velocity-Verlet integrator.

2.2. Model Configuration and Simulations

In the simulations described in this work, the model domain, $\mathcal{S}$, was rectangular, with length $L_x$, width $L_y=L_x/2$ and surface area $S=L_x{L}_y=\pi {\sum }_{i=1}^N{r}_i^2/A$. In isotropic-compression simulations, periodic boundary conditions were used in both x and y directions; in pure-shear simulations, only along the x-axis, with the grains along the lower model boundary defined as “frozen” (velocity set to zero throughout the simulation) and the grains along the upper model boundary moving with a prescribed velocity ${\mathbf{u}}_i=[u_b,0]$.

A complete list of the model parameters can be found in Table 1. The simulations were conducted for grains with a power-law (PL) GSD, with the mean grain radius $\overline{r}=4.0$ m and the slope of the distribution $\alpha =1.8$, a typical value observed in sea ice (see, e.g., [15,16]). The sample of $N=2·{10}^4$ grain radii was generated with a maximum-likelihood method (e.g., [26], chapter 6.5), which provides an estimate of the most probable value of the i-th element in the rank-ordered sample of finite size N from a given distribution. Thus, deviations from a power law in the tail of the GSD, resulting from the finite sample size, are properly accounted for.

Table 1. Physical and numerical model parameters used in the simulations. GSD, grain-size distributions.

Parameter	Symbol	Value	Units
Grain density	${\rho }_i$	910	kg/m${}^3$
Friction coefficient	$C_f$	1.025	kg/m${}^2$/s
Disk thickness	$h_i$	1.5	m
Mean grain radius	$\overline{r}$	4.0	m
Exponent of the power-law GSD	α	1.8	—
Elastic modulus	E	9.0$·{10}^9$	Pa
Poisson’s ratio	ν	0.33	—
Static yield criterion	μ	0.7	—
No. of grains	N	20,000	—
Speed at the upper boundary	$u_b$	0.2–1.0	m/s
Time step	$\Delta t$	$5·{10}^{-4}$	s

Table 1. Physical and numerical model parameters used in the simulations. GSD, grain-size distributions.

Parameter	Symbol	Value	Units
Grain density	${\rho }_i$	910	kg/m${}^3$
Friction coefficient	$C_f$	1.025	kg/m${}^2$/s
Disk thickness	$h_i$	1.5	m
Mean grain radius	$\overline{r}$	4.0	m
Exponent of the power-law GSD	α	1.8	—
Elastic modulus	E	9.0$·{10}^9$	Pa
Poisson’s ratio	ν	0.33	—
Static yield criterion	μ	0.7	—
No. of grains	N	20,000	—
Speed at the upper boundary	$u_b$	0.2–1.0	m/s
Time step	$\Delta t$	$5·{10}^{-4}$	s

In order to better illustrate the role of the extreme polydispersity in systems with a PL GSD, additional simulations were performed with a narrow, bidisperse (BD) GSD, corresponding to that used by Bi and colleagues [7,8], i.e., with the ratio of the radii of the coarse and fine fraction $r_1/r_2=1.16$ and the respective numbers of grains $n_1=0.2N$ and $n_2=0.8N$. The total number, N, and mean grain radius, $\overline{r}$, were the same as in the model setup with the PL GSD, resulting in $r_1=4.503$ m, $r_2=3.874$ m. In the remaining parts of the paper, all results and comments relate to the PL-GSD simulations unless clearly stated otherwise. Experiments with PL GSD with other values of α from the range $[1.5,2.0]$ produced very similar results and will not be discussed here.

The simulations were performed in two stages: (i) uniform, biaxial compression up to the jamming phase transition; and (ii) pure shear for a set of combinations of packing fraction A and strain rate ϵ values (with $ϵ=u_b/L_y$). The simulations of the second stage were initialized by sampling the results of the first stage at selected values of A and letting the system relax before applying shear strain.

3. Results and Discussion

The general model behavior in uniform-compression simulations is described in [21]. The jamming phase transition in the analyzed case occurs at $A_J\approx 0.918$. It is accompanied by a rapid increase of the internal pressure, p, the fraction of non-rattler grains (defined here as grains with at least two contacts) and the mean contact number, ${\eta }_c$, i.e., changes indicative of the percolation of the contact and force network. Additional simulations performed with different values of the GSD exponent, α, showed that, not surprisingly, the jamming packing fraction, $A_J$, increases with decreasing α (the wider the GSD, the denser the packing fraction attainable), but the course of the jamming transition (for example, the shape of the $p(A-A_J)$ curve) remains almost unaffected. This suggests that the results presented here are relevant for a wider range of model parameters than those actually used in the simulations.

The analysis below concentrates on the pure-shear model runs, with an emphasis on the role of polydispersity in the model behavior close to and at the jamming phase transition. All results have been obtained for packing fractions $A<A_J$, i.e., below the isotropic jamming point. Hence, the “jammed states” in the discussion below refer to regions of shear-jammed and fragile states on the jamming phase diagram proposed by Bi and colleagues [8]. Anticipating the further analysis of the results, the term “jammed state” used throughout the rest of the paper refers to states in which the largest contact network has percolated the whole system in both directions (and which are accompanied by certain global characteristics described further).

Figure 2. Snapshots of contact networks in the modeled system in the unjammed (left; $A=0.890$, $u_b=0.5$ m/s) and jammed (right; $A=0.905$, $u_b=1.0$ m/s) state. For each grain, i, a line is drawn from its center to the center of the neighboring grain, j, if $j\in C_i\left(t\right)$. Grains belonging to the ‘frozen’ and moving boundaries are not shown.

Figure 2. Snapshots of contact networks in the modeled system in the unjammed (left; $A=0.890$, $u_b=0.5$ m/s) and jammed (right; $A=0.905$, $u_b=1.0$ m/s) state. For each grain, i, a line is drawn from its center to the center of the neighboring grain, j, if $j\in C_i\left(t\right)$. Grains belonging to the ‘frozen’ and moving boundaries are not shown.

3.1. Shear Jamming: General Characteristics

The general behavior of the modeled system under pure shear deformation depends on the packing fraction, A, and the strain rate, ϵ [8]. At low A, the system remains in an unjammed state, in which the internal stress is generated via short, binary collisions between neighboring grains. Regions of jammed, more densely packed grains develop only locally (Figure 2); they are short-lived and disperse, due to interactions with the surrounding, more loosely packed regions. Hence, the internal stress level at the system scale remains very low, a few orders of magnitude lower than in the jammed states (Figure 3), when the force network between grains percolates the whole system (Figure 2) and the neighboring grains remain in contact for many seconds or even minutes (see further Section 3.2), i.e., periods of time up to a few orders of magnitude longer than the duration of a typical binary collision.

Figure 3. Time series of the average contact number, ${\eta }_c$ (a), force-network anisotropy ${\eta }_a$ (b), pressure p (c), shear stress τ (d) and the principal angle, ${\theta }_p$ (e), during simulations with: $A=0.908$ and $u_b=1.0$ m/s (blue); $A=0.905$ and $u_b=0.5$ m/s (black); $A=0.905$ and $u_b=0.2$ m/s (magenta).

Figure 3. Time series of the average contact number, ${\eta }_c$ (a), force-network anisotropy ${\eta }_a$ (b), pressure p (c), shear stress τ (d) and the principal angle, ${\theta }_p$ (e), during simulations with: $A=0.908$ and $u_b=1.0$ m/s (blue); $A=0.905$ and $u_b=0.5$ m/s (black); $A=0.905$ and $u_b=0.2$ m/s (magenta).

The large-scale system behavior can be described by means of the properties of the stress and fabric tensors: the pressure $p=({\sigma }_1+{\sigma }_2)/2$ and the shear stress $\tau =({\sigma }_2-{\sigma }_1)/2$ are calculated from the principal stresses, ${\sigma }_1$ and ${\sigma }_2$; the mean contact number ${\eta }_c={\lambda }_1+{\lambda }_2$ and the contact-network anisotropy ${\eta }_a=({\lambda }_2-{\lambda }_1)/{\eta }_c$ are calculated from the eigenvalues of the fabric tensor, ${\lambda }_1$ and ${\lambda }_2$ (see [8,21] for details). Both in the jammed and unjammed states, far from the jamming-transition point, those four large-scale system characteristics—p, τ, ${\eta }_c$ and ${\eta }_a$—remain relatively stable in time, and the system recovers fast from short rearrangement events that sporadically take place (Figure 3). In between those two extremes, the system undergoes rapid changes and shifts from unjammed to jammed states and vice versa (black lines in Figure 3). Between those two extremes, the force networks often have a fragile, “openwork” structure, with relatively large unjammed areas where the stress remains very low and with forces transmitted via long “strands” of approximately linearly aligned grains. As in the case of fragile states observed recently [8], those force networks may span the whole model domain in only one (compressive) direction, giving the material anisotropic strength in response to deformation, which manifests itself in high values of ${\eta }_a$ (see, also, Section 3.2). The present results suggest that, even in constant strain conditions, the fragile states are short-lived, at least in the range of A and ϵ combinations analyzed here.

Apart from the properties of the stress and contact-fabric tensors, a signature of jamming is also present in the grains’ velocity, both means and their anomalies. Let us define ${\mathbf{u}}_m(y,t)$ and ${\sigma }_m(y,t)$, as the mean and standard deviation, respectively, of the velocity of all grains that at time t have their y-coordinate within a certain small distance, δ, from y (i.e., that lie inside a stripe of length $L_x$ and width $2\delta$). Further, let $〈{\mathbf{u}}_m\left(y\right)〉$ and $〈{\sigma }_m\left(y\right)〉$ denote the time mean of ${\mathbf{u}}_m(y,t)$ and ${\sigma }_m(y,t)$ over the whole simulation time and ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$—the velocity anomaly of grain i, i.e., ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)={\mathbf{u}}_i(\mathbf{x},t)-〈{\mathbf{u}}_m\left(y\right)〉$.

Figure 4. Profiles of the average (a) and standard deviation (b) of the x-component of grain velocity in the function of the normalized y-distance ($y=0$ at the “frozen” boundary and $y=1$ at the moving boundary). Results obtained with $u_b=1.0$ m/s.

Figure 4. Profiles of the average (a) and standard deviation (b) of the x-component of grain velocity in the function of the normalized y-distance ($y=0$ at the “frozen” boundary and $y=1$ at the moving boundary). Results obtained with $u_b=1.0$ m/s.

The profiles of $〈{\mathbf{u}}_m\left(y\right)〉$ and $〈{\sigma }_m\left(y\right)〉$ are shown in Figure 4 for a range of A values corresponding to unjammed and jammed states. At low packing fractions, the motion of the grains is confined to the region close to the moving boundary, and a narrow zone of strong shear separates this region from the rest of the model domain, remaining almost at rest. To the contrary, jammed states are characterized by an almost constant velocity gradient $d〈{\mathbf{u}}_m\left(y\right)〉/dy$ and constant standard deviation of velocity $〈{\sigma }_m\left(y\right)〉$, independently on the distance from the moving boundary, i.e., the strain is distributed over the whole system.

In order to characterize the variability of velocity anomalies, it is convenient to define a measure analogous to entropy (as used in statistical mechanics), characterizing the spread of velocity anomalies of individual grains at a given time, t:

$$E\left(t\right)=-c\sum _{i=1}^n\left(p_i{log}_2{p}_i\right)$$

(5)

where n denotes the number of bins of the discrete pdfof ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$, $p_i$ is the probability density of the i-th bin and $c=1/{log}_{2}n$—a normalization constant, introduced so that the maximum value of $E=1$. In order to account for different ranges of ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$ in different model runs, the pdfs were estimated by dividing the range, $[q_{0.01},q_{0.99}]$, into $n=100$ bins of equal width, where $q_{0.01}$, $q_{0.99}$ denote the 1% and 99% quantiles of the data, respectively. Thus, E analyzed here reflects the shape of the pdfs within the respective inter-quantile range, and not their widths, which, as shown in Figure 4b, is much larger in the jammed than in the unjammed states.

Figure 5. Normalized entropy E of the anomalies, ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$, in the function of the grain packing fraction, A (a), and grain size (b). On each box, the central mark is the median, the edges of the box are the 25th and 75th percentiles and the whiskers extend to the most extreme data points not considered outliers. In (b), for the two selected values of A (0.890 and 0.905), the statistics are calculated three times: for all N grains and for the subsets of the 10% largest and 10% smallest grains, respectively. Results obtained with $u_b=1.0$ m/s.

Figure 5. Normalized entropy E of the anomalies, ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$, in the function of the grain packing fraction, A (a), and grain size (b). On each box, the central mark is the median, the edges of the box are the 25th and 75th percentiles and the whiskers extend to the most extreme data points not considered outliers. In (b), for the two selected values of A (0.890 and 0.905), the statistics are calculated three times: for all N grains and for the subsets of the 10% largest and 10% smallest grains, respectively. Results obtained with $u_b=1.0$ m/s.

As can be seen in Figure 5a, E increases with increasing packing fraction A. It has highest values, exceeding 0.85, and lowest time variability (see the boxes and whiskers in Figure 5) in shear-jammed states. In unjammed states, E, most of the time remains within the 0.65–0.7 range. Thus, the range of instantaneous velocity anomalies in jammed systems is significantly larger, even corrected for the width of the respective pdfs. On the other hand, jamming is associated with a transition from local to global correlations of ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$, as illustrated in Figure 6 and Figure 7, showing the linear correlation coefficient, C, between pairs of grains in two selected model runs (for two grains, i and j, C is a Pearson correlation coefficient between the x-components of ${\mathbf{u}}_i^{\prime }$ and ${\mathbf{u}}_j^{\prime }$ over time $t_c=100$ min). At low A, statistically significant correlation of velocity anomalies is observed only between grains within a small spatial distance from each other. At high A, the correlation remains high within the whole model domain. Those two facts—velocity anomalies correlated on the system-scale and high values of E—indicate that in a jammed state, the grains tend to have large velocity anomalies that are of the same sign.

Figure 6. Snapshots of the modeled system for $u_b=1.0$ m/s and the packing fraction $A=0.809$ (a) and $A=0.905$ (b), showing the linear correlation coefficient, C, between the velocity anomalies, ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$, of a selected grain (dark brown, $C=1$) and all other grains in the system. C was calculated for a period of time equal to 100 minutes. Grains belonging to the “frozen” and moving boundaries are not shown.

Figure 6. Snapshots of the modeled system for $u_b=1.0$ m/s and the packing fraction $A=0.809$ (a) and $A=0.905$ (b), showing the linear correlation coefficient, C, between the velocity anomalies, ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$, of a selected grain (dark brown, $C=1$) and all other grains in the system. C was calculated for a period of time equal to 100 minutes. Grains belonging to the “frozen” and moving boundaries are not shown.

Figure 7. The correlation coefficient, C, between the velocity anomalies, ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$, calculated for pairs of grains from a subset of the 10% largest (continuous lines) and 10% smallest (dashed lines) grains in the whole ensemble, in the function of the grain-grain distance. Results of simulations with $u_b=1.0$ m/s and $A=0.890$ (blue), $A=0.905$ (red).

Figure 7. The correlation coefficient, C, between the velocity anomalies, ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$, calculated for pairs of grains from a subset of the 10% largest (continuous lines) and 10% smallest (dashed lines) grains in the whole ensemble, in the function of the grain-grain distance. Results of simulations with $u_b=1.0$ m/s and $A=0.890$ (blue), $A=0.905$ (red).

3.2. The Role of Polydispersity

In systems with power-law GSD, the largest grains occupy a substantial part of the model domain (even with increasing system size N), and it is their locations and relative movement that have a deciding influence on the system as a whole. Sub-regions of the model domain that at a given time instance are filled with small grains can change their shape (and thus react to strain deformation) more easily than assemblies of large grains. In many respects, assemblies of very small grains act as a plastic, easily deformable ‘filler’ occupying empty spaces between very large grains. An analysis of animations illustrating the time evolution of the modeled system reveal that the rapid jamming and un-jamming events mentioned earlier (black curves in Figure 3) tend to be associated with reorganization of the positions of the largest grains. This observation seems confirmed by the fact that, at high packing fractions, the analyzed measures of the grains’ velocity anomalies, like the entropy, E, are strongly correlated to the global instantaneous pressure, p, and shear stress τ and that this correlation is higher for a subset of the largest grains than for the whole system. For example, in the model run with $A=0.905$ and $u_b=1$ m/s, the correlation of E with $log\left(p\right)$ equals 0.83 and 0.95 for, respectively, all and the subset of 10% of the largest grains.

Previous experiments with an earlier version of the model demonstrated that polydispersity plays an important role in many aspects of the dynamics of sea ice composed of floes with power-law size distribution, including the formation of clusters in response to wind [19,20]. Not surprisingly, polydispersity also influences the behavior of the sheared systems studied here. Many global characteristics of the system, including those analyzed above, have different values when they are calculated for a subset of the largest or smallest grains, revealing their different response to the forcing and interactions with neighboring grains. In particular, the entropy, E, of velocity anomalies of the largest grains in an ensemble is higher than the system average at all packing fractions analyzed, i.e., both in jammed and unjammed states (Figure 5b). The emergence of long-range correlations between velocity anomalies at the jamming transition, described in the previous section, takes place almost exclusively due to correlations between the largest grains in the system (Figure 7 and Figure 8). Similarly, at low A, the high values of C within clusters (0.5–0.6 on average) are observed only for pairs of the largest grains. Furthermore, whereas at low A, those values drop rapidly with increasing grain-grain distance, the rate of that decrease is much slower in jammed states (compare the continuous curves in Figure 7), resulting in a shift of the pdf of C towards larger values, representing statistically significant correlation (Figure 8b). To the contrary, the pdfs of C of the smallest grains (in this case, $r_i<1.87$ m) hardly change at the jamming transition, with most values of C remaining at a very low, statistically insignificant level.

Figure 8. pdfsof the correlation coefficient, C, between the velocity anomalies, ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$, calculated for pairs of grains from a subset of the 10% largest (blue) and 10% smallest (red) grains in the whole ensemble. Results of simulations with $u_b=1.0$ m/s and $A=0.890$ (a), $A=0.905$ (b).

Figure 8. pdfsof the correlation coefficient, C, between the velocity anomalies, ${\mathbf{u}}_i^{\prime }(\mathbf{x},t)$, calculated for pairs of grains from a subset of the 10% largest (blue) and 10% smallest (red) grains in the whole ensemble. Results of simulations with $u_b=1.0$ m/s and $A=0.890$ (a), $A=0.905$ (b).

Thus, the increase of the packing fraction, A, towards jamming is accompanied by the growth of clusters of coordinated motion of the relatively small subset of the largest grains (notably, the range of sizes of those grains is still very wide, between 6.5 and 180 m). It is worth noting that similar behavior has been described recently for sheared bidisperse granular systems, in which the dominant dynamical modes were associated with reorganization of grains within localized clusters [10]. Similarly, Weeks and colleagues [27] observed cooperative motion of particles within clusters in colloidal supercooled fluids, with the size of clusters rapidly increasing when the system approached the glass transition.

Figure 9. Selected properties of the contact networks in the modeled system for a number of packing fractions, A: number of contacts of individual non-rattler grains, $n_{c,i}$ (a); $n_{c,i}$ scaled with grain perimeter $2\pi r_i$ (b); percentage of the simulation time when individual grains were non-rattler grains (c); percentage of grains with at least three contacts (d); average contact number ${\eta }_c$ (e); and contact-network anisotropy ${\eta }_a$ (f). In (d–f), the elements of the box symbols are the same as in Figure 5; they reflect the time variability of the analyzed variables during the simulation.

Figure 9. Selected properties of the contact networks in the modeled system for a number of packing fractions, A: number of contacts of individual non-rattler grains, $n_{c,i}$ (a); $n_{c,i}$ scaled with grain perimeter $2\pi r_i$ (b); percentage of the simulation time when individual grains were non-rattler grains (c); percentage of grains with at least three contacts (d); average contact number ${\eta }_c$ (e); and contact-network anisotropy ${\eta }_a$ (f). In (d–f), the elements of the box symbols are the same as in Figure 5; they reflect the time variability of the analyzed variables during the simulation.

Due to obvious geometrical reasons, in strongly polydisperse materials, the number of contacts of individual grains strongly varies. Interestingly, the jamming transition (inferred from the size of the largest connected cluster) in the analyzed cases still takes place when the average contact number, ${\eta }_c$, exceeds the value of three (packing fraction $A=0.905$ in Figure 9e), similarly as in monodisperse and weakly polydisperse systems ([8] and Figure 10). For the large grains from the tail of the GSD, the contact number of individual grains, $n_{c,i}$, is an approximately linear function of their radius (and perimeter), with $n_{c,i}/\left(2\pi r_i\right)\approx 0.05$ in the jammed state (Figure 9a,b), e.g., $n_{c,i}\sim 30$ for $r_i=100$ m. A linear relationship for $n_c\left(r\right)$ has been obtained recently by Shaebani and colleagues [28] in simulations of 2D uniformly compressed, weakly polydisperse systems, in agreement with their mean-field solution. In our simulations, smaller grains often have just one or two neighbors (hence the points on the left side of Figure 9b tend to lie on the $r^{-1}$ curve), and importantly, it is their incorporation into the system-wide contact network that leads to its consolidation at the jamming transition: whereas the largest grains are non-rattler grains most of the time, even in unjammed states, the smallest grains switch at jamming from predominantly freely moving to predominantly non-rattler (Figure 9c; see also Figure 2). Consequently, jamming is associated with a decrease of the mean radius of grains forming the main contact network (not shown). On the other hand, it is the largest grains that build the stable skeleton of the global contact and force network, in the sense that the great majority of grains, predominantly from the left part of the GSD, does not participate in its formation. Even in a jammed state, only $\sim 35$% of grains have three or more contacts (Figure 9d), and even though ${\eta }_c>3$, the mean contact number for the 90% smaller grains (i.e., excluding the 10% largest) is smaller than three. The “sparse” character of the grain-grain contacts is clearly seen in Figure 11. In the jammed state at $A=0.905$, the force network percolates the whole model domain (see also Figure 2), but it has an “openwork” structure, with the largest grains incorporated into long force chains and irregular “cells”, surrounding unjammed regions, usually filled with very small grains. As already mentioned, the assemblies of the smallest grains act as a semi-plastic “filler”, adjusting its shape to the deforming cells of the main force network.

All those properties of the analyzed system are directly related to its extreme polydispersity. In the bidisperse reference case (Figure 10), the jamming transition is much more rapid in terms of the amount of grains that are incorporated into the global contact network: as soon as ${\eta }_c$ exceeds the value of three, roughly 80% of grains become non-rattler grains (Figure 10d,e). Moreover, the coarser and finer fractions contribute similarly to the contact network (compare Figure 10b,c), and the (very small) difference between the average number of neighbors of individual grains from the two fractions (Figure 10a) results simply from the difference between their perimeters.

Back to the PL-GSD simulation, it is also worth noticing that although the contact numbers of the largest grains are high independent of the packing fraction (Figure 9c and Figure 11), in unjammed states, most of those contacts do not form part of stable force strains, but reflect individual collisions as they ‘fight their way’ among smaller neighbors (in Figure 11 at $A=0.890$, most of the lines outgoing from the centers of the largest grains are black, i.e., they lead to grains with a number of contacts lower than three). Such contacts rarely survive more than a few seconds. Indeed, the exceedance probability of contact lifetime is at low A similar for all grain sizes (Figure 12) and only ∼10% of contacts survive for longer than one minute, as compared to ∼25% and 40% of contacts between the smallest and largest grains, respectively, observed in the jammed state.

Moreover, it must be remembered that the condition, $n_{c,i}\ge 3$, alone is not sufficient to stabilize the position of individual grains within the contact network. The arrangement of the contacts around the grain’s perimeter is important, as well. For a given grain, i, this arrangement is determined by a set of vectors, ${\mathbf{r}}_{ij}$, for all $j\in {\mathcal{C}}_i$ (see Section 2.1), which divide the grain into $n_{c,i}$ sectors. The central angle of the widest of those sectors (for brevity, we will call it the maximum contact angle, ${\alpha }_{max}$) provides a useful measure of the above-mentioned stability of the i-th grain within the force network. Obviously, ${\alpha }_{max}<{180}^{\circ }$, possible only if $n_{c,i}>2$, is necessary for stability; ${\alpha }_{max}<{120}^{\circ }$ is possible only if $n_{c,i}>3$.

Figure 10. Selected properties of the contact networks in the bidisperse (BD) model case for a number of packing fractions, A: mean number of contacts of individual non-rattler grains, $n_{c,i}$, from the finer and coarser fractions (a); percentage of the simulation time when individual grains from the finer (b) and coarser (c) fraction were non-rattler grains; percentage of grains with at least three contacts (d); average contact number ${\eta }_c$ (e); and contact-network anisotropy ${\eta }_a$ (f). In (b–f), the elements of the box symbols are the same as in Figure 5; they reflect the time variability of the analyzed variables during the simulation.

Figure 10. Selected properties of the contact networks in the bidisperse (BD) model case for a number of packing fractions, A: mean number of contacts of individual non-rattler grains, $n_{c,i}$, from the finer and coarser fractions (a); percentage of the simulation time when individual grains from the finer (b) and coarser (c) fraction were non-rattler grains; percentage of grains with at least three contacts (d); average contact number ${\eta }_c$ (e); and contact-network anisotropy ${\eta }_a$ (f). In (b–f), the elements of the box symbols are the same as in Figure 5; they reflect the time variability of the analyzed variables during the simulation.

As can be seen in Figure 11b, many non-rattler grains, i.e., those with $n_{c,i}>2$, have ${\alpha }_{max}>{180}^{\circ }$. In fact, ∼14% out of the ∼35% of grains building the percolated contact network (i.e., ∼5% of the total) have ${\alpha }_{max}>{180}^{\circ }$. For comparison, out of 82–83% of grains building the global network in the BD case, only ∼1% are unstable. Figure 13 shows the pdfs of ${\alpha }_{max}$ for the PL and BD cases in shear-jammed states. In the BD simulations, the pdfs have high peaks close to the 120${}^{\circ }$ value, indicating the prevailing—very stable—contact arrangement with $n_{c,i}=3$ and roughly uniform distribution of neighbors around the grain’s perimeter. Notably, the pdfs for the coarser and finer fractions have very similar shapes, with an exception of a small second peak close to 180${}^{\circ }$ for the finer fraction. To the contrary, in the PL simulations, the pdf of ${\alpha }_{max}$ of the smallest grains has a wide maximum shifted towards the unstable region and a long tail corresponding to individual, instantaneous collisions.

Figure 11. Zoomed fragments of the modeled system (top: $A=0.890$; bottom: $A=0.905$) corresponding to the situations shown in Figure 2. Color scale: number of contacts, $n_{c,i}$, of individual grains (dark blue: zero; light blue: one; yellow: two; brown: three or more). Red lines: forces between grains with $n_c\ge 3$; black lines: the remaining forces.

Figure 11. Zoomed fragments of the modeled system (top: $A=0.890$; bottom: $A=0.905$) corresponding to the situations shown in Figure 2. Color scale: number of contacts, $n_{c,i}$, of individual grains (dark blue: zero; light blue: one; yellow: two; brown: three or more). Red lines: forces between grains with $n_c\ge 3$; black lines: the remaining forces.

Finally, another very important property of force networks in sheared systems is their anisotropy, ${\eta }_a$. It has been identified as an order parameter for shear-jammed states, in that it is non-zero in such states and zero in isotropically jammed systems [8]. The values of ${\eta }_a$ obtained in this work, in the BD and, especially, PL simulations, are lower than those reported in [8], which may be related to the details of the contact-model formulation. Characteristically, for the PL case, the ${\eta }_a$ values are roughly 50% higher for the 10% largest grains than the ensemble mean (not shown), again underlying the special role of those grains in shaping the force network structure. The overall variability of ${\eta }_a\left(A\right)$ is, however, similar in both PL and BD simulations, i.e., a rapid drop of ${\eta }_a$ is observed at jamming, preceded by a slight increase when the jamming is approached from below (Figure 9f and Figure 10f). The strongest anisotropy was obtained in those model runs, in which the system underwent strong shifts between unjammed and jammed states (black curves in Figure 3).

Figure 12. Exceedance probability of the contact lifetime (in seconds) for two values of the packing fraction ($A=0.890$, dashed lines; $A=0.905$, continuous lines), calculated for contacts between the 10% largest (blue) and 10% smallest (red) grains.

Figure 12. Exceedance probability of the contact lifetime (in seconds) for two values of the packing fraction ($A=0.890$, dashed lines; $A=0.905$, continuous lines), calculated for contacts between the 10% largest (blue) and 10% smallest (red) grains.

Figure 13. pdfs of the maximum contact angle, ${\alpha }_{max}$, for two model runs: with bidisperse (BD; $A=0.816$) and power-law (PL; $A=0.905$) grain-size distribution. For the BD case, the pdfs are calculated separately for the coarser and finer grain fraction; for the PL case—for the subsets of 10% largest and 10% smallest grains. The vertical dotted and dashed lines mark the values ${\alpha }_{max}={120}^{\circ }$ and ${\alpha }_{max}={180}^{\circ }$, respectively (see the text for a description).

Figure 13. pdfs of the maximum contact angle, ${\alpha }_{max}$, for two model runs: with bidisperse (BD; $A=0.816$) and power-law (PL; $A=0.905$) grain-size distribution. For the BD case, the pdfs are calculated separately for the coarser and finer grain fraction; for the PL case—for the subsets of 10% largest and 10% smallest grains. The vertical dotted and dashed lines mark the values ${\alpha }_{max}={120}^{\circ }$ and ${\alpha }_{max}={180}^{\circ }$, respectively (see the text for a description).

4. Conclusions

The results of the present study suggest that many global characteristics of granular materials with very wide GSD, including those indicative of the jamming phase transition, are determined by the behavior of and interactions between a relatively small subset of the largest grains from the tail of the GSD. They build the “core” of the contact and force networks in the material and, consequently, react in a coordinated manner to strain deformation.

To summarize, the percolation of the contact networks in the analyzed system with PL GSD is associated with the following changes of the global system properties: (i) rapid increase of the entropy of grain velocity anomalies; (ii) emergence of large-scale correlation between the velocity anomalies of the largest grains; (iii) rapid increase of the mean contact number and the fraction of non-rattler grains, accompanied by a rapid decrease of the contact network anisotropy; and (iv) rapid increase of the contact lifetimes, especially between the largest grains. In comparison to less strongly polydisperse systems, the percolated contact networks are built of a significantly smaller subset of grains in stable positions, capable of sustaining non-zero strain.

In the context of sea ice (and presumably other real-world polydisperse granular material, as well), the behavior of the largest grains is relevant from a practical point of view: very likely, it is they that are subject to observation. Measuring equipment is usually mounted on the largest and thickest ice floes, providing stable and relatively safe observational platforms. Similarly, in the analysis of remote sensing (satellite or airborne) images of sea ice, it is the largest floes that can be easily identified and tracked. Thus, it is relevant to understand how the behavior of floes (or, more generally, grains) from the tail of the GSD is different from the behavior of the remaining grains and how it is related to the properties of the granular material as a whole. Even though the simulations described in this paper are highly idealized (no wind, ocean currents, etc.), in view of the fact of how little is known about the influence of granular effects on sea ice dynamics, they provide a starting point for further, more advanced studies.