Constitutive relations for compressible granular flow in the inertial regime

Granular flows occur in a wide range of situations of practical interest to industry, in our natural environment and in our everyday lives. This paper focuses on granular flow in the so-called inertial regime, when the rheology is independent of the very large particle stiffness. Such flows have been modelled with the $\unicode[STIX]{x1D707}(I),\unicode[STIX]{x1D6F7}(I)$ -rheology, which postulates that the bulk friction coefficient $\unicode[STIX]{x1D707}$ (i.e. the ratio of the shear stress to the pressure) and the solids volume fraction $\unicode[STIX]{x1D719}$ are functions of the inertial number $I$ only. Although the $\unicode[STIX]{x1D707}(I),\unicode[STIX]{x1D6F7}(I)$ -rheology has been validated in steady state against both experiments and discrete particle simulations in several different geometries, it has recently been shown that this theory is mathematically ill-posed in time-dependent problems. As a direct result, computations using this rheology may blow up exponentially, with a growth rate that tends to infinity as the discretization length tends to zero, as explicitly demonstrated in this paper for the first time. Such catastrophic instability due to ill-posedness is a common issue when developing new mathematical models and implies that either some important physics is missing or the model has not been properly formulated. In this paper an alternative to the $\unicode[STIX]{x1D707}(I),\unicode[STIX]{x1D6F7}(I)$ -rheology that does not suffer from such defects is proposed. In the framework of compressible $I$ -dependent rheology (CIDR), new constitutive laws for the inertial regime are introduced; these match the well-established $\unicode[STIX]{x1D707}(I)$ and $\unicode[STIX]{x1D6F7}(I)$ relations in the steady-state limit and at the same time are well-posed for all deformations and all packing densities. Time-dependent numerical solutions of the resultant equations are performed to demonstrate that the new inertial CIDR model leads to numerical convergence towards physically realistic solutions that are supported by discrete element method simulations.

Granular flows occur in a wide range of situations of practical interest to industry, in our natural environment and in our everyday lives. This paper focuses on granular flow in the so-called inertial regime, when the rheology is independent of the very large particle stiffness. Such flows have been modelled with the µ(I), Φ(I)-rheology, which postulates that the bulk friction coefficient µ (i.e. the ratio of the shear stress to the pressure) and the solids volume fraction φ are functions of the inertial number I only. Although the µ(I), Φ(I)-rheology has been validated in steady state against both experiments and discrete particle simulations in several different geometries, it has recently been shown that this theory is mathematically ill-posed in time-dependent problems. As a direct result, computations using this rheology may blow up exponentially, with a growth rate that tends to infinity as the discretization length tends to zero, as explicitly demonstrated in this paper for the first time. Such catastrophic instability due to ill-posedness is a common issue when developing new mathematical models and implies that either some important physics is missing or the model has not been properly formulated. In this paper an alternative to the µ(I), Φ(I)-rheology that does not suffer from such defects is proposed. In the framework of compressible I-dependent rheology (CIDR), new constitutive laws for the inertial regime are introduced; these match the well-established µ(I) and Φ(I) relations in the steady-state limit and at the same time are well-posed for all deformations and all packing densities. Time-dependent numerical solutions of the resultant equations are performed to demonstrate that the new inertial CIDR model leads to numerical convergence towards physically realistic solutions that are supported by discrete element method simulations.

Introduction
The original incompressible µ(I)-rheology was proposed (GDR MiDi 2004; to describe granular flow in which the particles are rigid and the solids volume fraction φ is constant and uniform. The key physical insight behind the theory was that, under these circumstances, the only non-dimensional groups are the bulk friction coefficient µ (the ratio τ /p of shear stress to pressure) and the inertial number where d is the grain diameter,γ is the shear rate, p is the pressure and ρ * is the intrinsic grain density. When forming a constitutive law from these groups, dimensionless analysis implies that one group must be a function of the other, i.e. τ p = µ(I).
In reality granular flow is compressible and the solids volume fraction φ may vary. With compressibility there are multiple non-dimensional groups, which greatly complicates the possible rheology. Nevertheless, GDR MiDi (2004), da Cruz et al. (2005) and Pouliquen et al. (2006) found that, in steady-state simulations and experiments, φ depended on just the inertial number where Φ is a monotonically decreasing function of I. Note that (1.3) implies Bagnold scaling (Bagnold 1954) for the pressure (i.e. it scales with the square of strain rate). Specifically, if Ψ (φ) is the inverse function of Φ(I), then (1.1) may be rearranged to yield (1.4) This is consistent with the discrete element method (DEM) simulations of Chialvo, Sun & Sundaresan (2012) in the inertial regime, in which the solids volume fraction φ is less than a critical volume fraction φ c , namely the jamming point (Liu & Nagel 1998). For φ > φ c , the stiffness of the particles becomes important and the rheology changes to either a 'quasi-static' or an 'intermediate' regime, which both depart from the Bagnold scaling. In general, the critical volume fraction φ c is dependent on the polydispersity of the grain-size distribution as well as the interparticle friction. These observations led to the compressible µ(I), Φ(I)-rheology (GDR MiDi 2004;Pouliquen et al. 2006;Forterre & Pouliquen 2008) in which the scalar relation (1.3), as well as (1.2), is assumed to hold even in non-steady situations. Although it might seem that compressibility would reduce the tendency to ill-posedness, in fact the compressible µ(I), Φ(I)-rheology is even more prone to ill-posedness in time-dependent calculations than incompressible µ(I)-rheology (Heyman et al. 2017). Such instability is perhaps not surprising since the equations are not always dissipative (see appendix B). To demonstrate this ill-posedness, § 3 presents computations for 928 D. G. Schaeffer and others a one-dimensional gravity-free shear cell in which the µ(I), Φ(I)-rheology blows up catastrophically once the mesh size is refined sufficiently, even though the initial conditions are just a small perturbation of the steady solution.
The purpose of this paper is to develop a viable alternative to µ(I), Φ(I)-rheology that preserves its successes. This alternative is based on the compressible I-dependent rheology (CIDR) introduced recently in . After briefly recalling the general formulation of CIDR, specific yield and dilatancy functions are introduced, which ensure that the theory is well-posed and reduces to µ(I), Φ(I)-rheology in the correct limits. Specifically, the new inertial CIDR equations recover both the µ(I) and Φ(I) relations (1.2)-(1.3) when the flow is isochoric (div u = 0), which in many geometries corresponds to steady flow. Sample computations with the one-dimensional inertial CIDR equations in a gravity-free shear cell converge to the steady-state solution, even for initial conditions that are a long way from steady state. Moreover, these computations agree well with transient DEM simulations of the same flows.
In § § 2 and 3 the µ(I), Φ(I)-rheology is described and computations that blow up because of ill-posedness are presented. Sections 4 and 5 introduce CIDR for the inertial regime and show that computations converge to steady state, and these computations are supported by DEM simulations described in § 6. The final two sections discuss a number of related issues and summarize the overall conclusions. Appendix A derives conditions for ill-posedness of the one-dimensional µ(I), Φ(I)-rheology in a gravity-free shear cell, appendix B relates the CIDR equations to the thermodynamic analysis of Goddard & Lee (2018) and appendix C describes the formulation of the DEM simulations. In a continuum formulation, dense granular flow is described by the solids fraction φ, the velocity vector u and the symmetric Cauchy stress tensor σ . In two dimensions (to which this paper is restricted) this formulation constitutes six scalar unknown functions of position and time. These satisfy governing equations including three conservation laws: one for conservation of mass (2.1) and two momentum conservation laws where ρ * is the intrinsic density of the grains and body forces have been neglected. In addition to the three conservation laws, three constitutive relations are needed to close this system. To write these constitutive laws, it is convenient to decompose the stress tensor into a pressure term p = −σ ii /2 plus a trace-free tensor τ , called the shear-stress tensor or the deviatoric stress tensor, such that (2. 3) The constitutive law for the incompressible µ(I)-rheology Forterre & Pouliquen 2008;Barker et al. 2015) is formulated in terms of the (total) strain-rate but for the compressible generalizations of the µ(I)-rheology (GDR MiDi 2004;Pouliquen et al. 2006;Forterre & Pouliquen 2008;Heyman et al. 2017;Goddard & Lee 2018) it is also useful to define the deviatoric strain-rate tensor where S ij is trace free. The alignment constitutive law, which is imposed in both µ(I), Φ(I)-rheology and CIDR, states that where the notation T = tr(T 2 )/2 denotes the second invariant of any symmetric tensor T . In particular, the eigenvectors of the deviatoric stress tensor and the deviatoric strain-rate tensor are parallel. This matrix equation is in fact equivalent to just one scalar equation and relies upon the implicit assumption that S = 0; i.e. that material is actually shearing. The other two constitutive laws in the µ(I), Φ(I)-rheology specify the deviatoric stress τ = τ and the solids fraction φ: (2.8) in terms of the pressure p and the inertial number where µ(I) is the bulk friction coefficient and d is the grain diameter. The constitutive laws (2.6)-(2.8) are referred to as µ(I), Φ(I)-rheology. Commonly used functional forms for µ(I) and Φ(I) in the constitutive laws (2.7) and (2.8) are where µ s , µ d , I 0 , φ c and a are constant material parameters (GDR MiDi 2004;Jop et al. 2006;Forterre & Pouliquen 2008;Trulsson et al. 2013). As shown in figure 1, these functions provide a good fit for the DEM simulations performed in this paper (see appendix C). Table 1 lists the best-fit parameters extracted from the steady-state DEM simulations. Note that figure 1 includes data for multiple cases with differing particle stiffness and system size, which shows that the parameters do not depend on the details of the DEM simulations. Incidentally the form Φ(I) = φ c − (φ c − φ min )/(I * /I + 1), where φ min and I * are constants, might be preferable to (2.11), which becomes negative for large I. However, the simpler form (2.11) is adequate for most practical purposes, because I rarely becomes large enough to drive Φ(I) negative.
Since (2.11) is a monotonically decreasing function, equation (2.8) can be inverted to write the inertial number as a function of the solids volume fraction which, for the specific function (2.11) used here, gives (2.13) From this it is also possible to determine an equation of state for the pressure by substituting (2.12) into the inertial number (2.9) and solving for p to give (2.14) which, for fixed φ, exhibits the Bagnold scaling in that the pressure scales with the square of strain rate as in (1.4). Note that this pressure tends to infinity as φ → φ c , behaviour that is discussed in § 7.2.
3. Catastrophic failure with the µ(I), Φ(I)-rheology As mentioned in the introduction, it follows from Heyman et al. (2017) that the dynamic equations of µ(I), Φ(I)-rheology for two-dimensional flow are always illposed. In this section it is demonstrated that ill-posedness can contaminate even the simplest of problems: flow in a two-dimensional shear cell that depends on only one spatial variable.

Equations for one-dimensional flow in a shear cell
In a planar shear cell granular material is confined from above and below by two parallel flat frictional walls, whose relative motion, in the absence of gravity, provides the only driving for the flow. This is a popular idealized geometry for the study of granular rheology, and the set-up can be either volume controlled, by fixing the wall separation distance H, or pressure controlled, by applying a normal force at the walls. In the following investigations H will be fixed and solutions will be restricted to those which are invariant in the shearing direction x and depend only on the vertical coordinate z and time. As such, this reduces the two-dimensional problem to one spatial dimension in the interval 0 < z < H. As the flow is taken to be compressible, both components of the velocity u = (u, w) may be non-zero and depend on z and t. For this special class of flows, the conservation laws (2.1)-(2.2) simplify to since all x-derivatives vanish. In addition, the deviatoric strain rate reduces to and its second invariant becomes It follows that in the shear cell the equation of state (2.14) reduces to which provides an explicit expression for the pressure. The alignment condition (2.6) implies that the relevant components of the deviatoric stress are which using the constitutive law (2.7) and equations (3.5) and (2.12) implies Equations (3.6) and (3.8) may then be substituted into the conservation laws (3.1)-(3.3) to obtain a system of three partial differential equations (PDEs) for φ, u and w that are first order in time and second order in space. For a given set of initial conditions φ(z, 0), u(z, 0), w(z, 0) the three PDEs must be solved subject to the boundary conditions that there is no slip at the walls where V 0 is the velocity of the top wall. It is easily verified that the steady linear shearing solution where φ 0 is the average initial solids volume fraction satisfies (3.1)-(3.3) given (3.6), (3.8) and the boundary conditions (3.9). It is therefore anticipated that any set of initial conditions will tend to this solution. 932 D. G. Schaeffer and others 3.2. Blow up and grid dependence with the µ(I), Φ(I)-rheology For one-dimensional flow in a shear cell, the physical variables are non-dimensionalized with the scalings where the hats indicate non-dimensional variables. Equations (3.1)-(3.3) with the relations (3.6) and (3.8) therefore reduce to the non-dimensional system where the grains size d, the grain density ρ * and the wall velocity V 0 scale completely out of the system. The conservation laws (3.13)-(3.15) are solved by the numerical method of lines (Schiesser 2012). This method discretizes the spatial derivatives in the PDEs which generates a system of coupled ordinary differential equations. These are then integrated forward in time using MATLAB's ODE15s routine. Two different methods are used to approximate the spatial derivatives: the first is a finite difference method using first-order differences, while the second method uses a Chebyshev spectral scheme (Canuto et al. 1988;Trefethen 2000).
The non-dimensional height of the cellĤ is chosen to equal 30 (i.e. the physical height of the cell is 30 grain diameters). The initial conditions are plotted in figure 2 and are identical to the non-dimensionalized steady solution except that the vertical velocity has a small smooth imperfection in the centre of the domain: where is a small parameter. It should be noted that this form does not satisfy the boundary conditions exactly but does within numerical precision. Provided is not too small, the ill-posedness condition (A 19), outlined in appendix A, is satisfied in the centre of the domain, as shown in figure 2(d), and it is therefore anticipated that numerical problems will develop. Confirming this anticipation, various snapshots of the non-dimensional vertical velocity are shown in figure 3. (The full time evolutions ofû(ẑ,t),ŵ(ẑ,t) and φ(ẑ,t) are shown in supplementary movies 1-3.) In the high-resolution plots (figure 3a,b), the growth of noise causes the numerical method to break down at the indicated time -at this point integration tolerances can no longer be met. Note that the misbehaviour originates in the centre of the domain, where the ill-posedness condition (A 19) is satisfied. Moreover, although the two solutions blow up at similar times, their spatial structure is completely different.
In contrast, a low-resolution (N z = 47) simulation using the finite difference scheme does not blow up and in fact decays towards the steady state, as shown in figures 3(c) and 4(a). (The same occurs for the low-resolution run with the Chebyshev spectral  scheme (see figure 4a).) This is exactly the behaviour that one might expect of a wellposed model, but here it is entirely due to the higher numerical diffusion in the lowresolution scheme. On grid refinement this diffusion is no longer sufficient to suppress the underlying ill-posedness of the equations, and instabilities develop. For all values of N z above a certain threshold (see figure 4b), the initial perturbation is amplified indefinitely, causing the method to fail. The timet * at which the Chebyshev spectral discretization fails is plotted in figure 4(b) as a function of the number of grid points N z . Higher spatial resolution computations resolve higher wavenumber instabilities that grow at a faster rate, leading to earlier failure.

Inertial compressible I-dependent rheology (iCIDR)
The CIDR is a general framework that retains the conservation laws (2.1)-(2.2) and the alignment condition (2.6), but it modifies the other two constitutive equations. The constitutive law for the deviatoric stress (2.7) is replaced by assuming that there is a yield condition such that, if material is deforming, then where Y is called the 'yield function'. In addition, it is assumed that the density evolves dynamically according to a flow rule that is reminiscent of critical state soil  mechanics (Schofield & Wroth 1968;Jackson 1983): where f is called the 'dilatancy function'.  showed that the CIDR equations are linearly well-posed provided the yield and dilatancy functions satisfy the following three conditions:  . The computations use the method of lines with a finite difference (dashed lines) and a Chebyshev (solids lines) discretization at both high (red lines) and low (blue lines) spatial resolution. The timet * = 0.01052 at which the Chebyshev scheme becomes non-convergent is indicated. (b) The final timet * before convergence failure is plotted for the Chebyshev spectral discretization as a function of the number of grid points N z . The same plot (not shown) for failure of the finite difference computation is qualitatively similar, but differs in detail.
A key result of this paper is the introduction of new constitutive functions which are motivated by the inertial regime. These are constructed to satisfy both the wellposedness conditions (4.3)-(4.5) and the observed asymptotic steady-state behaviour, known as the µ(I) and Φ(I) relations (1.2)-(1.3). These relationships are derived from isochoric (constant volume) flows in steady state (see figure 1 and appendix C). For the purpose of deriving Y and f they may be conveniently stated as follows: where Ψ (φ) is the inverse function of Φ(I). Even with these constraints there are infinitely many possible choices for the yield function and dilatancy function. What might be described as the simplest acceptable choice for these functions is now proposed. The starting point for this new theory is the relation Y = µ(Ψ (φ))p in (4.6). However, this cannot lead to a well-posed theory as ∂ I Y = 0 and thus (4.4) is not satisfied. Taking instead the simplest non-trivial I-dependence gives which reduces to the µ(I), Φ(I)-rheology when I = Ψ (φ). Figure 5 is a useful aid for comparing this formula with µ(I), Φ(I)-rheology. The corresponding dilatancy function is then found through substitution of (4.7) into the well-posedness PDE (4.3). Taking contributions from both the homogeneous and particular solutions of (4.3) and imposing (4.6) gives the expression This theory will be termed 'inertial CIDR' or iCIDR, since, as will be shown, it captures the Bagnold scaling. In addition to satisfying (4.3) and (4.6), it is readily verified that the iCIDR constitutive functions (4.7)-(4.8) satisfy the inequalities (4.4)-(4.5). As such, the dynamic equations which result from iCIDR are guaranteed to be linearly well-posed for all deformations and for all values of the solids volume fraction.
In µ(I), Φ(I)-rheology I is tied rigidly to φ through the I = Ψ (φ) relation (2.12), whereas for iCIDR I evolves through the flow rule (4.2) given the dilatancy function (4.8). Substituting (4.8) into (4.2) yields the quadratic equation which has the positive root (4.10) The right-hand side of (4.10) is a function of the solids volume fraction φ and the ratio div u/ S . Equation (2.9) can be used to determine an equation of state for the pressure (4.11) This equation differs from (2.14) only in that Ψ (φ) is replaced by I defined by (4.10). Since (4.10) depends on the velocity gradient only through the ratio div u/ S , which is unchanged by scaling the velocity gradient, (4.11) satisfies Bagnold scaling (Bagnold 1954). According to (4.9), if S → 0, then I tends to either zero or infinity, depending on the sign of div u. Hence, it seems that the iCIDR constitutive laws may break down if S = 0. However, the inertial number I can be bypassed by calculating an alternative expression for the pressure. Substituting (2.9) into (4.9) and solving for p gives (4.12) which is well defined for all deformation rates. Note that p is strictly positive unless S = 0 and div u 0, (4.13a,b) in which case the pressure is zero. Thus, in the absence of shear, there is no pressure if the grains are diverging from one another, but there is a finite, positive pressure if they are converging. Moreover, using the alignment condition (2.6), the yield function (4.7) and the definition of I, it follows that the deviatoric stress is given by which is also well defined for all deformation rates. Finally, note that if φ → φ c , then Ψ (φ) → 0 and the pressure tends to infinity (provided there is some straining).

One-dimensional flow in a shear cell with the iCIDR rheology
The response of the iCIDR rheology is now tested in the one-dimensional gravitationless shear cell, as studied in § 3.2 for the µ(I), Φ(I)-rheology. In this geometry the equation of state (4.12) implies the relations for the pressure

D. G. Schaeffer and others
and the deviatoric stresses (4.14) Substituting these expressions into the conservation laws (3.1)-(3.3) and using the scalings (3.12) implies the non-dimensional iCIDR equations are where the pressure is given bỹ (5.6) The pressure (5.6) can be substituted directly into (5.4)-(5.5) to produce a system of three equations for φ,û andŵ. Unlike the incompressible µ(I)-rheology and the incompressible Navier-Stokes equations, in which pressure must be determined globally as part of solving the equations, all terms in this system (5.4)-(5.5) are explicitly specified locally. This makes the development of numerical methods much simpler. Figure 6(a) and supplementary movies 4 and 5 show numerical solutions according to iCIDR for the one-dimensional gravity-free flow of § 3, specifically for initial conditions (3.16). Note that the solution converges smoothly to steady state, the same steady state as for the µ(I), Φ(I)-rheology (3.10), without any sign of the catastrophic resolution-dependent blow-up seen before. The simulations are computed using the method of lines with the finite difference discretization for both N z = 47 and 201 grid points. Note that the solutions in figure 6(a) lie directly on top of one another! This shows that these solutions are grid converged.
The fact that the iCIDR equations are mathematically well-posed suggests that they can handle larger perturbations from steady state than (3.16). An interesting test case is an initial condition with |∂ẑû| = 0 for at least one point. This is interesting because the function χ, defined in (A 18), is infinite when |∂ẑû| = 0, so the µ(I), Φ(I)-rheology is strongly unstable even at low resolution. Following this idea, the results of a simulation with the initial condition are plotted in figure 7. As figure 7(d) shows, with the parameters φ 0 = 0.78,â = 0.16 andb = 0.1, there is a large central region where χ > 0 and two points where χ is infinite. This problem is therefore very strongly ill-posed for the µ(I), Φ(I)rheology. Ill-posed behaviour has been verified numerically; indeed, for these initial conditions, even for the coarse grid with N z = 47, the solution blows up (not shown). By contrast, with the iCIDR equations (5.3)-(5.6) there is no catastrophic blow-up and the simulations smoothly evolve from the initial conditions towards the steady-state order of one time unit. The solids volume fraction starts out independent ofẑ, but it develops a non-trivial profile which changes sign and then decays much more slowly. The latter evolution is responsible for the local maximum of |ŵ| neart = 3 shown in figure 7(e).

Comparison of iCIDR with DEM simulations
The DEM calculations detailed in appendix C and figure 1, which were used to recover the steady µ(I) and Φ(I) relations, are also capable of verifying the time-dependent characteristics of the iCIDR solutions. Here new DEM simulations are initialized with the same procedure that was used to obtain the steady linear shearing solution (3.10), but then the velocity fields are replaced with the sinusoidally perturbed profiles (5.7) discussed in the previous section. As expected, this results in a decay from these applied fields back towards the steady solution. Figure 8 shows that the iCIDR solutions differ by less than 2.5 % relative error from the DEM simulations throughout the dynamics. Figure 8(a-c) shows that iCIDR captures the spatial variation of the three flow variables φ,û andŵ. Figure 8(d) indicates that it also captures complex details of the time evolution. Although more testing is needed, this agreement indicates that iCIDR correctly represents significant aspects of the rheology.
7. Discussion of related issues 7.1. Remarks on ill-posedness Issues of ill-posedness in µ(I)-rheology were first raised in Barker et al. (2015). It was shown that, although the dynamic equations derived from incompressible µ(I)rheology are mathematically well-posed for a large range of inertial numbers, the system is ill-posed when I is too high or too low. The following remarks may help reconcile this ill-posedness with the fact that problems of practical interest (e.g. Lagrée et al. 2011;Staron et al. 2012) have apparently been successfully solved numerically using µ(I)-rheology.
In the first place, ill-posedness effects may be masked in simulations performed on a low-resolution grid, as in § 4 of this paper. Specifically, numerical diffusion may be sufficient to suppress the instability. Ill-posedness may become apparent only through careful comparison of progressively mesh-refined simulations -see figure 4 of  and figure 21 of Martin et al. (2017); in these papers certain spurious flow features continue to become ever more finely scaled as the grid size gets smaller. It is sometimes suggested that low-resolution solutions of an ill-posed model might be sufficient for some practical purposes. However, such approaches are scientifically flawed, because the results rely on numerical diffusion to regularize the problem, which is dependent on both grid size and numerical scheme, and is often not known precisely. In our view it is far better to try to understand what physics is missing in the model, and only compute solutions when a well-posed theory has been formulated.
Other issues are the following. (i) In some problems, such as column collapse (Lagrée et al. 2011), the ill-posed region of parameter space may only be active for a short period of time. In such cases careful comparison of numerical results at different spatial resolutions, including some very fine grids, may be needed for non-convergence to become apparent (Barker et al. 2015;Martin et al. 2017). (ii) Ill-posedness may also be partially suppressed by attempts to remove the singularity in the viscosity at low strain rates. Many numerical codes do this by introducing an upper bound on the viscosity, which implies that the material reverts to a Newtonian fluid for slow flows. However, these procedures are ad hoc in nature, and there is no guarantee that ill-posedness is suppressed completely.
In this subsection attention has been focused on the consequences of ill-posedness that may be seen through computations. However, on the mathematical side, it 942 D. G. Schaeffer and others where = 0, cannot be solved on any non-zero interval {0 < t < t 0 } unless p is an even non-negative integer. The singularities of the initial conditions at x = nπ, n = 0, ±1, ±2, . . . (which are extremely mild if p is large) make such solutions impossible.
7.2. Behaviour at the boundary of and outside the inertial regime Regarding the boundary of the inertial regime, note in (2.14) and (5.11) that p tends to infinity as φ → φ c . In DEM simulations, Chialvo et al. (2012) found that this growth was cut off and blended into the pressure from the intermediate regime through the formula p blend = p inert p itm p inert + p itm , Thus the pressure tends to p itm as φ → φ c , where the limit pressure depends on the elastic modulus. Although it is beyond the scope of the present paper, it is anticipated that a more complete version of CIDR would remain valid across regime boundaries. Granular flow at densities φ > φ c (provided the strain rate is not too large) corresponds to what is called the quasi-static regime (Otsuki & Hayakawa 2009;Chialvo et al. 2012;Singh et al. 2015). In this regime stresses may remain non-zero even as the strain rate tends to zero, i.e. a static yield stress exists; the scale of these static stresses is set by k/d, where k is the spring constant in DEM simulations.
CIDR is a general theory that can model granular material outside, as well as inside, of the inertial regime. In particular, the version of CIDR in § 2(e) of , which was motivated by critical state soil mechanics (Schofield & Wroth 1968), has a non-zero static yield stress, as in the quasi-static regime in DEM simulations. (In that paper, to facilitate comparison with Silbert et al. (2001), the stress scale was chosen as ρ * gd, i.e. dependent on gravitational acceleration g. This was an unfortunate choice with no intrinsic significance. A more appropriate choice would have been k/d, where k is the spring constant in DEM simulations.) 7.3. Extensions to the theory Inertial CIDR is also able to incorporate non-monotonicity of the µ(I) function (DeGiuli & Wyart 2017), which is crucial for modelling hysteretic effects, such as coexisting static and moving regions, in depth-averaged avalanche models (Daerr & Douady 1999;Pouliquen & Forterre 2002;Edwards & Gray 2015;Edwards et al. 2017;Russell et al. 2019). Non-monotonic µ = µ(I) functions are problematic in the incompressible µ(I)-rheology, because they imply that the theory is always ill-posed in regions of decreasing friction (Barker et al. 2015). For iCIDR, however, having a non-monotonic µ = µ(I) function is not a problem, because it is formulated in terms of the solids volume fraction, i.e. µ = µ(Ψ (φ)), and so it does not affect the well-posedness conditions (4.3)-(4.5).
At present iCIDR is explicitly a local theory which cannot account for the observed role of fluctuations that inspired the non-local theories of Pouliquen & Forterre (2009), Kamrin & Koval (2012) and Bouzid et al. (2013). Inclusion of these effects is an important direction for future work. The iCIDR equations, introduced in this paper, provide a continuum model for fluid-like inertial flows of rigid spherical particles that lie below a critical volume fraction, above which the compressibility of the grains becomes important (Otsuki & Hayakawa 2009;Chialvo et al. 2012). One striking aspect of the iCIDR theory is its simplicity: it is a minimal extension of µ(I), Φ(I)-rheology with no extra variables, no extra parameters, no extra evolution equations beyond conservation of mass and momentum and no extra boundary conditions. While retaining the success of µ(I), Φ(I)-rheology for steady inertial flow, iCIDR is well-posed, thermodynamically sound (Onsager symmetric and dissipative) and agrees well with transient DEM simulations in a one-dimensional gravity-free shear cell.
Inertial CIDR is very well suited to numerical calculations because the pressure is defined by a local equation of state. This contrasts with incompressible theories, in which the pressure can only be found globally as part of solving the equations. The numerical simulations presented in this paper are therefore a particularly encouraging proof of concept and it is hoped that other existing numerical methods can similarly be modified in order to bring the advantages of the iCIDR formulation to a wide range of practical applications.
Although there are some technical complications, effectively the dynamic equations of µ(I), Φ(I)-rheology for two-dimensional flow are always ill-posed (Heyman et al. 2017). As in Barker et al. (2015), the ill-posedness has a directional character: plane-wave solutions (in all space) in certain directions suffer uncontrolled growth while those in other directions decay smoothly to uniform flow. If, as in § 3, solutions depending on only one spatial variable are sought, these one-dimensional equations may be either well-posed or ill-posed, depending on whether the one independent direction retained corresponds to one of the stable or unstable directions. In this appendix, conditions for the ill-posedness of the one-dimensional system (3.1)-(3.8) are derived.
The equations (3.1)-(3.3) together with (3.6) and (3.8) are linearized around a base flow (φ 0 , u 0 , w 0 ) in the dependent variables (φ, u, w) by introducing perturbations (φ,ȗ,w) of the form and then freezing the base flow coefficients. In the linearization the highest-order derivatives of the perturbed variables are retained in each equation, together with the convective derivatives. This leads to linearized equations of the form where a i and b ij are constant coefficients derived from the pressure (3.6) and the deviatoric stresses (3.8).
The linearized system (A 2) admits normal-mode solutions of the form in which ξ is the wavenumber, λ is the temporal growth rate andφ,ũ andw are constants. Substituting these into (A 2) reveals that λ must be an eigenvalue of the matrix The imaginary terms −iw 0 ξ on the diagonal shift the eigenvalues λ, but do not affect stability or well-posedness, which are governed by the real part of λ. The vertical component of the base state velocity w 0 is therefore set to zero in what follows. Denoting the bottom right 2 × 2 matrix of terms multiplying −ξ 2 as the eigenvalues are calculated by solving where the coefficients are To test for ill-posedness, the large-wavenumber limit ξ → ∞ is taken. Balancing the order of terms in (A 6), it is clear that λ(ξ ) ∼ ξ 2 , so the substitution λ = Λξ 2 is made. Then the terms of maximal order in (A 6), i.e. O(ξ 6 ), have the coefficient Λ 3 + (tr Γ )Λ 2 + (det Γ )Λ. Thus, to leading order, either Λ = 0, or Λ is a solution of i.e. one of the two roots To determine the signs of these roots, the coefficients in (A 2) are now evaluated. If the superscript 0 notation is dropped, the coefficients in (A 2) are then where all values are evaluated with the base-state fields From (A 10) and (A 5) it follows that The two roots in (A 9) are real since Since λ + ∼ Λ + ξ 2 as ξ → ∞, the system is linearly ill-posed if the larger root Λ + is positive, which occurs if either It may be seen from (A 12) and (A 13) that if tr Γ < 0 then det Γ < 0, so it suffices to consider only case (b). Thus, is a sufficient condition for ill-posedness. Conversely, suppose Appendix C. Details of the DEM calculations Two-dimensional DEM simulations (Cundall & Strack 1979) were performed in a shear-box geometry, both to confirm the well-established steady behaviour and to explore transient flows. The domain is a rectangle in (x, z) space, with all units non-dimensionalized with the scalings (3.12), such that 0 x <L and 0 ẑ <Ĥ. Boundary conditions and the system size are then designed to suppress confinement effects so that the volume can be taken to be a representative bulk element. Periodicity is enforced in thex direction and Lees-Edwards boundary conditions (Lees & Edwards 1972) are applied at the bottomẑ = 0 and topẑ =Ĥ of the domain. There is no gravity applied to the system and the only driving is provided by the difference in horizontal velocity between the top and bottom V 0 , as prescribed by the Lees-Edwards algorithm.
Details of the precise DEM simulation algorithm can be found in Otsuki & Hayakawa (2011) as an identical method is employed here. Normal forces f (n) between particles are calculated from a linear spring-dashpot arrangement with an associated spring constant k (n) and viscous dissipation constant η (n) . Tangential forces f (t) may either stick or slip depending on whether a Coulomb friction criterion, with particle friction constant µ p , is satisfied. Stick interactions are defined as those with |f (t) | < µ p |f (n) | and, like the normal force, are calculated from a linear spring-dashpot with parameters k (t) and η (t) . Interactions with greater computed tangential forces are labelled as slip events and the tangential force is truncated to |f (t) | = µ p |f (n) |. In this paper parameters are chosen which give k (n) /p > 10 4 so that calculations are in the rigid particle regime of da Cruz et al. (2005). Unless stated otherwise the values µ p = 0.4, k (n) = 10 4 and k (t) = 0.5k (n) are used. Particle interactions with these values are very short-lived so that results are invariant of the viscous dissipation. The tangential dashpot is therefore neglected (η (t) = 0) and η (n) = 4.2 is chosen, as in Silbert et al. (2001), so that the particles have a restitution coefficient of e = 0.9.
To select a mean packing fraction φ 0 , the system size along thex direction is determined as L =ĀN/(Hφ 0 ), whereĀ is the average grain area. The shear cell is then populated with N particles with density ρ * and mean diameter d and V 0 set to unity in order to match the non-dimensional iCIDR equations. In order to avoid crystallization effects, the individual particle diameters are chosen randomly from a discretized distribution. Here an even spread from 0.8d, 0.9d, d, 1.1d and 1.2d is taken so that the number of particles of each diameter is N/5. These particles, which are initially elastic but not frictional, are randomly distributed in the domain. This results in overlaps which would normally cause very large elastic forces, so firstly there is a period during which the total kinetic energy of all particles is scaled to a small constant value so that the system can reach an equilibrium state. The arrangement which results from this procedure has almost uniform packing density and very small overlaps. Then, the interaction algorithm is altered so that the particles are approximately rigid (very large k (n) ) and frictional. The true simulation begins after the velocity fields are prescribed.
The steady-state µ(I) and Φ(I) relations found in previous works (e.g. da Cruz et al. 2005) are first confirmed in order to obtain macroscopic rheological parameters. Both curves are derived from the same set of experiments in which the solids volume fractions φ 0 take values in the range 0.76-0.8. This range is expected to lie in the inertial regime and, due to the Φ(I) relation, each packing corresponds to a unique inertial number. Once the system has reached a steady state, the flow fields are coarse-grained. This is achieved by averaging in thex direction and then averaging within bins which discretizeẑ into boxes of height 2d. Each run is then repeated 10 times in order to calculate error estimates. The solids volume fraction φ and the two velocity componentsû andŵ are clearly defined in the DEM data and allow the inertial number to be readily calculated from its definition (1.1). Calculation of the bulk friction coefficient µ requires the macroscopic stress components to be defined. As in Silbert et al. (2001) this involves a sum over all particles α in the sampling volume: where r αβ is the centre-to-centre vector, f αβ is the total force between particle pairs and δû α = (û α −V 0ẑ /Ĥ,ŵ α ) is the velocity fluctuation of particle α. From the stress decomposition (2.3) these stress components are used to calculate the pressure p, deviatoric stress τ and hence µ across the domain. As these quantities are all found to be invariant ofẑ, the mean value is taken for the µ(I) data presented in figure 1(a). The volume fraction is also found to be constant at steady state so that the Φ(I) relation, plotted in figure 1(b), is simply verified using the mean packing density Φ = φ 0 .