Granular viscosity from plastic yield surfaces: The role of the deformation type in granular flows

Numerical simulations of granular flows with Navier–Stokes type models emerged in the last decade, challenging the well established depth-averaged models. The structure of these equations allows for extension to general rheologies based on complex and realistic constitutive models. Substantial effort has been put into describing the effect of the shear rate, i.e. the magnitude of the velocity gradient, on the shear stress. Here we analyse the effect of the deformation type. We apply the theories of Mohr–Coulomb and Matsuoka–Nakai to calculate the stresses under different deformation types and compare results to the theory of Drucker–Prager, which is formulated independently of the deformation type. This model is particularly relevant because it is the basis for many granular rheologies, such as the µ I ( ) –rheology. All models have been implemented into the open-source toolkit OpenFOAM ® for a practical application. We found that, within the context of these models, the deformation type has a large influence on the stress. However, for the geometries considered here, these differences are limited to specific zones which have little influence on the landslide kinematics. Finally we are able to give indicators on when the deformation type should be considered in modelling of landslides and when it can be neglected.


Introduction
Dense granular flows are substantial parts of many natural hazards, such as avalanches, landslides, debris flows and lahars.A constitutive description of granular materials is important for a deeper understanding and improved prediction of these processes.
The first models for granular materials stem from geotechnics and applications in the soil mechanics community, with the earliest examples of a mathematical description being in the 19th century, when Charles-Augustin de Coulomb formulated his famous friction law (see e.g.[1]), based on the rules found earlier by Guillaume Amontons (see e.g.[2]).It relates normal stress n and shear stress between two solids and defines the friction angle as = tan( ) .
n (1) This relation is limited to well defined sliding planes and not generally applicable to three-dimensional situations.Christian Otto Mohr generalized Coulomb's law by determining the decisive shear plane in a failing solid with arbitrary stress state.The resulting formulation is known as the Mohr-Coulomb (MC) failure criterion [3].The Mohr--Coulomb failure criterion describes the failure of brittle materials, such as concrete and granular material with good accuracy.In fact, it has been successfully applied to a wide range of problems, especially within elasto-plastic frameworks [4].Various extensions have been applied to the original idea, introducing strain hardening and softening or socalled caps, limiting the admissible pressure [5].Matsuoka and Nakai [6] proposed a smoother version of the Mohr-Coulomb criterion, the Matsuoka-Nakai (MN) failure criterion.It has gained a lot of popularity, as it improved numerical stability as well as physical accuracy.
All developments were merged into modern constitutive models for quasi-static granular materials, such as the hardening-soil-model [5], Severn-Trent-sand [7], SaniSand [8], Hypoplasticity [9,10] or Barodesy [11][12][13].All of the above models assume quasi-static conditions and describe the stresses at failure.The respective extension to dynamic cases has been of interest for a long time.Schaeffer [14] was the first to combine the two-dimensional Navier-Stokes Equations, with the quasi-static Mohr-Coulomb failure criterion.His work assumed that the Mohr-Coulomb failure stress is valid beyond failure and at any strain rate, an assumption which is consistent with ideal plasticity [4], where failure criterion and yield criterion are coinciding.In this sense we will use henceforth the term yield criterion when computing the stress in a flowing material.Schaeffer's theory matches Mohr-Coulomb only for plane strain (i.e. two dimensional deformations) and incompressible flows, i.e. isochoric shear.Although he found some problematic instabilities, his work was highly influential in the granular flow community and his approach was applied for both, plane strain and fully three-dimensional granular flows, first and foremost in chemical engineering [15,16].
Early models assumed that the magnitude of the stress tensor remains constant after failure and independent of the shear rate.This assumption fails to describe various phenomena from physical experiments, e.g.steady flows on inclined planes or roll waves [17][18][19][20] and a large amount of research has been attributed to describe the respective correlations.Most notable are the early works of Bagnold [21,22] and Voellmy [23], both combining Coulomb's friction law and an additional dynamic term.More recent developments have been achieved with kinetic theory [24,25] and with the so-called µ I ( )-rheology, relating the dimensionless friction coefficient = µ sin( ) to the dimensionless inertial number I [18].
The first generally applicable µ I ( )-rheology was introduced by Jop et al. [26], basically following the approach of Schaeffer [14] but with a shear rate dependent yield criterion.The instabilities found by Schaeffer [14] are partially present in the model of Jop et al. [26] and addressed, among others, by Barker et al. [27,28].
Although the success of respective models is impressive, we have to contemplate that Schaeffer's approach is limited to plane, incompressible deformation, i.e. isochoric shear.If his relation is used for arbitrary three-dimensional deformations, one gets what is commonly known as the Drucker-Prager (DP) yield criterion in solid mechanics [29].The differences between Mohr-Coulomb and Drucker-Prager are well known [30] and may be rather large, depending on the induced deformation [31,32].This yields, for example, remarkable errors in earth pressure calculations [33].Furthermore, it has been revealed that the Mohr-Coulomb yield criterion is fulfilled in discrete element simulations to a much better degree than the Drucker-Prager yield criterion [34].
These circumstances lead to a big gap between Mohr-Coulomb models, and the recent success of the µ I ( )-rheology.This gap is what we aim to close with this work.We will extend the approach of Schaeffer [14] and Jop et al. [26] to different yield criteria.This allows us to implement constitutive models based on Mohr-Coulomb and Matsuoka-Nakai alongside the classic relation of Schaeffer [14] (i.e.Drucker-Prager) into the CFD-toolkit OpenFOAM® [35].We are comparing the three relations by inducing three-dimensional deformations in granular flow simulations.We neglect the shear rate dependence of stresses in this work to focus solely on the effect of the yield criterion.However, the µ I ( )-scaling is perfectly compatible with all of the presented yield criteria and can be easily reintroduced.To keep computational demand to an acceptable level, we run axisymmetric cases, some of which have been validated with physical experiments [36].We are able to draw some conclusions and determine the circumstances for which the widely used Drucker-Prager relation differs from the traditional Mohr-Coulomb relation.

Method
The incompressible Navier-Stokes Equations are given as with velocity u, density , stress tensor , strain rate tensor 1 D, pressure = p T 1/3 tr( ) and gravitational acceleration g.In the following, we will assume that pore-pressure is negligibly small and hence that the pressure p is equal to the effective pressure p .
In order to combine granular rheologies with the Navier-Stokes Equations, it is required to calculate the deviatoric stress tensor T dev in terms of the known strain rate tensor, and the dynamic viscosity2 , or equivalently the kinematic viscosity = / (note the first term on the right-hand side of Eq. ( 3)).Therefore, without substantial modification of the underlying equations, all constitutive modelling is reduced to the scalar viscosity.This introduces some severe limitations, namely alignment between stress and strain rate tensor, which we will illuminate below.The hydrostatic part of the stress tensor, p I can not be established with a constitutive relation in incompressible flow models (as = D tr( ) 0) and is therefore calculated from the constraint of the continuity Eq. ( 2) with e.g. the PISO algorithm [37].

A tensorial description of deformation and stress
The set of all stress tensors, at which yielding occurs forms a surface in stress space, which we define as the yield surface (see Fig. 4).Admissible stress tensors of static material are located within the yield surface and stress tensors do not exist outside the yield surface [4].Vice-versa, a yield surface, in addition with the so called flow rule, is sufficient to describe a perfectly plastic material.As we will show later, the flow rule is required to determine the location on the yield surface and thus the stress tensor.
Many yield surfaces are defined in terms of stress invariants and we will shortly introduce the most important notations and relations.The first invariant is the trace of a tensor, (5) and the trace of the stress tensor is proportional to the pressure The second invariant is defined as and related to the norm of a tensor, which is defined here as 3 For the incompressible strain rate tensor we can thus write and the same is the case for the deviatoric stress tensor Finally, the third invariant is defined as and is, in combination with the second invariant, related to the type of deformation.The type of deformation can be described in terms of the Lode-angle [38], defined as symbol is reserved for the friction coefficient in granular flows and we will use instead. 3Note that often the Frobenius norm is used instead, defined as for symmetric tensors.
A physical interpretation of the Lode-angle is shown in Fig. 1.Deformation types can be reduced to three cases for incompressible flows: True triaxial compression is defined by one negative and two positive principal strain rates and characterized by a negative Lode-angle.Isochoric shear is the combination of plane strain (one principal strain rate is zero) and incompressibility, the Lode-angle in this case is zero.Finally, a positive Lode-angle indicates true triaxial extension, which is defined by two negative and one positive principal strain rates.Moreover, the two limiting cases shown in Fig. 3 are called triaxial compression and extension.Note that in reality, deformation types are much more complicated due to e.g.volumetric deformations.The Lode-angle as defined by Eq. ( 12) is limited to the interval °° [ 30 , 30 ], respectively.For an interpretation in principal stress space as deviatoric polar angle (Fig. 3), the Lode-angle is located in one of six equal sectors, depending on the order of the principal stresses (Haigh-Westergaard coordinates) [38].For isotropic constitutive models this plays no role and Eq. ( 12) can be used throughout all derivations.
The three invariants are sufficient to describe a tensor, except for its orientation in space (i.e.rigid body deformation).However, a rigid body deformation is not relevant for isotropic constitutive models.This makes invariants very convenient for defining yield criteria and constitutive models.

Stress tensor reconstruction
It is common to assume that the deviatoric stress and strain rate tensors are co-axial, meaning that their eigenvectors are parallel [14,39,40] where T i dev and D i are the eigenvectors of the respective tensor.This holds at least in steady, critical state [41] and allows for visualisation of the principal stress and strain rate tensors in the same principal-valuespace, which is for example utilized in Fig. 2.However, this assumption is not sufficient to determine the complete stress tensor.
A stricter assumption, which is sufficient to determine the complete stress tensor is alignment between deviatoric stress and strain rate tensor [26,40].Alignment of tensors implies co-axiality and an equal ratio between eigenvalues.This can be expressed in terms of principal stresses and strain rates as or more generally as This assumption is commonly applied in granular flows, among others by the µ I ( )-rheology [26].It allows for expression of the deviatoric stress tensor as the product of a scalar and the strain rate tensor and to establish a relation compatible with the Navier-Stokes Equations, Note that the viscosity is not constant and depends on the pressure and the strain rate (we will additionally introduce a dependence on the Lode-angle).Alignment is fundamental for many granular flow theories and its validity is indicated by good experimental agreement [17,26].
The validity of alignment can also be checked with models that include non-alignment between stress and strain tensor.We will use the constitutive soil model Barodesy [11,13] in here, as it has shown to provide a realistic relation between the directions of stresses and strain rates at critical state [42].Note that the intrinsic non-alignment of Barodesy makes it impossible to express it with a scalar viscosity, making it incompatible with the Navier-Stoskes-Equations.The stress in critical state follows from Barodesy as with The only constitutive parameter is usually expressed by the friction angle but in here used to match Barodesy to the Matsuoka-Nakai criterion, , where k MN is the constitutive parameter of the Matsuoka-Nakai cri- terion, see Eq. (34).Note that the factor 2 in Eq. ( 18) is related to the definition of the norm following Eq.( 8).Furthermore, plasticity theory allows us to classify the popular approach of alignment in a broader context.The associated flow rule, states that the plastic strain rate tensor (in here simply D, as we assume ideal plasticity) is oriented normal on the yield surface f T ( ).This re- lation follows the principle of maximal plastic dissipation [4] but the yield surface f T ( ) can be replaced with arbitrary plastic potentials f T ( ) in the flow rule (Eq.( 20)).In such a case one speaks of a non-associative flow rule.This procedure is applied often, as this matches soil behaviour better [43][44][45].
Inserting the definition for alignment (Eq.( 16)) into the flow rule (Eq.( 20)) yields which necessarily leads to a von Mises type surface for the plastic potential, with constant k vM .Conveniently, this flow rule is consistent with the continuity Eq. ( 2).Note, that associated flow rules for Drucker-Prager, Mohr-Coulomb and Matsuoka-Nakai predict volumetric strain, which is contradicting with Eq. ( 2) and experimental observations in critical state.This gives three methods to determine the direction of the stress tensor with respect to the strain rate tensor: Associated flow rule (20) and von Mises plastic potential (i.e.alignment), as well as sophisticated soil models, in here Barodesy.All three approaches are presented and compared in Fig. 2 for the Matsuoka-Nakai yield surface, as this matches Barodesy best [12].Results overlap in the deviatoric plane for triaxial extension and compression but differ substantially for shear.
The ad hoc assumption of alignment results in strain rate directions close to the associative flow rule and Barodesy, two physically reasonable models.In particular, the alignment assumption fits Barodesy and thus experimental behaviour, better than the associated flow rule.
In the following we will apply a von Mises plastic potential, as this allows simple implementation of various yield surfaces into the Navier-Stokes equations.Moreover, this approach guarantees a well defined stress tensor for all relevant yield surfaces and is compatible with incompressible flows.

Yield criteria
One of the simplest yield criteria, the von Mises yield surface has been introduced in Eq. (22).It represents basic plastic behaviour and is as such the basis for visco-plastic rheologies like Bingham or Herschel-Bulkley [46].However, it only takes into account the second stressinvariant and is thus not able to cover basic granular behaviour, i.e. pressure-dependent shear strength.
The Drucker-Prager yield criterion is the simplest criterion that takes the frictional character of granular materials into account [29].As such, it connects the pressure with the deviatoric part of the stress tensor, with the friction angle .Note that this is a special parametrisation for a cohesionless material in terms of the friction angle , such that Drucker-Prager and Mohr-Coulomb match for a Lode-angle = 0, i.e.
isochoric plane strain [14].This yield surface allows us, in combination with the von Mises flow rule or equivalently Eq. ( 16), to define the viscosity as It is also possible to model the friction coefficient = µ sin( ) as a function of the inertial number I, which leads to the µ I ( )-rheology [26], This shows that the µ I ( )-rheology can be classified as a Drucker-Prager yield surface with a von Mises plastic potential, within the framework of plasticity theory.
The Mohr-Coulomb criterion additionally takes into account the deformation type by incorporating both the major and minor principal stresses, sin( ) 0.
This relation can be expressed in terms of the Lode-angle and other invariants as [38 and the viscosity follows as Note that for = 0, i.e. isochoric shear, Eq. ( 28) reduces to Eq. ( 24).The Matsuoka-Nakai criterion is similar to the Mohr-Coulomb criterion, however with a continuous, smooth yield surface.It is defined as The first case in Eq. ( 30) corresponds to true triaxial extension, the second case to true triaxial compression and the last case to isochoric shear, where = I D ( ) 0

3
. The constant k MN is usually chosen such that Matsuoka-Nakai fits the Mohr-Coulomb criterion for triaxial extension and compression, i.e. for = ± °30 .Alternatively, k MN can be chosen such that it fits Mohr-Coulomb and Drucker-Prager for isochoric shear ( = 0), This is more appropriate for granular flows, as we will show later.Note that all three relations contain only a single constitutive parameter, the friction angle .

Stationary zones
As mentioned before, stresses lie only on the yield surface (in other words, the yield criterion is fulfilled = f T ( ) 0), if the respective mate- rial is flowing.Otherwise (e.g. in stable zones of the slide, levees and deposition) the stress has to lie within the yield surface.However, this is not included in the presented relations and the viscosity will reach a singularity in stationary zones where = D 0. To allow deviatoric stresses lower than the deviatoric stresses at yield, we can limit the viscosity to a predefined upper threshold.This way we can simulate quasi-stationary zones -zones with very high viscosity, that will deform very slowly.The viscosity threshold should be chosen high enough such that the deformation is negligible, but low enough to circumvent numerical problems.However, we want to emphasize that this approach is an extreme simplification of quasi-static granular material behaviour and results (e.g.stresses) in static zones are questionable.An elastoplastic model is more appropriate than the simple visco-plastic model applied in here.The time step in numerical simulations is chosen following the stability criterion (CFL condition) as given by Moukalled et al. [47, pp. 497].Note that the viscosity has to be taken into account in the calculation of the CFL number and that without a limitation of the viscosity, the time step would be very small.This method has been found to be sufficient for the here presented cases.

Element test
We want to consider three strain rate tensors, representing triaxial compression, extension and isochoric shear as shown in Fig. 1 and Table 1.The viscosities and stress tensors can be calculated explicitly in terms of the known relations ( 24), ( 28) and (30).The only parameter is the friction angle , which was set to °36.5 .Table 1 shows the respective stresses as absolute and relative values in relation to the Drucker-Prager yield surface.All models predict the same stresses for isochoric shear, since they have been calibrated for this case.In comparison to Drucker-Prager, Matsuoka-Nakai and Mohr-Coulomb predict higher deviatoric stresses for triaxial compression and lower deviatoric stresses for triaxial extension, which can be seen in Fig. 3 as well.The biggest difference, of 44%, can be observed between Drucker-Prager and Mohr-Coulomb for triaxial compression.To yield the high stresses of the Mohr-Coulomb model with Drucker-Prager, a friction angle of °59 would be required, highlighting the considerable discrepancy be- tween models.

Numerical experiments
For a practical application, the yielding criteria and the respective viscosity models have been implemented into the CFD-toolkit OpenFOAM.OpenFOAM provides a convenient interface for the implementation of custom viscosity models, which can then be used within various solvers [48].We chose the solver multiphaseInterFoam which allows simulation of complex geometries by tracking multiple phases with phase indicator functions.The first phase is representing dry sand and applies the granular viscosity model and a density of ) it has neglectable influence on the granular phase.Tracking of the motion of the interface is realized by transporting phase indicator functions s (sand) and a (air) with the velocity u, which is shared among all phases,  = .
i i i (37) This results in a simple model for granular flows with complex geometries and large strain.We apply a friction angle of = °36.5 in all simulations, similar to Savage et al. [49].The viscosity is truncated to the interval = [10 , 10 ] m s s 5 0 2 1 .The simplest and most common benchmarks for granular flow models are two-dimensional column collapses, flows on inclined planes and flows down chutes.However, they all lead to isochoric plane strain conditions and are therefore inappropriate for our investigations.As shown with the simple element test, all models will yield the same results in two-dimensional simulations.This behaviour has been confirmed with back-calculations of the experiments by Balmforth and Kerswell [50].Results are not shown here as they basically match previous results of e.g.Lagrée et al. [51] and Savage et al. [49].It follows that for a meaningful comparison of the proposed rheologies, we need to induce three-dimensional deformations.Axisymmetric simulations are an obvious choice for such a task, as they enforce threedimensional deformations while keeping computational expense comparable to a two-dimensional case.

Cylindrical granular collapse
The simplest axisymmetric experiment is the collapse of a granular cylinder (Fig. 5), which is presented in the following.Physical experiments of such have been conducted by Lube et al. [36].We choose to simulate two cases with different aspect ratios, .By exploiting the radial symmetry of the problem, we could apply a two-dimensional wedge-shaped mesh [35] with size × 0.5 m 0.2 m.The cell size was set to 1.67 mm and a mesh refinement study (cell sizes 2.5 mm, 1.67 mm, 1.25 mm) indicated that this resolution is sufficient for all practical purposes, with relative errors of less than 1%.The simulation duration was set to 0.5 s, which was sufficient for the duration of the collapse of about 0.35 s.The strict stability criterion following Moukalled et al. [47] and the high viscosity leads to very small time steps of approximately 10 s 6 , which is more restrictive than the traditional CFL type criterion.
Vertical sections of the granular pile three different times and for all rheologies are shown in Fig. 6.In addition, the experimental result of Lube et al. [36] for an aspect ratio of 0.54 is shown, scaled to match the finial pile height h and the final radius r .Fig. 7 highlights the Lode-angle and the strain rate D , as resulting in the simulation with Matsuoka-Nakai yield surface at = t 0.2 s.These fields are widely si- milar in all three simulations and thus not shown repetitively.The differences in stresses between Mohr-Coulomb, Matsuoka-Nakai and Drucker-Prager are highlighted in Fig. Cylindrical granular collapse: Only a small wedge with one cell across the wedge thickness, as highlighted in the figure, is simulated.The geometry is described by the aspect ratio r h / 0 0 .Fig. 6.Cylindrical granular collapse, aspect ratio 1/2, with various yield surfaces at = t 0.1 s (dotted line), = t 0.2 s (dot-dashed line) and = t 0.5 s (dashed line).Grid resolution is 1.67 mm.The colour marks the yield criteria and the black crosses mark the experimental final pile shape [36].Results are very similar, overlap almost entirely and differ by not more than 1%.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) .The same meshing strategy as before has been applied, however, with a computational domain of size × 0.4 m 0.3 m.Vertical sections of the granular pile are shown in Fig. 9, for three different times and all rheologies, together with the experiment of Lube et al. [36] for an aspect ratio of 1.81.The Lode-angle and the strain rate D from the simulation with Matsuoka-Nakai are shown in Fig. 10 for = t 0.1 s.The differences in stresses are highlighted in Fig. 11.Their maximum is similar as in the previous case.

Ring granular collapse
Cylinder collapses show solely negative Lode-angles = ° [ 30 , 0].The full range of the Lode-angle can be shown by the collapse of a ring, as shown in Fig. 12. Similarly to a cylinder, a ring allows to take advantage of the axisymmetry and only a vertical section has to be simulated.We choose an aspect ratio of = h r / 1 0 0 and an inner radius = r h 2.5 i 0 .Taller rings or rings with smaller inner radii collide at the centre, larger inner radii or smaller rings lead to lower Lode-angles and less visible effect.Indeed, for r i one approaches the two-dimensional granular collapse.Vertical sections of the ring at three time steps are shown in Fig. 13.Fig. 14 highlights the Lode-angle and the strain rate D .The differences in stresses are highlighted in Fig. 15.

Discussion and Conclusion
All conducted simulations show stable and smooth results.Mesh depended instabilities, as predicted by Schaeffer [14] and observed by others [52,53] do not show in any of our simulations.We conclude that the applied meshes are too coarse to resolve the small scale features that form the growing instabilities.We expect to see these instabilities in further grid refinements and in such a case the constant friction coefficient = µ sin( ) in Eqs. ( 24), (28), and (30) should be replaced with the regularised relation of Barker et al. [28].The modification of the yield criteria has no influence on the ill-posedness in two dimensions and isochoric plane strain conditions.Thus, all relations will be at least partially ill-posed in three dimensional cases and the presented yield surfaces will be no replacement for the regularisation of e.g.Barker et al. [28].We choose to keep the constant friction coefficient in favour simple equations and parameters.
The presented cases share some characteristics, which allows some conclusions on the general behaviour of investigated rheologies.Granular piles are collapsing while most of the deformation is located in a shear-band that encloses an angle of °35 -°40 with the horizontal plane, roughly matching the friction angle of the material.This corresponds to the failure mechanism as reported by Lube et al. [36].In fact, experiments of Lube et al. [36] correspond reasonably well to the numerical simulations, as shown in Figs. 6 and 9. Material below the shear band is mostly static, leading to a maximum pile inclination similar to the friction angle, roughly matching theoretical predictions [39].However, most regions are flatter due to the inertia of fast moving material.
Stresses in compressive zones differ substantially between the Mohr-Coulomb criterion and the Drucker-Prager criterion (up to 44%).Differences between the Matsuoka-Nakai and the Drucker-Prager criterion are smaller (up to 28%) but still relevant.However, the runout and the final shapes of the granular piles are little affected by the high variation of stresses.The highest variation in the pile shape is visible in the simulation of the high column with a difference of 5%.In other cases, the difference stays well below 1%.
At first glance, it seems unreasonable that a considerable increase of internal stresses does not affect the kinematics significantly.However, a closer investigation reveals that shear-bands overlap with regions where the Lode-angle is close to °0 .This is at least the case for the low collapses (see Figs. 16a and 16b).A Lode-angle close to zero means that stresses have to be similar, as yield surfaces intersect at this point (see Fig. 3).Regions of large strain rates and a non-zero Lode-angle are basically limited to the tall cylinder (see Fig. 16c), which also shows the highest differences in the kinematics.It is reasonable to assume that these shear bands control the kinematics and runout of the collapse.The dominance of these zones on the kinematic behaviour can be further investigated by plotting a histogram of the dissipated energy in terms of the Lode-angle, Fig. 17.All cases show a pronounced peak of dissipated energy at a Lode-angle of °0 .In fact, in the two low cases, almost all the energy is dissipated in the interval °° [ 5 , 5 ], leading to the similarities in kinematics.The granular cylinder with high aspect ratio differs from this behaviour, as a considerable amount of energy is  Finally we can estimate the effect of complex yielding criteria in real case landslides and avalanches.This can be done by introducing the non-dimensional variables, that are usually applied for landslide and avalanche models [54,1,19], see Fig. 18.These scaling laws are based on the assumption that the typical height of the slide H is much smaller than the typical length L, The non-dimensional coordinates and velocities (marked with a hat) follow as and furthermore the second and third invariants as (42) Introducing the dimensionless variables into the Lode-angle yields showing that is small since it contains and otherwise only dimensionless variables.For small angles we can make the approximation sin( ) and get This means that in shallow granular flows the Lode-angle is close to zero.The cylinder with an aspect ratio of 1/2 and the ring are ava- lanching in a relatively thin layer on top (see Figs. 7 and 14) and the shallowness assumption holds.The high cylinder on the other hand collapses completely and the shallowness assumption cannot be applied.Therefore one has to expect a Lode-angle far from zero and thus differences between yield criteria.

Summary and outlook
Using the approach of Schaeffer [14], almost arbitrary yield surfaces can be expressed as non-Newtonian viscosities and thus implemented into the incompressible Navier-Stokes Equations.We showed this by implementing three yield surfaces into the open source toolkit Open-FOAM, namely Drucker-Prager, Mohr-Coulomb and Matsuoka-Nakai.All three yield surfaces have been calibrated for isochoric shear.This means that they are equal for two-dimensional and similar for shear dominated flows.For other deformation types differences of up to 44% in internal stresses have to be expected.Numerical simulations of axisymmetric granular collapses revealed such deformation types and the respective differences.Deformation in mobilised zones was mainly characterized by shearing, in contrast to static zones, where triaxial compression was dominant.However, static zones are irrelevant for the kinematics, leading to a good agreement between yield surfaces in terms of runout.The only exception was the tall granular collapse with aspect ratio 2, where triaxial dominated the early stage of the rapid collapse.
A scaling analysis based on the shallowness assumption, H L / 1 and a Bagnold velocity profile, reveals that deformation in typical landslides and avalanches will be dominated by shearing, indicated by small Lode-angles 1 rad.However, in cases where the shallowness assumption does not hold, predictions of dynamics and kinematics will be different for the presented models.This is especially the case when dealing with obstacles, as the shear dominated flow pattern will be disturbed.Zones of triaxial compression and extension will emerge and deviatoric stresses will vary, leading to further flow pattern changes and possibly highly varying forces on the obstacles.There is strong evidence that Mohr-Coulomb and Matsuoka-Nakai might be the better choice here [33,31,32,34].However, a final conclusion can only be drawn with experiments and respective simulations.The extension of a Drucker-Prager yield criterion to a Mohr-Coulomb yield criterion is straight-forward and can be done by adding an additional factor, solely depending on the Lode-angle and the friction angle (see Eq. ( 28)).For the calculation of the runout and especially in depth-integrated models  with complex rheologies (e.g.[55]), the simpler Drucker-Prager model should be sufficient.The here presented methodology can be combined with the popular µ I ( )-rheology by introducing the velocity depended friction coefficient into the yield surfaces, µ I sin( ) ( ).Finally we want to note that all statements in this paper are limited to incompressible flows.In compressible flows, even two-dimensional deformations might be characterised with a Lode-angle unequal zero, making the here found similarities invalid.

Fig. 1 .Fig. 2 .
Fig. 1.Deformation types and the respective Lode-angle.The red solid shows the undeformed state, the blue solid shows the deformed state.Three deformation types are compatible with incompressible flows: Triaxial compression, isochoric shear and triaxial extension.True triaxial compression is a mix of triaxial compression and isochoric shearing, true triaxial extension a mix of triaxial extension and isochoric shearing with the respective range of Lodeangles.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (16) into Eq.(29) leads to a viscosity that represents this yield surface,

= 1430 kg m s 3 ,
whereas the second phase represents air.Because of its low viscosity ( = 1.48 •10 m s a 5 2 1 ) and density ( = kg m a 3

Fig. 3 .
Fig. 3.The yield surfaces in the deviatoric plane: Drucker-Prager (blue), Mohr-Coulomb (orange) and Matsuoka-Nakai (green).The angle enclosed with the horizontal axis is called Lode-angle and its values indicate the type of deformation: triaxial extension (TXE), triaxial compression (TXC) and isochoric shear (SHR).True triaxial extension (TTXE) and true triaxial compression (TTXC) are mixes of respective limiting cases and isochoric shear.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

.
The simulations with an aspect ratio of 1/2 have been realized with a cylinder of height =

8 ..2 m 0 Fig. 4 .
Fig. 4. Yield surfaces in three-dimensional principal stress space: von Mises (a), Drucker-Prager (b), Mohr-Coulomb (c) and Matsuoka-Nakai (d).The colour marks the Lodeangle .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Cylindrical granular collapse with the Matsuoka-Nakai yield criterion at = t 0.2 s.Grid resolution is 1.67 mm.The black line marks the free surface of the granular pile, the colour displays the Lode-angle (top) and the strain rate D (bottom).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .Fig. 9 .Fig. 10 .
Fig. 8. Cylindrical granular collapse with the Matsuoka-Nakai (top) and Mohr-Coulomb (bottom) yield criterion.The colour displays the ratio between the respective criterion and Drucker-Prager.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 11 .Fig. 12 .Fig. 13 .Fig. 14 .
Fig. 11.Cylindrical granular collapse with the Matsuoka-Nakai (top) and Mohr-Coulomb (bottom) yield criterion.The colour displays the ratio between the respective criterion and Drucker-Prager.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 15 .
Fig. 15.Ring granular collapse with the Matsuoka-Nakai (top) and Mohr-Coulomb (bottom) yield criterion.The colour displays the ratio between the respective criterion and Drucker-Prager.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) , are the slope-parallel coordinates (u v , the respective velo- cities) and z the slope-normal coordinate (w the respective velocity).Assuming continuous shearing along the whole flow depth (i.e.Bagnold-profile), the non-dimensional strain-rate tensor follows as =

Fig. 16 .Fig. 17 .
Fig. 16.Shearbands with a shear rate higher than 1 s 1 (blue) and a Lode-angle between °5 and °5 (red).Cylindrical collapse with = h r / 1/2 0 0 at = t 0.2 s (a), ring collapse at = t 0.2 s (b), cylindrical collapse with = h r / 2 0 0 at = t 0.2 s (c).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 18 .
Fig. 18.Typical height H and length L of a landslide as used in the scaling analysis of Savage and Hutter [54].

Table 1
Element test.
x D D / y D D / z D p T /