Revisiting the role of friction coefficients in granular collapses: confrontation of 3-D non-smooth simulations with experiments

Abstract

In this paper, transient granular flows are examined both numerically and experimentally. Simulations are performed using the continuous three-dimensional (3-D) granular model introduced in Daviet & Bertails-Descoubes (ACM Trans. Graph., vol. 35, no. 4, 2016b, p. 102), which represents the granular medium as an inelastic and dilatable continuum subject to the Drucker–Prager yield criterion in the dense regime. One notable feature of this numerical model is to resolve such a non-smooth rheology without any regularisation. We show that this non-smooth model, which relies on a constant friction coefficient, is able to reproduce with high fidelity various experimental granular collapses over inclined erodible beds, provided the friction coefficient is set to the avalanche angle – and not to the stop angle, as generally done. In order to better characterise the range of validity of the fully plastic rheology in the context of transient frictional flows, we further revisit scaling laws relating the shape of the final collapse deposit to the initial column aspect ratio, and accurately recover established power-law dependences up to aspect ratios of the order of 10. The influence of sidewall friction is then examined through experimental and simulated collapses with varying channel widths. The analysis offers a comprehensive framework for estimating the effective flow thickness in relation to the channel width, thereby challenging previously held assumptions regarding its estimation in the literature. Finally, we discuss the possibility to extend the constant coefficient model with a hysteretic model in order to refine the predictions of the early-stage dynamics of the collapse. This illustrates the potential effects of such phenomenology on transient flows, paving the way to more elaborate analysis.

1. Introduction and related work

Granular materials exhibit state transitions: under specific pressure and shear rates, the grain assembly may either behave as a fluid, a solid or a gas (Andreotti, Forterre & Pouliquen 2013). As an example, avalanches first flow over complex topographies as viscoplastic fluids and ultimately stop their course as a static solid deposit. Similarly, during granular impact cratering, the falling object ejects grains outwards, forming a gaseous entity for a brief instant. Performing efficient and reliable predictions of such flows remains challenging and raises several issues, from the fundamental underlying theory to applied industrial or geophysical processes.

In this context, two families of numerical approaches have been employed: the discrete element method (DEM), which consists in modelling the dynamics of each grain and its interactions; and the continuum-based method, which describes the granular material at a macroscopic scale, relying on rheological closures to capture the phase transitions and relate the state of stress to the material deformation. By modelling intergrain forces such as friction, adhesion or elasticity, DEM is the most ab initio approach to simulate the behaviour of granular materials, and has been used extensively to simulate various model experiments in unsteady (Staron & Hinch 2005; Lacaze, Phillips & Kerswell 2008) or steady (Silbert et al. 2001; Da Cruz et al. 2005; Azéma & Radjai 2014) regimes, yielding valuable insights about the relation between the microstructure and the macroscopic properties. However, DEM inherently suffers from high computational costs, which in practice restricts its use to small-scale systems. In contrast, continuum methods provide opportunities to simulate larger scale scenarios using depth-averaged (Naaim, Faug & Naaim-Bouvet 2003; Balmforth & Kerswell 2005; Moretti et al. 2012) or full two-dimensional (2-D) or three-dimensional (3-D) (Lagrée, Staron & Popinet 2011; Dunatunga & Kamrin 2015; Gaume et al. 2018) fluid models, factoring most of the microstructure complexity in macroscopic constitutive laws.

Many studies in the last two decades have therefore focused on the formulation of such rheological closures, using model experiments such as the plane shear or inclined plane set-ups. The main feature of a granular material is the presence of a pressure-dependent yield stress, characterised by an internal friction coefficient $\mu$, which sets the threshold on the maximal sustainable ratio of the shear stress over the pressure before flow, $\tau / p \leq \mu$. Experimental methods found in the literature to estimate the internal friction coefficient of a material are diverse, ranging from the use of the avalanche angle or the angle of repose of a wedge-type pile (Hutter & Koch 1991; Balmforth & Kerswell 2005), to triaxial tests (Ancey 2001; Adjemian & Evesque 2004).

Further analyses of dense granular flows (see GDR MiDi (2004) for a review) has led to the emergence of the $\mu (I)$ rheology introduced by Jop, Forterre & Pouliquen (2006) to describe granular matter in the liquid – flowing – regime. This viscoplastic rheology features a simple Drucker–Prager yield criterion, albeit with a non-constant friction coefficient depending on the non-dimensional inertial number $I = \dot {\varepsilon } d \sqrt {\rho /p}$, which naturally arises from dimensional analysis when considering strictly hard grains of size $d$ and density $\rho$ subjected to some pressure $p$, and strain rate $\dot {\varepsilon }$. As such, the $\mu (I)$ rheology is able to capture two crucial features of granular flows: (i) the existence of a minimal shear stress for the material to keep flowing, described by the stop friction coefficient $\mu _{stop}= \tan \theta _{stop}$ (with $\theta _{stop}$ the corresponding stop angle, also called repose angle) and (ii) the presence of shear-rate-dependent viscous dissipation at large inertial number $I > I_0 \sim 0.3$, precisely responsible for the existence of steady flows in a range of applied shear-rates.

While extensive testing of the $\mu (I)$ rheology has been performed in steady-flow contexts, the monotonic increase of the friction coefficient with $I$ is unable to describe the transition from rest to flow characterised by the existence of a static friction coefficient larger than the dynamic friction coefficient (Hutter & Koch 1991; Da Cruz et al. 2002; Pouliquen & Forterre 2002; Perrin et al. 2019). This friction gap results in a hysteresis effect in avalanches: when inclining a container initially filled with static grains, the material starts to flow at a higher start angle or avalanche angle $\theta _a$ (corresponding to a friction $\mu _a= \tan \theta _a$) than the stop angle.

As suggested in Pouliquen & Forterre (2002), the $\mu (I)$ law could be extended as a decreasing function of $I$ for small $I \in [0, I_*]$ in order to account for a higher start angle while recovering the classical $\mu (I)$ rheology at larger $I$. Exploring the $I<10^{-2}$ behaviour, DeGiuli & Wyart (2017) pointed out that $\mu (I)$ is indeed non-monotonic, exhibiting a decrease from $\mu (I=0)$ to $\mu ( I_* \sim 10^{-3}) = \mu _{stop}$ which amplitude $\varDelta _{\mu,hyst}$ is induced by endogenous acoustic noise depending on grain rigidity and appliedpressures.

Among the different model experiments available to study granular flows at low inertial number $I$, the simple collapse of a granular column onto a flat or inclined surface is one of the most studied yet distinctive cases (Lajeunesse, Mangeney-Castelnau & Vilotte 2004; Balmforth & Kerswell 2005; Lube et al. 2005; Lacaze & Kerswell 2009; Lagrée et al. 2011; Farin, Mangeney & Roche 2014; Ionescu et al. 2015). Closely reminiscent of natural avalanches or landslides, it challenges both the accuracy of constitutive laws in friction-dominated, unsteady and hysteretic regimes, and the sensitivity of granular flows to boundary conditions (through the presence of frictional sidewalls or a lifting gate for instance).

Granular collapse features were first studied experimentally (Boutreux & de Gennes 1997; Daerr & Douady 1999a; Lajeunesse et al. 2004; Balmforth & Kerswell 2005; Lube et al. 2005; Mangeney-Castelnau et al. 2005). In particular, Lajeunesse, Monnier & Homsy (2005), Lube et al. (2005) and Balmforth & Kerswell (2005) focused on the final deposit shape of collapses by studying the impact of the initial aspect ratio $a=H_0/L_0$ (i.e. the ratio between the initial height $H_0$ and length $L_0$ of the granular column) on the final height $H_f$ and run-out $L_f$, unravelling simple power-law relations between the normalised final heights or run-outs, and the aspect ratio $a$, however, with different regimes for the exponents, depending on the initial aspect ratio (with a regime transition around $a \approx 3$) or the channel width.

Following these first experimental studies, new numerical models demonstrated a good agreement between DEM simulations and the local $\mu (I)$-rheology for axisymmetric collapses on horizontal planes (Lacaze & Kerswell 2009), triggering a series of numerical studies based on continuous models to explore the relevance of the $\mu (I)$ viscoplastic rheology in the context of granular collapses (Lagrée et al. 2011; Dunatunga & Kamrin 2015; Ionescu et al. 2015; Martin et al. 2017; Chupin et al. 2021), and more generally transient flows. In particular, Lagrée et al. (2011) highlighted the ability of the $\mu (I)$-rheology to recover power-law dependences of the final shape on the initial aspect ratio of the column, albeit with a regime transition around $a \sim 7$, slightly higher than Lajeunesse et al. (2005) or Lube et al. (2005). Similar results were obtained by Dunatunga & Kamrin (2015) with a $\mu (I)$-rheology implemented in a material point method (MPM) algorithm.

Although these numerical studies have shown an ability to reproduce collapses with the $\mu (I)$-rheology, it is still, however, not established whether the well-known aforementioned power-law scaling relations originate from the viscoplastic $\mu (I)$-rheology or could be obtained from a purely plastic rheology. The performance of the $\mu (I)$-rheology for transient collapses was even questioned by Lagrée et al. (2011), whose observations indicate that the final deposits of granular collapses could also be recovered with a constant friction coefficient, but at the expense of adjusting the value of the coefficient for large initial aspect ratios. More generally, while purely plastic rheologies can undeniably not account for some well-known characteristic effects of granular flows (such as for instance the Bagnold velocity profiles observed in inclined plane geometries (Bagnold 1954; Silbert et al. 2001)), disentangling the respective roles of plastic friction and viscosity appears essential to devise future discriminatory transient flow experiments.

Beyond rheological matters, the sensitivity of granular collapses to frictional sidewalls has also raised concerns, largely disrupting the conclusions about the physical origins of observations (Balmforth & Kerswell 2005; Lajeunesse et al. 2005; Ionescu et al. 2015; Martin et al. 2017). Containing walls have been shown to play an important role in the analysis of steady inclined plane experiments, owing to the additional lateral frictional condition they involve, which can strongly affect the free-surface flow velocity (Jop, Forterre & Pouliquen 2005). In particular, sidewall friction is responsible for the existence of superstable static inclined piles (Taberlet et al. 2003) in narrow channels, and can substantially affect the apparent rheological parameters for channels of width $W$ below ${\sim }200$ grain diameters (Jop et al. 2005). In this context, Ionescu et al. (2015) have used the linear scaling law relating the incline flow thickness to the channel width (Savage 1979; Taberlet et al. 2003), originally designed in steady configurations, to support an increase of the effective 2-D friction coefficient and better match the 3-D collapse experiments, albeit without proper definition of a flow thickness in such transient set-ups. To our knowledge, this 2-D rescaling has, however, never been verified on granular collapses, and an advanced numerical study including direct account of sidewalls effects in 3-D simulations could provide valuable insight into the exact role of the walls in this unsteady configuration.

1.0.1. Contributions

In this paper, we perform 3-D numerical simulations at the continuum level, as well as experiments of granular column collapses over erodible beds with slopes varying from $0^\circ$ to $20^\circ$. We model the granular material as a dilatable, fully plastic continuum with a Drucker–Prager yield surface and non-associated flow rule. We solve the resulting non-smooth rheology without any regularisation by leveraging the Sand6 software (https://gitlab.inria.fr/elan-public-code/sand6), a numerical scheme introduced by Daviet & Bertails-Descoubes (2016a,b) and inspired from non-smooth optimisation algorithms originally developed in the context of rigid body contact dynamics (Moreau 1994). In such a framework, we can model granular collapses either as viscoplastic or as purely plastic, allowing us to investigate the role of viscosity as well as the relevance of a unique constant friction coefficient.

To this end, we perform numerical simulations of the inclined collapses using the $\mu (I)$-rheology or a constant friction coefficient, set either from the experimental avalanche or stop angles, and compare them with the experiments, faithfully accounting for the boundary conditions and the lifting door.

Our results confirm the weak impact of a $\mu (I)$-rheology in this particular context, as suggested by Ionescu et al. (2015), but also show that experimental granular collapses with various slopes, flume widths and aspect ratios, are better predicted by using the constant friction coefficient corresponding to the (higher) avalanche angle, as opposed to the usual (lower) stop angle provided by the Pouliquen (1999) or Pouliquen & Forterre (2002) experimental protocols.

In order to explore the domain in which this result applies, we perform simulations with a constant friction coefficient for various aspect ratios and show that we recover the various regimes in the power-law dependences of the final collapse heights and run-outs proposed by Balmforth & Kerswell (2005), Lajeunesse et al. (2005), Lube et al. (2005), Lagrée et al. (2011) and Dunatunga & Kamrin (2015), suggesting that a strain-rate dependant model is not necessary to recover the final run-out for a large range of aspect ratios.

The need for a higher avalanche friction coefficient $\mu _a$ instead of the stop friction coefficient $\mu _{stop}$ used in previous studies is then investigated in light of the results of Ionescu et al. (2015). Analysing the effect of wall friction on collapses in narrow flumes with variable width, we show that a linear rescaling of the 2-D friction coefficient as observed by Savage (1979), Taberlet et al. (2003) and Jop et al. (2005) and used by Ionescu et al. (2015) to simulate collapses in two-dimensions remains valid. However, this finding challenges previous assumptions in the literature that overestimated the effective flow thickness using the maximum flow thickness, thus attributing excessive importance to the influence of walls in the considered flume (Ionescu et al. 2015; Martin et al. 2017).

Finally, seeking a way to improve the correspondence between experiment and numerical scenarios at low inertial numbers $I$, we show that a simplified hysteresis rule built in light of previous knowledge (Pouliquen & Forterre 2002; DeGiuli & Wyart 2017) improves predictions in the incipient motion of the granular collapse, paving the way to more elaborate analysis.

Overall, our conclusions significantly deviate from those of previous works, which model granular collapses by using a $\mu (I)$-rheology based on the stop friction coefficient (Lagrée et al. 2011; Dunatunga & Kamrin 2015; Ionescu et al. 2015). Our results may suggest that unsteady collapses are more likely governed by local transitions from rest to flow, and highlight the irrelevance of viscous effects in the context of granular collapses, while steering further studies towards a better quantification of hysteretic laws in such configurations. Along the paper, we support our observations by validating carefully our simulator against experiments and/or prior research in a number of well-controlled scenarios.

The paper is organised as follows. We present in § 2 our modelling strategy combining an exact Drucker–Prager yield criterion with usual fluid conservation equations, along with our constraint-based approach to solve for the corresponding non-regularised discrete equations (cf. § 2.3). We show in § 2.4 that our numerical model is self-consistent using a numerical avalanche test. We then describe our different experimental set-ups for studying granular collapses and measuring macroscopic parameters in § 3, and validate our numerical model with a constant friction coefficient on model experiments in § 4, before showing its predictive potential on collapses of natural sediments. We further analyse experimentally and numerically the impact of the lateral walls – in combination with the initial column aspect ratio – during granular collapses in § 5. Finally, we discuss the possibility to extend the constant coefficient model with a hysteretic model in order to refine the predictions of the early-stage dynamics of the collapse in § 6.

The maintained C++ Sand6 software implementing the semi-implicit scheme to solve the continuum conservation equations with the Drucker–Prager rheology, used to study the role of the unique friction coefficient is available at https://gitlab.inria.fr/elan-public-code/sand6. Simulations of the experimental collapses were performed using a fork of the Sand6 code available at https://gitlab.com/groussea/sand6py which provides a Python binding of the C++ library along with analysis scripts.

2. A non-smooth numerical model for Drucker–Prager flows

In this section, we present our continuum dilatable Drucker–Prager plastic fluid model for granular flows, along with the non-smooth numerical discretisation used to robustly and efficiently simulate true-scale collapses. The numerical model is based on the original reformulation of the Drucker–Prager rheology as a conic constraint introduced in Daviet & Bertails-Descoubes (2016a,b), and implemented in the open-source Sand6 software.

In order to assess the physical relevance of this model in the context of transient granular flows, we numerically reproduce a classical avalanche experiment in two and three dimensions, and analyse the role of the Drucker–Prager friction coefficient in regards to the usual experimental avalanche and stop friction coefficients.

2.1. Continuum modelling

Focusing on the large deformation regimes, we model the granular material as a continuous yield-stress fluid and introduce the corresponding strain-rate tensor

(2.1)\begin{equation} \dot{\boldsymbol{\varepsilon}} = \tfrac{1}{2} ( \boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^{T} ),\end{equation}

with $\boldsymbol {u} \equiv \boldsymbol {u}(\boldsymbol {x},t)$ the velocity of the material in Eulerian coordinates. We also introduce the volume fraction field $\phi (\boldsymbol {x},t)$ and assume that all the grains composing the material have a constant density $\rho _g$, so that the density field is

(2.2)\begin{equation} \rho(\boldsymbol{x},t) = \rho_g \phi(\boldsymbol{x},t). \end{equation}

The mass conservation equation can then be written in terms of the volume fraction field as

(2.3)\begin{equation} \frac{\partial\phi}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} ( \phi \boldsymbol{u} ) = 0.\end{equation}

The momentum conservation equation is

(2.4)\begin{equation} \rho \left( \frac{\partial \boldsymbol{u}}{\partial t} + \left( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{u} \right) - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma} = \boldsymbol{f},\end{equation}

with $\boldsymbol {\sigma } \equiv \boldsymbol {\sigma }(\boldsymbol {x},t)$ is the Cauchy stress tensor and $\boldsymbol {f}$ denotes the external volumic forces.

Note that the volume fraction field $\phi$ takes by definition values between $0$ and some critical value $\phi _c$, which accounts for the maximal ‘packing’ fraction of the grains. As such, the granular material is not supposed uniformly dense, and the model can describe the different phases of granular matter. The varying volume fraction thus tightly couples the mass and momentum conservation equations (2.3) and (2.4) through the phase-dependent constitutive relation $\boldsymbol {\sigma }[\dot {\boldsymbol {\varepsilon }},\phi ]$ required to close the system.

2.2. The plastic Drucker–Prager rheology

The central feature of granular materials is the existence of a pressure-dependent yield stress, somehow reminiscent of the Amontons–Coulomb law for solid friction. As mentioned above, we consider the material as perfectly plastic, following the Drucker–Prager rheology proposed in Drucker & Prager (1952). Note that our choice for such perfectly plastic modelling naturally allows us to overcome the numerical stiffness inherent in the very small time scales resulting from granular elasticity. This is in contrast with the elastoplastic models of Dunatunga & Kamrin (2015), Mast et al. (2015) or Klár et al. (2016), which are theoretically able to account for more complex physics but require in practice artificially lower elastic stiffnesses in order to ensure numerical workability, thereby weakening the physical meaning of elastic modelling for hard granular material.

Introducing the usual isotropic-deviatoric decomposition of tensors

(2.5)\begin{equation} \forall\,{\boldsymbol{\tau}}, \quad \boldsymbol{\tau} = \frac{\mathrm{Tr}(\boldsymbol{\tau})}{\mathtt{d}} \,{\mathbb{I}_\mathtt{d}} + \mathrm{Dev}(\boldsymbol{\tau}) ,\end{equation}

where $\mathtt {d}$ is the dimension of the space; we assume that the dense state ($\phi = \phi _c$) is governed by the constitutive relation

(2.6a)$$\begin{gather} \lVert \mathrm{Dev}(\boldsymbol{\sigma}) \rVert \leq \mu p \quad \text{if Dev }(\dot{\boldsymbol{\varepsilon}}) = 0\ \textrm{(static regime)}, \end{gather}$$

(2.6b)$$\begin{gather}\mathrm{Dev}(\boldsymbol{\sigma}) = \mu p \frac{\mathrm{Dev}(\dot{\boldsymbol{\varepsilon}})}{\lVert \mathrm{Dev}(\dot{\boldsymbol{\varepsilon}})\rVert} \quad\text{if Dev }(\dot{\boldsymbol{\varepsilon}}) \neq 0\ \textrm{(flowing regime)}, \end{gather}$$

where we have defined the pressure $p = - {\mathrm {Tr}(\boldsymbol {\sigma })}/{\mathtt {d}}$ and the Frobenius tensor norm $\lVert \boldsymbol {\tau } \rVert \equiv \sqrt {{\boldsymbol {\tau } :\boldsymbol {\tau }}/{2}}$ associated with the natural inner product $\langle \boldsymbol {\tau }, \boldsymbol {\upsilon } \rangle = {\boldsymbol {\tau } : \boldsymbol {\upsilon }}/{2} = {\mathrm {Tr}(\boldsymbol {\tau }^T \boldsymbol {\upsilon })}/{2}$. The Drucker–Prager yield surface thus defines a second-order cone in the space of principal stresses (cf. (2.6a)), which leads to a simpler numerical treatment as compared with the hexagon-like Mohr–Coulomb model. To the best of our knowledge, the phenomenological differences between the Drucker–Prager or Mohr–Coulomb yield models are generally application dependent, and both surfaces provide very similar results for granular collapses (Rauter, Barker & Fellin 2020).

We also assume that for $\phi < \phi _c$, the material is in a disconnected stress-free state $\boldsymbol {\sigma } = 0$, as proposed in Narain, Golas & Lin (2010) and Dunatunga & Kamrin (2015). In our case, this can be naturally imposed through the Drucker–Prager rheology by simply constraining the pressure to vanish for $\phi < \phi _c$, which can be written concisely as a complementarity condition

(2.7)\begin{equation} 0 \leq \phi_{c} - \phi \perp p \geq 0,\end{equation}

or equivalently $\phi \leq \phi _c$, $p \geq 0$ and $(\phi _c - \phi )\,p = 0$. When the volume fraction is below the packing fraction $\phi _c$, the complementarity relation constrains the pressure to vanish, which imposes $\mathrm {Dev}(\boldsymbol {\sigma }) = 0$ and as a result $\boldsymbol {\sigma } = \mathrm {Dev}(\boldsymbol {\sigma }) - p {\mathbb {I}_\mathtt {d}} = 0$.

The resulting constitutive relation is a constrained multivalued and not everywhere differentiable functional. Solving for the whole system of (2.3) and (2.4) with (2.6a), (2.6b) and (2.7) is therefore numerically challenging, and is often handled through a regularisation of the rheology to reformulate the problem as a complex fluid as in Lagrée et al. (2011), Chauchat & Médale (2014) or Franci & Cremonesi (2019), or solved using a full elastoplastic model with explicit or implicit return-mapping projections as in Dunatunga & Kamrin (2015), Mast et al. (2015) or Klár et al. (2016).

The former approach suffers from viscous artefacts such as effective creeping flows, which can affect some physical observations. The latter approach is able to recover accurately the yield condition at the expense of small time-steps imposed by the small elastic time scale of hard granular material, thereby strongly increasing the computational costs, and intrinsically opposes recompaction after plastic expansion, leading to eventual volume gains.

One can also mention the approach of Ionescu et al. (2015) which is based on an augmented Lagrangian formulation to solve for the corresponding variational inequality, and could therefore in theory deal with the non-differentiability of the rheology but requires in practice a viscous regularisation to ensure viable convergence of the iterative fixed point algorithm.

In contrast, our numerical approach is entirely based on non-smooth optimisation tools, and exploits in particular the similarities between the Drucker–Prager yield criterion and the Amontons–Coulomb friction law to leverage recent developments in the field of contact dynamics.

2.3. Numerical method

As mentioned above, our simulation framework is based on the 3-D non-smooth numerical model introduced in Daviet & Bertails-Descoubes (2016b), which rewrites the plastic Drucker–Prager equations as a non-smooth root-finding problem, and leverages efficient Gauss–Seidel algorithms originally developed in the context of contact dynamics (Moreau 1994; Jean 1999) to solve for the resulting nonlinear equation. Details about this approach can be found in Daviet & Bertails-Descoubes (2016a,b), and are briefly described here for the sake of completeness.

Introducing the material derivative

(2.8)\begin{equation} \frac{{\rm D}\bullet}{{\rm D}t} \equiv \frac{\partial \bullet}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \bullet , \end{equation}

we can rewrite the mass conservation equation (2.3) as

(2.9)\begin{equation} \frac{{\rm D}\phi}{{\rm D}t} + \phi \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0 , \end{equation}

which can be discretised in time with a first-order Euler scheme,

(2.10)\begin{equation} \phi^{(n+1)} = \phi^{(n)} - \delta t \phi^{(n)} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} , \end{equation}

where $\phi ^{(n+1)} \equiv \phi (\boldsymbol {x}(t+\delta t), t+\delta t)$ and $\phi ^{(n)} \equiv \phi (\boldsymbol {x}(t), t) = \phi (t)$. This allows us to linearise the complementarity relation (2.7) by enforcing the $\phi \leq \phi _c$ constraint on the Lagrangian transported volume fraction, namely

(2.11)\begin{equation} 0 \leq \phi(t) \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} + \frac{\phi_{c} - \phi(t)}{\delta t} \perp p \geq 0 , \end{equation}

which is now linear in $\boldsymbol {u}$.

Defining the auxiliary tensor fields

(2.12)\begin{equation} \left. {\begin{array}{@{}c@{}} \boldsymbol{\lambda} \equiv{-} \boldsymbol{\sigma},\\ \boldsymbol{\gamma} \equiv \phi(t) \dot{\boldsymbol{\varepsilon}} + \dfrac{\phi_{c} - \phi(t)}{\mathtt{d} \delta t} {\mathbb{I}_\mathtt{d}} , \end{array}} \right\} \end{equation}

the whole rheology can be written as

(2.13)\begin{equation} \mathcal{DP}(\mu) \equiv \left\{ \begin{array}{@{}ll@{}} \lVert \mathrm{Dev}(\boldsymbol{\lambda}) \rVert \leq \mu \dfrac{\mathrm{Tr}(\boldsymbol{\lambda})}{\mathtt{d}} & \textrm{if} \ \mathrm{Dev}(\boldsymbol{\gamma}) = 0 \textrm{ (static regime)}\\ \begin{aligned} & \mathrm{Dev}(\boldsymbol{\lambda}) ={-}\mu \dfrac{\mathrm{Tr}(\boldsymbol{\lambda})}{\mathtt{d}} \dfrac{\mathrm{Dev}(\boldsymbol{\gamma})}{\lVert\mathrm{Dev}(\boldsymbol{\gamma})\rVert} \\ & \quad 0 \leq \mathrm{Tr}(\boldsymbol{\gamma}) \perp \mathrm{Tr}(\boldsymbol{\lambda}) \geq 0\end{aligned} & \textrm{if} \ \mathrm{Dev}(\boldsymbol{\gamma}) \neq 0\ \textrm{(flowing regime)} \end{array}\right. . \end{equation}

As shown in Daviet & Bertails-Descoubes (2016a,b), this rheological law can be recast as a normal cone inclusion in $\mathbb {R}^{s_\mathtt {d}}$ with $s_\mathtt {d} \equiv {\mathtt {d}(\mathtt {d}+1)}/{2}$ the dimension of $S(\mathtt {d})$ the space of symmetric $\mathtt {d} \times \mathtt {d}$ matrices. It is thus equivalent to a root-finding problem on a generalised Fischer–Burmeister non-smooth function introduced in Fukushima, Luo & Tseng (2002) and modified in Daviet, Bertails-Descoubes & Boissieux (2011) to handle non-symmetric and non-associated second-order cone complementarity problems. We denote here $f_{MFB}$ as this modified Fischer–Burmeister function. Additional details regarding the definition of $f_{MFB}$ are given in § A.1. We then have

(2.14)\begin{equation} (\boldsymbol{\gamma},\boldsymbol{\lambda}) \in \mathcal{DP}(\mu) \iff f_{MFB}(\boldsymbol{\gamma},\boldsymbol{\lambda}) = 0. \end{equation}

The conservation equations are then discretised in space using a hybrid MPM–finite-element method described in § A.2, finally giving the following algebraic problem:

(2.15)\begin{equation} \begin{gathered} \textrm{find } \boldsymbol{\gamma}, \boldsymbol{\lambda} \in S(\mathtt{d}) \times S(\mathtt{d}) \textrm{ subject~to}\\ \left\{\begin{array}{@{}l@{}} \boldsymbol{\gamma} = \boldsymbol{W} \boldsymbol{\lambda} + \boldsymbol{b}\\ (\boldsymbol{\gamma},\boldsymbol{\lambda}) \in \mathcal{DP}(\mu) \end{array},\right. \end{gathered} \end{equation}

which is solved by minimising the functional

(2.16)\begin{equation} \mathcal{F}[\boldsymbol{\lambda}] \equiv \tfrac{1}{2} \lVert\, f_{MFB} (\boldsymbol{W} \boldsymbol{\lambda} + \boldsymbol{b},\boldsymbol{\lambda}) \rVert^2, \end{equation}

using a Gauss–Seidel iterative procedure with a generalised Newton algorithm to solve the local problems. The overall algorithm is summarised in § A.3.

Note that the constraint $(\boldsymbol {\gamma },\boldsymbol {\lambda }) \in \mathcal {DP}(\mu )$ in (2.15) actually denotes a vector concatenation of local constraints: for each finite-element interpolation node $\boldsymbol {x}_i$, we require that the local auxiliary stress and strain rate tensors satisfy the Drucker–Prager yield constraint, i.e. $\forall \, i, (\boldsymbol {\gamma }(\boldsymbol {x}_i),\boldsymbol {\lambda }(\boldsymbol {x}_i)) \in \mathcal {DP}(\mu )$. The yield constraint is thus only imposed strongly at interpolation points, but not at the continuous nor integral level. In practice, due to the strong nonlinearity of the constraint and the MPM interpolation, this can induce small but noticeable constraint violations, which manifest themselves in particular as potential volume losses over the course of long simulations (note that the mass is, however, accurately conserved as imposed by (2.3)). In the following simulations, we have carefully adjusted the numerical parameters to ensure that the total relative volume losses always stay below $3\,\%$.

2.4. Drucker–Prager friction coefficient and yield transition

As shown above, our model solves for a perfectly plastic Drucker–Prager rheology, without regularising the yield criterion. As such, it features a single friction coefficient $\mu$ which defines the yield surface, and thereby the transition between the static and flowing regimes.

To assess the numerical behaviour of our model close to the yield point, we consider a simple avalanche experiment: a 2-D or 3-D dense bed of granular material with height $h$ in a closed box under gravity is quasistatically inclined, and we measure for each bed inclination $\theta$ the corresponding stabilised granular free-surface $\varphi$.

This simple ‘inclined chute’ configuration usually involves two characteristic angles: the start – or avalanche – angle $\theta _{a}$, which is the maximal possible angle attainable with a static granular bed; and the stop angle $\theta _s$, which is the minimal angle needed to sustain flow in the same granular bed (cf. Savage 1979; Daerr & Douady 1999b; Pouliquen & Forterre 2002; Artoni, Santomaso & Canu 2011). In our model, however, we only parameterise the yield transition using a single friction coefficient. As such, we do not expect to recover the full hysteretic phenomenology of such experiments, but rather seek to characterise our Drucker–Prager rheology from a macroscopic point of view, and to provide a way to set the friction coefficient $\mu$ in the context of granular collapses.

We consider in practice a 2-D (or 3-D) box of dimensions $1~\textrm {m} \times 0.3~\textrm {m}$ ($\times 0.1~\textrm {m}$), discretised with a Cartesian mesh with resolution $250 \times 150$ (resp. $100 \times 30 \times 10$), filled with a 0.1 m-high bed of granular material with slip boundary conditions on the sides and stick boundary conditions at the bottom of the box. To minimise the effect of boundary conditions, we measure the angle of the free surface between $x = 0.3$ m and $0.7$ m using a $\tan ^{-1}$ on the end-point coordinates of a line fit. The simulation time-step is chosen as $\delta t = 1$ ms and we increase $\theta$ by steps of $0.5^\circ$ ensuring equilibrium is always reached between each change of inclination. Note that the simulated box is actually kept fixed and we vary the inclination $\theta$ by rotating the external gravity force.

Figure 1 shows the results for the free surface angle $\varphi$, once the surface has stabilised, as a function of $\theta$ for a simulation with $\mu =0.44$. As we increase the bed inclination $\theta$, we observe that the free surface first remains aligned with the bed ($\varphi = \theta$); the material is jammed and remains in static equilibrium. After some critical inclination $\theta _a$ (the so-called avalanche angle), the material starts to flow before stopping. If we again increase quasistatically the bed inclination, the same happens again, and the angle of the stabilised free surface basically takes a constant value, corresponding to the so-called stop angle $\varphi _s = \theta _s$ for $\theta > \theta _a$.

Figure 1. Inclination $\varphi$ of the stabilised granular bed free surface as a function of the box inclination $\theta$ for the avalanche numerical set-up obtained in our 2-D and 3-D simulations. (a) Simulation snapshots of the 2-D granular bed for $\theta =23^\circ$ and $27^\circ$. (b) Plot of $\varphi$ versus $\theta$; the vertical and horizontal dashed lines indicate the friction angle used in the simulation $\tan ^{-1}(0.44) \approx 23.75^{\circ }$.

As expected from our simple non-hysteretic Drucker–Prager model, we observe that the avalanche and stop angles are indistinguishable. Furthermore, as shown in figure 1(b), both angles correspond precisely to the friction angle prescribed by the friction coefficient $\mu$ of the Drucker–Prager law, as $\theta _a = \theta _s = \theta _{\mu } = \tan ^{-1}(\mu )$.

On the one hand, these two observations demonstrate that our numerical model provides a consistent discretisation of the Drucker–Prager law, and gives us confidence in using the corresponding Sand6 code for further experiments. On the other hand, this elementary numerical experiment makes it clear that modelling granular material with a simple plastic Drucker–Prager rheology requires us to be careful in the choice of the – constant – friction coefficient.

In the following, we show that our numerical model, which features a single friction parameter, is sufficient to predict the macroscopic flowing of a real granular collapse, provided the friction coefficient is chosen as the avalanche angle $\theta _a$ of the real granular material. This observation, which departs from the usual choice for setting the friction parameter (the stop angle), is supported by complementary validations and comparisons of our model throughout the paper, and eventually discussed.

3. Experimental collapses and measure of macroscopic parameters

In order to validate our non-smooth simulator in transient flow configurations, we have conducted various experiments of granular column collapses in a 6 cm-wide channel. Our experiments are performed on various inclinations of the channel bed and using two different materials, which we describe in the following. To relate our physical experiments to our numerical simulations, we also need to measure appropriately, from the experiments, all the – macroscopic – parameters at play in the continuum model: in particular the yield friction coefficient $\mu$ of the Drucker–Prager rheology, but also the friction coefficient between the granular medium and the lateral walls. These measurement processes are discussed in the next section.

3.1. Materials

We consider two different materials: 0.5 mm glass beads (Sigmund Linder – SiliBeads – 45015) and 2.7 mm natural granules. Their aspects are shown in figure 2, and their physical parameters – which we have experimentally measured – are summarised in table 1.

Figure 2. Photograph of the two granular materials used for our validation experiments: glass beads on the left, and natural irregular granules on the right. Multicoloured material improves feature detection, hence velocimetry performance.

Table 1. Measured material and rheological parameters for the two granular materials used in our experiments: $d$ is the median grain diameter, $\rho _c$ the density of the material in the dense state, $\mu _a$ the avalanche friction coefficient, $\mu _{stop}$ the stop friction coefficient and $\mu _w$ the estimated friction coefficient between the granular material and the lateral walls. Note that the density $\rho _c$ is measured on the granular media in its dense state and not on the composing grain, and corresponds as such to a direct measurement of $\rho _c = \rho _g \phi _c$.

Our 0.5 mm glass beads share the same geometric and physical features as the glass beads used in Pouliquen (1999), Pouliquen & Forterre (2002) and Jop et al. (2005). The latter references represent seminal experimental works for calibrating the $\mu (I)$ constitutive law of Jop et al. (2006), relying on a steady granular flow in a wide channel with various inclinations. In § 3.3, we verify experimentally that our beads behave in accordance with the results reported in the aforementioned literature, when subjected to a similar setting. We note that our beads are also similar to those used for the granular collapses studied in Mangeney-Castelnau et al. (2005) and Farin et al. (2014), whose experimental results are exploited to test modelling assumptions formulated by Ionescu et al. (2015) and later by Martin et al. (2017). Aligning our choice on this material thus makes it possible to check such assumptions and better transfer insights. For this reason we have used this well-calibrated model material preferably for extensive and systematic validation, as reported in § 4.2.

In order to test the robustness of our model on natural materials – generally exhibiting larger friction coefficients than spherical glass beads, we have also used coarse natural granules (2.7 mm irregular grains). The extension of our validation study to the collapse of such natural granules is shown in § 4.2.

3.2. Collapse set-up

The experimental apparatus is inspired by Balmforth & Kerswell (2005) and sketched in figure 3. An initial granular pile with aspect ratio $a=H_0/L_0\sim 0.5$, where $H_0$ is the depth and $L_0$ is the length along the channel, is confined in a $W = 6$ cm-wide flume with a lifting gate preventing the pile collapsing. The experiment consists in opening the pneumatic lifting gate to trigger the granular collapse. Similarly to the second set of experiments of Farin et al. (2014), the base of the box is covered with a horizontal erodible granular bed of ${\sim }2$ cm height made of the same material.

Figure 3. (a) Sketch of our granular column-collapse experimental configuration showing the initial column dimension ($L_0$,$H_0$) and final deposit state (run-out $L_f$ and final upslope height $H_f$). Collapses are triggered by a pneumatic lifting gate. (b) Photograph of the final deposit.

Free surface profiles as well as velocities are collected using a high-speed camera from the transparent sidewalls with an interframing time of 5.3 ms. We estimate the local optical flow of good feature to track (Shi 1994; Miozzi, Jacob & Olivieri 2008), i.e. features having optimal contrasts to be tracked (enhanced by the presence of multicoloured materials as shown in figure 2). We interpolate feature velocities on a mesh to produce the velocity fields using the opyf Python package. The reader may refer to Annex 1 of Rousseau & Ancey (2020) for details on the velocimetry procedure.

Figure 4 illustrates the methodological data used to compare experiments and simulations in this article: figure 4(a,c) shows snapshots of the recorded experimental collapse for the initial and some intermediate state, coloured with the velocity magnitude as measured by velocimetry, while figure 4(b,d) shows the corresponding snapshots of the simulation, where the background medium and velocity are displayed by the material points quantities. We superpose on all these snapshots the postprocessed profile height and static–flowing transition lines determined as the $0.5$ contour on the volume fraction $\phi$ and $0.01\ \textrm {m}\ \textrm {s}^{-1}$ contour on the velocity, respectively. Figure 4(e) illustrates the direct comparison of the height profiles in the final state.

Figure 4. Capturing the free surface and velocity field from experimental and simulated granular collapses at different instants. Camera snapshots are shown in (a,c) and (e) where experimental velocities are extracted using image velocimetry and shown on the snapshots foreground. The static–flowing transition corresponds to the 0.01 m s$^{-1}$ contour of the experimental and numerical velocity fields. The simulated material points are shown in (b,d) and coloured according to their velocity magnitude. The black rectangle in (b) has dimensions $5 \delta _x \times 5 \delta _z$ with $\delta _x$ and $\delta _z$ the respective horizontal and vertical MPM resolutions. Data correspond to the B05 collapse (0.5 mm beads – $i=5^\circ$). See movie 1 in the online Supplementary material and movies are available at https://doi.org/10.1017/jfm.2023.835 for a visualisation of the experimental and simulated evolution in time.

3.3. Measurement of macroscopic parameters

As mentioned in § 2.4, the flowing threshold in our Drucker–Prager constitutive law is determined by a single friction coefficient $\mu = \tan \theta$. Yet, there is still no clear consensus regarding which of the avalanche angle $\theta _a$ and the stop angle $\theta _s$ should be used as the yield angle. In order to discriminate between these two options, we first need to measure both parameters for the granular media used in our experiments. Since the role of sidewalls is also investigated, we also evaluate the friction coefficient between the granular medium and the lateral walls $\mu _w$.

3.3.1. Measurement of the avalanche friction coefficient

Our first experiment for measuring the avalanche angle $\theta _a$, is based on the configuration used by Balmforth & Kerswell (2005) for measuring the ‘internal’ angle of friction – a terminology which is then equivalent to our avalanche angle from an experimental point of view. Similarly to the numerical set-up presented in § 2.4, we consider a 20 cm-wide flume filled with the granular material and inclined the flat granular bed until motion downslope begins (see figure 5a). Our second experiment consists in pouring slowly granular material upon a conical heap (see figure 5b). In both experiments, we measure $\theta _a$ as the critical angle corresponding to the maximal slope of the material, just before an avalanche occurs and causes the slope to decrease (see movie 2 in the online Supplementary material and movies for an illustration of the process). As emphasised by Balmforth & Kerswell (2005) or Russell et al. (2019), there is a large spread in the $\theta _{a}$ measurements owing to the sensitivity of the protocol to small perturbations.

Figure 5. The avalanche angle $\theta _a$, used to set the yield friction coefficient $\mu$ of the Drucker–Prager model as $\mu = \mu _a = \tan \theta _a$ using either the large channel inclination test (a), or the conical heap test (b). For our 0.5 mm glass beads depicted here, both protocols converged to a same average measurement value ($\tan \theta _a = \tan 23.7^\circ = 0.44 \pm 0.03$).

For the 0.5 mm glass beads material, no matter the experimental protocol used (inclined channel or conical heap), we found the avalanche angle $\theta _{a} = 23.7^\circ \pm 1$, which corresponds to the friction coefficient $\mu _a = \tan \theta _{a} \approx 0.44 \pm 0.03$. This value is consistent with those made in the literature for similar materials, such as the 0.7 mm or 0.8 mm glass beads of, respectively, Farin et al. (2014) and Balmforth & Kerswell (2005) which give $\tan \theta _a = \tan 25^\circ = 0.47 \pm 0.01$ and $\tan \theta _a = \tan 24.5^\circ = 0.46 \pm 0.04$. For the 2.7 mm granules, we have only conducted the experiment with the first protocol, and found $\theta _{a} = 36.8^\circ \pm 5$ ($\mu _a = 0.75 \pm 0.1$).

3.3.2. Measurement of the stop friction coefficient

As mentioned above, another feature of granular flows is the existence of a minimal angle necessary to sustain motion in an already flowing material. In the case of a granular layer with constant height flowing on rough inclines, this stop angle depends on the height $h$ as a function $\varTheta (h)$ (Pouliquen 1999; Pouliquen & Forterre 2002), which naturally converges towards a constant value as the height increases, namely $\theta _s = \varTheta _s(h \to \infty )$.

This phenomenology has led to nowadays-well-established protocols to measure $\theta _s$ in a steady granular flow, in which either the bed with steady flow of known height is progressively lowered until flow ceases (Pouliquen 1999), or the initially static bed is inclined up until it starts to flow, and then stops (Pouliquen & Forterre 2002).

We have replicated the steady inclined bed protocol experimentally with our 0.5 mm beads and 2.7 mm granules, and have measured the resulting stop angle as in Pouliquen & Forterre (2002) using a 3 m-long, 10 cm-wide channel ($W/d = 200$) for beads and a 2 m-long, 12 cm-wide channel for granules ($W/d = 44$).

While our channels are narrower than the one used in Pouliquen & Forterre (2002), the maximal characteristic flowing heights $h$ are systematically small with respect to the channel $W$, with $h/W \sim 0.06$ for both beads and granules. Following Jop et al. (2005), sidewall friction can thus result in a small overestimation of the stop friction coefficient of order $\mu _w h/W \sim 0.02$ in both cases.

Our results for beads are shown in figure 6(a), together with the experimental data and corresponding fits of Pouliquen (1999), Pouliquen & Forterre (2002), Forterre & Pouliquen (2003) and Farin et al. (2014). We observe that despite some spreading in all the collected data, which might be attributed to the different experimental protocols used (Pouliquen & Forterre (2002) versus Pouliquen (1999)) and/or to the slightly different beads used (Pouliquen & Forterre (2002) versus Farin et al. (2014)), we obtain a good agreement with Pouliquen & Forterre (2002), in particular regarding the $h/d_p \to +\infty$ asymptote, which corresponds to the $I \to 0$ stop angle.

Figure 6. (a) Variation of the stop angles $\varTheta _s$ as a function of the non-dimensional thickness $h/d$ for our glass beads (B) of diameter $d=0.5$ mm, measured using the protocol described in Pouliquen & Forterre (2002). The plain curve corresponds to the best fit of our data with the function $\tan (\varTheta _s(h)) = \mu _1 + (\mu _2 - \mu _1)/(1+h/L)$, obtained for $\mu _1 = 0.37$, $\mu _2 = 0.54$ and $L = 1.16$ mm. We also report the data and fits of Pouliquen & Forterre (2002), Pouliquen (1999), Forterre & Pouliquen (2003) and Farin et al. (2014) for the sake of comparison. Curves do not all collapse perfectly, which may be explained by slight differences between experimental conditions, protocols or materials. Overall, our data look consistent with the previous studies performed with 0.5 mm glass beads. (b) Same measurements made for our granules (G) with a best fit obtained for $\mu _1 = 0.65$, $\mu _2 = 0.95$ and $L = 3.1$ mm.

Indeed, we retrieve a stop friction coefficient $\mu _{stop} = \tan (\varTheta _s(\infty )) = \mu _1 = 0.37$ which is very close to the values for the same type of beads, $\mu _{stop} = 0.38$ (Pouliquen 1999; Pouliquen & Forterre 2002; Jop et al. 2006). Accounting for the $10$ cm width of our channel, we thus deduce that our actual stop friction coefficient is slightly less than $0.37$. In practice, for our comparisons we shall stick to $\mu _{stop} = 0.38$ (see § 4.3).

Following the same measurement protocol with our granules, we find a stop friction coefficient of $\mu _{stop} = 0.65 \pm 0.1$ for this natural sediment (see table 1 and figure 6b).

These results confirm that our beads behave fairly similarly to the ones from Pouliquen & Forterre (2002), which enables direct comparison of our results with other works.

It is noteworthy that whatever the material used, the estimated stop friction coefficient ($\mu _{stop} = 0.38$ and $\mu _{stop} = 0.65$ for the beads and granules, respectively) is significantly lower compared with the measured avalanche coefficient ($\mu _a = 0.44$ and $\mu _a = 0.75$, respectively). In § 4, we find that using $\mu = \mu _a$ in our simulations yields more realistic collapses than taking $\mu = \mu _{stop}$.

3.3.3. Measurement of the friction coefficient between the granular medium and the lateral walls

Following Balmforth & Kerswell (2005) and Hutter & Koch (1991), we estimate the friction coefficient $\mu _w$ between the glass walls and the granular materials by measuring the critical angle above which a block of particles held together within a rigid light plastic cylinder begins to slide on an inclined glass surface. The cylinder is 4 cm in diameter and 6 cm high. We obtain $\mu _w = \tan (13^\circ ) = 0.23$ for the $0.5$ mm glass beads and $\mu _w = \tan (17^\circ ) = 0.3$ for the $2.7$ mm granules. These values are reported in table 1 along with the corresponding estimated standard deviations.

4. Validation: numerics versus experiments

To evaluate the performance of our simulation approach and assess the validity of the Drucker–Prager rheological law – featuring only one friction coefficient – in the case of transient collapse flows, we perform direct 3-D simulations of the collapse set-ups presented in § 3.2, with a full account of frictional interactions with the encasing walls and the lifting door.

We first compare the simulation results with the experiments for the glass beads (B), in order to validate our numerical model with a controlled material, and then highlight its robustness by considering the case of natural irregular granules (G) (cf. figure 2).

4.1. Simulation parameters

The rheological parameters used in our simulations are derived from independent experimental measurements performed with the same materials. As discussed in § 3.3.1, the Drucker–Prager friction coefficients are determined using the respective avalanche angles, measured consistently with both the channel and heap set-ups. To reduce uncertainties on the actual granular dense packing fraction, the density is directly measured for the granular media in the dense state, and thus corresponds to the critical density $\rho _c = \rho _g \phi _c$. This amounts to rescaling the volume fraction field $\phi$ by its dense critical value $\phi _c$, so as to ensure $\rho _c ({\phi }/{\phi _c}) = \rho _g \phi = \rho$. As such, we take in the simulations $\phi _c = 1$. The values for the material parameters and the friction coefficients are summarised in table 1.

Gravity is set to its standard value $g = 9.81$ m s$^{-2}$, and aligned with a rotated downward unit vector to account for the various respective inclinations of the collapse beds.

The $152 \times 15 \times 6$ cm$^3$ simulation domain is discretised with a $152 \times 30 \times 6$ mesh of regular rectangular cuboids, which gives a $1$ cm length and width resolution, and a $0.5$ cm resolution for the height, in order to better capture the free surface. Figure 4(b) illustrates the resolution of the background MPM-grid: the black rectangle has dimensions $5 \delta _x \times 5 \delta _z$ where $\delta _x$ and $\delta _z$ are the horizontal and vertical resolutions, respectively.

We impose frictional boundary conditions with the lateral walls and the left upstream wall – the bottom wall being covered with a $2$ cm bed of erodible granular material to ensure the same bulk friction conditions. As mentioned in § 3.2, the experimental collapses are initiated by the lifting of a pneumatic gate, which can interact frictionally with the granular material and affect the flow in the initial phase. We therefore directly account for the lifting of the gate by introducing a rigid plane with imposed motion mimicking the gate. This ‘numerical’ door also interacts frictionally with the material, with the same friction coefficient $\mu _w$ as the lateral walls. Figure 7 illustrates the numerical set-up with the lifting of the gate. Note that while the upward lifting motion of the door influences the early dynamics, frictional interaction between the material and the door plays a negligible role on the flow and the final deposit, as discussed in § B.2.

Figure 7. The 3-D MPM simulations of the B15 granular column collapses over a $15^\circ$-tilted bed at three different times. The motion due to the lifting gate is particularly visible at $t=0.2\ \textrm{s}$. A 3-D animated visualisation is provided in the Supplementary material and movies (cf. movie 3).

The simulations are run with approximately $216$ MPM particles per mesh cell, with random initial positions, and a time-step of $8.3 \times 10^{-4}\ \textrm {s}$.

4.2. Validation on various materials for different inclinations

Figure 8 shows a comparison between profiles for the collapse experiments and simulations performed on our 0.5 mm glass beads for five bed inclinations, ranging from ${0}^{\circ }$ to $20^{\circ }$. Profiles are plotted at three different times, to highlight the behaviour of the collapse in the starting, flowing and stopped phases. Although numerical profiles presented a negligible difference along the flume width ($\kern0.7pt y$-axis), they were extracted next to the sidewall position to faithfully compare with measurements. The right-hand column shows the final state of the deposit, when the mass has stopped moving, i.e. when we detect no velocity above $0.01\ \textrm {m}\ \textrm {s}^{-1}$ in the whole domain, with corresponding respective experimental and simulation final – stop – times $t_{f,exp}$ and $t_{f,sim}$ which naturally depends on the bed inclination, as summarised in table 2. Note that the differences observed for the final time $t_f$ between the experiment and the simulation are mainly caused by the local nature of the ‘rest’ criterion we use (no velocity above $0.01\ \textrm {m}\ \textrm {s}^{-1}$ in the whole domain), which artificially delays the final time of the experiment due to marginal grain movement at the free surface. We also plot the experimental and simulation static–flowing transition lines determined as the $0.01\ \textrm {m}\ \textrm {s}^{-1}$ contour of the velocity field in order to highlight the bulk dynamics of the collapse.

Figure 8. Comparison between collapse experiments with the 0.5 mm glass beads and 3-D simulations for various bed inclinations ranging from $0^{\circ }$ (a) to $20^{\circ }$ (e) (B00, B05, B10, B15 and B20 runs). The profiles are extracted next to the sidewall position. We compare both the thickness profiles (solid pink line for the experiment, dash–dotted black line for the simulation) and the static–flowing transition (pink dashed line for the experiment, black dotted line for the simulation). The velocity heatmap in the background is computed from the simulation. In the simulations, the constant friction $\mu$ was set to $\mu _a = 0.44$, and the wall friction to $\mu _w = 0.23$.

Table 2. Metrics of the experimental and numerical runs for the beads. The experimental and simulation rest times ($t_{f,{exp}}$ and $t_{f,{sim}}$) correspond to the times when we detect no velocity above 0.01 m s$^{-1}$ in the domain. Here $H_0$ is the initial pile height and $a$ the initial collapse aspect ratio; $H_{f,{exp}}$ and $H_{f,{sim}}$ are the experimental and simulated final pile height while $L_{f,{exp}}-L_0$ and $L_{f,{mod}}-L_0$ are the experimental and simulated run-out distances measured from the gate position ($L_{f}$ is defined as the final run-out distance from the left collapse wall). The run-out distance is defined as the length of the continuous deposit extent which height is above a 2 mm threshold, following Balmforth & Kerswell (2005). Note that the missing final experimental run-out distance for the B20 collapse is due to restrictions on the field of view of the camera.

We first observe an excellent agreement between the experimental and computed final thickness profiles for all the inclinations, highlighting the ability of our numerical method to capture threshold effects and large strains. The $t=0.2$ s snapshots also highlight the role of simulating the lifting of the door, which influences the initial dynamics, as also observed by Ionescu et al. (2015).

The dynamics of the collapses at early steps also appears to be well described by the simulation for small inclinations. Above $15^{\circ }$ inclination, we observe some difference in the thickness profiles, with in particular a systematic under-estimation of the profile height near the left-hand wall, which was also observed by Martin et al. (2017). The static–flowing transition lines follow the same overall behaviour, with a fair agreement for small inclinations, degrading for larger ones.

We also show the position and velocity of the granular front as a function of time in figure 9. We can first observe that our simulations accurately capture the stopping time of the front, with a rather good prediction of the front velocity for each collapse inclination. Note that the front always stops before the final stopping time $t_f$ of the whole flow.

Figure 9. Front position (a) and velocity (b) for the experimental and simulated granular collapses, for bed inclinations from $0^\circ$ (top) to $15^\circ$. Note that the run at $20^\circ$ inclination was left out since our experimental recording set-up could not frame the run-out during the whole collapse. The 3-D simulations were performed with the constant friction $\mu = 0.44$, and the wall friction $\mu _w = 0.23$.

Overall, while the early dynamics of the collapses seem to depend on the inclination, these results validate our simulation approach, and in particular support the use of a simple Drucker–Prager law with the constant avalanche friction coefficient to capture final deposit profiles.

The robustness of these results has been evaluated using a more natural and irregular material (denoted granules above, cf. figure 2 and table 1). The corresponding profiles are presented in figure 10, and the final pile heights and run-out values in table 3. We again observe a good agreement between the experiments and simulations.

Figure 10. Comparison between collapse experiments with the natural granules and 3-D simulations for various bed inclinations ranging from $0^\circ$ (a) to $15^\circ$ (d) (G00, G05, G10 and G15 runs). The profiles are extracted next to the sidewall position. We compare both the thickness profiles (solid pink line for the experiment, dash–dotted black line for the simulation) and the static–flowing transition (pink dashed line for the experiment, black dotted line for the simulation). The velocity heatmap in the background is computed from the simulation. In the simulations, the constant friction $\mu$ was set to $\mu _a = 0.75$, and the wall friction to $\mu _w=0.3$. See movie 4 in the Supplementary material and movies for an animated evolution of the G15 simulated collapse against experiment.

Table 3. Metrics of the experimental and numerical runs for the granules. The letter $G$ indicates granules. Please refer to the definition of the quantities in table 2.

The small discrepancies observed close to the lifting gate at the beginning of the flow suggest that the frictional interaction of the material with the door is higher. We also notice, as was also the case for the glass beads, that the material flows less in the early dynamics, suggesting that more subtle rheological effects are at play in this stage. However, again this does not impact the final run-outs and profiles.

We should stress that the measurement protocols for the avalanche and stop friction coefficients are subjected to important uncertainties (as reported in table 1) in the case of such gritty material, due to the inherently large inclination angles required to reach the yield and stop transitions, and to the increased impact of boundary effects in the resulting very inclined channels. As such, we do not try to discriminate accurately against the use of either the avalanche or stop friction coefficient, but rather interestingly illustrate the reliability of the plain Drucker–Prager rheology in the context of transient collapse flows of natural materials.

4.3. Comparison with the $\mu (I)$ rheology

The $\mu (I)$ rheology, initially proposed by Jop et al. (2006) in the context of steady inertial flows, introduces a local dependence of the Drucker–Prager coefficient on the inertial number

(4.1)\begin{equation} I = \frac{d \lVert \dot{\varepsilon} \rVert}{\sqrt{p/\rho}},\end{equation}

with $d$ the particle diameter, $\dot {\varepsilon }$ the strain-rate tensor introduced in (2.1), $p$ the pressure and $\rho$ the material density. This inertial number characterises the ratio between the inertial time scale $d \sqrt {\rho / p}$ and the deformation time scale $1/\lVert \dot {\varepsilon } \rVert$.

The $\mu (I)$ rheology then writes as

(4.2)\begin{equation} \mu(I) = \mu_{stop} + \frac{\mu_2 - \mu_{stop}}{1 + \dfrac{I_0}{I}},\end{equation}

where $\mu _{stop}$, $\mu _2$ and $I_0$ are material-dependent coefficients which can be calibrated using the steady inclined plane experiment as described in Jop et al. (2005).

This model has recently been used with success for transient collapse prediction in several studies (Lagrée et al. 2011; Ionescu et al. 2015; Martin et al. 2017).

The $\mu (I)$ rheology equation (4.2) can be straightforwardly implemented in our simulation framework by adding an explicit friction coefficient update in the yield criterion and solving for the Drucker–Prager law with this new coefficient as before. As discussed in § 3.3, we have performed the steady inclined flow protocol with our $0.5$ mm beads, which are similar to the one from Jop et al. (2005), and extracted the $\mu _{stop}$ and $\mu _2$ coefficients from the fit of $\varTheta _s(h/d)$ (cf. figure 6). As expected, we retrieve values close to the ones from Jop et al. (2005); for the sake of comparison, and owing to the materials similarities, we have not determined the characteristic inertial number $I_0$, and have chosen to use the parameters from Jop et al. (2005) to perform our $\mu (I)$ simulations: $\mu _{stop} = 0.38$, $\mu _2 = 0.64$ and $I_0 = 0.279$.

Figure 11 shows a comparison of the experimental and computed collapse profiles with the simple Drucker–Prager ($\mu _a=0.44$) and the $\mu (I) (\mu _{stop}=0.38)$ rheologies. For completeness, we also show the numerical results obtained with a constant $\mu _{stop}$ friction coefficient – the lower bound of $\mu (I)$ – in figure 12. For the sake of readability, we focus here on the results for the $15^\circ$-inclined configuration; the corresponding comparison figures are provided for all inclinations in § B.1, and support the sameobservations.

Figure 11. Comparison of 3-D simulations performed with $\mu = \mu _a = 0.44$ (black lines) and $\mu = \mu (I)$ (green lines) for the $15^\circ$ bead collapse (B15) at $t=0.2$ s, $0.6$ s and $t_f$. Experimental curves (solid pink line) are provided for reference, and we show both the free-surface lines (resp. plain and dash–dotted) and the static–flowing transition contours (resp. dashed and dotted). The results for the other collapse inclinations can be found in § B.1.

Figure 12. Comparison of 3-D simulations performed with $\mu = \mu _a = 0.44$ (black lines), $\mu =\mu (I)$ and $\mu = \mu _{stop} = 0.38$ (blue lines) for the $15^\circ$ bead collapse (B15) at $t=0.2$ s, $0.6$ s and $t_f$. Experimental curves (solid pink line) are provided for reference, and we show both the free-surface lines (resp. plain and dash–dotted) and the static–flowing transition contours (resp. dashed and dotted). The results for the other collapse inclinations can be found in § B.1.

We observe that both the $\mu (I)$ and the constant $\mu _{stop}$ rheologies fail to capture accurately the free surface and the static–flowing transition region, in both the dynamical and final resting phases. More precisely, they behave very similarly and both underestimate the flowing threshold, leading to an increased flowability and broaderdeposits.

The similarity observed between the $\mu (I)$ and constant $\mu _{stop}$ rheologies is coherent with the systematic low measured inertial numbers $I$, as shown in figure 13: except for marginal surface points, the inertial numbers are of the order of at most $I \sim 10^{-2}$ during all the collapse, which is much lower than the characteristic value $I_0 = 0.279$ of the $\mu (I)$ rheology. The additional viscous dissipation introduced by the $\mu (I)$ law compared with a constant $\mu _{stop}$ is therefore negligible for such collapses, as was also observed by Valette et al. (2019) for horizontal column collapses with a regularised $\mu (I)$rheology.

Figure 13. Inertial number $I$ during the simulated $15^\circ$ bead collapse (B15) with $\mu =\mu _a=0.44$: (a) evolution of $I$ as a function of time. Median, $30$- and $70$-quantiles are estimated in the regions where $I \geq 5 \times 10^{-4}$ thereby filtering out static parts of the collapse; (b) instantaneous inertial number heatmap at $t_{I_{70,max}} = 1.1$ s, i.e. the time for which the $70$-quantile of $I$ is maximal. The $I$ values remain significantly lower than $I_0 \sim 0.3$ during our collapses explaining why the increase of the friction coefficient from $\mu _{stop}$ to $\mu _2$ in the $\mu (I)$ rheology equation (4.2) has a negligible effect. (See movie 5 in the Supplementary material and movies for an animated evolution of the inertial number heatmap.)

The necessary use of a friction coefficient larger than the ‘stop’ friction coefficient $\mu _{stop}$ measured from steady inclined plane experiments to quantitatively match experiments was also discussed by Lagrée et al. (2011), or Ionescu et al. (2015) and Martin et al. (2017). The last two explain their increased friction coefficient by the role of the lateral walls, an effect that we shall, however, discard below (cf. § 5).

4.4. Experimental validation of wall friction effect

The validation of the numerical method is ended by evaluating the simulation of the lateral walls as such boundary conditions are crucial when modelling 3-D configurations (Jop et al. 2005; Lajeunesse et al. 2005; Ionescu et al. 2015; Martin et al. 2017). For that purpose, we perform dedicated collapse experiments with two different channel widths, and compare the corresponding height profiles with the simulated ones. For simplicity, the granular collapses are performed in a rough horizontal – $0^\circ$ – channel, with the 0.5 mm glass beads described above. We consider two different channel widths: $W = 1$ cm and $W = 4$ cm, which thus correspond to $W/d = 20$ and $80$, respectively. Experimental conditions are similar to the aforementioned $0^\circ$ collapse, but with initial column dimensions of $H_0 = L_0 = 12$ cm, giving an aspect ratio $a=1$, which allows us to make the effects of the walls more visible on the final upslope height (cf. Balmforth & Kerswell (2005) or Zhang et al. (2021)), thereby providing a simple measure to compare the effect of lateral friction, which we shall use to determine the corresponding 2-D effective friction in § 5.2. The corresponding numerical collapses are carried out consistently, using the constant bulk friction coefficient $\mu = 0.44$ and a friction coefficient between the material and the walls $\mu _w = 0.23$.

Figure 14 shows the rest-state height profiles for the simulations and experiments, with two independent experimental runs for each width in order to ensure reliability of the results. We observe that the simulated profiles compare noticeably well with the experimental ones, which underlines the validity of our numerical method, and indicates that a plain Drucker–Prager rheology, with – independently measured – friction coefficients $\mu$ and $\mu _w$ can accurately account for confinement effects in such transient flows. Note that the good agreement between experiments and simulation observed in the narrow case $W/d = 20$ also suggests that continuum modelling remains valid in this range, at least for such collapse flows.

Figure 14. Influence of the sidewalls: rest-state height profiles of simulated and experimental collapses for $a=1$ and two different widths ($W=1$ cm and $W=4$ cm) at $0^\circ$ slope. Two experimental replicates (A and B) are compared for each width run. The green solid line delineates the initial column with size $0.12 \times 0.12\ \textrm {cm}^2$. Numerical collapses are performed with a constant internal friction coefficient $\mu =0.44$ and a friction coefficient between the glass walls and the material $\mu _w=0.23$.

5. Revisiting features of granular collapses with 3-D simulations in light of a constant friction coefficient

As demonstrated in § 4, the plain Drucker–Prager rheology, with a constant friction coefficient set from the avalanche angle of the material, can accurately predict experimental column collapses on inclined beds, for a large range of inclinations. Our validated 3-D non-smooth numerical model therefore enables us to explore and revisit phenomenological features of transient collapses, which have been previously interpreted in the framework of a steady viscoplastic rheology. While viscous rheological contributions – as introduced by the $\mu (I)$ rheology – are undeniably crucial to represent inertial granular flows, notably to explain the existence of steady inclined flows, the range of applicability of the fully plastic Drucker–Prager rheology has, to the best of our knowledge, never been thoroughly framed, in particular in the context of transient flows at low strain-rates. It thus appears essential to elucidate the different ingredients responsible for the various features of granular collapses, in order to disentangle the respective roles of frictional or viscous effects, and clarify the links between steady and transient geometries in 3-D configurations. In the following, we first investigate this question by referring to established scaling laws relating the final pile run-out and height to the initial column aspect ratio (Balmforth & Kerswell 2005; Lajeunesse et al. 2005; Lube et al. 2005; Lagrée et al. 2011; Dunatunga & Kamrin 2015). We then propose a comprehensive analysis of the impact of sidewalls in granular collapses, in light of the steady law proposed by Savage (1979), Taberlet et al. (2003) and Jop et al. (2005) and transposed to the transient case by Ionescu et al. (2015).

Our results notably highlight the weak role played by viscosity in transient granular collapses by demonstrating that established scaling laws are well accounted for by the plain Drucker–Prager rheology, and show that the effect of sidewall friction, while possessing a similar scaling behaviour than in the steady case, has much less impact on transient collapses than previously assumed in Ionescu et al. (2015) and Martin et al. (2017).

5.1. Aspect ratio scaling laws in wide and narrow configurations

We perform simulations in wide and narrow channels using different aspect ratios for the initial granular column, and compare the final run-out and upslope height with previously proposed scaling laws. The normalised rest state run-out $L_{f}$ and upslope height $H_{f}$ are shown in figure 15 for $40$ simulated collapses in both wide and narrow configurations with aspect ratios $a = H_0 / L_0$ varying logarithmically from $1$ to $20$. The wide configuration is obtained from 3-D collapses by setting the wall friction to zero (which is equivalent to the 2-D collapse situation). The narrow situation corresponds to a $1$ cm-wide channel with a wall friction coefficient set to $\mu _w=0.23$. The granular material considered is again the glass beads, with a bulk friction coefficient of $\mu =0.44$.

Figure 15. (a) Inverse normalised final upslope height ($H_0/H_{f}$) and (b) normalised run-out ($L_{f}-L_0)/L_0$ as a function of the aspect ratio of the initial column $a$. Here $L_{f}$ is obtained using the position the MPM particle that goes the furthest at the final state.

The power-law fits obtained for our data (summarised in table 4) exhibit scaling behaviours similar to the previous experimental and numerical studies from Lajeunesse et al. (2004), Balmforth & Kerswell (2005), Lagrée et al. (2011) and Dunatunga & Kamrin (2015), and in particular feature the change of regime between the short ($a \lesssim 7$), and high ($a > 7$) columns in the wide configuration. Note that we also observe the theoretical scaling exponents from the shallow 2-D model of Kerswell (2005) and Balmforth & Kerswell (2005), which predicts $H_{0}/H_{f} \sim a^{0.69}$ in our aspect ratio range for the wide configuration and $H_{0}/H_{f} \sim a^{0.5}$ for the narrow one. In comparison with the work of Lagrée et al. (2011) and Dunatunga & Kamrin (2015), our results thus not only confirm the robustness of the power-law scalings observed by Lajeunesse et al. (2005), Balmforth & Kerswell (2005) and Lube et al. (2005), but also highlight that this scaling behaviour of collapse flows is well described by the plain Drucker–Prager rheology with a constant friction coefficient, for a large range of aspect ratios – up to $a\sim 10$.

Table 4. Power-law scaling laws for the run-out and upslope height of the final state of the collapse depending on the initial aspect ratio $a$.

5.2. Role of sidewall friction

A linear relation between the flow thickness and the channel width has been reported by Savage (1979), Taberlet et al. (2003) and Jop et al. (2005) in the context of steady inclined granular flows, stating that the presence of friction with the lateral walls can be accounted for by rescaling the friction coefficient as

(5.1)\begin{equation} \mu_{eq} = \mu + \mu_w \frac{h}{W},\end{equation}

where $W$ is the flow width, $h$ the steady flow thickness and $\mu _w$ the friction coefficient with the walls. This empirical law, which can be interpreted from force balance principles in the steady flow configuration, has also been used in the context of transient granular collapses by Ionescu et al. (2015) to support an increase of the effective 2-D friction coefficient, albeit without proper definition of a flow thickness in such framework. In the following, we use our validated numerical simulator to explore the relevance of such linear scaling in the context of unsteady flows.

5.2.1. Impact of the channel width in transient collapses

We investigate the role of sidewalls friction on collapse dynamics by considering the $15^\circ$-inclined granular collapse case introduced in § 3.2, and perform several 3-D simulations with various channel widths ranging from $W=1$ cm to $W=20$ cm. The results for transient and final states are shown in figure 16. While the early collapse dynamics, outlined by the static–flowing transition contours, seems impacted for narrow channels, the influence of the width on the free surface height appears negligible for $W \geq 6$ cm, and wall friction only starts to play a significant role on the final deposit for very narrow channels $W \lesssim 1$ cm.

Figure 16. Influence of the sidewalls: comparison of 3-D simulations of the $15^\circ$ bead collapse with several channel widths, and a wall friction coefficient $\mu _w = 0.23$.

Note that this weak impact of the channel width on the collapse rest state was also observed in the recent DEM simulations of Zhang et al. (2021), which report negligible impact of the width for $W \gtrsim 20d = 5$ cm for the collapse of granular columns similar to ours.

5.2.2. Equivalent 2-D friction coefficient

To further investigate the linear relation (see (5.1)) in the context of unsteady flows, we perform multiple simulations of 3-D collapses with different channel widths, and determine for each 3-D simulation (with a given width $W$) the best equivalent friction coefficient able to reproduce the corresponding 3-D collapse with a 2-D – unconfined – simulation. Comparison to determine this 2-D equivalent coefficient $\mu _{2\textrm {D},eq}$ is performed on the rest state upslope height $H_{f}$, and as such does admittedly not fully allow us to replace the 3-D simulation by a 2-D one, but still provides insight into the role of lateral confinement and friction.

We consider horizontal and $15^\circ$-inclined collapses in order to test the robustness of the scaling. As mentioned above, we use an initial column aspect ratio of $a=1$ to increase the effect of wall friction on $H_{f}$. In fact, as already observed by for example Balmforth & Kerswell (2005) or Zhang et al. (2021), the upslope part of the column remains static for aspect ratios $a$ below $0.5$ in the horizontal case, so that the final upslope height cannot help assess the role of lateral confinement for such stocky collapses.

Figure 17 gives the results for channel widths ranging from $W=20$ cm down to ${W=0.4}$ cm. Note that we use the initial column height $H_0$ as a length scale, and consider the non-dimensional width $W/H_0$, since the flow thickness $h$ is not relevant in our transient case. To compare our data with the linear relation (5.1), we plot $\mu _{2\textrm {D},eq}$ as a function of the inverse non-dimensional width $H_0/W$.

Figure 17. The 2-D equivalent friction coefficient determined to best match the corresponding 3-D final free-surface, as a function of the channel width, for the $0^\circ$ and $15^\circ$ beads collapses. Our results seem compatible with a linear dependency, which we illustrate by the dashed fit. We also report the linear scaling (pink dotted line) used by Ionescu et al. (2015).

For the sake of comparison, we also show the linear scaling hypothesis used by Ionescu et al. (2015) to obtain the 2-D friction coefficient corresponding to the 3-D $\mu (I)$ law with an effective flowing thickness at $h \sim 0.05$ m estimated using the maximum flowing thickness for collapses with $H_0 \sim 0.1$ m (similar to our collapses). This hypothesis was used to account for the confinement effects in collapses within 2-D simulations, and typically led to a rescaling of $\mu _1$ from $0.38$ to $0.38 + \mu _w 0.05/W \simeq 0.48$ for $10$ cm-wide collapses.

We can, however, observe that our 3-D simulations do not support this hypothesis, and exhibit a significantly weaker impact of the lateral walls, as was also suggested in § 5.2.1. Our equivalent 2-D friction coefficient nevertheless still appears linearly dependent on $H_0/W$, with fitted coefficients corresponding to

(5.2)\begin{equation} \mu_{2{\rm D},eq} = 0.44 + 1.2 \times 10^{{-}2} \frac{H_0}{W}. \end{equation}

Interpreted from (5.1), it would correspond to an effective flow height $h \simeq 1.2 \times 10^{-2} H_0 / \mu _w \simeq 6.3$ mm. This value, which is approximately one order of magnitude lower than the maximum flowing thickness used by Ionescu et al. (2015) to estimate an effective flow thickness, illustrates that special care is required to estimate an effective flow thickness in order to deduce the 2-D equivalent friction coefficient.

6. May a hysteresis phenomenology explain the collapse onset behaviour?

Our non-smooth granular model can accurately predict the material profiles of granular collapses using a single friction coefficient inferred from measured avalanche angles. While generally consistent with experimental measures, the simulated flow dynamics, however, exhibits significant differences in the early stages of the collapse, as illustrated by the static–flowing transition contours in figure 8. This early deviation, also observed using the $\mu (I)$ rheology by Martin et al. (2017), suggests that more complex rheological effects are at play during the onset of the flow, which naturally involves low velocity and hysteresis phenomenology (Pouliquen & Forterre 2002). As a first step to highlight the role of a static friction coefficient and a potential hysteresis between the solid and flowing phases, we implement a non-constant friction coefficient law, somehow reminiscent of the well-known static–dynamic transition in solid friction,

(6.1)\begin{equation} \mu_{hyst} (I) = \left\{ \begin{array}{@{}ll@{}} \mu + {\rm \Delta}\mu_{hyst} \left( 1 - \dfrac{I}{I_*} \right) & \textrm{if } I < I_{*}\\ \mu & \textrm{if } I \geq I_{*} \end{array}\right., \end{equation}

and illustrated in figure 18.

Figure 18. Friction law as a function of the inertial number $I$. The dash–dotted red curve represents a non-monotonic friction law as expected for granular material (DeGiuli, McElwaine & Wyart 2016). The blue solid line is the simplified linear law implemented with a hysteresis gap ${\rm \Delta} \mu _{hyst}$.

The friction coefficient $\mu (I=0)$, which describes the static–flowing transition, is then $\mu + {\rm \Delta} \mu _{hyst}$, and $I_{*}$ characterises the transition from ‘static’ to ‘dynamic’ friction. While simpler than the non-monotonic law discussed by DeGiuli & Wyart (2017), our two-valued law still allows for hysteretic instability due to the decrease of the friction coefficient between $I=0$ and $I_*$, and is chosen in this context to illustrate the potential effects of such phenomenology on transient collapses, paving the way to more elaborate analyses. Note that the increase of the effective friction coefficient at larger $I$ due to collisional dissipation, as predicted by the $\mu (I)$ rheology is also not accounted for in (6.1) as we consider only low inertial number flows.

Figure 19 shows the results for the $15^\circ$-inclined collapse, using the hysteretic law (6.1) with $\mu =\mu _{stop}=0.38$, ${\rm \Delta} \mu _{hyst}=0.16$ and a transition inertial number $I_{*} = 4 \times 10^{-2}$, consistent with the study of DeGiuli & Wyart (2017). Despite the simplicity of the law, we observe that the increase of the friction coefficient for very low inertial numbers, characterised by the static friction coefficient $\mu + {\rm \Delta} \mu _{hyst} = 0.54$, definitely improves the prediction of the flow dynamics at early stages, while only weakly affecting the final rest state.

Figure 19. Hysteresis effect on the $15^\circ$-collapse (B15): simulation results obtained with the hysteretic law equation (6.1) against experimental profiles and static–flowing transition contours. The non-hysteretic simulation results from figure 8 are recalled for reference, and the hysteretic simulation corresponds to $\mu = \mu _{stop} = 0.38$, ${\rm \Delta} \mu _{hyst} = 0.16$ and $I_*= 4 \times 10^{-2}$. See movie 6 in the Supplementary material and movies for a visualisation of the collapse in time.

We should stress that the corresponding increase of the ‘static’ friction coefficient (up to $\mu + {\rm \Delta} \mu _{hyst} = 0.38 + 0.16 = 0.54$) provided here is only phenomenological, and does not correspond to an independent experimental measurement. However, it suggests that the actual rest-to-flow yield transition, characterised by $\mu (I = 0)$, does not directly correspond to the experimental avalanche friction coefficient $\mu _a$, which measurement is highly sensitive to the mechanical noise and could thus incorporate nucleation effects (DeGiuli & Wyart 2017; de Geus & Wyart 2022) close to the free-surface, where the pressure conditions can be impacted by preparation effects and small irregularities. Within our approach, the measured $\mu _a$ would thus correspond to some small but non-zero $I$, somewhere between $0$ and $I_*$, more representative of the characteristic collapse inertial numbers than the smallest ones typically used to extract $\mu _{stop}$ from steady inclined flows, which are of order $I \lesssim I_0 = 0.279$ (Pouliquen 1999; Jop et al. 2005). Note that setting a higher $\mu (I=0)$ than $\mu _a$ also appears consistent with the higher friction coefficient values deduced from uniaxial or triaxial compression tests on similar micrometric glass bead materials (Ancey 2001; Adjemian & Evesque 2004; Cui et al. 2017) where friction coefficient values range from $\tan 26.5^\circ \sim 0.5$ to $\tan 30^\circ \sim 0.58$, much larger than $\mu _a=0.44$, illustrating the sensitivity of avalanche onsets to pressure conditions. Despite its simplicity, this ad hoc hysteretic model thus supports the prominence of low inertial number effects in granular collapses, which appear mostly driven by solid–liquid transitions in the transient case, where the slope is not sufficient to sustain steady flow.

7. Conclusion

Our 2-D and 3-D non-smooth numerical model can faithfully simulate granular collapses, and quantitatively predict the final deposits for a wide range of bed inclinations, channel widths and column aspect ratios using a fully plastic model with a plain Drucker–Prager rheology. In contrast to previous numerical investigations advocating more complex rheologies (Lagrée et al. 2011; Dunatunga & Kamrin 2015; Mast et al. 2015), this suggests that transient granular flows are mostly driven by transitions from rest to flow, and that the final stable states can accurately be described by a unique constant bulk friction coefficient.

Comparisons with experimental collapses furthermore support the use of a friction coefficient corresponding to the avalanche angle, as opposed to the stop angle measured from steady experiments. This observation is not restricted to our numerical study, and was for example also noted by Ionescu et al. (2015), albeit not interpreted in this way: while the correction initially attributed to the role of the lateral walls cannot hold in the light of the study provided in § 5, the effective friction coefficient they use turns out to coincide with the avalanche angle measured on their experiment.

Owing to the low inertial numbers involved in the friction-dominated collapse flows, the $\mu (I)$ rheology, which accounts for additional viscous collisional dissipation, is not able to quantitatively improve flow predictions. Furthermore, adjustment of the parameters of the $\mu (I)$ law using steady inclined flow experiments, in particular the stop friction coefficient $\mu _{stop}$, appears unable to accurately predict experimental collapse profiles, and systematically underestimates frictional dissipation within the material.

Studying the effect of friction on the lateral walls on 3-D simulated collapses in narrow flumes with various widths, we have shown that a linear rescaling of the bulk friction coefficient for predicting 2-D collapses is still valid, but with an effective flow thickness much smaller than the maximum flow thickness used in Ionescu et al. (2015).

As a first step to explore low-velocity extensions of the rheology, we have implemented a simple hysteretic law inspired by the work of Pouliquen & Forterre (2002) and DeGiuli & Wyart (2017). This two-valued rheology, reminiscent of the static and dynamic regimes of Coulomb friction, significantly improves the prediction of the collapse dynamics, especially at the onset of the flow, while only weakly affecting the final state and run-out. Again, this points out the crucial role played by static–flowing transitions in such transient configurations, which cannot be reproduced quantitatively by parameters measured in steady flow experiments. The extraction of solid–liquid friction coefficients is, however, very sensitive to mechanical noise and preparation effects, as the decreasing dependence of friction at low inertial numbers is a large source of instability, which can prevent access to the true $I=0$ limit in free-surface geometries. We should also mention that the solid–liquid transition is strongly affected by non-local effects, as mechanical interactions at the grain scale are precisely responsible for nucleation or strengthening effects (Kamrin 2019; Perrin et al. 2021), and improving our non-smooth numerical model to account for a non-local rheology would definitely provide additional insights into the role of spatial inhomogeneities in transient flows (Mowlavi & Kamrin 2021).

Our non-smooth solver provides good predictions of complex 3-D granular flows involving high granular deformations and interactions with frictional surfaces. However, the full potential of Sand6 supporting cohesion and dynamic interaction with complex objects remains to be fully exploited. A next step would be to simulate full 3-D steady flows with well controlled frictional boundary conditions and compare them to experiments in order to accurately bridge the gap between the constant friction rheology and the $\mu (I)$ rheology.