Analytic inclined flow rule for determining granular rheology under strong non-local effects

Abstract

This study theoretically establishes a flow rule for a granular flow down a rough inclined plane, capable of determining granular rheology in the presence of strong non-local effects resulting from grain cooperativity. To describe the non-local rheology, a Landau–Ginzburg model is formulated in terms of the fluidisation parameter represented by the granular inertial number. The exact solutions of the inertial-number field are solved and provide physical insights into the evolution of the internal rheology and the flow arrest process controlled by the flow height. Through asymptotic analysis in the regime dominated by strong non-locality, the exact solutions are further reduced to yield an analytical flow rule for the mean flow velocity. A comparison between the prediction of the flow rule and experimental data from the literature for sand grains determines the underlying rheology law and the relevant rheological parameters. Thus, the proposed flow rule serves as an effective tool for inferring granular rheology from strongly non-local inclined flow data, surpassing the limitations of the classical flow rule deduced from the local rheology framework.

1. Introduction

Flows of granular materials encompass a wide range of geophysical phenomena, such as debris flows, landslides and avalanches, and find applications in various industries, such as metallurgy, pharmaceuticals and agriculture. Despite their ubiquity, granular flows are often challenging to predict due to complex inter-particle frictional and collisional interactions, rendering the flows to exhibit solid- and fluid-like behaviours (Campbell 2006; Andreotti, Forterre & Pouliquen 2013). To understand such a biphasic feature and to simulate the flow dynamics in relevant applications, a unified model able to describe granular rheology in various flow conditions is desired. Much effort has been made to the so-called local $\mu (I)$-rheology (GDR-Midi 2004; da Cruz et al. 2005; Jop, Forterre & Pouliquen 2006), based on dimensional analysis for homogeneous shear flow to relate an effective frictional coefficient, $\mu =\tau /P$, defined as the ratio between the shear stress $\tau$ and the pressure $P$, to an inertial number

(1.1)\begin{equation} I=\frac{\dot{\gamma} D}{\sqrt{P/\rho}}, \end{equation}

defined as the ratio of the time scales for macroscopic deformation and microscopic grain rearrangement, where $\dot {\gamma }$ is the shear rate, $D$ is the mean grain diameter and $\rho$ is the intrinsic density of solid grains. The inertial number serves as a granular state indicator from the quasi-static, dense inertial, to collisional regime as its magnitude increases. A crucial concern arises regarding the determination of the $\mu (I)$ function, which has been extensively studied through experiments and numerical simulations across various shear flow configurations (GDR-Midi 2004; da Cruz et al. 2005; Jop, Forterre & Pouliquen 2005; Koval et al. 2009; Govender 2016; Barker, Zhu & Sun 2022). One typical method is to infer the $\mu (I)$ function from an empirical flow rule for grains flowing down rough inclined surfaces (Pouliquen 1999; Silbert, Landry & Grest 2003; GDR-Midi 2004). The flow rule involves a relation for the depth-averaged mean velocity, $\bar {u}$, to the flow height, $h$, and to the slope angle of inclination, $\theta$. For spherical grains such as glass beads, the flow rule is observed to follow a well-defined relation

(1.2)\begin{equation} \frac{\bar{u}}{\sqrt{gh}}=\beta \frac{h}{h_{stop}(\theta)}, \end{equation}

known as the Pouliquen flow rule (Pouliquen 1999). Here, $g$ is the acceleration of gravity, $\beta$ is an empirical constant, and $h_{stop}(\theta )$ is a critical height as a function of $\theta$ below which a flow stops. The significance of the Pouliquen flow rule lies in its analytical connection to the local $\mu (I)$ rheology. In the configuration of inclined-plane flow, the rheology predicts the Bagnold velocity profile and the resulting mean velocity relation is consistent with (1.2) (GDR-Midi 2004; Gray & Edwards 2014). Hence, the measured relation (1.2) contains embedded information about the $\mu (I)$ function and the related rheological parameters.

Nevertheless, numerous investigations have revealed that granular rheology as well as its mean flow characteristics can be significantly affected by non-local effects arising from cooperative motions among contacting grains (Pouliquen & Forterre 2009; Reddy, Forterre & Pouliquen 2011; Bouzid et al. 2015; Zhang & Kamrin 2017). Such non-local effects render the rheology at a specific location dependent on flow conditions over the entire domain, thereby failing the local rheology law. In the case of spherical grains such as glass beads, non-local effects are locally significant close to flow arrest, causing the mean velocity to continuously decrease to zero as $h$ is reduced to $h_{stop}$, instead of converging to the finite magnitude $\beta$ as predicted by the Pouliquen flow rule (1.2) (Deboeuf et al. 2006). On the other hand, for irregular grains with rough and faceted shapes such as sand, a pronounced deviation from the Pouliquen flow rule is observed across a wide range of flow heights (Forterre & Pouliquen 2003; Félix & Thomas 2004; Aranson et al. 2008; Malloggi, Andreotti & Clément 2015), indicating that the angularity of grain shape induces a significant non-local effect. While the deviation can be described by a modified Pouliquen flow rule, $\bar {u}/\sqrt {gh}=\gamma +\beta (h/h_{stop})$, with $\gamma$ being an offset (Forterre & Pouliquen 2003; Félix & Thomas 2004; GDR-Midi 2004; Malloggi et al. 2015), or by another form of the flow rule derived from a modified kinetic theory of dense granular gases (Jenkins 2006; Börzsönyi & Ecke 2007), their connections to the underlying $\mu (I)$ rheology is unclear.

A promising method incorporating non-local effects into the $\mu (I)$ rheological framework is the Landau–Ginzburg approach, originally developed in the context of thermodynamic phase transitions (Aranson & Tsimring 2002). This approach introduces an order parameter, $\lambda$, as a field variable to quantify the degree of material fluidisation at a specific location. Its evolution in space and time is governed by a reaction–diffusion equation, known as the Landau–Ginzburg equation, which includes source terms associated with the local rheology and diffusive terms that describe non-locality resulting from cooperative motions. Within the framework, Aranson & Tsimring (2002) first proposed a partially fluidised theory, which postulates an order parameter to characterise the ratio of solid to fluid stress and a polynomial source term to describe the local rheology. Subsequently, a granular fluidity theory was developed by adapting a non-local fluidity theory for emulsions to the granular Landau–Ginzburg framework (Kamrin & Koval 2012; Kamrin & Henann 2015). More recently, Lee & Yang (2017) treated the inertial number as an order parameter and used scaling arguments to formulate the Landau–Ginzburg equation. These models have demonstrated their ability to capture many significant features of granular inclined flows, including creep-flow behaviour (Aranson & Tsimring 2002; Kamrin & Henann 2015; Lee & Yang 2017), the Bagnold profile (Kamrin & Henann 2015; Lee & Yang 2017) and hysteresis of the flow thresholds (Aranson & Tsimring 2002; Aranson et al. 2008; Lee & Yang 2017; Mowlavi & Kamrin 2021).

Despite these efforts, an analytical solution for the inclined flow rule derived from granular non-local rheology is still missing. The difficulties lie in the mathematical complexity arising from the introduction of non-local mechanisms and the uncertainty of a unified formulation of the Landau–Ginzburg model. Yet such a solution holds significance, as the irregular shapes exhibited by most grains encountered in reality give rise to pronounced non-local effects that challenge the applicability of the Pouliquen flow rule. Hence, this study focuses on granular materials with strong non-local effects and aims to derive an analytic inclined flow rule from the Landau–Ginzburg framework. As will be demonstrated, this is accomplished by utilising the inertial-number-based Landau–Ginzburg approach (Lee & Yang 2017), which enables the derivation of exact solutions for the order parameter and asymptotic solutions for the mean flow velocity in the regime dominated by non-locality. The derived flow rule is shown to well capture literature data for sand and provide insights into the rheology of sandy materials.

2. Theoretical analysis

2.1. Landau–Ginzburg granular rheology model

By treating granular material as a bi-phasic system undergoing a solid–fluid phase transition, a Landau–Ginzburg equation describing the evolution of the fluidisation parameter, $\lambda$, is formulated as

(2.1)\begin{equation} t_0 \frac{{\rm D} \lambda}{{\rm D} t}=l^2 \nabla^2 \lambda +r(\mu)\lambda -B \lambda^{1+b}. \end{equation}

Here, $t_0$ represents a relaxation time, $l$ denotes a microscopic length associated with non-local transport and $B>0$ and $b>0$ are rheological parameters. The function $r(\mu )$ controls the occurrence of granular phase transitions at a critical effective frictional coefficient, $\mu _c$, requiring that $r(\mu )>0$ for $\mu >\mu _c$ and $r(\mu )<0$ for $\mu <\mu _c$. The functional form of $r(\mu )$ will be determined by comparing the model solution with empirical data. Note that, in the partially fluidised theory (Aranson & Tsimring 2002; Aranson et al. 2008) and the granular fluidity theory (Kamrin & Henann 2015; Mowlavi & Kamrin 2021), the source terms are assumed a priori to present the form of an analytic expansion with integer power exponents. However, Lee & Yang (2017) have demonstrated that the power exponent in the dissipative term, $1+b$, is crucial for accurately reproducing the observed rheology, so here, $b$ is treated as a free parameter to be determined. Note also that an additional bistability term may be introduced in (2.1) to account for hysteresis of flow onset and arrest (Aranson & Tsimring 2002; Lee & Yang 2017; Mowlavi & Kamrin 2021), which is, however, beyond the scope of this study focusing on the flow-to-arrest process.

Three potential candidates for the physical definition of $\lambda$ are discussed as follows. In the partially fluidised theory (Aranson & Tsimring 2002; Aranson et al. 2008), $\lambda$ conceptually denotes the concentration of granular solid phase to fluid phase, requiring an additional empirical relation to relate it to the flow variables, thereby introducing further complexity. In the granular fluidity theory (Kamrin & Koval 2012; Kamrin & Henann 2015; Mowlavi & Kamrin 2021), $\lambda$ is defined as a granular fluidity parameter, the inverse of effective viscosity rescaled by pressure, $\dot {\gamma }/\mu$. This definition introduces dimensional aspects and makes the coefficient $B$ a flow variable (Kamrin & Henann 2015), posing challenges for analytical analysis. For the sake of simplicity, here $\lambda$ is chosen as the inertial number, i.e. $\lambda = I$, which is dimensionless and defined directly by the flow variables (Lee & Yang 2017). With this choice, in a steady homogeneous shear flow, (2.1) yields a local rheology relation given by

(2.2)\begin{equation} I_{loc}(\mu)=\left( \frac{r(\mu)}{B} \right)^{1/b}. \end{equation}

2.2. Application to inclined-plane flow

Consider a steady granular flow on a rough inclined plane, where the flow height $h$ is uniform and the slope angle $\theta$ is measured from the horizon. A Cartesian coordinate is set at the bottom somewhere upstream, where $x$ points down the inclined plane and $z$ is the direction normal to the plane. Assuming the material is incompressible with a constant solid volume fraction $\phi$, the shear stress $\tau$ and the pressure $P$ can be derived from the momentum balance equations at steady state, given by

(2.3)$$\begin{gather} \tau = \rho \phi g \sin\theta (h-z), \end{gather}$$

(2.4)$$\begin{gather}P = \rho \phi g \cos\theta (h-z), \end{gather}$$

which leads to the effective frictional coefficient

(2.5)\begin{equation} \mu = \tan \theta. \end{equation}

In this configuration, the Landau–Ginzburg equation (2.1) reduces to

(2.6)\begin{equation} 0=l^2 \frac{{\rm d}^2 I}{{\rm d} z^2} +r(\theta)I -B I^{1+b}. \end{equation}

Multiplied by ${\rm d}I/{\rm d}z$, the (2.6) can be manipulated to yield

(2.7)\begin{equation} 0=\frac{{\rm d}}{{\rm d}z}\bigg[ \frac{l^2}{2} \left(\frac{{\rm d} I}{{\rm d}z}\right)^2 +\frac{r(\theta)}{2}I^2-\frac{B}{2+b} I^{2+b} \bigg], \end{equation}

or

(2.8)\begin{equation} \frac{l^2}{2}\left(\frac{{\rm d} I}{{\rm d} z}\right)^2 +\frac{r(\theta)}{2}I^2-\frac{B}{2+b} I^{2+b} =\textrm{const}. \end{equation}

To proceed with the analysis, the boundary conditions for the order parameter at the free surface and bottom are required. Since no fluidisation enters or leaves the free surface, the surface condition is typically assumed to be flux free (Aranson & Tsimring 2002; Kamrin & Henann 2015; Lee & Yang 2017)

(2.9)\begin{equation} \frac{{\rm d} I}{{\rm d} z} =0,\quad \textrm{at}\ z=h. \end{equation}

For irregular grains flowing down a rough plane considered here, the bottom can be fairly assumed to be highly dissipative, indicating a pure solid phase

(2.10)\begin{equation} I =0,\quad \textrm{at}\ z=0. \end{equation}

This bottom condition is supported by direct experimental observation of sand flows on a rough inclined plane in which a thin jammed layer forms at the bottom below the flowing layer (Aranson et al. 2008; Malloggi et al. 2015). Now, with the surface condition, (2.8) reduces to

(2.11)\begin{equation} \frac{l^2}{2}\left(\frac{{\rm d} I}{{\rm d} z}\right)^2 +\frac{r(\theta)}{2}I^2-\frac{B}{2+b} I^{2+b} = \frac{r(\theta)}{2}I_s^2-\frac{B}{2+b} I_s^{2+b}, \end{equation}

where $I_s$ denotes the inertial number at the free surface. In terms of the following normalised variables:

(2.12a–c)\begin{equation} I^*=\frac{I}{I_s},\quad z^*=\frac{z}{h},\quad \epsilon = \frac{I_s}{I_{loc}(\theta)}, \end{equation}

where $I_{loc}$ is given by (2.2), (2.11) can be rewritten as

(2.13)\begin{equation} \frac{{\rm d} I^*}{{\rm d}z^*}= \frac{\rm \pi}{2} \kappa \sqrt{ 1-I^{*2} -\frac{2}{2+b}\epsilon ^{b} ( 1 - I^{*2+b})}, \end{equation}

where

(2.14)\begin{equation} \kappa = \frac{2\sqrt{B}I_{loc}(\theta)^{b/2} h}{{\rm \pi} l}. \end{equation}

With the bottom condition $I^*=0$ at $z^*=0$, (2.13) gives rise to an integral equation

(2.15)\begin{equation} \frac{\rm \pi}{2}\kappa z^* = \int^{I^*}_0 \frac{{\rm d}X}{\sqrt{ 1-X^{2} -\dfrac{2}{2+b}\epsilon ^{b} ( 1 - X^{2+b})}}. \end{equation}

Substitution of the surface condition $I^*=1$ at $z^*=1$ into (2.15) gives

(2.16)\begin{equation} \frac{\rm \pi}{2}\kappa = \int^{1}_0 \frac{{\rm d}X}{\sqrt{ 1-X^{2} -\dfrac{2}{2+b}\epsilon ^{b} ( 1 - X^{2+b})}}. \end{equation}

The above two equations suggest $I^*=I^*(z^*,\kappa,b)$ and $\epsilon =\epsilon (\kappa,b)$. With (2.14) and (2.16), the stopping height $h_{stop}$ can be determined by prescribing the flow-arrest condition $\epsilon =0$, given by

(2.17)\begin{equation} \kappa_{\epsilon=0}=\frac{2\sqrt{B}I_{loc}(\theta)^{b/2} h_{stop}}{{\rm \pi} l} =1. \end{equation}

By using (2.2) to replace $I_{loc}$ in (2.17), the stopping height is expressed in terms of $r(\theta )$ as

(2.18)\begin{equation} h_{stop}=\frac{{\rm \pi} l}{2\sqrt{r(\theta)}}. \end{equation}

As a result, the relation (2.14) can be rewritten as

(2.19)\begin{equation} \kappa = \frac{h}{h_{stop}(\theta)}, \end{equation}

which reveals that $\kappa$ is the ratio of the flow height to the stopping height.

2.3. Determination of $b$

The parameter $b$ is determined by comparing the mean velocity relation with the phenomenological observation. The mean velocity is evaluated by integrating the local velocity over the depth divided by the flow height, $\bar {u} = (1/h)\int ^{h}_0 u(z)\, {\rm d} z$. By exploiting integration by parts and assuming a no-slip condition at the bottom, $u=0$ at $z=0$,

(2.20)\begin{equation} \bar{u} = \left.\frac{1}{h} (z-h) u\right|^h_0 -\frac{1}{h} \int^{h}_0 (z-h) \frac{{\rm d}u}{{\rm d} z}\,{\rm d}z=\frac{1}{h} \int^{h}_0 (h-z) \frac{{\rm d}u}{{\rm d} z}\,{\rm d}z. \end{equation}

By replacing ${\rm d}u/{\rm d}z$ with $I$ using (1.1)

(2.21)\begin{equation} \bar{u} =\frac{\sqrt{\phi g \cos\theta }}{Dh} \int^{h}_0 I (h-z)^{3/2} \, {\rm d}z. \end{equation}

In terms of the normalised quantities (2.12a–c), the above equation can be rewritten as

(2.22)\begin{equation} Fr = \sqrt{\phi \cos\theta} \frac{h}{D} I_{loc}(\theta) \epsilon(\kappa,b) \int^{1}_0 I^*(z^*,\kappa,b) (1-z^*)^{3/2} \, {\rm d} z^*, \end{equation}

where $Fr\equiv \bar {u}/\sqrt {gh}$ is the Froude number. Note that the $\sqrt {\cos \theta }$ term can be treated as a constant since it varies slightly within the steady flow range (Jop et al. 2005). With the help of (2.2), (2.18) and (2.19), (2.22) can be manipulated to be

(2.23)\begin{equation} Fr = \sqrt{\phi\cos\theta}\bigg(\frac{2\sqrt{B}D}{{\rm \pi} l}\bigg)^{-{2}/{b}} \left[\frac{D}{h_{stop}(\theta)}\right]^{{2}/{b}-1} \kappa \epsilon(\kappa,b) Q(\kappa,b), \end{equation}

where

(2.24)\begin{equation} Q=\int^{1}_0 I^*(z^*,\kappa,b) (1-z^*)^{3/2} \, {\rm d} z^*. \end{equation}

There has been substantial experimental and numerical evidence showing that the Froude number is one to one with the height ratio $h/h_{stop}$ over a wide range of $h/h_{stop}$ (Pouliquen 1999; Silbert et al. 2003; Deboeuf et al. 2006; Aranson et al. 2008; Malloggi et al. 2015). To be consistent with this observation,

(2.25)\begin{equation} b=2, \end{equation}

has to be selected such that $Fr$ depends on $\kappa$ only

(2.26)\begin{equation} Fr =\alpha_0 \kappa \epsilon (\kappa ) Q(\kappa), \end{equation}

with

(2.27)\begin{equation} \alpha_0 = \frac{{\rm \pi} l}{2D} \sqrt{\frac{\phi\cos\theta}{B}}. \end{equation}

The above analysis significantly demonstrates that the choice of $b=2$, leading to the cubic term in the Landau–Ginzburg equation (2.1), is crucial to reproducing the observed mean flow characteristic. This explains why the common formulation with a quadratic dissipative term fails to yield the one-to-one $Fr-\kappa$ relation (Kamrin & Henann 2015). Note that the previous study by Lee & Yang (2017) employed scaling analysis to derive the same result (2.25) to reproduce the Bagnold scaling in the dense inertial regime. In contrast, the mean velocity analysis adopted here is not limited to a specific flow regime, indicating that the cubic term in the Landau–Ginzburg equation is a general result valid across the quasi-static and the dense inertial regime.

2.4. Exact solutions of the internal flow profiles

With $b=2$, (2.15) and (2.16) reduce to

(2.28)\begin{equation} \frac{\rm \pi}{2}\kappa z^* = \int^{I^*}_0 \frac{{\rm d}X}{\sqrt{ 1-X^{2}} \sqrt{1-\frac{1}{2}\epsilon ^{2} ( 1 + X^{2})}}, \end{equation}

and

(2.29)\begin{equation} \frac{\rm \pi}{2}\kappa = \int^{1}_0 \frac{{\rm d}X}{\sqrt{ 1-X^{2}}\sqrt{1 -\frac{1}{2}\epsilon ^{2} ( 1 + X^{2})}}, \end{equation}

respectively, which can be further manipulated to yield

(2.30)\begin{equation} \frac{\rm \pi}{2}\sqrt{1-\frac{1}{2}\epsilon^2}\kappa = K\left(\sqrt{\frac{\epsilon^2}{2-\epsilon^2}}\right), \end{equation}

and

(2.31)\begin{equation} z^* K\left(\sqrt{\frac{\epsilon^2}{2-\epsilon^2}}\right) = F\left(I^*;\sqrt{\frac{\epsilon^2}{2-\epsilon^2}}\right). \end{equation}

Here, $K(\sqrt {{\epsilon ^2}/({2-\epsilon ^2})})$ and $F(I^*;\sqrt {{\epsilon ^2}/({2-\epsilon ^2})})$ are the complete and incomplete elliptic integrals of the first kind, respectively (Abramowitz & Stegun 1964). Equation (2.30) indicates that $\epsilon$ is a function of $\kappa$ only. On the other hand, (2.31) describes the cross-depth profiles of the normalised inertial number $I^*$, whose shapes depend on $\epsilon$ and hence $\kappa$.

The exact solutions (2.30) and (2.31) reveal that the internal flow rheology evolves with the height ratio $\kappa$. Figure 1 plots the profile of $\epsilon$ vs $\kappa$. It shows that, as $\kappa$ increases from unity, $I_s$ monotonically increases from zero and eventually saturates to $I_{loc}$. Figure 2(a) displays the cross-depth profiles of $I/I_s$ at various $\kappa$. For a very thick flow as $\kappa \gg 1$, the profile is nearly uniform except at the very bottom where a boundary layer with large gradients of $I$ forms to conform to the solid-phase bottom. According to the $\epsilon -\kappa$ profile shown in figure 1, the values of $I$ in the uniform region approximate $I_{loc}$. This indicates that very thick flows fall into the regime dominated by local rheology. On the other hand, for a thin flow with $\kappa \gtrsim 1$, the profiles of $I$ are non-uniform across the layer. This case corresponds to the regime dominated by non-local effects, where the solid-phase bottom strongly correlates with the rheology throughout the depth. The transition between the local and the non-local features of the two regimes can be further illustrated by plotting the cross-depth profiles of $I/I_{loc}$ at a particular slope angle $\theta$ in figure 2(b). This shows that as $\kappa$ is increased, the profiles gradually collapse onto the uniform distribution $I/I_{loc}=1$, while when $\kappa$ is decreased to unity, the profile gradually vanishes. The above results provide physical insights into how the flow height controls the flow-to-no-flow transition: as the flow becomes thinner, the non-local transport capability is progressively enhanced, causing the bottom solid phase to spread and diminish fluidisation over the entire domain.

Figure 1. Profile of $\epsilon$ vs $\kappa$ given by (2.30).

Figure 2. Depth profiles of (a) the inertial number $I$ normalised by the surface magnitude $I_s$, (b) the inertial number $I$ rescaled by the local rheology $I_{loc}$ at $\theta =32.2^\circ$, where the functional form and the parameters in $I_{loc}$ are obtained in § 3, and (c) downslope velocity $u$ normalised by the surface magnitude; (a–c) $\kappa =1.004,1.111,1.987,3.736,5.81,13.07$. Dashed line: the Bagnold profile given by (2.33).

Note that similar transitions in the inertial number profiles have been predicted by the numerical solutions of the $I$-based Landau–Ginzburg equation with a hysteretic bistability term reported in Lee & Yang (2017). This suggests that the hysteresis effect does not significantly influence the non-local flow behaviour. However, the hysteresis effect can result in sudden-jump transitions between flow and arrest (Pouliquen & Forterre 2002; Lee & Yang 2017; Edwards et al. 2019; Mowlavi & Kamrin 2021), which cannot be captured by the non-hysteresis exact solutions (2.30) and (2.31) which predict the transition to be smooth.

The rheological features of the locality- and the non-locality-dominated regime is also reflected in the shapes of the velocity profiles. With the no-slip bottom, the velocity is evaluated by $u=\int ^{z}_0 ({\rm d}u/{\rm d}z) \, {\rm d}\tilde {z}$. By replacing ${\rm d}u/{\rm d}z$ with $I$ and using the normalised quantities (2.12a–c), the velocity is derived to be

(2.32)\begin{equation} u=\frac{\sqrt{\phi g h^3 \cos\theta }}{D} I_{loc}(\theta) \int^{z^*}_0 I^*(\tilde{z},\kappa) (1-\tilde{z})^{1/2} \, {\rm d} \tilde{z}. \end{equation}

Figure 2(c) displays the normalised velocity profiles given by (2.30), (2.31) and (2.32) for the values of $\kappa$ corresponding to those in figure 2(a). For the large values of $\kappa$ in the locality-dominated regime, the velocity profiles collapse onto the Bagnold profile (obtained by setting $I^*=1$ in (2.32)), i.e.

(2.33)\begin{equation} u_{Bag} =\frac{2}{3}\frac{\sqrt{\phi g h^3 \cos\theta }}{D} I_{loc}(\theta) [1-(1-z^*)^{3/2} ], \end{equation}

plotted as the dashed line in figure 2(c). In the non-locality-dominated regime with $\kappa$ close to unity, the profile shapes become noticeably concave near the bottom, characterising creep-flow behaviours. The predicted shape transition is qualitatively consistent with the simulation findings of inclined-plane flows of polyhedral irregular grains (Azéma et al. 2012). Similar theoretical predictions of the velocity shapes have been reported in the previous Landau–Ginzburg modelling studies (Kamrin & Henann 2015; Lee & Yang 2017).

Note that, in the locality-dominated regime with a vanishing boundary layer as $\epsilon \rightarrow 1$ and $I^*\rightarrow 1$, the Pouliquen flow rule is obtained from (2.26), given by

(2.34)\begin{equation} \frac{\bar{u}}{\sqrt{gh}} =\frac{2}{5}\alpha_0 \frac{h}{h_{stop}(\theta)}. \end{equation}

This demonstrates that the flow with a highly dissipative bottom considered here can be captured by the Pouliquen flow rule when it is sufficiently thick. However, this situation takes place when $\kappa \gg 10$, according to the fact that the boundary layer thickness of $I$ is predicted to be approximately one tenth of the flow height at $\kappa \approx 10$ as shown in figure 1(b). This result agrees with the observation that the Pouliquen flow rule does not emerge within the typical experimental ranges for sandy flows, $1<\kappa <20$ (GDR-Midi 2004; Deboeuf et al. 2006; Börzsönyi & Ecke 2007; Aranson et al. 2008; Malloggi et al. 2015), which should fall in the non-locality-dominated regime.

2.5. Asymptotic solution for the flow rule

In the non-locality-dominated regime as $\kappa \rightarrow 1$, the asymptotic solution of the mean flow rule (2.26) is sought as follows. In this limit, since $\epsilon ^2 \rightarrow 0$ as $\kappa \rightarrow 1$ according to (2.30), (2.28) can be asymptotically expanded in terms of $\epsilon ^2$, resulting in

(2.35)\begin{equation} \frac{\rm \pi}{2}\kappa z^* = \int^{I^*}_0 \frac{{\rm d}X}{\sqrt{1-X^{2}}}+\frac{\epsilon^2}{4} \int^{I^*}_0 \frac{1+X^2}{\sqrt{ 1-X^{2}}}\, {\rm d} X +O(\epsilon^4). \end{equation}

After evaluating the integrals, the equation becomes

(2.36)\begin{equation} \frac{\rm \pi}{2}\kappa z^* = \sin^{-1}(I^*)+\frac{\epsilon^2}{8}[3 \sin^{-1}(I^*)-I^*\sqrt{1-I^{*2}}] +O(\epsilon^4). \end{equation}

By inserting the free-surface condition $I^*=1$ at $z^*=1$ into (2.36), an asymptotic relation for $\kappa (\epsilon )$ is then obtained, given by

(2.37)\begin{equation} \kappa -1 = \tfrac{3}{8} \epsilon^2 +O(\epsilon^4), \end{equation}

which indicates an asymptotic relation for $\epsilon$ in terms of $\kappa -1$

(2.38)\begin{equation} \epsilon^2 = \tfrac{8}{3} (\kappa -1) + O[(\kappa-1)^2]. \end{equation}

To evaluate the integral relation (2.24) for $Q$, the explicit solution of $I^*$ is derived as follows. First, (2.36) is rewritten as

(2.39)\begin{equation} I^*=\sin\left( \frac{4 {\rm \pi}}{8+3\epsilon^2} \kappa z^* + \frac{\epsilon^2}{8+3\epsilon^2} I^*\sqrt{1-I^{*2}}+O(\epsilon^4) \right). \end{equation}

By replacing $\epsilon ^2$ with $\kappa -1$ given by the relation (2.38), the above equation reduces to

(2.40)\begin{equation} I^*=\sin\left( \frac{\rm \pi}{2} z^* + \frac{1}{3} \frac{\kappa-1}{\kappa }I^*\sqrt{1-I^{*2}}+ O[(\kappa-1)^2]\right), \end{equation}

which can be expanded in terms of $\kappa -1$ to yield

(2.41)\begin{equation} I^*= \sin\left( \frac{\rm \pi}{2} z^*\right) + \frac{1}{3} (\kappa-1) \cos\left(\frac{\rm \pi}{2}z^*\right) I^*\sqrt{1-I^{*2}}+ O[(\kappa-1)^2]. \end{equation}

According to the leading-order behaviour $I^*\sim \sin ({\rm \pi} z^*/2)$, (2.41) can then be asymptotically reduced to

(2.42)\begin{equation} I^*=\sin\left( \frac{\rm \pi}{2} z^*\right) + \frac{1}{3} (\kappa-1) \sin\left( \frac{\rm \pi}{2} z^*\right) \cos^2\left(\frac{\rm \pi}{2}z^*\right)+ O[(\kappa-1)^2]. \end{equation}

Substitution of (2.42) into (2.24) gives

(2.43)\begin{align} Q &= \int^{1}_0 \sin\left( \frac{\rm \pi}{2} z^*\right) (1-z^*)^{3/2} \, {\rm d} z^* \nonumber\\ & \quad + \frac{1}{3} (\kappa-1) \int^{1}_0 \sin\left( \frac{\rm \pi}{2} z^*\right) \cos^2\left(\frac{\rm \pi}{2}z^*\right) (1-z^*)^{3/2} \, {\rm d} z^*+O[(\kappa-1)^2]. \end{align}

After evaluating the integrals, $Q$ is derived to be

(2.44)\begin{equation} Q \approx 0.1625+0.03227 (\kappa -1) + O[(\kappa-1)^2]. \end{equation}

Substitution of (2.38) and (2.44) into (2.26) yields an asymptotic relation for $Fr$ in terms of $\kappa -1$, given by

(2.45)\begin{equation} Fr = \alpha_0 C_0 (\kappa -1)^{1/2} [1+ C_1 (\kappa -1) +O[(\kappa-1)^2] ], \end{equation}

where $C_0=0.265$ and $C_1=0.407$. Finally, the analytical solution of the flow rule (2.45) to the $(\kappa -1)^{3/2}$-order can be explicitly expressed in terms of $u/\sqrt {gh}$ and $h/h_{stop}$ as

(2.46)\begin{equation} \frac{\bar{u}}{\sqrt{gh}} = \alpha_0 C_0 \left(\frac{h}{h_{stop}(\theta)} -1\right)^{1/2} \left[1+ C_1 \left(\frac{h}{h_{stop}(\theta)} -1\right) \right]. \end{equation}

Note that $\alpha _0$ is to be determined by fitting of the solution (2.46) to empirical data.

3. Comparison with empirical data

The analytical solutions are compared with experimental measurements of inclined-plane flow for sand reported in Malloggi et al. (2015). In the experiment, the material is Fontainebleau sand of a narrow size distribution around $D=312\pm 60\ \mathrm {\mu }{\rm m}$, and the width of the plane was around 1200 times larger than the grain size, so the lateral wall effects can be neglected, justifying the use of the current two-dimensional solutions. Figure 3 shows the experimental measurement of the scaled stopping height $h_{stop}/D$ vs $\tan \theta$ (open circles), which is well fitted by the function (solid line)

(3.1)\begin{equation} \frac{h_{stop}}{D}= A \frac{\mu_2-\tan\theta}{\tan\theta-\mu_c}, \end{equation}

where the fit parameters are $\mu _c=0.559$, $\mu _2=1.423$ and $A=0.73$. By comparing the fit function (3.1) with (2.18), $l/D=(2/{\rm \pi} )A=0.465$ is obtained and the functional form of $r(\mu )$ is determined in replacement of $\tan \theta$ with $\mu$, given by

(3.2)\begin{equation} r(\mu)=\begin{cases} \left(\displaystyle\frac{\mu-\mu_c}{\mu_2-\mu}\right)^2 & \text{for $\mu>\mu_c$}, \\ -\left(\displaystyle\frac{\mu-\mu_c}{\mu_2-\mu}\right)^{2} & \text{for $\mu<\mu_c$}. \end{cases} \end{equation}

According to (2.2), the local $\mu (I)$ rheology for sand is obtained as

(3.3)\begin{equation} \mu=\mu_c+(\mu_2-\mu_c)\frac{I}{I_0+I}, \end{equation}

where $I_0=B^{-1/2}$ is a characteristic inertial number. Note that this form is identical to the typical $\mu (I)$ function used for glass beads determined by the Pouliquen flow rule (Jop et al. 2005). Yet it should be noted that the form of the $\mu (I)$ function varies if a different fit function other than (3.1) is used.

Figure 3. The stopping height $h_{stop}$ vs the tangent of slope angle $\theta$. Symbols: experimental data in Malloggi et al. (2015). Solid line: the fit function (3.1).

Figure 4 shows the experimental measurement (open circles) and the analytical prediction (2.46) (solid line) for $u/\sqrt {gh}$ vs $h/h_{stop}$, where the $\sqrt {\cos \theta }$ term in the relation for $\alpha _0$ (2.27) is calculated with the average angle $\theta = 42.05^\circ$ of the minimum $tan^{-1}(\mu _c)$ and the maximum $tan^{-1}(\mu _2)$. As can be seen, the flow rule (2.45) quantitatively captures the experimental data over the experimental range with the fitted value of $\alpha _0=0.4742$. Using the relation (2.27) and taking $\phi =0.6$ for dense flows as commonly adopted (Jop et al. 2005), along with the fitted values of $\alpha _0$ and $l/D$, the parameter $B=0.248$ and hence $I_0=2.025$ in the local rheology (2.2) is finally determined.

Figure 4. The Froude number $\bar {u}/\sqrt {gh}$ vs the height ratio $h/h_{stop}$. Solid line: analytical solution (2.46). Dashed line: numerical solution solved from (2.26), (2.30) and (2.31). Symbols: experimental data in Malloggi et al. (2015), where the legend shows the corresponding symbols to different slope angles. Inset: the relative error of the predicted Froude number between the analytical and numerical solutions vs $h/h_{stop}$.

To assess the validity of the asymptotic approximations, the analytical solution (2.45) is compared with the numerical solution solved from (2.26), (2.30) and (2.31) using the material parameters given by the forgoing fitting. The deviation between the two solutions is quantified by the relative error of the Froude number, $|(Fr_{num}-Fr_{analy})/Fr_{num}|$, as a function of $h/h_{stop}$. As shown in figure 4, the analytical solution (solid line) and the numerical solution (dashed line) show a good quantitative agreement with slight deviations due presumably to the neglected higher-order effects in the current analysis. Nonetheless, the inset of figure 4 reveals that the relative error is below $7\,\%$ over the investigated range, justifying the asymptotic approximations.

Figure 5 compares the obtained local $\mu (I)$ rheology for sand (solid line) and that for glass beads (dashed line) proposed in Jop et al. (2005) with the parameters $\mu _c=0.382$, $\mu _2=0.644$ and $I_0=0.279$. The result shows that the $\mu$ values for sand are higher than those for glass beads across the entire range of $I$. This difference can be attributed to the angularity of sand grains which enhances bulk frictional resistance. Similar enhancements in bulk friction due to angularity effects have been reported for non-spherical grains in other flow configurations (Mandal & Khakhar 2016; Salerno et al. 2018; Fazelpour, Tang & Daniels 2022). To further validate the obtained rheological law, it is compared with the data of contact dynamic simulation from Azéma et al. (2012) for inclined-plane flows of polyhedral grains, which share similar geometric characteristics with sand grains. The simulation evaluated $\mu$ and $I$ by extracting the internal flow data far from the surface and bottom to avoid non-locality near the boundaries. Figure 5 shows that the magnitude of $\mu$ in the theoretical prediction for sand and the simulation results for polyhedra (symbols) are comparable, with slight deviations likely stemming from differing grain properties and angularity.

Figure 5. The local $\mu (I)$ rheology for sand (sold line) and glass beads (dashed line). Symbols: simulation data for polyhedral grains in Azéma et al. (2012).

Furthermore, figure 5 shows that $\mu$ saturates to its maximum value $\mu _2$ at larger $I$ for sand than for glass beads, arising from the larger characteristic $I_0$ for the former. The saturation has been known to indicate the transition from the dense inertial to the collisional regime (Forterre & Pouliquen 2008), suggesting that generating a collisional flow of sand requires imposing a higher deformation rate than glass beads as other flow conditions are fixed. This finding aligns with the simulation results of Azéma et al. (2012), showing flows of polyhedral grains enter the collisional regime at larger $I$ than flows of spheres. It may result from an angularity effect that limits the degree of freedom for irregular grains to undergo collisional interactions.

4. Conclusion

In this study, an inclined flow rule has been analytically established to determine the rheology of granular materials with strong non-local effects. The Landau–Ginzburg model based on the inertial-number approach has been formulated to describe the non-local rheology. By analysing the mean velocity property, an important finding is that the dissipative term in the Landau–Ginzburg equation must take the form of a cubic term to account for the one-to-one relationship between $Fr$ and $h/h_{stop}$ observed in extensive experiments and simulations. It is emphasised that this result is valid for all flow regimes, generalising the scaling analysis in Lee & Yang (2017) in which the same result is obtained but limited to the dense inertial regime. The resulting model allows for obtaining of the exact solutions of the inertial number, expressed in the forms of the complete and incomplete elliptic integrals of the first kind. The solutions have been shown to account for how the internal flow evolves from the locality- to the non-locality-dominated regime and leads to the flow arrest when the flow height is reduced. Through asymptotic expansions in the non-locality-dominated regime, these solutions are reduced to yield the analytical flow rule model (2.46). By comparing the flow rule with the literature experimental data on sand, the control function $r(\mu )$ and the model parameters $l$, $B$ and $\alpha _0$ are determined. These results establish the local $\mu (I)$ rheological law for sand, showing consistently higher $\mu$ than the magnitudes of glass beads due presumably to the angularity effect. It is worth noting that the obtained $\mu (I)$ relation can also be employed to develop other types of non-local models based on the inertial-number rheology (Bouzid et al. 2013; Kamrin & Henann 2015). Hence, the proposed flow rule (2.46) provides an effective approach for determining granular rheology with significant non-local effects, which cannot be achieved using the classical Pouliquen flow rule limited to the locality-dominated regime.

The current analysis can be extended to the case of glass beads by considering a more sophisticated bottom flow condition beyond the solid-phase assumption. As mentioned in the introduction, glass bead flows revert to the Pouliquen flow rule once the flow height slightly exceeds the stopping height (Pouliquen 1999; Silbert et al. 2003; Deboeuf et al. 2006). This indicates that the non-locality-dominated regime for glass beads is rather narrow, implying that the bottom flow drastically transitions from a solid jammed state to a Bagnold shearing state. To capture this phenomenon and the resulting inclined flow dynamics, further investigation is warranted to develop an effective bottom boundary condition describing the jammed-to-shearing transition.