Crossover from quasi-static to dense flow regime in compressed frictional granular media

Figure 2 shows the macroscopic response of a granular sample loaded using ${\delta \sigma_1^{t_r}=1.10^{-6} \sigma_3}$. The inertial number I quantifies the ratio between inertial forces and imposed forces and is defined as $I = \dot{\epsilon}_1 \sqrt{\overline{m}/\sigma_3}$ [16], where $\overline{m}$ is the average grain mass. Initially, I is of the order of $10^{-6}\text{-}10^{-5}$ (see fig. 2(a)). As τ increases towards $2\sigma_3$, while undergoing fluctuations associated to plastic events, I remains lower than $10^{-4}$, which is often considered an upper bound for quasi-static conditions [16,21]. Hence, in the range $0 \le \tau < 2\sigma_3$, the sample undergoes quasi-static deformation. At values of τ larger than $\tau_c \approx 2\sigma_3$, I increases by several orders of magnitudes, reaching values of about $10^{-3}\text{-}10^{-2}$. This drastic change in I-values defines the macroscopic instability and the transition to a dense flow regime, where inertial forces become significant. Figure 2(b) shows the surface change of the granular assembly $\Delta S/S_0=\frac{(h_x h_y)_{1}-(h_x h_y)_{0}}{(h_x h_y)_{0}}$, where 0 and 1 denote the initial and current configuration. Sample contraction is observed in the early stages of mechanical testing and essentially results from contacts elasticity, as it would not be observed in the limit of rigid grains [22]. Sample dilation irreversibly occurs after the maximum contraction and the dense flow transition observed at $\tau_c\approx 2\sigma_3$ occurs at packing fractions larger than the initial one.

Zoom In Zoom Out Reset image size

The critical shear stress value $\tau_c \approx 2\sigma_3$ materializes the softening of the granular assembly, as it corresponds to a macroscopic friction coefficient of $\mu_{\textit{macro}}\approx 0.5$ that is significantly smaller than the friction $\mu_{\textit{micro}}=1$ considered at the grain scale. While macroscopic hardening can be associated to the non-zero macroscopic friction obtained for an assembly of frictionless circular particles submitted to a similar compressional loading [3], the assemblies of frictional particles herein considered show a reverse behaviour, characterized by macroscopic softening. In the present configuration, collective effects and structural rearrangements generate weak zones whose spatial extent that are much larger than the particle size grows towards the macroscopic instability. Inelastic deformation and flow concentrate within these weak zones, resulting in a macroscopic softening. The quantification of these effects towards the macroscopic instability through the analysis of stress and strain fields is the backbone of this manuscript.

We first characterize the stress field by quantifying the spatial extent of stress concentration regions using a coarse-graining analysis [23,24]. An average shear stress rate $\langle\dot{\tau}\rangle$ (see footnote ¹ ) is computed at different stages of mechanical testing (cf. color dots in fig. 2) over a time window T and over a broad range of spatial scales L, ranging from the micro-scale corresponding to mesh size (see fig. 3 (left)), to the macro-scale corresponding to sample size. Subsystems of the granular assembly are selected by defining square boxes of size W (see fig. 3 (right)). Grains are considered to belong to a given subsystem if their mass centers are included in the box of size W. For a given box of size W, the actual area spanned by the mesh associated to the selected grains assembly slightly differs from the area of the box W². This difference is due to the irregularity of the contour of the mesh at the boundary of the domain. Thus, for a set of $N_{\textit{box}}$ boxes of size W, we define the average scale L as $L = \frac{1}{N_{\textit{box}}} \sum_{k=1}^{N_{\textit{box}}} L_k$, where L_k is the scale associated to the box number k, computed as the square root of the sum of mesh element surfaces. The smallest subsystem of average size $L = 1.1 D_{\textit{max}}$ shows a maximum standard deviation of L_k-values of $0.2 D_{\textit{max}}$ for the Voronoi triangulation used to compute stresses and of $0.13 D_{\textit{max}}$ for the Delaunay triangulation used to compute strains (see fig. 3). No overlap of scales L_k occurs through successive values of W. For a given assembly of grains included in the box number k, the stress tensor is computed over all contacts carried by each grain using [25], i.e. writing

$$\begin{equation} \sigma_{ij}^k = \frac{1}{L_k^2} \sum_{g_k} \sum_{c_k} (r_i^{c_k} - r_i^{g_k}) f_j^{g_k c_k}, \end{equation} \tag{ 1 }$$

where f_j^g_kc_k is the j-th component of the contact force exerted on grain g_k at contact c_k, r_i^c_k is the i-th component of the position vector of c_k, and r_i^g_k is the i-th component of the position vector of the mass center of grain g_k. The stress rate tensor is obtained by differentiating the stress tensor components calculated at two successive configurations. These successive configurations are selected in time and separated one from the other by a time interval $T = \sqrt{N_g} \times 100 \times t_r$ that corresponds to the travel time of elastic waves through the whole sample [17].

Zoom In Zoom Out Reset image size

Figure 4 (top) shows the results of the multiscale analysis performed on the stress rate field from the early stages of biaxial testing up to $\tau_c \approx 2\sigma_3$. Curves are selected with respect to the control parameter Δ defined as $\Delta = \frac{\tau_c - \tau}{\tau_c}$. At the early stages of macroscopic deformation, i.e. at large values of Δ, and at small spatial scales, $\langle\dot \tau\rangle$ decreases as L increases. Moreover, for L-values larger than a crossover scale $l^*_\tau$, $\langle\dot \tau\rangle$ remains constant. We conclude that associated shear stress rate fields are heterogeneous for $l \ll l^*_\tau$ and homogeneous for $l\gg l^*_\tau$. Associated fields provided on the right-hand side of fig. 4 illustrate the connection between $l^*_\tau$ and the size of spatial zones of stress concentration. For $\Delta = 0.42$, the stress rate field shows a roughly homogeneous pattern at large spatial scales, i.e. scales larger than $l^*_\tau \approx 3 D_{\textit{max}}$ in that case, while the region of stress localization observed on the top-left part of the snapshot exhibits a spatial extent of the order of $l^*_\tau \approx 3 D_{\textit{max}}$ that corresponds to the crossover scale of the associated curve. We interpret $l^*_\tau$ as the correlation length [24]. As macroscopic deformation proceeds, $l^*_\tau$ grows until reaching the size of the system at $\tau_c$. At that point, a power law scaling $\langle\dot \tau\rangle \sim L^{-\rho_\tau}$ is observed over the entire scale range, with $\rho_\tau = 0.38$. We verified using various sample sizes that the cut-off remaining on the scaling at $\Delta \rightarrow 0$ is a finite-size effect [17]. At $\tau = \tau_c$, large and strongly localized structures characterize the shear stress rate field (see snapshot at $\Delta = 0.005$). These scalings suggest a progressive structuring of the stress rate field as approaching the transition to the dense flow regime. We report an associated divergence of the correlation length $l^*_\tau$ by performing a collapse analysis (inset of fig. 4 (top)) and formulating that $l^* = l^*_\tau$ diverges as

$$\begin{equation} l^* \sim \frac{L_s^{\delta}}{\Delta^{\nu} L_s^{\delta} + C}, \end{equation} \tag{ 2 }$$

where L_s is the square root of the sample area and $\nu = \nu_\tau = 1.3 \pm 0.1$ is the exponent of divergence. The parameters δ and C characterize the finite-size effect [17]. The divergence of $l_\tau^* \sim \Delta^{-\nu_\tau}$ at $\Delta \rightarrow 0$ is confirmed from a similar analysis performed on other moments $\langle\dot \tau^q\rangle$ of the shear rate [17], revealing the multi-fractality of the shear rate field at the critical point. These particular features of the shear stress rate field are only observed at the specific timescale $t = T$, that corresponds to the travel time of elastic waves. At increasing timescales, i.e. for t > T, the departure from the power law scaling observed in fig. 4 (bottom) denotes a progressive homogenization of the shear stress rate field, materialized by the decrease of $l^*_\tau$ as t increases (see also the associated snapshots). We interpret the multi-scale properties of the shear stress rate field observed at $t = T$ to be associated to stress concentration events occurring simultaneously to the propagation of elastic information throughout the sample, which operates by means of elastic waves. The loss of scaling properties reported beyond this time, i.e. for t > T, implies that the multi-scale localization structures observed at $t = T$ have little memory, limited to the travel time of an elastic wave. In other words, the loss of scaling properties at t > T is explained by the superposition of several stress concentration events operating at $t = T$, but uncorrelated in time one with the successive other.

Zoom In Zoom Out Reset image size

We then investigate whether similar features characterize the shear strain field. A Delaunay triangulation is performed on the grain centers of the non-rattler-grains set (fig. 3). The partial derivatives are computed at the mesh scale as $\epsilon_{ij} = 1/2 (\partial u_i/\partial x_j + \partial u_j/\partial x_i)$, where $(u_1,u_2)$ and $(x_1,x_2)$ are, respectively, the incremental displacements and spatial coordinates of grain centers. By averaging the partial derivatives, an average value of the shear strain rate is computed at increasing spatial scales L using the coarse-graining analysis exposed previously. Using the constant timescale t = T to compute displacement rates, the correlation length associated to the shear strain rate field does not diverge as approaching the transition to the dense flow regime (not shown). Thus, the structure of the strain rate field does not form simultaneously in time with the stress rate field, which would be the case if one considered only the elastic component of the strain. Here, for $N_g=10000$, approximately 10⁴ stress increments are prescribed during the propagation time of an elastic wave throughout the sample. Hence, a multitude of contacts, in our case about 5% of the whole contact network, are sliding although elastic information does not have time to travel across the entire sample. However, a progressive structuring of the shear strain field is observed when computing the scaling of the incremental shear strain $\langle\delta \gamma\rangle$ for constant macroscopic deformation windows $\delta \epsilon_1=\delta\epsilon_p=1.10^{-5}$. The divergence of the correlation length $l_\gamma^*$ as approaching the transition to the dense flow regime is reported on the incremental shear strain field in fig. 5. This divergence is similar to the one observed on the shear stress rate field, as we find $\nu_\gamma=\nu_\tau=1.3$ (eq. (2)). For larger macroscopic deformation windows $\delta \epsilon_1>\delta\epsilon_p$, the deformation field no longer exhibits multi-scale properties and $l_\gamma^*$ continuously decreases as $\delta \epsilon_1$ increases. This observation is in agreement with the study of [14], who reported a narrowing of the fluctuation velocity distributions at increasing timescales. As the material softens when τ increases (see fig. 2(a)), the time of integration associated to the constant deformation window $\delta \epsilon_1=\delta \epsilon_p$ decreases when approaching the dense flow transition. This time can be either smaller or larger than the travel time T of an elastic wave, depending on the imposed loading rate $\delta \sigma_1^{t_r}$. All the presented scalings are indifferently obtained for loading rate values $\delta \sigma_1^{t_r}$ smaller than $\delta \sigma_1^{t_r}=1.10^{-6}$. For larger values of $\delta \sigma_1^{t_r}$, the divergence of the correlation length is no longer observed [17]. In that case, the mechanical behavior of the granular assembly is related to a dense flow regime from the initial stages of deformation, attesting that the scalings reported in this study are associated to the transition from quasi-static to dense flow regimes of deformation. Interestingly, $\delta \epsilon_1=\delta \epsilon_p$ remains equal to $1.10^{-5}$ whatever the loading rate considered, suggesting that a given, constant, amount of plastic activity has to operate in order to observe multi-scale properties within the incremental shear strain field.

Zoom In Zoom Out Reset image size

To better understand the role of plastic activity in setting the characteristic value $\delta\epsilon_p$, we perform a four-point dynamic susceptibility $\chi_4$ analysis on the inter-particles contact network [26,27]. From an "initial" contact configuration selected at an axial deformation $\epsilon_1^{\textit{init}}$, we compute the self-overlap order parameter $Q_{\epsilon_1^{\textit{init}}}(\delta \epsilon_1)=\frac{1}{N_c}\sum_{i=1}^{N_c}w_i$, where N_c is the number of the non-sliding contacts of the initial configuration and w_i is a step-function cutoff that equals 1 if no sliding event occurred at contact i over the deformation window $\epsilon_1^{\textit{init}}\rightarrow \epsilon_1^{\textit{init}}+\delta \epsilon_1$, and 0 otherwise. The first two moments $Q(\delta \epsilon_1)=\langle Q_{\epsilon_1^{\textit{init}}}(\delta \epsilon_1) \rangle$ and $\chi_4(\delta \epsilon_1)=N_c[\langle Q_{\epsilon_1^{\textit{init}}}(\delta \epsilon_1)^2 \rangle-\langle Q_{\epsilon_1^{\textit{init}}}(\delta \epsilon_1)\rangle^2]$ of $Q_{\epsilon_1^{\textit{init}}}(\delta \epsilon_1)$ (calculated from sample-to-sample fluctuations) are computed for initial configurations selected at the previously used Δ-values, i.e. within the quasi-static region $(\tau < \tau_c)$. This analysis evaluates the number and the associated spatial heterogeneity of sliding events that occur as the axial deformation increases. As no significant variations of $Q(\delta \epsilon_1)$ and $\chi_4(\delta \epsilon_1)$ are observed with respect to Δ, fig. 6 shows the average functions $\langle Q(\delta \epsilon_1)\rangle$ and $\langle\chi_4(\delta \epsilon_1)\rangle$ computed over all the values of Δ. Major and minor force networks, defined by selecting contact forces larger and smaller than the average force, are considered separately. By construction, $\langle Q(\delta \epsilon_1)\rangle$ is initially equal to 1. As $\delta \epsilon_1$ increases, $\langle Q(\delta \epsilon_1)\rangle$ decreases but never reaches 0, as about 35% (respectively, 55%) of the contacts of the minor (respectively, major) network never slide. These non-sliding contacts associated to the remaining rigid bodies illustrate the localization of permanent deformation. At low values of $\delta \epsilon_1$, $\langle\chi_4(\delta \epsilon_1)\rangle$ increases as $\langle\chi_4(\delta \epsilon_1)\rangle \sim \delta \epsilon_1^\beta$ with $\beta = 1.8$. This increase results from the nucleation of spatially correlated events of contact sliding. As $\delta \epsilon_1$ exceeds the threshold value $\delta \epsilon_p = 1.10^{-5}$, $\langle\chi_4(\delta \epsilon_1)\rangle$ saturates. Thus, at $\delta \epsilon_p = 1.10^{-5}$, i.e. when multi-scale properties are observed on the shear deformation field (see fig. 5), all the spatially correlated contacts located close to the Coulomb criterion, i.e. susceptible to slide, have been destabilized. At this stage, only 4% (respectively, 3%) of contacts have slid in the minor (respectively major) network.

Zoom In Zoom Out Reset image size