Well-posed continuum equations for granular flow with compressibility and μ(I)-rheology.

Continuum modelling of granular flow has been plagued with the issue of ill-posed dynamic equations for a long time. Equations for incompressible, two-dimensional flow based on the Coulomb friction law are ill-posed regardless of the deformation, whereas the rate-dependent μ(I)-rheology is ill-posed when the non-dimensional inertial number I is too high or too low. Here, incorporating ideas from critical-state soil mechanics, we derive conditions for well-posedness of partial differential equations that combine compressibility with I-dependent rheology. When the I-dependence comes from a specific friction coefficient μ(I), our results show that, with compressibility, the equations are well-posed for all deformation rates provided that μ(I) satisfies certain minimal, physically natural, inequalities.

Introduction

Much effort has been devoted to formulating constitutive laws for continuum models of granular materials [1-5]. However, the lack of acceptable dynamic theories, i.e. well-posed equations in the sense of Joseph & Saut [6], for granular flow has severely hampered progress in modelling many geophysical and industrial problems. In the simplest class of models, flow is described by partial differential equations (PDEs) for the density, the velocity vector and the stress tensor; conceptually, such models are hardly more complicated than the Navier–Stokes equations. The equations represent conservation laws for mass and momentum coupled to constitutive equations to close the system. However, despite the appeal of their simplicity, they have been plagued with ill-posedness, i.e. small perturbations grow at an unbounded rate in the limit that their wavelength tends to zero [6]. Such behaviour is clearly unphysical. However, the immediate practical implication of ill-posedness is that numerical computations either blow-up, even at finite resolution, or do not converge to a well-defined solution as the grid is refined, i.e. the numerical results are grid dependent [7-10].

The first model of this type [2,11,12] specifies constitutive laws that represent a tensorial generalization of the work of de Coulomb [13] on earthwork fortifications. In the language of plasticity theory, it is a rate-independent, rigid/perfectly plastic model with a yield condition based on friction between the grains. However, it was shown to be ill-posed in all two-dimensional contexts and all realistic three-dimensional contexts [2]. Critical-state soil mechanics (CSSM) [1] is a sophisticated elaboration of Coulomb behaviour that allows for compressibility. It also suffers from ill-posedness, depending on the degree of consolidation. This ill-posedness is much less severe than for a Coulomb material [3,12], but is still physically unsatisfactory and introduces potential issues for the numerical simulation of transient granular flows. More recently, the μ(I)-rheology [4,5,14] introduces a modest amount of rate dependence into (incompressible) Coulomb behaviour through the non-dimensional inertial number, which is proportional to the shear rate and inversely proportional to the square root of the pressure. As shown in Barker et al. [9], this theory leads to well-posed (two dimensional) equations in a significant region of state space, but it is ill-posed at both low and high inertial numbers.

This paper presents an analysis of constitutive equations that extend the incompressible μ(I)-rheology of Jop et al. [5] to compressible deformations, through combination with CSSM. The main result is that, in two dimensions, compressible I-dependent equations can be made well-posed for all densities, for all stress states and for all deformation rates. In other words, to obtain well-posedness, Coulomb behaviour is modified by including only two natural, fairly small, perturbations of the theory, namely compressibility and rate dependence. Following this very general treatment, which has implications for many existing formulations [15-17], we elucidate our findings with an illustrative model that includes physically motivated features and reduces to the μ(I)-rheology in the incompressible limit. This has the advantage that it retains the conceptual simplicity of the original theory. Although we consider only two-dimensional flow, it should be noted that in numerous cases it has been found that flow in two dimensions is more prone to ill-posedness than in three [2,3,18]. Thus, we anticipate that the corresponding three-dimensional equations including these effects will also be well-posed.

Currently, a wide range of new constitutive laws for granular materials are being developed including the μ(I)-rheology [4,5], elasto-plastic formulations [19,20], non-local rheologies [21-24], kinetic theory [25], as well as Cosserat [26], micro-structural [27] and hypoplastic theories [28]. Enormous progress has been made over the past decade and there is the realistic and exciting prospect that practical granular flows, which span the solid-like, liquid-like and gaseous regimes, may shortly be described by continuum models. In this paper, we seek to understand one of the conceptually simplest formulations that leads to well-posed equations.

In §2, we introduce the equations to be studied and formulate our well-posedness result for them. This theorem is proved in §§3 and 4. In §5, we solve the new equations for steady, uniform chute flow. In two appendices, we summarize key ideas from CSSM and survey topics regarding ill-posed partial differential equations.

Governing equations

Dense granular flow is described by the solids volume fraction ϕ, the velocity vector and the stress tensor . In two dimensions, this constitutes six scalar unknowns that are spatially and temporally dependent. These are governed by conservation laws plus constitutive relations. Conservation of mass gives the scalar equation and conservation of momentum gives the vector equation where ρ* is the constant intrinsic grain density and g is the acceleration due to gravity. Closure of these equations is achieved through three constitutive relations.

The Coulomb constitutive model

For a Coulomb material, which is assumed to be incompressible, the first constitutive relation states that ϕ is a constant. This then reduces (2.1) to the For the next constitutive relation, the stress tensor is decomposed into a pressure term (where p=−σ/2) plus a trace-free tensor called the deviatoric stress. The second relation is then the where μ is a constant and for any tensor the norm is defined by This yield condition expresses the idea that a granular material cannot deform unless the shear stress is sufficient to overcome friction.[1] The third constitutive relation requires that the eigenvectors of the deviatoric stress tensor and the deviatoric strain-rate tensor[2] are aligned (see figure 1 for motivation), which may be written In words, the above equation may be interpreted as asserting that in the space of trace-free symmetric 2×2 matrices, which is two-dimensional, and are parallel. Thus, this matrix equation entails only one scalar relation. For reference below we record that

Figure 1.

(a) Illustrative stress eigenvectors; along the major axis the stress eigenvalue is −(p+∥∥), with the minus sign indicating compression. (b) A possible material deformation that is consistent with the stress field in (a). (Online version in colour.)

It is customary [2,5] to process these equations by expressing the deviatoric stress in terms of p and the strain rate as follows: where we have invoked (2.5) and (2.8). We may substitute (2.10) into (2.2) to obtain and the resulting equation, together with (2.3), gives three equations for pressure p and velocity . In form, at least, these equations resemble the incompressible Navier–Stokes equation. However, in two dimensions (as considered here) they are always ill-posed [2].

Incompressible μ(I)-rheology

Work described by the Groupement De Recherche Milieux Divisés [4] has significantly improved the Coulomb model by including some rate dependence (in the sense of plasticity [29]) in the yield condition while making no changes in the incompressible flow rule (2.3) and the alignment condition (2.8). Specifically, a wide range of experiments are captured by replacing the constant μ in (2.5) by an increasing function μ(I) of the inertial number, where d is the particle diameter. The expression where μ1, μ2 and I0 are constants with μ2>μ1, is a frequently used form [30]. Below we shall assume that

The modified yield condition changes (2.11) to read The effect of this seemingly small perturbation is profound. Unlike for Coulomb material, equations (2.15) and (2.3) are linearly well-posed for I, μ satisfying For the specific μ(I)-curve (2.13) this inequality covers a significant range of inertial numbers, specifically when the deformation rate is neither too small nor too large relative to the pressure. Outside of this range, the maximal-order linear stability analysis and numerical simulations show that perturbations grow exponentially with growth rates tending to infinity as their wavelength is reduced [9]. This behaviour is the hallmark of ill-posedness and leads to unphysical numerical solutions that strongly depend on the grid resolution used.

Compressibility and I-dependent rheology

We refer to CSSM (see appendix A) for guidance in introducing compressibility into the rheology. Thus, we make no change in the alignment condition (2.8); we assume ϕ-dependence in the yield condition, and we allow for volumetric changes by introducing a new function f(p,ϕ,I) and modifying the flow rule to To get well-posed equations, our analysis (see §3) shows that the yield condition and the flow-rule functions must be related by the equation[3] and satisfy the inequalities and

We may now state our main result, the well-posedness theorem for the system (2.1), (2.2), (2.8), (2.17), (2.18), which we call the CIDR equations. (Mnemonic: compressible I-dependent rheology.)

Theorem 2.1

Under hypotheses (2.19) and (2.20), the CIDR system is linearly well-posed.

The term linearly well-posed is defined in appendix B, and the result is proved in §§3 and 4. Table 1 contains a comparison between the conditions for linear well-posedness for the different constitutive models that have been discussed here.

Table 1.

Table of criteria for linear well-posedness.

model	Coulomb	incomp. μ(I)	CSSM	CIDR
conditions	always ill-posed	(2.16)	always ill-posed	(2.19), (2.20)

Remark 2.2

The I-dependence in these equations need not relate to a friction coefficient μ(I). In §2e, we connect the equations to μ(I)-rheology.

Derivation of evolution equations

To place the equations in a larger continuum-mechanics context, we show that the CIDR equations of motion can be rewritten as a system of three evolution equations for the velocity and the solids fraction ϕ. In form, these equations are analogous to the Navier–Stokes equations for a viscous, compressible fluid. We make no use of this form of the equations in our proof of well-posedness.

We want to eliminate stresses from the equations of motion. To this end, we propose to solve for the mean stress p using the flow rule (2.18), which we rewrite as Note that f(p,ϕ,I) depends on p both directly in its first argument and indirectly through in its third argument. However, which by assumption (2.20b) is non-zero. Thus, we may apply the implicit function theorem to (2.21) to solve p=P(∇,ϕ).[4] Given this, we may define and substitute into conservation of momentum to obtain an equation This equation, along with (2.1), gives a system of three evolution equations for the velocity and the solids fraction ϕ. It is possible that previous formulations of compressible μ(I) equations [15-17] may be seen as CIDR equations with specific constitutive laws specified. In this paper, we choose to elucidate the well-posedness result with more generic choices of f and Y in order to remain impartial.

Connection to μ(I)-rheology

Without making any attempt to be general, we illustrate one example of how μ(I)-rheology may be included in constitutive relations of the form (2.17), (2.18). Motivated by equation (.3) in appendix A, we make the ansatz and In these equations, it is worth emphasizing that p,ϕ,I are treated as independent variables, not to be confused with the dependence of I on p in the previous subsection. The function C(ϕ) is an increasing function of ϕ. As ϕ varies (with I fixed) the yield loci ∥∥=Y (p,ϕ,I) derived from (2.24a) form a nested family of convex curves in stress space (figure 2b). Observe from (2.18) that deformation without volumetric strain is possible if f(p,ϕ,I)=0; i.e. for (2.24b), if p/C(ϕ)=β(I)/2. Substituting this formula into (2.17) and using (2.24a), we derive for such isochoric deformation to be possible. Thus, to recover the yield condition ∥∥=μ(I)p of the μ(I)-rheology, let us require that

Figure 2.

(a) An example curve for the function C(ϕ) with a minimum solids volume fraction ϕmin and a vertical asymptote at ϕmax. (b) Nested yield surfaces of the form (2.24) for a fixed value of I with differing solids volume fractions. (The solid blue line, the dashed arrows and the labels A and B refer to a discussion of CSSM in appendix A.)

Lemma 2.3

Equations (2.19) and (2.25) imply that and

Proof.

Substituting the relations (2.24) into (2.19), and using (2.25) to eliminate β, we derive the linear ordinary differential equation for α=α(I): Solving this linear equation for α(I), with an integrating factor, we obtain from which the formula (2.26) follows after integrating the second integral by parts. Finally, substituting this formula for α(I) into (2.25), we obtain the formula (2.27) for β(I). ■

Lemma 2.4

The yield condition and flow-rule function (2.24a,b) that follow from (2.26), (2.27) verify hypotheses (2.19) and (2.20), provided μ(I) satisfies (2.14).

Of course (2.19) is satisfied because this equation was imposed in deriving (2.26), (2.27).

Differentiating (2.24b), we see that ∂f(p,ϕ,I)=−2/C(ϕ)<0. To calculate ∂f(p,ϕ,I), we first reparametrize the integral in (2.27) to obtain . Then By (2.14), μ′′(I)<0, so μ′(sI)>μ′(I) for 0

For inequality (2.20a), we reparametrize the integral (2.26) and differentiate to obtain as desired. ■

The incompressible limit

Based on an analogy with CSSM, let us suppose that C(ϕ) is a sensitive function of ϕ, say of the form where g is the constant of gravitational acceleration, b is a non-dimensional parameter and the factor ρ⋆gd gives C the dimensions of pressure. The non-dimensional function has an argument that is dependent on the minimum solids fraction ϕmin for sustained stress transmission between grains (random loose packing[5]) and ϕmax, which is the maximum packing fraction that can be attained. Typically, Δϕ=ϕmax−ϕmin is small. For definiteness we may take as in figure 2. Note that C(ϕ) diverges as ϕ→ϕmax; thus, (2.30) requires that ϕ is confined to a narrow range, In physical terms, the maximum solids fraction ϕmax represents the jamming threshold. We call the limit Δϕ→0 incompressible because, as may be seen from (2.32), the density of the material becomes essentially constant.

Lemma 2.5

As Δϕ→0, the CIDR equations reduce to the equations of incompressible μ(I)-rheology, (2.3), (2.15).

We process the CIDR equations, which have the six unknowns ϕ, u and σ, as follows. First, we reduce to five unknowns—ϕ, u, p and τ=∥∥—by recalling the definition (2.4) and the alignment condition (2.8) to write Next we use the yield condition to eliminate ϕ, reducing this number to four. Specifically, substituting (2.24a) into (2.17), we write the yield condition Solving (2.33) for ϕ we obtain where the dependence on ∇ comes from the fact that Substitution of this formula into the conservation laws (2.1), (2.2) yields the equations and Finally, we show that the flow rule (2.18) may be rewritten To see this, we combine (2.24b) with (2.25) to conclude and then substitute the relation α(I)=τ/p+p/C(ϕ) derived by manipulating (2.33). Thus, the system (2.35), (2.36) governs the evolution of the four unknowns u, p and τ.

Now we claim that if C(ϕ) has the form (2.30), then (2.35), (2.36) is a singular perturbation of (2.3), (2.15). It follows from (2.30) that (2.34) has the expansion Substituting (2.37) into the continuity equation (2.35a), we find If Δϕ=0, then this equation reduces to div =0, although this is of course a highly singular limit. Thus, if Δϕ=0, the left-hand side of (2.36) vanishes, so this equation simplifies to the yield condition τ=μ(I)p, and substitution into (2.35b) yields (2.15). This proves the lemma. ■

Proofs, part I: linearization

An alternative formulation of the alignment condition

It is convenient to study the linearized equations with a reformulated alignment condition that describes stress in terms of eigenvectors of, rather than entries of, the stress tensor. Since defined by (2.4) has trace zero, it has eigenvalues[6] ±∥∥. Taking ψ as the angle that the eigenvector with eigenvalue −∥∥ makes with the x1-axis gives which may be verified by checking that is an eigenvector of this matrix with eigenvalue −∥∥. Thus, the stress tensor σ is completely specified by the three scalars p, ∥∥ and ψ.

Focusing on the first rows of the strain-rate tensor (2.9) and of (3.1), we extract from the matrix equation (2.8) the vector equation where k=−2∥∥<0. Since and lie in the two-dimensional space of trace-free, symmetric matrices, (3.2) is equivalent to (2.8). It follows from (3.2) that In point of fact, this equation is slightly weaker than the alignment condition since (3.3) is consistent with the possibility that k>0 in (3.2); to rule out the latter possibility we impose the supplemental inequality[7] that

The calculation

Substitution of the stress tensor (3.1) into the momentum balance equations (2.2) allows for the full set of equations to be written as This system has five scalar unknowns, =(u1,u2,ϕ,p,ψ). In (3.5a), (3.5b), τ is a mnemonically suggestive abbreviation for the yield function Y (p,ϕ,I) in (2.17), and in (3.5d), a repetition of (2.18), the function f depends on arguments (p,ϕ,I) that are not written explicitly.

As in appendix B, to linearize the equations we substitute a perturbation of a base solution , say into the equations, retain only terms that are linear in the perturbation and freeze the coefficients at an arbitrary point . It is convenient to temporarily drop most terms not of maximal order and estimate their effect in a calculation at the end of the argument. For example, this construction applied to (3.5c) yields the constant-coefficient, linear equation where and . Lower-order terms and in the full linearization of (3.5c) have been dropped in (3.7).

In expanding the fully nonlinear factor ∥∥ in (3.5d), we may take advantage of the rotational invariance of the equations to arrange that ψ*=0; i.e. we may calculate in a rotated coordinate system for which, at , the x1-axis is the maximal stress axis. Then by the alignment condition (3.3), the base-state deviatoric strain-rate tensor is diagonal at and by (3.4), in the 1,1-position of this matrix, . This corresponds to non-zero compression along the major stress axis, as illustrated in figure 3. Now and where the approximation follows from the expansion if A>0 and |X|,|Y |≪A. Thus, as given in table 2, the (local) linearization of ∥∥ equals .

Figure 3.

An example of a base-state velocity field for the strain-rate tensor (3.8) with and . (Online version in colour.)

Table 2.

List of maximal-order linearizations of terms in (3.5a)–(3.5e), to assist in deriving (3.10a)–(3.10e). (only in this table).

term in (3.5a)–(3.5e)	contribution to (3.10a)–(3.10e)
∥D∥	−D^_11
I	−I^∗∥D^∗∥D^_11−I^∗2p^∗p^
∂_j[τ cos⁡(2ψ)]	(∂_pτ)^∗∂_jp^+(∂_ϕτ)^∗∂_jϕ^
	+ (∂_Iτ)^∗{−I^∗∥D^∗∥∂_jD^_11−I^∗2p^∗∂_jp^}
∂_j[τ sin⁡(2ψ)]	2τ^∗∂_jψ^
f∥D∥	−f^∗D^_11+∥D^∗∥(∂_pf)^∗p^+∥D^∗∥(∂_ϕf)^∗ϕ^
	+ ∥D^∗∥(∂_If)^∗{−I^∗∥D^∗∥D^_11−I^∗2p^∗p^}

In (3.5d), the function f contains p, ϕ and I as implicit arguments. As reflected in the table, the dependence on p and ϕ contributes zeroth-order terms in these variables to the linearization.

In (3.5a), (3.5b), τ also depends on p, ϕ and I, and the terms involving τ are differentiated; hence new issues arise in linearizing them. For example, by the chain rule, Since ψ*=0, the full linearization of, say, the first term here equals , a term given in the table, plus lower-order terms All of these terms, as well as numerous other analogous terms in the full linearization of (3.5a) that are not of maximal order, have been dropped in (3.10a)–(3.10e).

Putting all the pieces together, we obtain the linearization[8] of the system (3.5a)–(3.5e) where and Observe that, by hypothesis (2.19), B=C, a fact that we use in (4.3) and below.

Proofs, part II: calculation of growth rates

The eigenvalue problem

We now look for exponential solutions of (3.10a)–(3.10e), where is a 5-vector of scalars, =(ξ1,ξ2) is a vector wavenumber, 〈,〉 indicates the inner product and λ is the growth rate. The function (4.1) is a solution of (3.10a)–(3.10e) iff satisfies the generalized eigenvalue problem where , and On the right side of (4.2), the modified eigenvalue parameter is λ+i〈*,〉 because

Equation (4.2) is a generalized eigenvalue problem because , the matrix of coefficients of time-derivative terms, is not invertible. To extract an ordinary eigenvalue problem, we decompose into blocks where and 12, 21 and 22 fill out the rest of the matrix. Defining and , we rewrite (4.2) as The zero entries in the last two rows of mean that so we can solve for Substitution of into (4.7) then reduces this problem[9] to the ordinary 3×3 eigenvalue problem, where 11 is the 3×3 block in the upper left of .

We decompose the 3×3 matrix in (4.9) into smaller blocks, where we calculate as the contribution of 11, as the contribution of , which is symmetric, and

Estimation of the eigenvalues

We claim that the growth-rate eigenvalues (4.10) satisfy By compactness, it suffices to prove that Since only the real parts of eigenvalues matter, we may drop the term i〈*,〉 in (4.10) and verify (4.14) for the eigenvalue problem[10] where we write for the matrix in (4.10) and we shorten the notation by dropping the subscript 1 on . For large , it is instructive to use perturbation theory to compare the eigenvalues (4.15) with the eigenvalues , where

Lemma 4.1

Provided ≠0, the 2×2 matrix + is positive definite.

Since and are symmetric, it suffices to show that the trace and determinant of + are positive. According to (2.20), A>0 and Γ>0, from which it follows immediately that tr(+)>0.

Regarding the determinant, for any 2×2 matrices where accounts for the cross terms. For the specific matrices (4.11) and (4.12), , and

This proves the lemma. ■

Remark 4.2

It is noteworthy that except for the two directions Effectively, this calculation rederives the result of Pitman & Schaeffer [12] that the equations of CSSM, even without I-dependence, are well-posed for all directions except possibly those defined by (4.22).

It follows from lemma 4.1 that has two eigenvalues, say Λ1,Λ2, where Λ1,Λ2<0 and is homogeneous of degree 2 in . Since is an -perturbation of 0, two of the growth-rate eigenvalues of (4.15) satisfy both of which are negative in the limit ; i.e. they are bounded above by zero in this limit. The third growth rate is given by The first term on the extreme right is the ratio of two quartics, the denominator being non-zero, so it is bounded, and the perturbation decays at infinity. This verifies (4.14) for all three eigenvalues derived from (3.10a) to (3.10e).

It remains to consider the effect of the lower-order terms that were neglected in (3.10a)–(3.10e). Inclusion of these terms would lead, after a calculation as above, to an eigenvalue problem (4.15) for a perturbed matrix As above, two of the eigenvalues of this matrix are negative and , and invoking the determinant shows that the third is bounded. This verifies (4.14) for eigenvalues of the full linearization of (3.5a)–(3.5e) and hence shows that the system is linearly well-posed. It is, therefore, expected that numerical solutions of the full two-dimensional nonlinear transient equations will not exhibit the exponential blow-up of perturbations or the dependence on grid resolution seen in the incompressible μ(I) equations [9].

Chute flow in compressible I-dependent rheology

Let us recall the steady-uniform ‘Bagnold’ solution [4,31] to the (incompressible) μ(I)-rheology [4,5] for chute flow. Assuming that μ(I) is given by (2.13), these solutions exist for inclination angles ζ between the maximum and minimum angles and , respectively. Letting Oxz be a Cartesian coordinate system, with the x-axis pointing downslope and the z-axis being the upward pointing normal, the Bagnold solution is where h is the flow depth, ϕ⋆ is a constant solids volume fraction and u is the downslope velocity. At a fixed inclination, the inertial number is equal to the positive constant There is strong experimental and discrete element method (DEM) evidence [4,5,31,32] for both the lithostatic pressure distribution and the three-halves power in the dependence of the velocity on the thickness in (5.1). Fortunately, the CIDR solution is close to this, as we now show.

Assume that all variables depend only on z and only the x-component of the velocity is non-zero, say u=u(z). Motivated by Bagnold flow, the scalings are used to non-dimensionalize the variables and the following ordinary differential equations (ODEs) are derived for and below.

Lemma 5.1

Under the above assumptions, the non-dimensional pressure and where the pressure is zero at the free surface

Steady uniform solutions to the CIDR model are shown in figure 4 for χ=0.5, 2, 5, 20, 100 and with ϕmax=0.6 and ϕmin=0.5. The nature of the solution is controlled by the non-dimensional parameter χ, which is larger for thicker flows. This parameter is also inversely proportional to β(I), which by (5.2) and (2.27) is a weakly increasing function of the inclination angle as shown in figure 5. For large χ, the concentration graphed in figure 4a is close to ϕmax over a significant proportion of the flow depth, with a narrow boundary layer near the free surface, where the concentration decreases to ϕmin. In the limit as , the solution tends to the red line in figure 4a, which corresponds to Bagnold flow with ϕ⋆=ϕmax. As χ is decreased, the surface boundary layer becomes thicker and the flow becomes progressively more dilute, tending to the blue line in figure 4a as χ→0, which corresponds to Bagnold flow with ϕ⋆=ϕmin. In figure 4b,c, the non-dimensional pressure and the velocity are shown, which are also bounded between the red and blue lines corresponding to Bagnold flow with ϕ⋆ equal to ϕmax and ϕmin, respectively. As well as the profiles being bounded by two non-dimensional Bagnold solutions, the steady-uniform CIDR solutions also have exactly the same scaling properties (5.3), on the flow density ρ⋆, gravity g, the flow depth h, the chute inclination ζ, the particle diameter d and the inertial number I, as the Bagnold solution (5.1). The steady-uniform CIDR solutions are, therefore, very closely related to the classical Bagnold solution and are almost indistinguishable for large χ.

Figure 4.

Solutions for (a) the solids volume fraction ϕ, (b) the non-dimensional pressure and (c) the non-dimensional velocity as a function of the non-dimensional thickness for five different values of the non-dimensional parameter χ=0.5, 2, 5, 20, 100 (black lines). The red line corresponds to and Bagnold flow with ϕ⋆=ϕmax, while the blue line corresponds to χ→0 and Bagnold flow with ϕ⋆=ϕmin. In all cases, ϕmin=0.5 and ϕmax=0.6.

Figure 5.

The function β(I) as a function of the inclination angle ζ for the case ζ1=28°, ζ2=36° and I0=0.279. (Online version in colour.)

Proof of the lemma.

The form of the assumed solution implies that D=D=0 and D=(1/2) du/dz>0; in particular, div =0. It follows from the alignment condition (2.8) that σ=σ=−p and σ=τ, where τ=||||. Thus, the downslope and normal momentum balances (2.2) reduce to and We may apply the chain rule to (5.7) and (5.8) to derive , and, since tractions at the surface vanish, On the other hand, since div =0, the flow rule (2.18) implies f(p,ϕ,I)=0, and, given the ansatz (2.24b), we conclude that From (5.10), we have p/C(ϕ)=β(I)/2, so by (2.24b) and (2.25) Comparing (5.9) with (5.11), it follows that , and hence the inertial number I is equal to the constant I everywhere. The value of I is given by the formula (5.2), which is derived by inverting (2.13) with . Since the inertial number is constant and C(ϕ) satisfies (2.30), equation (5.10) can be inverted and then non-dimensionalized, using the scalings (5.3), to give the solids volume fraction in (5.6). Using the same scalings (5.3), the normal momentum component (5.8) yields the ODE (5.4), which can be solved for the pressure by integrating down from the free surface using (5.6). The concentration then follows from (5.6). From the definition of the inertial number (2.12) and that 2∥∥=|du/dz|, the downslope velocity satisfies the ODE and the scalings (5.3) then give the ODE (5.5), which, since is known, can be integrated upwards from the no-slip condition at the base. ■

Conclusion and discussion

In this paper, we have analysed a generalization of the μ(I)-rheology that allows for changes in the granular solids volume fraction. The equations of motion, the CIDR model, are found to be linearly well-posed when the constitutive laws satisfy certain criteria. We indicate (in §2e) how the specific I-dependence of the μ(I)-rheology can be fitted into the theory, and we have shown that the inclusion of compressibility removes the ill-posedness of the incompressible rheology [9].

CIDR equations for steady, uniform flow down an inclined chute are solved in §5, where we observe that the solution is comparable to the classical Bagnold solution under conditions similar to the many experimental and DEM results [4,5,31,32] for the velocity profile in steady chute flows. Crucially, the equations are well-posed at all inertial numbers, which corresponds to the full range of inclination angles for steady chute flow, in contrast with the limited range of angles in which the incompressible theory is well-posed [9]. At the same time, the new theory captures the dilatant behaviour of over-consolidated granular material as the structure of the equations is motivated by CSSM.

Incidentally, the chute flow calculation relates to Forterre & Pouliquen [33], who modify incompressible μ(I)-rheology by postulating a constitutive law in which ϕ is explicitly specified as a function of the inertial number, say Φ(I). In the CIDR model, there is no constitutive law relating these two variables; nevertheless, dependence of ϕ on I is implicit in the solution. It is also quite possible that a more subtle dependence of volume fraction on inertial number will be needed to reproduce some phenomena reported recently in experiments and numerical simulations, while retaining the property of well-posedness.

Our primary viewpoint in this paper has been to regard the CIDR model as modifying μ(I)-rheology with compressibility. However, it is equally valid to regard CIDR as modifying CSSM with rate dependence. To recapitulate, our result shows well-posed equations result from such modification provided that the yield locus and flow rule satisfy (2.19) and (2.20). Unlike in the incompressible μ(I)-rheology and the rate-independent CSSM equations, the CIDR equations are linearly well-posed for all deformations and for perturbations in all directions in Fourier space.

An important next step will be to specify constitutive laws, satisfying the general conditions for linear well-posedness of this paper, and formulated to accurately match the available experimental results and discrete numerical simulations for granular flows such as two-dimensional steady flow in a Couette geometry and time-dependent chute flow. Then the CIDR model and well-posedness result can be tested with fully two-dimensional nonlinear transient numerical computations of these flows.