Non-isochoric stable granular models taking into account fluidisation by pore gas pressure

Abstract

In this paper, we study non-isochoric models for mixtures of solid particles, at high volume concentration, and a gas. One of the motivations for this work concerns geophysics and more particularly the pyroclastic density currents, whose dense basal parts are precisely mixtures of pyroclasts and lithic fragments and air. They are extremely destructive phenomena, capable of devastating urbanised areas, and are known to propagate over long distances, even over almost flat topography. Fluidisation of these dense granular flows by pore gas pressure is one response that could explain this behaviour and must therefore be taken into account in the models. Starting from a solid–gas mixing model and invoking the compressibility of the gas through a law of state, we rewrite the conservation of mass equation of the gas phase into an equation of the pore gas pressure whose net effect is to reduce the friction between the particles. The momentum equation of the solid phase is completed by generic constitutive laws, specified as in Schaeffer et al. (J. Fluid Mech., vol. 874, 2019, pp. 926–951) by a yield function and a dilatancy function. Therefore, the divergence of the velocity field, which reflects the ability of the granular flow to expand or compress, depends on the volume fraction, pressure, strain rate and inertial number. In addition, we require the dilatancy function to describe the rate of volume change of the granular material near an isochoric equilibrium state, i.e. at constant volume. This property ensures that the volume fraction, which is the solution to the conservation of mass equation, is positive and finite at all times. We also require that the non-isochoric fluidised model be linearly stable and dissipate energy (over time). In this theoretical framework, we derive the dilatancy models corresponding to classical rheologies such as Drucker–Prager and $\mu (I)$ (with or without expansion effects). The main result of this work is to show that it is possible to obtain non-isochoric and fluidised granular models satisfying all the properties necessary to correctly account for the physics of granular flows and which are well posed, at least linearly stable.

1. Introduction

Dense granular flows are present in many research areas, e.g. chemistry, geophysics, biology, engineering and mathematics, in industrial applications, e.g. pharmaceutical production processes, the food industry and construction engineering, and in nature. In the latter case, they constitute a major source of potential danger to human life, buildings and infrastructure in inhabited areas, for example in the case of landslides, which may be caused by seismic activity, soil instability or volcanic eruptions. Landslides, which are a potential source of tsunamis, can be subaerial or submarine and are therefore characterised at first glance by a granular medium submerged in air or water. Pyroclastic density currents (PDCs), as described in Druitt (2007) for instance, whose dense basal parts are mixtures of solid particles (pyroclasts and lithic fragments) and air, are one of the most significant hazards of volcanic eruptions. The basal flows of PDCs behave like a fluid and can travel long distances from the source of the eruption, sometimes more than 100 km. Although understanding the mechanisms responsible for this particular behaviour of concentrated PDCs is one of the main scientific questions related to volcanic processes, modelling the interactions between the dense granular medium and the interstitial gas (air) remains a real challenge.

Because of the large amount of material involved in volcanic processes, the small size of the particles (between $10$ and $100$ microns) in laboratory experiments and the relatively high volume concentration, of the order of $40$ % for moderately expanded flow to $60$ % for dense flow, it seems preferable to consider the granular material and the fluid containing it as a continuous medium. In their pioneering work on solid (particle) and fluid mixtures, Anderson & Jackson (1967) derived a system of coupled equations based on the principles of conservation of mass and momentum. A drag force reflects the solid–fluid interaction in each phase. In the particular case of a gas, by using among other things the fact that the viscosity and density of the gas phase are low, the system can be simplified, resulting in a Darcy-type law relating the velocities of the two phases. An equation relating the density and pressure is then deduced from the mass conservation equation of the gas phase. For compressible gases, the system is closed with a state law, which allows the density to be eliminated as an unknown, and an equation satisfied by the pressure of the interstitial gas is obtained. The resulting fluidised gas–particle mixing model is then formed by the conservation of mass and momentum equations for the solid phase, supplemented by the pore gas pressure equation. The latter intervenes through its gradient in the conservation of momentum equation. The net effect of the fluidisation is therefore to decrease the solid pressure. We thus obtain a system of coupled equations which models the granular flow while taking into account the interstitial gas through its pressure effects. However, a rheology for the solid phase still needs to be specified in order to have a complete and closed system. Two main approaches, which provide constitutive laws for continuous granular models, have been developed over the last two decades. One is based on kinetic theory, initiated by Garzó & Dufty (1999) and extended by Berzi and Jenkins among others – see Berzi, Jenkins & Richard (2020) for recent advances – while the other is based on phenomenological relationships between characteristic quantities of granular flows measured in laboratory experiments. The latter is studied in this article.

A granular medium flows only if the stress exceeds a threshold, the yield stress, otherwise it does not deform and behaves like a solid. The most advanced model, and certainly the one most used for simulating granular flows as a continuum and that accounts for this peculiar behaviour, is the $\mu (I)$-rheology proposed by GDR MiDi (2004), Da Cruz et al. (2005) and Jop, Forterre & Pouliquen (2006). The underlying yield criterion of the $\mu (I)$-model is of the Drucker–Prager type, i.e. the internal friction coefficient of the granular material is proportional to the solid pressure. Recall that the pore gas pressure acts in the granular momentum conservation equation by decreasing the solid pressure, so the friction between particles is reduced and the granular flow is fluidised, allowing it to travel greater distances. The question of a possible fluidisation effect of the pore gas pressure has been addressed through laboratory experiments on the collapse of dense granular columns fluidised by air injection from below in Roche (2012). Indeed, it has been shown in Roche (2012) that fluidised columns flow over distances twice as long as non-fluidised columns.

Since the work of Jop et al. (2006), the $\mu (I)$-rheology has been at the centre of much research both for its contribution to the production of numerical results (let us quote for example Lagrée, Staron & Popinet 2011; Ionescu et al. 2015; Martin et al. 2017) and for theoretical questions related to the well posedness of granular models based on this constitutive law (see Barker et al. 2015, 2017; Barker & Gray 2017; Goddard & Lee 2017; Schaeffer et al. 2019). For dense granular flows with a volume concentration $\phi$ close to the packing limit, of the order of $60$ % per unit volume, the variation in $\phi$ can reasonably be neglected. In this case, the mass conservation equation implies that the granular flow is incompressible (isochoric), neglecting to model the packing and dilatant effects of the granular material when sheared. This approach, because of its apparent simplicity, is attractive and has been used to conduct numerical simulations for granular column collapse (see for instance Lagrée et al. 2011; Ionescu et al. 2015). Whilst successful results have been obtained for predicting the profile of the granular mass during its collapse as well as for estimating the velocity of the front, the incompressible $\mu (I)$-model is ill posed, in the sense of Hadamard, as demonstrated by Barker et al. (2015). Indeed, at low and high inertial numbers, small perturbations grow at an exponential rate in the high-wavenumber limit. The numerical solutions, when they do not blow up, depend on the grid with bands of high gradients appearing in the strain rate and pressure, as shown in Martin et al. (2017). Note that similar instabilities have been observed with a viscous Drucker–Prager model (see for example Martin et al. 2017; Chupin et al. 2021). A simple way to circumvent the ill posedness of the incompressible $\mu (I)$-model is to regularise the constitutive law as proposed in Barker & Gray (2017) and implemented in Gesenhues et al. (2019).

Taking into account the expansion of the granular material, when it is sheared, allows us to regularise the $\mu (I)$-rheology and to obtain a linearly stable model (as shown in Barker et al. 2017; Heyman et al. 2017; Schaeffer et al. 2019). In Barker et al. (2017), the authors introduced a yield function and a dilatancy function, both of which depend on the volume fraction, the (solid) pressure and the inertial number. The divergence of the velocity field is assumed to be proportional to the dilatancy function and the strain rate, which makes local variations of the volume of the granular medium possible. Also, the yield condition specifies the deviatoric stress in terms of a yield function. Note that, unlike the $\mu (I){-}\varPhi (I)$ rheology which provides a state law defining the pressure in terms of the volume fraction, here the volume fraction evolves with the flow according to the conservation of mass equation, and the rate of volume change is specified according to the yield function. As shown in Barker et al. (2017), if the yield and dilatancy functions satisfy three conditions, namely one equation and two inequalities, then the resulting compressible $\mu (I)$-model is linearly well posed.

Without imposing any further constraints, it is possible to find yield and dilatancy functions that satisfy these stability conditions. However, they are in this case obtained by purely mathematical arguments and do not take into account the physics of granular flows. On the other hand, Pailha & Pouliquen (2009) combined the $\mu (I)$-rheology with a dilatancy model based on the critical state theory proposed by Roux & Radjaï (1998). This theory introduces a dilatation angle which reflects the need for a granular material to expand, by increasing the volume it occupies, when sheared. Unfortunately, the resulting constitutive laws, rewritten in the terms of the theory developed in Barker et al. (2017), do not satisfy the stability conditions, so the model may not be well posed.

The objective of this work is to propose non-isochoric granular models, i.e. with local variations of the volume, which take into account the effects of fluidisation by the pore gas pressure. The constitutive laws are written in terms of two functions as in the linear stability theory developed by Barker et al. (2017) (see also Schaeffer et al. 2019). We require that these models be linearly stable, be compatible with the dilatancy model of Roux & Radjaï (1998) and dissipate energy over time and that the volume fraction, the solution of the conservation of mass equation, be positive and bounded from above for any time. By compatibility with the Roux and Radjaï dilatancy model, we mean that at equilibrium the divergence of the velocity field must be zero and that the rate of volume change depends on the deviation of the volume fraction from the equilibrium state. It should be noted that the requirement that the energy of the system must be dissipated over time is motivated by two reasons, one numerical and one theoretical. Indeed, if the energy is dissipated, one can hope to prove that the model is well posed, i.e. that a solution exists and can be unique. Moreover, in order to develop stable numerical schemes, i.e. with bounded solutions, it is more than desirable that the continuous model be dissipative. In this framework, we derive the dilatancy models obtained for specific choices of classical rheologies, such as Drucker–Prager and $\mu (I)$. This approach allows us to derive non-isochoric fluidised granular models with the above-mentioned properties.

The paper is organised as follows. In § 2, a solid–gas mixing model derived from Anderson and Jackson's equations is described. The fluidisation of granular flows by compressible gases, in the case of a general state law, is studied in § 3. The special case of a perfect gas (such as air) is also discussed. In § 4, the fluidised model is completed by generic constitutive laws: the yield condition and the dilatancy law are defined by introducing functions similar to those used in Barker et al. (2017). Section 5 is devoted to the study of the main properties of this generic fluidised and non-isochoric granular model, namely energy dissipation, compatibility of the dilatancy law with equilibrium conditions, linear stability and volume fraction bounds. Finally, within this theoretical framework, dilatancy laws corresponding to classical rheologies, such as Drucker–Prager and $\mu (I)$, are derived in § 6. The resulting models of fluidised and non-isochoric granular flows satisfy all the aforementioned properties.

2. Solid–gas two-phase model

2.1. Governing equations

We consider a mixture of solid particles, i.e. a granular phase of (constant) density $\rho _s$, and a gas whose (variable) density is denoted by $\rho _f$. If $\phi$ denotes the local volume fraction of the particles within the mixture, then the mass conservation for both constituents is written as

(2.1)$$\begin{gather} \partial_t (\phi \rho_s) + \mathrm{div} ( \phi \rho_s {\boldsymbol u} ) = 0, \end{gather}$$

(2.2)$$\begin{gather}\partial_t ((1-\phi) \rho_f) + \mathrm{div} ( (1-\phi) \rho_f {\boldsymbol u}_f ) = 0, \end{gather}$$

where ${\boldsymbol u}$ and ${\boldsymbol u}_f$ correspond respectively to the velocities of the granular phase and the gas phase.

The conservation of momentum equations for the two phases involve the forces between the two components. These equations are derived from Jackson's book (Jackson 2000). A detailed explanation of each of the terms involved is also given in Pitman & Le (2005, Appendix A). The resulting system of equations is written as

(2.3)$$\begin{gather} \phi \rho_s ( \partial_t {\boldsymbol u} + {\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol u} ) = \phi \rho_s {\boldsymbol g} - \boldsymbol{\nabla} p + \mathrm{div} \,\boldsymbol{\tau} - \phi \boldsymbol{\nabla} p_f + \boldsymbol F_{{drag}}, \end{gather}$$

(2.4)$$\begin{gather}(1-\phi) \rho_f ( \partial_t {\boldsymbol u}_f + {\boldsymbol u}_f \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol u}_f ) = (1-\phi) \rho_f {\boldsymbol g} - (1-\phi) \boldsymbol{\nabla} p_f - \boldsymbol F_{{drag}}, \end{gather}$$

where $p$ and $p_f$ correspond to the pressures within each of the phases, i.e. the granular phase and the gas phase, respectively. Note that the pore pressure $p_f$ acts equally on both phases while the solid pressure $p$, also called the effective pressure, which represents normal chain forces between particles, is present only in the momentum equation of the granular phase. The vector ${\boldsymbol g}$ corresponds to the gravity force. The tensor $\boldsymbol {\tau }$ expresses the extra stresses associated with the granular phase (whereas it is assumed that the only stresses associated with the gas phase are due to pressure). An explicit expression for $\boldsymbol {\tau }$ as a function of $\phi$, $p$ and $\boldsymbol {\nabla } {\boldsymbol u}$ will be specified later in the section on the rheology of the medium. Finally, the drag force, denoted by $\boldsymbol F_{{drag}}$, acts in the direction of the relative velocity ${\boldsymbol u}_f-{\boldsymbol u}$ and also depends on the particle concentration. The usual models – see for instance Ergun (1952), Wen (1966) and more recently Gidaspow (1994) and Beetstra, van der Hoef & Kuipers (2007) – can be written as

(2.5)\begin{equation} \boldsymbol F_{{drag}} = \tilde\beta(\phi,|{\boldsymbol u}_f-{\boldsymbol u}|)({\boldsymbol u}_f-{\boldsymbol u}). \end{equation}

Using the arguments of Bouchut et al. (2015), for small values of $|{\boldsymbol u}_f-{\boldsymbol u}|$ this force is proportional to the relative velocity, so we can write

(2.6)\begin{equation} \boldsymbol F_{{drag}} = \beta(\phi)({\boldsymbol u}_f-{\boldsymbol u}), \end{equation}

where $\beta$, which depends only on $\phi$, is called the drag coefficient. For instance, the following expressions are obtained:

(2.7)\begin{equation} \beta(\phi) = \frac{150 \eta_f \phi^2}{d^2 (1-\phi)} \end{equation}

with the model in Gidaspow (1994) and

(2.8)\begin{equation} \beta(\phi) = \frac{18 \eta_f (1-\phi)^2\phi}{d^2 (1-\phi)}\left(10\frac{\phi}{(1-\phi)^2}+(1-\phi)^2\left(1+\frac 32\sqrt{\phi}\right)\right) \end{equation}

with the model given in Beetstra et al. (2007). In the above relations, $d$ is the diameter of the grains and $\eta _f$ the viscosity of the gas.

The work of Anderson, Sundaresan & Jackson (1995, p. 331) – see also Pailha & Pouliquen (2009, (3.12)) – compares the different contributions to the fluid momentum conservation. The authors introduce approximations associated with the smallness of density of a typical gas. Following their work, (2.4) reduces to

(2.9)\begin{equation} \beta(\phi) ({\boldsymbol u}_f - {\boldsymbol u}) =- (1-\phi) \boldsymbol{\nabla} p_f.\end{equation}

This equation can be seen as a Darcy law: it allows us to express the fluid velocity ${\boldsymbol u}_f$ in terms of the solid one ${\boldsymbol u}$, the gradient of the fluid pressure $p_f$ and the volume solid fraction $\phi$. Using (2.6) and (2.9), equations (2.1), (2.2) and (2.3) now read

(2.10)$$\begin{gather} \partial_t \phi + \mathrm{div} ( \phi {\boldsymbol u} ) = 0, \end{gather}$$

(2.11)$$\begin{gather}\partial_t ((1-\phi) \rho_f) + \mathrm{div} ( (1-\phi) \rho_f {\boldsymbol u} ) = \mathrm{div} ( \kappa(\phi) \rho_f \boldsymbol{\nabla} p_f), \end{gather}$$

(2.12)$$\begin{gather}\phi \rho_s ( \partial_t {\boldsymbol u} + {\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol u} ) = \phi \rho_s {\boldsymbol g} - \boldsymbol{\nabla} p + \mathrm{div} \,\boldsymbol{\tau} - \boldsymbol{\nabla} p_f. \end{gather}$$

where $\kappa (\phi ) = (1-\phi )^2/\beta (\phi )$. Note that by using the expression (2.7), the well-known Carman–Kozeny relationship is obtained (see Carman 1937 or Carman 1997 and Kozeny 1927) namely $\kappa (\phi ) = d^2 (1-\phi )^3/(150\eta _f \phi ^2)$.

2.2. Energy estimate

One of the key points we wish to emphasise in this article is that the proposed model is energetically consistent, i.e. it has an energy that decreases over time (in the absence of external forces, such as gravity forces). The energy estimates associated with this type of flow are generally obtained by taking the scalar product of the conservation of momentum equation (2.12) and the velocity ${\boldsymbol u}$, then integrating with respect to the spatial variable. Here, we deduce that

(2.13) \begin{align} &\int \phi \partial_t \left(\frac{1}{2} \rho_s |{\boldsymbol u}|^2 \right) + \int \phi {\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla}\left(\frac{1}{2}\rho_s |{\boldsymbol u}|^2 \right)\nonumber\\ & \quad = \int \phi \rho_s {\boldsymbol g} \boldsymbol{\cdot} {\boldsymbol u} - \int \boldsymbol{\nabla} p \boldsymbol{\cdot}{\boldsymbol u} + \int (\mathrm{div} \,\boldsymbol{\tau})\boldsymbol{\cdot} {\boldsymbol u} - \int \boldsymbol{\nabla} p_f\boldsymbol{\cdot} {\boldsymbol u}. \end{align}

Note that, throughout this document, integrations with respect to the space variable are denoted by $\int$. They are performed over a domain $\varOmega \subset \mathbb {R}^2$ on which we will assume that there is no exchange with the outside: the velocities and normal stresses are assumed to be zero on the boundary $\partial \varOmega$ so that there will never be any boundary terms due to the various integrations by parts.

Multiplying (2.10) by $\frac {1}{2}\rho _s |{\boldsymbol u}|^2$ and integrating, we also obtain

(2.14)\begin{equation} \int \partial_t \phi \left(\frac{1}{2}\rho_s |{\boldsymbol u}|^2 \right) + \int \mathrm{div}(\phi {\boldsymbol u}) \left(\frac{1}{2}\rho_s |{\boldsymbol u}|^2 \right) = 0. \end{equation}

The sum of the last two equalities (2.13) and (2.14), combined with integrations by parts, gives the following estimate:

(2.15)\begin{equation} \frac{{\rm d}}{{\rm d} t} \left( \frac{1}{2} \int \phi \rho_s |{\boldsymbol u}|^2 \right) = \int \phi \rho_s {\boldsymbol g} \boldsymbol{\cdot} {\boldsymbol u} + \underbrace{\int p_f \,\mathrm{div} \, {\boldsymbol u}}_{A} - \underbrace{\int \boldsymbol{\tau} : \boldsymbol{\nabla} {\boldsymbol u}}_{B} + \underbrace{\int p \,\mathrm{div} \,{\boldsymbol u}}_{C}. \end{equation}

Remark 2.1 Note that it is possible to include other, more realistic boundary conditions in this study. For instance, if a non-zero (given) velocity is imposed on part of the boundary, then additional source terms such as those related to gravity forces will have to be taken into account on the right-hand side of (2.15). These terms can increase the energy of the system. In order to take into account the friction of the granular material with a boundary, we may enforce the following conditions:

(2.16)\begin{equation} {\boldsymbol u}_n = 0 \quad \text{and} \quad \boldsymbol{\tau}_{\!t} =- \tan(\alpha_b) (-\boldsymbol{\tau}_{\!n})^+ \frac{\boldsymbol{u}_t}{|\boldsymbol{u}_t|}, \end{equation}

mixing the normal and tangential components ${\boldsymbol u}_t$ and ${\boldsymbol u}_n$ of the velocity ${\boldsymbol u}$ and the normal and tangential components $\boldsymbol {\tau }_t$ and $\boldsymbol {\tau }_n$ of the normal stress $\boldsymbol {\tau }\boldsymbol{\cdot} {\boldsymbol n}$, with $\boldsymbol {n}$ corresponding to the outward unit normal at the boundary. The coefficient $\alpha _b$ depends on the friction angle of the solid particles with the boundary material. With such boundary conditions, additional terms appear during integration by parts in order to deduce (2.15), but these are dissipative terms, i.e. they tend to reduce the energy of the system. For the sake of completeness, the estimate (2.15) becomes

(2.17)\begin{equation} \frac{{\rm d}}{{\rm d} t} \left( \int \phi \frac{\rho_s |{\boldsymbol u}|^2}{2} \right) + \int_{\partial\varOmega} \tan(\alpha_b) (-\boldsymbol{\tau}_{\!n})^+ |\boldsymbol{u}_t| = \int \phi \rho_s {\boldsymbol g} \boldsymbol{\cdot} {\boldsymbol u} + A - B + C. \end{equation}

So far, the three-equation model (2.10)–(2.12) has six unknowns, namely $\phi$, ${\boldsymbol u}$, $\rho _f$, $p_f$, $p$ and $\boldsymbol {\tau }$. In order to close the system, (2.10)–(2.12) must be completed by three closing relations. The terms $A$, $B$ and $C$ of the energy equation are intrinsically linked to the choice of closure laws. The term $A$ depends on the pore gas pressure $p_f$ and thus reflects the fluidisation of the solid phase by the presence of the gas phase. We will see in the next section how an equation for $p_f$ can be derived from (2.11) by accounting for the gas compressibility through a (generic) state law. A second closure relation is given by the constitutive law specifying the deviatoric stress tensor $\boldsymbol {\tau }$ as a function of $\phi$, $\boldsymbol {\nabla } {\boldsymbol u}$ and $p$. The rheology must be physically consistent but must also allow control of the term $B$ in (2.15). Finally, the last closure relation will specify the divergence of the velocity field ${\boldsymbol u}$, which expresses the local change in volume and thus governs the expansion or compression of the flow when it is sheared. Note that control of the terms $A$ and $C$ depends on $\mathrm {div} \, {\boldsymbol u}$.

3. Gas compressibility and fluidisation model

The fluidisation phenomenon of the granular flow is mainly due to the fact that the gas trapped between the particles is compressible. In other words, its density $\rho _f$ is not constant. Depending on the nature of the gas considered, a state law relating the pressure $p_f$ and the density can be imposed, namely

(3.1)\begin{equation} p_f = Q(\rho_f), \end{equation}

where $Q$ is a differentiable function.

In order to obtain an energy estimate and to compensate for the term $A$ in (2.15), we take inspiration from the methods used for the study of compressible fluids; see for instance Lions (1998). For a general state law of the form (3.1), (2.11) becomes

(3.2)\begin{equation} \partial_t ((1-\phi) {\rho}_f) + \mathrm{div} ( (1-\phi) {\rho}_f {\boldsymbol u} ) = \mathrm{div} ( \kappa(\phi) {\rho}_f {Q}'({\rho}_f) \boldsymbol{\nabla} {\rho}_f ). \end{equation}

First, we write (3.2) in non-conservative form

(3.3)\begin{equation} (1-\phi) ( \partial_t {\rho}_f + {\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\rho}_f ) + {\rho}_f \,\mathrm{div} \, {\boldsymbol u} = \mathrm{div} ( \kappa(\phi) {\rho}_f {Q}'({\rho}_f) \boldsymbol{\nabla} {\rho}_f ). \end{equation}

By multiplying the last equation by ${H}'({\rho }_f)$ where ${H}:\mathbb {R} \rightarrow \mathbb {R}$ is any smooth function, we obtain

(3.4)\begin{equation} (1-\phi) ( \partial_t {H}({\rho}_f) + {\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla} {H}({\rho}_f) ) + {\rho}_f {H}'({\rho}_f)\, \mathrm{div} \, {\boldsymbol u} = {H}'({\rho}_f)\, \mathrm{div} ( \kappa(\phi) {\rho}_f {Q}'({\rho}_f) \boldsymbol{\nabla} {\rho}_f ), \end{equation}

which can be rewritten in conservative form as

(3.5)\begin{align} &\partial_t ((1-\phi){H}({\rho}_f)) + \mathrm{div}( (1-\phi){H}({\rho}_f) {\boldsymbol u}) \nonumber\\ &\quad + ({\rho}_f {H}'({\rho}_f) - {H}({\rho}_f))\, \mathrm{div} \, {\boldsymbol u} = {H}'({\rho}_f)\, \mathrm{div} ( \kappa(\phi) {\rho}_f {Q}'({\rho}_f) \boldsymbol{\nabla} {\rho}_f ). \end{align}

Noting that the term $A$ that we wish to control is $A = \int Q(\rho _f)\, \mathrm {div} \, {\boldsymbol u}$, it suffices to choose ${H}$ such that $x{H}'(x)-{H}(x) = {Q}(x)$ and to integrate with respect to the space variable, which leads to

(3.6)\begin{equation} -\!A = \frac{{\rm d}}{{\rm d} t} \left( \int (1-\phi) {H}({\rho}_f) \right) - \int {H}'({\rho}_f)\, \mathrm{div} ( \kappa(\phi) {\rho}_f {Q}'({\rho}_f) \boldsymbol{\nabla} {\rho}_f ). \end{equation}

Since $x{H}''(x)={Q}'(x)$, an integration by parts allows us to rewrite the last term as

(3.7)\begin{equation} -\!A = \frac{{\rm d}}{{\rm d} t} \left( \int (1-\phi)H(\rho_f) \right) + \int \kappa(\phi) |\boldsymbol{\nabla} p_f|^2. \end{equation}

Remark 3.1 As previously, in all integrations by parts, the boundary terms are zero. In the latter case, this cancellation comes from the assumption of zero normal velocities at the boundary and from Darcy's law (2.9). More precisely, on the boundary $\partial \varOmega$, we have

(3.8)\begin{equation} \kappa(\phi) {\rho}_f {Q}'({\rho}_f) \boldsymbol{\nabla} {\rho}_f\boldsymbol{\cdot} \boldsymbol{n} =-(1-\phi) {\rho}_f ({\boldsymbol u}_f-{\boldsymbol u})\boldsymbol{\cdot} \boldsymbol{n} = 0. \end{equation}

In practice, to determine the function $H$ involved in the energy estimate from the function $Q$ specifying the state law, we integrate the differential equation $x{H}'(x)-{H}(x) = {Q}(x)$. Indeed, we find

(3.9)\begin{equation} H(x) = x\int_{c_1}^x \frac{Q(\zeta)}{\zeta^2}\, \mathrm d \zeta + c_2. \end{equation}

The reader is referred to Lions (1998, p. 36), where similar calculations are conducted. It can be noted here that the choice of integration constants $c_1$ and $c_2$ has no influence on the energy. Indeed, if we add a constant $c_2$ to $H$ then the quantity $A$ defined above by the equality (3.7) will be increased by $c_2 ({{\rm d}}/{{\rm d} t}) \int (1-\phi )$, which is zero owing to the mass conservation equation (2.10). In the same way, if we add a linear term $c_1 x$ to $H$ then $A$ will be increased by $c_1 ({{\rm d}}/{{\rm d} t})( \int (1-\phi )\rho _f )$, which is zero owing to (3.2) and the remark above.

In the particular case of an ideal gas, the use of an affine relationship between pressure and density makes it possible, as before, to derive (see supplementary material S.1 available at https://doi.org/10.1017/jfm.2023.1010) from (2.11) a ‘diffusion’ equation describing the evolution of $p_f$, namely

(3.10)\begin{equation} \partial_t ((1-\phi) p_f) + \mathrm{div} ( (1-\phi) p_f {\boldsymbol u} ) + p_{{atm}} \,\mathrm{div} \, {\boldsymbol u} = p_{{atm}} \,\mathrm{div} ( \kappa(\phi) \boldsymbol{\nabla} p_f ), \end{equation}

where $p_{{atm}}=1.013 \times 10^{5} \ \mathrm {Pa}$ is the atmospheric pressure. In the remainder of this article, this equation will be used to describe the evolution of the pore gas pressure, although a more general model can be chosen (i.e. for a general gas). Note that a similar equation (see supplementary material S.1) is used in Goren et al. (2010, (7)) and in McNamara, Flekkøy & Måløy (2000, (7)), corrected in Anghel et al. (2006).

4. Generic model for rheology and dilatation

The rheology of granular media is complex, and most classical models are known to be ill posed when granular flow is assumed to be incompressible; see Barker et al. (2015) and Martin et al. (2017). In Barker et al. (2017) – see also Schaeffer et al. (2019) and Barker et al. (2023) – the authors have shown that taking into account the dilation of the granular medium allows us, in the two-dimensional case, to regularise these models and to remedy these instabilities. More precisely, the constitutive and the dilatancy laws, prescribing respectively the extra stress and the rate of volume change, must be defined in a concordant way.

In this section, we will use the same notation as in Schaeffer et al. (2019). We define the deviatoric strain-rate tensor by

(4.1)\begin{equation} {\boldsymbol S} = \frac{\boldsymbol{\nabla} {\boldsymbol u} + (\boldsymbol{\nabla} {\boldsymbol u})^\mathrm{T}}{2} - \frac{1}{2}(\mathrm{div} \, {\boldsymbol u}) \, {\boldsymbol Id}, \end{equation}

which is a symmetric and traceless tensor, and its second invariant $|{\boldsymbol S}|$ by $|{\boldsymbol S}|^2 = \frac {1}{2}{\boldsymbol S}:{\boldsymbol S} = \frac {1}{2} \sum _{i,j}S_{ij}S_{ij}$. We also introduce the inertial number (see GDR MiDi 2004; Jop et al. 2006)

(4.2)\begin{equation} I = \frac{{d} |{\boldsymbol S}|}{\sqrt{p/\rho_s}}, \end{equation}

where $d$ is the grain diameter of the granular medium under consideration (note that in Barker et al. (2023), the $\mu (J)$-rheology is used to describe granular material immersed in water; the $\mu (J)$-rheology applies to granular flows with a low Stokes number ($St=\rho _s d^2|{\boldsymbol S}|/\eta _f$), and because the viscosity of air is approximately $50$ times smaller than that of water, the granular flows we consider have a higher Stokes number; the interstitial gas is taken into account by the pressure equation (3.10)).

4.1. Granular rheology

The rheology of a granular flow can be described in a fairly general way by the following relationship (in Barker et al. (2017), Schaeffer et al. (2019) and Barker et al. (2023), $Y(\phi, p, I)$ is used instead of $Z(\phi, I)p$; the latter form seems appropriate, since in all the models studied subsequently, the stress is proportional to the pressure and this also makes it relatively easy to write the stability conditions as discussed below):

(4.3)\begin{equation} \boldsymbol{\tau} = Z(\phi, I) p \frac{{\boldsymbol S}}{|{\boldsymbol S}|}. \end{equation}

To be rigorous, the relation (4.3), having no sense when ${\boldsymbol S}$ vanishes, should be rewritten as

(4.4)\begin{equation} \left.\begin{array}{c@{}} \boldsymbol{\tau} : {\boldsymbol S} = 2\, Z(\phi, I) p |{\boldsymbol S}|, \\ |\boldsymbol{\tau}| \leq Z(\phi, I) p, \quad \text{where $\boldsymbol{\tau}$ is symmetric and traceless} \end{array}\right\}\!. \end{equation}

This corresponds to the usual threshold rheology in granular media: the value of $Z(\phi, I)$ is related to a threshold at which the flow starts to deform. One of the fundamental points is therefore to make $Z(\phi, I)$ explicit. Currently, only empirical measurements allow us to have access to this threshold, and we will see in § 6 several examples of such laws obtained by fitting experimental measurements.

Remark 4.1 It is important to note that the constitutive law (4.3) excludes the case of purely viscous rheologies, of the form $\boldsymbol {\tau } = \eta (\phi,|{\boldsymbol S}|) {\boldsymbol S}$. The addition of such a term to (4.3) does not change the results and will be also discussed in Remark 5.2.

4.2. Dilatancy model

The quantity $\mathrm {div}\,{\boldsymbol u}$ expresses the evolution of the elementary volumes of fluid under the action of a flow moving at the velocity ${\boldsymbol u}$. Imposing a zero-divergence condition is therefore equivalent to imposing that, locally, elementary volumes do not vary. If we want to take into account the expansion and compression of the medium, we need to impose an additional law which specifies how the divergence of the velocity field depends on certain characteristic quantities of the flow. One way of expressing the effects of dilation (see supplementary material S.2) is to impose a relation of the form $\mathrm {div} \, {\boldsymbol u} = 2 |{\boldsymbol S}| \,f(\phi, I)$. By adding this expansion law, the constitutive law (4.3) and the equation for the pressure of the interstitial gas (3.10) to the equations for the conservation of mass and momentum (respectively (2.10) and (2.12)), we have a complete system whose unknowns are the volume fraction $\phi$, the granular velocity ${\boldsymbol u}$, the pressure $p$ and the pore gas pressure $p_f$, namely

(4.5)$$\begin{gather} \partial_t \phi + \mathrm{div} ( \phi {\boldsymbol u} ) = 0, \end{gather}$$

(4.6)$$\begin{gather}\partial_t ((1-\phi) p_f) + \mathrm{div} ( (1-\phi) p_f {\boldsymbol u} ) + p_{{atm}}\, \mathrm{div} \, {\boldsymbol u} = p_{{atm}} \,\mathrm{div} ( \kappa(\phi) \boldsymbol{\nabla} p_f ), \end{gather}$$

(4.7)$$\begin{gather}\phi \rho_s ( \partial_t {\boldsymbol u} + {\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol u} ) = \phi \rho_s {\boldsymbol g} - \boldsymbol{\nabla} p + \mathrm{div} \left( Z(\phi, I) p \frac{{\boldsymbol S}}{|{\boldsymbol S}|} \right) - \boldsymbol{\nabla} p_f, \end{gather}$$

(4.8)$$\begin{gather}\mathrm{div} \, {\boldsymbol u} = 2 |{\boldsymbol S}|\, f(\phi, I). \end{gather}$$

5. Main properties of the generic model

5.1. Energy dissipation

As announced in § 2 (see (2.15)), the quantity $\mathcal {E}_s = \frac {1}{2}\int \phi \rho _s |{\boldsymbol u}|^2$ satisfies the following equation:

(5.1)\begin{equation} \frac{{\rm d}\mathcal{E}_s}{{\rm d} t} = A-B+C + \int \phi \rho_s {\boldsymbol g} \boldsymbol{\cdot} {\boldsymbol u}. \end{equation}

Moreover, we have seen in the previous section (see (3.7)) that the term $A,$ owing to the presence of the interstitial gas, satisfies

(5.2)\begin{equation} -\!A = \frac{{\rm d}\mathcal{E}_f}{{\rm d} t} + \mathcal{D}_f, \end{equation}

where ${\mathcal {E}_f} = \int (1-\phi )H(\rho _f)$ and $\mathcal {D}_f = \int \kappa (\phi ) |\boldsymbol {\nabla } p_f|^2$ so that the energy equation (2.15) becomes

(5.3)\begin{equation} \frac{{\rm d} {\mathcal{E}_s}}{{\rm d} t} +\frac{{\rm d} {\mathcal{E}_f}}{{\rm d} t} + {{\mathcal{D}_f}} + B - C = \int \phi \rho_s {\boldsymbol g} \boldsymbol{\cdot} {\boldsymbol u}. \end{equation}

Let us now examine the effects of the rheology and the dilatancy law on the energy of the system (4.5)–(4.8). Replacing the deviatoric stress tensor $\boldsymbol {\tau }$ and the divergence of the velocity field by the relations (4.3) and (4.8) in the terms $B$ and $C$ appearing in (2.15), we obtain

(5.4)\begin{equation} B-C = \int \boldsymbol{\tau} : \boldsymbol{\nabla} {\boldsymbol u} - p\, \mathrm{div} \, {\boldsymbol u} = 2\int (Z-f)p|{\boldsymbol S}|. \end{equation}

Substituting this relation into (2.15), we deduce that the total energy defined by $\mathcal {E} = {\mathcal {E}_s + \mathcal {E}_f}$ satisfies the following equation:

(5.5)\begin{equation} \frac{{\rm d} \mathcal{E}}{{\rm d} t} + \mathcal{D} = \int \phi \rho_s {\boldsymbol g} \boldsymbol{\cdot} {\boldsymbol u}, \end{equation}

where $\mathcal {D} = \mathcal {D}_f + \mathcal {D}_s$ is a dissipation rate composed of the fluid dissipation $\mathcal {D}_f$ appearing due to the fluid pressure and the solid dissipation $\mathcal {D}_s$ defined by $\mathcal {D}_s = 2\int (Z-f) p |{\boldsymbol S}|$. Indeed, if $Z-f\ge 0$, the total energy of the fluidised granular model decreases over time in the absence of any external force, such as a gravitational force. The positivity of the dissipation rate, which can be obtained by the condition

(5.6)\begin{equation} Z \geq f, \end{equation}

is one of the main properties of the fluidised granular models proposed in the following, and can now be stated.

Theorem 5.1 Under the dissipation condition (5.6) and without external force, the model (4.5)–(4.8) has a decreasing energy.

Remark 5.2 The addition of a viscous term, as suggested in Remark 4.1, does not affect this result. On the contrary, adding a viscous contribution $\eta (\phi,|{\boldsymbol S}|){\boldsymbol S}$ to the stress enhances the dissipation rate by adding the non-negative term $\mathcal {D}_{{dis}} = 2\int \eta (\phi, |{\boldsymbol S}|)|{\boldsymbol S}|^2$ to the left-hand side of the energy equation (5.5).

5.2. Consistency of the dilatation law at equilibrium

According to the arguments of Pailha & Pouliquen (2009) and those of Roux & Radjaï (1998), the rate of volume change, which is proportional to the function $f$, is related to a deviation from an equilibrium steady state of the granular medium, which is isochoric, i.e. with constant volume. It has been experimentally observed that at this equilibrium, the volume fraction $\phi$ depends linearly on the inertial number $I$ such that we have

(5.7)\begin{equation} \phi_{{eq}}(I) = \phi_{{max}} - \Delta\phi I. \end{equation}

(Other formulations similar to (5.7) exist. In Schaeffer et al. (2019), the authors suggest $\phi _{{eq}}^{{sch}}(I) = \phi _{{max}} - \Delta \phi /(1+ 1/I),$ arguing that this law would prevent $\phi$ from becoming negative for large values of $I$. We will show below (see Theorem 5.5) that the $\phi$ solution of (4.5) always remains positive if the equilibrium conditions (5.10) are satisfied, which makes the expression $\phi _{{eq}}^{{sch}}(I)$ useless in the present study.)

The parameters $\phi _{{max}}$ and $\Delta \phi$ are obtained by fitting (5.7) with experimental measurements (e.g. $\phi _{{max}}=0.6$ and $\Delta \phi =0.2$ in Forterre & Pouliquen 2008). Note that discrete element method simulations confirm this result (see Schaeffer et al. (2019) for example). The relation (5.7) can be used to determine the volume fraction at equilibrium by setting $\phi =\phi _{{eq}}(I)$. This empirical law can be transformed (see for instance Heyman et al. 2017) into a model for the effective pressure $p$, leading to the $\mu (I)$–$\phi (I)$ rheology. However, the latter is always ill posed in the two-dimensional case as proved by Heyman et al. (2017) (see also Schaeffer et al. 2019). To overcome this difficulty, a dilatancy function $f$ (see (4.8)) specifying the rate of volume change was introduced into the compressible models of Barker et al. (2017) and then Schaeffer et al. (2019). The empirical law (5.7) is not prescribed in the fluidised granular model (4.5)–(4.8), nor in former models. The approach we propose below is to impose conditions on $f$ so that the dilation law (4.8) is compatible with the physics of dense granular media. Note that as (5.7) is empirical, it is possible to use experimental measurements to refine it and improve its accuracy in certain cases; see Robinson, Holland & Fullard (2023) and Breard et al. (2022). In the present article, we will always use (5.7), but any other reasonable choice could be considered.

In Roux & Radjaï (1998), the authors proposed a dilatancy model inspired by critical state mechanics that relates the rate of volume change to the deviation from equilibrium, i.e.

(5.8)\begin{equation} \mathrm{div} \, {\boldsymbol u}=2a|{\boldsymbol S}|(\phi-\phi_{{eq}}(I)), \end{equation}

where $a$ is a constant. This relation expresses how two layers of beads confined at a given pressure should expand when subjected to a constant shear rate. When a deformation occurs ($|{\boldsymbol S}|>0$), we deduce from (5.8) that the material dilates ($\mathrm {div} \, {\boldsymbol u}>0$) if $\phi >\phi _{{eq}}(I)$ while it contracts ($\mathrm {div} \, {\boldsymbol u}<0$) if $\phi <\phi _{{eq}}(I)$. The granular flow is isochoric ($\mathrm {div} \, {\boldsymbol u}=0$) when there is no deformation (i.e. $|{\boldsymbol S}|=0$) or when $\phi =\phi _{{eq}}(I)$. This behaviour is essential and must be reproduced by the generic dilatancy law (4.8). Note that (5.8) is found for $f(\phi,I)=a(\phi -\phi _{{eq}}(I))$ in (4.8).

From the empirical law (5.7), the inertial number at equilibrium $I_{{eq}}(\phi )$ can be defined by

(5.9)\begin{equation} I_{{eq}}(\phi) = \frac{\phi_{{max}}-\phi}{\Delta\phi}. \end{equation}

For the rest (see § 6), defining the equilibrium conditions in terms of the deviation of the inertial number $I$ from the equilibrium state via $I_{{eq}}(\phi )$, instead of $\phi _{{eq}}(I)$, makes things easier. We now state those conditions.

Theorem 5.3 Under the equilibrium conditions

(5.10)\begin{equation} \left. \begin{array}{c@{}} f(\phi,I)=0 \quad \text{if $I=I_{{eq}}(\phi)$}, \\ f(\phi,I)>0 \quad \text{if $I>I_{{eq}}(\phi)$}, \\ f(\phi,I)<0 \quad \text{if $I< I_{{eq}}(\phi)$}, \end{array}\right\} \end{equation}

the model (4.5)–(4.8) is consistent with the physics presented in Roux & Radjaï (1998).

The conditions (5.10) ensure that the dilatancy function $f$, when it deviates from the isochoric equilibrium position, is consistent with the Roux and Radjaï model. Indeed, at equilibrium, the granular flow is isochoric. If $I>I_{{eq}}(\phi )$, the granular medium expands, whereas it compresses if $I< I_{{eq}}(\phi )$. Figure 1 illustrates this behaviour.

Figure 1. Change in volume of a granular medium near an isochoric equilibrium state characterised by $I=I_{{eq}}(\phi ),$ or equivalently $\phi =\phi _{{eq}}(I)$.

5.3. Linear stability of the model

According to Barker et al. (2017), the granular model consisting of (4.5), (4.7) with $p_f=0$ and (4.8) is linearly stable as soon as the functions $Z$ and $f$ satisfy the following properties:

(5.11)$$\begin{gather} Z - \tfrac{1}{2}I\partial_I Z = f+ I\partial_I \,f, \end{gather}$$

(5.12)$$\begin{gather}Z + I\partial_I Z \geq 0, \end{gather}$$

(5.13)$$\begin{gather}{\partial_I \,f > 0}. \end{gather}$$

We shall see – the proof is detailed in the supplementary material S.3 – that accounting for the presence of the pore gas in the granular medium, as achieved in § 3 by supplementing the mass and momentum conservation equations of the solid particles with (4.6) and adding the term $-\boldsymbol {\nabla } p_f$ to the right-hand side of (4.7), does not change the stability result stated in Barker et al. (2017). More precisely, we have the following result.

Theorem 5.4 Under the conditions (5.11), (5.12) and (5.13), the model (4.5)–(4.8) is linearly stable.

While one might expect the presence of interstitial gas to improve the stability by reducing friction between particles, this result is similar to that obtained by Barker et al. (2017) for the non-fluidised model. The presence of the pore gas pressure in the solid-phase momentum equation tends to lower or even cancel out the effective pressure and so the friction between particles. In this case, we might expect the system to be well posed. On the other hand, when the interstitial pressure diffused by (4.6) decreases, the effective pressure increases and there is no reason why the stability of the fluidised model should be improved.

5.4. Volume fraction bounds

In order for the model to make sense, it is fundamental to ensure that the volume fraction $\phi$ always remains positive and does not exceed the maximum value $\phi _{{max}}$ previously introduced. We now state one of the main results of this paper.

Theorem 5.5 Suppose that the initial condition satisfies $0 \leq \phi |_{t=0} \leq \phi _{{max}}$ and consider any smooth solution to (4.5)–(4.8). We have $\phi \geq 0.$ Moreover, if the equilibrium conditions (5.10) are satisfied, then we have $\phi \leq \phi _{{max}}.$

The technical proof of this result is detailed in Appendix A. Although to show the positiveness of the volume fraction we only need to assume that the velocity field is regular, the equilibrium conditions (5.10) are essential to show that $\phi$ remains bounded above at all times.

6. Physical examples

The objective of this section is to complete the generic fluidised (4.5)–(4.8) model presented in § 4, with constitutive laws specified by the yield function $Z$ and the dilatancy function $f$. The fundamental question is how to determine $Z$ and $f$ in such a way that the physics of granular flows is properly described and the mathematical properties introduced in the previous section are verified, specifically that the energy of the system must be dissipated, the volume fraction must be positive and bounded above, and the model must be linearly stable. According to Theorems 5.1–5.5 these properties hold if the dissipation condition (5.6), the equilibrium conditions (5.10) and the stability conditions (5.11)–(5.13) are satisfied.

Note that the assumptions of Theorem 5.3 require that the sign of the function $f$ correspond to a deviation of the granular flow from an isochoric equilibrium, i.e. of constant volume. The expansion or contraction of the granular material near this equilibrium state is then taken into account in the model. The conditions (5.10) must therefore be fulfilled for both physical and mathematical reasons.

The question of the choice of rheology, i.e. the choice of the function $Z$, must now be addressed. This choice is guided by the state of knowledge on the physics of granular materials. The model classically used is the Drucker–Prager model, which is a multidimensional form of the Mohr–Coulomb law (see for example Forterre & Pouliquen 2008). Later, the $\mu (I)$-rheology emerged thanks to the work of the Groupement de Recherche Milieux Divisés (GDR MiDi 2004) and became established. In practice, these constitutive laws defining the deviatoric stress tensor are obtained by fitting experimental measurements and can therefore be used as a starting point for the construction of a model.

Our methodology is as follows: once the yield function $Z$ is prescribed according to the criteria discussed above, the stability condition (5.11), which can be reduced to a differential equation, is used to determine $f$. At this stage, $f$ is known up to a constant which is independent of the inertial number but may depend on the volume fraction $\phi$. This constant will be fixed so that at equilibrium the divergence of the velocity field equals zero, in agreement with the critical state theory proposed in Roux & Radjaï (1998). The last step will be to make sure that the other conditions (5.6), (5.10), (5.12) and (5.13) are fulfilled. Following this approach, we propose hereafter several fluidised and non-isochoric granular models.

Finally, it should be mentioned that relatively simple models should be preferred as long as a numerical scheme can be designed and implemented. However, this is beyond the scope of this paper and will be the subject of future work.

6.1. Non-isochoric models with classical rheologies

The $\mu (I)$-rheology, which is almost unanimously accepted, provides an expression for the friction coefficient in terms of the inertial number, namely

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

where the function $\mu$ depends on three parameters, $I_0>0$ and $\mu _2>\mu _1>0$. (Following Forterre & Pouliquen (2008, p. 10), typical values of the constants obtained for mono-dispersed glass beads are $\mu _1 = \tan (21^\circ )$, $\mu _2 = \tan (33^\circ )$ and $I_0 = 0.3$. Note that these values are obtained experimentally, and that some authors specify that $I_0$ can depend on $\phi$; see for instance Gray & Edwards (2014, (2.12)) or Jop, Forterre & Pouliquen (2005, Appendix A).) By replacing $Z(\phi,I)$ by $\mu (I)$ in (4.3), we can reformulate the deviatoric stress tensor as in Ionescu et al. (2015), showing that the $\mu (I)$-rheology is a viscoplastic rheology with a Drucker–Prager plasticity criterion, with a yield stress equal to $\mu _1$ and a viscosity that is variable in space and time. It then becomes attractive to replace in numerical simulations the variable viscosity of the $\mu (I)$-rheology, which covers a wide range of values, with a constant viscosity, which amounts to using a viscous Drucker–Prager rheology. This approach has been successfully used to reproduce granular column collapse experiments (see Ionescu et al. 2015; Chupin et al. 2021). For this reason, we also study the latter rheology in what follows. In this case, the function $Z$ is defined by $Z(\phi,I)=\sin (\delta )$ – see Andreotti, Forterre & Pouliquen (2012, table 4.1) – where $\delta$ is the internal angle of friction (note that this rheology is most often written using the tangent function instead of the sine function, but with a different angle; see Andreotti et al. (2012, p. 144) for the details). Before going further in the study of these constitutive laws, we state a general result.

Proposition 6.1 Let $Z$ be a function depending on the volume fraction and on the inertial number. The function $f$ given by

(6.2)\begin{equation} f(\phi,I) = \frac{1}{2} \left(\frac{I_{{eq}}(\phi)}{I} Z(\phi,I_{{eq}}(\phi)) - Z(\phi,I) \right) + \frac{3}{2I}\int_{I_{{eq}}(\phi)}^I Z(\phi,J)\, {\rm d} J \end{equation}

is the solution of (5.11) vanishing at $I=I_{{eq}}(\phi )$.

It can be easily checked that the function $f$ so obtained for the Drucker–Prager model, namely with $Z(\phi,I)=\sin (\delta )$ (see table 1), satisfies the equilibrium condition (5.10) and the stability condition (5.13). Also, the function $Z$ satisfies (5.12). As $Z(\phi,I)-f(\phi,I)=2\sin (\delta )(\phi -\phi _{{max}})/I$ is non-negative, since $\phi \le \phi _{{max}}$ (see Theorem 5.1), the condition (5.6) is satisfied too.

Table 1. Yield and dilatancy functions, i.e. $Z(\phi,I)$ and $f(\phi,I)$, ensuring that the stability condition (5.11) is satisfied.

By using (4.2) and the expression for $f$, the dilatancy law (4.8) can be rewritten as

(6.3)\begin{equation} \mathrm{div} \, {\boldsymbol u} = 2 \sin(\delta) |{\boldsymbol S}| - \frac{2 \sin(\delta)}{d \sqrt{\rho_s}} I_{{eq}}(\phi) \sqrt{p}. \end{equation}

While the Drucker–Prager rheology is known to be linearly unstable for incompressible granular flows, as shown in Martin et al. (2017), a stable model is obtained when this dilatancy law is used. From a physical point of view, this relationship highlights a competition between shear forces, which tend to increase the volume (solid particles separate when the material is sheared), and pressure forces, which compress the flow. Figure 2 illustrates these effects on the granular medium.

Remark 6.2 Using (5.9) and (5.7), the dilatancy law (6.3) can be rewritten as

(6.4)\begin{equation} \mathrm{div}{\boldsymbol u}=2a \vert{\boldsymbol S}\vert( \phi - \phi_{{eq}}(I) ), \end{equation}

with $a={\sin (\delta )}/({\Delta \phi I}).$ With this formulation, the dilatation law obtained for the Drucker–Prager rheology is similar to that proposed in Roux & Radjaï (1998) (see (5.8)).

Figure 2. Volume change under shear (left) and normal (right) stresses.

For the $\mu (I)$-rheology, the function $Z$, defined by $Z(\phi,I)=\mu (I)$, satisfies the second stability condition (5.12). The function $f$ derived from $Z$ using (6.2) is listed in table 1. By noting that

(6.5)\begin{equation} \partial_I f(\phi,I) = \frac{\mu_2-\mu_1}{2I_0} H'\left(\frac{I}{I_0}\right) + \frac{I_{{eq}}(\phi)}{I^2} \left(\mu_2+ \frac{\mu_2-\mu_1}{2} H\left(\frac{I}{I_0}\right)\right) \end{equation}

and observing that $H'(x)>0$ for $x>0$ and $H(0)=-2$, we have $\partial _I f(\phi,I)>0$ so that (5.13) is satisfied. Moreover, we have

(6.6)\begin{align} Z(\phi,I)-f(\phi,I) &= \frac{3(\mu_2-\mu_1)}2\left(\frac{\ln(1+I/I_0)}{I/I_0}-\frac 1{1+I/I_0}\right) \nonumber\\ &\quad +\frac{I_{{eq}}(\phi)}{I}\left(\mu_2+\frac{\mu_2-\mu_1}2 H\left(\frac{I_{{eq}}(\phi)}{I_0}\right)\right)\!, \end{align}

showing that $Z(\phi,I)-f(\phi,I)$ is non-negative. This ensures that the energy is dissipated (see (5.6)). We easily deduce that $f$ also satisfies (5.10). In summary, the fluidised and non-isochoric granular models, based on the Drucker–Prager model and on the $\mu (I)$-rheology, that is, with the functions $Z$ and $f$ specified in the two first lines of table 1, have all the required properties.

6.2. Rheologies with dilatancy effects

In order to take into account the variations of the volume fraction of the granular medium in the rheology, Wood (1990) proposed to introduce a dilatancy angle $\psi$. The latter is a measure of the ratio between the relative vertical and horizontal displacements between two layers of grains when they are sheared. It can be positive (expansion) or negative (contraction). The change in volume (see Andreotti et al. 2012) is expressed in terms of the dilatancy angle by $\mathrm {div} \, {\boldsymbol u} = 2|{\boldsymbol S}| \sin (\psi )$, which can be approximated by

(6.7)\begin{equation} \mathrm{div} \, {\boldsymbol u} = 2|{\boldsymbol S}| \psi,\end{equation}

as long as $\psi$ remains small. By comparing (6.7) and (4.8), we obtain $f(\phi,I)=\psi (\phi,I)$ so that the smallness assumption on $\psi$ is transferred to $f$. Recalling (5.10), we can deduce that $f$ will be small near the equilibrium state, which translates into $I\approx I_{{eq}}(\phi )$ in terms of the inertial number.

Intuitively, we can guess that a close packing that has to dilate in order to deform ($\psi > 0$) has a friction coefficient larger than that of a loose packing that will undergo compaction ($\psi < 0$). The granular viscosity is then an increasing function of $\psi$, cancelling out when $\psi =0.$ In the case of the Drucker–Prager rheology, Andreotti (see Andreotti et al. 2012) proposed the constitutive law $\boldsymbol {\tau } = \sin (\delta +\psi )\, p {\boldsymbol S}/|{\boldsymbol S}|$, which can be approximated by

(6.8)\begin{equation} \boldsymbol{\tau} = ( \sin (\delta) + \cos (\delta) \psi ) p \frac{{\boldsymbol S}}{|{\boldsymbol S}|}, \end{equation}

if $\psi$ remains small. (Curiously, it seems that the development used in Andreotti et al. (2012) (corrected in the next French edition, but repeated in Robinson et al. 2023) is not correct; the coefficient $\cos (\delta )$ had been omitted. For their work this does not change anything, since this constant can be incorporated into the multiplicative coefficient in front of $\psi$.) Similarly, the dilatancy effects can be taken into account in the $\mu (I)$-rheology by modifying the stress $\boldsymbol {\tau }$ (see Robinson et al. 2023) as follows:

(6.9)\begin{equation} \boldsymbol{\tau} = ( \mu(I) + \psi ) p \frac{{\boldsymbol S}}{|{\boldsymbol S}|}. \end{equation}

These laws imply finding functions $Z$ and $f$ such that

(6.10)\begin{equation} Z(\phi,I)=\sin (\delta) + \cos (\delta) f(\phi,I) \end{equation}

for the Drucker–Prager rheology and

(6.11)\begin{equation} Z(\phi,I)=\mu(I) + f(\phi,I) \end{equation}

for the $\mu (I)$-rheology. Using (6.10), the condition (5.11) becomes $(2+\cos (\delta ))I\partial _I f + 2(1-\cos (\delta )) f= 2 \sin (\delta ).$ The solution of this differential equation vanishing at $I_{{eq}}(\phi )$ is provided in the two last lines of table 1. Note that, for common materials such as glass beads or sand, the angle of friction $\delta$ is of the order of $30^\circ$ so that the coefficient $\beta$ (see table 1) is of the order of $0.1$. Now, using (6.11), the condition (5.11) becomes

(6.12)\begin{equation} f(\phi,I)= \frac 23 \int_{I_{{eq}}(\phi)}^I \frac{\mu(J)}J \, {\rm d} J - \frac 13 (\mu(I)-\mu(I_{{eq}})). \end{equation}

By substituting (6.1), integrating and enforcing $f(\phi,I_{{eq}}(\phi ))=0$, we find the expression for $f$ which is written in table 1. It can be easily checked that (5.13) is satisfied for both rheologies.

The condition (5.6) is obviously satisfied for the $\mu (I)$-rheology. For the Drucker–Prager model, we have $Z(\phi,I)-f(\phi,I) = \sin (\delta )\, (I_{{eq}}(\phi )/I)^\beta$ so that $Z-f$ is non-negative, which ensures that the energy is dissipated. We can easily show that, for both rheologies, (5.12) is not true for all values of $I$ but is satisfied as soon as $I\approx I_{{eq}}(\phi )$. Therefore, the fluidised and non-isochoric granular models obtained with the yield and dilatancy functions derived in this section are stable near the equilibrium state.

7. Concluding remarks

In this paper, we have studied non-isochoric fluidised granular models, i.e. models which take into account the fluidisation effects of a compressible interstitial gas and local volume changes. The motivation is to model mixtures of particles, with a high concentration (between $40$ % and $60$ % in volume), and a gas, in particular air. Pyroclastic density flows, which are frequent and can have devastating effects, fit into this framework. They are a major hazard in volcanic eruptions because of the great distances they can travel, up to 100 km in some cases, at high speed, even on gentle slopes. The starting point of this study is the fluid–solid mixing model of Anderson and Jackson, which is simplified here because of the physical characteristics of the interstitial fluid, which is a gas. The compressibility of the latter allows us to transform the mass conservation equation of the gas phase into an equation for the pressure. The effect of the pressure of the interstitial gas is to reduce the friction between the particles and thus make it possible for the flow to accelerate and travel greater distances. This phenomenon has been observed in the laboratory (see Roche 2012), showing that columns of particles (glass beads) travel twice as far when fluidised with an air stream injected from below. In order to close the equations for the solid phase, the constitutive laws, written as in Schaeffer et al. (2019), are specified in terms of a yield function and a dilatancy function, both of which depend on the volume fraction, inertial number and (solid) pressure. The resulting fluidised granular model is non-isochoric as it allows for volume change, i.e. the granular flow can expand when sheared. Moreover, we require that the model be linearly stable, that it dissipate energy (over time), that it be compatible with the dilatancy model proposed by Roux & Radjaï (1998) and that the volume fraction, which is the solution of the mass conservation equation, be positive and bounded at all times. Working within this theoretical framework, we have studied dilatancy laws that are compatible with classical rheologies, i.e. Drucker–Prager and $\mu (I)$, with or without taking into account the effects of dilation in the yield condition. The main objective of this work was to derive fluidised and non-isochoric granular models that satisfy the aforementioned mathematical properties and take into account the known physics of granular materials, in particular in terms of the rheology and deformation of the medium when sheared. The next steps will be to prove that these models are well posed, i.e. that they admit solutions that can be unique, and to propose stable numerical schemes. Then, numerical studies will have to be carried out to prove the validity of the models, in particular by comparing the results of the simulations with laboratory experiments such as those described in Roche (2012). These are the outlines of our future work on the subject.