Effect of particle stiffness on contact dynamics and rheology in a dense granular flow

Dense granular flows have been well described by the Bagnold rheology, even when the particles are in the multibody contact regime and the coordination number is greater than 1. This is surprising, because the Bagnold law should be applicable only in the instantaneous collision regime, where the time between collisions is much larger than the period of a collision. Here, the effect of particle stiffness on rheology is examined. It is found that there is a rheological threshold between a particle stiffness of ${10}^4\text{–}{10}^5$ for the linear contact model and ${10}^5\text{–}{10}^6$ for the Hertzian contact model above which Bagnold rheology (stress proportional to square of the strain rate) is valid and below which there is a power-law rheology, where all components of the stress and the granular temperature are proportional to a power of the strain rate that is less then 2. The system is in the multibody contact regime at the rheological threshold. However, the contact energy per particle is less than the kinetic energy per particle above the rheological threshold, and it becomes larger than the kinetic energy per particle below the rheological threshold. The distribution functions for the interparticle forces and contact energies are also analyzed. The distribution functions are invariant with height, but they do depend on the contact model. The contact energy distribution functions are well fitted by Gamma distributions. There is a transition in the shape of the distribution function as the particle stiffness is decreased from ${10}^7$ to ${10}^6$ for the linear model and ${10}^8$ to ${10}^7$ for the Hertzian model, when the contact number exceeds 1. Thus, the transition in the distribution function correlates to the contact regime threshold from the binary to multibody contact regime, and is clearly different from the rheological threshold. An order-disorder transition has recently been reported in dense granular flows. The Bagnold rheology applies for both the ordered and disordered states, even though the rheological constants differ by orders of magnitude. The effect of particle stiffness on the order-disorder transition is examined here. It is found that when the particle stiffness is above the rheological threshold, there is an order-disorder transition as the base roughness is increased. The order-disorder transition disappears after the crossover to the soft-particle regime when the particle stiffness is decreased below the rheological threshold, indicating that the transition is a hard-particle phenomenon.

Figure 1

The variation in the volume fraction $\varphi$ [(a) and (c)] and the coordination number $Z$ [(b) and (d)] with height $z$ for different values of the particle stiffness $k_n$ for the linear contact model [(a) and (b)] and the Hertzian contact model [(c) and (d)]. The parameters for the different symbols are provided in Table 1.

Figure 2

The variation of the strain rate with height for different values of the particle stiffness $k_n$ for the linear contact model (a) and the Hertzian contact model (b). The parameters for the different symbols are provided in Table 1.

Figure 3

The variation of the normal stresses ${\sigma }_{xx}$ and ${\sigma }_{zz}$ as a function of $(d{u}_x/dz)$ for different values of the particle stiffness $k_n$ for the linear contact model (a) and the Hertzian contact model (b). The dashed lines are power-law fits [Eq. (2)], the dotted line on the left shows ${\sigma }_{ij}\propto {(d{u}_x/dz)}^2$, and the dotted line on the right shows ${\sigma }_{ij}\propto (d{u}_x/dz)$. The parameters for the different symbols are provided in Table 1, and the values of ${\eta }_{\sigma }$ are shown as a function of the particle stiffness later in Fig. 5.

Figure 4

The variation of the kinetic energy per particle as a function of $(d{u}_x/dz)$ for different values of the particle stiffness $k_n$ for the linear contact model (a) and the Hertzian contact model (b). The dashed lines show the fits to Eq. (4), the dotted line on the left shows $K\propto {(d{u}_x/dz)}^2$, and the dotted line on the right shows $K\propto {(d{u}_x/dz)}^{3/2}$. The parameters for the different symbols are provided in Table 1, and the values of ${\eta }_T$ are shown as a function of the particle stiffness later in Fig. 5.

Figure 5

The height-averaged coordination number $Z$ ($\circ$), the exponents ${\eta }_{\sigma }$ in Eq. (2) $\left(▵\right)$, ${\eta }_T$ in Eq. (4) ($▿$), ${\eta }_{N2}$ in Eq. (3) ($◃$), and ${\eta }_{\gamma }$ in Eq. (1) ($▹$), as a function of the particle stiffness $k_n$ for the linear contact model (a) and Hertzian contact model (b). The dotted horizontal lines show the exponent values 0.5, 1.0, 1.5, and 2.0 for reference.

Figure 6

The variation of the magnitude of the average resultant force on a particle $F_R$ defined in Eq. (5) for the linear contact model (a) and the Hertzian contact model (b). The dashed lines show linear fits for the data. The parameter values for the different symbols are shown in Table 1.

Figure 7

The height-averaged values of $F_R$ defined in Eq. (5) ($\circ$) and $F_M$ defined in Eq. (6) ($▵$) referenced to the left $y$ axis on a log scale, and the ratio $(F_M/F_R)(\diamond )$ referenced to the right $y$ axis on a linear scale, as a function of the particle stiffness for the linear contact model (a) and Hertzian contact model (b).

Figure 8

The variation of the kinetic energy per particle $T$ (open symbols) and contact energy per particle $E_c$ (filled symbols) as a function of height. The parameters for the different symbols are shown in Table 1. The dashed lines show the power-law fits [Eq. (8)].

Figure 9

The exponent ${\eta }_{Tz}$ in Eq. (8) $\left(▵\right)$, ${\eta }_{Ez}$ in Eq. (8) ($\circ$), and the exponent ${\eta }_{ET}$ in Eq. (9) ($▿$). The predicted exponent ${\eta }_{TE}$, from the scaling relationship between the kinetic and contact energies, ${\eta }_{ET}=1.5$ for the linear contact model and ${\eta }_{ET}=1.4$ for the Hertzian contact model, is shown by the dotted lines.

Figure 10

The height-averaged kinetic energy per particle ($\circ$), contact energy per particle $\left(▵\right)$, and the ratio of the contact and kinetic energies $▿$ a function of the particle stiffness for the linear contact model (a) and the Hertzian contact model (b). The dotted line has a slope $-1/2$ in panel (a) and $-2/5$ in panel (b).

Figure 11

The distribution function for the contact force magnitude at a height of 6 particle diameters from the base ($\circ$), 26 particle diameters from the base ($▵$), 46 particle diameters from the base ($▿$), and 62 particle diameters from the base ($\diamond$) for the linear contact model with particle stiffness $k_n={10}^3$ (a) and ${10}^8$ (b), and for the Hertzian particle contact model for particle stiffness ${10}^4$ (c) and ${10}^8$ (d).

Figure 12

The distribution function for the contact energy at a height of 6 particle diameters from the base ($\circ$), 26 particle diameters from the base ($▵$), 46 particle diameters from the base ($▿$), and 62 particle diameters from the base ($\diamond$) for the linear contact model with particle stiffness $k_n={10}^3$ (a) and ${10}^8$ (b), and for the Hertzian particle contact model for particle stiffness ${10}^4$ (c) and ${10}^8$ (d).

Figure 13

The height-averaged distribution function for the force magnitude (a) and the contact energy (b) for the linear contact model for particle stiffness $k_n=5\times {10}^2$ ($\circ$), $k_n={10}^6$ ($▵$), $k_n={10}^7$ ($▿$), and $k_n={10}^8$ ($\diamond$). In panel (a), the dotted line is the exponential distribution. In panel (b), the dashed line is the Gamma distribution $P_{\mathrm{\Gamma }}(0.7,(E_c/{\overline{E}}_c))$, and the dotted line is the Gamma distribution $P_{\mathrm{\Gamma }}(0.4,(E_c/{\overline{E}}_c))$.

Figure 14

The height-averaged distribution function for the force magnitude (a) and the contact energy (b) for the Hertzian contact model for particle stiffness $k_n=2\times {10}^3$ ($\circ$), $k_n={10}^7$ ($▵$), $k_n={10}^8$ ($▿$), and $k_n={10}^{10}$ ($\diamond$). In panel (a), the dotted line is the exponential distribution. In panel (b), the dashed line is the Gamma distribution $P_{\mathrm{\Gamma }}(0.6,(E_c/{\overline{E}}_c))$, and the dotted line is the Gamma distribution $P_{\mathrm{\Gamma }}(0.55,(E_c/{\overline{E}}_c))$.

Figure 15

The variation with height of the number density (a), the $q_6$ order parameter (b), and the mean velocity (c) for the linear model with dimensionless particle stiffness $k_n={10}^7$, and for a hexagonally ordered frozen-particle base with base particle diameter 0.5 $[\circ$ in (b) and (c)], 1.34 [thin line in (a), $▵$ in (b) and (c)], 1.35 [thick line in (a), $▿$ in (b) and (c)], and 2.00 $[\diamond$ in (b) and (c)].

Figure 16

The variation with height of the number density [(a) and (b)], the $q_6$ order parameter [(c) and (d)], and the mean velocity [(e) and (f)] for the linear model [(a), (c), and (e)] and the Hertzian model [(b), (d), and (f)]. In panels (a), (c), and (e) particle stiffness $k_n={10}^7[\circ$ in (c) and (e)], $k_n=6\times {10}^4$ [thin line in (a), $▵$ in (c) and (e)], $k_n=2\times {10}^4$ [thick line in (a), $▿$ in (c) and (e)], and $5\times {10}^2[\diamond$ in (c) and (e)]. In panels (b), (d), and (f) particle stiffness $k_n={10}^8[\circ$ in (d) and (f)], $k_n=6\times {10}^5$ [thin line in (b), $▵$ in (d) and (f)], $k_n=2\times {10}^5$ [thick line in (b), $▿$ in (d) and (f)], and $2\times {10}^3[\diamond$ in (d) and (f)].