On the shearbanding instability of arbitrarily inelastic granular shear flows

On the shearbanding instability of arbitrarily inelastic granular shear flows


One prototypical instability in granular flows is the shearbanding 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 shearbanding 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 shearbanding 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 onset of the shearbanding instability.


revtex4-1Repair the float

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, shearbanding, etc. Jaeger et al. (1996); Ottino and Khakhar (2000); Aranson and Tsimring (2006); Alam et al. (2009). Among all the non-uniform structures exhibited by granular flows, the shearbanding in granular flows—in particular—has received tremendous attention mainly due to its analogy with the shearbanding in soft matters, e.g. foams, emulsions, colloidal suspensions, etc. Schall and van Hecke (2010). The shearbanding 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 and Richman (1985a); Sela and Goldhirsch (1998); Alam and Luding (2003a); Khain and Meerson (2006); Shukla and 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 and Glasser (2004); Liao et al. (2010); Murdoch et al. (2013)) as well as theoretical studies (see e.g. Alam and Nott (1998); Conway and Glasser (2004); Conway et al. (2006); Khain and Meerson (2006); Conway et al. (2006); Khain (2007); Shukla and Alam (2009); Khain (2009); Shukla and Alam (2011b, a); Berzi and Jenkins (2015)) have confirmed the shearbanding phenomenon in granular shear flows. The theoretical studies based on linear and nonlinear stability analyses Alam and Nott (1998); Shukla and Alam (2009, 2011b, 2011a) show that a granular shear flow admits stationary and traveling instabilities leading to clustering and shearbanding—with the latter being the focus of the present work. In the phenomenon of shearbanding, 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 shearbanding instabilities have also been observed in experiments Mandl et al. (1977); Veje et al. (1999); Conway and 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 and Briels (2008), which provide comprehensive details of the shearbanding phenomenon.

To investigate such a vast varieties 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. Similar 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ó and Dufty (1999); Sela and Goldhirsch (1998); Brilliantov and Pöschel (2004); Gupta (2011). Alternatively, these constitutive relations can also be derived from the moment equations, see e.g. Jenkins and Richman (1985b, a); Kremer and Marques Jr. (2011); Garzó (2013); Gupta and 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 and 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 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 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 () 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 shearbanding 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 kept the other transport coefficients same as those in the original JR model, and showed that some of their models underpin the shearbanding while the others do not. Moreover, they also established that the onset of shearbanding instabilities could be predicted with the knowledge of the pressure and shear viscosity of the system. This result indicates that the emergence of the shearbanding 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 shearbanding, Khain and Meerson (2006) substantiated that the experimentally observed shearbanding instabilities could not be perceived with the usual constitutive relations for dense granular gases; however, shearbands, 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 shearbanding instability predictions of both Alam et al. (2008) and Khain and 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ó and 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ó and Dufty (1999) were generalized to arbitrary number of dimensions by Lutsko (2005). We shall refer to the constitutive relations obtained in Ref. Garzó and 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 in the latter while only the collisional dissipation depends on 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 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 shearbanding instability 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 shearbanding instability in arbitrarily inelastic dense granular flows has never been addressed explicitly using the granular hydrodynamics. In this context, the goal of the present work is to analyze the shearbanding instability in arbitrarily inelastic dense granular flows of hard disks by using the GDL constitutive relations. By means of the linear stability analysis of a uniform shear flow (USF), the onset of the shearbanding instability is predicted. In contrast to previous studies Khain and 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.

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 can be described by the mass, momentum and energy balance equations, which read Brilliantov and Pöschel (2004); Garzó (2019)


Here, is the mass density with being the material density and being the volume fraction of grains; is the coarse-grained velocity with , and being its components in the -, - and -directions, respectively; is the granular temperature; is the pressure tensor, is the (granular) heat flux; is the collisional dissipation due to inelastic collisions among grains; and 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: , and . These unknowns are typically expressed in terms of the hydrodynamic variables , , and their spatial gradients by means of the Chapman–Enskog expansion, see e.g. Goldshtein and Shapiro (1995); Brey et al. (1998); Sela and Goldhirsch (1998); Garzó and Dufty (1999); Brilliantov and Pöschel (2004); Gupta (2011). To first order in spatial gradients of the hydrodynamic variables, these unknowns are expressed as Goldshtein and Shapiro (1995); Sela and Goldhirsch (1998)


where ; is the identity tensor; and is the deviatoric strain rate tensor; the quantities , , , , , and 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:


Here all the ’s are the dimensionless functions of the volume fraction and restitution coefficient only. It is worthwhile to note that the values of ’s are different for different models for the constitutive relations. From the GDL model Garzó and Dufty (1999); Lutsko (2005), ’s in the case of hard-disk flows read


where with being the pair correlation function and


The pair correlation function is adopted from the work of Torquato (1995)


which approaches to the random close packing fraction while achieving the freezing point . Following Garzó and Dufty (1999) and Lutsko (2005), a transport coefficient can be decomposed in its kinetic and collisional shares:

where is the representation of the transport coefficients (5), and the superscripts and denote its kinetic and collisional shares, respectively. In order to do so, the dimensionless functions, , are decomposed as , where . Owing to this decomposition, two limiting cases, namely the kinetic or dilute limit ( and ) and the collisional or dense limit, of the GDL model (6)–(7) can be distinguished. In the kinetic or dilute limit, the quadratic terms of density can be neglected whereas in the collisional or dense limit, the terms of the linear terms of density can be neglected. The two limiting cases of the GDL model are as follows.

Kinetic limit


Collisional limit


Iii Non-dimensionalization and base state

We shall investigate a granular plane shear flow of hard disks driven by two oppositely moving parallel plates at , with a speed of each plate along the streamwise () direction; the plates are separated by a height or gap along the -direction and hence the overall shear rate is , see figure 1.

Figure 1: Schematic of the plane Couette flow of granular hard disks.

The velocity and temperature boundary conditions for the problem are taken as the no-slip and zero heat flux on the plates, respectively. In other words, the boundary conditions read


For the purpose of non-dimensionalization, we choose as a reference length, as a reference velocity, the inverse of the shear rate as a reference time and 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,


where is referred to as the dimensionless height between the parallel plates or the Couette gap. With as the scaling for the density, the volume fraction also denotes dimensionless density. With these scales, governing equations (1) along with the constitutive relations (5) in the dimensionless form read




Boundary conditions (11) in the dimensionless form read


We shall be investigating the shearbanding 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. -) direction. Such instability arises from the perturbations having no variation in the streamwise (i.e. -) direction. The streamwise independent (i.e. ) governing equations (13) in the dimensionless form read


iii.1 Base state: Uniform shear flow

The basic flow, whose stability is to be analyzed here, is assumed to be a steady , fully developed , plane shear flow of the following form

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

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


where and are the integration constants. The above set of equations with no-slip and zero heat flux boundary conditions (15) admits the following base state solution: where the superscript ‘0’ denotes the base state solutions, and are constants, and and . 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. , thus leading to the USF. The base state density , the Couette gap and the restitution coefficient are the control parameters of the problem.

Iv Linear stability analysis

iv.1 Perturbation equations

For the linear stability analysis of the USF, the field variables are decomposed into the base state solution plus the perturbation from the base state solution as follows:


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 (18) into governing equations (16) and retaining only the linear terms of the perturbed field variables, one obtains the governing equations for the perturbed field variables:


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


where is the vector of perturbed fields and is the matrix of linear differential operators given by


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


iv.2 Analytical solutions

We perform the standard linear stability analysis on the USF by assuming a normal mode solution of the form , where and is the complex frequency whose real part represents growth or decay rate of the perturbations and imaginary part denotes the oscillation of the perturbation. Substituting this normal mode solution into the linearized perturbation equations (20)–(22), we get the following matrix eigenvalue problem:


where . It has been verified that eigenvalue problem (23) has an analytical solution in terms of sine and cosine functions Alam and Nott (1998), as


where are the mode numbers and are the constant amplitude of the normal mode solution. For instance, is the fundamental mode and is the second harmonic of the normal mode solution, etc. With solution (24), problem (23) simplifies to


where and


For the nontrivial solutions of (25) . This condition is the dispersion relation and an be written as




Here the coefficients 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 (27) is a fourth degree polynomial in with real coefficients and, therefore, there are three possibilities for four roots of (27): (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 with , one can find the analytical expressions of the eigenvalues (the roots of (27)), as follows. Let the frequency is represented by


where are unknown coefficients. Substituting this ansatz along with (28) into (27), and comparing each power of on both sides of the resulting equations, one obtains algebraic equations, which are solved for the unknowns ’s in (29). Exploiting these values, one obtains four roots form (27), which are given by


where the subscripts and represent the real and imaginary parts, respectively, of the roots and


It is evident from (30)–(32) that in the limit of large , dispersion relation (27) has two real roots and a complex conjugate pair of roots . In the limit of large , it is verified numerically that and are always negative resulting into the least stable mode as . Figure 2 illustrates the four eigenvalues for large obtained through the asymptotic analysis of dispersion relation (27) (solid line) and those obtained by solving (25) numerically (symbols) for , and . It can be seen from the figure that the eigenvalues from both the methods are in excellent agreement for large .

Figure 2: The eigenvalues of (25) for large for parameters , and . The solid lines denote the eigenvalues from (30)–(32), which were obtained through the asymptotic analysis of (27), and the symbols delineate those obtained by solving matrix eigenvalue problem (25) numerically.

V Results

The linear stability of the USF of granular material has been studied previously for the case of nearly elastic particles Alam and 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ó and Dufty (1999); Lutsko (2005). In order to perform the linear stability analysis, eigenvalue problem (23) is solved numerically using the Chebyshev spectral collocation method. In addition, eigenvalue problem (25) analytically, which is obtained using the exact solutions (24), is also solved. It is verified that the eigenvalues obtained numerically using the Chebyshev spectral collocation method and analytically using the exact solutions (24) are found in an excellent agreement (figure not shown for brevity). Therefore in the present analysis, the eigenvalues are computed by solving (25). All the computations have been performed in MATLAB. For analyzing the stability results, we define the least stable eigenvalue as one of the eigenvalues whose real part is maximum for a fixed mode number, and the dominant eigenvalue is one of the least stable eigenvalues whose real part is maximum over all the mode numbers , i.e.


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, , are shown in the -plane for the inelastic particles with the restitution coefficient , see figure 3(a). The real part of the dominant mode is positive (i.e ) inside the zero contour and negative (i.e ) 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 is always zero. This can also be seen from the asymptotic analysis presented in  Sec. IV.3 that the least stable eigenvalue is always real. Therefore the shearbanding instabilities in granular shear flow are due to the stationary waves.

Figure 3: (a) Contours of the positive, negative and zero growth rates in the -plane, and (b) the variation of the growth rate with the Couette gap for . The restitution coefficient is set to .

We have also examined that the lowest mode number, i.e. , 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 Couette gap and a critical mean density above which the USF is unstable and below which it is stable. That is to say, the USF becomes unstable if either or . In addition, it can also be seen from figure 3(a) that the USF remains stable for all densities and for all Couette gaps below an onset value of the mean density, say (see Sec. V.1 for more details). Note that, the onset mean density depends only on the restitution coefficient. In particular, for the flow is stable when for all values of .

In order to get more insight, the variation of the dominant growth rate with the Couette gap 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 to , where is a positive integer. It is seen that the dominant mode number increases with the Couette gap, i.e.  mode is the first one to loose stability and remains unstable with increasing Couette gap until it crosses mode (first kink in figure 3(b)), thereafter mode becomes dominant mode until it crosses the next mode (second kink in figure 3(b)), and so on. For parameter values shown in figure 3(b), the USF becomes unstable for .

Figure 4: Eigenfunctions (24) for three values of the Couette gap (first row), (second row) and (third row) corresponding to , and , respectively. The other parameter values are the same as those in figure 3(b).

Figure 4 illustrates the density, temperature and velocity eigenfunctions for three values of the Couette gap, , , , with other parameters being the same as in figure 3(b). The eigenfunctions displayed in the first row correspond to mode for which that vanish once at in the flow domain, therefore the density and temperature solutions attain their mean values i.e.  once in the flow domain. The second row of figure 4 corresponds to mode which gives eigenfunctions that vanish twice along the flow domain at . Similarly, the third row corresponds to mode having the eigenfunctions that vanish thrice at . Consequently, the density and temperature solutions attain its mean values once, twice and thrice in the flow domain for , and , respectively. By analyzing the corresponding velocity components (second column in figure 4), it is evident that while the horizontal component of the velocity eigenfunction varies significantly, its transverse component varies only slightly.

Figure 5 shows the neutral stability curve for the case when only the kinetic contribution to the transport coefficients (9) (i.e. for the dilute limit) is considered. It is evident from the figure that the USF is always stable in the dilute limit (), which is also known from the previous studies Alam et al. (2008); Conway and Glasser (2004); Conway et al. (2006). Therefore the dilute model (the dilute limit of the GDL model) (9) is able to capture the stability of USF in the dilute limit ().

Figure 5: Neutral stability curve in the -plane using dilute constitutive relations (9) of the GDL model. The restitution coefficient is .

While the dilute model gives correct prediction, the dense model (the dense limit of the GDL model) (10), which is expected to be valid in collisional limit, fails to capture the shearbanding instabilities appearing in dense granular shear flows Alam and Luding (2003b); Conway and Glasser (2004); Conway et al. (2006) correctly (phase diagram is not shown for brevity). The dense model predicts that the USF is stable in the dense limit whereas the full GDL model (see figure 3) and MD simulations Conway and Glasser (2004); Conway et al. (2006) predict the shearbanding instability for the dense USF. In addition, the kinetic contribution is small in the dense limit and therefore both models (the full GDL model and its dense limit) are expected to support shearbanding however this is not true. By analyzing three variants of the GDL model viz. the full model (6), the kinetic limit (9) and collisional limit (10), one can conclude that the choice of constitutive relations plays an important role in determining the shearbanding instability. For the correct instability predictions, both kinetic as well as collisional mechanisms are significant, and therefore none of them can be discarded throughout the flow regimes. In the following subsections, we focus on the effect of the inelasticity on the shearbanding instability through the full GDL model.

v.1 Effect of the restitution coefficient and Couette gap on the shearbanding: a global criterion

Figure 6: (a) Neutral stability curves in the -plane for increasing restitution coefficient ; the flow is stable outside each neutral stability curve and unstable inside it. (b) Variation of with where the flow is unstable for for each .

Figure 6(a) illustrates the neutral stability curves in the -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. 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 Couette gap for which the USF unstable decreases with increasing restitution coefficient. Indeed, there exists a set of critical parameters () 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 6(a) asymptotically (as ) 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.


where the USF is stable for and vice versa. The onset mean density (36) as a function of the restitution coefficient is plotted in figure 6(b). It is clear that the onset density decreases with increasing the restitution coefficient, which implies that the USF tends to lose at a lower density in the elastic limit than that in the inelastic case. Note that the above conclusion holds for very large Couette gaps.

Figure 7: (a) Neutral stability curves in the -plane for increasing Couette gap , , , and (b) variation of with (main panel), and with (inset), where the solid line and circles represent the values from relation (39) and those extracted from the neutral stability curves in figure 6, respectively.

It can also be noticed from figure 6(a) that there exists an onset value of the Couette gap—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 Couette gap as the minimum of all critical Couette gaps:


Note that and depend on both the mean density and restitution coefficient whereas and are only function of the restitution coefficient and hence of the inelasticity. It is evident from figure 6(a) that the value of increases with increasing the restitution coefficient implying that the USF is more stable as increases which qualitatively matches with the simulation results of Conway and Glasser (2004); Conway et al. (2006). Owing to the functional dependence (see (37)), the shearbanding in granular shear flow starts to appear at small Couette gaps 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 Couette gaps.

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


It can also be seen from the inset of figure 7(b) that increases with increasing the Couette gap.

Although the neutral stability curves shown in figure 6 and 7 depict the overall behavior of the instability, it still remains to understand how the onset parameters and vary and to discern if there exists any relation relating these parameters. We shall now seek the onset of the shearbanding instability in terms of these onset parameters. As mentioned earlier, the USF is unstable for all