Rheology of three dimensional granular chute flows at large inertial numbers

We perform DEM simulations for inelastic spheres having restitution coefficient $e_{n}=0.5$ spanning a wide range of inclination angles $\theta=22^{\circ}$ to $38^{\circ}$ for the initial flowing layer thickness $L_{y}=40d$. The average kinetic energy of the particles at these inclinations are shown in Fig. 3(a) and Fig. 3(b). The average value of the kinetic energy keeps increasing with time and eventually becomes constant, indicating that the system is able to achieve a steady state at these inclinations. The average kinetic energy of the particles, shown in Fig. 3(a) and Fig. 3(b), becomes constant for all the inclination angles considered. This confirms that steady, fully developed flows are possible over a larger range of inclinations for $e_{n}=0.5$. The average kinetic energy of the particles at steady state is plotted with the inclination angle in Fig. 3(c). As expected, the kinetic energy increases with inclination. A linear increase in the average kinetic energy is observed with the inclination angle $\theta$. Fig. 3(d) shows the time taken to achieve $95\%$ of steady state kinetic energy, denoted as $T_{95}$, for different inclinations. The time taken to reach a steady state increases with the inclination angle and a sharp change is observable after $\theta=32^{\circ}$. Note that the kinetic energy in the $y$ axis is non-dimensionalized using $mgd$ where $m$ is the mass, $d$ is the particle diameter and $g$ is the acceleration due to gravity.

From the DEM simulations data obtained from two different flowing layer thicknesses $L_{y}=40d$ and $L_{y}=20d$, we report the variation of the effective friction coefficient $\mu$, solids fraction $\phi$, the ratio of first normal stress difference to pressure $N_{1}/P$ and the ratio of second normal stress difference to pressure $N_{2}/P$ with inertial number $I$ for $e_{n}=0.5$ in Fig. 4. The black circles correspond to the DEM simulations starting from a settled layer height of $L_{y}=40d$ and the red circles correspond to the DEM simulations starting from a settled layer height of $L_{y}=20d$. The data for the two cases do not differ significantly from each other and can be described using a single curve. The black solid line represents the best fit to the data accounting for both layer heights. The effective friction coefficient obtained from the simulations (shown using black and red symbols in Fig. 4(a)) increases up to $I\approx 1.8$ and shows a saturating behavior at large inertial numbers. There is a slight decrease in the effective friction coefficient $\mu$ at large inertial numbers $I$, however, the decrease is not prominent, and the inclusion of $\mu-I$ data for inclination angles $\theta>38^{\circ}$ is needed to confirm the decrease in the effective friction coefficient at large inertial numbers. Given the sharp increase in the computational time for higher inclinations, we proceed with the analysis of the existing data.

Fitting the 4-parameter MK model using Eq. (6) to the simulation data shown in Fig. 4(a) indicates that the MK model (shown using solid lines) is able to describe the $\mu-I$ variation very well. Note that the simulation data for the inclination angle $\theta=22^{\circ}-38^{\circ}$ is considered to obtain the model parameters. The MK model parameters are reported in Table 1. The variation of solids fraction is plotted with the inertial number in Fig. 4(b). The solid line shown in Fig. 4(b) is obtained by fitting the power law variation of $\phi$ with $I$ (Eq. (7)) to the simulation data. The dilatancy law model parameters are reported in Table 2. In Fig. 4(c), we report the variation of the ratio of the first normal stress difference to the pressure ratio $N_{1}/P$ with the inertial number $I$. The variation for the first normal stress difference to the pressure ratio $N_{1}/P$ with the inertial number $I$ is described well by a linear relation of the form $N_{1}/P=f_{1}(I)=AI+B$. The first normal stress difference law model parameters are reported in Table 3. We also report the variation of the ratio of the second normal stress difference to the pressure ratio $N_{2}/P$ with the inertial number $I$. The variation for the second normal stress difference to the pressure ratio $N_{2}/P$ with the inertial number $I$ is described well by a quadratic relation of the form $N_{2}/P=f_{2}(I)=FI^{2}+GI+H$. The second normal stress difference law model parameters are reported in Table 4. A comparison of Fig. 4(c) and Fig. 4(d) shows that while $N_{1}/P$ is zero at low inertial numbers, $N_{2}/P$ is found to be positive even for the quasistatic systems with $I\sim 0$. In addition, $N_{2}/P>N_{1}/P$ for $I\leq 1.25$. For $I>1.25$, however, $N_{1}/P$ becomes larger than $N_{2}/P$. This transition of $N_{2}/P-N_{1}/P$ from positive to negative at $I\sim 1.25$ needs to be explored further.

Using the rheological model parameters given in Table 1-4, we predict various flow properties of interest. These theoretical predictions are reported as solid lines in Figs. 5 and 6. The velocity profiles predicted using the theory are found to be in excellent agreement with the simulation results for inclination angles ranging from $\theta=22^{\circ}-32^{\circ}$ (see Fig. 5(a)). The slip velocity, which appears to be small for these inclinations is taken from the simulation data. The predicted values of the solids fraction $\phi$ are also found to be in very good agreement for inclination angles $\theta=22^{\circ}$ to $\theta=32^{\circ}$ as shown in Fig. 5(b). The theoretical inertial numbers are also found to be in good agreement with the inertial numbers obtained from DEM simulations (see Fig.6(c)).

Figs.6(a-d) shows the comparison between DEM simulations and the theoretical predictions for the pressure $P$, shear stress $\tau_{yx}$, first normal stress difference $N_{1}$, and the second normal stress difference $N_{2}$ for the range of inclinations $\theta=22^{\circ}$ to $\theta=32^{\circ}$. The theoretical predictions using the rheological parameters are in very good agreement with the DEM simulation results. Comparison of $N_{1}$ and $N_{2}$ shows that $N_{2}$ is higher than $N_{1}$ for these range of inclinations.

We next use the rheological model parameters given in Table 1-4 to predict various flow properties at large inclination angles ranging from $\theta=34^{\circ}-38^{\circ}$. These theoretical predictions are reported as solid lines in Figs. 7 and 8. The velocity profiles predicted using the theory are found to be in good agreement with the simulation results even at large inclination angles ranging from $\theta=34^{\circ}-38^{\circ}$ (see Fig. 7(a)). We use the slip velocity obtained from DEM simulations in the theoretical predictions. The predicted values of the solids fraction $\phi$ are also found to be in very good agreement for inclination angles $\theta=34^{\circ}$ to $\theta=38^{\circ}$ as shown in Fig. 7(b). In addition, the theoretical inertial numbers are also found to be in good agreement with the inertial numbers obtained from DEM simulations (see Fig.8(c)) except the largest inclination angle $\theta=38^{\circ}$ where the theory predicts the inertial number $I=1.76$ in comparison to the inertial number $I=1.9$ obtained from DEM simulations. Figs. 8(a)-(d) show the comparison between DEM simulations and the theoretical predictions for the pressure $P$, shear stress $\tau_{yx}$, first normal stress difference $N_{1}$, and the second normal stress difference $N_{2}$ for the large inclinations angles ranging from $\theta=34^{\circ}$ to $\theta=38^{\circ}$. The theoretical predictions using the rheological parameters are in very good agreement with the DEM simulation results.

Many industrial granular flows in conveyors and chutes have flowing layer thickness in the range $10d-20d$. It is well known that the inertial number $I$ and solids fraction $\phi$ are nearly constant in the bulk region of the flowing layer of chute flow for any given inclination angle $\theta$. The bulk region is defined as the region that is away from the base of the chute and the free surface. The bulk region is very important in order to estimate the rheological properties. Lack of bulk region in the case of low flowing layer thickness will lead to fluctuations of flow properties in the bulk and will lead to significant errors in estimating the rheological parameters. Due to the lack of bulk region in the case of low flowing layer thickness $L_{y}=20d$, we have performed DEM simulations starting from a settled flowing layer height $L_{y}=40d$. We will now use the rheological model parameters given in Table 1-4 obtained from the flowing layer thickness $L_{y}=40d$ to predict various flow properties for a wider range of inclination angles ranging from $\theta=22^{\circ}-38^{\circ}$ for low flowing layer thickness $L_{y}=20d$ which is the initial settled flowing layer thickness.

These theoretical predictions are reported as solid lines in Figs. 9 and 10. The velocity profiles predicted using the theory are found to be in good agreement with the simulation results for inclination angles ranging from $\theta=22^{\circ}-32^{\circ}$ (see Fig. 9(a)). Figs. 9(b) and Fig. 9(c) show that the predicted values of the solids fraction $\phi$ and inertial number $I$ are also found to be in very good agreement for inclination angles $\theta=22^{\circ}$ to $\theta=32^{\circ}$. Figs.10(a-d) shows the comparison between DEM simulations and the theoretical predictions for the pressure $P$, shear stress $\tau_{yx}$, first normal stress difference $N_{1}$, and the second normal stress difference $N_{2}$ for inclinations angles ranging from $\theta=22^{\circ}$ to $\theta=32^{\circ}$. Again, the theoretical predictions using the rheological parameters are in very good agreement with the DEM simulation results. Next, the theoretical predictions for the various flow properties of interest have been compared with the DEM simulation results for higher inclination angles in the range of $\theta=34^{\circ}$ to $\theta=38^{\circ}$. The predictions from the rheological parameters obtained using the simulations data of $L_{y}=40d$ layer are in very good agreement with the DEM simulation results for low flowing layer thickness $L_{y}=20d$ suggesting that rheological parameters obtained from large flowing layer thickness data can be used to predict the flow properties for low flowing layer thickness as well (shown in Figs.11-12).

Figs. 13(a)-13(d) summarise the results over the entire range of inclinations investigated in this study and show the variation of the inertial number $I_{bulk}$, solids fraction $\phi_{bulk}$, height of the bulk region of the flowing layer $h_{bulk}$ and the slip velocity at the base $v_{slip}$ as a function of inclination angle $\theta$ for restitution coefficient $e_{n}=0.5$. The slope of the curve for the lower angles differs significantly from that for the higher angles. This slope contrast is evident in the case of slip velocity shown in Fig. 13 for different flow regimes: dense flows with negligible slip velocity at lower angles and rapid, dilute flows with large slip velocities at higher angles. The steady-state solids fraction $\phi$ is found to decrease with the increase in the inclination angle $\theta$ and can be seen in Fig. 13(b). The decrease in the solids fraction is followed by the dilation of the bulk flowing layer thickness $h_{bulk}=h_{initial}\times\phi_{max}/\phi_{steady}$ which increases with the inclination angle. This increase in the bulk flowing layer thickness with the inclination angle is shown in Figs. 13(c). The slip velocity at the base of the chute is shown in Figs. 13(d) which increases with the increase in inclination angle and changes its slope gradually beyond inclination angle $\theta=30^{\circ}$.

In this section, we will discuss the effect of the first and the second normal stress differences on the flow properties of granular materials for the range of inclination angle $\theta=22^{\circ}$ to $\theta=32^{\circ}$. Figure 14(a) shows the theoretical predictions of the inertial number using the JFP model with and without considering the effect of the first and the second normal stress differences i.e. by assuming $f_{1}(I)=N_{1}/P=0$ and $f_{2}(I)=N_{2}/P=0$ in the latter case. Symbols represent the DEM simulation results. The dashed lines correspond to the theoretical predictions of the inertial number $I$ using the JFP model where the normal stress differences have been ignored whereas the solid lines correspond to the theoretical predictions of the inertial number $I$ using the JFP model where the normal stress differences have been considered. It is evident from figure 14(a) that the theoretical inertial number starts to deviate from the DEM simulation results as we increase the inclination angle $\theta$. Figure 14(b) shows the theoretical predictions of the inertial number using the modified rheological model with and without considering the effect of the first and the second normal state differences. The latter case assumes $f_{1}(I)=N_{1}/P=0$ and $f_{2}(I)=N_{2}/P=0$. As before, symbols represent the DEM simulation results. The dashed lines correspond to the theoretical predictions of the inertial number $I$ using the modified rheological model where the normal stress differences have been ignored whereas the solid lines correspond to the theoretical predictions of the inertial number $I$ using the modified rheological model where the normal stress differences have been considered. The deviation in the theoretical predictions of the inertial number $I$ from the DEM simulation results are observed in the case of the modified rheologial model as well and this deviation is more significant at higher inclinations. These results suggest that the role of the normal stress differences becomes crucial not only for large inclinations but also for small inclinations as well and one must consider these effects in order to accurately predict the flow properties. Neglecting the effect of the normal stress differences will lead to significant errors in the theoretical predictions.

The steady-state oscillations in the bulk layer height $h_{bulk}$ are shown in Fig. 15(a) over a time period of $200$ time units starting from a time instant $t_{0}$. Oscillations in the center of mass $y_{com}$ are shown in Fig. 15(b). These height measurements are done after the flow has achieved steady kinetic energy (hence $t_{0}$ varies with the inclination angle $\theta$). The oscillations in the bulk layer height $h_{bulk}$ as well as center of mass $y_{com}$ keep increasing with the inclination angle. While the difference between maximum and minimum bulk height is less than the particle diameter for inclinations $\theta\leq 30^{\circ}$, this variation becomes $\Delta h_{bulk}\approx 4d$ for $\theta=38^{\circ}$ case (Fig. 15(a)). These oscillations indicate that the role of density and height variations become important and cannot be ignored for granular flows at inertial numbers comparable to unity.

In order to investigate the time-periodic behavior of the layer at high inertial numbers, we performed a Fast Fourier Transform (FFT) analysis of the time series data of the center of mass position at steady state. Figure 16 shows the amplitude spectrum of the center of mass for different inclinations. The amplitude spectrum for inclinations $\theta=22^{\circ}$ to $\theta=28^{\circ}$ shows nearly uniform distribution for all frequencies as shown in Fig. 16(a). At higher inclinations, a dominant frequency with a clear peak starts to appear and is shown in Fig. 16(b). The occurrence of a dominant frequency in the amplitude spectrum confirms that the variation of layer height occurs with a characteristic time period at these large inclinations. The peak frequency keeps moving to lower values as the inclination angle increases, indicating that the time period of oscillations increases with the inclination angle. Figure 16(b) also shows that the amplitude of oscillations at the $\theta=38^{\circ}$ is nearly an order of magnitude higher compared to those at $\theta=30^{\circ}$.

FFT analysis of the bulk layer height, center of mass, and kinetic energy data are shown for three different inclinations i.e. $\theta=30^{\circ}$, $\theta=34^{\circ}$, and $\theta=38^{\circ}$ as shown in Fig. 17(a-c). The black solid line corresponds to the center of mass position, the red solid line corresponds to the bulk flowing layer thickness, and the black solid line corresponds to the kinetic energy. We find that the maximum amplitude of the spectrum is observed at the same frequency $f=0.05$ for different properties at the inclination $\theta=38^{\circ}$. For $\theta=34^{\circ}$ and $\theta=30^{\circ}$, the corresponding frequencies are $f=0.09$ and $f=0.13$ respectively. Hence we conclude that the oscillations at high inclinations affect these properties of the flow in a nearly identical fashion.

We also investigate the role of the simulation box size on the dominant frequency of oscillation and its corresponding amplitude. For this, we have considered three different simulation box size of base area $15d\times 15d$, $20d\times 20d$ and $30d\times 30d$ where $d$ is the diameter of the . The initial height of the simulation setup at $t=0$ for all the three cases is same i.e. $40d$. Fig. 18 shows the FFT analysis of the center of mass data at $\theta=38^{\circ}$ for three different box sizes. We find that the maximum amplitude of the spectrum is observed at the same frequency $f=0.05$ at the inclination $\theta=38^{\circ}$ for three different simulation box size confirming that the amplitude and frequency of the oscillation is independent of the simulation box size.