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 (Reference Drucker and Prager1952). 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 (Reference Dunatunga and Kamrin2015), Mast et al. (Reference Mast, Arduino, Mackenzie-Helnwein and Miller2015) or Klár et al. (Reference Klár, Gast, Pradhana, Fu, Schroeder, Jiang and Teran2016), 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 Reference Rauter, Barker and Fellin2020).

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 (Reference Narain, Golas and Lin2010) and Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015). 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. (Reference Lagrée, Staron and Popinet2011), Chauchat & Médale (Reference Chauchat and Médale2014) or Franci & Cremonesi (Reference Franci and Cremonesi2019), or solved using a full elastoplastic model with explicit or implicit return-mapping projections as in Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015), Mast et al. (Reference Mast, Arduino, Mackenzie-Helnwein and Miller2015) or Klár et al. (Reference Klár, Gast, Pradhana, Fu, Schroeder, Jiang and Teran2016).

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. (Reference Ionescu, Mangeney, Bouchut and Roche2015) 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 (Reference Daviet and Bertails-Descoubes2016b), 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 Reference Moreau1994; Jean Reference Jean1999) to solve for the resulting nonlinear equation. Details about this approach can be found in Daviet & Bertails-Descoubes (Reference Daviet and Bertails-Descoubes2016a,Reference Daviet and Bertails-Descoubesb), 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 (Reference Daviet and Bertails-Descoubes2016a,Reference Daviet and Bertails-Descoubesb), 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 (Reference Fukushima, Luo and Tseng2002) and modified in Daviet, Bertails-Descoubes & Boissieux (Reference Daviet, Bertails-Descoubes and Boissieux2011) 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 Reference Savage1979; Daerr & Douady Reference Daerr and Douady1999b; Pouliquen & Forterre Reference Pouliquen and Forterre2002; Artoni, Santomaso & Canu Reference Artoni, Santomaso and Canu2011). 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.2. Collapse set-up

The experimental apparatus is inspired by Balmforth & Kerswell (Reference Balmforth and Kerswell2005) 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. (Reference Farin, Mangeney and Roche2014), 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 Reference Shi1994; Miozzi, Jacob & Olivieri Reference Miozzi, Jacob and Olivieri2008), 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 (Reference Rousseau and Ancey2020) 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 (Reference Balmforth and Kerswell2005) 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 (Reference Balmforth and Kerswell2005) or Russell et al. (Reference Russell, Johnson, Edwards, Viroulet, Rocha and Gray2019), 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. (Reference Farin, Mangeney and Roche2014) and Balmforth & Kerswell (Reference Balmforth and Kerswell2005) 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 Reference Pouliquen1999; Pouliquen & Forterre Reference Pouliquen and Forterre2002), 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 Reference Pouliquen1999), or the initially static bed is inclined up until it starts to flow, and then stops (Pouliquen & Forterre Reference Pouliquen and Forterre2002).

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 (Reference Pouliquen and Forterre2002) 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 (Reference Pouliquen and Forterre2002), 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. (Reference Jop, Forterre and Pouliquen2005), 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 (Reference Pouliquen1999), Pouliquen & Forterre (Reference Pouliquen and Forterre2002), Forterre & Pouliquen (Reference Forterre and Pouliquen2003) and Farin et al. (Reference Farin, Mangeney and Roche2014). We observe that despite some spreading in all the collected data, which might be attributed to the different experimental protocols used (Pouliquen & Forterre (Reference Pouliquen and Forterre2002) versus Pouliquen (Reference Pouliquen1999)) and/or to the slightly different beads used (Pouliquen & Forterre (Reference Pouliquen and Forterre2002) versus Farin et al. (Reference Farin, Mangeney and Roche2014)), we obtain a good agreement with Pouliquen & Forterre (Reference Pouliquen and Forterre2002), 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 (Reference Pouliquen and Forterre2002). 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 (Reference Pouliquen and Forterre2002), Pouliquen (Reference Pouliquen1999), Forterre & Pouliquen (Reference Forterre and Pouliquen2003) and Farin et al. (Reference Farin, Mangeney and Roche2014) 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 Reference Pouliquen1999; Pouliquen & Forterre Reference Pouliquen and Forterre2002; Jop et al. Reference Jop, Forterre and Pouliquen2006). 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 (Reference Pouliquen and Forterre2002), 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}$.

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 (Reference Balmforth and Kerswell2005). 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. (Reference Ionescu, Mangeney, Bouchut and Roche2015).

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. (Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017). 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. (Reference Jop, Forterre and Pouliquen2006) 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. (Reference Jop, Forterre and Pouliquen2005).

This model has recently been used with success for transient collapse prediction in several studies (Lagrée et al. Reference Lagrée, Staron and Popinet2011; Ionescu et al. Reference Ionescu, Mangeney, Bouchut and Roche2015; Martin et al. Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017).

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. (Reference Jop, Forterre and Pouliquen2005), 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. (Reference Jop, Forterre and Pouliquen2005); 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. (Reference Jop, Forterre and Pouliquen2005) 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. (Reference Valette, Riber, Sardo, Castellani, Costes, Vriend and Hachem2019) 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. (Reference Lagrée, Staron and Popinet2011), or Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015) and Martin et al. (Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017). 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. Reference Jop, Forterre and Pouliquen2005; Lajeunesse et al. Reference Lajeunesse, Monnier and Homsy2005; Ionescu et al. Reference Ionescu, Mangeney, Bouchut and Roche2015; Martin et al. Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017). 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 (Reference Balmforth and Kerswell2005) or Zhang et al. (Reference Zhang, Su, Lei and Chen2021)), 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.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. (Reference Lajeunesse, Mangeney-Castelnau and Vilotte2004), Balmforth & Kerswell (Reference Balmforth and Kerswell2005), Lagrée et al. (Reference Lagrée, Staron and Popinet2011) and Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015), 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 (Reference Kerswell2005) and Balmforth & Kerswell (Reference Balmforth and Kerswell2005), 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. (Reference Lagrée, Staron and Popinet2011) and Dunatunga & Kamrin (Reference Dunatunga and Kamrin2015), our results thus not only confirm the robustness of the power-law scalings observed by Lajeunesse et al. (Reference Lajeunesse, Monnier and Homsy2005), Balmforth & Kerswell (Reference Balmforth and Kerswell2005) and Lube et al. (Reference Lube, Huppert, Sparks and Freundt2005), 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 (Reference Savage1979), Taberlet et al. (Reference Taberlet, Richard, Valance, Losert, Pasini, Jenkins and Delannay2003) and Jop et al. (Reference Jop, Forterre and Pouliquen2005) 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. (Reference Ionescu, Mangeney, Bouchut and Roche2015) 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. (Reference Zhang, Su, Lei and Chen2021), 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 (Reference Balmforth and Kerswell2005) or Zhang et al. (Reference Zhang, Su, Lei and Chen2021), 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. (Reference Ionescu, Mangeney, Bouchut and Roche2015).

For the sake of comparison, we also show the linear scaling hypothesis used by Ionescu et al. (Reference Ionescu, Mangeney, Bouchut and Roche2015) 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. (Reference Ionescu, Mangeney, Bouchut and Roche2015) 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.