A critical state μ(I)-rheology model for cohesive granular flows

Abstract

The dynamic behaviour of granular media can be observed widely in nature and in many industrial processes. Yet, the modelling of such media remains challenging as they may act with solid-like and fluid-like properties depending on the rate of the flow and can display a varying apparent friction, cohesion and compressibility. Over the last two decades, the $\mu (I)$-rheology has become well established for modelling granular liquids in a fluid mechanics framework where the apparent friction $\mu$ depends on the inertial number $I$. In the geo-mechanics community, modelling the deformation of granular solids typically relies on concepts from critical state soil mechanics. Along the lines of recent attempts to combine critical state and the $\mu (I)$-rheology, we develop a continuum model based on modified cam-clay in an elastoplastic framework which recovers the $\mu (I)$-rheology under flow. This model permits a treatment of plastic compressibility in systems with or without cohesion, where the cohesion is assumed to be the result of persistent inter-granular attractive forces. Implemented in a two- and three-dimensional material point method, it allows for the trivial treatment of the free surface. The proposed model approximately reproduces analytical solutions of steady-state cohesionless flow and is further compared with previous cohesive and cohesionless experiments. In particular, satisfactory agreements with several experiments of granular collapse are demonstrated, albeit with shear bands which can affect the smoothness of the surface. Finally, the model is able to qualitatively reproduce the multiple steady-state solutions of granular flow recently observed in experiments of flow over obstacles.

1. Introduction

Understanding the mechanical behaviour of granular media is essential in various pharmaceutical and manufacturing processes as well as in geophysical events like landslides, debris flows and snow avalanches (Mangeney et al. 2007; Persson, Alderborn & Frenning 2011; Lucas, Mangeney & Ampuero 2014; Steinkogler et al. 2015; Delannay et al. 2017). As the number of grains in such flows are usually too high to be modelled in a discrete manner, accurate continuum models are needed to efficiently simulate these large-scale natural phenomena and contribute to improving, e.g. hazard assessment and mitigation in mountainous regions. Significant progress has been made, in particular over the last decades, with improved constitutive and rheological models as well as new numerical schemes. Despite this progress, it remains challenging to achieve general models that can describe the diverse behaviour of granular systems. The diversity stems from the varying apparent friction and solid fraction, as well as the effect of possible inter-granular cohesive forces, or even elastic forces in the slow flow regime (Campbell 2002, 2005; Jop, Forterre & Pouliquen 2006; Pouliquen et al. 2006; Pirulli & Mangeney 2008; Rao & Nott 2008; Rognon et al. 2008b; Moretti et al. 2012; Andreotti, Forterre & Pouliquen 2013; Radjai, Roux & Daouadji 2017; Moretti et al. 2020; Vo et al. 2020; Mandal, Nicolas & Pouliquen 2020; Gans et al. 2023; Poulain et al. 2023).

A major contribution to the understanding of dense granular flows came with the collection of experimental data and discrete simulations of granular systems under simple shear, in particular by GDR MiDi (2004). In various configurations, they showed that the shear stress is proportional to the normal stress with a factor of proportionality, a friction $\mu = \mu (I)$, depending on a dimensionless quantity called the inertial number $I$. Small inertial numbers describe quasi-static flows where the macroscopic deformation is slow compared with the microscopic rearrangements, while large inertial numbers refer to rapid flows. The $I$-dependence of the friction was further generalized to tensorial form by Jop et al. (2006) by considering an analogous relation between the shear stress tensor and the isotropic pressure. Now known as the $\mu (I)$-rheology, it has become a well-established continuum theory for describing dense granular flows. Barker et al. (2015) showed that the system of the incompressible $\mu (I)$-Navier–Stokes equations is well-posed unless the inertial number becomes too high or too small, and Barker & Gray (2017) later suggested partial regularization by altering the functional form of the $\mu (I)$-curve. Such ill-posedness was also observed numerically by Barker et al. (2015), Barker & Gray (2017) and Martin et al. (2017).

With regard to numerical implementations, Lagrée, Staron & Popinet (2011) was the first to incorporate the $\mu (I)$-rheology equations in a two-dimensional finite volume Navier–Stokes model, and they applied it to simulate granular collapse as a continuum. Formally, the $\mu (I)$-rheology can be interpreted as a viscoplastic constitutive law involving a Drucker–Prager yield criterion being dependent on pressure and strain rate (Ionescu et al. 2015; Martin et al. 2017). In such a viscoplastic model, Ionescu et al. (2015) simulated granular collapse using an arbitrary Lagrangian–Eulerian finite element formulation of the Navier–Stokes equations. Relying on a two-dimensional smoothed particle hydrodynamics (SPH) model and a three-dimensional particle finite element method (PFEM), Chambon et al. (2011) and Franci & Cremonesi (2019), respectively, demonstrated the implementation of the $\mu (I)$-rheology in such particle-based schemes, thus enabling the trivial handling of the free surface which can pose difficulties in mesh-based schemes. As a purely mesh-free scheme, SPH requires a computationally expensive neighbour search for each time step, and handling boundary conditions can be challenging (Vacondio et al. 2021). Similarly, PFEM can also become expensive due to the need of creating a new mesh every time the finite element mesh becomes too distorted, although efficient remeshing algorithms may improve this (Cremonesi et al. 2020). Moreover, classical fluid solvers have issues in capturing material disconnections and with representing fully static zones, which become important features in the fluid-to-solid transition. It has been shown that regularization methods with a sufficiently small regularization parameter can reproduce the static–flowing transition (Frigaard & Nouar 2005; Lusso et al. 2017) as well as augmented Lagrangian methods (Martin et al. 2017). In a different approach, Dunatunga & Kamrin (2015) considered an $I$-dependent elastoplastic model with a Drucker–Prager yield in a two-dimensional material point method (MPM).

The ill-posedness of the $\mu (I)$-rheology for very small and large values of inertial numbers is an indication that important physics are missing in the model. In particular, for small inertial numbers, i.e. slow flows, concepts from solid mechanics like elasticity and shear banding become important (Barker et al. 2015). While the $\mu (I)$-rheology has been successful in modelling granular liquids, critical state soil mechanics (CSSM) (Schofield & Wroth 1968) is well established in the geo-mechanics community for describing the, potentially large, deformation of granular media as elastoplastic solids. Initially conceived to model soils and clay, this theory is based on the existence of a critical state, i.e. as a frictional fluid state where shear deformations can continue without any further change of volume, pressure or shear stress. All such critical states must reside along a so-called critical state line (CSL) in the space given by the pressure and shear stress. Among CSSM models, the modified cam-clay (MCC) model first formulated by Roscoe & Burland (1968) is arguably one of the most used in numerical implementations, giving an ellipsoidal yield surface in principal stress space. The CSSM models have been successfully applied to modelling the solid to fluid transition in snow avalanches (Gaume et al. 2018, 2019; Trottet et al. 2022). They have also been used to simulate the avalanche dynamics over complex topography (Li et al. 2021, 2024; Cicoira et al. 2022). However, these CSSM models are based on a rate-independent theory, which can therefore not capture the $I$-dependent behaviour of granular systems in their liquid state.

The success of CSSM in capturing snow avalanches is, arguably, due to its ability to consider cohesion as well as plastic compaction and dilatancy. Cohesion refers to the presence of an attractive force between the grains, with the attractive force being either persistent or not. In fact, cohesion is key in understanding snow avalanches. The extensive snow experiments of Rognon et al. (2008a) showed that flowing dry snow behaves as dense granular flows with varying grain sizes, suggesting that this variation is a result of aggregates forming and size segregation within the flow, thus implying an important role of the cohesive forces in the snow. The cohesion has varying physical origins. In snow, it is a result of sintering of ice crystals or due to liquid water creating bonds by capillarity between snow particles (Szabo & Schneebeli 2007; Steinkogler et al. 2015; Ligneau, Sovilla & Gaume 2022). In fine powders, it may also be of electrostatic or magnetic origin (Castellanos 2005).

In an effort to describe the rheology of cohesive granular flows, Rognon et al. (2006, 2008b) suggested on the basis of two-dimensional discrete simulations to make the coefficients in the linear $\mu (I)$-relation proposed by da Cruz et al. (2005) dependent on a dimensionless number giving the ratio between the maximum inter-granular attractive force to the average normal force. Similar conclusions were made by Khamseh, Roux & Chevoir (2015) and Badetti et al. (2018) based on simulations and experiments of wet granular materials, showing that the apparent friction should depend both on the inertial number and the dimensionless cohesion number. Based on discrete simulations in both two dimensions and three dimensions, Berger et al. (2015) and Vo et al. (2020), respectively, proposed instead a new generalized inertial number dependent on the same dimensionless cohesion number.

Pertaining to the plastic compressibility of granular systems, in CSSM this is realized in a simple manner through a hardening law of the yield criterion. Specifically, in CSSM it is typically assumed that the inverse solid volume fraction decreases with the logarithm of pressure. While compressibility is not considered in the $\mu (I)$-rheology, a relation between the solid volume fraction $\phi$ and the inertial number was shown by GDR MiDi (2004), da Cruz et al. (2005) and Pouliquen et al. (2006). This prompted the compressible $\mu (I), \phi (I)$-rheology, with a monotonically decreasing, typically linear, function $\phi (I)$ for the solid fraction. Unfortunately, this $\mu (I), \phi (I)$-rheology is mathematically ill-posed in time-dependent problems (Heyman et al. 2017; Schaeffer et al. 2019).

The critical behaviour of granular flow was recently demonstrated in an extensive set of three-dimensional discrete simulations of granular systems by Wang, Li & Liu (2022), showing that, under shear loading, a critical state is reached at sufficient shear for all shear rates, and that the exact critical state depends on the rate in a way that can be described by the inertial number. Previous depth-averaged models have introduced the critical state concept into the $\mu (I)$-rheology, with the first attempt made by Pailha & Pouliquen (2009) focusing on granular materials immersed in a secondary fluid. This study was augmented by Bouchut et al. (2016) and Garres-Díaz et al. (2020) in a multi-layer model accounting for in-depth heterogeneities and later presented for the fully dry case by Bouchut et al. (2021) in order to provide a compressible $\mu (I)$-rheology model, all studies successfully compared with laboratory experiments. While all of these previous models relied on depth-averaged equations in one or more layers based on shallow water assumptions, they also relied on the critical state theory of Roux & Radjai (1998). In this viscoplastic approach, the dynamics is governed by a frictional parameter which can evolve depending on the deviation of the solid fraction from a prescribed critical solid fraction, thus different from the elastoplastic theory of CSSM. In an alternative approach, Rauter (2021) supplements the otherwise linearly decreasing $\phi (I)$ in a Navier–Stokes-type model with the critical state model of Johnson & Jackson (1987) and Johnson, Nott & Jackson (1990), giving an expression for the pressure as a function of the shear rate and solid fraction.

The idea of combining CSSM and the $\mu (I)$-rheology was conceived by Barker et al. (2017) and Schaeffer et al. (2019), albeit in a purely viscoplastic form, resulting in equations they termed compressible $I$-dependent rheology. While not based on the elastoplastic framework of CSSM, compressibility was considered through the evolution of a (closed) yield surface and a carefully chosen plastic flow rule. Their motivation came from the specific desire of making the rheology well-posed for all inertial numbers through the introduction of compressibility. While succeeding in proving the well-posedness of their model subject only to reasonable restrictions, this conceptual and analytical study did not offer a complete numerical realization beyond simple one-dimensional demonstrations, nor an experimental comparison or validation.

In an effort to capture the behaviour of granular media through both their solid and fluid states in a continuum model, it seems natural to seek a model that combines CSSM, which readily incorporates cohesion and compressibility, with the $\mu (I)$-rheology, allowing the inertial number dependence in flows. Now, guided by advances in particle-based numerical schemes and concepts from computational elastoplasticity, we are able to formulate and implement an elasto-viscoplastic MCC model that follows the $\mu (I)$-rheology under flow conditions. This allows the modelling of cohesive and compressible granular materials as elastoplastic solids which, in the fluid transition, recovers the expected $I$-dependent behaviour. In short, this is accomplished by making the CSL dependent on the inertial number and by relying on MPM which, due to its hybrid Eulerian–Lagrangian nature, has the ability to handle large deformations as the material transitions from solid-like to fluid-like, including material separation and reconsolidation. While MPM also offers the advantage of not requiring any special numerical treatment of the free surface, it comes with increasing computational costs (compared with classical finite element methods) as well as the need for a careful treatment of numerical dissipation. The type of cohesion treated in the model is one that does not cease to exist, e.g. in order to consider the situation of permanent attractive inter-granular force. However, as will be discussed, other types of cohesion could be considered in the proposed general framework that would make the model more suitable for landslides and debris flows. As an elastoplastic model implemented in MPM, the proposed model is inspired by the work of Dunatunga & Kamrin (2015), although differing both in the elastic and plastic equations. In addition to considering cohesion and compressibility through the critical state, the proposed model relies on advancements from the MPM community which improve stability and reduce numerical dissipation. In this paper, the model is demonstrated for collapse and flow problems where the magnitude of the elastic deformations is negligible compared with the (rate-dependent) plastic deformations, and the general solid behaviour and elastic nature of granular media will not be treated in detail here.

2. Theory

The granular material is here considered as a continuum governed by the conservation laws for mass and momentum. Conservation of mass is given by

(2.1)\begin{equation} \frac{{\rm D} \rho}{{\rm D} t} + \rho(\boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{\cdot} \boldsymbol{v}) = 0, \end{equation}

and conservation of momentum under gravity states that

(2.2)\begin{equation} \rho \frac{{\rm D} \boldsymbol{v}}{{\rm D} t} - \boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{\cdot} \boldsymbol{\sigma}_0 - \rho \boldsymbol{g} = 0 , \end{equation}

where $\rho = \rho (\boldsymbol {x},t)$ is the mass density, $\boldsymbol {v} = \boldsymbol {v}(\boldsymbol {x},t)$ is the velocity, $\boldsymbol {\sigma }_0 = \boldsymbol {\sigma }_0(\boldsymbol {x},t)$ is the symmetric Cauchy stress tensor, $\boldsymbol {g}$ is the constant gravitational acceleration and ${\rm D}/{\rm D}t$ denotes the material time derivative. Without further assumptions, the system of equations given by these conservation laws is not closed. In this section, constitutive laws relating the deformation to the stress are outlined in order to close the system. In particular, the granular material will be treated as an elastoplastic continuum where the plastic, i.e. irreversible, deformations are highly rate dependent.

2.1. Definitions of deformation

With the deformation of a body over a time $t$ described by a one-to-one mapping from a material (initial) coordinate $\boldsymbol {X}$ to a spatial (current) coordinate $\boldsymbol {x}$, the deformation gradient $\boldsymbol {F}$ is defined with $F_{ij} = \partial x_i / \partial X_j$ and $\det (\boldsymbol {F}) > 0$. The singular values of this second-order invertible tensor are the principal stretches $\lambda _\alpha > 0$, $\alpha =1,\ldots,\dim$, where $\dim =1,2,3$ denotes the number of spatial dimensions. From the deformation gradient, the Hencky strain tensor is defined as

(2.3)\begin{equation} \boldsymbol{\varepsilon}= \frac{1}{2} \ln \boldsymbol{F} \boldsymbol{F}^T = \sum_{\alpha=1}^{\dim} \ln \lambda_\alpha\ \boldsymbol{n}_\alpha \otimes \boldsymbol{n}_\alpha, \end{equation}

where the unit vectors $\boldsymbol {n}_\alpha$ define the spatial principal directions. Moreover, the velocity gradient tensor can be written as

(2.4)\begin{equation} \boldsymbol{L} = \boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{v} = \dot{\boldsymbol{F}} \boldsymbol{F}^{{-}1}, \end{equation}

where $\dot {\boldsymbol {F}}$ denotes the material time derivative of the deformation gradient. The trace of the Hencky strain, $\varepsilon _V = {\rm tr}\,\boldsymbol {\varepsilon } = \ln ( \det (\boldsymbol {F}) )$, is a measure of the volumetric deformation, with material time derivative $\dot {{\varepsilon }}_V = {\rm tr}(\boldsymbol {L})$.

2.2. Elasticity and yield surface

To account for both elastic, reversible, and plastic, irreversible, deformations, it is common to consider the multiplicative decomposition of the deformation gradient

(2.5)\begin{equation} \boldsymbol{F} = \boldsymbol{F}^E \boldsymbol{F}^P, \end{equation}

where $\boldsymbol {F}^E$ refers to the deformation arising from the elastic forces, while $\boldsymbol {F}^P$ is the permanent, plastic, component. We seek a law that relates the stresses (due to elastic forces) to the elastic component of the Hencky strain

(2.6)\begin{equation} \boldsymbol{\varepsilon}^E = \tfrac{1}{2} \ln ( \boldsymbol{F}^E (\boldsymbol{F}^E)^T ) {,} \end{equation}

through a so-called strain energy density function $\varPsi$, i.e. a hyperelastic law. In this work, a frame-indifferent isotropic hyperelastic model known as the Hencky model is adopted, where

(2.7)\begin{equation} \varPsi(\boldsymbol{\varepsilon}^E) = \tfrac{1}{2} \lambda ({\rm tr}\,\boldsymbol{\varepsilon}^E)^2 + G \,{\rm tr}((\boldsymbol{\varepsilon}^E)^2), \end{equation}

which results in the Kirchhoff stress $\boldsymbol {\sigma } = \det (\boldsymbol {F}) \boldsymbol {\sigma }_0$ being related to the Hencky strain through

(2.8)\begin{equation} \boldsymbol{\sigma} = \frac{\partial \varPsi(\boldsymbol{\varepsilon}^E)}{\partial \boldsymbol{\varepsilon}^E} = \lambda ({\rm tr}\,\boldsymbol{\varepsilon}^E) \boldsymbol{I} + 2 G \boldsymbol{\varepsilon}^E , \end{equation}

where $\boldsymbol {I}$ denotes the second-order identity tensor, and $\lambda$ and $G$ are the two Lamé parameters which may be related to Young's modulus $E$ and Poisson's ratio $\nu$. Equation (2.8) may also be written

(2.9)\begin{equation} \boldsymbol{\sigma} = \mathcal{C} : \boldsymbol{\varepsilon}^E, \end{equation}

with the fourth-order isotropic stiffness matrix $\mathcal {C}$ well known from linear elasticity with elements $\mathcal {C}_{ijkl} = \lambda \delta _{ij}\delta _{kl} + G(\delta _{ik}\delta _{jl} + \delta _{il}\delta _{jk})$. The properties of this model are elegantly demonstrated in Xiao & Chen (2002), showing that, with this model, the Kirchhoff stress and the Hencky strain constitute a work–conjugate pair.

It is instructive to relate Hencky's elasticity model to hypoelastic models where a widely used relation between an objective stress rate $\mathring {\boldsymbol {\sigma }}$ and the elastic rate-of-deformation $\boldsymbol {D}^E = \frac {1}{2}(\boldsymbol {L}^E+(\boldsymbol {L}^E)^T)$ can be written as

(2.10)\begin{equation} \mathring{\boldsymbol{\sigma}} = \mathcal{C} : \boldsymbol{D}^E {.} \end{equation}

Choosing the Jaumann stress rate, this is similar to the rate equation used in the elastoplastic $\mu (I)$-rheology model of Dunatunga & Kamrin (2015) except they considered the Cauchy stress. Equation (2.10) should be exactly integrable to give an elastic relation. If this so-called self-consistency criterion is to be fulfilled, then it can be shown that the logarithmic rate of the Kirchhoff stress is the only choice among all corotational rates (Xiao, Bruhns & Meyers 1999). To go even further, the logarithmic rate is argued by Xiao, Bruhns & Meyers (2000) as the only choice among corotational and non-corotational rates in the context of elastoplasticity. Interestingly, the integration of (2.10) with the logarithmic rate assuming a stress-free initial state yields exactly Hencky's elasticity model (Bruhns, Xiao & Meyers 1999; Xiao & Chen 2002).

The onset of plastic deformations is characterized by a yield function $y = y(\boldsymbol {\sigma })$. Conventionally, $y \leqslant 0$ defines admissible states, and the yield surface $y=0$ defines the limit where deformations are no longer elastic. Thus, plastic deformations only occur on the yield surface. Here, the yield function will be expressed in terms of the stress invariants

(2.11a,b)\begin{equation} p ={-}\frac{1}{\dim} \,{\rm tr}\,{\boldsymbol{\sigma}} \quad\textrm{and}\quad q = \frac{1}{\sqrt{2} } || \boldsymbol{\tau} || , \end{equation}

where $p$ is termed the isotropic pressure and $q$ the equivalent shear stress. The latter was defined from the deviatoric stress tensor

(2.12)\begin{equation} \boldsymbol{\tau} = \operatorname{dev}(\boldsymbol{\sigma}) = \boldsymbol{\sigma} - \frac{1}{\dim}\,{\rm tr}(\boldsymbol{\sigma})\boldsymbol{I}, \end{equation}

using the notation $|| \boldsymbol {\tau } || = \sqrt {\boldsymbol {\tau } : \boldsymbol {\tau }}$.

While yielding of many granular media, including sand and systems of glass beads, follows a Drucker–Prager criterion, other granular media, e.g. snow, are best described by a yield function providing a closed yield surface in stress space (Reiweger, Gaume & Schweizer 2015; Ritter, Gaume & Löwe 2020; Blatny et al. 2021). To this end, we consider the yield criterion given by the (cohesive) MCC model

(2.13)\begin{equation} y(\boldsymbol{\sigma}) = y(p,q) = q^2 - \mu^2 (p+\beta p_c) (p_c-p) , \end{equation}

where $p_c \geqslant 0$ represents an isotropic compressive strength, $\beta \geqslant 0$ is a dimensionless measure of cohesion and $\mu >0$ defines the slope of the CSL along which all critical states reside. Alternatively, one may say that all stress states along the CSL are on a Drucker–Prager yield surface $q=\mu (p+\beta p_c)$. The MCC yield surface is illustrated in figure 1.

Figure 1. The MCC yield surface $y(p,q)=0$ from (2.13) in the space of the stress invariants $p$ and $q$ defined in (2.11a,b). Elastic states satisfy $y<0$. The dashed line represents the CSL. The dash-dotted line marks the value of $p$ at the apex of the ellipse, separating states undergoing plastic compaction and dilation as a result of the chosen associative plastic flow rule.

In this model, $-\beta p_c$ represents the isotropic tensile strength. The case $\beta = 0$ refers to the cohesionless case, and in this case the material cannot sustain tensile stress states. For $\beta > 0$, the system will always have a non-zero tensile strength (unless $p_c = 0$, i.e. the stress-free state). As such, cohesion will always be present, which may not be the case in, e.g. landslides. For granular systems in which the cohesion is better described by breakable bonds, $\beta$ may easily be set to zero after yield, however, such a situation will not be considered here.

In order to complete the description of the model, a plastic flow rule will be outlined in § 2.3, consequently determining the onset and magnitude of plastic dilation and compaction. In § 2.4, a hardening law is introduced that governs the evolution of $p_c$ and thus controls the, possibly extensive, variations of solid fraction. Finally, § 2.5 introduces dependence on the inertial number.

2.3. Plastic flow rule

Under the assumption that the rate of work per unit undeformed volume can be additively decomposed into an elastic and plastic part, it follows that

(2.14)\begin{equation} \boldsymbol{\sigma} : \boldsymbol{L} = \boldsymbol{\sigma} : \boldsymbol{L}^E + \boldsymbol{\sigma} : \boldsymbol{L}^P, \end{equation}

where an elastic and a plastic part of the velocity gradient were defined. Following conventions, the plastic part $\boldsymbol {L}^P$ is chosen as the symmetric tensor

(2.15)\begin{equation} \boldsymbol{L}^P = \dot{\gamma} \frac{\partial g / \partial \boldsymbol{\sigma} }{||\partial g / \partial \boldsymbol{\sigma}||} ,\end{equation}

where the equivalent plastic strain $\dot {\gamma } \geqslant 0$ is defined as $\dot {\gamma } = || \boldsymbol {L}^P ||$ and $g$ is the so-called plastic potential. In MCC models, the typical choice is an associative plastic flow rule, i.e. $g=y$ (the plastic potential is associated with the yield function). In other words, the plastic rate of deformation coincides with the gradient of the yield surface. An associative flow rule is a result of the principle of maximum plastic dissipation, in particular it can be shown that maximizing $\boldsymbol {\sigma } : \boldsymbol {L}^P$ subject to $y \leqslant 0$ leads to (2.15) where $\dot {\gamma }$ can be interpreted as a Lagrange multiplier (Simo 1992; Bonet & Wood 2008).

With $g = y(\boldsymbol {\sigma }) = y(p,q)$, (2.15) can be written

(2.16)\begin{equation} \boldsymbol{L}^P = \frac{1}{\sqrt{2}} \dot{\gamma} \frac{\partial y / \partial q }{||\partial y / \partial \boldsymbol{\sigma}||} \frac{ \boldsymbol{\tau}}{|| \boldsymbol{\tau} ||} - \frac{1}{\dim} \dot{\gamma} \frac{\partial y / \partial p}{||\partial y / \partial \boldsymbol{\sigma}||} \boldsymbol{I} {,} \end{equation}

promoting the definition of the equivalent plastic shear strain rate as

(2.17)\begin{equation} \dot{\gamma}_{_S} = \dot{\gamma} \frac{\partial y / \partial q }{||\partial y / \partial \boldsymbol{\sigma}||} = \sqrt{2} || \operatorname{dev}(\boldsymbol{L}^P) || , \end{equation}

and the equivalent plastic volumetric strain rate as

(2.18)\begin{equation} \dot{\gamma}_{_V} ={-} \dot{\gamma} \frac{\partial y / \partial p}{||\partial y / \partial \boldsymbol{\sigma}||} = {\rm tr}(\boldsymbol{L}^P) , \end{equation}

such that the plastic volumetric Hencky strain is $\varepsilon _V^P = \ln ( \det (\boldsymbol {F}^P) ) = \int \dot{\gamma}_{_V} \,{\rm d}t$. The associative plastic flow rule implies a plastic volume increase for stress states with $p < p_{cs}$ and a volume decrease for $p > p_{cs}$, where $p_{cs} = \frac {1}{2}p_c(1-\beta )$ denotes the pressure at the apex of the yield surface (marked by the dash-dotted vertical line in figure 1). Note that the CSL always goes through the apex, so there should indeed be no volume change for stress states at critical state.

At this point, it can be instructive to relate the current model to the critical state model of Roux & Radjai (1998) which has been used in the depth-averaged (single- or multi-layer) models of Pailha & Pouliquen (2009), Bouchut et al. (2016, 2021) and Garres-Díaz et al. (2020). Assuming for simplicity the cohesionless case, the stress states in that model are forced to be on the (Drucker–Prager) surface $q=\tan (\delta + \psi ) p$, where $\tan (\delta )=\mu$ is the critical friction coefficient and $\psi$ is a so-called dilatancy angle which satisfies $\psi > 0$ under dilation and $\psi < 0$ under compaction. If one imagines $q=\tan (\delta ) p$ intersecting an MCC-ellipsoid at its apex where $p = p_{cs}$, it follows that $q=\tan (\delta + \psi ) p$ must intersect the same ellipsoid at $p < p_{cs}$ under dilation and at $p > p_{cs}$ under compaction.

The choice of plastic flow rule (or plastic potential) is fundamental to the possibility and extent to which the material can undergo volumetric deformations, and accordingly to how critical state is reached. In Drucker–Prager/Mohr–Coulomb models, an associative flow rule would imply every plastic deformation being dilating. The use of non-associative flow rules in such models prevents this. In addition, the singularities in the Drucker–Prager/Mohr–Coulomb yield surfaces require special treatment which often may induce unwanted dilation, although they can be corrected with various strategies (Dunatunga & Kamrin 2015; Pradhana 2017).

2.4. Hardening law

In soil mechanics, the isotropic compressive strength $p_c$ is typically increased with plastic compaction (i.e. hardening – yield surface expands) and decreased with plastic dilation (i.e. softening – yield surface shrinks). As such, $p_c = p_c(\varepsilon _V^P)$ which entails that $p_c$ will not change at critical state. Here, a similar relation will be assumed, specifically

(2.19)\begin{equation} p_c(\varepsilon_V^P) = p_c^0 {e}^{-\xi \varepsilon_V^P}, \end{equation}

where $p_c^0 \geqslant 0$ is an initial compressive strength before any plastic compaction/dilation has occurred, and $\xi \geqslant 0$ is a dimensionless parameter.

This hardening law can be motivated from the observed behaviour of soils (Wood 1991), cohesive and cohesionless powders (Walker 1923; Poquillon et al. 2002; Castellanos 2005) and porous brittle media such as snow (Blatny, Löwe & Gaume 2023), where the inverse solid volume fraction behaves as

(2.20)\begin{equation} \frac{1}{\phi} = \frac{1}{\phi_0} - \varLambda \ln(p_c/p_c^0), \end{equation}

where $\varLambda >0$ is called the plasticity index and $\phi _0$ is the initial solid fraction. In (2.20), the compressive strength is $p_c = p_c^0$ when the solid fraction is $\phi =\phi _0$, with $p_c^0$ and $\phi _0$ being two independent parameters. In reality, $p_c^0$ and $\phi _0$ may be related, with typically an increasing $p_c^0$ for an increasing $\phi _0$. However, in this study, systems with varying $\phi _0$ will not be investigated.

Equation (2.20) can be rearranged as

(2.21)\begin{equation} \frac{\phi_0}{\phi} - 1 ={-} \phi_0 \varLambda \ln(p_c/p_c^0). \end{equation}

Here, the left-hand side may be approximated by the volumetric plastic Hencky strain

(2.22)\begin{equation} \varepsilon_V^P = \det(\boldsymbol{F}^P) \stackrel{(i)}{\approx} \det(\boldsymbol{F}) = \ln(\phi_0/\phi) \stackrel{(ii)}{\approx} \frac{\phi_0}{\phi} - 1, \end{equation}

under the assumption that (i) the elastic volumetric deformation is negligible compared with the plastic volumetric deformation and that (ii) the volumetric deformation is small. This allows (2.21) to be written as

(2.23)\begin{equation} p_c = p_c^0 {e}^{-\varepsilon_V^P / (\varLambda \phi_0) }, \end{equation}

which is exactly (2.19) with $\xi = (\phi _0 \varLambda )^{-1}$. Note that the volumetric plastic Hencky strain $\varepsilon _V^P$ is zero before any deformation occurred (at which point $\phi =\phi _0$ and $p_c = p_c^0$).

2.5. Inertial number dependence

In the $\mu (I)$-rheology, the shear stress is made proportional to the isotropic pressure with a factor of proportionality $\mu$ depending on a dimensionless number known as the inertial number $I$. Assuming the deformation is dominantly plastic, the inertial number can be written

(2.24)\begin{equation} I=\frac{\dot{\gamma}_{_S} d}{\sqrt{p /\rho_*}}, \end{equation}

where $d$ denotes the diameter of the grains and $\rho _*$ is the constant intrinsic grain density.

In the MCC model outlined in §§ 2.2–2.4, it is only at critical state that $\mu$ represents the factor of proportionality between the equivalent shear stress and the pressure. As most stress states in an unobstructed smooth flow would gather close to the CSL, the $\mu (I)$-rheology would be approximately recovered in such flows by letting $\mu$ depend on $I$ in the same way as in the $\mu (I)$-rheology. A typical expression for $\mu (I)$ is such that it is bound by a lower value $\mu _1$ and approaching an upper value $\mu _2$ asymptotically with increasing inertial number, in particular

(2.25)\begin{equation} \mu(I) = \mu_1 + \frac{\mu_2-\mu_1}{\dfrac{I_0}{I}+1}, \end{equation}

from Jop, Forterre & Pouliquen (2005), where $I_0 > 0$ is a dimensionless constant.

Evidently, (2.24) is not well defined for negative pressures, which is possible in the cohesive case, i.e. when $\beta >0$. As a remedy, a new pressure $\bar {p} = p + \beta p_c$ is introduced, where the cohesive stress is added to the pressure, leading to the definition of a cohesion-dependent inertial number

(2.26)\begin{equation} I=\frac{\dot{\gamma}_{_S} d}{\sqrt{\bar{p} /\rho_*}} {.} \end{equation}

With $\beta p_c>0$ here representing the magnitude of the isotropic cohesive stress (tensile strength), this cohesive inertial number was proposed by Berger et al. (2015) and Vo et al. (2020). They demonstrated, through discrete simulations, that this inertial number improved the scaling behaviour for granular flows where cohesion was added to the system through an (unbreakable) attractive inter-granular force. This is in contrast to other studies (Rognon et al. 2006, 2008b; Khamseh et al. 2015; Badetti et al. 2018) which have proposed making the apparent friction dependent on a separate dimensionless cohesion number in addition to the (non-cohesive) inertial number. Note that the non-cohesive inertial number, (2.24), is related to its cohesive counterpart, (2.26), through a multiplication of $\sqrt {1+\beta p_c/p}$. With this inertial number, (2.25) becomes

(2.27)\begin{equation} \mu(I) = \mu_1 + \frac{\mu_2-\mu_1}{\omega\dfrac{\sqrt{\bar{p}}}{\dot{\gamma}_{_S}}+1}, \end{equation}

where, for convenience, the constant $\omega = {I_0}/({d\sqrt {\rho _*}})$ was introduced. The evolution of $\mu$ with $I$ when changing $\omega$ is sketched in figure 2.

Figure 2. Evolution of $\mu (I)$ from (2.27) as the parameter $\omega$ is increased.

In a nutshell, the proposed model consists of the elastic law given by (2.8), the yield function given by (2.13) with $\mu$ being the function defined in (2.27), the plastic flow rule given by (2.15) and the hardening law given by (2.19). A summary of the model parameters is presented in table 1, where it should be noted that $I_0$, $d$ and $\rho _*$ appear in the proposed model through $\omega$ only.

Table 1. Model parameters.

Formally, this elasto-viscoplastic model can be written as a Perzyna overstress model (Perzyna 1963, 1966), in particular

(2.28)\begin{equation} {\dot{\gamma}_{_S}} = \begin{cases} \dfrac{1}{\eta} \left( q-q_y(p) \right), & \textrm{if } q > q_y(p) \\ 0, & \textrm{otherwise}, \end{cases} \end{equation}

with a $p$-dependent yield shear stress

(2.29)\begin{equation} q_y(p) = \mu_1 \sqrt{ (p_c-p) \bar{p}} \geqslant 0, \end{equation}

and an $I$-dependent viscosity

(2.30)\begin{equation} \eta = \frac{ (\mu(I)-\mu_1) \sqrt{ (p_c-p) \bar{p}} }{ \dot{\gamma}_{_S} } = \frac{(\mu_2-\mu_1) \sqrt{ (p_c-p) \bar{p}} }{\dot{\gamma}_{_S} + \omega \sqrt{\bar{p}}} \geqslant 0 {.} \end{equation}

The Perzyna yield surface given by $q_y(p)$ differs from the previous concept of a yield surface as stress states outside the Perzyna yield surface are admissible.

2.6. An operator splitting scheme

Solving for the stress (and deformation) can be accomplished through an elastic predictor–plastic corrector scheme (Simo 1992; de Souza Neto, Perić & Owen 2008). While the details of this operator splitting scheme are presented in Appendix A, the concept is simple: assume first the deformation is purely elastic, and if the yield condition is violated, correct for any plastic deformation through the plastic flow rule. In particular, a trial stress state $(p^{trial}, q^{trial})$ obtained elastically from the previous stress state $(p^n,q^n)$ at time $t^n$ gives the new stress state $(p^{n+1}, q^{n+1})$ after a plastic correction

(2.31)$$\begin{gather} p^{n+1} = p^{trial} + K \dot{\gamma}_{_V} \Delta t \end{gather}$$

(2.32)$$\begin{gather}q^{n+1} = q^{trial} - G \dot{\gamma}_{_S} \Delta t, \end{gather}$$

where $\Delta t = t^{n+1} - t^n$. This operation is known as a return mapping.

Solving the return mapping must be done iteratively, the convergence of which is difficult to guarantee. As an approximation, since the hardening law $p_c = p_c(\varepsilon _V^P)$ is only dependent on the update of $p$, the return mapping can be performed in the view of Perzyna, assuming $\mu = \mu _1$ is constant (cf. (2.29)), producing the desired $p^{n+1}$. With $p^{n+1}$ now known, $q^{n+1}$ can be found by inserting (2.32) into (2.13) and using (2.27), yielding a quadratic equation for $\dot{\gamma}_{_S}$

(2.33)\begin{equation} a \dot{\gamma}^2_{_S} + b \dot{\gamma}_{_S} + c = 0, \end{equation}

where

(2.34)\begin{equation} \left.\begin{gathered} a = G \Delta t \\ b = (\mu_2-\mu_1) \sqrt{ (p_c-p^{n+1}) \bar{p}^{n+1} } + G \Delta t \omega \sqrt{\bar{p}^{n+1}} - (q^{trial} - q_y(p^{n+1})) \\ c ={-}\omega\sqrt{\bar{p}^{n+1}} (q^{trial} - q_y(p^{n+1})) \end{gathered}\right\}. \end{equation}

Since $a>0$ and $c<0$, exactly one positive solution for $\dot{\gamma}_{_S}$ is guaranteed, and the new stress state $(p^{n+1}, q^{n+1})$ is found. Notably, there are no issues with zero pressure or strain rate. The return mapping is illustrated in figure 3.

Figure 3. The return mapping in the view of Perzyna's overstress approach. A trial state $(p^{trial}, q^{trial})$ is projected to the stress state $(p^{n+1},q^{n+1})$ through an intermediate step $(p^{n+1}, q_y(p^{n+1}))$. Note that, since here $p^{trial} > p_{cs}$, volumetric compaction occurred and the yield surface therefore expanded (hardened) from the initial surface coloured grey to the final surface coloured black.

3. Numerical solver

Here, MPM is adopted as a numerical scheme for approximating solutions to (the weak form of) (2.2) subject to the elasto-viscoplastic constitutive model. In short, MPM consists of tracking moving particles (or ‘material points’) which store associated material properties such as mass, velocity and deformation gradient, while also adopting a background Eulerian grid facilitating the momentum computations. This is sketched in figure 4. Although it is commonly stated that the original formulation of MPM dates back to Sulsky, Chen & Schreyer (1994), MPM can be considered an extension of particle-in-cell (PIC) methods (Harlow 1964), notably the fluid-implicit particle (FLIP) method (Brackbill, Kothe & Ruppel 1988). It was first applied to the flow of granular media by Wieckowski, Youn & Yeon (1999) in the context of silo discharge, assuming simply a cohesionless constant Drucker–Prager yield criterion.

Figure 4. Illustrating sketch of the MPM discretization (grid and particles) of a flow on a surface subject to gravity. In this study, the domain is initialized with four to eight MPM particles per grid cell.

The material domain is discretized into $N_p$ particles with constant and equal mass $m_p$, initially equal volume $V_p$ and initially zero velocity $\boldsymbol {v}_p$. With a fixed Eulerian grid with uniform grid spacing $\Delta x$, we consider $C^1$-continuous interpolation functions $N_{i}(\boldsymbol {x}_p) \geqslant 0$ between a grid node $i$ and a particle at a coordinate $\boldsymbol {x}_p$, in particular

(3.1)\begin{equation} N_i(\boldsymbol{x}_p) = \sum_{\alpha=0}^{\dim} N \left( \frac{x_{p \alpha}-x_{i \alpha}}{\Delta x} \right) , \end{equation}

where $\alpha$ denotes the component and with $N( {\cdot } )$ given by the quadratic B-spline

(3.2)\begin{equation} N(x) =\begin{cases} \frac{3}{4} - |x|^2, & \textrm{if } |x| < \frac{1}{2} \\[6pt] \frac{1}{2}(\frac{3}{2} - |x|)^2, & \textrm{if } \frac{1}{2} \leqslant |x| < \frac{3}{2} \\[6pt] 0, & \textrm{otherwise.} \end{cases} \end{equation}

Such B-spline MPM circumvents the issue of cell crossing instabilities (Steffen, Kirby & Berzins 2008) associated with the original MPM which relied on linear hat functions (e.g. used in the elastoplastic $\mu (I)$-rheology model of Dunatunga & Kamrin 2015). Not considering the boundary term, the temporal and spatial discretization of the weak form of the momentum conservation equation can be written on the Eulerian grid nodes $i$ as

(3.3)\begin{equation} \frac{ m^n_{i} \boldsymbol{v}^{n+1}_{i} - m^n_{i} \boldsymbol{v}^{n}_{i} }{\Delta t} ={-} \sum_p V^0_{p} \boldsymbol{\sigma}(\boldsymbol{\varepsilon}^{E,n}) \boldsymbol{\nabla} N_{i}(\boldsymbol{x}^n_p) + m_i^n \boldsymbol{g} ,\end{equation}

where $V^0_p$ is the initial particle volume and the nodal mass

(3.4)\begin{equation} m_i^n = \sum_p m_p N_i(\boldsymbol{x}_p^n)\end{equation}

is obtained under the common assumption of mass lumping. The MPM discretization from the weak form can be thought of as the equivalent (updated Lagrangian) finite element discretization where now the quadrature points are the material points which can move in space. Note that conservation of mass is automatically fulfilled as the total mass of the system $\sum _{p} m_p$ is the sum of all particle masses, which individually will remain constant in time.

The detailed steps of the numerical scheme are outlined in Algorithm 1. At each time step of the scheme, particle momentum is transferred to the grid on which the momentum update, (3.3), is performed. Boundary conditions are applied on the grid velocities, before the new grid velocities are transferred back to the particles. As such, boundary conditions are treated in a similar way as in finite element methods, avoiding expensive neighbour particle searches needed in, e.g. SPH. Before commencing the next time step, the particles’ stress states are checked against a yield criterion, and a return mapping is performed if necessary. Adaptive time stepping in this explicit scheme is accomplished by respecting the Courant–Friedrichs–Lewy (CFL) condition and the elastic wave speed condition (i.e. the elastic wave should not travel more than one grid cell in one time step). Accordingly, we use the constants $C_{CFL} = C_{el} = 0.5$, where the elastic wave constant $C_{el}$ is defined analogously to the CFL constant $C_{CFL}$.

Algorithm 1: Symplectic B-spline AFLIP MPM

It should be emphasized that there are several methods for transferring the momentum/velocities between the grid nodes and particles, and the chosen method can greatly influence the numerical dissipation of the scheme. With the time-dependent constitutive model proposed here, it is crucial to minimize this dissipation in order to capture the dynamics of the deformation correctly. Earlier versions of MPM used PIC- or FLIP-style transfers, i.e.

(3.5)\begin{equation} m^n_{i} \boldsymbol{v}^n_i = \sum_{p} m_p \boldsymbol{v}^n_p N_{i}(\boldsymbol{x}^n_p) \end{equation}

for the particle-to-grid momentum transfer and either (PIC)

(3.6)\begin{equation} \boldsymbol{v}_p^{n+1} = \sum_i \boldsymbol{v}_i^{n+1} N_{i}(\boldsymbol{x}_p^n), \end{equation}

or (FLIP)

(3.7)\begin{equation} \boldsymbol{v}_p^{n+1} = \boldsymbol{v}_p^n + \sum_{i} (\boldsymbol{v}_i^{n+1}-\boldsymbol{v}_i^n) N_i(\boldsymbol{x}_p^n) ,\end{equation}

for the grid-to-particle transfer. The PIC scheme is highly dissipative and does not preserve angular momentum in grid to particle, which could lead to rotational artefacts. While FLIP (as used, e.g. in the elastoplastic $\mu (I)$-rheology model of Dunatunga & Kamrin 2015) partly improves the conservation of angular momentum, it also reduces dissipation by only transferring the velocity increments, however, at the cost of more pronounced ‘ringing instability’ (Love & Sulsky 2006). Later, Jiang et al. (2015) proposed affine PIC (APIC) as an improvement of PIC, conserving angular momentum in both transfers, consequently resolving rotations correctly. The APIC particle-to-grid transfer is given by

(3.8)\begin{equation} m^n_{i} \boldsymbol{v}^n_i = \sum_p m_p ( \boldsymbol{v}_p^n + \boldsymbol{C}_p^n (\boldsymbol{x}_i^n - \boldsymbol{x}_p^n) ) N_{i}(\boldsymbol{x}^n_p), \end{equation}

with a ‘velocity derivative’ tensor $\boldsymbol {C}_p = \boldsymbol {B}_p \boldsymbol {D}_p^{-1}$ where $\boldsymbol {D}_p$ is analogous to an inertia tensor and can be shown to be a constant diagonal tensor given by $\frac {1}{4} (\Delta x)^2 \boldsymbol {I}$ in the case of a uniform background grid with quadratic B-splines (Jiang et al. 2015; Jiang, Schroeder & Teran 2017), and $\boldsymbol {B}_p$ is given by

(3.9)\begin{equation} \boldsymbol{B}_p^{n+1} = \sum_i \boldsymbol{v}^{n+1}_i (\boldsymbol{x}^n_i - \boldsymbol{x}^n_p)^T N_{i}(\boldsymbol{x}_p^n) . \end{equation}

In an effort to reduce dissipation further, APIC was recently generalized to AFLIP by Fei et al. (2021) by using the grid-to-particle scheme of FLIP, i.e. (3.7), instead of that of PIC, i.e. (3.6). Because of the low-dissipation property of AFLIP, this scheme is adopted in this work. Like in Fei et al. (2021), a parameter $\alpha _f$ is introduced which represents a PIC/FLIP ratio, with $\alpha _f=0$ meaning pure APIC and $\alpha _f=1$ meaning pure AFLIP. In particular, $\alpha _f=0.95$ is used to introduce as little numerical dissipation as possible. More details and examples of the low-dissipation property of AFLIP are presented in the paper of Fei et al. (2021), and the reader is also referred to the excellent and recent reviews by de Vaucorbeil et al. (2020), Sołowski et al. (2021) and Nguyen, de Vaucorbeil & Bordas (2023) for further details on MPM.

4. Results

In this section, various granular flow and collapse simulations using the critical state $\mu (I)$-rheology model are presented. In § 4.1, simulation results from flow on inclined planes are compared with approximate analytical solutions of the velocity profile. Next, comparisons with previous experiments of granular collapse problems are made in § 4.2, and cohesive media are considered in § 4.3. Finally, the ability of the model to capture the two steady-state solutions of granular flow over a smooth two-dimensional bump is demonstrated in § 4.4. In all simulations, the Young's modulus $E$, Poisson's ratio $\nu$ and hardening factor $\xi$ are kept constant with values shown in table 2. Furthermore, the initial compressive strength $p_c^0$ is taken to be 100 Pa and the parameter $\beta$ is set to zero, except in § 4.3 where cohesive simulations are presented. Appendix B presents further reasoning behind these constitutive parameter choices.

Table 2. Elastic and plastic model parameters used throughout this study, except in Appendix B where a sensitivity analysis of all parameters is presented.

4.1. Analytical solutions

Following the analytical treatment of GDR MiDi (2004), a cohesionless two-dimensional flow of height $h$ as sketched in figure 5 is considered. Under the assumption $\sigma _{xx} = \sigma _{zz}$ justified in da Cruz et al. (2005) and Andreotti et al. (2013), $q = \rho g (h-z) \sin \theta$ and $p = \rho g (h-z) \cos \theta$. Assuming all stress states lie on the CSL results in $\mu (I) = q/p = \tan \theta$. Combined with (2.26) and (2.27), using $\beta =0$ and assuming all deformation is plastic, one obtains

(4.1)\begin{equation} \frac{{\rm d} v_x(z)}{{\rm d}z} = \dot{\gamma}_{_S} = \omega \frac{\tan \theta - \mu_1}{\mu_2 - \tan \theta} \sqrt{\rho g (h-z)\cos \theta} {.} \end{equation}

Integrating over $z$, assuming a no-slip bottom surface, yields

(4.2)\begin{equation} v_x(z) = v_{max} ( 1 - (1-z/h)^{3/2} ), \end{equation}

where

(4.3)\begin{equation} v_{max} = v_x(h) = \frac{2}{3} \omega \frac{\tan \theta - \mu_1}{\mu_2 - \tan \theta} \sqrt{\rho g h^3 \cos \theta}\end{equation}

is the velocity at the surface of the flow.

Figure 5. Two-dimensional flow on a no-slip surface showing a Bagnold velocity profile. While (a) is a sketch of the set-up and flow, (b) shows the time evolution of the flow profile at a fixed $x$-position $2.5$ m downstream for a simulation with slope angle $\theta = 24^\circ$. Here, the black dashed line is the analytical solution from (4.2) where no fitting to the measured data has been made.

To see how the proposed model compares with the approximate analytical solution given by (4.2) and (4.3), two-dimensional simulations are performed assuming cohesionless spherical glass bead parameters given by Jop et al. (2005) based on the experiments by Forterre & Pouliquen (2003), i.e. $I_0=0.279$, $\mu _1=\tan (20.9^\circ )$ and $\mu _2=\tan (32.8^\circ )$, with initial density $\rho ^0=1550\ \textrm {kg}\ \textrm {m}^{-3}$ and intrinsic density $\rho _*=2500\ \textrm {kg}\ \textrm {m}^{-3}$ as summarized in table 3. In these simulations, a mass of 280 kg (per unit width) is being slowly fed onto a no-slip inclined surface and through a gate restricting the maximum height of the flow. The bead diameter is $d=2.5$ mm here, and the grid size $\Delta x$ is taken to be 1.2 mm. As expected, for slope angles $\theta$ smaller than $\tan ^{-1}(\mu _1) = 20.9^\circ$, the flow eventually stops, while for $\theta$ between $\tan ^{-1}(\mu _1)$ and $\tan ^{-1}(\mu _2)$, i.e. between $20.9^\circ$ and $32.8^\circ$, a steady state forms after some time. For the particular case of $\theta = 24^\circ$ at a point 2.5 m downstream, the plot in figure 5 shows how the velocity profile converges towards the Bagnold velocity profile given by (4.2) after a few seconds. The surface velocity measured for all steady-state flow cases is shown in figure 6(a), showing the agreement with the prediction by (4.3). As the flow transitions to a steady state, figures 6(b) and 6(c) show that $\mu$ approaches the expected constant value. Since the analytical solution assumes all stress states lie on the CSL, a slight deviation from the analytical solution is expected. The minor observed deviation can also be attributed to the non-idealized set-up where the particles are fed through a gate without a fixed flow height, resulting in a not perfectly smooth flow and constant flow height, and the plotted quantities are determined through an averaging scheme over a predefined region.

Table 3. Parameters used in § 4.1.

Figure 6. Simulations on varying slope angles $\theta$. In (a), the maximum velocity $v_{max}$ is measured along the surface of the flow and the analytic curve is (4.3) with $h = 10-0.3\theta$ cm fitted from the measured data. The error bars represent one standard deviation as the surface velocity was found by averaging over the spatial extent of the flow, neglecting only the area close to the gate and the area close to the flow front. In (b,c), the evolution of the average $\mu$, denoted $\langle \mu \rangle$, obtained from a slice of thickness 0.02 m in each simulation is shown. The latter is plotted against a scaled time, where $t_f$ denotes the time at which all mass has passed through the gate. The two dashed horizontal lines marks the lower and upper bounds $\mu _1$ and $\mu _2$, respectively.

In the steady-state flow simulations presented here, the solid fraction $\phi = \rho /\rho _*$ settles at an approximately temporally and spatially constant value for each slope angle. The hardening factor $\xi$ affects the compressibility of the system, here it is chosen large enough to only induce small changes in the solid fraction of the system. In particular, the solid fraction, initially being 0.62, varies between 0.60 and 0.63. The lowest fractions are measured for the highest slope angles where the flow is more rapid, qualitatively similar to the $\phi (I)$-relation suggested by GDR MiDi (2004), da Cruz et al. (2005) and Pouliquen et al. (2006). This is shown in figure 7, which was obtained by averaging the solid fraction in a slice in the middle of the flow. Initially, at $t/t_f<0.1$, the erratic behaviour of the measured average solid fraction is also attributed to the unstable nature of the first particles coming through the gate combined with the lack of particles needed for a stable average.

Figure 7. Evolution of the solid volume fraction $\phi$ with time as a steady-state flow is reached on different inclinations. As in figure 6, $t_f$ denotes the time at which all mass has passed through the gate.

Studying here the continuous (highly dynamic) flow where the deformations are dominantly plastic, the results presented are largely independent of the elastic moduli. Thus, whether $E$ is 1 MPa or 1 GPa, the results are indistinguishable. The same can be said about the exact onset of the plastic deformation as long as it occurs at sufficiently small stress values, in particular it does not matter in the presented results whether $p_c^0$ is 10 Pa or 100 Pa. This is demonstrated in Appendix B.

4.2. Granular collapse

Here, the problem of cohesionless granular collapse on horizontal and tilted surfaces is addressed. For the horizontal case, the experiments of Xu et al. (2017) are considered, where glass beads initially confined between two gates separated 8 cm apart were released by opening the gates such that the beads could flow between two parallel transparent glass sheets spaced 2 cm apart. The initial height of the system was 5 cm. Three-dimensional simulations are performed with the spherical glass bead parameters of § 4.1, using the known bead diameter $d = 2$ mm and intrinsic density $\rho _* = 2500\ \textrm {kg}\ \textrm {m}^{-3}$ reported from the experiments, as well as $\Delta x = 0.7$ mm ($N_p = 1.7$M). While the initial bulk density $\rho ^0$ was measured to vary between 1450 and $1625\ \textrm {kg}\ \textrm {m}^{-3}$, in the simulations it is for simplicity assumed that $\rho ^0 = 1550\ \textrm {kg}\ \textrm {m}^{-3}$, close to the average experimental value. The parameters are summarized in table 4. The walls are assumed frictionless and the ground no-slip, inferring the gate lift velocity $0.45\ \textrm {m}\ \textrm {s}^{-1}$ from the experimental photos. Figure 8 shows the results alongside the reported experimental observations, indicating an overall better agreement between the model and the experiments at time $t=0.42$ s (close to the rest state) than at the intermediate time. At the final time, note that the experimental data indicate that some granules spread out rather far from the main part of the system, resulting in long tails which were not captured by the simulations. This is believed to be due to the no-slip boundary condition at the bottom surface which in the experiments was a steel plate. Nevertheless, neglecting the tails, the repose angle of the main part of the collapsed system was reproduced remarkably well. The elastoplastic nature of the model naturally leads to the simulated granular collapse exhibiting localized bands of increased plastic (shear) strain. Shear bands are observed in granular flows (e.g. Mueth et al. 2000; Fenistein & van Hecke 2003; Jop et al. 2006; Mandal, Nicolas & Pouliquen 2021) and will be further discussed at the end of this section as well as in Appendices B and C. Regarding computational time, reaching the last state shown in figure 8 took almost 27 h on 8 cores (Intel Xeon Gold 6136 CPU 3.00 GHz, 187 GB RAM, OpenMP parallelization). In general, the computational time is, in addition to the degree of discretization, dependent on the Young's modulus and density as these parameters control the adaptive time stepping of the explicit scheme.

Figure 8. Granular collapse on a horizontal plane. (a) A middle slice from the simulations is shown and compared with the experiments from Xu et al. (2017) at different times. Note that the plotted slice is arbitrary as the front remains more or less uniform along the width as can be seen in the corresponding three-dimensional visualizations in (b). Here, the material is coloured according to the horizontal speed.

Table 4. Parameters used to simulate the experiment by Xu et al. (2017).

Having studied granular collapse on a horizontal surface, we now consider collapse on tilted surfaces. In the dam break experiments of Mangeney et al. (2010), cohesionless glass beads with diameter $d=0.7 \pm 0.1$ mm and intrinsic density $\rho _* = 2500\ \textrm {kg}\ \textrm {m}^{-3}$ were released on inclined rough surfaces. The system had an initial density $\rho ^0 = 1550\ \textrm {kg}\ \textrm {m}^{-3}$, initial height $h_0=14$ cm and initial length 20 cm along the inclined flow direction, enclosed by Plexiglas sidewalls 10 cm apart. These experiments are here simulated in three dimensions (using $\Delta x = 2.4$ mm, $N_p = 1.6$M) with a frictionless back wall preventing movement in the upflow direction, assuming a no-slip bottom surface as well as frictionless sidewalls. With the authors reporting the measured repose and avalanche angle as $23.5 \pm 0.5^\circ$ and $25.5 \pm 0.5^\circ$, respectively, $\mu _1$ is thus chosen here as the inverse tangent of the mean of these two angles. The remaining parameters, i.e. $I_0$ and $\mu _2$, are chosen as the typical spherical glass bead parameters as used for the collapse on the horizontal plane and in § 4.1. The parameters are summarized in table 5. With these assumptions, the final equilibrium state from the experiments is reproduced in the simulations. This can be seen in figure 9 for three different inclination angles $\theta = 10^\circ$, $16^\circ$ and $22^\circ$. The three-dimensional simulations required approximately 18 computational hours per real second on 8 cores on the same computer architecture as given above. While the agreement in figure 9 appears excellent, it should be emphasized that the simulations were conducted under the assumption of frictionless sidewalls. The sidewalls may influence the final deposit as demonstrated by Martin et al. (2017). With the rough bed surface being covered by the same beads glued to the surface, the no-slip assumption at the base is also questionable. In particular, Ionescu et al. (2015) argue that there is slip near the front.

Table 5. Parameters used to simulate the experiments by Mangeney et al. (2010). These parameters were also used in § 4.4 to simulate flow over a smooth obstacle.

Figure 9. Final rest state of the flow of a dam break on planes with three different inclinations. Simulations are compared with experiments from Mangeney et al. (2010).

Comparisons with the dam break experiments at intermediate times are shown in figure 10, which also highlights the static–fluid transition. This figure demonstrates that the simulations give very similar results compared with the $\mu (I)$-rheology model of Martin et al. (2017) without wall friction. Martin et al. (2017) also conducted simulations with wall friction, naturally resulting in a better match to the experimental data.

Figure 10. Dam break on two different inclinations, $\theta =10^\circ$ (a,c) and $\theta =22^\circ$ (b,d), at two intermediate times $t$ prior to reaching the final equilibrium position. The simulations are coloured in grey scale to highlight the static–flowing transition: dark grey is static and light grey is flowing. Experimental data from Mangeney et al. (2010) as well as previous simulations by Martin et al. (2017) with and without wall friction are included for comparison.

The three inclination angles in the dam break simulations are all below $\tan ^{-1}(\mu _1) = 24.5^\circ$, and thus the flow will come to rest, i.e. the front position converges to a fixed value. Figure 11(a) shows the overall capability of the model in capturing the evolution of the downslope front position $x_f$ in the dam break simulations. Comparing the front velocities, figure 11(b), highlights the not perfectly smooth front evolution of the simulated collapse on the largest inclination. On the two smaller inclinations, while still not achieving a perfectly smooth evolution, the front velocity remains very close to the experimental measurements. Figure 11(c) shows the evolution of the solid volume fraction $\phi$ with time, indicating a smaller final density with increasing inclination angle. After an initial increase in solid fraction (for $t<0.1$ s) due to consolidation under gravity, the solid fraction decreases gradually during the collapse. In the experiments, the system had already settled under gravity and the experimental data is therefore not expected to feature the initial jump in solid fraction observed in the simulations. Although the experimental data are scattered, especially for the largest inclination, the data seem to indicate an overall decrease in solid fraction during the collapse. This is qualitatively the same behaviour shown in the simulations after the initial consolidation.

Figure 11. (a) Front position $x_f$, (b) velocity $v_f$ and (c) solid fraction $\phi$ of the flow following the dam break on inclined planes. In (c), note that this is the solid volume fraction of the complete system and not a local measurement. Experimental data from Mangeney et al. (2010) and Martin et al. (2017).

A three-dimensional rendering of the $\theta =22^\circ$ dam break is shown in figure 12(a) visualizing the downslope velocity. While the final equilibrium state features a smooth surface, a non-smooth surface can be observed at intermediate times, even in the cohesionless case. This was not reported in the experiments. Visualizing instead the equivalent plastic shear strain rate in figure 13, it is immediately clear that the surface bumps are correlated with the occurrence of shear bands, i.e. thin zones of localized plastic shear strain. This is related to the initial state of the granular system and the consolidation that immediately follows due to gravity, thus producing states with large compressive strengths $p_c$ due to the imposed hardening law. In Appendix B, it is demonstrated that adjusting the plastic parameters, especially $\xi$, the gravity-induced initial compaction can be significantly reduced, and consequently it is possible to partly smooth out the surface bumps observed in these granular collapse simulations. Strain localization is expected in elastoplastic modelling featuring strain softening, and therefore the proposed model will fail in accurately capturing localization patterns without a rigorous regularization technique. Nevertheless, this is shown in Appendix C to not influence the general dynamics of the collapse and flow.

Figure 12. Dam break on slope angle $\theta =22^\circ$ of (a) cohesionless and (b) cohesive glass beads. The full time evolution of these two simulations is shown in supplementary movies 1 and 2 available at https://doi.org/10.1017/jfm.2024.643, respectively.

Figure 13. Cohesionless dam break on slope angle $\theta =22^\circ$ coloured according to the equivalent plastic shear strain rate $\dot{\gamma}_{_S}$.

4.3. The case of cohesion

While the above experiments and simulations were conducted with cohesionless glass beads, it is instructive to consider the case of cohesive beads. To this end, the dam break case on the slope angle $\theta =22^\circ$ is studied with $\beta >0$. In figure 12, snapshots of the dam break at different times are presented for the case $\beta = 0$ and $\beta = 0.3$, allowing us to visually observe the difference in the evolution of the flow when cohesion is introduced. In particular, the extension of the flow is reduced and lumps of material adhere to each other at the surface of the flow. The two simulations presented in figure 12 are also shown in supplementary movies 1 and 2, respectively, at $\frac {1}{3}$ speed. Figure 14 shows how the evolution of the front position decreases with increasing $\beta$.

Figure 14. Front position $x_f$ of flow on inclined plane $\theta =22^\circ$ with increasing cohesion $\beta = 0$, 0.1, 0.2, 0.3.

The recent publication of Gans et al. (2023) presented detailed experiments of granular collapse with cohesive glass beads produced using a novel technique by Gans, Pouliquen & Nicolas (2020) of coating the beads with a thin layer of polyborosiloxane. Since the coating provides a stable cohesion in the form of a constant-in-time inter-granular adhesive force, these experiments can be considered suitable for assessing the proposed model, which, in its current form, assumes a persistent (unbreakable) cohesion. The experiments were conducted in a 15.4 cm wide channel having a rough bed constructed by gluing identical grains as the flowing grains, the granular system having an equal initial height $h_0$ and length of 8.9 cm. Two of their experiments, one with slightly more cohesive beads than the other, are presented in figure 15, where we call the least cohesive case Exp. A and the most cohesive case Exp. B. In the case of Exp. A, the surface of the system is plotted at two different times. The figure also includes results from two-dimensional simulations where the spherical glass bead parameters as presented in § 4.1 were utilized besides from those reported from the experimental measurements, in particular $d=0.8$ mm, $\rho _* = 2600\ \textrm {kg}\ \textrm {m}^{-3}$ and $\rho ^0 = 1508\ \textrm {kg}\ \textrm {m}^{-3}$. These parameters are summarized in table 6. As in the previous section, a no-slip ground was assumed. In the presented simulations the cohesion through $\beta p_c^0$ is tuned to match the experiments, here using $p_c^0=500$ Pa. As before, it is clear that runout is decreased with increasing cohesion. In particular, using $\beta =0.2$ and $\beta =0.3$ can be seen to give a front position evolution that aligns with the experiments, albeit not perfectly. As in the previous section, it should be mentioned that the assumption of a no-slip base and frictionless walls may not be fully accurate. Moreover, the previous comments made about strain localization in the cohesionless case also applies here, and the reader is again referred to Appendices B and C. More importantly, since the cohesion in the proposed model is tuned to match the experiments, it is pertinent to emphasize that the presented simulations do not constitute a rigorous validation of the model.

Figure 15. Cohesive granular collapse simulations compared with the experiments from Gans et al. (2023) on a horizontal rough surface. (a,b) Snapshots of the simulation of Exp. A using $\beta =0.2$ at two different times where the markers indicate the surface measured in the experiments. (c) Evolution of the front position $x_f$ with time in both experiments are shown together with simulations with $\beta =0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3$.

Table 6. Parameters used to simulate the experiments by Gans et al. (2023) in figure 15.

Having investigated the effect of cohesion on granular collapse, one can also consider more prolonged flows. From discrete simulations, it has been demonstrated that cohesive granular flows can reach a steady state where the introduction of cohesive forces between the grains leads to the occurrence of a plug flow near the free surface (Rognon et al. 2008b; Mandal et al. 2020). It appears possible for the proposed critical state $\mu (I)$-rheology model to qualitatively capture this behaviour, as presented in figure 16. Here, on an inclined slope of $\theta = 32^\circ$ with a no-slip base, steady-state flows are reached for various values of $\beta$. With increasing $\beta$ the occurrence of a plug flow is observed, and above a critical $\beta$ the flow stops. In these simulations, the same set-up and parameters as in § 4.1 were used (including $p_c^0 = 100$ Pa), except with a twice as large gate height and $\mu _2 = \tan (40^\circ ) > \tan (\theta )$ to ensure steady-state flow. These parameters are listed in table 7. A comparison with the results of the discrete simulations of Mandal et al. (2020) reveals clear differences, both in terms of the shape of the velocity profile under the plug and in terms of the relative change in velocity magnitude as the cohesion is increased. While the qualitative results from the cohesive simulations presented in this section are promising and in line with expectations, a robust validation of the cohesive model and its implementation remains to be addressed.

Figure 16. Steady-state velocity profiles for increasing cohesion parameter $\beta$ on a no-slip surface tilted at $\theta = 32^\circ$. The velocity profiles are measured at the same position $x=0.8$ m from the gate which had a height of $40d$.

Table 7. Parameters used to simulate cohesive steady-state flow (figure 16).

4.4. Multiple solutions for flow over a smooth bump

In small-scale experiments of granular flow over a smooth obstacle, Viroulet et al. (2017) showed that two very different steady-state regimes can occur. In particular, either a ‘jet’ flow detaching from the apex of the bump or a ‘shock’ forming upstream of the bump is observed, the latter being the case if a sufficient amount of erodible particles are placed in front of the bump. This initially static mass reduces the energy of the impacting flow and accordingly also the flow velocity, generating an increased flow thickness before the bump. Curious as to whether the proposed critical state $\mu (I)$-rheology model can capture both steady-state solutions, simulations are performed assuming beads with the same diameter and on the same terrain/bump as used in the experiments. Along the flow direction $x$, the elevation of the smooth symmetric bump is given by (in metres)

(4.4)\begin{equation} b(x) = 0.0475 \,{\rm sech}\left(\frac{x-0.43}{0.04} \right), \end{equation}

as sketched in figure 17. At $x=0$, i.e. 43 cm from the centre of the bump, the beads flow through a gate of height 1.5 cm. With the bump being two-dimensional and the flow being confined in a narrow chute, the problem is here simulated in two spatial dimensions. Cohesionless spherical glass beads with identical parameters as in the dam break simulations of § 4.2 are assumed, i.e. the parameters presented in table 5. In the experiments of Viroulet et al., the beads show slippage along the base, prompting the introduction of a base friction in the simulation, for simplicity as a Coulomb friction law approximated on the velocities. In particular, introducing a base friction parameter $\mu _b$, the velocity $\boldsymbol {v}$ relative to the boundary is obtained as

(4.5)\begin{equation} \boldsymbol{v} = \begin{cases} \boldsymbol{v}^*_{T} - \mu_b ||\boldsymbol{v}^*_{N}||_{L_2} \dfrac{\boldsymbol{v}^*_{T}}{||\boldsymbol{v}^*_{T}||_{L_2}} & \textrm{if}\; ||\boldsymbol{v}^*_{T}||_{L_2} > \mu_b ||\boldsymbol{v}^*_{N}||_{L_2} {,} \\ \boldsymbol{0} & \text{otherwise}, \end{cases} \end{equation}

Figure 17. Sketch of terrain in experiments of Viroulet et al. (2017).

where $\boldsymbol {v}^*_{T}$ and $\boldsymbol {v}^*_{N}$ refer to the tangential and normal velocity components, respectively, of $\boldsymbol {v}^*$, i.e. the relative velocity before the boundary condition is considered. Equation (4.5) is only used if $\boldsymbol {v}^*_N$ suggests penetration of the boundary, otherwise $\boldsymbol {v} = \boldsymbol {v}^*$. Such a boundary condition (4.5), which has been extensively used in MPM flow modelling before (Gaume et al. 2018; Li et al. 2020; Cicoira et al. 2022) and which is further validated in Appendix E, is motivated from having to apply boundary conditions on the grid where the velocities are readily available, unlike the stresses. It is reported from the experiments that, with no initial mass in front of the bump, the beads detach from the apex of the bump for slope angles around approximately $35^\circ$, which justifies roughly $\mu _b=0.5$. As in the experiments, the problem is here simulated with and without an initial mass (here 0.13 kg) in front of the bump. Figures 18 and 19 show that the simulations at a slope angle $\theta =39^\circ$ capture both of the two different steady-state regimes. In the case with initial mass in front of the bump, figure 18 shows the ‘shock’ forming upstream, demonstrating an overall good agreement with the experimental observations. The complete time evolution in this simulation can be seen in supplementary movie 3. Without initial mass in front of the bump, figure 19 shows the detaching ‘jet’ steady-state flow. In this case, the experimental images indicate a significant vertical spreading of the particles, a situation which is not observed in the simulations. This is believed to be due to the sidewalls in the experiments, where the wall friction reduces the speed of the particles closer to the walls. Consequently, these particles will experience a shorter trajectory over the bump compared with the particles in the middle of the flow. As the experimental images are taken from the side through the transparent walls, the different trajectories are not visually distinguishable. Three-dimensional simulations with sidewall friction would allow for a more proper comparison with experiments, and a more extensive parameter tuning could improve the agreement with experimental trajectories. Using a grid size of $\Delta x = 0.6$ mm and $N_p = 0.65$M, the two-dimensional simulations presented here required approximately 8 computational hours per second on 8 cores on the same computer architecture as presented in § 4.2. Supplementary movie 4 shows the complete time evolution in the jet flow simulation. In this movie, it can be noted how the flow slightly pulses in response to small waves that are propagating down from the inflow. This is here due to the not perfectly smooth feeding of particles through the gate, but could also be a physical result of the flow breaking down into roll waves on sufficiently long chutes (Barker & Gray 2017; Viroulet et al. 2018). Nevertheless, it is not observed in the experiments.

Figure 18. Flow over bump with initial mass and $\theta =39^\circ$. Experiments of Viroulet et al. (2017) (left) and simulations (right) at increasing time from top to bottom ($t = 0, 0.4, 0.7, 1, 1.5, 4$ s). The full time evolution is shown in supplementary movie 3. Note the colour bar only applies to the simulations, not the experiments.

Figure 19. Flow over bump without initial mass and $\theta =39^\circ$. Experiments of Viroulet et al. (2017) (left) and simulations (right) at increasing time from top to bottom ($t = 0.3, 0.6, 0.9, 4$ s). The full time evolution is shown in supplementary movie 4. Note the colour bar only applies to the simulations, not the experiments.

5. Discussion and conclusions

This paper has presented a continuum model based on the combination of the $\mu (I)$-rheology with the critical state concept from geo-mechanics in a way that captures both solid and fluid properties of granular media in a single framework. Incorporating the elastoplastic CSSM theory, the proposed model is different from the previously proposed critical state $\mu (I)$-rheology models for dry granular flows (Barker et al. 2017; Schaeffer et al. 2019; Bouchut et al. 2021). While these previous models are either implemented in depth-averaged schemes (in one or multiple layers) or realized only for one-dimensional problems, the elasto-viscoplastic model of this paper is implemented in a B-spline MPM framework allowing full three-dimensional simulations. Accounting for cohesion and compressibility through CSSM, the proposed model is more general than the two-dimensional MPM-implemented $\mu (I)$-rheology of Dunatunga & Kamrin (2015). The proposed model approximately reproduces the expected flow profiles on inclined slopes, and with the introduction of a smooth bump along the slope, the model is able to qualitatively provide the two steady-state regimes observed in the experiments of Viroulet et al. (2017). While for cohesionless glass beads a Drucker–Prager yield may be more appropriate, for other granular materials, in particular snow, closed yield surfaces (such as the MCC yield surface) are more appropriate (Reiweger et al. 2015; Ritter et al. 2020; Blatny et al. 2021). Snow also undergoes significant volumetric deformations under typical loading conditions compared with cohesionless glass beads. The MCC model already facilitates an uncomplicated treatment of such compaction and dilation.

The simulations of granular collapse of spherical glass beads displayed shear bands affecting the smoothness of the surface at intermediate times before reaching the final equilibrium position. This was attributed to the initial state of the granular system and the immediate consolidation occurring as a result of the onset of gravity, thus leading to a sudden increase in $p_c$. As the gravity-induced compaction can be suppressed through adjustments of the plastic parameters of the hardening law, Appendix B demonstrated that the surface during the collapse can be smoothed. While such shear band-induced surface patterns were not reported in the spherical glass bead collapse experiments, they are generally observed in experiments involving dry granular media (Gudehus & Nübel 2004; Torres-Serra, Romero & Rodríguez-Ferran 2020). Nevertheless, neglecting the surface effects, §§ 4.2 and 4.3 showed that the proposed CSSM-based $\mu (I)$-rheology model is able to capture the observed dynamics and final deposit in granular collapse experiments of spherical glass beads on both flat and inclined surfaces reasonably well.

While Lacaze & Kerswell (2009) showed, using discrete simulations, that the $\mu (I)$-rheology holds across the whole flow of such problems, others have presented numerical experiments that seem to indicate that the $\mu (I)$-rheology is not needed to capture the specific problem of granular collapse. Specifically, the recent paper of Rousseau et al. (2022) claims that only a single friction parameter is needed to describe granular collapse and further argues that the $\mu (I)$-rheology is an over-complication for modelling such problems. This is analogous to the previous claims of Ionescu et al. (2015) and Martin et al. (2017) who, implementing the $\mu (I)$-rheology in a Drucker–Prager viscoplastic approach with variable viscosity, showed that the results did not differ significantly if the viscosity was assumed constant. Similarly, Crosta, Imposimato & Roddeman (2009) showed that granular collapse can be accurately modelled with a constant friction angle in a Mohr–Coulomb yield rule. Along these lines, it seems reasonable to assume that, with appropriately chosen constitutive parameters, a rate-independent CSSM (i.e. with a constant, $I$-independent, slope of the CSL) may also be adequate in modelling the granular collapse problem. This is shown and discussed in Appendix D.

Moreover, it should be mentioned that the wall friction may influence the results. In the simulations of granular collapse presented here (on various inclinations), all sidewalls were assumed frictionless. While investigating the influence of wall friction in granular flow was not within the scope of this study, such simulations could be accomplished by employing the frictional boundary conditions outlined in § 4.4. The effect of sidewalls during granular collapse has been analysed, e.g. numerically by Martin et al. (2017), showing that introducing friction may help in better capturing the experimental results. Along the same lines, slippage of the material along the base may also influence the results. The base was assumed no-slip in all granular collapse simulations presented here, and it was highlighted as a potential reason for the short runout compared with the experimental observations by Xu et al. (2017) of granular collapse on the horizontal surface. In particular, Ionescu et al. (2015) reported a 10 % increase in runout going from a no-slip base to a basal friction of 0.48. The reader is referred to these references for further insight into the influence of wall and base friction on granular collapse. The flow over bump problem presented in § 4.4 was the only set-up in this work where a non-zero finite base friction was employed. Following the MPM community, this was accomplished by applying boundary conditions on the grid velocities, proving an approximation to the Coulomb friction law which was further validated in Appendix E. While one could foresee imposing boundary conditions on the stresses, this would come at an increasing computational cost as the stresses are not readily available on the grid.

The cohesion assumed in the proposed model is a form of cohesion that is always present (for $\beta$, $p_c > 0$). In order to apply the model to flows where this may not be an appropriate description of cohesion, the model could be altered to remove cohesion, i.e. setting $\beta$ to zero, after yield or another pre-defined criterion. With increasing cohesion, it was found that steady-state flow can still be observed, with increasing $\beta$ leading to plug flow near the surface, qualitatively similar to the observations in the discrete simulations of Rognon et al. (2008b) and Mandal et al. (2020) of cohesive granular flows. An exact quantitative comparison between the proposed model and previous discrete simulations remains challenging as this entails relating $\beta$ to the prescribed inter-particle force parameter used in the discrete simulations. In this study, the cohesion parameter $\beta$ imposes a translation of the CSL in principal stress space. This is analogous to the usual treatment of cohesion in viscoplastic fluids, such as Bingham fluids, where cohesion can be accounted for by altering the yield stress. For example, Abramian, Staron & Lagrée (2020) and Gans et al. (2023) considered cohesion similarly by altering a simple $I$-dependent Mohr–Coulomb yield criterion in their volume-of-fluid incompressible Navier–Stokes scheme. Comparing with discrete simulations, Abramian et al. (2020) concluded that such approach can accurately capture the effective cohesion induced by cohesive discrete particles. In fact, they proposed explicit relations between their (macroscopic) Mohr–Coulomb cohesion parameter and the (microscopic) inter-particle cohesive force parameter. Such a study could also be performed in the context of the model proposed here, in order to establish a link between $\beta$ (and $p_c^0$) and the inter-particle attractive force. On the contrary, the numerical studies by Mandal et al. (2020, 2021) suggest that, while the inter-particle attractive force controls the initiation of the flow, the cohesive flow dynamics is better described by an effective attractive force that also takes into account the stiffness and inelasticity of the grains.

In the critical state $\mu (I)$-rheology model presented here, the granular material can undergo both elastic and plastic compaction and dilation. Unless the initial compressive strength is chosen to be very large, the elastic volumetric change is negligible and all volumetric deformation can be considered plastic. In § 2.4, this assumption was used to derive a hardening law that recovers the logarithmic relationship between the inverse solid volume fraction and the isotropic compressive strength, which is a relationship commonly used for the (slow) deformation of soils and powders (Walker 1923; Poquillon et al. 2002; Castellanos 2005). However, as mentioned in the Introduction, for granular media undergoing flow, a $\phi (I)$-relation has been reported, relating the solid volume fraction to the inertial number. In particular, above a critical inertial number, the solid fraction typically decreases linearly with this number. A hardening law that recovers the $\phi (I)$-relation may be possible, although this was not explored in this work, where we opted for a $\phi (p)$-relation. However, as briefly mentioned in §§ 4.1 and 4.2, this did induce a decreasing solid volume fraction with flow speed, qualitatively similar to the $\phi (I)$-relation. Formulating an accurate unified hardening law for both solid and fluid granular media, i.e. a law reproducing $\phi (I)$ in the high inertial number limit and $\phi (p)$ in the low inertial number limit, would require extensive further work.

The critical state $\mu (I)$-rheology assumed in this paper is a local theory in the sense it does not explicitly account for the observed situation where stress or velocity fluctuations caused by rearrangement at one location trigger rearrangements elsewhere in the system (Gaume, Chambon & Naaim 2011, 2020). This is what led to the non-local theories of Pouliquen & Forterre (2009) and Kamrin & Koval (2012). In the future, the model presented in this paper could be extended to introduce non-local effects in order to reproduce, e.g. the previous experimental observations of Pouliquen (1999), Pouliquen & Forterre (2002) and GDR MiDi (2004) indicating that there is a minimum thickness required to observe a steady-state flow at a given inclination angle. This minimal thickness decreases with inclination angle. In the local $\mu (I)$-rheology this cannot be reproduced without possibly considering a threshold on the flow speed. Very recently, Dunatunga & Kamrin (2022) and Haeri & Skonieczny (2022) included non-local effects in their elasto-viscoplastic Drucker–Prager $\mu (I)$-rheology model.

The flow and collapse simulations highlighted in this paper provide a small glimpse at what is feasible with the presented theoretical and computational framework. In the future, it may be applied to various granular dynamics problems. For example, the numerical implementation in MPM makes this framework particularly suited to study erosion (Li et al. 2022, 2024; Vicari et al. 2022). While this may be accomplished well with other methods (Lusso et al. 2017; Fernández-Nieto et al. 2018), in MPM the terrain can simply be filled with the same or different material as that of the incoming flow, with no extra coupling or changes to Algorithm 1 needed. Other phenomena that could be studied include the formation of (free-surface) roll waves which may develop in granular flow on rough inclined planes (Forterre & Pouliquen 2003; Barker & Gray 2017; Viroulet et al. 2018). Implementing the model in other hybrid/mesh-free numerical solvers such as SPH or PFEM should also allow such studies. In addition, future studies could investigate the capability of the model, or extensions of the model, to capture levee formation, self-channelization (Félix & Thomas 2004; Mangeney et al. 2007; Rocha, Johnson & Gray 2019; Edwards et al. 2023) and finger formation (Woodhouse et al. 2012; Baker, Johnson & Gray 2016) observed, e.g. in debris flows, snow avalanches and pyroclastic flows (Jessop et al. 2012). One could also foresee coupling the model with particle size segregation in order to capture the segregation phenomena that occur in poly-disperse granular assemblies (Rognon et al. 2007; Barker et al. 2021; Maguire et al. 2024). While an increased interest in understanding the effect of cohesion in granular flows has been witnessed (Rognon et al. 2006, 2008b; Berger et al. 2015; Khamseh et al. 2015; Badetti et al. 2018; Abramian et al. 2020; Vo et al. 2020; Mandal et al. 2020, 2021; Macaulay & Rognon 2021; Dong et al. 2023; Gans et al. 2023), this topic has also recently become increasingly relevant in the snow and avalanche community (Rognon et al. 2008a; Bartelt et al. 2015; Li et al. 2020). This is particularly motivated by the notion that climate change induces an increase in the proportion of wet cohesive snow avalanches compared with dry ones (Ballesteros-Cánovas et al. 2018; Součková et al. 2022). Hence, we expect the proposed approach to have important impacts in this field, in particular for the evaluation of avalanche impact pressures, runout distances and more generally risk management and mitigation (Ortner et al. 2023).