Laurent Chupin

Thierry Dubois

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

In this paper, we study non-isochoric models for mixtures of solid particles, at high volume concentration, and a gas. One of the motivations of this work concerns geophysics and more particularly the pyroclastic density currents which are precisely mixtures of pyroclast 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 gas-solid 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 on 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 (2019, Journal of Fluid Mechanics, 874, 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 is linearly stable and dissipates energy (over time). In this theoretical framework, we derive the dilatancy models corresponding to classical rheologies such as Drucker-Prager and µ(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 being well-posed, at least linearly stable.

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, food industry, construction engineering, as well as 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 characterized at first glance by a granular medium submerged in air or water. Pyroclastic density currents (PDCs), see [START_REF] Druitt | Pyroclastic density currents[END_REF] for instance, which are mixtures of solid particles (pyroclasts and lithic fragments) and air, are one of the most significant hazards of volcanic eruptions. PDCs are often characterized by the presence of a dense basal flow, which behaves 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. Due to the large amount of material carried out in volcanic processes, the small size of the particles (between 10 and 100 microns) in laboratory experiments and, the relatively high volume concentration, in 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 and Jackson [START_REF] Anderson | Fluid mechanical description of fluidized beds[END_REF] 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 verified by the pressure of the interstitial gas is obtained. The resulting fluidized 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. 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 most used for simulating granular flows as a continuum and that accounts for this peculiar behaviour is the µ(I)-rheology proposed by Jop et al [START_REF] Jop | A constitutive law for dense granular flows[END_REF]. The underlying yield criterion of the µ(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 fluidized, allowing it to travel greater distances. The question of a possible fluidisation effect of the pore gas pressure has been addressed through laboratory experiments of the collapse of dense granular columns fluidized by air injection from below in [START_REF] Roche | Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective[END_REF]. Indeed, it has been shown in [START_REF] Roche | Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective[END_REF] that fluidized columns flow over distances twice as long as non-fluidised columns. Since the work of Jop et al [START_REF] Jop | A constitutive law for dense granular flows[END_REF], the µ(I)-rheology has been at the centre of many researches both for its contribution in the production of numerical results, let us quote for example [START_REF] Lagrée | The granular column collapse as a continuum: Validity of a two-dimensional Navier-Stokes model with a µ(I)-rheology[END_REF][START_REF] Ionescu | Viscoplastic modeling of granular column collapse with pressure-dependent rheology[END_REF][START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I)-rheology and lateral wall effects[END_REF], and for theoretical questions related to the well-posedness of granular models based on this constitutive law [START_REF] Barker | Well-posed and ill-posed behaviour of the µ(I)-rheology for granular flow[END_REF][START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF][START_REF] Barker | Partial regularisation of the incompressible µ(I)-rheology for granular flow[END_REF][START_REF] Goddard | On the stability of the µ(I)-rheology for granular flow[END_REF][START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF]. For dense granular flows with a volume concentration ϕ close to the packing limit, on the order of 60% per unit volume, the variation in ϕ 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 produce numerical simulations for granular column collapse, see for instance [START_REF] Lagrée | The granular column collapse as a continuum: Validity of a two-dimensional Navier-Stokes model with a µ(I)-rheology[END_REF][START_REF] Ionescu | Viscoplastic modeling of granular column collapse with pressure-dependent rheology[END_REF]. 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 µ(I)-model is ill-posed, in the sense of Hadamard, as demonstrated by Barker et al. in [START_REF] Barker | Well-posed and ill-posed behaviour of the µ(I)-rheology for granular flow[END_REF]. 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 on the strain rate and pressure, as shown in [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I)-rheology and lateral wall effects[END_REF]. Note that similar instabilities have been observed with a viscous Drucker-Prager model, see for example [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I)-rheology and lateral wall effects[END_REF][START_REF] Chupin | Pressure-dependent threshold in a granular flow: Numerical modeling and experimental validation[END_REF]. A simple way to circumvent the ill-posedness of the incompressible µ(I)model is to regularise the constitutive law as proposed in [START_REF] Barker | Partial regularisation of the incompressible µ(I)-rheology for granular flow[END_REF] and implemented in [START_REF] Gesenhues | Finite element simulation of complex dense granular flows using a well-posed regularization of the µ(I)-rheology[END_REF]. Taking into account the expansion of the granular material, when it is sheared, allows to regularise the µ(I)rheology and to obtain a linearly stable model, as shown in [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF][START_REF] Heyman | Compressibility regularizes the µ(I)-rheology for dense granular flows[END_REF][START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF]. In [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF], Barker et al 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 µ(I) -Φ(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 [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF], if the yield and dilatancy functions satisfy three conditions, namely one equation and two inequalities, then the resulting compressible µ(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 argumentation and do not take into account the physics of granular flows. On the other hand, Pailha and Pouliquen in [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] have proposed a dilatancy model based on the critical state theory proposed by Roux and Radjaï (see [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF]). 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 [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF], do not satisfy the stability conditions so that 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 [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF][START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF]. We require that these models are linearly stable, that they are compatible with the dilatancy model of Roux and Radjaï [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF], that they dissipate energy over time, and that the volume fraction, solution of the conservation of mass equation, is for any time positive and bounded from above. 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 that it 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 is dissipative. In this framework, we derive the dilatancy models obtained for specific choices of classical rheologies, such as Drucker-Prager and µ(I). This approach allows us to derive non-isochoric fluidised granular models with the above mentioned properties. The paper is organised as follows. In Section 2, a solid-gas mixing model derived from the 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 Section 3. The special case of a perfect gas (such as air) is also discussed. In Section 4, the fluidized model is completed by generic constitutive laws: the yield condition and the dilatancy law are defined by introducing functions similar to those used in [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF]. 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 µ(I), are derived in Section 6. The resulting models of fluidised and non-isochoric granular flows satisfy all the aforementioned properties.

2 Solid-gas two-phase model

Governing equations

We consider a mixture of solid particles, i.e. a granular phase of (constant) density ρ s , and a gas whose (variable) density is noted ρ f . If ϕ denotes the local volume fraction of the particles within the mixture, then the mass conservation for both constituents is written

∂ t (ϕρ s ) + div(ϕρ s u) = 0, (1)
∂ t ((1 -ϕ)ρ f ) + div((1 -ϕ)ρ f u f ) = 0, (2)
where u and u f correspond respectively to the velocity 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 [START_REF] Jackson | The dynamics of fluidized particles[END_REF]. A detailed explanation of each of the terms involved is also given by Pitman and Le in [32, Appendix A]. This system of equations is written as

ϕρ s ∂ t u + u • ∇u = ϕρ s g -∇p + div τ -ϕ∇p f + β(ϕ)(u f -u), (3)
(1 -ϕ)ρ f ∂ t u f + u f • ∇u f = (1 -ϕ)ρ f g -(1 -ϕ)∇p f -β(ϕ)(u f -u), (4)
where p and p f correspond to the pressures within each of the phases, i.e. the granular phase and the gas phase respectively. The tensor τ 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 of τ as a function of ϕ, p and ∇u will be specified later in the section on the rheology of the medium. Finally, β is a drag coefficient which depends on ϕ and which, according to Andreotti et al. [3, p.290], can be written

β(ϕ) = 150 η f ϕ 2 d 2 (1 -ϕ) ,
where d is the diameter of the grain and η f the viscosity of the gas.

Following the work of Anderson et al [1, p.331] (see also [31, equation (3.12)]) which compare the different contributions in the fluid momentum conservation, we could introduced approximations associated with the smallness of density of a typical gas. The equation (4) reduces to

β(ϕ)(u f -u) = -(1 -ϕ)∇p f . (5)
This equation can be seen as a Darcy law: it allows us to express the fluid velocity u f with respect to the solid one u, the gradient of the fluid pressure p f and the volume solid fraction ϕ. The equations ( 1), ( 2) and (3) now read

∂ t ϕ + div(ϕ u) = 0, (6)
∂ t ((1 -ϕ)ρ f ) + div((1 -ϕ)ρ f u) = div(κ(ϕ)ρ f ∇p f ), (7)
ϕρ s ∂ t u + u • ∇u = ϕρ s g -∇p + div τ -∇p f . ( 8
)
where

κ(ϕ) = (1 -ϕ) 2 β(ϕ) = d 2 (1 -ϕ) 3 150 η f ϕ 2 .
Note that this coefficient can be related to the permeability of the material via the well-known Carman-Kozeny relationship, see [START_REF] Carman | Fluid flow through a granular bed[END_REF] or [START_REF] Carman | Fluid flow through granular beds[END_REF], and [START_REF] Kozeny | Uber kapillare leitung des wassers im boden-aufstieg, versickerung und anwendung auf die bewasserung, sitzungsberichte der akademie der wissenschaften wien[END_REF].

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 gravity forces). The energy estimates associated with this type of flow are generally obtained by performing the scalar product of the conservation of momentum equation ( 8) by the velocity u, then integrating with respect to the spatial variable. Here, we deduce

ϕ∂ t ρ s |u| 2 2 + ϕu • ∇ ρ s |u| 2 2 = ϕρ s g • u -∇p • u + (div τ ) • u -∇p f • u. (9)
Note that throughout this document, the integrations with respect to the space variable are noted . They are performed on a domain Ω ⊂ 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 ∂Ω so that there will never be any boundary terms due to the various integrations by parts.

Multiplying then (6) by 1 2 ρ s |u| 2 and integrating, we also obtain

∂ t ϕ ρ s |u| 2 2 + div(ϕu) ρ s |u| 2 2 = 0. ( 10
)
The sum of the last two equality ( 9) and [START_REF] Carman | Fluid flow through a granular bed[END_REF], combined with integrations by parts, gives the following estimate:

d dt ϕ ρ s |u| 2 2 = ϕρ s g • u + p f div u A -τ : ∇u B + p div u C . ( 11
)
So far, the three-equation model ( 6)-( 8) has six unknowns, namely ϕ, u, ρ f , p f , p and τ . In order to close the system, ( 6)-( 8) must be completed by three closing relations. The terms A, B and C of the energy equation are closely related 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 equation ( 7) 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 τ as a function of ϕ, ∇u and p. The rheology must be physically consistent but must also allow the control of the term B in [START_REF] Carman | Fluid flow through granular beds[END_REF]. Finally, the last closure relation will specify the divergence of the velocity field u, which expresses the local change in volume and thus governs the expansion or compression of the flow when it is sheared. Note that the control of the terms A and C depends on div u.

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 ρ f is not constant. Depending on the nature of the gas considered, a state law relating pressure p f and density can be imposed, namely

p f = Q(ρ f ), ( 12
)
where Q is a differentiable function. For instance, for an ideal gas, an affine relationship between pressure and density will be imposed (see Subsection 3.2):

Q(ρ f ) = p atm 1 - ρ f ρ 0 f , (13)
where p atm = 1.013 × 10 5 Pa is the atmospheric pressure, and ρ 0 f = 1 kg.m -3 corresponds to the density of air at atmospheric pressure.

Case of a general gas

In order to obtain an energy estimate and to compensate the term A in the global estimate [START_REF] Carman | Fluid flow through granular beds[END_REF], we take inspiration from the methods used for the study of compressible fluids, see for instance [START_REF] Lions | Compressible Models[END_REF]. For a general state law of the form ( 12), the equation ( 7) becomes

∂ t ((1 -ϕ)ρ f ) + div((1 -ϕ)ρ f u) = div(κ(ϕ)ρ f Q ′ (ρ f )∇ρ f ). ( 14
)
First, we write the equation ( 14) in non-conservative form:

(1 -ϕ) ∂ t ρ f + u • ∇ρ f ) + ρ f div u = div κ(ϕ)ρ f Q ′ (ρ f )∇ρ f .
By multiplying the last equation by H ′ (ρ f ) where H : R → R is any smooth function, we obtain

(1 -ϕ) ∂ t H(ρ f ) + u • ∇H(ρ f )) + ρ f H ′ (ρ f )div u = H ′ (ρ f )div κ(ϕ)ρ f Q ′ (ρ f )∇ρ f ,
which can be rewritten in conservative form as

∂ t ((1 -ϕ)H(ρ f )) + div((1 -ϕ)H(ρ f )u) + (ρ f H ′ (ρ f ) -H(ρ f ))div u = H ′ (ρ f )div κ(ϕ)ρ f Q ′ (ρ f )∇ρ f . (15)
Noting that the term A that we wish to control is

A = Q(ρ f )div u, it suffices to choose H such that xH ′ (x) -H(x) = Q(x)
and to integrate with respect to the space variable, which leads to

-A = d dt (1 -ϕ)H(ρ f ) -H ′ (ρ f )div κ(ϕ)ρ f Q ′ (ρ f )∇ρ f . Since xH ′′ (x) = Q ′ (x)
, an integration by parts allows to rewrite the last term as

-A = d dt (1 -ϕ)H(ρ f ) + κ(ϕ)|∇p f | 2 . ( 16
)
Remark 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 [START_REF] Barker | Well-posedness and ill-posedness of single-phase models for suspensions[END_REF].

More precisely, on the boundary ∂Ω, we have

κ(ϕ)ρ f Q ′ (ρ f )∇ρ f .n = -(1 -ϕ)ρ f (u f -u).n = 0,
where n corresponds to the outward unit normal at the boundary.

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 xH ′ (x) -H(x) = Q(x). Indeed, we find

H(x) = x x c1 Q(ζ) ζ 2 dζ + c 2 . ( 17
)
The reader is referred to [26, 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 appearing in equality ( 16) will be increased by c 2 d dt

(1 -ϕ) which is zero due to the mass conservation equation ( 6). In the same way, if we add a linear term c 1 x to H then the quantity A appearing in equality ( 16) will be increased by c 1 d dt (1 -ϕ)ρ f which is zero due to ( 14) and the remark above.

The particular case of an ideal gas

As mentioned at the beginning of this section, if the gas under consideration is air, it is reasonable to consider it as an ideal gas and to impose an affine relation between p f and ρ f . More precisely, the relation ( 13) can be written as (see [START_REF] Goren | Pore pressure evolution in deforming granular material: A general formulation and the infinitely stiff approximation[END_REF]):

ρ f = ρ 0 f 1 + p f p atm . ( 18
)
It is then possible, for instance using H = Q in [START_REF] Gesenhues | Finite element simulation of complex dense granular flows using a well-posed regularization of the µ(I)-rheology[END_REF], to derive from [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF] an equation describing the evolution of the pressure p f , namely

∂ t ((1 -ϕ)p f ) + div((1 -ϕ)p f u) + p atm div u = div(κ(ϕ)(p atm + p f )∇p f ).
Considering that the pore gas pressure in the granular medium is negligible compared to the atmospheric pressure, that is p f ≪ p atm , it is reasonable to approximate p f + p atm by p f in the right-hand side of the last equation. As a consequence, we obtain the following pressure "diffusion" equation

∂ t ((1 -ϕ)p f ) + div((1 -ϕ)p f u) + p atm div u = p atm div(κ(ϕ)∇p f ). (19)
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).

Remark 2. This convection-diffusion equation for the pore gas pressure is frequently used, sometimes in slightly different forms. Thus, in [17, equation [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF]] or in [28, equation (7)] corrected in [START_REF] Anghel | Erratum: Grains and gas flow: Molecular dynamics with hydrodynamic interactions[END_REF], the authors use the non-conservative form of (19), namely

(1 -ϕ) ∂ t p f + u • ∇p f + (p atm + p f )div u = p atm div(κ(ϕ)∇p f ),
and then use the approximation p atm + p f ≈ p atm . A similar equation is often used in the incompressible case, i.e. when ϕ is constant and div u = 0. It then reduces to a "classical" convection/diffusion equation, and even to a diffusion equation if the transport term is not taken into account. This is for example the case in [START_REF] Montserrat | Pore fluid pressure diffusion in defluidizing granular columns[END_REF][START_REF] Roche | Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective[END_REF].

As in the case of a general state law, it is possible to estimate the term A in the energy equation [START_REF] Carman | Fluid flow through granular beds[END_REF]. More precisely, if the state law is given by ( 18), we obtain H by using [START_REF] Goren | Pore pressure evolution in deforming granular material: A general formulation and the infinitely stiff approximation[END_REF], namely

H(ρ f ) = p atm ρ f ρ 0 f ln ρ f ρ 0 f -1 = (p atm + p f ) ln 1 + p f p atm -1 ,
so that we get

-A = d dt (1 -ϕ)(p atm + p f ) ln 1 + p f p atm -1 + κ(ϕ)|∇p f | 2 .

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 (isochoric). More precisely, the Drucker-Prager and µ(I) rheologies are linearly unstable in some regimes, see [27, Appendix C] for the viscous Drucker-Prager model and [START_REF] Barker | Well-posed and ill-posed behaviour of the µ(I)-rheology for granular flow[END_REF] for the µ(I)-rheology.

In several recent works, see [START_REF] Barker | Well-posed and ill-posed behaviour of the µ(I)-rheology for granular flow[END_REF][START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF][START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF][START_REF] Barker | Well-posedness and ill-posedness of single-phase models for suspensions[END_REF], the authors have shown that taking into account the dilation of the granular medium allows, in the two-dimensional case, to regularise these models and to remedy these instabilities. More precisely, it is shown that the rheology, i.e. the expression of the additional stress τ , and the dilatancy law, i.e. the expression of div u, should be defined in a concordant way.

In this section, we will use the same notations as in [START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF]. We define the deviatoric strain-rate tensor by

S = ∇u + (∇u) T 2 - 1 2 (div u) Id,
which is a symmetric and trace-less tensor. We also introduce the inertial number (see [START_REF] Andreotti | Les milieux granulaires-entre fluide et solide: Entre fluide et solide[END_REF]):

I = d |S| p/ρ s , (20)
where d represents the grain diameter of the granular medium under consideration, and where the matrix norm is defined by the second invariant of any symmetric tensor, namely

|S| 2 = 1 2 S : S = 1 2 i,j S ij S ij .
Note that in a recent paper by Barker et al [START_REF] Barker | Well-posedness and ill-posedness of single-phase models for suspensions[END_REF], the authors use the µ(J)-rheology, instead of µ(I), to describe fluidized granular flow, i.e. granular material immersed in water. The dimensionless number J, defined by

J = η f |S| p ,
is used for granular flows with a low Stokes number (St = ρsd 2 η f |S|) and is therefore well suited for granular flows in liquids. In the present study, we focus on granular flows with a relatively high Stokes number, which is the case in air with a viscosity about 50 times smaller than that of water, while taking into account the effect of the interstitial gas through the pressure equation [START_REF] Heyman | Compressibility regularizes the µ(I)-rheology for dense granular flows[END_REF].

Granular rheology

The rheology of a granular flow can then be described in a fairly general way by a relationship of the form

τ = Z(ϕ, I)p S |S| . (21)
To be rigorous, the relation [START_REF] Jackson | The dynamics of fluidized particles[END_REF], having no sense when S vanishes, should be rewritten as

τ : S = 2 Z(ϕ, I)p|S|,
|τ | ≤ Z(ϕ, I)p where τ is symmetric and trace-less.

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

Remark 3. In [START_REF] Barker | Well-posed and ill-posed behaviour of the µ(I)-rheology for granular flow[END_REF][START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF][START_REF] Barker | Well-posedness and ill-posedness of single-phase models for suspensions[END_REF][START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF], the authors write Y (ϕ, p, I) instead of Z(ϕ, I)p. The writing proposed here seems to be relevant since in all models investigated later the stress is proportional to the pressure. Moreover, it allows the conditions that will be proposed later to be written relatively simply.

It is important to note that the constitutive law (21) excludes the case of purely viscous rheologies, which write τ = η(ϕ, |S|)S. The addition of such a term to [START_REF] Jackson | The dynamics of fluidized particles[END_REF] does not change the announced results and will be discussed regularly in the following, see Remarks 4 and 7.

Dilatancy model

The quantity div u expresses the evolution of the elementary volumes of fluid under the action of a flow moving at the velocity 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/compression phenomena, we must therefore impose an additional law which specifies how the divergence of the velocity field depends on certain quantities which characterise the flow. One way of expressing the effects of dilation is to impose a relation of the form1 div u = 2|S|f (ϕ, p, I).

With the constitutive law [START_REF] Jackson | The dynamics of fluidized particles[END_REF], the dilatancy law [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF] and the pore gas pressure equation [START_REF] Heyman | Compressibility regularizes the µ(I)-rheology for dense granular flows[END_REF] added to the conservation of mass and momentum equations (respectively ( 6) and ( 8)), we have a complete system whose unknowns are the volume fraction ϕ, the granular velocity u and pressure p, and the pore gas pressure p f , namely

∂ t ϕ + div(ϕ u) = 0, (23)
∂ t ((1 -ϕ)p f ) + div((1 -ϕ)p f u) + p atm div u = p atm div(κ(ϕ)∇p f ), ( 24
)
ϕρ s ∂ t u + u • ∇u = ϕρ s g -∇p + div Z(ϕ, I)p S |S| -∇p f , (25)
div u = 2|S|f (ϕ, p, I).

Note that this system can also be written in non-conservative form. Indeed, introducing the notation d t = ∂ t + u • ∇ for the convective derivative, equations ( 23)-( 26) rewrite

d t ϕ + ϕ div u = 0, (1 -ϕ)d t p f + (p atm + p f )div u = p atm div(κ(ϕ)∇p f ), ϕρ s d t u = ϕρ s g -∇p + div Z(ϕ, I)p S |S| -∇p f , div u = 2|S|f (ϕ, p, I).
5 Main properties of the generic model

Energy dissipation

As announced in the Section 2 (see equation ( 11)), the quantity E 0 (ϕ, u) = 1 2 ϕρ s |u| 2 satisfies the following

equation d dt E 0 (ϕ, u) = ϕρ s g • u + A -B + C.
Moreover, we have seen in the previous section (see equation ( 16)), that the term A, due to the presence of the interstitial gas, satisfies

-A = d dt E 1 (ϕ, p f ) + κ(ϕ)|∇p f | 2 ,
where E 1 (ϕ, p f ) = (1 -ϕ)H(ρ f ), so that the energy equation [START_REF] Carman | Fluid flow through granular beds[END_REF] rewrites

d dt E 0 (ϕ, u) + d dt E 1 (ϕ, p f ) + κ(ϕ)|∇p f | 2 + B -C = ϕρ s g • u.
Let us now examine the effects of the rheology and the dilatancy law, specified by ( 21) and ( 22) respectively, on the energy of the system ( 23)- [START_REF] Lions | Compressible Models[END_REF]. Replacing the deviatoric stress tensor τ and the divergence of the velocity field by the relations ( 21) and ( 22) in the terms B and C appearing in the energy equation [START_REF] Carman | Fluid flow through granular beds[END_REF], we obtain

B -C = τ : ∇u -p div u = 2 (Z -f )p|S|.
Reporting this relation in [START_REF] Carman | Fluid flow through granular beds[END_REF], we deduce that the total energy defined by E = E 0 (ϕ, u) + E 1 (ϕ, p f ) satisfies the following equation

dE dt + D = ϕρ s g • u,
where

D = κ(ϕ)|∇p f | 2 + 2 (Z(ϕ, I) -f (ϕ, p, I))p|S|
is a dissipation rate. Indeed, if Z(ϕ, I) -f (ϕ, p, I) ≥ 0, the total energy of the fluidised granular model decreases over time in the absence of any external force, such as gravitational force. The positivity of the dissipation rate, which is expressed by the condition

Dissipation condition: Z ≥ f, (27)
is one of the main properties of the fluidised granular models proposed in the following, and can now be stated.

Theorem 1. Under the dissipation condition [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I)-rheology and lateral wall effects[END_REF] and without external force, the model ( 23)-( 26) has a decreasing energy.

Remark 4. The addition of a viscous term, as suggested in Remark 3, does not affect this result. On the contrary, adding a viscous contribution η(ϕ, |S|)S to the stress induces the term

D dis = 2 η(ϕ, |S|)|S| 2 > 0.
in the left-hand side of the energy equation. In this case, the dissipation rate is enhanced.

Consistency of the dilatation law at equilibrium

According to the arguments of Pailha and Pouliquen in [31, section 3.3] and those of Roux and Radjaï in [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF], 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 ϕ depends linearly on the inertial number I and that an equilibrium relation of the following form

ϕ eq (I) = ϕ max -∆ϕ I (28)
can be found. The parameters ϕ max and ∆ϕ are obtained by fitting [START_REF] Mcnamara | Grains and gas flow: Molecular dynamics with hydrodynamic interactions[END_REF] with experimental measurements (e.g. ϕ max = 0.6 and ∆ϕ = 0.2 in [START_REF] Forterre | Flows of dense granular media[END_REF]). Note that DEM (Discrete Element Methods) simulations confirm this result (see [START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF] for example). The relation [START_REF] Mcnamara | Grains and gas flow: Molecular dynamics with hydrodynamic interactions[END_REF] can be used to determine the volume fraction at equilibrium by setting ϕ = ϕ eq (I). From this empirical law and the definition of the inertial number [START_REF] Ionescu | Viscoplastic modeling of granular column collapse with pressure-dependent rheology[END_REF], a state equation for the granular flow can easily be derived (see for instance [START_REF] Heyman | Compressibility regularizes the µ(I)-rheology for dense granular flows[END_REF]), which expresses the pressure as

p = ρ s d∆ϕ|S| ϕ max -ϕ 2 .
It is therefore attractive to supplement the conservation equations of mass and momentum with the µ(I)rheology and the above state law, so that a closed compressible model is obtained. However, such system based on the µ(I), ϕ(I)-rheology is always ill-posed in the two-dimensional case as proved by Heyman et al in [START_REF] Heyman | Compressibility regularizes the µ(I)-rheology for dense granular flows[END_REF] (see also [START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF]). Barker et al in [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF] followed by Schaeffer et al in [START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF] developed another approach: a dilatancy function f specifies the rate of volume change through the equation [START_REF] Lions | Compressible Models[END_REF].

Although the empirical law ( 28) is not prescribed in the fluidised granular model ( 23)- [START_REF] Lions | Compressible Models[END_REF], nor in the compressible models of Barker et al [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF] and Schaeffer et al [START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF], it must be used to impose conditions on the function f so that the dilatancy law ( 26) is consistent with the physics of dense granular medium. This is what we aim to do in what follows by introducing equilibrium conditions on f . Remark 5. There exists other formulations approaching this relationship between I and ϕ eq . For instance, in [36, p.929], the authors suggest ϕ sch eq (I) = ϕ max -∆ϕ 1 + 1/I , arguing that this law would prevent ϕ from becoming negative for large values of I. In practice, we will show in Subsection 5.4 that, independently of this law, the volume fraction solution of the mass conservation equation always remains positive, which makes the more complex law proposed in [START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF] useless. Nevertheless, as this law of equilibrium is empirical, experiments make it possible to refine the relationship and improve its accuracy in some cases. Hence, in [33, p.3], the authors propose

ϕ rob eq (I) = ϕ max -A I a , div u > 0 I > I eq (ϕ) div u = 0 I = I eq (ϕ) div u < 0 I < I eq (ϕ)
Figure 1: Change in volume of a granular medium near an isochoric equilibrium state characterized by I = I eq (ϕ), or equivalently ϕ = ϕ eq (I).

with A = 0.1305 and a = 0.8156, whereas in [9, p.10], the following law

ϕ bre eq (I) = ϕ max 1 + I .
is suggested. However, in the present article, we will always use the relation (28) but any other reasonable choice could be considered.

In [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF], Roux and Radjaï proposed a dilatancy model inspired by critical state mechanics that relates the rate of volume change to the deviation from equilibrium, i.e.

div u = 2a|S|(ϕ -ϕ eq (I)) ( 29
)
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. Written in the form [START_REF] Midi | On dense granular flows[END_REF], it is clear, that when a deformation occurs (i.e |S| > 0), if the volume fraction exceeds the equilibrium fraction, then div u is positive, so the material dilates. On the other hand, if ϕ < ϕ eq (I), we have div u < 0 which induces a contraction of the granular medium. The granular flow is isochoric (div u = 0) when there is no deformation (i.e |S| = 0). This behavior is essential and must be reproduced by the generic dilatancy law [START_REF] Lions | Compressible Models[END_REF]. Note that the Roux and Radjaï's model [START_REF] Midi | On dense granular flows[END_REF] fits in the general framework proposed by Barker et al [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF] and used here, with the function f defined as: f (ϕ, p, |S|) = a(ϕ -ϕ eq (I)). The empirical law (28) can be equivalently written

I eq (ϕ) = ϕ max -ϕ ∆ϕ , (30)
which defines the inertial number at the equilibrium by I = I eq (ϕ). As we will see in Section 6, it is preferable, when we place ourselves in the framework of the linear stability theory developed by Barker et al in [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF], to work with this alternative form and to define the equilibrium conditions in terms of deviation of the inertial number I from the equilibrium state characterized by I eq (ϕ). We now enonce those conditions which ensure that the dilatancy function f is consistent with the Roux and Radjaï's model [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF].

Theorem 2. Under the following equilibrium conditions Equilibrium conditions:

    
f (ϕ, p, I eq (ϕ)) = 0, f (ϕ, p, I) > 0 for I > I eq (ϕ), f (ϕ, p, I) < 0 for I < I eq (ϕ), [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] the model ( 23)-( 26) is consistent with the physics described by Roux and Radjaï in [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF].

The conditions [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] describe the same behaviour of the granular material, when it deviates from the isochoric equilibrium position, as the Roux and Radjaï dilatancy model. Indeed, at equilibrium we have I = I eq (ϕ)

so that to the divergence of the velocity field equals zero. If I > I eq (ϕ) then the rate of volume change is positive, resulting in an expansion of the granular flow. On the other hand, when the inertial number is less than the equilibrium value, the divergence of the velocity field is negative, resulting in compression of the material. Figure 1 illustrates this behaviour. Note that the non-isochoric and fluidized granular models derived in Section 6 will all satisfy the equilibrium conditions (31).

Linear stability of the model

According to the work of Barker et al [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF], we know that the granular model ( 23), [START_REF] Lagrée | The granular column collapse as a continuum: Validity of a two-dimensional Navier-Stokes model with a µ(I)-rheology[END_REF] with p f = 0 and ( 26) is linearly stable as soon as the functions Z and f satisfy the three following properties First stability condition:

Z - I 2 ∂ I Z = f + I∂ I f. (32)
Second stability condition:

Z + I∂ I Z ≥ 0. ( 33
)
Third stability condition:

∂ p f - I 2p ∂ I f < 0. ( 34
)
The purpose of this section is to show that accounting for the presence of the pore gas in the granular medium as achieved in Section 3 by supplementing the mass and momentum conservation equations of the solid particles with the equation ( 24) and adding the term -∇p f in the right-hand side of (25), does not change the stability result stated in [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF], as shown below.

Theorem 3. Under the conditions (32), ( 33) and (34), the model ( 23)-( 26) is linearly stable.

Proof. To prove this result, we adopt the ideas of the proof made in [7, Section 3] or in [5, section 3]. We must estimate the effects of the additional equation describing the evolution of the pore gas pressure, namely [START_REF] Kozeny | Uber kapillare leitung des wassers im boden-aufstieg, versickerung und anwendung auf die bewasserung, sitzungsberichte der akademie der wissenschaften wien[END_REF], and its coupling with the momentum conservation equation through the gradient term -∇p f in the right-hand side of [START_REF] Lagrée | The granular column collapse as a continuum: Validity of a two-dimensional Navier-Stokes model with a µ(I)-rheology[END_REF].

Let us consider a solution V 0 = (ϕ 0 , u 0 , p 0 , p 0 f ) of the system of equations ( 23)- [START_REF] Lions | Compressible Models[END_REF]. The first step is to linearize this system around V 0 by looking for a perturbed solution in the form V = V 0 + Ṽ. As in [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF], we retain only the terms that are linear in the perturbation V and neglect most of the terms that are not of maximal order. As an example, the linearized version of equation ( 24) describing the evolution of p f writes

∂ t p f = -c∆ p f where c = p atm κ(ϕ 0 ) 1 -ϕ 0 . ( 35
)
In the next step, the coefficients in the resulting linear system are frozen and we look for exponential solutions V(t, x) = e iξ•x+λt V, in order to obtain an eigenvalue problem that can be written as λ V = M V.

By specifying the unknown p f , i.e. by decomposing V = U, p f , the matrix M takes the following form

M =       N 0 iξ 0 0 0 0 c|ξ| 2      
, where N exactly corresponds to the matrix obtained when the fluidisation by the pore gas pressure is not modeled. The term iξ comes from the pressure gradient ∇p f involved in the momentum equation ( 25) whereas the last line of the matrix M corresponds to the equation [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF]. In [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF], it is proved that the conditions (32), ( 33) and [START_REF] Roche | Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective[END_REF] imply that all eigenvalues of the matrix N are positive. Since c|ξ| 2 ≥ 0, these conditions are also sufficient in the present case, i.e. with the additional pressure p f . Remark 6. The second condition [START_REF] Robinson | Complex behavior in compressible nonisochoric granular flows[END_REF] is not exactly the same as the one given in [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF], which is stronger since it requires a strict inequality ∂ I Z > 0. Looking in details at the proof in [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF], it is clear that the result remains valid if I∂ I Z + Z ≥ 0. Indeed, this condition is used to show that the trace of a matrix, namely

Tr(M + N) = I∂ I Y 2∥D ⋆ ∥ + Y 2∥D ⋆ ∥ |ξ| 2 + 1 Γ (1 + B) 2 ξ 2 1 + (1 -B) 2 ξ 2 2 ,
is strictly positive. Note the above equation is written with the notations used in [7, Proof of lemma 4.1]. Since Γ > 0, according to [START_REF] Roche | Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective[END_REF], the second term in the right-hand side of the above equality is always strictly positive. Therefore, the condition I∂ I Y + Y ≥ 0 is sufficient to ensure the positiveness of the trace. Note that the conditions for the linear stability of the model established in [START_REF] Barker | Well-posed continuum equations for granular flow with compressibility and µ(I)-rheology[END_REF] are only sufficient conditions. There is no evidence that they are optimal.

Volume fraction bounds

In order for the model to make sense, it is fundamental to ensure that the volume fraction ϕ always remains positive, and does not exceed the maximum value ϕ max previously introduced. We are now stating one of the main results of this paper.

Theorem 4. Suppose that the initial condition satisfies 0 ≤ ϕ| t=0 ≤ ϕ max and consider any smooth solution to (23)-( 26). We have ϕ ≥ 0.

Moreover if the consistency conditions [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] are satisfied then we have

ϕ ≤ ϕ max .
Proof. The proof is based on the following observation: If (ϕ, u) are regular and satisfy the equation ∂ t ϕ + div(ϕu) = 0 then for any function β ∈ C 1 (R, R), we have

∂ t β(ϕ) + div(β(ϕ)u) + (ϕβ ′ (ϕ) -β(ϕ)) div u = 0. ( 36
)
Step 1: Positiveness of ϕ. In order to show that ϕ remains positive as soon as ϕ| t=0 ≥ 0, we will use the above equation with the functions β ε defined, for all ε > 0, by

β ε (ϕ) =        0 if ϕ ≥ ε, -ϕ if ϕ ≤ -ε, 1 4ε (ϕ -ε) 2 if |ϕ| ≤ ε.
It is easy to show that β ε is of class C 1 and that

ϕβ ′ ε (ϕ) -β ε (ϕ) =    0 if |ϕ| ≥ ε, 1 4ε (ϕ 2 -ε 2 ) if |ϕ| ≤ ε.
In particular, we have

|ϕβ ′ ε (ϕ) -β ε (ϕ)| ≤ ε 4
so that after integrating the equation ( 36), we obtain

d dt β ε (ϕ) ≤ ε 4 |div u|.
By passing to the limit when ε tends to 0, we deduce that

d dt β(ϕ) ≤ 0 where β(ϕ) = 0 si ϕ ≥ 0, -ϕ si ϕ ≤ 0.
The condition ϕ| t=0 ≥ 0 means that β(ϕ| t=0 ) = 0. Since β(ϕ) decreases with time, and is by definition always positive, we deduce that, for any time, β(ϕ) = 0 which means that ϕ ≥ 0.

Step 2: Upper bound on ϕ. In general, there is no reason why a solution ϕ of the equation ∂ t ϕ + div(ϕu) = 0 should remain bounded if the divergence of the field u is not zero. Here, we will use the following additional information which comes from the equilibrium conditions [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF]. The equation ( 26) specifying the divergence of the velocity field implies:

ϕ ≥ ϕ max =⇒ div u ≥ 0. ( 37
)
Indeed, if ϕ ≥ ϕ max then ϕ ≥ ϕ max -∆ϕ I which is equivalent to I ≥ I eq (ϕ). According to the condition [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] we have, in this case, f ≥ 0 which implies div u ≥ 0. Next, we use the same argument as for showing the positiveness of ϕ but with the function β defined by

β(ϕ) = 0 if ϕ ≤ ϕ max , (ϕ -ϕ max ) 2 if ϕ ≥ ϕ max .
This function is of class C 1 and we have

ϕβ ′ (ϕ) -β(ϕ) = 0 if ϕ ≤ ϕ max , ϕ 2 -ϕ 2 max if ϕ ≥ ϕ max .
Using this expression and the implication (37), we notice that (ϕβ ′ (ϕ)-β(ϕ))div u ≥ 0 so that the integration with respect to the space variables of the equation ( 36) implies

d dt β(ϕ) ≤ 0.
Thus, if we have ϕ| t=0 ≤ ϕ max then β(ϕ| t=0 ) = 0 and β(ϕ) = 0 so that the result follows, namely ϕ ≤ ϕ max .

Physical examples

The objective of this section is to complete the generic fluidised ( 23)-( 26) model presented in Section 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, precisely the energy of the system must be dissipated, the volume fraction must be positive and upper bounded, and the model must be linearly stable. These last three properties are verified if the assumptions of Theorems 1, 3 and 4 are fulfilled, which is the case if the dissipation condition [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I)-rheology and lateral wall effects[END_REF], the stability conditions (32)-( 34) and the equilibrium consistency condition (31) are satisfied. Note that due to the Theorem 2, this last condition requires that the sign of the function f corresponds 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 function f must therefore fulfil the condition (31) 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 [START_REF] Forterre | Flows of dense granular media[END_REF]). Later, the µ(I)-rheology emerged thanks to the work of the Groupement de Recherche Milieux Divisés [START_REF] Midi | On dense granular flows[END_REF] 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 [START_REF] Pitman | A two-fluid model for avalanche and debris flows[END_REF], which can be reduced to a differential equation to the function f , is used to determine f . At this stage, f is known up to a constant which is independent on the inertial number but may depend on the volume fraction ϕ and the pressure p. 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 by Roux and Radjaï in [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF]. The last step will be to make sure that the other conditions ( 27), ( 33), ( 34) and (37) are fulfilled. Following this approach, we derive hereafter several fluidised and non-isochoric granular models. Before going into the details of the choice of the constitutive law and the determination of the rate of volume change, let us mention that the equation ( 32) is a linear in Z and f . Thus, any linear combination of pairs of functions (Z, f ) verifying [START_REF] Pitman | A two-fluid model for avalanche and debris flows[END_REF] will also be a solution of [START_REF] Pitman | A two-fluid model for avalanche and debris flows[END_REF]. Moreover, as no derivative with respect to the variable ϕ appears in the calculation of f (see the condition (32)), the coefficients of such linear combinations may depend on ϕ. This observation will facilitate the developments presented in what follows.

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.

Stabilisation of the Drucker-Prager model using a divergence condition

Granular media have been studied extensively (see for instance [START_REF] Forterre | Flows of dense granular media[END_REF] and the references therein) and several rheologies have been proposed. It is therefore useful to have in mind the usual forms that the constitutive law, defining the function Z, can take in the literature. One of the first choices comes from solid mechanics and is known as Drucker-Prager model. Following the book of Andreotti et al [3, Table 4.1] the deviatoric stress tensor is described by

τ = sin(δ) p S |S| ,
where δ corresponds to 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 [3, p.144] for the details).

In this case, the function Z previously introduced reads Z(ϕ, I) = sin(δ).

The stability condition [START_REF] Pitman | A two-fluid model for avalanche and debris flows[END_REF] implies that ∂ I (If (ϕ, p, I)) = sin(δ). Integrating this equation, we deduce that there exists a function K independent of I such that If (ϕ, p, I) = sin(δ) I + K(ϕ, p). In order to determine the function K, we invoke the consistency condition at equilibrium [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] and we obtain

f (ϕ, p, I) = sin(δ) 1 - I eq (ϕ) I . (38)
Noting that I eq (ϕ) and I are positive, it immediately follows that Z ≥ f so that the dissipation condition [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I)-rheology and lateral wall effects[END_REF] holds. The equilibrium condition [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] and the first stability condition [START_REF] Pitman | A two-fluid model for avalanche and debris flows[END_REF] are also fulfilled by construction of f from Z. The two other stability conditions [START_REF] Robinson | Complex behavior in compressible nonisochoric granular flows[END_REF] and [START_REF] Roche | Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective[END_REF], which in this particular case reads Z +I∂ I Z ≥ 0 and ∂ I f > 0, are obvious. By using the definition (20) of the inertial number I and the expression of f found above (38), the dilatancy law ( 22) can be rewritten

div u = 2 sin(δ)|S| - 2 sin(δ) d √ ρ s I eq (ϕ) √ p.
From a physical point of view, this relationship highlights a competition between shear effects which tend to increase the volume (solid particles separate when the material is sheared) and pressure effects which compress the flow. Figure 2 illustrates these effects on the granular medium.

Using the definition of the equilibrium via I eq (ϕ) or via ϕ eq (I), see [START_REF] Montserrat | Pore fluid pressure diffusion in defluidizing granular columns[END_REF] or [START_REF] Mcnamara | Grains and gas flow: Molecular dynamics with hydrodynamic interactions[END_REF], it is possible to rewrite the function f obtained in (38) as follows f (ϕ, p, I) = sin(δ) ∆ϕ I ϕ -ϕ eq (I) .

Written in this form, the dilatancy model obtained is close to the one proposed by Roux and Radjaï [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF] for which we have f (ϕ, p, I) = a ϕ -ϕ eq (I) .

We are now able to write a fluidized non-isochoric granular model based on the Drucker-Prager rheology, which is linearly stable and physically consistent, namely

∂ t ϕ + div(ϕ u) = 0, ∂ t ((1 -ϕ)p f ) + div((1 -ϕ)p f u) + p atm div u = p atm div(κ(ϕ)∇p f ), ϕρ s ∂ t u + u • ∇u = ϕρ s g -∇p + div sin(δ)p S |S| -∇p f , div u = 2 sin(δ)|S| -2λ sin(δ)(ϕ max -ϕ) √ p, (40)
where λ = 1 ∆ϕd √ ρs . Moreover, as we have previously seen in the Subsection 5.1, the energy for this model is dissipated over time, namely

dE dt + D = ϕρ s g • u,
where

D = 2λ sin(δ) (ϕ max -ϕ)p √ p + κ(ϕ)|∇p f | 2 .
Note that this rheology is known to be linearly unstable when associated with the incompressibility condition div u = 0, see [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I)-rheology and lateral wall effects[END_REF]Appendix C]. Here, by coupling it with the dilatancy law (40), we obtain a model satisfying the stability conditions ( 32), ( 33) and (34). where the function µ depends on three parameters I 0 > 0, µ 2 > µ 1 > 0 and is given by

µ(I) = µ 1 + µ 2 -µ 1 1 + I 0 /I . (41)
Following [14, p.10] or [3, p.226], typical values of the constants obtained for mono-dispersed glass beads are µ 1 = tan(21 o ), µ 2 = tan (33 o ) and I 0 = 0.3. Note that when I tends to 0, the Drucker-Prager model is recovered, which means that µ 1 = sin(δ) must be related to the internal friction. It may also be noted that these values are obtained experimentally, and that some authors specify that I 0 can depend on ϕ (see for instance [18, equation (2.12)] or [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF]Appendix A]).

Let us now see if it is possible to derive a complete and stable model with this choice of rheology. First of all, the condition (33) is satisfied since µ > 0 and µ ′ > 0. In order to ensure that the condition (32) is satisfied, the function f must be a solution of the differential equation

∂ I (If ) = µ(I) - I 2 µ ′ (I).
By integrating this ordinary differential equation, with respect to I, and using the condition f (ϕ, p, I eq (ϕ)) = 0, we find f (ϕ, p, I) = F (I) -I eq (ϕ) I F (I eq (ϕ)),

where F (I) = 3 2I M (I) -1 2 µ(I), M being the primitive of µ which vanishes at I = 0. Using the explicit expression (41), we derive

F (I) = µ 2 + µ 2 -µ 1 2 H I I 0 where H(x) = 1 1 + x - 3 ln(1 + x) x .
Finally, imposing the third stability condition [START_REF] Roche | Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective[END_REF] amounts to writing

F ′ (I) + I eq (ϕ) I 2 F (I eq (ϕ)) > 0. ( 42
)
F ′ (I) + I eq (ϕ)

I 2 F (I eq (ϕ)) = µ 2 -µ 1 2I 0 H ′ I I 0 + I eq (ϕ) I 2 µ 2 + µ 2 -µ 1 2 H I I 0 > 0,
which means that the condition (34) (equivalent to (42)) is always satisfied.

In the case of the µ(I)-rheology, the energy of the model is also dissipated, as

D = 2 p|S|(µ(I) -F (I)) + 2 d √ ρ s I eq (ϕ)F (I eq (ϕ))p √ p is positive since µ(I) -F (I) = 3(µ 2 -µ 1 ) 2 ln(1 + I/I 0 ) I/I 0 - 1 1 + I/I 0 > 0.
As a consequence, the fluidised and non-isochoric granular model, based on the µ(I)-rheology, satisfies all the required properties and reads

∂ t ϕ + div(ϕ u) = 0, ∂ t ((1 -ϕ)p f ) + div((1 -ϕ)p f u) + p atm div u = p atm div(κ(ϕ)∇p f ), ϕρ s ∂ t u + u • ∇u = ϕρ s g -∇p + div µ(I)p S |S| -∇p f , div u = 2F (I)|S| -2 F (ϕ) √ p,
where the functions F and F are defined by

F (I) = µ 2 + µ 2 -µ 1 2 1 1 + I/I 0 - 3 ln(1 + I/I 0 ) I/I 0 , F (ϕ) = 1 d √ ρ s
I eq (ϕ)F (I eq (ϕ)) with I eq (ϕ) = ϕ max -ϕ ∆ϕ .

Of course, when I goes to 0, we recover the Drucker-Prager model, presented in the previous subsection, with µ 1 = sin(δ) since lim

Stabilisation for a generic rheology using a divergence condition

Since the stability conditions (32), ( 33) and [START_REF] Roche | Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective[END_REF] do not involve a derivative with respect to the volume fraction ϕ, it is also possible to propose a law of the µ(I)-type taking into account the dependence with respect to ϕ. This means that, prescribing Z(ϕ, I), integrating the stability condition (32) and using the equilibrium condition (31), we find

f (ϕ, p, I) = W (ϕ, I) -I eq (ϕ)W (ϕ, I eq (ϕ)) 1 I , ( 43
)
where W is defined from Z by

W (ϕ, I) = 3 2I I I1 Z(ϕ, J) dJ - Z(ϕ, I) 2 .
Note that the choice of the constant I 1 has no influence on the value of the function f since replacing W by W + 3I1 2I does not change f . The choice of I 1 can be used in such a way that the integral makes sense, for example if Z admits a singularity at a point which one wishes to ignore. By performing an expansion in the ϕ variable in the vicinity of the equilibrium ϕ = ϕ eq (I), it is possible to rewrite f as

f (ϕ, p, I) = W (ϕ eq (I), I) ∆ϕ I ϕ -ϕ eq (I) + O ϕ -ϕ eq (I) 2 , Z(ϕ, I) f (ϕ, p, I) Contribution to dissipation D 1 1 - I eq (ϕ) I 2I eq (ϕ) d √ ρ s p √ p I 1 4 I - I eq (ϕ) 2 I 3 4 I + 1 4 I eq (ϕ) 2 d √ ρ s p √ p I n , (n ̸ = -1) 2 -n 2(n + 1) I n - I eq (ϕ) n+1 I 3n 2(n + 1) I n + 2 -n 2(n + 1) I eq (ϕ) n+1 d √ ρ s p √ p
Table 1: Dilatancy functions f (ϕ, p, I) for yield functions Z(ϕ, p, I) given as powers of the inertial number.

The last column shows the corresponding contribution to the dissipation rate D, in the energy equation.

which again is similar to the model (39) and therefore takes the same form as the one proposed in [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF] with a = W (ϕeq(I),I) ∆ϕ I . In the Table 1, we show the values of f obtained by the general formula (43) when polynomial functions in the variable I are chosen for the function Z. Because of the linearity of the Z → f correspondence, as indicated in the introduction to the Section 6, it is possible to write Z as linear combinations of each row of the table with coefficients depending on ϕ.

Remark 7. The case Z(ϕ, I) = I 2 , corresponding to the case n = 2 in the Table 1, is particular. Indeed, the corresponding dilatancy function is f = 0. Using the definition of the inertial number and the deviatoric stress tensor, such a choice of function Z corresponds to the following rheology

τ = Z(ϕ, I)p S |S| = d 2 ρ s |S|S,
which is purely viscous. More generally, viscous rheology may be written as

τ = η(ϕ, |S|)S. (44)
Such rheologies are not treated in the present study as noted in the Remarks 3 and 4. However, the stability condition equivalent to (32) in [7, p.5], shows that f = 0 is compatible with (44). It will therefore be possible to add any viscous contribution to the stress without changing the stability conditions. Note also that the case where Z = I n with n > 2 produces a negative dissipation term. It is therefore likely that the choice of rheology containing such terms is not dissipative, and should therefore be excluded.

Rheologies with dilatancy effects

The Ducker-Prager rheology with dilatancy effect

In order to take into account the variations of the volume fraction of the granular medium in the rheology, Wood [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF] proposed to modify the Drucker-Prager model by introducing what is called a dilatation angle ψ. This angle ψ 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). Intuitively, we can guess that a close packing that has to dilate in order to deform (ψ > 0) has a friction coefficient larger than that of a loose packing that will undergo compaction (ψ < 0). The granular viscosity is then an increasing function of ψ, and must be equal to sin(δ) when the dilatancy angle is zero. More precisely, the resulting laws, proposed in [3, p.150-151], are written as In this two-dimensional case, the definition of ψ via tan(ψ) = dY dX coincides with sin(ψ) = div u 2|S| .

τ = sin(δ + ψ) p S |S| div u = 2|S| sin(ψ). ( 45
The definition of the dilatancy angle ψ in the case of a 3-dimensional flow is different. In [3, p.150-151], the authors define ψ by the formula (45) by replacing the constant 2 by 3, while explaining that this definition no longer coincides with that given in the case of simple shear. Nevertheless, it is possible to do the same reasoning as in dimension 2 as shown below.

dX ∂ z u dY ∂ z v dZ ∂ z w = div u ℓ ψ
Assuming that u = (u(z), v(z), w(z)). We have div u = ∂ z w and

S =      -1 3 ∂ z w 0 1 2 ∂ z u 0 -1 3 ∂ z w 1 2 ∂ z v 1 2 ∂ z u 1 2 ∂ z v
When the dilatancy angle is small, the above relation is a close approximation to (45). The coefficient 2 remains valid in the three-dimensional case.

If the dilatancy angle ψ remains small, we can use the following approximations2 τ ≈ sin(δ) + cos(δ)ψ p S |S| div u ≈ 2|S|ψ.

In terms of the notations previously introduced, these laws corresponds to the following functions Z(ϕ, I) = sin(δ) + cos(δ)ψ ( 46)

f (ϕ, p, I) = ψ. (47)
In a general case, the dilatancy angle ψ can depend on the inertial number I and the volume fraction ϕ. This angle must be zero at equilibrium, i.e. ψ(ϕ eq (I), I) = 0. The usual closure proposed by Roux and Radjaï [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF] consists in using the first term of the expansion near the equilibrium state, namely ψ(ϕ, I) ≈ a(ϕ -ϕ eq (I)) where a = ∂ψ ∂ϕ (ϕ eq (I), I).

(

) 48
But it is clear that the choice (48), which implies the values of Z via (46) and f via (47) does not ensure that the stability conditions are fulfilled. In particular, the relation (32) linking Z and f is not satisfied. The question is now to determine whether it is possible to close the system, i.e. to propose an expression for the expansion angle ψ, so that the stability conditions are satisfied. Using (46) and (47), we can rewrite the condition (32) as

For common materials such as glass beads or sand, the angle of friction δ is of the order of 30 o . The power β is then of the order of 0.1. By using this dilatancy angle, the equations ( 46) and (47) give an expression for Z and f , and thus a complete model which fulfills, by construction, the stability condition [START_REF] Pitman | A two-fluid model for avalanche and debris flows[END_REF]. It is not difficult to verify that the other two stability conditions, namely ( 33) and ( 34), are satisfied.

Remark 9. For dilatancy angle ψ given by (50), it is therefore possible to provide a value for the "constant" introduced in the Roux and Radjaï's model [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF]. Indeed, using the equation (48), we find a = 2 sin(δ) 2 + cos(δ)

1 ∆ϕ I .
Finally, the fluidized and non-isochoric granular model based on the Drucker-Prager rheology incorporating the dilatancy effects writes

∂ t ϕ + div(ϕ u) = 0, (51)

Figure 2 :

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

6. 2

 2 Stabilisation of the µ(I)-rheology using a divergence conditionIn the class of consitutive laws for granular media, one of them is almost unanimously accepted: this is the µ(I)-rheology. It is basically a kind of generalisation of the Drucker-Prager model in which the effective viscosity depends on the inertial number I. With the notations introduced in this article, this model corresponds to Z(ϕ, I) = µ(I),

I→0 µ(I) = µ 1 F

 1 and limI→0 (I) = µ 1 .

) Remark 8 .= 1 2 -

 82 As for the friction angle, some authors use the tangent function instead of the sine function. Specifically, the dilatancy angle is the angle of motion relative to the horizontal arising from displacement, with dY = tan(ψ)dX where dY and dX are the vertical and horizontal displacements. This definition of ψ is specific to planar shear but can be generalised by the equation (45). The diagram below represents the dilatancy angle for a two-dimensional flow whose velocity field depends only on the vertical variable y. dX ∂ y u dY ∂ y v = div u 2 | S | ψ In a two-dimensional case and assuming that the velocity field is written u = (u(y), v(y)), we have div u = ∂ y v and S ∂ y v ∂ y u ∂ y u ∂ y v so that 2|S| = (∂ y u) 2 + (∂ y v) 2 .

. 2 .

 2 Hence 2|S| = (∂ z u) 2 + (∂ z v) 2 + 4 3 (∂ z w) 2 so that the length ℓ (see figure) writes ℓ 2 = 4|S| 2 -1 3 (div u)It is then possible to express the divergence of the velocity field as follows

( 2 +

 2 cos(δ))I∂ I ψ + 2(1 -cos(δ))ψ = 2 sin(δ). (49)By specifying that at the equilibrium, we must have ψ(ϕ, I eq (ϕ)) = 0, we solve (49) and we deduce that the dilatancy angle must have the following expression

∂ t (( 1 - 2 . 6 . 4 . 2

 12642 ϕ)p f ) + div((1 -ϕ)p f u) + p atm div u = p atm div(κ(ϕ)∇p f ),ϕρ s ∂ t u + u • ∇u = ϕρ s g -∇p + div sin(δ) + cos(δ)ψ(ϕ, I) p S |S| -∇p f , div u = 2|S|ψ(ϕ, I),(52)where ψ(ϕ, I) is given by (50). This model dissipates energy over time since the condition[START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I)-rheology and lateral wall effects[END_REF] (that is Z ≥ f ) is satisfied. Indeed, we haveZ -f = sin(δ) I eq (ϕ) I β .so that the energy equation rewritesdE dt + D = ϕρ s g • u,where D = 2 sin(δ) I eq (ϕ) I β p|S| + κ(ϕ)|∇p f | The µ(I)-rheology with dilatancy effect As with the Drucker-Prager model, the dilatancy effects can be taken into account in the µ(I)-rheology. In [33], Robinson et al proposed to modify the stress τ and the divergence condition as follows τ = µ(I) + ψ p S |S| div u = 2|S|ψ.

We will see in Remark 8 that this choice can be motivated by introducing a dilatancy angle.

Curiously, it seems that the development used in[3, p.133] (corrected in the next french edition, but repeated in[START_REF] Robinson | Complex behavior in compressible nonisochoric granular flows[END_REF]) is not correct; the coefficient cos(δ) being omitted. For their works, this does not change anything since this constant can be incorporated into the value of a hereafter.

Acknowledgments

This is contribution no. XXX of the ClerVolc program of the International Research Center for Disaster Sciences and Sustainable Development of the University of Clermont Auvergne. Declaration of Interests. The authors report no conflict of interest.

Using the same method as previously, we can determine the dilatancy angle ψ so that the model is stable. In that case, the condition [START_REF] Pitman | A two-fluid model for avalanche and debris flows[END_REF] becomes

By knowing the expression of the function µ (see (41)), it is easy to obtain the solution ψ which vanishes for I = I eq (ϕ), namely

with

The complete model is then very similar to the model ( 51)-( 52), the dilatancy angle being given by the relation (53). In this case, since Z -f = µ, the energy dissipation is given by

Note also that in this case, the dilatancy law (52) can be approximated near the equilibrium by the Roux and Radjaï's model [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF] div u = 2a(ϕ -ϕ eq (I))|S| where the "constant" a is given by a = G ′ (I) ∆ϕ .

Concluding remarks

In this paper, we have studied non-isochoric fluidized granular models, i.e. 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 very devastating, 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 [START_REF] Roche | Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective[END_REF]) 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 [START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF], 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 fluidized granular models is non-isochoric as it allows for volume change, namely the granular flow can expand when sheared. Moreover, we impose that the model be linearly stable, that it dissipate energy (over time), that it be compatible with the Roux and Radjaï [START_REF] Roux | Texture-dependent rigid-plastic behavior[END_REF] dilatancy model, and that the volume fraction, which is the solution of the mass conservation equation, be positive and bounded at all times. Working with this theoretical framework, we have studied dilatancy laws that are compatible with classical rheologies, i.e. Drucker-Prager and µ(I), with or without taking into account the effects of dilation in the yield condition. The main objective of this work was to derive fluidized 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 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 their validity, in particular by comparing the results of the simulations with laboratory experiments such as those described in [START_REF] Roche | Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective[END_REF]. These are the outlines of our future work on the subject.