Shear-banding instability in arbitrarily inelastic granular shear flows

TL;DR

This paper analyzes the shear-banding instability in inelastic granular flows using linear stability analysis of hydrodynamic equations, highlighting the importance of constitutive relations and identifying criteria for instability onset.

Contribution

It provides a detailed linear stability analysis of shear-banding in inelastic granular flows, emphasizing the role of constitutive relations and deriving criteria for instability onset.

Findings

Constitutive relations significantly influence shear-banding predictions.

Two global criteria for the onset of shear-banding instability are identified.

The study explores effects of restitution coefficient, channel width, and density.

Abstract

One prototypical instability in granular flows is the shear-banding instability, in which a uniform granular shear flow breaks into alternating bands of dense and dilute clusters of particles having low and high shear (shear stress or shear rate), respectively. In this work, the shear-banding instability in an arbitrarily inelastic granular shear flow is analyzed through the linear stability analysis of granular hydrodynamic equations closed with Navier--Stokes-level constitutive relations. It is shown that the choice of appropriate constitutive relations plays an important role in predicting the shear-banding instability. A parametric study is carried out to study the effect of the restitution coefficient, channel width and mean density. Two global criteria relating the control parameters are found for the onset of the shear-banding instability.

Equations

$$\begin{array}{rl}\left(\frac{\partial }{\partial \bar{t}}+\bar{\mathbf{U}}\cdot \bar{\nabla }\right)\bar{\rho }+\bar{\rho }\bar{\nabla }\cdot \bar{\mathbf{U}}& =0,\\ \bar{\rho }\left(\frac{\partial }{\partial \bar{t}}+\bar{\mathbf{U}}\cdot \bar{\nabla }\right)\bar{\mathbf{U}}+\bar{\nabla }\cdot \bar{\mathbf{P}}& =0,\\ \bar{\rho }\left(\frac{\partial }{\partial \bar{t}}+\bar{\mathbf{U}}\cdot \bar{\nabla }\right)\bar{T}+\bar{\mathbf{P}}:\bar{\nabla }\bar{\mathbf{U}}+\bar{\nabla }\cdot \bar{\mathbf{q}}& =-\bar{\mathcal{D}}.\end{array}\}$$

$$\bar{\mathbf{P}}$$

$$\bar{\mathbf{q}}$$

$$\bar{\mathcal{D}}$$

$$\begin{array}{rlrl}\bar{p}(\varphi ,\bar{T},e)& =f_1(\varphi ,e){\rho }_p\bar{T},& & \\ \bar{\eta }(\varphi ,\bar{T},e)& =f_2(\varphi ,e){\rho }_p{d}_p\sqrt{\bar{T}},& \bar{\gamma }(\varphi ,\bar{T},e)& =f_3(\varphi ,e){\rho }_p{d}_p\sqrt{\bar{T}},\\ \bar{\kappa }(\varphi ,\bar{T},e)& =f_4(\varphi ,e){\rho }_p{d}_p\sqrt{\bar{T}},& \bar{\mu }(\varphi ,\bar{T},e)& =f_{4h}(\varphi ,e){\rho }_p{d}_p\bar{T}\sqrt{\bar{T}},\\ {\bar{\mathcal{D}}}_0(\varphi ,\bar{T},e)& =f_5(\varphi ,e)\frac{{\rho }_p}{d_p}\bar{T}\sqrt{\bar{T}},& {\bar{\mathcal{D}}}_1(\varphi ,\bar{T},e)& =f_{5u}(\varphi ,e){\rho }_p\bar{T}.\end{array}\}$$

$$\displaystyle\left.\begin{aligned} f_{1}(\phi,e)&=\phi\,\big{[}1+(1+e)\,G(\phi)\big{]},\\ f_{3}(\phi,e)&=\frac{1}{\sqrt{\pi}}(1+e)\left(1-\frac{c}{32}\right)\phi\,G(\phi),\\ f_{2}(\phi,e)&=\frac{\sqrt{\pi}}{8}\frac{\left[1-\frac{1}{4}(1+e)(1-3e)\,G(\phi)\right]\left[1+\frac{1}{2}(1+e)\,G(\phi)\right]}{\nu_{\eta}^{\ast}-\frac{1}{2}\zeta_{0}^{\ast}}+\frac{1}{2}f_{3}(\phi),\\ f_{4}(\phi,e)&=\frac{\sqrt{\pi}}{2}\left[1+\frac{3}{4}(1+e)\,G(\phi)\right]\kappa_{k}^{\ast}+\frac{1}{\sqrt{\pi}}(1+e)\left(1+\frac{7c}{32}\right)\phi\,G(\phi),\\ f_{4h}(\phi,e)&=\frac{\sqrt{\pi}}{2\phi}\left[1+\frac{3}{4}(1+e)\,G(\phi)\right]\mu_{k}^{\ast},\\ f_{5}(\phi,e)&=\frac{4}{\sqrt{\pi}}(1-e^{2})\left(1+\frac{3c}{32}\right)\phi\,G(\phi),\\ f_{5u}(\phi,e)&=\frac{3}{2}(1-e^{2})\left[\frac{3}{32}\frac{\frac{1}{8}\omega^{\ast}-\frac{c}{3}(1+e)(1-3e)}{\nu_{\zeta}^{\ast}-\frac{3}{4}(1-e^{2})}-1\right]\phi\,G(\phi),\end{aligned}\right\}$$

$$c$$

$$\chi (\varphi )=\{\begin{array}{ll}\frac{1-7\varphi /16}{(1-\varphi {)}^2}& \text{mbox}for0\le \varphi <{\varphi }_f,\\ \frac{1-7{\varphi }_f/16}{(1-{\varphi }_f{)}^2}\left(\frac{{\varphi }_c-{\varphi }_f}{{\varphi }_c-\varphi }\right)& \text{mbox}for{\varphi }_f\le \varphi \le {\varphi }_c.\end{array}$$

$$\displaystyle\left.\begin{aligned} \zeta_{0}^{\ast}&=\frac{1}{2}(1-e^{2})\left(1+\frac{3c}{32}\right)\,\chi(\phi),\\ \nu_{\eta}^{\ast}&=\frac{1}{8}(1+e)(7-3e)\left(1+\frac{7c}{32}\right)\,\chi(\phi),\\ \kappa_{k}^{\ast}&=\frac{1+c+\frac{3}{8}(1+e)^{2}\left[2e-1+\frac{c}{2}(1+e)\right]G(\phi)}{2\big{(}\nu_{\kappa}^{\ast}-2\zeta_{0}^{\ast}\big{)}},\\ \mu_{k}^{\ast}&=\frac{\zeta_{0}^{\ast}\kappa_{k}^{\ast}(1+\phi\,\partial_{\phi}{\ln{\chi}})+\frac{c}{4}+\frac{3}{8}(1+e)\left(1+\frac{1}{2}\phi\,\partial_{\phi}{\ln{\chi}}\right)\left[e(e-1)+\frac{c}{12}(14-3e+3e^{2})\right]G(\phi)}{\nu_{\kappa}^{\ast}-\frac{3}{2}\zeta_{0}^{\ast}},\\ \nu_{\kappa}^{\ast}&=\frac{1}{4}(1+e)\left[1+\frac{15}{4}(1-e)+\frac{365-273e}{128}c\right]\,\chi(\phi),\\ \omega^{\ast}&=(1+e)\left[(1-e^{2})(5e-1)+\frac{c}{12}(41-69e+3e^{2}-15e^{3})\right],\\ \nu_{\zeta}^{\ast}&=\frac{1}{192}(1+e)(185-153e+30e^{2}-30e^{3}).\end{aligned}\right\}$$

$$\begin{array}{rl}{f}_1^k(\varphi ,e)& =\varphi ,\\ f_3^k(\varphi ,e)& =0,\\ f_2^k(\varphi ,e)& =\frac{\sqrt{\pi }}{8}\frac{1-\frac{1}{4}(1+e)(1-3e)G(\varphi )}{\left({\nu }_{\eta }^{\ast }-\frac{1}{2}{\zeta }_0^{\ast }\right)},\\ f_4^k(\varphi ,e)& =\frac{\sqrt{\pi }}{2}{\kappa }_k^{\ast }\\ f_{4h}^k(\varphi ,e)& =\frac{\sqrt{\pi }}{2\varphi }{\mu }_k^{\ast },\end{array}\}$$

$$\begin{array}{rl}{f}_1^c(\varphi ,e)& =(1+e)\varphi G(\varphi ),\\ f_3^c(\varphi ,e)& =\frac{1}{\sqrt{\pi }}(1+e)\left(1-\frac{c}{32}\right)\varphi G(\varphi ),\\ f_2^c(\varphi ,e)& =\frac{1}{2}(1+e)G(\varphi )f_2^k+\frac{1}{2}{f}_3(\varphi ),\\ f_4^c(\varphi ,e)& =\frac{3\sqrt{\pi }}{8}(1+e)G(\varphi ){\kappa }_k^{\ast }+\frac{1}{\sqrt{\pi }}(1+e)\left(1+\frac{7c}{32}\right)\varphi G(\varphi ),\\ f_{4h}^c(\varphi ,e)& =\frac{3\sqrt{\pi }}{8}(1+e)\chi (\varphi ){\mu }_k^{\ast }.\end{array}\}$$

$$f_5(\varphi \to 0,e)=\frac{4}{\sqrt{\pi }}(1-e^2)\left(1+\frac{3c}{32}\right){\varphi }^2\text{and}{f}_{5u}(\varphi \to 0,e)=0.$$

$$\bar{u}\left(\pm \frac{h}{2}\right)=\pm \frac{U_w}{2}\text{and}{\frac{\partial \bar{T}}{\partial \bar{y}}|}_{\bar{y}=\pm h/2}=0.$$

$$\nabla =h\bar{\nabla },\mathbf{U}=\frac{1}{U_w}\bar{\mathbf{U}},t=\frac{U_w}{h}\bar{t},\rho =\frac{\bar{\rho }}{{\rho }_p},T=\frac{H^2}{U_w^2}\bar{T},$$

$$\displaystyle\left.\begin{aligned} \left(\frac{\partial}{\partial t}+\bm{U}\cdot\bm{\nabla}\right)\phi&=-\phi\,\bm{\nabla}\cdot\bm{U},\\ \phi\left(\frac{\partial}{\partial t}+\bm{U}\cdot\bm{\nabla}\right)\bm{U}&=\frac{1}{H^{2}}\bigg{[}-\bm{\nabla}p+\bm{\nabla}\left(\lambda\,\bm{\nabla}\cdot\bm{U}\right)+\bm{\nabla}\cdot\Big{\{}\eta\,\Big{(}\bm{\nabla}\bm{U}+(\bm{\nabla}\bm{U})^{\mathsf{T}}\Big{)}\Big{\}}\bigg{]},\\ \frac{d}{2}\,\phi\left(\frac{\partial}{\partial t}+\bm{U}\cdot\bm{\nabla}\right)T&=\frac{1}{H^{2}}\bm{\nabla}\cdot\left(\kappa\,\bm{\nabla}T+\mu\,\bm{\nabla}\phi\right)-p\,\bm{\nabla}\cdot\bm{U}+\lambda\,(\bm{\nabla}\cdot\bm{U})^{2}\\ &\quad+\eta\,\Big{\{}\bm{\nabla}\bm{U}+(\bm{\nabla}\bm{U})^{\mathsf{T}}\Big{\}}:\bm{\nabla}\bm{U}-\mathcal{D}_{0}-\mathcal{D}_{1}\,\bm{\nabla}\cdot\bm{U},\end{aligned}\right\}$$

$$\begin{array}{rlrl}p(\varphi ,T,e)& =f_1(\varphi ,e)T,& \eta (\varphi ,T,e)& =f_2(\varphi ,e)\sqrt{T},\\ \gamma (\varphi ,T,e)& =f_3(\varphi ,e)\sqrt{T},& \lambda (\varphi ,T,e)& =\gamma -\frac{2}{d}\eta ,\\ \kappa (\varphi ,T,e)& =f_4(\varphi ,e)\sqrt{T},& \mu (\varphi ,T,e)& =f_{4h}(\varphi ,e)T\sqrt{T},\\ {\mathcal{D}}_0(\varphi ,T,e)& =f_5(\varphi ,e)T\sqrt{T},& {\mathcal{D}}_1(\varphi ,T,e)& =f_{5u}(\varphi ,e)T.\end{array}\}$$

$$u\left(\pm \frac{1}{2}\right)=\pm \frac{1}{2}\text{mbox}and{\frac{\partial T}{\partial y}|}_{y=\pm 1/2}=0.$$

$$\begin{array}{rl}\frac{\partial \varphi }{\partial t}+v\frac{\partial \varphi }{\partial y}+\varphi \frac{\partial v}{\partial y}& =0,\\ \varphi \left(\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial y}\right)& =\frac{1}{H^2}\frac{\partial }{\partial y}\left(\eta \frac{\partial u}{\partial y}\right),\\ \varphi \left(\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial y}\right)& =\frac{1}{H^2}\left[-\frac{\partial p}{\partial y}+\frac{\partial }{\partial y}\left(\lambda \frac{\partial v}{\partial y}\right)+2\frac{\partial }{\partial y}\left(\eta \frac{\partial v}{\partial y}\right)\right],\\ \varphi \left(\frac{\partial T}{\partial t}+v\frac{\partial T}{\partial y}\right)& =\frac{1}{H^2}\frac{\partial }{\partial y}\left(\kappa \frac{\partial T}{\partial y}+\mu \frac{\partial \varphi }{\partial y}\right)-p\frac{\partial v}{\partial y}+\lambda {\left(\frac{\partial v}{\partial y}\right)}^2+\eta {\left(\frac{\partial u}{\partial y}\right)}^2+2\eta {\left(\frac{\partial v}{\partial y}\right)}^2-{\mathcal{D}}_0-{\mathcal{D}}_1\frac{\partial v}{\partial y}.\end{array}\}$$

$$\varphi =\varphi (y),u=u(y),v=0,T=T(y).$$

$$\frac{\mathrm{d}}{\mathrm{d}y}\left(\eta \frac{\mathrm{d}u}{\mathrm{d}y}\right)=0,-\frac{\mathrm{d}p}{\mathrm{d}y}=0,\frac{1}{H^2}\frac{\mathrm{d}}{\mathrm{d}y}\left(\kappa \frac{\mathrm{d}T}{\mathrm{d}y}+\mu \frac{\mathrm{d}\varphi }{\mathrm{d}y}\right)+\eta {\left(\frac{\mathrm{d}u}{\mathrm{d}y}\right)}^2-{\mathcal{D}}_0=0.$$

$$\begin{array}{c}{f}_2(\varphi ,e)\sqrt{T}\frac{\partial u}{\partial y}=c_1,f_1(\varphi ,e)T=c_2,\\ \frac{1}{H^2}\frac{\mathrm{d}}{\mathrm{d}y}\left(f_4(\varphi ,e)\sqrt{T}\frac{\mathrm{d}T}{\mathrm{d}y}+f_{4h}(\varphi ,e)T\sqrt{T}\frac{\mathrm{d}\varphi }{\mathrm{d}y}\right)+f_2(\varphi ,e)\sqrt{T}{\left(\frac{\mathrm{d}u}{\mathrm{d}y}\right)}^2-f_5(\varphi ,e)T\sqrt{T}=0,\end{array}\}$$

$${\varphi }^0(y)={\varphi }_0,u^0(y)=y,v^0(y)=0,T^0(y)=f_2^0/f_5^0=T_0,$$

$$\begin{array}{r}\varphi (t,y)={\varphi }^0(y)+{\varphi }^{\prime }(t,y),u(t,y)=u^0(y)+u^{\prime }(t,y),v(t,y)=v^0(y)+v^{\prime }(t,y),T(t,y)=T^0(y)+T^{\prime }(t,y).\end{array}$$

$$\begin{array}{rl}\frac{\partial {\varphi }^{\prime }}{\partial t}& =-{\varphi }_0\frac{\partial v^{\prime }}{\partial y},\\ \frac{\partial u^{\prime }}{\partial t}& =\frac{{\eta }_{\varphi }^0}{{\varphi }_0{H}^2}\frac{\partial {\varphi }^{\prime }}{\partial y}+\frac{{\eta }^0}{{\varphi }_0{H}^2}\frac{{\partial }^2{u}^{\prime }}{\partial y^2}-v^{\prime }+\frac{{\eta }_T^0}{{\varphi }_0{H}^2}\frac{\partial T^{\prime }}{\partial y},\\ \frac{\partial v^{\prime }}{\partial t}& =-\frac{p_{\varphi }^0}{{\varphi }_0{H}^2}\frac{\partial {\varphi }^{\prime }}{\partial y}+\frac{{\lambda }^0+2{\eta }^0}{{\varphi }_0{H}^2}\frac{{\partial }^2{v}^{\prime }}{\partial y^2}-\frac{p_T^0}{{\varphi }_0{H}^2}\frac{\partial T^{\prime }}{\partial y},\\ \frac{\partial T^{\prime }}{\partial t}& =\frac{1}{{\varphi }_0}\left(\frac{{\mu }^0}{H^2}\frac{{\partial }^2}{\partial y^2}+{\eta }_{\varphi }^0-{\mathcal{D}}_{0,\varphi }^0\right){\varphi }^{\prime }+\frac{2{\eta }^0}{{\varphi }_0}\frac{\partial u^{\prime }}{\partial y}-\frac{p^0+{\mathcal{D}}_1^0}{{\varphi }_0}\frac{\partial v^{\prime }}{\partial y}+\frac{1}{{\varphi }_0}\left(\frac{{\kappa }^0}{H^2}\frac{{\partial }^2}{\partial y^2}+{\eta }_T^0-{\mathcal{D}}_{0,T}^0\right)T^{\prime },\end{array}\}$$

$$\frac{\partial \mathbf{X}}{\partial t}=\mathcal{L}\mathbf{X},$$

$$\mathcal{L}=\left[\begin{array}{cccc}0& 0& -{\varphi }_0\frac{\partial }{\partial y}& 0\\ \frac{{\eta }_{\varphi }^0}{{\varphi }_0{H}^2}\frac{\partial }{\partial y}& \frac{{\eta }^0}{{\varphi }_0{H}^2}\frac{{\partial }^2}{\partial y^2}& -1& \frac{{\eta }_T^0}{{\varphi }_0{H}^2}\frac{\partial }{\partial y}\\ -\frac{p_{\varphi }^0}{{\varphi }_0{H}^2}\frac{\partial }{\partial y}& 0& \frac{{\lambda }^0+2{\eta }^0}{{\varphi }_0{H}^2}\frac{{\partial }^2}{\partial y^2}& -\frac{p_T^0}{{\varphi }_0{H}^2}\frac{\partial }{\partial y}\\ \frac{1}{{\varphi }_0}\left(\frac{{\mu }^0}{H^2}\frac{{\partial }^2}{\partial y^2}+{\eta }_{\varphi }^0-{\mathcal{D}}_{0,\varphi }^0\right)& \frac{2{\eta }^0}{{\varphi }_0}\frac{\partial }{\partial y}& -\frac{p^0+{\mathcal{D}}_1^0}{{\varphi }_0}\frac{\partial }{\partial y}& \frac{1}{{\varphi }_0}\left(\frac{{\kappa }^0}{H^2}\frac{{\partial }^2}{\partial y^2}+{\eta }_T^0-{\mathcal{D}}_{0,T}^0\right)\end{array}\right].$$

$$u^{\prime }\left(\pm \frac{1}{2}\right)={\frac{\partial T^{\prime }}{\partial y}|}_{y=\pm 1/2}=0.$$

$$\hat{L}\hat{\mathbf{X}}=\omega \hat{\mathbf{X}},\hat{u}\left(\pm \frac{1}{2}\right)=\frac{\mathrm{d}\hat{T}}{\mathrm{d}y}=0,$$

$$(\hat{\varphi },\hat{T}{)}^{\mathsf{T}}$$

$$L_1{\mathbf{X}}_1=\omega {\mathbf{X}}_1,$$

$$L_1=\left[\begin{array}{cccc}0& 0& -\pi \beta {\varphi }_0& 0\\ -\frac{\pi \beta {\eta }_{\varphi }^0}{{\varphi }_0{H}^2}& -\frac{{\pi }^2{\beta }^2{\eta }^0}{{\varphi }_0{H}^2}& -1& -\frac{\pi \beta {\eta }_T^0}{{\varphi }_0{H}^2}\\ \frac{\pi \beta p_{\varphi }^0}{{\varphi }_0{H}^2}& 0& -\frac{{\pi }^2{\beta }^2({\lambda }^0+2{\eta }^0)}{{\varphi }_0{H}^2}& \frac{\pi \beta p_T^0}{{\varphi }_0{H}^2}\\ \frac{1}{{\varphi }_0}\left(-\frac{{\pi }^2{\beta }^2{\mu }^0}{H^2}+{\eta }_{\varphi }^0-{\mathcal{D}}_{0,\varphi }^0\right)& \frac{2\pi \beta {\eta }^0}{{\varphi }_0}& -\frac{\pi \beta (p^0+{\mathcal{D}}_1^0)}{{\varphi }_0}& \frac{1}{{\varphi }_0}\left(-\frac{{\pi }^2{\beta }^2{\kappa }^0}{H^2}+{\eta }_T^0-{\mathcal{D}}_{0,T}^0\right)\end{array}\right].$$

Full text

I Introduction

Granular materials in the so-called rapid flow regime Campbell (1990); Goldhirsch (2003) exhibit various non-uniform structures, such as granular vortices, density waves, clustering, shear banding, etc. Jaeger et al. (1996); Ottino & Khakhar (2000); Aranson & Tsimring (2006); Alam et al. (2009). Among all the non-uniform structures exhibited by granular flows, the shear banding in granular flows—in particular—has received tremendous attention mainly due to its analogy with the shear banding in soft matters, e.g. foams, emulsions, colloidal suspensions, etc. Schall & van Hecke (2010). The shear banding is manifested even in a granular shear flow, one of the simplest types of flows, which serves as a prototype for the rheology and pattern formation Lun et al. (1984); Jenkins & Richman (1985a); Sela & Goldhirsch (1998); Alam & Luding (2003a); Khain & Meerson (2006); Shukla & Alam (2009, 2011a, 2011b). Several experimental (see e.g. Mandl et al. (1977); Veje et al. (1999); Mueth et al. (2000); Viggiani et al. (2001); Bocquet et al. (2001); Conway & Glasser (2004); Liao et al. (2010); Murdoch et al. (2013)) as well as theoretical studies (see e.g. Alam & Nott (1998); Conway & Glasser (2004); Conway et al. (2006); Khain & Meerson (2006); Conway et al. (2006); Khain (2007); Shukla & Alam (2009); Khain (2009); Shukla & Alam (2011b, a); Berzi & Jenkins (2015)) have confirmed the shear-banding phenomenon in granular shear flows. The theoretical studies based on linear and nonlinear stability analyses Alam & Nott (1998); Shukla & Alam (2009, 2011b, 2011a) show that a granular shear flow admits stationary and traveling instabilities leading to clustering and shear banding—with the latter being the focus of the present work. In the phenomenon of the shear banding, a homogeneous flow transforms into an inhomogeneous flow characterized by the coexisting bands of different rheological properties due to the shearing motion. Shearing along the streamwise direction renders inhomogeneities in the other two directions—commonly referred to as the gradient (or transverse) direction and the vorticity (or spanwise) direction. Consequently, shearing along the streamwise direction leads to two different banding instabilities, namely the gradient banding instability and the vorticity banding instability. In the former, bands of high and low shear rate form along the gradient direction while bands of high and low shear stress appear along the vorticity direction in the latter. The localized high and low shear (shear rate or shear stress) regions correspond to low and high density regions, respectively. Both types of the shear-banding instabilities have also been observed in experiments Mandl et al. (1977); Veje et al. (1999); Conway & Glasser (2004); Mueth et al. (2000); Conway et al. (2006) conducted by shearing a granular material in a shear cell wherein shearing remains localized within the narrow regions leaving behind bands of unsheared regions. For more details, the reader is referred to the articles Fielding (2007); Dhont & Briels (2008), which provide comprehensive details of the shear-banding phenomenon.

To investigate such a vast variety of phenomena exhibited by granular flows through hydrodynamic theory is one of the most challenging tasks among the granular community. The hydrodynamic modeling of granular fluids is more involved in comparison to regular fluids since interactions among granular particles are inherently inelastic; this very nature of granular materials poses some undesirable complexities, e.g. microscopic irreversibility, lack of scale separation, mesoscopic flow behaviour, strong nonlinearities in the momentum and energy balance equations Campbell (1990); Goldhirsch (2003). Notwithstanding, hydrodynamic models for dilute/dense granular flows can be derived from the (inelastic) Boltzmann/Enskog–Boltzmann equation within the framework of kinetic theory. Similarly to the well-established Navier–Stokes and Fourier (NSF) equations for regular fluids, a hydrodynamic model for granular fluids consists of the mass, momentum and energy balance equations, and the pressure tensor and heat flux appear as additional unknowns in the momentum and energy balance equations. However, unlike the NSF equations for regular fluids, the energy balance equation in the case of granular fluids contains an additional term, referred to as the collisional dissipation, that accounts for the energy loss due to inelastic collisions among granular particles, and the constitutive relations for the pressure tensor and heat flux are, in general, quite different from the Navier–Stokes’ and Fourier’s laws for regular fluids. For granular fluids, the constitutive relations for the pressure tensor, heat flux and collisional dissipation are typically derived from the Boltzmann equation (in the dilute case) or from the Enskog–Boltzmann equation (in the dense case) by means of the Chapman–Enskog expansion at first-order of expansion, see e.g. Brey et al. (1998); Garzó & Dufty (1999); Sela & Goldhirsch (1998); Brilliantov & Pöschel (2004); Gupta (2011). Alternatively, these constitutive relations can also be derived from the moment equations, see e.g. Jenkins & Richman (1985b, a); Kremer & Marques Jr. (2011); Garzó (2013); Gupta & Shukla (2017); Gupta et al. (2018). It is worthwhile to note that determining the constitutive relations from the moment equations is not only much simpler than that by the Chapman–Enskog expansion performed on the full Boltzmann equation but can also yield more accurate constitutive relations on considering more moments Gupta et al. (2018). The mass, momentum and energy balance equations for granular flows closed with the first-order constitutive relations are referred to as the granular NSF equations. The validity of the granular NSF equations—even for rapid granular flows Campbell (1990); Goldhirsch (2003)—is subjected to the conditions under which the constitutive relations for the pressure tensor, heat flux and collisional dissipation are derived. Hence the constitutive relations involved in the granular NSF equations ought to be chosen carefully, especially while dealing with dense granular flows for which the “molecular chaos” assumption is inadequate.

The papers Jenkins & Richman (1985b, a) by Jenkins and Richman (JR) may be regarded as the pioneering works which derive the balance equations and associated constitutive relations for dense granular gases of identical rough circular disks and spheres, respectively, based on the revised Enskog theory. It is important to note that the restitution coefficient $e$ in the JR model enters only the energy balance equation through the collisional dissipation while the pressure tensor and heat flux do not depend on $e$ due to the approximations made in their theories. Consequently, the NSF transport coefficients—that appear in the constitutive relations for the pressure tensor and heat flux—from the JR model are essentially the same as those for regular (elastic) fluids. Furthermore, the JR model was derived for nearly elastic ($e\approx 1$) granular fluids. Despite these limitations, the JR model has been widely exploited since its derivation and has been validated against particle simulations Bizon et al. (1999) as well as against experiments Rericha et al. (2001) for nearly elastic granular flows. The JR model has also been employed in analyzing various instabilities in granular flows. For instance, Alam et al. (2008) investigated the shear-banding instability in two-dimensional dilute and dense granular shear flows using different variants of the JR model. In Alam et al. (2008), the authors essentially considered different models for the global equation of state (i.e. for pressure) and for the shear viscosity (with and without a viscosity divergence term in its expression) while keeping the other transport coefficients the same as those in the original JR model, and they showed that some of their models underpin the shear banding while the others do not. Moreover, they also established that the onset of the shear-banding instabilities could be predicted with the knowledge of the pressure and shear viscosity of the system. This result indicates that the emergence of the shear banding is linked to the transport properties and, hence, to the choice of constitutive relations. Therefore, the selection of constitutive relations is crucial in describing instabilities and patterns in granular flows.

In another study on the the shear banding, Khain & Meerson (2006) substantiated that the experimentally observed shear-banding instabilities could not be perceived with the usual constitutive relations for dense granular gases; however, shear bands, miraculously, appear on slightly changing the coefficient of the shear viscosity since the shear viscosity diverges at a lower density than the other transport properties. Note that the shear-banding instability predictions of both Alam et al. (2008) and Khain & Meerson (2006) are valid only for the nearly elastic (or quasi-elastic) granular flows. However, in order to understand instability induced patterns in arbitrarily inelastic granular flows correctly, one must utilize proper constitutive relations which incorporate all the microscopic features of granular materials Campbell (1990); Goldhirsch (2003); Garzó (2019).

To overcome the “nearly elastic” limitation of the JR model, Garzó & Dufty (1999) performed the Chapman–Enskog expansion on the Enskog–Boltzmann equation and, in contrast to the JR model, obtained the restitution coefficient-dependent NSF-level constitutive relations for dense granular gases of hard spheres. Subsequently, the results of Garzó & Dufty (1999) were generalized to an arbitrary dimension $d$ by Lutsko (2005). We shall refer to the constitutive relations obtained in Refs. Garzó & Dufty (1999) and Lutsko (2005) by the GDL model. The main differences between the JR and GDL models are as follows: (i) the heat flux in the latter contains an additional term proportional to the density gradient that vanishes identically for elastic particles; this term is absent in the former and (ii) all the transport coefficients depend on the restitution coefficient $e$ in the latter while only the collisional dissipation depends on $e$ in the former. As a consequence, the GDL model is not limited to nearly elastic granular fluids. Recently, Almazán et al. (2013) conducted a comparative study of the Faraday instability in granular flows through the JR and GDL models and through the event-driven simulations, and also concluded that the choice of appropriate constitutive relations is crucial for analyzing granular patterns.

There have been several theoretical studies on the shear-banding instability in a granular uniform shear flow (USF) for quasi-elastic particles, and thus all of them are valid only for nearly elastic granular fluids. Nonetheless, to the best of authors’ knowledge, the shear-banding instability in arbitrarily inelastic dense granular flows has never been addressed theoretically. One of the reasons—among others—is that the granular NSF equations are valid strictly for processes involving small spatial gradients, and moreover since the (reduced) shear rate in the granular USF is inversely proportional to the collision frequency, small spatial gradients (the validity region of the NSF equations) would again mean large restitution coefficient (or nearly elastic particles) Lutsko (2006); Garzó (2006). Furthermore, the granular USF state is inherently anisotropic and hence requires the generalized transport coefficients for predicting the instabilities in the USF correctly Garzó (2006). The generalized transport coefficients depend on the shear rate and have tensorial form. For dilute granular flows, the generalized transport coefficients were independently obtained by Lutsko (2006) and Garzó (2006), and have been utilized by Garzó (2006) to analyze the linear stability of the USF. From his linear stability analysis, Garzó (2006) showed significant discrepancies between the results obtained with the generalized transport coefficients and usual NSF transport coefficients; nonetheless, the unavailability of the numerical/experimental results thwarted quantitative comparisons in Garzó (2006). On the other hand, the derivation of such generalized transport coefficients in the case of dense granular flows is extremely involved and is beyond the scope of the present paper. Hence the NSF equations along with the GDL model Garzó & Dufty (1999); Lutsko (2005) for transport coefficients, which is the best known NSF-level hydrodynamic model for describing dense granular flows, can be considered as an intermediate way of investigating the linear stability of the USF.

In this context, the goal of the present work is to analyze the shear-banding instability in arbitrarily inelastic dense granular flows of hard disks by using the GDL constitutive relations. By means of the linear stability analysis of the USF, the onset of the shear-banding instability is predicted. In contrast to previous studies Khain & Meerson (2006); Alam et al. (2008), the present stability results are valid for dilute-to-dense arbitrarily inelastic particles with the only assumption that the restitution coefficient is constant, i.e. it does not depend on the impact velocity. It is emphasized that the focus of the present work is to determine the control parameters for the onset of the shear-banding instability but not the shape and location of a shear band, for which a theory, based on the principle of minimum energy dissipation, developed e.g. in Unger et al. (2004); Török et al. (2007) may be needed.

The rest of the paper is organized as follows. The problem description and the governing equations are presented in Sec. II. The non-dimensionalization and base state flow whose stability is to be investigated are demonstrated in Sec. III. The linear stability of the base state flow is analyzed in Sec. IV. The results and discussion are elucidated in Sec. V. The paper ends with conclusions and outlook in Sec. VI.

II Granular hydrodynamic equations

A granular flow of mono-disperse smooth identical inelastic hard disks of diameter $d_{p}$ can be described by the mass, momentum and energy balance equations, which read Brilliantov & Pöschel (2004); Garzó (2019)

[TABLE]

Here, $\bar{\rho}=\rho_{p}\,\phi$ is the mass density with $\rho_{p}$ being the material density and $\phi$ being the volume fraction of grains; $\bar{\bm{U}}=(\bar{u},\bar{v},\bar{w})$ is the coarse-grained velocity with $\bar{u}$, $\bar{v}$ and $\bar{w}$ being its components in the $\bar{x}$-, $\bar{y}$- and $\bar{z}$-directions, respectively; $\bar{T}$ is the granular temperature; $\bar{\bm{P}}$ is the pressure tensor; $\bar{\bm{q}}$ is the (granular) heat flux; $\bar{\mathcal{D}}$ is the collisional dissipation due to inelastic collisions among grains; and $d$ denotes the dimension of the problem which takes value two for hard-disk flows and three for hard-sphere flows.

Clearly, system (1) is not closed due to the presence of the additional unknowns: $\bar{\bm{P}}$, $\bar{\bm{q}}$ and $\bar{\mathcal{D}}$. These unknowns are typically expressed in terms of the hydrodynamic variables $\bar{\rho}$, $\bar{\bm{U}}$, $\bar{T}$ and their spatial gradients by means of the Chapman–Enskog expansion, see e.g. Goldshtein & Shapiro (1995); Brey et al. (1998); Sela & Goldhirsch (1998); Garzó & Dufty (1999); Brilliantov & Pöschel (2004); Gupta (2011). To first order in spatial gradients of the hydrodynamic variables, these unknowns are expressed as Goldshtein & Shapiro (1995); Sela & Goldhirsch (1998)

[TABLE]

where $\bar{\lambda}=\bar{\gamma}-\bar{\eta}$; $\bar{\bm{I}}$ is the identity tensor; and \bar{\bm{S}}=\frac{1}{2}\big{[}\bar{\bm{\nabla}}\bar{\bm{U}}+(\bar{\bm{\nabla}}\bar{\bm{U}})^{\mathsf{T}}\big{]}-(\bar{\bm{\nabla}}\cdot\bar{\bm{U}})\,\bar{\bm{I}} is the deviatoric strain rate tensor; the quantities $\bar{p}$, $\bar{\eta}$, $\bar{\gamma}$, $\bar{\kappa}$, $\bar{\mu}$, $\bar{\mathcal{D}}_{0}$ and $\bar{\mathcal{D}}_{1}$ are the pressure, shear viscosity, bulk viscosity, pseudo-thermal conductivity, Dufour-like coefficient (which identically vanishes for ordinary fluids), zeroth- and first-order contributions to the collisional dissipation, respectively, and are given in the form of constitutive relations:

[TABLE]

Here all the $f_{i}$’s are the dimensionless functions of the volume fraction $\phi$ and restitution coefficient $e$ only. It is worthwhile to note that the values of $f_{i}$’s are different for different models for the constitutive relations. From the GDL model, $f_{i}$’s in the case of hard-disk flows ($d=2$) read Lutsko (2005); Garzó (2019)

[TABLE]

where $c$, the fourth cumulant estimating the lowest order correction to the Gaussian distribution function, is given by van Noije & Ernst (1998)

[TABLE]

and $G(\phi)=\phi\,\chi(\phi)$ with $\chi(\phi)$ being the pair correlation function adopted from Torquato (1995):

[TABLE]

In (8), $\phi_{f}$ and $\phi_{c}$ are the freezing packing fraction and the random close-packing fraction, respectively. For hard-disk flows ($d=2$), $\phi_{f}\approx 0.69$ and $\phi_{c}\approx 0.82$ Torquato (1995). The other variables in (6) for hard-disk flows ($d=2$) are given by Garzó (2012)

[TABLE]

The quantities in (9) emanate from the so-called modified Sonine approximation introduced by Garzó et al. (2007), and are presented in a more coherent form—for an arbitrary dimension $d$—by Garzó (2012) (see also the textbook Garzó (2019)). All the quantities, except $\mu_{k}^{\ast}$, in (9) are also given in Almazán et al. (2013) for hard-disk flows ($d=2$). The $\mu_{k}^{\ast}$ in the present work is twice of that of  Almazán et al. (2013) but is the same as that in Garzó (2012) for $d=2$ in order to keep the standard form of the reduced Dufour-like coefficient ($\mu^{\ast}=n\mu/(\kappa_{0}T)$) given, e.g., in Brey et al. (1998); Garzó (2012); Gupta et al. (2018); Garzó (2019).

Notably, the physical properties in a granular flow are transported via two mechanisms, namely the kinetic and collisional. The former is attributed to streaming, i.e. to the movement of particles from one place to another, while the latter to collisions among grains. Needless to say, the former is dominant in dilute flows while the latter in dense flows. Owing to these transport mechanisms, the pressure tensor $\bar{\bm{P}}$ and heat flux $\bar{\bm{q}}$ for a (dense) granular flow can be decomposed into their kinetic and collisional contributions, i.e. $\bar{\bm{P}}=\bar{\bm{P}}^{k}+\bar{\bm{P}}^{c}$ and $\bar{\bm{q}}=\bar{\bm{q}}^{k}+\bar{\bm{q}}^{c}$, where the superscripts ‘$k$’ and ‘$c$’ denote the kinetic and collisional contributions, respectively, see e.g. Garzó & Dufty (1999); Lutsko (2005); Garzó et al. (2007) and references therein. Accordingly, the pressure $\bar{p}$ and the transport coefficients $\bar{\eta},\bar{\gamma},\bar{\kappa},\bar{\mu}$ and, hence the dimensionless functions $f_{i}$’s for $i\in\{1,2,3,4,4h\}$, can be decomposed into their kinetic and collisional parts, i.e. $f_{i}=f_{i}^{k}+f_{i}^{c}$. From the expressions of $\bar{\bm{P}}^{k}$, $\bar{\bm{P}}^{c}$, $\bar{\bm{q}}^{k}$ and $\bar{\bm{q}}^{c}$ given in Lutsko (2005); Garzó (2012), it is straightforward to determine the kinetic and collisional parts of $f_{i}$’s, which read

[TABLE]

and

[TABLE]

Indeed, as expected, the collisional contributions to the pressure tensor ($\bar{\bm{P}}^{c}$) and heat flux ($\bar{\bm{q}}^{c}$) for dilute granular flows ($\phi\to 0$) vanish and hence the pressure tensor and heat flux for dilute granular flows are given by $\bar{\bm{P}}=\bar{\bm{P}}^{k}$ and $\bar{\bm{q}}=\bar{\bm{q}}^{k}$ Garzó (2013, 2019). Consequently, $f_{i}^{k}$’s in (10) can be referred to as the dilute limit of $f_{i}$’s for $i\in\{1,2,3,4,4h\}$. The collisional dissipation $\bar{\mathcal{D}}$ in a (dilute or dense) granular flow, on the other hand, is attributed only to inelastic collisions among grains. Interestingly, the first-order contribution (in spatial gradients) to the collisional dissipation, $\bar{\mathcal{D}}_{1}$, is zero for dilute granular flows Garzó (2013, 2019), and the zeroth-order contribution (in spatial gradients) to the collisional dissipation, $\bar{\mathcal{D}}_{0}$, for dilute granular flows can be obtained from its dense counterpart given, e.g., in Lutsko (2005); Garzó (2012, 2013, 2019) by taking $\chi(\phi\to 0)=1$. Consequently, the dimensionless functions $f_{5}$ and $f_{5u}$ for dilute granular hard-disk flow ($d=2$) are given by

[TABLE]

In what follows, the constitutive relations (5) obtained using $f_{i}\approx f_{i}^{k}$ for $i\in\{1,2,3,4,4h\}$ with $f_{i}^{k}$ from (10) and $f_{5},f_{5u}$ from (12) will be referred to as the dilute limit of the GDL model. Similarly, the constitutive relations (5) obtained using $f_{i}\approx f_{i}^{c}$ for $i\in\{1,2,3,4,4h\}$ with $f_{i}^{c}$ from (11) and $f_{5},f_{5u}$ from (6) will be referred to as the collisional limit of the GDL model.

III Non-dimensionalization and base state

We shall investigate a plane shear flow of granular hard disks confined in a two-dimensional channel of width $h$. The flow is driven by the two oppositely moving walls (of the channel) placed at $\bar{y}=\pm h/2$ with a speed $U_{w}/2$ along the streamwise ($\bar{x}$) direction; hence the overall shear rate is $U_{w}/h$, see figure 1.

It is worthwhile to note that the velocity slip at a rigid wall is inherently present in granular flows. This has also been verified through experiments and computer simulations Mueth et al. (2000); Conway & Glasser (2004); Conway et al. (2006). The velocity slip, in turn, leads to the generation of the pseudo-thermal energy at the wall that competes with the energy lost due to particle-wall collisions. Consequently, energy flux at the wall is, in general, non zero. Nevertheless, for simplicity, we choose the no-slip and zero heat flux (adiabatic) boundary conditions for the velocity and temperature, respectively. These idealized boundary conditions have also been used in previous studies, e.g. in Wang et al. (1996); Alam & Nott (1998); Alam (2006); Alam et al. (2008); Khain (2007), pertaining to the stability of granular shear flows. The no-slip and adiabatic boundary conditions for the problem under consideration read

[TABLE]

For the purpose of non-dimensionalization, we choose $h$ as a reference length, $U_{w}$ as a reference velocity, the inverse of the shear rate $h/U_{w}$ as a reference time and $\rho_{p}$ as the reference mass density. In the following, the quantities with overbar are dimensional and their bare counterparts are dimensionless. Let us define the dimensionless quantities without bars as follows,

[TABLE]

where $H=h/d$ is referred to as the dimensionless gap between the walls or the channel width. With $\rho_{p}$ as the scaling for the density, the volume fraction $\phi$ also denotes dimensionless density. With these scales, governing equations (1) along with the constitutive relations (5) in the dimensionless form read

[TABLE]

where

[TABLE]

Boundary conditions (13) in the dimensionless form read

[TABLE]

We shall be investigating the shear banding in a two-dimensional plane shear flow. Specifically, we focus on the gradient banding due to which bands of dense and dilute regions form along the gradient (i.e. $y$-) direction. Such instability arises from the perturbations having no variation in the streamwise (i.e. $x$-) direction. The streamwise independent (i.e. $\partial(\cdot)/\partial x=0$) governing equations (15) in the dimensionless form read

[TABLE]

III.1 Base state: Uniform shear flow

The basic flow, whose stability is to be analyzed here, is assumed to be a steady $(\partial(\cdot)/\partial t\equiv 0)$, fully developed $(\partial(\cdot)/\partial x\equiv 0)$, plane shear flow of the following form

[TABLE]

Hence, the mass balance equation (18)1 is identically satisfied for this flow, while the remaining equations in (18) reduce to

[TABLE]

Using (16), the above equations, respectively, lead to

[TABLE]

where $c_{1}$ and $c_{2}$ are the integration constants. The above set of equations with no-slip and zero heat flux boundary conditions (17) admits the following base state solution:

[TABLE]

where the superscript ‘0’ denotes the base state solutions, $\phi_{0}$ and $T_{0}$ are constants, and $f_{2}^{0}=f_{2}(\phi^{0},e)$ and $f_{5}^{0}=f_{5}(\phi^{0},e)$. Note that the base flow velocity is linear with constant density and constant temperature. Such a base state gives constant or uniform shear rate, i.e. $u^{0}_{y}=1$, thus leading to the USF. The base state density $\phi_{0}$, the channel width $H$ and the restitution coefficient $e$ are the control parameters for the problem.

IV Linear stability analysis

IV.1 Perturbation equations

For the linear stability analysis of the USF, the field variables $(\phi,u,v,T)$ are decomposed into the base state solution plus the perturbation from the base state solution as follows:

[TABLE]

Here, the field variables with prime denote the perturbations from their respective base states. These perturbations are assumed to be small so that the linear theory remains valid. Substituting the field variables from (25) into governing equations (18) and retaining only the linear terms of the perturbed field variables, one obtains the governing equations for the perturbed field variables:

[TABLE]

where the subscripts ‘$\phi$’ and ‘$T$’ represent the partial derivatives with respect to $\phi$ and $T$, respectively, and the superscript zero represents the variables calculated at the base state. The above equations can be written as a matrix system

[TABLE]

where $\bm{X}=(\phi^{\prime},u^{\prime},v^{\prime},T^{\prime})^{\mathsf{T}}$ is the vector of perturbed fields and $\mathcal{L}$ is the matrix of linear differential operators given by

[TABLE]

From (17), the boundary conditions for the perturbed field variables read

[TABLE]

IV.2 Analytical solutions

We perform the standard linear stability analysis on the USF by assuming a normal mode solution of the form $\bm{X}(t,y)=\hat{\bm{X}}(y)\exp(\omega t)$, where $\hat{\bm{X}}=(\hat{\phi},\hat{u},\hat{v},\hat{T})^{\mathsf{T}}$ and $\omega=\omega_{r}+\mathsf{i}\,\omega_{i}$ is the complex frequency whose real part $\omega_{r}$ represents growth or decay rate of the perturbations and imaginary part $\omega_{i}$ denotes the oscillation of the perturbation, with $\mathsf{i}$ being the imaginary unit. Substituting this normal mode solution into the linearized perturbation equations (27)–(29), we get the following matrix eigenvalue problem:

[TABLE]

where $\hat{L}=\mathcal{L}\left(\partial/\partial y\rightarrow{\rm d}/{\rm d}y,\partial^{2}/\partial y^{2}\rightarrow{\rm d}^{2}/{\rm d}y^{2}\right)$. It has been verified that eigenvalue problem (30) has an analytical solution in terms of sine and cosine functions Alam & Nott (1998), as

[TABLE]

where $\beta=1,2,3,\dots$ are the mode numbers and $(\phi_{1},u_{1},v_{1},T_{1})^{\mathsf{T}}$ is the constant amplitude of the normal mode solution. For instance, $\beta=1$ is the fundamental mode and $\beta=2$ is the second harmonic of the normal mode solution, etc. With solution (31), problem (30)1 simplifies to

[TABLE]

where $\bm{X}_{1}=({\phi}_{1},{u}_{1},{v}_{1},{T}_{1})^{\mathsf{T}}$ and

[TABLE]

For the nontrivial solutions of (32) ${\rm det}(L_{1}-\omega I)=0$, where $I$ is the identity matrix of size 4. This condition is the dispersion relation and can be written as

[TABLE]

where

[TABLE]

Here the coefficients $a_{ij}$ are the functions of the transport coefficients evaluated at the base state. However, their explicit expressions are relegated to appendix A for better readability. Dispersion relation (34) is a fourth-degree polynomial in $\omega$ with real coefficients and, therefore, there are three possibilities for four roots of (34): (i) all roots are real, (ii) two complex conjugate pairs of roots and (iii) two real roots and one complex conjugate pair of roots.

IV.3 Asymptotic analysis

With the help of the classical asymptotic analysis in powers of $H^{-1}$ with $H^{-1}\to 0$, one can find the analytical expressions of the eigenvalues (the roots of (34)), as follows. Let the frequency $\omega$ be represented by

[TABLE]

where $\omega_{0},\omega_{1},\omega_{2},\dots$ are unknown coefficients. Substituting this ansatz along with (35) into (34), and comparing each power of $H$ on both sides of the resulting equations, one obtains algebraic equations, which are solved for the unknowns $\omega_{i}$’s in (36). Exploiting these values, one obtains four roots $\omega=\omega^{(1,2,3,4)}$ from (34), which are given by

[TABLE]

where the subscripts ‘$r$’ and ‘$i$’ represent the real and imaginary parts, respectively, of the roots and

[TABLE]

It is evident from (37)–(39) that in the limit of large $H$, dispersion relation (34) has two real roots $\omega^{(1,2)}$ and a complex conjugate pair of roots $\omega^{(3,4)}$. In the limit of large $H$, it is verified numerically that $\omega^{(2)}$ and $\omega_{r}^{(3,4)}$ are always negative resulting into the least stable mode as $\omega^{(1)}$. Figure 2 illustrates the four eigenvalues for large $H$ obtained through the asymptotic analysis of dispersion relation (34) (solid line) and those obtained by solving (32) numerically (symbols) for $\phi^{0}=0.6$, $e=0.5$ and $\beta=1$. It can be seen from the figure that the eigenvalues from both the methods are in excellent agreement for large $H$.

V Results

The linear stability of the USF of granular materials has been studied previously for the case of nearly elastic particles Alam & Nott (1998); Alam (2005); Alam et al. (2005); Alam (2006); Alam et al. (2008). In this section, we analyze the linear stability of the USF of arbitrary inelastic granular hard disks by appropriately choosing the GDL constitutive relations pertaining to dense granular flows Garzó & Dufty (1999); Lutsko (2005). In order to perform the linear stability analysis, eigenvalue problem (30) is solved numerically using the Chebyshev spectral collocation method. In addition, eigenvalue problem (32), which is obtained using the exact solutions (31), is also solved analytically. It is verified that the eigenvalues obtained numerically using the Chebyshev spectral collocation method and analytically using the exact solutions (31) are found in an excellent agreement (figure is not shown for brevity). Therefore in the present analysis, the eigenvalues are computed by solving (32). All the computations have been performed in Matlab®. For analyzing the stability results, we define the least stable eigenvalue $\omega^{l}$ as one of the eigenvalues whose real part is maximum for a fixed mode number, and the dominant eigenvalue $\omega^{d}$ is one of the least stable eigenvalues whose real part is maximum over all the mode numbers $\beta$, i.e.

[TABLE]

In order to avoid the mode number and grid dependencies on the stability predictions, the dominant modes are analyzed in the present work. The contours of positive, negative and zero dominant growth rates, $\omega^{d}_{r}$, are shown in the $(H,\phi^{0})$-plane for the inelastic particles with the restitution coefficient $e=0.6$, see figure 3(a). The real part of the dominant mode is positive (i.e. $\omega^{d}_{r}>0$) inside the zero contour and negative (i.e. $\omega^{d}_{r}<0$) outside the zero contour, therefore the flow is unstable inside the zero contour and stable outside. It is also verified that the instability depicted in figure 3(a) is due to the stationary waves since the imaginary part of the complex frequency $\omega^{d}$ is always zero. This can also be seen from the asymptotic analysis presented in  Sec. IV.3 that the least stable eigenvalue $\omega^{l}$ is always real. Therefore the shear-banding instabilities in a granular shear flow are due to the stationary waves.

We have also examined that the lowest mode number, i.e. $\beta=1$, is the first one to become unstable for a fixed restitution coefficient. This is similar to the classical Rayleigh–Bénard convection in which the first mode is the dominant mode Chandrasekhar (1961). Figure 3(a) also illustrates that for a fixed restitution coefficient there exist a critical channel width $H_{c}(\phi^{0})$ and a critical mean density $\phi^{0}_{c}(H)$ above which the USF is unstable and below which it is stable. That is to say, the USF becomes unstable if either $H>H_{c}(\phi^{0})$ or $\phi^{0}>\phi_{c}^{0}(H)$. In addition, it can also be seen from figure 3(a) that the USF remains stable for all densities and for all channel widths below an onset value of the mean density, say $\phi^{0}_{onset}$ (see Sec. V.1 for more details). Note that the onset mean density $\phi^{0}_{onset}$ depends only on the restitution coefficient. In particular, for $e=0.6$, the flow is stable when $\phi^{0}<\phi^{0}_{onset}\approx 0.38$ for all values of $H$. Figure 3(a) also reveals that the USF is stable in the dilute limit ($\phi^{0}\to 0$).

In order to get more insight, the variation of the dominant growth rate with the channel width for fixed values of the mean density and restitution coefficient is illustrated in figure 3(b). The kinks (crests) in the figure correspond to the eigenmode crossing from the mode number $j$ to $j+1$, where $j=1,2,3,\dots$ is a positive integer. It is seen that the dominant mode number increases with the channel width, i.e. $\beta=1$ mode is the first one to lose stability and remains unstable with increasing channel width until it crosses $\beta=2$ mode (first kink in figure 3(b)), thereafter $\beta=2$ mode becomes dominant mode until it crosses the next $\beta=3$ mode (second kink in figure 3(b)), and so on. For parameter values shown in figure 3(b), the USF becomes unstable for $H>H_{c}\approx 10.446$.

Figure 4 illustrates the density, temperature and velocity eigenfunctions for three values of the channel width, $H=20$, $40$, $50$, with other parameters being the same as those in figure 3(b). The eigenfunctions displayed in the first row correspond to $\beta=1$ mode for which $(\hat{\phi},\hat{T})\propto\cos(\pi y\pm\pi/2)=\sin(\pi y)$ that vanish once at $y=0$ in the flow domain, therefore the density and temperature solutions attain their mean values $\phi^{0}$ and $T^{0}$ once in the flow domain. The second row of figure 4 corresponds to $\beta=2$ mode which gives eigenfunctions $(\hat{\phi},\hat{T})\propto\cos(2y\pi\pm\pi)=\cos(2\pi y)$ that vanish twice along the flow domain at $y=\pm 1/4$. Similarly, the third row corresponds to $\beta=3$ mode having the eigenfunctions $(\hat{\phi},\hat{T})\propto\cos(3y\pi\pm 3\pi/2)=\pm\sin(3\pi y)$ that vanish thrice at $y=0,~{}\pm 1/3$. Consequently, the density and temperature solutions attain its mean values once, twice and thrice in the flow domain for $H=20$, $40$ and $60$, respectively. By analyzing the corresponding velocity components (second column in figure 4), it is evident that while the horizontal component of the velocity eigenfunction $\hat{u}$ varies significantly, its transverse component $\hat{v}$ varies only slightly.

Figure 5 exhibits the neutral stability curve obtained with the dilute limit of the GDL model (i.e. by using $f_{i}\approx f_{i}^{k}$ for $i\in\{1,2,3,4,4h\}$ with $f_{i}^{k}$ from (10) and $f_{5},f_{5u}$ from (12) in the constitutive relations (5)). It is evident from the figure that the USF is always stable in the dilute limit ($\phi^{0}\to 0$), which is also known from the previous studies Alam et al. (2008); Conway & Glasser (2004); Conway et al. (2006). Therefore the dilute limit of the GDL model is able to capture the stability of dilute granular shear flows. On the other hand, although the collisional mechanism is dominant over the kinetic mechanism in dense granular flows, the collisional limit of the GDL model (the constitutive relations (5) obtained using $f_{i}\approx f_{i}^{c}$ for $i\in\{1,2,3,4,4h\}$ with $f_{i}^{c}$ from (11) and $f_{5},f_{5u}$ from (6)) alone fails to capture the shear-banding instability in dense granular shear flows Alam & Luding (2003b); Conway & Glasser (2004); Conway et al. (2006) correctly (phase diagram is not shown for brevity): the collisional limit of the GDL model predicts that the USF is stable for high densities whereas the full GDL model (i.e. the constitutive relations with (6)) (see figure 3) and molecular dynamics simulations Alam & Luding (2003b); Conway & Glasser (2004); Conway et al. (2006) predict the shear-banding instability for dense granular shear flows. By analyzing three variants of the GDL model viz. the full model, its dilute limit and its collisional limit, one can conclude that the choice of constitutive relations plays an important role in determining the shear-banding instability. For correct instability predictions in dense granular flows, both the kinetic and collisional mechanisms are important and hence none of them should be discarded. Therefore, in the following, we focus on the effect of the inelasticity on the shear-banding instability through the full GDL model.

Figure 6 compares qualitatively a typical result from molecular dynamics simulation Conway & Glasser (2004) with that from the present linear stability theory for parameter values $e=0.8$, $H=80$ and $\phi^{0}=0.6$. The left panel in the figure shows two parallel, high-density regions located on both sides of the centerline Conway & Glasser (2004). The right panel in the figure exhibits the density eigenfunction $\hat{\phi}$ for the same parameter values and for $\beta=2$. Similarly to the molecular dynamics simulation result Conway & Glasser (2004), the density eigenfunction shows minimum density at the centerline ($y=0$) of the channel and higher densities on both sides of the centerline. For these parameter values, the corresponding growth rate is positive showing the instability of the USF.

V.1

Effect of the restitution coefficient and channel width on the shear banding: a global criterion

Figure 7(a) illustrates the neutral stability curves in the $(H,\phi^{0})$-plane for various restitution coefficients ranging from moderately inelastic to quasi-elastic limit. For a fixed restitution coefficient, the USF is unstable inside each contour and stable outside. Similarly to figure 3(a), figure 7(a) also shows that the USF is stable in the dilute limit ($\phi^{0}\to 0$) for all values of the restitution coefficient, which is in agreement with previous studies, e.g. Alam & Nott (1998); Garzó (2006); Alam et al. (2008). We see that the neutral stability curve shifts towards right as the restitution coefficient increases thereby leading to more stable region behind it. In other words, the range of the channel width for which the USF is unstable decreases with increasing the restitution coefficient. Indeed, there exists a set of critical parameters ($\phi^{0}_{c},H_{c}$) corresponding to the boundary of a neutral stability curve outside which the USF becomes unstable.

As discussed earlier, the lower branch of each neutral stability curve in figure 7(a) asymptotically (as $H\rightarrow\infty$) approaches to a minimum mean density below which the flow is always stable. We define this minimum critical density as the onset density, i.e.

[TABLE]

where the USF is stable for $\phi^{0}<\phi^{0}_{onset}$ and vice versa. The onset mean density (43) as a function of the restitution coefficient is plotted in figure 7(b). It is clear that the onset density decreases with increasing the restitution coefficient, which implies that the USF tends to lose stability at a lower density in the elastic limit than that in the inelastic case. Note that the above conclusion holds for very large channel widths.

It can also be noticed from figure 7(a) that there exists an onset value of the channel width—associated with the nose of each neutral stability curve—below which the USF is always stable for all mean densities and above which it is unstable. Let us define this onset value of the channel width as the minimum of all critical channel widths:

[TABLE]

Note that $\phi_{c}^{0}$ and $H_{c}$ depend on both the mean density and restitution coefficient whereas $\phi_{onset}^{0}$ and $H_{onset}$ are only functions of the restitution coefficient and hence of the inelasticity. It is evident from figure 7(a) that the value of $H_{onset}$ increases with increasing the restitution coefficient implying that the USF is more stable as $e$ increases which qualitatively matches with the simulation results of Conway & Glasser (2004); Conway et al. (2006). Owing to the functional dependence (see (44)), the shear banding in granular shear flow starts to appear at small channel widths when particles are more inelastic as compared to the quasi-elastic ones. In contrast, the granular shear flow remains uniform for large values of the channel widths.

In order to get further insight, the neutral stability curves in the $(e,\phi^{0})$-plane are shown for three values of the channel width in figure 8(a). The flow is stable outside (towards right) of each contour and unstable inside the bounded region of the curve (towards left). As the channel width increases, the neutral stability curve shifts towards right such that it covers more unstable region in the $(e,\phi^{0})$-plane, i.e. the flow becomes more unstable with the increasing channel width. Therefore for a fixed density, the range of the restitution coefficient for which the USF is unstable increases with increasing channel width. Figure 8(a) also depicts that there is an onset restitution coefficient below which the shear-banding instability persists and above which the USF remains stable. We define this onset restitution coefficient as:

[TABLE]

It can also be seen from the inset of figure 8(b) that $e_{onset}$ increases with increasing the channel width.

Although the neutral stability curves shown in figures 7 and 8 depict the overall behavior of the instability, it still remains to understand how the onset parameters $H_{onset}(e)$ and $e_{onset}(H)$ vary and to discern if there exists any relation relating these parameters. We shall now seek the onset of the shear-banding instability in terms of these onset parameters. As mentioned earlier, the USF is unstable for all $H>H_{onset}$ and for all $e<e_{onset}$, and these onset values correspond to the nose of each neutral stability curve, see figure 8(a). The onset parameters $H_{onset}(e)$ (main panel) and $e_{onset}(H)$ (inset) are shown for some points by circles in figure 8(b). By curve fitting, one can find a functional relationship between the onset parameters $e_{onset}$ and $H_{onset}$, which reads

[TABLE]

where $\alpha=115$, $\gamma=17$ and $\delta=12.86$. The solid line in figure 8(b) represents the values obtained from relation (46). It is important to note that (46) is a global criterion as this does not depend on the spatial positions.

Let us now analyze the effect of the mean density on the shear-banding instability. Figure 9 illustrates the neutral stability curves in the $(H,e)$-plane for various mean densities. For a fixed mean density, the USF is unstable inside (towards right of) each of the neutral stability curves. The neutral stability curves for less dense flows look markedly different from those of moderately-to-highly dense flows. It is also seen that the unstable region increases with increasing the mean density, therefore the shear-banding instabilities are more prone to the dense flows, in general. However, for densities $\phi^{0}>\phi_{f}\approx 0.69$, the instability region decreases with increasing densities, see figure 9.

V.2

Effect of the restitution coefficient and density on the shear banding: another global criterion

As discussed above, the shear-banding instability is stationary which implies that the least stable mode, and hence the dominant mode, is real. Therefore the onset of the shear-banding instability, where the growth rate is zero, can be determined analytically from dispersion relation (34) by substituting $\omega=0$ which gives $a_{0}=0$. Using the expression of $a_{0}$ from (35) (see appendix A), one obtains

[TABLE]

where

[TABLE]

Thus at the onset of the shear-banding instability, equality (47) must be satisfied. Note that the left-hand side of (47) is always positive and therefore

[TABLE]

because $\xi_{2}^{0}>0$. Following Alam (2006); Alam et al. (2008), the condition $\xi_{1}^{0}>0$ is equivalent to

[TABLE]

which must be satisfied at the onset of instability. Thus condition (50) leads to a necessary criterion for the shear-banding instability. It is worth noticing that the term $\sqrt{f_{2}^{0}f_{5}^{0}}/f_{1}^{0}$ is none other than the ratio of the shear stress to the pressure of the USF. This further allows us to express the dynamic friction coefficient, the ratio of sliding force $F_{s}$ to normal force $F_{n}$, as

[TABLE]

Here we have used the definition of the base state granular temperature $T^{0}=f_{2}^{0}/f_{5}^{0}$ and of the base state shear rate $u^{0}_{y}=1$, see Sec. III.1. From (50) and (51), it can be concluded that the existence of the shear-banding instability requires the dynamic friction coefficient $C_{f}$ to be an increasing function of mean density $\phi^{0}$, i.e. $\partial C_{f}/\partial\phi^{0}>0$. In other words, the condition for the onset of the shear-banding instability is given by

[TABLE]

Note that condition (52) is also a global criterion for the onset of the shear-banding instability as it also does not depend on the spatial positions.

Figure 10 shows the variation of the dynamic friction coefficient $C_{f}$ and its gradient with respect to the mean density $\phi^{0}$ for various restitution coefficients. It is seen from figure 10(a) that the dynamic friction coefficient $C_{f}$ varies non-monotonically—it first decreases, attains a minimum, and increases thereafter with increasing density. A value of the mean density $\phi^{0}$ where $C_{f}$ attains its minimum is the same as the onset density $\phi_{onset}^{0}$ in figure 8(a) (for corresponding $e$), which was extracted from the neutral stability curve. Clearly, at $\phi^{0}=\phi_{onset}^{0}$, the slope $\partial C_{f}/\partial\phi^{0}$ is zero, see figure 10(b).

Similarly to figure 7(b), figure 11 also depicts the variation of the onset mean density $\phi_{onset}^{0}$ with the restitution coefficient calculated using criterion (52) (solid line). Clearly, the onset mean density decreases monotonically with increasing $e$, which implies that the onset of the shear-banding instability in the elastic limit ($e\to 1$) occurs at slightly lower mean density than in the case of restitution coefficient $e<1$. By curve fitting, one can find a functional relationship between $\phi_{onset}^{0}$ and $e$, which reads

[TABLE]

where $\alpha_{1}=-0.1350$, $\alpha_{2}=0.2144$, $\alpha_{3}=-0.2208$ and $\alpha_{4}=0.4768$. This implies that criterion (52) for the onset of the shear-banding instability is equivalent to relation (53).

VI Conclusions and outlook

In this paper, we have extended a previous study Alam et al. (2008) and investigated the shear-banding instability in an arbitrary inelastic granular shear flow. In particular, we have analyzed the stability of the USF using granular hydrodynamic equations closed with the Navier–Stokes-level constitutive relations proposed by Garzó & Dufty (1999) and Lutsko (2005). Two limiting cases of the full GDL model, namely the kinetic limit, which is valid for dilute flows, and the collisional limit have also been discussed. It has been shown that the kinetic limit of the GDL model is able to give correct prediction about the stability of dilute granular shear flows. Surprisingly, the collisional limit alone fails to give meaningful prediction on the instability of dense granular shear flows. Thus we have shown that both the kinetic and collisional mechanisms are important to capture the shear-banding instability in dense granular shear flows correctly. Furthermore, by analyzing the neutral stability curves for various parameters in different planes, we have found that the USF is always stable in the dilute limit and therefore non-uniformity in terms of shear band does not appear in dilute granular flows. In contrast to dilute granular shear flows, moderately-to-highly dense granular shear flows become unstable once the control parameters exceed their critical values.

The influence of the restitution coefficient, channel width and mean density on the shear banding has also been explored. It has been found that the USF is more unstable with: (i) decreasing the restitution coefficient (or increasing the inelasticity) for any fixed mean density and fixed channel width, (ii) increasing the channel width for any fixed mean density and fixed restitution coefficient, and (iii) increasing the mean density (in general) for any fixed channel width and fixed restitution coefficient. These findings all together lead to the fact that the shear banding is more prone to dense flows of highly inelastic particles confined in channels having large channel widths. This fact agrees qualitatively with the simulation results of Conway & Glasser (2004); Conway et al. (2006), who showed that a pseudo one-dimensional cluster appears with increasing the restitution coefficient or with increasing the channel width while fixing the other parameters, and therefore the intensity of the clustering instability can be controlled with these parameters.

In the present study, the onset values of the density and channel width as a function of the restitution coefficient have been assessed. For this, we have defined three parameters: (i) $\phi^{0}_{onset}(e)$, the smallest mean density below which the USF is always stable for all channel widths and above which it is unstable, (ii) $H_{onset}(e)$, the critical channel width below which the flow is stable for all mean densities and above which it becomes unstable and (iii) $e_{onset}(H)$, the critical restitution coefficient above which the flow is always stable for all mean densities and below which it is unstable. Furthermore, a relation between $H_{onset}$ and $e_{onset}$ has been obtained in a power-law form. From this relation, one can easily find the onset channel width for the shear-banding instability at a fixed restitution coefficient and vice versa. This is a global criterion (as it does not depend on spatial locations) for the onset of the shear banding in terms of $e$ and $H$—the shear banding appears for $H>H_{onset}(e)$ or for $e<e_{onset}(H)$.

It has been found that the onset of the shear-banding instability in granular USF is tied to the increasing dynamic friction coefficient, $C_{f}(\phi^{0},e)=\eta^{0}u_{y}^{0}/p^{0}$ with $\eta^{0}$ and $p^{0}$ being the shear viscosity and pressure, respectively, in the USF. In particular, the USF breaks into the dense and dilute regions of low and high shear (shear stress or shear rate) along the gradient direction when the dynamic friction coefficient increases with the mean density. In other words, the USF cannot sustain higher friction with increasing the mean density, and therefore rearranges to a non-uniform shear-banded state of lower dynamic friction. For a fixed restitution coefficient, we have found that the gradient of the dynamic friction coefficient with respect to the mean density takes values from negative to positive and crosses zero at the onset mean density $\phi_{onset}^{0}$. The onset mean density has been found to be a monotonically decreasing function of the restitution coefficient satisfying the cubic polynomial relation (53). Consequently, the USF of nearly elastic particles reaches to the onset of the shear-banding instability at a lower mean density as compared to that of relatively more inelastic particles, which is in contrast to Alam et al. (2008) as the constitutive models employed in Alam et al. (2008) are valid only for nearly elastic particles. Furthermore, the onset of the shear banding has been found to follow two global criteria relating (i) the channel width and restitution coefficient and (ii) the mean density and restitution coefficient (or the shear viscosity and pressure).

Within the framework of dense granular flows with arbitrary inelasticity, the present work provides the control parameters for the onset of the shear-banding instability. Functional relationships relating the onset parameters have been obtained that enable us to predict the existence of shear bands in a granular shear flow by merely knowing the control parameters. This could also be of much interest from an experimental point of view as the onset of the shear-banding instability becomes completely known from the present work in terms of control parameters. It is important, however, to note that the findings of the present paper are based on the NSF equations closed with the GDL model that neglects the anisotropy of the USF. To incorporate the effect of anisotropy of the USF, the generalized transport coefficients for dense granular flows are required. The linear stability analysis of the USF with the generalized transport coefficients will be a topic of future research. Notwithstanding, the present work paves the way for simulations and experiments on granular shear flows of arbitrarily inelastic particles.

Acknowledgment

The authors acknowledge the anonymous reviewers for their valuable suggestions, which significantly improved the paper. PS acknowledges financial support from IIT Madras through the grant MAT/16-17/671/NFSC/PRIY. VKG gratefully acknowledges the financial supports through the “MATRICS” project MTR/2017/000693 funded by the SERB, India and that through the Commonwealth Rutherford Fellowship.

Appendix A Coefficients in dispersion relation (34)

The coefficients $a_{ij}$ in (35) that enters dispersion relation (34) are as follows.

[TABLE]