3.1 Base state and non-dimensionalisation

Appropriate non-dimensionalisation of our equations is motivated by examining the base state of the system, which is easily obtained by seeking a unidirectional steady flow $\boldsymbol{u}=(u_{0}(z),0)$ of constant depth $h=h_{0}$ and solid fraction $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{0}(z)$ (Barker et al. Reference Barker, Schaeffer, Shearer and Gray2017). As this flow is incompressible, we have $f=0$ , and hence it follows from conservation of momentum and (2.5) that $\unicode[STIX]{x03C4}=p\tan \unicode[STIX]{x1D703}=\unicode[STIX]{x1D707}(I)p$ . Thus

(3.3) $$\begin{eqnarray}\displaystyle I=I_{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D707}^{-1}(\tan \unicode[STIX]{x1D703}) & & \displaystyle\end{eqnarray}$$

is constant throughout the flow. To ensure that $I$ is positive and bounded, the inclination must exist in the range $\tan ^{-1}\unicode[STIX]{x1D707}_{1}<\unicode[STIX]{x1D703}<\tan ^{-1}\unicode[STIX]{x1D707}_{2}$ , i.e. the range for which steady uniform chute flow is obtained by Forterre & Pouliquen (Reference Forterre and Pouliquen2003). Negative and unbounded $I$ represent transitions to the quasi-static and gaseous regimes, respectively. From conservation of momentum and the definition of the inertial number, we have that the base state pressure and velocity profiles (assuming that $u$ increases with $z$ (Gray & Edwards Reference Gray and Edwards2014)) are determined by the ordinary differential equations

(3.4a,b ) $$\begin{eqnarray}\displaystyle \frac{\text{d}p_{0}}{\text{d}z}=-\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{0}g\cos \unicode[STIX]{x1D703},\quad p_{0}(h_{0})=0,\quad \frac{\text{d}u_{0}}{\text{d}z}=\frac{I_{\unicode[STIX]{x1D703}}}{d}\sqrt{\frac{p_{0}}{\unicode[STIX]{x1D70C}_{s}}},\quad u_{0}(0)=0, & & \displaystyle\end{eqnarray}$$

where from (2.8) and (2.10) $\unicode[STIX]{x1D719}_{0}$ is given by

(3.5a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{0}=\unicode[STIX]{x1D719}_{min}+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D719}\left(\displaystyle \frac{2p_{0}b}{\unicode[STIX]{x1D6FD}(I_{\unicode[STIX]{x1D703}})\unicode[STIX]{x1D70C}_{s}gd}\right)}{1+\displaystyle \frac{2p_{0}b}{\unicode[STIX]{x1D6FD}(I_{\unicode[STIX]{x1D703}})\unicode[STIX]{x1D70C}_{s}gd}},\quad \unicode[STIX]{x0394}\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{max}-\unicode[STIX]{x1D719}_{min}. & & \displaystyle\end{eqnarray}$$

Thus (3.4) suggests we non-dimensionalise the model using the natural scales

(3.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle x,z,h\sim h_{0},\quad t\sim \frac{h_{0}}{U},\quad \boldsymbol{u}\sim U=\frac{I_{\unicode[STIX]{x1D703}}\sqrt{gh_{0}^{3}\cos \unicode[STIX]{x1D703}}}{d}, & \displaystyle\end{eqnarray}$$

(3.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle p,\unicode[STIX]{x03C4}_{ij},Y\sim \unicode[STIX]{x1D70C}_{s}gh_{0}\cos \unicode[STIX]{x1D703},\quad \unicode[STIX]{x1D60B}_{ij}\sim \frac{U}{h_{0}}. & \displaystyle\end{eqnarray}$$

(3.7a,b ) $$\begin{eqnarray}\displaystyle \frac{\text{d}p_{0}}{\text{d}z}=-\unicode[STIX]{x1D719}_{0},\quad p_{0}(1)=0,\quad \frac{\text{d}u_{0}}{\text{d}z}=\sqrt{p_{0}},\quad u_{0}(0)=0, & & \displaystyle\end{eqnarray}$$

while $\unicode[STIX]{x1D719}_{0}$ is given by

(3.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{0}=\unicode[STIX]{x1D719}_{min}+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D719}(2p_{0}\unicode[STIX]{x1D712})}{1+2p_{0}\unicode[STIX]{x1D712}},\quad \unicode[STIX]{x1D712}=\frac{bh_{0}\cos \unicode[STIX]{x1D703}}{\unicode[STIX]{x1D6FD}(I_{\unicode[STIX]{x1D703}})d}, & & \displaystyle\end{eqnarray}$$

and $\unicode[STIX]{x1D712}$ is a non-dimensional parameter. Note that one recovers a uniform volume fraction profile, as predicted from the incompressible rheology, in the limit $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}\rightarrow 0$ . The ordinary differential equation system (3.7) can readily be solved numerically (using ode45 in Matlab) to give the base state profiles illustrated in figure 2 for various values of $\unicode[STIX]{x1D712}$ . We recover an approximately lithostatic pressure distribution (figure 2 a), while figure 2(b) demonstrates the Bagnold velocity profile for granular chute flow (MiDi Reference Midi2004).

Figure 2.

Non-dimensional base state (a) pressure, (b) velocity and (c) solid fraction profiles as determined by (3.7) and (3.8) for parameters as in table 1.

As discussed by Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017), the parameter $\unicode[STIX]{x1D712}$ essentially controls the variation in the base state solid fraction profile over a given $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ range, and hence its mean value (figure 2 c). Experimental evidence illustrates that for steady uniform flow the solid fraction is approximately constant over the interior of the flow before rapidly decreasing near the free surface (see § 3), and thus this suggests the use of a large value for $\unicode[STIX]{x1D712}$ . While this parameter depends both on the ratio of the base state flow thickness to particle size $h_{0}/d$ ( ${\sim}10$ ) and inclination $\unicode[STIX]{x1D703}$ from (3.8), it is also proportional to the largely unconstrained parameter $b$ . In practice, it may be possible to determine this free parameter using suitable base state flow data (e.g. see § 4.2). In the case of inclined plane flow, provided knowledge of $h_{0}$ , $\unicode[STIX]{x1D703}$ and $d$ , one could fit the base state volume fraction profile (3.8) to available experimental data in order to yield an estimate for $b$ . However, in practice we anticipate that, in a similar manner to the parameters characterising the $\unicode[STIX]{x1D707}(I)$ phenomenological law (2.9c ), this parameter varies with the material properties of the granular medium and bed roughness conditions.

In the absence of additional information, we assume $\unicode[STIX]{x1D712}$ is a relatively large free parameter (e.g.  $25\leqslant \unicode[STIX]{x1D712}\leqslant 100$ ). Note that at various points throughout this work, results are illustrated for the representative value $\unicode[STIX]{x1D712}=75$ given that, in this case, (i) the base state $\unicode[STIX]{x1D719}$ profile mimics that found experimentally and (ii) we find $\overline{\unicode[STIX]{x1D719}_{0}}\approx 0.606$ , in agreement with the average solid fraction of $0.59\pm 0.03$ for steady chute flows as obtained by Pouliquen (Reference Pouliquen1999). Note that taking $\unicode[STIX]{x1D712}\rightarrow \infty$ introduces a discontinuity in the base state solid fraction profile given that $\unicode[STIX]{x1D719}_{0}$ must reduce to $\unicode[STIX]{x1D719}_{min}$ at the free surface (from (3.8)). This boundary layer will naturally introduce numerical difficulties when discretising $\unicode[STIX]{x1D719}_{0}$ (and derivatives thereof) in the solution of the generalised eigenvalue problem outlined in § 3.3. Thus we avoid taking $\unicode[STIX]{x1D712}$ too large ( ${>}100$ say) to circumvent such issues.

Finally, using the scaling (3.6), the non-dimensional CIDR system is given by

(3.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}\boldsymbol{u})=0, & \displaystyle\end{eqnarray}$$

(3.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle F^{2}\unicode[STIX]{x1D719}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}\right)=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\unicode[STIX]{x03C4}\frac{\unicode[STIX]{x1D63F}}{\Vert \unicode[STIX]{x1D63F}\Vert }\right)+\unicode[STIX]{x1D719}(\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D703}},-1), & \displaystyle\end{eqnarray}$$

(3.10a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=2f(p,\unicode[STIX]{x1D719},I)\Vert \unicode[STIX]{x1D63F}\Vert ,\quad \unicode[STIX]{x03C4}=Y(p,\unicode[STIX]{x1D719},I), & & \displaystyle\end{eqnarray}$$

where

(3.11a-c ) $$\begin{eqnarray}\displaystyle Y=\unicode[STIX]{x1D6FC}(I)p-\frac{\unicode[STIX]{x1D6FD}(I_{\unicode[STIX]{x1D703}})\unicode[STIX]{x1D712}p^{2}}{{\hat{C}}(\unicode[STIX]{x1D713})},\quad f=\unicode[STIX]{x1D6FD}(I)-\frac{2\unicode[STIX]{x1D6FD}(I_{\unicode[STIX]{x1D703}})\unicode[STIX]{x1D712}p}{{\hat{C}}(\unicode[STIX]{x1D713})},\quad I=\frac{2I_{\unicode[STIX]{x1D703}}\Vert \unicode[STIX]{x1D63F}\Vert }{\sqrt{p}}, & & \displaystyle\end{eqnarray}$$

with $I_{\unicode[STIX]{x1D703}}$ given by (3.3). Here we have defined the Froude number $F$ and constant $\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D703}}$ as

(3.12a,b ) $$\begin{eqnarray}\displaystyle F=\frac{U}{\sqrt{gh_{0}\cos \unicode[STIX]{x1D703}}},\quad \unicode[STIX]{x1D707}_{\unicode[STIX]{x1D703}}=\tan \unicode[STIX]{x1D703}. & & \displaystyle\end{eqnarray}$$

In addition, the surface boundary conditions (3.1) and (3.2) become

(3.13) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\displaystyle ph_{x}-\unicode[STIX]{x03C4}_{11}h_{x}+\unicode[STIX]{x03C4}_{12}=0,\quad -\unicode[STIX]{x03C4}_{12}h_{x}-p+\unicode[STIX]{x03C4}_{22}=0,\\ \displaystyle v=h_{t}+uh_{x},\end{array}\right\}\quad \text{at}~z=h(x,t), & & \displaystyle\end{eqnarray}$$

respectively, supplemented with the no-slip condition $\boldsymbol{u}=0$ at $z=0$ .

3.2 Linear stability analysis

The linear stability of the base state outlined in § 3.1 is analysed using a standard normal mode analysis (e.g. Drazin & Reid Reference Drazin and Reid1981, chap. 4). We consider perturbations of the form

(3.14) $$\begin{eqnarray}\displaystyle [u,v,\unicode[STIX]{x1D719},p,h]=[u_{0}(z),0,\unicode[STIX]{x1D719}_{0}(z),p_{0}(z),1]+[\tilde{u} ,\tilde{v},\tilde{\unicode[STIX]{x1D719}},\tilde{p},\tilde{h}](x,z,t), & & \displaystyle\end{eqnarray}$$

where perturbations are represented by the tilde variables. We then substitute these expansions into the governing equations (3.9), (3.10) and (3.13) and linearise the resulting system. After some tedious but nonetheless relatively straightforward algebra, we obtain

(3.15a ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D719}}_{t}+u_{0}\tilde{\unicode[STIX]{x1D719}}_{x}+\unicode[STIX]{x1D719}_{0}\tilde{u} _{x}+(\unicode[STIX]{x1D719}_{0}\tilde{v})_{z}=0, & \displaystyle\end{eqnarray}$$

(3.15b ) $$\begin{eqnarray}\displaystyle & \displaystyle F^{2}\unicode[STIX]{x1D719}_{0}(\tilde{u} _{t}+u_{0}\tilde{u} _{x}+u_{0}^{\prime }\tilde{v})=-\tilde{p}_{x}+\unicode[STIX]{x1D702}_{0}(2\tilde{u} _{x}/3-\tilde{v}_{z}/3)_{x}+[A_{p}\tilde{p}+A_{\unicode[STIX]{x1D719}}\tilde{\unicode[STIX]{x1D719}}+A_{u}(\tilde{u} _{z}+\tilde{v}_{x})]_{z}+\tilde{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D703}}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.15c ) $$\begin{eqnarray}\displaystyle & \displaystyle F^{2}\unicode[STIX]{x1D719}_{0}(\tilde{v}_{t}+u_{0}\tilde{v}_{x})=-\tilde{p}_{z}+(\unicode[STIX]{x1D702}_{0}[2\tilde{v}_{z}/3-\tilde{u} _{x}/3])_{z}+A_{p}\tilde{p}_{x}+A_{\unicode[STIX]{x1D719}}\tilde{\unicode[STIX]{x1D719}}_{x}+A_{u}(\tilde{u} _{z}+\tilde{v}_{x})_{x}-\tilde{\unicode[STIX]{x1D719}}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.15d ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{u} _{x}+\tilde{v}_{z}=C_{p}\tilde{p}+C_{\unicode[STIX]{x1D719}}\tilde{\unicode[STIX]{x1D719}}+C_{u}(\tilde{u} _{z}+\tilde{v}_{x}). & \displaystyle\end{eqnarray}$$

$z$

$u_{0}^{\prime }$

$(~)_{x,z,t}$

(3.16a ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{u} =\tilde{v}=0\quad \text{at}~z=0, & \displaystyle\end{eqnarray}$$

(3.16b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}l@{}}\displaystyle \unicode[STIX]{x1D707}_{\unicode[STIX]{x1D703}}p_{0}^{\prime }\tilde{h}+A_{p}\tilde{p}+A_{\unicode[STIX]{x1D719}}\tilde{\unicode[STIX]{x1D719}}+A_{u}(\tilde{u} _{z}+\tilde{v}_{x})=0,\\ \displaystyle -p_{0}^{\prime }\tilde{h}-\tilde{p}+\unicode[STIX]{x1D702}_{0}(2\tilde{v}_{z}/3-\tilde{u} _{x}/3)=0,\\ \displaystyle \tilde{v}=\tilde{h}_{t}+u_{0}\tilde{h}_{x},\end{array}\right\}\quad \text{at}~z=1. & \displaystyle\end{eqnarray}$$

$A$

$C$

$\unicode[STIX]{x1D702}_{0}$

$z$

$\tilde{h}$

$A_{u}$

$\unicode[STIX]{x1D702}_{0}$

$\tilde{u} _{z}$

$\tilde{v}_{z}$

$A_{u}\tilde{u} _{z}$

$\unicode[STIX]{x1D702}_{0}\tilde{v}_{z}$

$\widetilde{\unicode[STIX]{x1D6E4}}$

$\widetilde{\unicode[STIX]{x1D6F9}}$

(3.17a,b ) $$\begin{eqnarray}\displaystyle \widetilde{\unicode[STIX]{x1D6E4}}=A_{u}(\tilde{u} _{z}+\tilde{v}_{x}),\quad \widetilde{\unicode[STIX]{x1D6F9}}=\unicode[STIX]{x1D702}_{0}\left(\frac{2\tilde{v}_{z}}{3}-\frac{\tilde{u} _{x}}{3}\right), & & \displaystyle\end{eqnarray}$$

so that the boundary conditions are now regular in these variables. One should note, however, that (3.16b ) differs from the equivalent linearised boundary condition in the IIDR model (see Forterre Reference Forterre2006, (2.12)) due to the presence of the additional $\widetilde{\unicode[STIX]{x1D6F9}}$ stress term. Incompressibility implies that this stress term should vanish at the free surface; however, by consideration of the linearised flow rule (3.15d ), it appears there is no reason a priori to assume $\widetilde{\unicode[STIX]{x1D6F9}}=0$ in the compressible model. As such, we have chosen to include this term in the boundary conditions.

Note that (3.15) is a variable-coefficient linear system whose coefficients are independent of both time and the downstream direction. Therefore, following standard practice in the literature (see e.g. Drazin & Reid (Reference Drazin and Reid1981, chap. 4) and Chandrasekhar (Reference Chandrasekhar1981, chap. 1)), we seek normal mode solution

(3.18) $$\begin{eqnarray}\displaystyle [\tilde{u} ,\tilde{v},\tilde{\unicode[STIX]{x1D719}},\tilde{p}]=\text{e}^{\text{i}kx+\unicode[STIX]{x1D70E}t}[\hat{u} ,\hat{v},\hat{\unicode[STIX]{x1D719}},\hat{p}](z), & & \displaystyle\end{eqnarray}$$

where $\hat{u} =\hat{u} (z)$ , and similarly so for the other eigenfunctions (Malik Reference Malik1990; Alam & Nott Reference Alam and Nott1998; Forterre Reference Forterre2006). While a perturbation solution in terms of Kelvin modes may be suitable in the case of unbounded flow (Alam & Nott Reference Alam and Nott1997; Goddard & Lee Reference Goddard and Lee2017), the imposition of vertical boundary conditions in our case restricts the perturbation solution to the form (3.18) (Wang, Jackson & Sundaresan Reference Wang, Jackson and Sundaresan1996). Here $k\in \mathbb{R}$ is the downstream wavenumber and $c=\unicode[STIX]{x1D70E}/k$ , where $\unicode[STIX]{x1D70E}\in \mathbb{C}$ , is the wave speed. Thus for a given wavenumber the linear evolution of the system is governed by the growth rate $\text{Re}(\unicode[STIX]{x1D70E})$ and angular frequency $-\text{Im}(\unicode[STIX]{x1D70E})$ . This is referred to as a temporal stability analysis, which contrasts with a spatial stability analysis in which a pulsation $\unicode[STIX]{x1D70E}=-\text{i}\unicode[STIX]{x1D714}$ for $\unicode[STIX]{x1D714}\in \mathbb{R}$ is specified and $k\in \mathbb{C}$ is the eigenvalue of interest (see Schmid & Henningson Reference Schmid and Henningson2001, chap. 7). A spatial stability analysis is carried out in Forterre (Reference Forterre2006) and Gray & Edwards (Reference Gray and Edwards2014) in order to compare their results directly with the spatial stability data of Forterre & Pouliquen (Reference Forterre and Pouliquen2003); however, we avoid this approach given that the eigenvalue $k$ appears quadratically in the stability problem, which is typically more difficult to handle numerically (Malik Reference Malik1990; Varsakelis & Papalexandris Reference Varsakelis and Papalexandris2016). In any case, we discuss in § 3.4 how the temporal analysis can be used to compare against the spatial stability data.

Substitution of (3.18) into the linearised equations yields a system of four linear differential–algebraic equations (see (A 3) and (A 4) in appendix A) for the velocity, pressure and solid fraction perturbations, forming a generalised eigenvalue problem in $\unicode[STIX]{x1D70E}$ . Note that in the incompressible $\unicode[STIX]{x1D707}(I)$ limit, i.e. by taking $f=0$ , $\unicode[STIX]{x1D719}_{0}=\text{const.}$ , $\tilde{\unicode[STIX]{x1D719}}=0$ and $Y=\unicode[STIX]{x1D707}(I)p$ , equations (3.15) and (3.16) reduce to the system presented in Forterre (Reference Forterre2006), as expected.

3.3 Numerical method for stability problem

Following standard practice in the hydrodynamic stability of compressible flows (e.g. Malik Reference Malik1990; Malik, Alam & Dey Reference Malik, Alam and Dey2006), we implement a Chebyshev collocation method (see Trefethen Reference Trefethen2000; Boyd Reference Boyd2001, chap. 6) to solve the stability problem (A 3) subject to the boundary conditions (A 4). We propose a series expansion of a given eigenfunction $Q(z)$ in terms of the Chebyshev polynomials $T_{n}(\unicode[STIX]{x1D709})=\cos (n\cos ^{-1}\unicode[STIX]{x1D709})$ , i.e. 

(3.19) $$\begin{eqnarray}\displaystyle Q(z)=\mathop{\sum }_{n=0}^{N}a_{n}^{Q}T_{n}(\unicode[STIX]{x1D709}), & & \displaystyle\end{eqnarray}$$

where $N$ is the degree of the expansion, $a_{n}^{Q}$ are the expansion coefficients for eigenfunction $Q$ and we define the change of variables $\unicode[STIX]{x1D709}=2z-1$ to map the physical domain $z\in [0,1]$ to the Chebyshev domain $\unicode[STIX]{x1D709}\in [-1,1]$ . The governing equations are then collocated at the Gauss–Lobatto points

(3.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D709}_{i}^{{\mathcal{G}}{\mathcal{L}}}=\cos \left(\frac{\unicode[STIX]{x03C0}i}{N}\right),\quad i=0,\ldots ,N, & & \displaystyle\end{eqnarray}$$

i.e. the extrema of the polynomial $T_{n}$ , while evaluation at the boundary nodes is replaced by the relevant boundary conditions.

However, discretisation of the stability problem in this manner requires the use of artificial boundary conditions for the solid fraction perturbation. It is well known that the numerical solution of generalised eigenvalue problems with artificial boundary conditions can suffer from spurious (i.e. unphysical) eigenvalues (e.g. see Alam & Nott (Reference Alam and Nott1998) and Boyd (Reference Boyd2001, chap. 7)). To avoid the specification of such boundary conditions, we follow the approach of Müller & Kleiser (Reference Müller and Kleiser2008) by implementing a staggered collocation method. In particular, the momentum equations are evaluated on the Gauss–Lobatto points $\unicode[STIX]{x1D709}_{i}^{{\mathcal{G}}{\mathcal{L}}}$ , while conservation of mass, the flow rule and the definitions of the auxiliary variables (3.17) are collocated on the Gauss points

(3.21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D709}_{i}^{{\mathcal{G}}}=\cos \left[\frac{(2i+1)\unicode[STIX]{x03C0}}{2N}\right],\quad i=0,\ldots ,N-1. & & \displaystyle\end{eqnarray}$$

Thus the velocity eigenfunctions $(\hat{u} ,\hat{v})$ are expanded as $N+1$ term series, while the other eigenfunctions $(\hat{\unicode[STIX]{x1D719}},\hat{p},\hat{\unicode[STIX]{x1D6E4}},\hat{\unicode[STIX]{x1D6F9}})$ are of one degree less. Evaluation of the Chebyshev polynomials (and derivatives thereof) at the Gauss–Lobatto and Gauss points is carried out using the formulae detailed in appendix A.2 of Boyd (Reference Boyd2001).

Upon discretisation of the stability problem in this manner, we obtain the system

(3.22) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63C}\boldsymbol{Q}=\unicode[STIX]{x1D63D}\unicode[STIX]{x1D70E}\boldsymbol{Q}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D63C}$ is a $(6N+2)^{2}$ dense matrix, $\boldsymbol{Q}$ is the vector of unknown coefficients and $\unicode[STIX]{x1D63D}$ is a rank-deficient matrix, which results in the system (3.22) possessing a large number of infinite eigenvalues. Thus, we adopt the approach of Moyles & Wetton (Reference Moyles and Wetton2015) and consider instead

(3.23) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63D}\boldsymbol{Q}=\frac{1}{\unicode[STIX]{x1D70E}}\unicode[STIX]{x1D63C}\boldsymbol{Q}=\unicode[STIX]{x1D63C}\unicode[STIX]{x1D706}\boldsymbol{Q},\quad \unicode[STIX]{x1D706}=\frac{1}{\unicode[STIX]{x1D70E}}\;\Longrightarrow \;\unicode[STIX]{x1D70E}=\frac{\bar{\unicode[STIX]{x1D706}}}{|\unicode[STIX]{x1D706}|^{2}}, & & \displaystyle\end{eqnarray}$$

which is both more computationally efficient and accurate to solve. The eigenvalues $\unicode[STIX]{x1D706}$ and eigenvectors $\boldsymbol{Q}$ can then be obtained using the standard eig solver in Matlab. For linear stability we are only interested in the value of $\unicode[STIX]{x1D70E}$ with maximal real part, and henceforth we assume that $\unicode[STIX]{x1D70E}$ represents this maximal value.

3.4 Incompressible limit

Our stability analysis of the CIDR model represents a generalisation of the IIDR analysis presented in Forterre (Reference Forterre2006), and thus it is natural to consider an incompressible limit of our eigenvalue problem as a means of validation of the numerical method. Given that variability in the solid fraction is characterised by the quantity $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{max}-\unicode[STIX]{x1D719}_{min}$ , we can consider an incompressible limit by taking $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}\rightarrow 0$ . Indeed, Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) demonstrate that the CIDR model suitably reduces to its incompressible equivalent in this limit. To make a sensible comparison between our results and those of Forterre (Reference Forterre2006), we must take into account a number of important differences between the two studies. The first simply relates to the choice of scales and the subsequent definition of the Froude number $F$ , which in this work is given by (3.12), i.e.  $F=I_{\unicode[STIX]{x1D703}}h_{0}/d$ . However since the solid fraction is assumed constant in Forterre (Reference Forterre2006), which we denote $\unicode[STIX]{x1D719}_{Bulk}$ , the authors incorporate it into their scaling of the model. As such, we find that their Froude number, ${\mathcal{F}}$ , relates to ours via

(3.24) $$\begin{eqnarray}\displaystyle {\mathcal{F}}=\frac{2\sqrt{\unicode[STIX]{x1D719}_{Bulk}}}{5}F, & & \displaystyle\end{eqnarray}$$

where, for simplicity, we assume that $\unicode[STIX]{x1D719}_{Bulk}$ is given by the mean value of our base state solid fraction profile, i.e.  $\unicode[STIX]{x1D719}_{Bulk}=\overline{\unicode[STIX]{x1D719}_{0}}$ .

The second key difference relates to the spatial stability analysis in which perturbations to the base state are proportional to $\text{e}^{\text{i}k^{\ast }x-\text{i}\unicode[STIX]{x1D714}^{\ast }t}$ , where $k^{\ast }\in \mathbb{C}$ , $\unicode[STIX]{x1D714}^{\ast }\in \mathbb{R}$ and asterisked quantities refer to those used by Forterre (Reference Forterre2006). Thus the growth rate and wave speed are given by $\unicode[STIX]{x1D70E}^{\ast }=-\text{Im}(k^{\ast })$ and $c^{\ast }=\unicode[STIX]{x1D714}^{\ast }/\text{Re}(k^{\ast })$ , respectively. In order to compare our results against the spatial dispersion relations, we perform an equivalent analysis by seeding the eigenvalue solver with imaginary values for the wavenumber and searching over complex $k$ space. The set of eigenvalues which have real part zero can thus be interpreted as values for the perturbation frequency $\unicode[STIX]{x1D714}=-\text{Im}(\unicode[STIX]{x1D70E})$ . Hence from this set we can obtain values for the growth rate $-\text{Im}(k)$ and wave speed $\unicode[STIX]{x1D714}/\text{Re}(k)$ , allowing for direct comparison with the results of Forterre (Reference Forterre2006).

Figure 3.

Spatial dispersion relations for the dimensionless (a) growth rate and (b) wave speed as a function of the perturbation frequency. The solid line indicates the predicted profile from the solution of the full compressible problem with $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=10^{-3}$ , while the markers represent the equivalent prediction from the IIDR model of Forterre (Reference Forterre2006). Parameter values as in table 1 with $\unicode[STIX]{x1D719}_{min}=0.649$ , $N=100$ , $\unicode[STIX]{x1D703}=29^{\circ }$ , ${\mathcal{F}}=1.02$ and $\unicode[STIX]{x1D712}=75$ .

Such a comparison is illustrated in figure 3 where the leading growth rate and corresponding wave speed are plotted as a function of the frequency $\unicode[STIX]{x1D714}$ . Note that the time scales in each study differ by a factor of $5/(2\sqrt{\unicode[STIX]{x1D719}_{Bulk}})$ , and hence to make a practical comparison between the two, we rescale our computed eigenvalue (and thus $\unicode[STIX]{x1D714}$ ) by this factor. For $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=10^{-3}$ we observe excellent agreement between the spatial dispersion relation computed using our numerical method and that given by Forterre (Reference Forterre2006) (see figure 2 therein), with a slight discrepancy between the two at larger values of $\unicode[STIX]{x1D714}$ (corresponding to larger $\text{Re}(k)$ ). Such differences at large $k$ are anticipated given that the incompressible rheology is linearly ill-posed. The overall agreement in figure 3 provides strong evidence for the veracity of the solver. Further validation by comparison with a long-wavelength asymptotic analysis of the stability problem is presented in the supplementary material available at https://doi.org/10.1017/jfm.2019.43.

4.1 Trulsson et al. (Reference Trulsson, Bouzid, Claudin and Andreotti2013) model

While conservation of mass and momentum (2.1)–(2.2) form the basis for both the Trulsson et al. (Reference Trulsson, Bouzid, Claudin and Andreotti2013) and CIDR models, several key differences exist in their constitutive relations. Firstly, in lieu of the flow rule (2.7) the authors invoke a linear dilatancy law observed in steady plane shear (Da Cruz et al. Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005; Pouliquen et al. Reference Pouliquen, Cassar, Jop, Forterre and Nicolas2006; Trulsson et al. Reference Trulsson, Andreotti and Claudin2012), relating the volume fraction to the inertial number, i.e. 

(4.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{c}-BI, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{c}$ is the critical volume fraction at the jamming point and $B$ is some constant. By combining this relation with the definition of the inertial number (2.6), the authors express the pressure $p$ as a function of the volume fraction and strain rate, i.e. 

(4.3a,b ) $$\begin{eqnarray}\displaystyle p=\unicode[STIX]{x1D70C}_{s}d^{2}\Vert \unicode[STIX]{x1D738}\Vert ^{2}\mathfrak{f}(\unicode[STIX]{x1D719}),\quad \mathfrak{f}(\unicode[STIX]{x1D719})=\left(\frac{B}{\unicode[STIX]{x1D719}_{c}-\unicode[STIX]{x1D719}}\right)^{2},\quad \unicode[STIX]{x03B3}_{ij}=\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x03B3}_{ij}$ is the strain rate tensor. Note that this relation is obtained under the assumption that the flow is incompressible, i.e.  $\unicode[STIX]{x1D60B}_{ij}=(1/2)\unicode[STIX]{x03B3}_{ij}$ , evidently violated in unsteady flow. A more appropriate generalisation, in relation to the CIDR model, would be to use the norm of the deviatoric strain rate tensor $\Vert \unicode[STIX]{x1D63F}\Vert$ instead of $\Vert \unicode[STIX]{x1D738}\Vert$ in (4.3). In addition, using dimensional arguments, the authors propose that the stress tensor can be decomposed as

(4.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x03C3}_{ij}=-p\unicode[STIX]{x03B4}_{ij}+\frac{\unicode[STIX]{x1D707}(\unicode[STIX]{x1D719})p}{\Vert \unicode[STIX]{x1D738}\Vert }\unicode[STIX]{x03B3}_{ij}+\frac{2g(\unicode[STIX]{x1D719})p}{\Vert \unicode[STIX]{x1D738}\Vert }(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u})\unicode[STIX]{x03B4}_{ij}, & & \displaystyle\end{eqnarray}$$

where $g(\unicode[STIX]{x1D719})$ is some unknown function and $\unicode[STIX]{x1D707}(\unicode[STIX]{x1D719})$ is that given by (2.9c ), under the proviso that $I=I(\unicode[STIX]{x1D719})$ from (4.2). In practice the authors assume $g(\unicode[STIX]{x1D719})$ constant, although it is largely unconstrained. While (4.4) bears some similarity to the deviatoric stress–strain rate relation (2.3)–(2.4), there are a number of crucial differences. Indeed, the two can only be cast in a similar form under the modifications that (i) the strain rate norm is replaced by the deviatoric strain rate norm and (ii) one assumes $g(\unicode[STIX]{x1D719})=-(1/3)\unicode[STIX]{x1D707}(\unicode[STIX]{x1D719})$ . While a comparison between the merits and limitations of the Trulsson et al. (Reference Trulsson, Bouzid, Claudin and Andreotti2013) and CIDR models is beyond the scope of this work, we note that the former (i) does not investigate the well-posedness of the model, (ii) lacks exposition of the stress tensor formulation and (iii) invokes an additional phenomenological law in (4.2), as opposed to the more physically motivated flow rule (2.7).

Specification of the rheological and constitutive parameters in the CIDR model for this problem can be achieved by making use of those detailed in the model above. In particular, given that both models utilise the same functional form for $\unicode[STIX]{x1D707}(I)$ , we assume equivalent values for $\unicode[STIX]{x1D707}_{1,2}$ and $I_{0}$ as determined by Trulsson et al. (Reference Trulsson, Andreotti and Claudin2012) (see table 3). Note that $\unicode[STIX]{x1D719}_{c}$ in (4.2) is precisely the maximum volume fraction $\unicode[STIX]{x1D719}_{max}$ in (2.10), while we can also take $B=\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ (e.g. see Pouliquen et al. Reference Pouliquen, Cassar, Jop, Forterre and Nicolas2006; Forterre & Pouliquen Reference Forterre and Pouliquen2008) to yield a suitable estimate for $\unicode[STIX]{x1D719}_{min}$ .

Table 3.

Parameter values assumed for the case of plane shear flow in the absence of gravity (see Trulsson et al. Reference Trulsson, Andreotti and Claudin2012, Reference Trulsson, Bouzid, Claudin and Andreotti2013).

4.2 Base state

In order to compare the CIDR predictions to the plane shear data available, it is more convenient to present our analysis using the dimensional model (2.1)–(2.10). The base state of the system again follows from assuming a unidirectional flow of the form $\boldsymbol{u}=(u_{0}(z),0)$ with constant volume fraction $\unicode[STIX]{x1D719}_{0}$ , whose value is imposed in the discrete element simulations. From conservation of momentum in the downstream and normal directions, we have that the base state pressure, $p_{0}$ , and shear stress, $Y_{0}$ , are constant. As the base state flow is incompressible, it follows from (2.8) and (2.9) that $I=I_{c}=\unicode[STIX]{x1D707}^{-1}(Y_{0}/p_{0})$ is also constant. Thus, assuming a positive velocity gradient, the base state downstream velocity follows from the definition of the inertial number (2.6), i.e. 

(4.5a,b ) $$\begin{eqnarray}\displaystyle \frac{\text{d}u_{0}}{\text{d}z}=\frac{I_{c}}{d}\sqrt{\frac{p_{0}}{\unicode[STIX]{x1D70C}_{s}}}\;\Longrightarrow \;u_{0}=\dot{\unicode[STIX]{x1D6FE}}z,\quad \dot{\unicode[STIX]{x1D6FE}}=\frac{I_{c}}{d}\sqrt{\frac{p_{0}}{\unicode[STIX]{x1D70C}_{s}}}, & & \displaystyle\end{eqnarray}$$

where the stationary plate is located at $z=0$ and the velocity gradient $\dot{\unicode[STIX]{x1D6FE}}$ is constant. Finally, given that $f=0$ for this flow, it follows from (2.8) that the base state must satisfy

(4.6) $$\begin{eqnarray}\displaystyle C(\unicode[STIX]{x1D719}_{0})=\frac{2p_{0}}{\unicode[STIX]{x1D6FD}(I_{c})}. & & \displaystyle\end{eqnarray}$$

As gravity is omitted in this problem, it is more suitable to scale $C(\unicode[STIX]{x1D719})$ in (2.10) with a characteristic shear stress rather than the hydrostatic scaling used previously (see § 3). A balance of the inertial and pressure terms in (2.2) suggests ${\mathcal{P}}=\unicode[STIX]{x1D70C}_{s}U^{2}$ as an appropriate scale, assuming that $\boldsymbol{u}$ scales with the wall velocity, where $U=\dot{\unicode[STIX]{x1D6FE}}H$ follows from (4.5). Note that this pressure scale follows that previously utilised in Goddard & Lee (Reference Goddard and Lee2017) upon the rescaling $(d/H)^{2}$ . Indeed it may be more natural to use such a scale, which is independent of external accelerations, so as to avoid problems regarding frame-indifference of the model, as discussed in § 3. In this case, note that one can use (4.6) to determine the free parameter $b$ in (2.10) as

(4.7) $$\begin{eqnarray}\displaystyle b=\left(\frac{H}{d}\right)^{2}\frac{\unicode[STIX]{x1D6FD}(I_{c})I_{c}^{2}}{2}{\hat{C}}(\unicode[STIX]{x1D719}_{0}), & & \displaystyle\end{eqnarray}$$

given knowledge of the geometry and the base state volume fraction and inertial number. Thus provided $b$ can be determined for the problem at hand, no free parameter exists.

4.3 Linear stability analysis

We now perform a temporal stability analysis analogous to that presented in § 3.2. In this case, however, the calculation is significantly simplified given that (i) we follow Trulsson et al. (Reference Trulsson, Bouzid, Claudin and Andreotti2013) by assuming $x$ -independence for all variables and (ii) the base state pressure, volume fraction and velocity gradient are all constant. As such, perturbations to the base state can be expressed in terms of plane wave solutions, i.e. 

(4.8) $$\begin{eqnarray}\displaystyle [u,v,\unicode[STIX]{x1D719},p]=[u_{0}(z),0,\unicode[STIX]{x1D719}_{0},p_{0}]+\text{e}^{\text{i}Kz+\text{i}\unicode[STIX]{x1D714}t}[\hat{u} ,\hat{v},\hat{\unicode[STIX]{x1D719}},\hat{p}], & & \displaystyle\end{eqnarray}$$

where $K$ is the transverse wavenumber and we use $\unicode[STIX]{x1D70E}=\text{i}\unicode[STIX]{x1D714}$ as the eigenvalue to ensure consistency with Trulsson et al. (Reference Trulsson, Bouzid, Claudin and Andreotti2013). Upon substitution of (4.8), linearisation of (2.1), (2.2) and (2.7) yields the (dimensional) system

(4.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}\unicode[STIX]{x1D714}\hat{\unicode[STIX]{x1D719}}+\text{i}\unicode[STIX]{x1D719}_{0}K\hat{v}=0, & \displaystyle\end{eqnarray}$$

(4.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{0}(\text{i}\unicode[STIX]{x1D714}\hat{u} +\dot{\unicode[STIX]{x1D6FE}}\hat{v})=\text{i}K(D_{p}\hat{p}+D_{\unicode[STIX]{x1D719}}\hat{\unicode[STIX]{x1D719}}+\text{i}KD_{u}\hat{u} ), & \displaystyle\end{eqnarray}$$

(4.9c ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}_{0}\unicode[STIX]{x1D714}\hat{v}=-\text{i}K\hat{p}+\unicode[STIX]{x1D702}_{c}(\text{i}K)^{2}\hat{v}, & \displaystyle\end{eqnarray}$$

(4.9d ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{i}K\hat{v}=E_{p}\hat{p}+E_{\unicode[STIX]{x1D719}}\hat{\unicode[STIX]{x1D719}}+\text{i}KE_{u}\hat{u} , & \displaystyle\end{eqnarray}$$

$D$

$E$

$\unicode[STIX]{x1D702}_{c}$

$\unicode[STIX]{x1D714}$

$K$