2.1 Constitutive expressions and closure relations

In order to close (2.1)–(2.3), suitable expressions for the terms $\hat{p}^{v}$ , $\hat{{\bf\tau}}$ , $\hat{{\it\beta}}$ and $\hat{{\it\Gamma}}$ have to be employed. For the irreversible isotropic pressure $\hat{p}^{v}$ and the deviator $\hat{{\bf\tau}}$ we opt for the constitutive expressions derived from the ${\it\mu}(I)$ -rheology. This is a phenomenological rheology law that combines a fluid-like, rate-dependent approach with a yield criterion. It is consistent with dimensional analysis and its predictions compare favourably against experimental measurements in a multitude of configurations; see, for example, MiDi (Reference MiDi2004), Jop et al. (Reference Jop, Forterre and Pouliquen2006) as well as the more recent review article of Forterre & Pouliquen (Reference Forterre and Pouliquen2008). Herein, we employ the three-dimensional extension of the ${\it\mu}(I)$ -rheology law of Jop et al. (Reference Jop, Forterre and Pouliquen2006) in the volume fraction representation. Then, the constitutive relations for $\hat{p}^{v}$ and $\hat{{\bf\tau}}$ read,

(2.5a,b ) $$\begin{eqnarray}\displaystyle \hat{{\it\tau}}=\frac{\hat{{\it\mu}}({\it\phi})\,|\hat{\dot{{\bf\gamma}}}|\,\hat{d}_{p}^{2}\hat{{\it\rho}}}{I^{2}}\hat{\dot{{\bf\gamma}}},\quad \hat{p}^{v}=\frac{|\hat{\dot{{\bf\gamma}}}|^{2}\,\hat{d}_{p}^{2}\hat{{\it\rho}}}{I^{2}}, & & \displaystyle\end{eqnarray}$$

(2.6a,b ) $$\begin{eqnarray}\displaystyle I=\frac{{\it\phi}_{max}-{\it\phi}}{{\it\phi}_{max}-{\it\phi}_{min}},\quad {\it\mu}({\it\phi})={\it\mu}_{1}+\frac{{\it\mu}_{2}}{\displaystyle \frac{I_{0}}{I}+1}. & & \displaystyle\end{eqnarray}$$

In the above equations, $\dot{{\bf\gamma}}=(\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\text{T}})-(2/3)\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}$ is the deviator of the strain-rate tensor whilst $|\hat{\dot{{\bf\gamma}}}|=\sqrt{(\hat{\dot{{\it\gamma}}}_{ij}\hat{\dot{{\it\gamma}}}_{ij})/2}$ designates its second invariant (written in the Einstein summation convention). Also, $I=|\hat{\dot{{\bf\gamma}}}|\hat{d}_{p}/\sqrt{(\hat{p}+\hat{p}^{v})/\hat{{\it\rho}}}$ is the inertial number, with $\hat{d}_{p}$ standing for the particle diameter. As illustrated by da Cruz et al. (Reference da Cruz, Prochnow, Roux and Chevoir2005), $I$ describes the ratio between the macroscopic time scale, induced by deformation due to shear, to the inertial time scale associated with the pressure force. Further, ${\it\phi}_{min}$ and ${\it\phi}_{max}$ stand for the minimum and maximum particle concentration whereas ${\it\mu}_{1}$ , ${\it\mu}_{2}$ and $I_{0}$ are constants whose value is acquired though experimental measurements.

According to (2.5) and (2.6), both $\hat{p}^{v}$ and $\hat{{\bf\tau}}$ diverge as ${\it\phi}\rightarrow {\it\phi}_{max}$ . This divergence describes the jamming transition that grains experience upon attaining their maximum concentration. In other words, in the volume fraction representation, the ${\it\mu}(I)$ -rheology gives rise to Krieger–Dougherty type laws for both the normal and shear viscosities. It is also interesting to observe that $I$ is a decreasing function of ${\it\phi}$ and is maximized at ${\it\phi}={\it\phi}_{min}$ with maximum value $I({\it\phi}_{min})=1$ . This inverse correlation describes the fact that an increased concentration should translate into larger frictional effects. Moreover, the value of unity is also not coincidental. The ${\it\mu}(I)$ -rheology predicts a shear-thickening behaviour through (i) the inverse dependence of the shear stresses on the inertial number $I$ , and (ii) the nonlinear dependence of the shear stresses on the shear rate. Therefore, when $I=1$ , shear thickening can only happen through the second term whilst for $I>1$ , the shear-thickening behaviour is not guaranteed since there is a competition between (i) and (ii). In this respect, the value $I=1$ is the borderline value upon which the behaviour of the ${\it\mu}(I)$ -rheology changes qualitatively.

At this point, a comment on the appropriateness of the employed rheology law is relevant. Strictly speaking, the ${\it\mu}(I)$ -rheology and its incarnations have been developed for and measured in steady granular flows. Consequently, the presumed validity of this law for transient flows constitutes more of an ad hoc postulate than a documented fact. Nevertheless, we stress that the ${\it\mu}(I)$ -rheology has been already exercised in the description of transient flows with a notable success. For instance, Lagrée, Staron & Popinet (Reference Lagrée, Staron and Popinet2011) combined the ${\it\mu}(I)$ -rheology with the two-phase Navier–Stokes equations for the investigation of a granular column collapse. The resulting numerical predictions accorded reasonably well with experimental data. Thus, on the premise that the limitations are acknowledged, the utilization of the ${\it\mu}(I)$ -rheology for the flows of interest can be justified.

It is also important to note that, as mentioned above, $\hat{p}^{v}$ is in general decomposed into a bulk viscous pressure and a shear-induced pressure. The ${\it\mu}(I)$ -rheology provides expressions only for the latter component since it based on the assumption that the flow is isochoric. A more complete constitutive expression for $\hat{p}^{v}$ would read,

(2.7) $$\begin{eqnarray}\displaystyle \hat{p}^{v}=\frac{|\hat{\dot{{\bf\gamma}}}|^{2}\hat{d}_{p}^{2}\hat{{\it\rho}}}{I^{2}}-\hat{{\it\zeta}}({\it\phi},|\hat{\dot{{\bf\gamma}}}|)\hat{\boldsymbol{{\rm\nabla}}}\boldsymbol{\cdot }\hat{\boldsymbol{u}}, & & \displaystyle\end{eqnarray}$$

where $\hat{{\it\zeta}}({\it\phi},|\hat{\dot{{\bf\gamma}}}|)$ is a bulk viscosity coefficient. However, since we are not acquainted with any systematic measurements for $\hat{{\it\zeta}}({\it\phi},|\hat{\dot{{\bf\gamma}}}|)$ , we will assume that, for the flows of interest,

(2.8) $$\begin{eqnarray}\displaystyle \hat{{\it\zeta}}({\it\phi},|\hat{\dot{{\bf\gamma}}}|)\hat{\boldsymbol{{\rm\nabla}}}\boldsymbol{\cdot }\hat{\boldsymbol{u}}\ll \min \left\{\frac{|\hat{\dot{{\bf\gamma}}}|^{2}\hat{d}_{p}^{2}\hat{{\it\rho}}}{I^{2}},\frac{\hat{{\it\mu}}({\it\phi})|\hat{\dot{{\bf\gamma}}}|\hat{d}_{p}^{2}\hat{{\it\rho}}}{I^{2}}\hat{\dot{{\bf\gamma}}}\right\}, & & \displaystyle\end{eqnarray}$$

which practically amounts to invoking a Stokes-type hypothesis.

For the intergranular stress $\hat{{\it\beta}}$ , we adopt the empirical correlation of Powers, Stewart & Krier (Reference Powers, Stewart and Krier1989) which is based on the experimental measurements of Elban & Chiarito (Reference Elban and Chiarito1986). This correlation assumes the following form,

(2.9) $$\begin{eqnarray}\displaystyle \hat{{\it\beta}}=\hat{k}_{1}\left(\frac{{\it\phi}}{2-{\it\phi}}\right)^{2}\log \left(\frac{1}{1-{\it\phi}}\right). & & \displaystyle\end{eqnarray}$$

Here, $\hat{k}_{1}$ is a strictly positive, material dependent, constant. Equation (2.9) is clearly a monotonically increasing, convex function of the volume fraction. As mentioned in Powers et al. (Reference Powers, Stewart and Krier1989), the monotonicity reflects the fact that, as the particle concentration and thus the intergranular stress increase, an increased hydrostatic pressure is required to offset $\hat{{\it\beta}}$ ; see also the earlier discussion of Carroll & Holt (Reference Carroll and Holt1972). Convexity, on the other hand, is welcomed from the mathematical point of view because it ameliorates technical difficulties associated with the establishment of the existence of solutions, (Varsakelis & Papalexandris Reference Varsakelis and Papalexandris2014a ).

As regards the equilibrated stress coefficient $\hat{{\it\Gamma}}$ , we utilize a ${\it\phi}$ -weighted variation of the functional relation proposed by Passman et al. (Reference Passman, Nunziato and Bailey1986),

(2.10) $$\begin{eqnarray}\displaystyle \hat{{\it\Gamma}}=\frac{\hat{k}_{2}{\it\phi}}{({\it\phi}_{max}-{\it\phi})^{2}}, & & \displaystyle\end{eqnarray}$$

where $\hat{k}$ is a material-dependent, positive constant. Equation (2.10) represents a form of Krieger–Dougherty type relation and, as such, diverges to infinity as ${\it\phi}\rightarrow {\it\phi}_{max}$ but additionally converges to zero as ${\it\phi}\rightarrow 0$ . Passman et al. (Reference Passman, Nunziato and Bailey1986), who based their selection on the experimental measurements of Savage (Reference Savage1979), argued that the divergence at the maximum concentration is consistent with the behaviour of (2.5) and (2.6) at this limit.

The last quantity that needs to be determined is the compaction viscosity $\hat{{\it\nu}}_{c}$ . This is somewhat problematic because, as remarked by Powers, Stewart & Krier (Reference Powers, Stewart and Krier1990) and much later by Lowe & Greenaway (Reference Lowe and Greenaway2005), reliable empirical correlations for $\hat{{\it\nu}}_{c}$ are not available in the literature. In the absence of any further guidance, we adopt the heuristic analysis of Baer & Nunziato (Reference Baer and Nunziato1986), and assume that $\hat{{\it\nu}}_{c}$ can be approximated by a constant.

2.2 Non-dimensionalisation and base-flow profile

In a plane, unidirectional Couette flow, the granular material is placed between two infinite parallel planes, in the $\hat{x}$ – ${\hat{y}}$ Cartesian plane, that are kept at a constant distance ${\hat{H}}$ (Couette gap). Further, the line ${\hat{y}}=0$ is assumed to coincide with the lower plane. The upper plane is presumed to move at a constant velocity $\hat{U} _{w}>0$ in the $\hat{x}$ -direction.

In the context of a linear stability analysis, it is advantageous to utilize the non-dimensional form of the governing equations. For this, we first determine the (dimensional) reference values that will be denoted with the subscript ‘ $r$ ’.

As is typically the case in the study of Couette flows, the reference velocity $\hat{u} _{r}$ is taken equal to the wall velocity $\hat{U} _{w}$ . However, due to the explicit presence of the particle diameter $\hat{d}_{p}$ in the expressions for the deviator and the irreversible isotropic pressure, there is an ambiguity in the choice of the reference length. Herein, we employ the Couette gap ${\hat{H}}$ as the reference length $\hat{L}_{r}$ . With these choices, velocity vectors, length vectors, and time are non-dimensionalised as follows,

(2.11a−c ) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{i}={\displaystyle \frac{\hat{\boldsymbol{u}}_{i}}{\hat{u} _{r}}},\quad \boldsymbol{x}={\displaystyle \frac{\hat{\boldsymbol{x}}}{\hat{L}_{r}}},\quad t=\hat{t}_{r}{\displaystyle \frac{\hat{u} _{r}}{\hat{L}_{r}}}. & & \displaystyle\end{eqnarray}$$

The reference values for the granular density and the pressure $\hat{p}$ are taken equal to $\hat{{\it\rho}}_{r}=\hat{{\it\rho}}$ , $\hat{p}_{r}=\hat{{\it\rho}}_{r}\hat{u} _{r}^{2}$ hence the non-dimensional density and pressure read,

(2.12a,b ) $$\begin{eqnarray}\displaystyle p_{i}={\displaystyle \frac{\hat{p}}{\hat{p}_{r}}},\quad {\it\rho}={\displaystyle \frac{\hat{{\it\rho}}}{\hat{{\it\rho}}_{r}}}=1. & & \displaystyle\end{eqnarray}$$

As regards $\hat{{\it\beta}}$ and $\hat{{\it\Gamma}}$ , based on the order-of-magnitude analysis of Varsakelis & Papalexandris (Reference Varsakelis and Papalexandris2011), we non-dimensionalise them with respect to $\hat{{\it\rho}}_{r}\hat{u} _{r}^{2}$ and $\hat{{\it\rho}}_{r}\hat{u} _{r}\hat{L}_{r}^{2}$ , respectively. Then,

(2.13a,b ) $$\begin{eqnarray}\displaystyle {\it\beta}=\frac{\hat{{\it\beta}}}{\hat{{\it\rho}}_{r}\hat{u} _{r}^{2}},\quad {\it\Gamma}=\frac{\hat{{\it\Gamma}}}{\hat{{\it\rho}}_{r}\hat{u} _{r}\hat{L}_{r}^{2}}. & & \displaystyle\end{eqnarray}$$

Subsequently, by introducing the expressions for $\hat{p}^{v}$ , $\hat{{\bf\tau}}$ $\hat{{\it\beta}}$ and ${\it\Gamma}$ to (2.1)–(2.3), and by performing the above non-dimensionalisation, we arrive at the non-dimensional form of the governing equations which reads

(2.14) $$\begin{eqnarray}\displaystyle \frac{\partial {\it\phi}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\phi}\boldsymbol{u})=0, & & \displaystyle\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle \frac{\partial {\it\phi}\boldsymbol{u}}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\phi}\boldsymbol{u}\boldsymbol{u})+\boldsymbol{{\rm\nabla}}(p{\it\phi}) & = & \displaystyle -\boldsymbol{{\rm\nabla}}\left(\frac{|\dot{{\bf\gamma}}|^{2}{d_{p}}^{2}}{I^{2}}{\it\phi}\right)+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\frac{{\it\mu}({\it\phi})|\dot{{\bf\gamma}}|{d_{p}}^{2}}{I^{2}}\dot{{\bf\gamma}}{\it\phi}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\frac{k_{2}{\it\phi}}{({\it\phi}_{max}-{\it\phi})^{2}}\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{{\rm\nabla}}{\it\phi}\right),\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle \frac{\partial {\it\phi}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}=\frac{1}{Re_{c}}\left(p-k_{1}{\it\phi}^{2}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\frac{k_{2}{\it\phi}}{({\it\phi}_{max}-{\it\phi})^{2}}\boldsymbol{{\rm\nabla}}{\it\phi}\right)\right). & & \displaystyle\end{eqnarray}$$

Equations (2.14)–(2.16) have four dimensionless groups, namely $d_{p}$ , $k_{1}$ , $k_{2}$ and $Re_{c}$ . The quantity $d_{p}$ describes the ratio between the particle diameter and the height of the domain. We remark that the validity of the continuum hypothesis requires that $d_{p}\sim O(10^{-k})$ , for some $k\in \mathbb{Z}$ , however, the exact value of $k$ is problem dependent, (Hutter & Rajagopal Reference Hutter and Rajagopal1994). Further, $k_{1}$ and $k_{2}$ measure the ratio between compaction to inertia effects and between intergranular surface variation to inertia effects. Finally, $Re_{c}$ is the ‘compaction Reynolds number’ which measures the ratio of viscous forces, due to compaction, to inertia.

Next, we advert to the choice of the base-flow profile. We invoke the time-independent, quasi-parallel flow assumption and search for profiles of the following form,

(2.17a−d ) $$\begin{eqnarray}\displaystyle u^{0}=u^{0}(y),\quad v^{0}=v^{0}(y),\quad p^{0}=p^{0}(y),\quad {\it\phi}^{0}={\it\phi}^{0}(y). & & \displaystyle\end{eqnarray}$$

Note that, although it is sufficient to impose that the base-flow is one-dimensional, the quasi-parallel flow assumption is more general. According to both experimental measurements and numerical predictions, expounded in Forterre & Pouliquen (Reference Forterre and Pouliquen2008), a base-flow profile that conforms with our postulates reads,

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\phi}^{0}(y)={\it\phi}^{0}=\text{constant}, & \displaystyle\end{eqnarray}$$

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle u^{0}(y)=y, & \displaystyle\end{eqnarray}$$

(2.20) $$\begin{eqnarray}\displaystyle & \displaystyle v^{0}(y)=0, & \displaystyle\end{eqnarray}$$

(2.21) $$\begin{eqnarray}\displaystyle & \displaystyle p^{0}(y)=k_{1}\left(\frac{{\it\phi}^{0}}{2-{\it\phi}^{0}}\right)^{2}\log \left(\frac{1}{1-{\it\phi}^{0}}\right). & \displaystyle\end{eqnarray}$$

It is straightforward to verify that the above relations constitute solutions to the non-dimensionalised (2.14)–(2.16).

Finally, the governing equations have to be endowed with suitable boundary conditions which, due to the quasi-parallel flow postulate, are required only for the upper and the lower planes. For the sake of simplicity, and keeping in mind the restrictions that arise from such a choice, a no-slip boundary condition is imposed on the granular velocity at both planes $u(0)=(0),v(0)=0$ and $u(1)=1,v(1)=0$ .

For the volume fraction, we first note that, in weakly non-local theories that involve gradients as internal variables such as the one that (2.1)–(2.3) have been based upon, the prescription of boundary conditions cannot be circumvented, unless the flow under consideration enjoys multiple symmetries, (Massoudi Reference Massoudi2007). For the problem at hand, the presence of the second-order differential operator in the right-hand side of the compaction equation in conjunction with the quasi-parallel flow postulate, imply that two boundary conditions are needed for the volume fraction as well. Following Wang & Hutter (Reference Wang and Hutter1999a ,Reference Wang and Hutter b ), a Dirichlet boundary condition is assigned to the volume fraction, i.e. ${\it\phi}(0)={\it\phi}(1)={\it\phi}^{0}$ . A plausible rationale of physical origin in favour of this choice is offered by Massoudi & Mehrabadi (Reference Massoudi and Mehrabadi2001) and Massoudi & Phuoc (Reference Massoudi and Phuoc2005).

As regards the pressure $p$ , a comment is in order. For general flowing conditions, e.g. transient, multi-dimensional flows, boundary conditions for the pressure are inherited from the momentum equation (2.14) as dictated by the Helmholtz decomposition and Ladyzhenskaya’s decomposition theorem, (Ladyzhenskaya & Solonnikov Reference Ladyzhenskaya and Solonnikov1978; Lions Reference Lions1996; Chorin & Marsden Reference Chorin and Marsden2000). For such cases, the pressure $p$ is computed as the solution to an appropriate Neumann problem by taking the divergence of the momentum equation and upon combination with the continuity equation, Varsakelis & Papalexandris (Reference Varsakelis and Papalexandris2014b ), Varsakelis, Monsorno & Papalexandris (Reference Varsakelis, Monsorno and Papalexandris2015). This is actually the exact analogue for the computation of the pressure field in the variable density Navier–Stokes equations. Importantly, as shown in Ladyzhenskaya (Reference Ladyzhenskaya1969) for the Navier–Stokes equations and in Varsakelis & Papalexandris (Reference Varsakelis and Papalexandris2010) for the equations at hand, imposing additional boundary conditions to the pressure field, results in an overdetermined problem. For the particular case at hand, however, the prescription of boundary conditions is trivial since $p$ is uniform throughout the domain; assuming that $p$ is continuous at the boundary, then a zero Neumann as well as the Dirichlet condition $p(0)$ or $p(1)=k_{1}({\it\phi}^{0}/(2-{\it\phi}^{0}))^{2}\log (1/(1-{\it\phi}^{0}))$ work equally well.

Additionally, we assume that the base-flow profiles $\boldsymbol{u}^{0}$ and ${\it\phi}^{0}$ satisfy the same boundary conditions as $\boldsymbol{u}$ and ${\it\phi}$ , respectively,

(2.22a,b ) $$\begin{eqnarray}\displaystyle u(0)=u^{0}(0)=0,\quad u(1)=u^{0}(1)=1, & & \displaystyle\end{eqnarray}$$

(2.23a,b ) $$\begin{eqnarray}\displaystyle v(0)=v^{0}(0)=0,\quad v(1)=v^{0}(1)=0, & & \displaystyle\end{eqnarray}$$

(2.24a,b ) $$\begin{eqnarray}\displaystyle {\it\phi}(0)={\it\phi}^{0}(0)={\it\phi}^{0},\quad {\it\phi}(1)={\it\phi}^{0}(1)={\it\phi}^{0}. & & \displaystyle\end{eqnarray}$$

3.1 Numerical method

Due to the complexity of the generalized eigenvalue problem (3.6), the analytical computation of the spectrum and the eigenfunctions is a formidable task. For this reason, we have resorted to the numerical approximation of these quantities. The voluminous literature on the numerical aspects of the Orr–Sommerfeld equation, arising in the study of the linear stability of the Navier–Stokes equations, corroborates that the accurate computation of the spectra and the eigenfunctions requires a judiciously chosen numerical method, (Orszag Reference Orszag1971; Huang & Sloan Reference Huang and Sloan1994; Dongarra, Straughanb & Walker Reference Dongarra, Straughanb and Walker1996). Thus, since the eigenvalue problem (3.6) is markedly more complex than the one induced by the Orr–Sommerfeld equation, its numerical approximation requires an elaborate numerical scheme as well.

In the present study, the generalized eigenvalue problem (3.6), along with the boundary conditions (3.7) and (3.8), is approximated via a Chebyshev collocation method. A detailed description of this algorithm is given in Varsakelis & Papalexandris (Reference Varsakelis and Papalexandris2015b ) and herein only its flowchart is presented.

1. (i) Primitive algebraic functions, i.e. the components of the base flow and the disturbances, are approximated via Chebyshev interpolation or equivalently by expansion in series of Chebyshev polynomials.

2. (ii) Algebraic differential operators are computed via applying a pseudo-spectral collocation method to the Chebyshev representation of the various quantities.

3. (iii) Boundary conditions are evaluated with the help of ghost cells to maintain the structure of the approximation throughout the computational domain.

4. (iv) The resulting discretized eigenvalue problem is computed via the QR-decomposition. The computation of the eigenvalues and the corresponding eigenvectors is performed in a series of adapted (denser) grids so that spurious modes are excluded and higher accuracy is achieved.

5. (v) All Chebyshev approximations require a very strict convergence criterion, close to machine precision.

All computations have been carried with the Chebfun computational suite, (Hale & Trefethen Reference Hale and Trefethen2014). The robustness and accuracy of Chebfun have been assessed in a multitude of cases, including the computation of the spectrum of the Orr–Sommerfeld operator.

4.1 Parametric study with respect to ${\hat{H}}$ and ${\it\phi}^{0}$

Figure 2.

Selected growth rates ${\it\sigma}$ plotted against the Couette gap ratio ${\hat{H}}/{\hat{H}}_{ref}$ , for various initial concentrations ${\it\phi}^{0}$ and for ${\it\alpha}=1$ and $\hat{U} _{w}=\hat{U} _{w,ref}=1$ . The flow remains linearly stable for $1\leqslant {\hat{H}}/{\hat{H}}_{ref}\leqslant 30$ and $0.4\leqslant {\it\phi}^{0}\leqslant 0.55$ . Further, according to our predictions, shear flows at higher Couette gaps and with lower concentrations correspond to higher growth rates.

For the scopes of this parametric study we modify ${\hat{H}}$ and ${\it\phi}^{0}$ while the remaining variables have maintained their values with respect to the reference configuration. Figure 2 shows representative growth rates ${\it\sigma}$ versus the Couette gap ratio ${\hat{H}}/{\hat{H}}_{ref}$ for different initial concentrations and for ${\it\alpha}=1$ , $\hat{U} _{w}=\hat{U} _{w,ref}=1$ . For all values of these parameters considered herein, namely $1\leqslant {\hat{H}}/{\hat{H}}_{ref}\leqslant 30$ , the flow is predicted to be linearly stable.

As asserted by figure 2, the Couette gap plays a destabilizing role in the flows of interest. Indeed, it correlates positively with the growth rate, however, its rate of change experiences a significant decrease. Therefore, shear flows confined in very wide channels, and at the limit unbounded ones, might be susceptible to instabilities. On the other hand, the magnitude of the initial concentration is found to have a stabilizing effect. This is not surprising since, as the concentration increases, frictional forces become dominant and render relative motion between the grains more difficult. As a consequence, disturbances in the form of travelling waves that conform with the governing equations at hand correspond to increasingly smaller growth rates and, therefore, decay faster.

The combined effect of modifying both the Couette gap and the initial concentration does not yield any additional information about the interplay of these two parameters. Indeed, the simultaneous modification of both parameters amounts to simply shifting the curves depicted in figure 2 accordingly. Nonetheless, although the magnitude of the growth rate is predicted to be a complicated function of ${\hat{H}}$ and ${\it\phi}^{0}$ , its sign remains invariant.

Previous stability analyses that have reported the ( ${\hat{H}},{\it\phi}^{0}$ ) diagram have also affirmed the destabilizing and stabilizing role of the Couette gap and the initial concentration, respectively, (Wang et al. Reference Wang, Jackson and Sundaresan1996; Alam & Nott Reference Alam and Nott1998; Nott et al. Reference Nott, Alam, Agrawal, Jackson and Sundaresan1999). Most notably, however, these studies have reported the existence of unstable modes. We now examine whether this disparity is expected and justified.

As mentioned in the introduction, extant studies in the linear stability of the flows of interest have focused on either dilute or rapid granular flows. These flows belong to the collisional regime, where collisions between grains dominate over frictional effects. In turn, this is reflected in the functional form of the governing equations and in particular in the constitutive expressions for the normal and the shear stresses that act on the granular conglomerate. By contrast, the focus of the present study pertains precisely to the regime where frictional effects are important. In the context of the employed theory, these effects are accounted for by embodying appropriate viscous terms to the granular material’s stress tensor. In other words, even if we ignore compaction effects, the differences between the stress tensor employed in kinetic theoretical models and herein are large enough to lead to different stability diagrams.

Being dissimilar, however, does not imply that the two stability diagrams are in contradiction to each other. Rather, it reflects the fact that the dynamics of granular materials in the aforementioned two regimes bear strong differences as well. Consequently, reconciling the findings of the present study with those obtained from kinetic theoretical models requires a ‘hybrid’ theory that is predicated on continuum mechanics but additionally incorporates aspects of kinetic theory of granular gases in a consistent manner. Nevertheless, to the best of our knowledge, a widely accepted theory with the above characteristics has yet to appear in the literature.

4.2 Parametric study with respect to $\hat{U} _{w}$ and ${\it\phi}^{0}$

For this parametric study, we have kept the value of the Couette gap fixed and have modified the wall velocity and the initial concentration. More specifically, we have considered values $1\leqslant \hat{U} /\hat{U} _{w,ref}\leqslant 100$ and $0.4\leqslant {\it\phi}^{0}\leqslant 0.55$ . The results of this parametric study are presented in figure 3, that portrays representative growth rates ${\it\sigma}$ versus the wall velocity ratio $\hat{U} /\hat{U} _{w,ref}$ for various initial concentrations. Two conclusions can be directly drawn. First, the flow remains linearly stable for all values of $\hat{U} /\hat{U} _{w,ref}$ and ${\it\phi}^{0}$ considered herein, and second, the effects of increasing the initial concentration are stabilizing, in line with what was observed in the previous parametric study; clearly, the latter conclusion is nothing but a sanity test for the self-consistency of the acquired results.

Figure 3.

Representative growth rates ${\it\sigma}$ plotted against the wall velocity ratio $\hat{U} _{w}/\hat{U} _{w,ref}$ , for various initial concentrations and for ${\it\alpha}=1$ , ${\hat{H}}/{\hat{H}}_{ref}=1$ . The flow remains linearly stable for $1\leqslant \hat{U} _{w}/\hat{U} _{w,ref}\leqslant 100$ and $0.4\leqslant {\it\phi}^{0}\leqslant 0.55$ . Increasing $\hat{U} _{w}/\hat{U} _{w,ref}$ initially translates into lower growth rates which is subsequently followed by an abrupt increase and the presence of a local maximum. The post-maximum behaviour of ${\it\sigma}$ switches from monotonically increasing to decreasing as ${\it\phi}^{0}$ increases.

By contrast, the role of $\hat{U} _{w}$ and its interplay with ${\it\phi}^{0}$ are substantially more intricate. According to figure 3, the curves of the growth rates comprise three segments. In the first segment, increasing $\hat{U} /\hat{U} _{w,ref}$ translates into gradually decaying growth rates. This segment spans from the start of the horizontal axis until the point where a minimum value has been attained in the cavity. Once this minimum value is crossed, the sign of the correlation changes and this marks the onset of the second segment. Here, ${\it\sigma}$ experiences a rapid increase and reaches a local maximum that designates the border of the second segment. The third segment extends from this local maximum until the end of the horizontal axis. The shape of ${\it\sigma}$ in this segment depends explicitly on ${\it\phi}^{0}$ and signifies the competition between ${\it\phi}^{0}$ and $\hat{U} _{w}$ . For lower values of ${\it\phi}^{0}$ , ${\it\sigma}$ correlates positively with $\hat{U} /\hat{U} _{w,ref}$ suggesting that the increasing magnitude of the applied shear, induced by the wall velocity, eventually acts as a destabilizing factor. As ${\it\phi}^{0}$ increases, this correlation appears to progressively change sign. However, it should be stressed that, according to our numerical experiments, the negative correlation does not hold for arbitrary high wall velocities. As a matter of fact, our results conclusively show that for each concentration ${\it\phi}^{0}$ there exists a threshold wall velocity upon which ${\it\sigma}$ becomes a monotonically increasing function of $\hat{U} /\hat{U} _{w,ref}$ . In other words, increasing the concentration amounts to translating the borderline from where the destabilizing effects of the applied shear gradually counterbalance the frictional forces, accordingly.

One of the anonymous referees of this paper remarked that the complex, non-monotonous dependence of the growth rate ${\it\sigma}$ on $\hat{U} /\hat{U} _{w,ref}$ , depicted in figure 3, can be partially attributed to the interplay between the shear-thickening behaviour that the ${\it\mu}(I)$ -rheology predicts and the destabilizing role of the applied shear. Indeed, as $\hat{U} /\hat{U} _{w,ref}$ increases so do (i) the magnitude of the frictional forces and hence of the (stabilizing) dissipation, and (ii) the strength of the applied shear. With these in mind, we observe that, in the first of the aforementioned three segments that the growth rate curves consist of, increasing $\hat{U} /\hat{U} _{w,ref}$ renders the flow more stable; hence in this regime the stabilizing effect of dissipation dominates. Nevertheless, as $\hat{U} /\hat{U} _{w,ref}$ increases further, the (still increasing) dissipation is progressively offset by the destabilizing effect of the (also increasing) applied shear. Although this offset is oscillatory, hinting that secondary concentration-dependent mechanisms also influence the flow, our numerical predictions suggest that the shear-thickening behaviour induced by the ${\it\mu}(I)$ -rheology can only transiently stabilize the flows of interest.

4.3 Parametric study with respect to ${\hat{H}}$ and $\hat{U} _{w}$

This parametric study concerns the effects of the combined modification of ${\hat{H}}$ and $\hat{U} _{w}$ on the stability properties of flow. Figure 4(a) depicts the predicted stability diagram for $1\leqslant {\hat{H}}/{\hat{H}}_{ref}\leqslant 30$ and $1\leqslant \hat{U} /\hat{U} _{w,ref}\leqslant 100$ while the initial concentration is ${\it\phi}^{0}=0.55$ . It can be readily inferred that, for values of $\hat{U} _{w}/\hat{U} _{w,ref}$ and ${\hat{H}}/{\hat{H}}_{ref}$ inside the domain defined by the neutral stability curve, the flow becomes linearly unstable. Accordingly, the maximum growth rate is ${\it\sigma}_{max}=0.0054179$ and corresponds to a wall velocity $\hat{U} _{w}=0.67~\text{m}~\text{s}^{-1}$ and Couette gap ${\hat{H}}=3~\text{m}$ . We remark that ${\it\sigma}_{max}$ corresponds to the eigenvalue with the largest real part and not to the first eigenvalue with positive real part. In this respect, ${\it\sigma}_{max}$ is associated to a cutoff frequency, similar to the one obtained by Forterre (Reference Forterre2006), i.e. the spectrum is not unbounded as is the case of Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015). Another interesting finding concerns the magnitude of ${\it\sigma}_{max}$ as a function of the initial concentration ${\it\phi}^{0}$ . According to our results, in the unstable regime, ${\it\sigma}_{max}$ is negatively correlated to ${\it\phi}^{0}$ suggesting a reversal in the behaviour depicted in figures 2 and 3.

The prediction of an instability regime for values of ${\hat{H}}/{\hat{H}}_{ref}$ , $\hat{U} _{w}/\hat{U} _{w,ref}$ and ${\it\phi}^{0}$ that have been employed in the previous parametric studies is not at odds with the results of these studies. Rather, it reflects the aforementioned fact that, due to the particular dependence of the dimensionless groups $d_{p}$ , $k_{1}$ , $k_{2}$ and $Re_{c}$ on $\hat{L}_{r}={\hat{H}}$ and $\boldsymbol{U}_{r}=\hat{U} _{w}$ , changes of ${\hat{H}}$ cannot be represented by changes in $\hat{U} _{w}$ . Thus, as confirmed by our numerical predictions, a combined modification of these parameters is necessary for the appearance of an unstable mode. This is further confirmed by figure 4(b,c) that illustrate the same stability diagram but for different concentrations. From these figures, we can additionally infer that modifying the concentration shifts the instability regime along the main diagonal of the axes, accordingly.

The stability diagrams portrayed in figure 4(a–c) assert that for small values of the Couette gap and of the wall velocity, the flow is linearly stable. As far as the Couette gap is concerned, this finding is routinely recovered in stability analyses based on kinetic theoretical models. However, as regards the wall velocity, this appears to be a novel, consistent result. The scarcity of results in this direction may be attributed to the dimensionless groups that kinetic theoretical models admit. In general, they are markedly different from the ones discussed in the previous section and do not incorporate the reference velocity.

Subsequently, we shift our focus to the manifestation of the detected instability and its effects on the granular concentration and velocity fields. From the various unstable modes that we acquired, we opt to study the case corresponding to the maximum growth rate ${\it\sigma}_{max}=0.0048$ ( $\hat{U} _{w}/\hat{U} _{w,ref}=67$ , ${\hat{H}}/{\hat{H}}_{ref}=30$ , ${\it\phi}^{0}=0.55$ ). Figure 5 shows contour plots of the granular concentration ${\it\phi}(x,y)$ ,

(4.1) $$\begin{eqnarray}\displaystyle {\it\phi}(x,y)={\it\phi}^{0}+{\it\phi}^{1}(x,y,t)=0.55+\text{Re}(\overline{{\it\phi}}(y)\,\text{e}^{\text{i}(x-{\it\sigma}_{max}t)}), & & \displaystyle\end{eqnarray}$$

at $t=0.1$ . Superimposed in this figure is a vector plot of the corresponding (normalized) velocity field so that the motion of granular material can be pictured.

Figure 4.

Contour plot of growth rates ${\it\sigma}$ for $1\leqslant \hat{U} _{w}/\hat{U} _{w,ref}\leqslant 100$ and $1\leqslant {\hat{H}}/{\hat{H}}_{ref}\leqslant 30$ . (a) ${\it\phi}^{0}=0.55$ , (b) ${\it\phi}^{0}=0.53$ and (c) ${\it\phi}^{0}=0.5$ . The wavenumber is ${\it\alpha}=1$ . The thick line (–) designates the neutral stability curve where ${\it\sigma}=0$ and the filled diamond ( $\blacklozenge$ ) the maximum growth rate for each case. For values of $\hat{U} _{w}/\hat{U} _{w,ref}$ and ${\hat{H}}/{\hat{H}}_{ref}$ inside the region defined by the neutral stability curve, the flow is predicted to be linearly unstable. For case (a) the maximum growth rate is predicted to be ${\it\sigma}=0.0054179$ and corresponds to the point $\hat{U} _{w}/\hat{U} _{w,ref}=67$ , ${\hat{H}}/{\hat{H}}_{ref}=30$ . As the concentration decreases, the instability regime is shifted downwards along the diagonal of the axes asserting that offsetting the stabilizing effect of the frictional forces requires a lower wall velocity and Couette gap.

Figure 5.

Contour plots of the granular concentration at $t=0.1$ for $\hat{U} _{w}/\hat{U} _{w,ref}=67$ , ${\hat{H}}/{\hat{H}}_{ref}=30$ and ${\it\sigma}=0.0048$ . The concentration is ${\it\phi}^{0}=0.55$ . Superimposed is a vector plot of the velocity field. The predicted instability manifests itself through the formation of two granular bulbs of high and low concentration that move away from and towards the upper wall, respectively. In turn, this results in particle migration.

Figure 6.

Contour plots of the granular vorticity at $t=0.1$ for $\hat{U} _{w}/\hat{U} _{w,ref}=67$ , ${\hat{H}}/{\hat{H}}_{ref}=30$ and ${\it\sigma}=0.0048$ . The concentration is ${\it\phi}^{0}=0.55$ . As a consequence of the predicted instability, two counter-rotating vortices have appeared that are stretched in the streamwise direction by a factor of $3.5$ . The relative position with respect to the granular bulbs of figure 5 indicates a form of preferential concentration.

The initial uniform concentration ${\it\phi}^{0}=0.55$ has been substantially modified and two granular bulbs of high and low concentration have emerged, with the former preceding the latter in the streamwise direction. These bulbs arise due to the interaction between the disturbance ${\it\phi}^{1}(x,y,t)$ and the base-flow concentration ${\it\phi}^{0}$ , i.e. they constitute the manifestation of the predicted instability. We remark that the second bulb constitutes a mirror image of the first one, stemming directly from the fact that we have employed a normal mode analysis. The same observation applies to the velocity and vorticity fields as well. Thus, it suffices to adhere to the examination of the first (high concentration) bulb and the associated velocity and vorticity structures.

The motion of the first bulb takes place along curves of clockwise orientation and shows a vertical displacement of particles from the upper plane towards the middle of the domain. In other words, the instability gives rise to a shear dilatancy that, in turn, induces particle migration.

At this point it is useful to note that the prediction of bulbs in granular shear flows is not without precedence. Indeed, Conway & Glasser (Reference Conway and Glasser2004) and Conway, Liu & Glasser (Reference Conway, Liu and Glasser2006), who studied granular shear flows via molecular dynamics simulations, documented the formation and propagation of elongated and oscillatory structures of high particle concentration located in the middle of the domain. These structures are qualitatively similar to those predicted by our linear stability analysis; cf. Figure 5, and provide further evidence on the appropriateness of continuum mechanical models to describe the primary properties of the flows of interest.

Further insight to the properties of the predicted instability can be obtained upon examination of the vorticity field of the fluid. Figure 6 depicts contours of the granular vorticity at $t=0.1$ . We observe that, following the onset of the instability, two counter-rotating vortices are formed that are stretched in the streamwise direction by a factor of, approximately, 3.5. More important, however, is the positioning of these vortices relative to the bulbs. In particular, the maximum concentration of the first bulb coincides with the area of minimum vorticity and vice-versa. Therefore, the manifestation of the instability through particle migration is also accompanied by a form of preferential concentration.

4.4 The role of the intergranular stress and the equilibrated stress coefficient

As mentioned above, the exact form of the compaction equation for granular flows remains a subject of debate. Moreover, the experimental validation or refutation of existing forms is a cumbersome task due to the multitude of assumptions that the derivation of these equations is based upon and the presence of terms for which experimental correlations are difficult to be acquired. It is, therefore, appropriate to comment on the predictions of the stability analysis concerning the effect of the intergranular stress ${\it\beta}$ and the equilibrated stress coefficient ${\it\Gamma}$ on the flow properties.

Summing up, the results of the stability analysis assert that large Couette gaps, high wall velocities and low concentrations are destabilizing mechanisms whilst the converse is also true. On the other hand, by virtue of (2.9), (2.10) and (2.13), we can infer that as ${\hat{H}}$ , $\hat{U} _{w}$ increase and ${\it\phi}$ decreases, then the magnitude of ${\it\beta}$ and ${\it\Gamma}$ decrease. In other words, as the configuration of the flow approaches the region of instability, the strength of both the intergranular stresses and of the tendency of the granular material to suppress inhomogeneities are accordingly reduced. By contrast, the exact opposite scenario is observed following the decrease of ${\hat{H}}$ , $\hat{U} _{w}$ and the increase of ${\it\phi}$ . Collectively, these observations suggest that the intergranular stress ${\it\beta}$ and the equilibrated stress coefficient ${\it\Gamma}$ play a stabilizing role in the flows of interest.