Granular materials and disordered systems are commonly encountered in natural systems and industrial processes, such as landslides, pyroclastic flows, fresh concrete, and pharmaceutical and chemical particulates [1]–[4]. In recent decades, extensive investigations have been conducted to describe the rheology[5], [6], compaction [7], [8], and jamming transition [9], [10] of granular materials. The proposal of the inertial number \(I_c\) and the viscous number \(I_v\) suggests that the dynamics of granular flows is governed by the combination of three time scales: \(t_i = d_p/\sqrt{\sigma_n/\rho_p}\), \(t_v = \eta_f/\sigma_n\), and \(t_M = 1/\dot{\gamma}\), where \(t_i\), \(t_v\), and \(t_M\) are inertial, viscous, and macroscopic time scales, \(d_p\) and \(\rho_p\) are the particle diameter and density, \(\sigma_n\) is the pressure, and \(\eta_f\) denotes the fluid viscosity if the system is submerged or fully saturated.

However, the current rheological framework cannot fully capture the influence of either particle shapes or size distributions. The grain size distribution (GSD) plays an important role in the rheology of granular systems, but previous studies often provide seemingly contradictory opinions, where some studies suggested that the existence of small particles initiates a lubrication effect that helps the system increase its mobility [11], while others believed that a wider GSD induces a denser packing fraction that increases the shear strength of granular systems and makes the system easier to transition to a jammed state [12].

[13] studied the rheology of bidispersed granular flows using the discrete element model (DEM) [13], [14] and suggested that the pressure of a granular system is determined not only by its size distribution, but also by the effective free volume per particle, which was adopted from the jammed packing of monosized particles. They stated that the dimensionless pressure, \((1/I_c)^2\), scales as a function of the coordination number, \(Z_c\). [15] studied the rheological behavior of polydispersed granular systems, which concluded that such systems still behave similarly to monodispersed systems if the right length scale is chosen, but the solid fraction can be quite different due to polydispersity. Recently, [11] also studied the behavior of granular systems with fractal distributions, and found that the shear strength of confined flow decreases at large fractal dimensions, which can be further linked to help decipher the mechanisms within the highly localized basal zones and the high mobility of general geological avalanches.

In this Letter, we focus on disordered physical systems with various continuous GSDs, where sheared granular materials with three different cumulative distribution functions (CDF) are examined using DEM. We show that, unlike the length scale proposed in [15], a simply defined “average" particle size, which has a clear physical root in the contact probability, already unifies the rheological behavior of dry granular systems with a wide range of different particle size distributions. Further analyses are devoted to better understand the critical phenomena when a granular system with different GSDs and frictional properties is approaching the”jamming" transition.

Methodology and simulation setup In this work, the classical DEM is used to simulate sheared granular systems, where the particles are spherical and follow a linear spring-dashpot contact model. The normal contact force \(F_{ij}^n\) is calculated as a Hookean contact with energy damping. The magnitude of the tangential contact force \(|k_t\delta_{ij}^t|\) cannot exceed its friction limit, \(\mu_p|F_{ij}^n|\), where \(\mu_p\) is the particle frictional coefficient and for most cases we set \(\mu_p=0.3\). The system has no rolling resistance. The granular system is located between two sawteeth plates and is sheared with constant pressure and shear rate. For different simulations, we vary the pressure \(\sigma\) from 50 to 1000 Pa and the shear rate \(\dot{\gamma}\) from 0.01 to 1.0 s\(^{-1}\), and measure the mean value of the key variables when the sheared system is already in steady state, shown in Appendix 1, where its solid fraction \(\phi_s\), pressure \(\sigma\), shear stress \(\tau\), and velocity profile become stable. We can obtain the effective frictional coefficient \(\mu_{\mathrm{eff}} = \tau/\sigma\), the effective viscosity \(\eta_{\mathrm{eff}} = \tau/\dot{\gamma}\), and the average inertial number \(I_c\) from each simulation for further investigation. The plot of a shear cell is shown in the inset of FIG. 1 (d). The particle size distributions used in this work, i.e., Fuller distribution and two power-law distributions, are elaborated in detail in Appendix 2.

Rheological behaviors According to the classical \(\mu(I)\) rheology of dry granular systems, both \(\mu_{\mathrm{eff}}=\tau/\sigma\) and \(\phi_s\) should scale with the inertial number \(I_{ca} = \dot{\gamma}d_p/\sqrt{\sigma/\rho_p}\), where \(d_p\) is the particle diameter and also a characteristic length scale. [6] proposed that the length scale plays an important role in determining the rheological behavior of granular systems, and we further argue that this length scale should be adjusted to reflect changes in particle size distributions. FIGs. 1 (a-b) demonstrate the \(\mu_{\mathrm{eff}}\sim I_{ca}\) relationship for systems with GSDs of D1, D2, and D3, which are clarified in Appendix 2, whereas FIG. 1 (c) shows the \(\phi_s\sim I_{ca}\) relationship, where \(I_{ca}\) represents the conventional inertial number that uses the average particle diameter of the system as the characteristic length.

FIGs. 1 (a-c) shows that, even though the \(\mu_{\mathrm{eff}}\sim I_{ca}\) and \(\phi_s\sim I_{ca}\) relationships collapse for systems with any specific GSD, changing the particle gradation leads to deviations in both \(\mu_{\mathrm{eff}}\) and \(\phi_s\). FIGs. 1 (a) and (b) show that systems with different particle gradations result in different transitional \(I_{ca}\), which denotes a transition from quasi-static to inertial flow. The change in this transitional \(I_{ca}\) implies that the characteristic length of the system may also change due to changes of GSDs. [13] reported that the average particle size should be calculated based on each particle volume. Thus, we define the characteristic length of the system as \(d_{aw} = [\sum_{i=1}^{N_p}( d_{pi}V_{pi})] / (\sum_{j=1}^{N_p}V_{pj})\) and particles with different volumes are assumed to have different weight in calculating this “average" particle size. Thus, a new inertial number can be defined as \[\begin{align} I_{cw} = \dot{\gamma}d_{aw}(\rho_p/\sigma)^{1/2}\;. \end{align} \label{eq:Icw}\tag{1}\]

After adopting the new inertial number, as shown in FIGs. 1 (d,e), the relationships between \(\mu_{\mathrm{eff}}\) and \(I_{cw}\) collapse onto a single master curve, indicating that the characteristic length, \(d_{aw}\), reflects some physical and topological nature of the sheared granular systems, and choosing an appropriate length scale helps unify the frictional rheology of dry granular systems with any size distribution. However, we admit that \(I_{cw}\) fails to unify the \(\phi_s\sim I_{cw}\) relationship, although FIG. 1 (f) looks much better than FIG. 1 (c). For systems with different GSDs, their transitional \(I_{cw}\) seems to be unified, and the only difference lies in the \(\phi_s-\)axis, i.e., moving the \(\phi_s\sim I_{cw}\) relationships upward or downward can also unify this relation. For a series of simulations with the same GSD, the relationship between \(\phi_s\) and \(I_{cw}\) can be written as \(\phi_s = \phi_c/(1+\alpha_{\phi}I_{cw})\), where \(\alpha_{\phi} \approx 0.32\). Thus, if we change the vertical coordinate into \(\phi_s/\phi_c\), where \(\phi_c\) is the critical solid fraction for each GSD, we can obtain a master curve of \(\phi_s/\phi_c = 1/(1+\alpha_{\phi}I_{cw})\), shown as the blue dashed curve in the inset of FIG. 1 (f).

Critical length scale and critical solid fraction Our choices of characteristic length and critical solid fraction have deep physical and geometrical reasons, even though they seem to be deliberately fitted. On one hand, \(d_{aw}\) can be related to the coordination number of each polydispersed granular system. The characteristic length should be linked to particle volumes because they are directly related to the number of contacts. Intuitively, a particle with larger size has a greater chance to contact with more particles, indicating that such a particle plays a more important role in bearing strong force networks and determining the macroscopic behavior of the whole system. To validate this, we calculated another length scale, \(d_{ac}\), based on contact numbers of each particle, where \(d_{ac} = [\sum_i(d_{pi}N_{ci})]/(\sum_j N_{cj})\) and \(N_{ci}\) is the contact number of particle \(i\), and plotted its relationship with \(d_{aw}\) in the inset of FIG. 1 (b). This shows that \(d_{aw}\) and \(d_{ac}\) are linearly related, which verifies that \(d_{aw}\), which can be easily obtained from GSD, has a deep root in contact statistics.

\(\phi_c\) varies with the GSD and frictional properties of particles. It is clearly neither the jamming solid fraction, \(\phi_{J}\) [10], nor the random close packing fraction, \(\phi_{RCP}\) [16]. However, [16] and his follow-up research on the random close packing of the polydisperse granular system [17] encourage us to argue that \(\phi_c\) is a function of \(\phi_{RCP}\) and \(\mu_p\). Since we kept \(\mu_p\) constant, \(\phi_c\) is predominantly controlled by \(\phi_{RCP}\), which should be related to the GSD of a granular system. Thus, we could obtain a functional relationship between \(\phi_c\) and the GSD.

We note that \(\phi_c\) for each GSD should not be obtained from fitting of the \(\phi_s\sim I_{cw}\) relationship, since the functional form of the \(\phi_s\sim I_{cw}\) relationship is based on phenomena, not physics. Fortunately, [16] derived the relationship between effective viscosity \(\eta_{\mathrm{eff}}=\tau/\dot{\gamma}\) and \(\phi_s\) from first principles, which yields \(\eta_{\mathrm{eff}}\sim |\phi_c - \phi_s|^{-1.3}\). Thus, if we plot \({\eta_{\mathrm{eff}}}^{-1/1.3}\) against \(|\phi_c - \phi_s|\), we should be able to obtain a linear relationship. In FIG. 2, we plot \(\eta_{\mathrm{eff}}\) against \(\phi_s\) for all simulations and find that \(\eta_{\mathrm{eff}}\) always diverges at relatively high solid fractions (which are different for systems with different GSDs). The solid fraction, where a granular system is approaching infinite viscosity, indicates a critical solid fraction, \(\phi_c\). For each GSD, we can obtain a specific \(\phi_c\), and \({\eta_{\mathrm{eff}}}^{-1/1.3}\) and \(|\phi_c - \phi_s|\) are linearly related, as shown in the inset of FIG. 2.

Our above analyses indicate that the length scale \(d_{aw}\), which implies contact statistics, helps to unify the \(\mu_{\mathrm{eff}}\sim I_{cw}\) relationship but leaves a problem of \(\phi_c\) for us to solve. This implies that, for systems with different GSDs but exactly the same \(d_{aw}\), their stress-strain rate relationships behave similarly, but their solid fractions differ from each other. Since stresses in granular materials are strongly related to force networks (especially the strong force network), granular systems with different GSDs may have different amounts of rattler particles, which do not alter the strong force network, but occupy different amounts of volumes, resulting in different solid fractions.

We define particles with zero contact or only one contact as rattlers, while particles with more than one contact as non-rattlers. While calculating the contact number of each particle, \(N_c\), and the coordination number, \(Z_c\equiv\langle N_c\rangle\), of each simulation, we can also obtain the volumetric percentage of rattlers, \(V_{p01}/V_{p}\), and the coordination number of non-rattlers, \(Z_{c\geq 2}\). FIG. 3 (a) plots the histograms for a set of simulations with GSD of D3 and \(\eta_3 = 1.0\), where the legend provides the information of \(I_{cw}\). Increasing \(I_{cw}\) leads to a decrease in the number of particles with large contact numbers and increases the chances of finding particles with only one contact, indicating that increasing the inertial number results in more binary collisions. With such a wide range of particle sizes, a sheared granular system naturally has many “small" particles, which do not play a role in the force network. FIG. 3 (b) indicates that, for some cases, \(Z_c\) is even below 1.0, but for each set of simulations with the same GSD, the data collapse onto a \(Z_c\sim\phi_s\) master curve. However, it is difficult to figure out how changing the GSD alters \(Z_c\).

Neglecting particles with only zero and one contact number and changing the vertical axis into the average coordination number of particles with more than one contact, \(Z_{c\geq 2}\), in FIG. 3 (c), the relationship \(Z_{c\geq 2}\sim \phi_s\) shows power-law features with systematic changes introduced by changing the GSD. The inset of FIG. 3 (c) shows that dividing \(\phi_s\) by \(\phi_c\) helps to get a unified \(Z_{c\geq2}\sim\phi_s/\phi_c\) relationship, which is similar to the \(\eta_{\mathrm{eff}}\sim\phi_s/\phi_c\) relationship. We argue that it is mostly particles with \(N_c\leq 2\) that determine the effective viscosity of the system. All of these imply that changing the GSD leads to a different percentage of particles that participate in bearing forces in granular systems. However, FIG. 3 (d) contradicts our expectation and shows that changing the GSD does not change the volume percentage of particles much, as \(N_c\leq 1\).

We further plot the relationship between the normalized granular temperature \(T_g/(\dot{\gamma}d_{aw})^2\) and \(\phi_s\) in FIG. 3 (e), where \(T_g \equiv \langle (\vec{u} - \langle \vec{u} \rangle)^2 \rangle\) and \(\vec{u}\) is the particle velocity. Again, the critical solid fraction \(\phi_c\), which is heavily influenced by the GSD, plays a key role in the scaling of \(T_g/(\dot{\gamma}d_{aw})^2\sim\phi_s\), and \(T_g/(\dot{\gamma}d_{aw})^2\) scales with \(\sim |\phi_c-\phi_s|^{-0.6}\). We argue that \(\phi_c\), as a critical solid fraction, marks the transition from a system with finite viscosity to that with infinite viscosity. Physically, it has similar characteristics as random close packing fractions \(\phi_{RCP}\), which are affected by both \(\mu_p\) and GSD.

[17] and [18] suggested that \(\phi_{RCP}\) can be expressed as a function of the standard deviation of GSD, while [12] argued that they could draw a unified relation between \(\phi_{RCP}\) and a combination of skewness and polydispersity of GSD with the data obtained from simulations. However, we find that both methods do not work for our system. We adopted the idea of [12] and calculated both the skewness \(S_v\) and the polydispersity \(\delta_v\) of the normalized grain volume distribution (GVD), instead of the GSD, so that \[\begin{align} \tilde{V}_p^i = \left[V_p^i - \min(V_p)\right]/\left[ \max(V_p)-\min(V_p) \right]\;, \\ S_v = \langle \Delta\tilde{V}_p^3 \rangle/\left(\langle \Delta\tilde{V}_p^2 \rangle\right)^{3/2}\;,\;\delta_v = \sqrt{\langle \Delta\tilde{V}_p^2 \rangle}/\langle \tilde{V}_p \rangle\;, \end{align}\] where \(\tilde{V}_p\) is the normalized particle volume, \(\Delta\tilde{V}_p = \tilde{V}_p - \langle \tilde{V}_p\rangle\), \(\langle \tilde{V}_p^n \rangle = \int \tilde{V}_p^n\cdot P(\tilde{V}_p) d\tilde{V}_p\), and \(\langle \Delta\tilde{V}_p^{(n)} \rangle = \int \Delta\tilde{V}_p^n\cdot P(\tilde{V}_p) d\tilde{V}_p\). To obtain a universal relationship, we varied both the size ranges ([1, 10] cm, [2, 10] cm, and [3.16, 10] cm) and the particle frictional coefficient (\(\mu_p =\)​0.0001, 0.1, 0.3 and 0.9). With a combined dimensionless index of the GVD, \(I_{GVD}=\delta_v + S_v\delta_v^2\), we obtained an empirical relationship between \(\phi_c\) and \(I_{GVD}\) so that \[\label{eq:phic} \begin{align} \phi_c =& 1-(1-\phi_{c0})/\left[1 + \beta\ln{(1 + I_{GVD}/I_{g0})} \right]\;, \end{align}\tag{2}\] where \(\beta\approx 0.087\) and \(I_{g0} = 4.5\) are fitting constants, \(\phi_{c0}\) is the base \(\phi_c\) and varies with \(\mu_p\), and the functional form of Eq. 2 is inspired by the time evolution of the solid fraction of granular system under tapping compactions in [19]. FIG. 3 (f) suggested that we only need to adjust \(\phi_{c0}\) to get a reasonable \(\phi_c\sim I_{GVD}\) relationship for systems with different \(\mu_p\), and Eq. 2 , although being empirical, introduces a logarithmic dependence that ensures the eventual “flattening" of \(\phi_c\) with increasing \(I_{GVD}\). Ultimately, \(\phi_c\) should flatten out to a”maximum" packing fraction that is lower than or equal to 1 when \(I_{GVD}\rightarrow +\infty\), as also suggested in previous publications [12], [18], [20]. In Appendix 3, we plot the \(\phi_{RCP}\sim I_{GVD}\) relationship for data extracted from previous publications to show that Eq. 2 also works for predicting the random close packing fractions.

This Letter focuses on the rheology of granular systems with different grain size distributions and identifies a characteristic length scale that can be associated with the contact probability. We find that varying the GSD leads to the change in criticality behavior, i.e. the critical solid fraction \(\phi_c\) that marks the transition from a finite \(\eta_{\mathrm{eff}}\) to an infinite one. The analyses of \(Z_c\) and \(T_g\) verify the importance of determining \(\phi_c\) based on the GVD, and we have introduced a dimensionless index, \(I_{GVD}\), to successfully acquire a unified solution for granular systems with different particle gradations, which has the potential to be adopted to other athermal disordered systems with various polydispersities.

1 Simulations and stress calculations↩︎

The discrete element method (DEM) solves the momentum balance equation for both translational and rotational motions. The governing equations of DEM can be written in the form of the Newton’s second law, where the calculation of contact forces is of vital importance. In each simulation, the governing equations are integrated using a Verlet method [21], which is commonly used in molecular dynamics. Particles are spherical and follow a linear spring-dashpot contact model, where the normal contact force \(F_{ij}^n\) is calculated as a Hookean contact with energy damping, so that \(F_{ij}^n = -k_n\delta_{ij}^n-c_n\dot{\delta}_{ij}^n\), where \(k_n\) is the contact stiffness, \(c_n\) is the damping coefficient, and \(\delta_{ij}^n\) is the contact overlap. \(F_{ij}^n\) has a viscous term for energy dissipation, where \(c_n\) is calculated based on the coefficient of restitution, \(e_n\), so that \[\begin{align} c_n = 2 \xi_n \sqrt{m_{\mathrm{eff}}k_n}\;,\xi_n = -\ln{e_n} / \sqrt{\pi^2+(\ln{e_n})^2}\;, \end{align}\] where \(m_{\mathrm{eff}} = (1/m_i+1/m_j)^{-1}\) is the reduced mass and we set \(e_n=0.5\). The magnitude of the tangential contact force \(|k_t\delta_{ij}^t|\) cannot exceed its frictional limit, \(\mu_p|F_{ij}^n|\), where \(\mu_p\) is the friction coefficient of particles. The granular system is located within and sheared by two sawteeth plates at constant pressure and shear rate, as shown in the inset of FIG. 4 (b). At each timestep, the sheared granular system reaches a steady state when \(\phi_s\), \(\tau\), and \(\sigma_n\) become stable. Meanwhile, instead of measuring stresses from the upper and lower plates, we calculated the average stress tensor for each output timestep. The average stress tensor \(\sigma_{ij}\) can be obtained by \[\begin{align} \sigma_{ij} = (1/V)\sum_c^{N_{ct}}{f_c^i l_c^j}\;, \end{align}\] where \(V\) is the volume of the granular system, \(f_c^i\) is the \(i\)th component of a contact force vector (a summation of normal and tangential contact forces of the \(c\)th contact), and \(l_c^j\) is the \(j\)th component of a contact branch vector that links centroids of two contacting particles. FIG. 4 (a) presents typical time-lapse evolutions of \(\phi_s\), \(\tau\), and \(\sigma_n\), showing a clear steady state for determining the representative average values. Since we utilized sawteeth plates as upper and lower boundaries, the velocity profile across height is always linear [shown in FIG. 4 (b)], indicating a constant shear rate.

2 Particle size distributions↩︎

In this letter, we use two different types of particle size distributions, i.e., Fuller distribution (D1) and power-law distributions (D2 and D3). The cumulative distribution function (CDF) of D1 can be written as \[\label{eq:D1} \begin{align} C_1(d_p) = \left[ (d_p - d_{\min})/(d_{\max}-d_{\min}) \right]^{\eta_1}\;. \end{align}\tag{3}\] Then, a random particle diameter that follows this CDF can be obtained by imposing a uniformly distributed random number within \([0, 1]\), \(\mathrm{Rand}(0,1)\), on the left-hand side of Eq. 3 so that \[\begin{align} d_p = d_{\min} + [\mathrm{Rand}(0,1)]^{1/\eta_1}\cdot(d_{\max}-d_{\min}). \end{align}\] Considering the possible slopes of the PDF of power-law distributions, we take two different functional forms: \[\begin{align} &C_2(d_p) = (d_p/d_{\max})^{\eta_2}\;, \\ &d_p = d_{\max}\cdot[\mathrm{Rand}(0,1)]^{1/\eta_2}\;,\\ &C_3(d_p) = 1-(d_p/d_{\min})^{-\eta_3}\;,\\ &d_p = d_{\min}\cdot[1-\mathrm{Rand}(0,1)]^{-1/\eta_3}\;. \end{align}\] where we impose the \(d_p\in [d_{\min}, d_{\max}]\) condition every time we generate a random particle diameter. Figures 5 (a-f) show both the PDF and CDF of systems with different GSDs, and FIGs. 5 (g, h) show the rheological properties of granular systems with a different size range, which indicate that changing the particle size range would not change the rheological behavior, especially when we consider \(I_{cw}\) and eliminate the effect of \(\phi_c\).

3 \(\phi_{RCP}\) extracted from previous publications↩︎

We mentioned that previous research has focused on predicting the random close packing solid fraction, \(\phi_{RCP}\), which is closely related to \(\phi_c\) in our present work. However, we make clear that \(\phi_{RCP}\) is different from \(\phi_c\). The studies of \(\phi_{RCP}\) are mostly focused on frictionless systems. Furthermore, DEM simulations, where particles are in loading conditions, allow for the existence of overlaps among contacting particles, which will result in relatively larger solid fractions. To validate the \(\phi_c\sim I_{GVD}\) relationship, we turn to published works to also link \(\phi_{RCP}\) to \(I_{GVD}\). In FIG. 6, we plot the data extracted from both [22] with lognormal size distribution and [18] with \(\Gamma\) distribution. We converted their particle size distributions into particle volume distributions and obtained the corresponding statistical parameters to calculate \(I_{GVD}\). FIG. 6 shows that our empirical relationship of Eq. 2 also works for predicting \(\phi_{RCP}\) for frictionless spheres with a single fitting parameter of \(\phi_{c0}\approx 0.645\).

References↩︎