2.1. Physical set-up, control parameters and assumptions for the theoretical analysis

Figure 1 shows the problem under consideration: a quasi-two-dimensional box of length $L_{x}$ and depth $D$ (of a few particle diameters) with open top (the height of the box is large enough that no particles will escape under shaking) which is partially filled with spheres of diameter $d$ having a material density of ${\it\rho}_{p}$ and mass $m={\it\rho}_{p}{\rm\pi}d^{3}/6$ . This box of particles is vibrated vertically with a sine wave:

where $a_{s}$ and $f_{s}$ are the amplitude and frequency of shaking, respectively. The energy injected per particle over one shaking period ( ${\it\tau}_{s}=1/f_{s}$ ) is

where $v_{s}=\text{d}y_{s}/\text{d}t=2{\rm\pi}f_{s}a_{s}\cos (2{\rm\pi}f_{s}t)$ is the translation velocity of the shaker.

Figure 1. Schematic diagram of a harmonically shaken box of particles. This is an ‘open-top’ quasi-2D box of horizontal length $L_{x}$ and a depth $D$ of a few particle diameters.

Four dimensionless control parameters are needed to describe this system: (i) the particle loading parameter or the number of particle layers at rest,

where $N_{p}=(L_{x}/d)\times (D/d)\times (h_{0}/d)$ is the total number of particles, with $h_{0}=Fd$ being the height of the granular bed at rest; (ii) the shaking strength,

which is a ratio of the kinetic energy injected to the system ( $T_{p}N_{p}$ ) and the potential energy associated with all particles ( $mgdN_{p}$ ); (iii) the normal restitution coefficient $e$ and (iv) the length of the box $L_{x}/d$ . The shaking strength (2.4) can also be re-expressed as

where

is the dimensionless shaking acceleration, called the shaking intensity, and $g$ is the acceleration due to gravity.

For the stability analysis below, the quasi-2D box in figure 1 is assumed to be two-dimensional (i.e. to have no variations along the depth) with periodic boundary conditions along the $x$ -direction. We further assume that the granular Leidenfrost state (see § 2.4) can be characterized by a ‘quasi-steady’ state, wherein we specify a constant temperature (2.2) at the base plate which is calculated over one shaking period. (The temperature perturbation being zero at the base follows from the above mean-field approximation, see § 3.1.) A detailed energy balance at the base can yield accurate boundary conditions (Hui et al. Reference Hui, Haff, Ungar and Jackson1984; Jenkins & Richman Reference Jenkins and Richman1986; Rao & Nott Reference Rao and Nott2008) at the expense of making the present nonlinear analysis complicated. As it is the first nonlinear analysis of granular convection phenomena, the present work is rather focused on identifying the nonlinear equilibrium solutions by employing the hydrodynamic model (§ 2.2 and § 2.3) and boundary conditions with identical assumptions to those in the previous linear analysis of the same problem (Eshuis et al. Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010, Reference Eshuis, van der Weele, Alam, van Gerner, van der Höf, Kuipers, Luding, van der Meer and Lohse2013). The effect of different forms of base boundary conditions on the present results can be taken up in a future work.

2.2. Balance equations and constitutive relations

We use a Navier–Stokes-order hydrodynamic model for which the balance equations for mass, momentum and granular energy are (Jenkins & Richman Reference Jenkins and Richman1985a ; Sela & Goldhirsch Reference Sela and Goldhirsch1998):

where $n(\boldsymbol{x},t)$ is the number density of particles, $\boldsymbol{u}(\boldsymbol{x},t)=\langle \boldsymbol{c}\rangle$ is the coarse-grained velocity and $T(\boldsymbol{x},t)=\langle mC^{2}/2\rangle$ is the granular temperature, which is defined as the mean-square fluctuation energy, with $\boldsymbol{C}=(\boldsymbol{c}-\boldsymbol{u})$ being the peculiar velocity of particles and $\boldsymbol{ c}$ the instantaneous particle velocity; $\boldsymbol{g}=(0,-g)$ is the gravitational acceleration and $m$ is the mass of a particle. The flux terms are the stress tensor, ${\it\bf\Sigma}$ , and the granular heat flux, $\boldsymbol{q}$ . In (2.9), ${\it\bf\Sigma}:\boldsymbol{{\rm\nabla}}\boldsymbol{u}$ represents the production of fluctuation energy due to the coupling between ${\it\bf\Sigma}$ and the velocity gradient $\boldsymbol{{\rm\nabla}}\boldsymbol{u}$ , and $\mathscr{D}$ is the rate of dissipation of granular energy per unit volume. To close the above balance equations, the stress tensor is taken to be of Newtonian form and the heat flux follows the Fourier law: where $\unicode[STIX]{x1D644}$ is the identity tensor and

is the deviatoric part of the deformation rate tensor. Here, $p$ , ${\it\mu}$ , ${\it\zeta}$ and ${\it\kappa}$ are the pressure, shear viscosity, bulk viscosity and thermal conductivity, respectively. In many recent works, instead of the standard Fourier law (2.11), a generalized Fourier law is adopted in which a density-gradient-dependent term ( $=-{\it\kappa}_{h}\boldsymbol{{\rm\nabla}}n$ , where ${\it\kappa}_{h}$ is the Dufour coefficient) is considered (Sela & Goldhirsch Reference Sela and Goldhirsch1998; Brilliantov & Pöschel Reference Brilliantov and Pöschel2004), although a recent work (Almazan et al. Reference Almazan, Carrillo, Saluena, Garzo and Pöschel2013) on the dynamics of a vibrated bed has confirmed that the incorporation of this term can lead to a temperature maximum far away from the bottom plate in contradiction to experimental and simulation results. It may be noted that this Dufour term is of order $O({\it\epsilon}^{3})$ (where ${\it\epsilon}=\sqrt{(1-e^{2})}$ is the degree of inelasticity), as demonstrated by Sela & Goldhirsch (Reference Sela and Goldhirsch1998), and hence does not belong to the Navier–Stokes order $O({\it\epsilon}^{2})$ .

As long as we restrict our consideration to the Navier–Stokes-order description, the balance equations (2.7)–(2.9) and flux relations (2.10)–(2.11) are general and no assumptions have been made up to this point. Next, we need constitutive relations for $p$ , ${\it\mu}$ , ${\it\zeta}$ , ${\it\kappa}$ and $\mathscr{D}$ for which different models are available (Goldhirsch Reference Goldhirsch2003). In this paper, we have chosen the following constitutive relations used by Eshuis et al. (Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010, Reference Eshuis, van der Weele, Alam, van Gerner, van der Höf, Kuipers, Luding, van der Meer and Lohse2013):

where $n_{c}=2/\sqrt{3}d^{2}$ is the number density of a triangular-packed crystal, $e$ is the coefficient of restitution and $l$ is the mean free path (Grossman, Zhou & Ben-Naim Reference Grossman, Zhou and Ben-Naim1997):

with ${\it\alpha}_{1}=1-\sqrt{3/8}$ . It should be noted that the expression of shear viscosity is taken to be proportional to the thermal conductivity, with the coefficient of proportionality being dubbed the ‘effective’ Prandtl number ( $\mathit{Pr}$ ); the numerical value of $Pr=1.7$ was used by Eshuis et al. (Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010), which yielded good agreement between the linear stability predictions and experiments. For a dense gas, the Prandtl number ( $\mathit{Pr}$ ) is no longer a constant but should be a function of the density; the effect of such a density-dependent Prandtl number is probed in § 6.1. We have set the bulk viscosity to zero ( ${\it\zeta}=0$ ), and its impact on stability results is assessed in § 6.1. Of course, the constitutive model (2.13) is based on heuristic arguments (Haff Reference Haff1983; Grossman et al. Reference Grossman, Zhou and Ben-Naim1997) but is expected to be valid for nearly elastic particles ( $e\sim 1$ ) for a wide range of density, as we discuss below.

Grossman et al. (Reference Grossman, Zhou and Ben-Naim1997) derived the expressions for $p$ , ${\it\kappa}$ and $\mathscr{D}$ by employing free-volume arguments in the dense limit and subsequently using interpolations with the respective constitutive relations in the dilute limit. These interpolations include two order-one fitting parameters ${\it\gamma}_{c}$ and ${\it\alpha}$ that were obtained by matching the theoretical profiles for density and granular temperature in a thermally driven granular gas between two walls (with zero gravity) with those obtained from molecular dynamics simulations. The predicted density profile had a large density variation ranging from a dilute granular gas near one wall to a high-density cluster at the other end. Later on, the same constitutive relations were used by Meerson et al. (Reference Meerson, Pöschel and Bromberg2003) and Eshuis et al. (Reference Eshuis, van der Weele, van der Meer and Lohse2005) to predict the granular Leidenfrost state, which also consists of an extreme density contrast, with a dense cluster floating over a dilute gas under vertical vibration with gravity. Moreover, the linear stability analysis of the granular Leidenfrost state with (2.13) yielded quantitative agreement for the onset of convection instability with both experiments and simulations with just one fitting parameter ( $Pr$ ), and hence the constitutive relations (2.13) are likely to be valid for a wide range of density spanning from the dilute to the dense regimes. This motivated us to continue with the same constitutive model to analyse the nonlinear stability of the granular Leidenfrost state.

Given that the simplified constitutive model (2.13) works well for linear stability predictions, the goal of the present paper is to ascertain the nonlinear saturation of the same linear instability modes when one takes into account finite-amplitude perturbations. To establish the sensitivity of our nonlinear stability predictions (§ 5) on the adopted constitutive relations, we will assess the effects of (i) the bulk viscosity ( ${\it\zeta}\neq 0$ ), (ii) the viscous dissipation or the shear work and (iii) a density-dependent Prandtl number in § 6.1 of this paper. Moreover, the different boundary conditions at the free surface and their possible effects are discussed in § 6.2. A detailed parametric study with a different set of constitutive relations (Jenkins & Richman Reference Jenkins and Richman1985a ; Sela & Goldhirsch Reference Sela and Goldhirsch1998; Garzo, Santos & Montanero Reference Garzo, Santos and Montanero2007) is, however, beyond the scope of the present work. In this regard, we would like to remind the readers that the different forms of constitutive relations (at Navier–Stokes order) do not qualitatively affect the onset of instabilities, as demonstrated previously for various types of granular flow: (i) plane shear flow (Alam & Nott Reference Alam and Nott1998; Alam et al. Reference Alam, Arakeri, Nott, Goddard and Herrmann2005; Gayen & Alam Reference Gayen and Alam2006; Alam, Shukla & Luding Reference Alam, Shukla and Luding2008; Alam & Shukla Reference Alam and Shukla2012), (ii) Poiseuille flow (Alam, Chikkadi & Gupta Reference Alam, Chikkadi and Gupta2009) and (iii) inclined Chute flow (Forterre & Pouliquen Reference Forterre and Pouliquen2002; Woodhouse & Hogg Reference Woodhouse and Hogg2010). Going beyond the Navier–Stokes order, it might be worthwhile to carry out a comparative study for both linear and nonlinear stability of granular convection using a non-Newtonian constitutive model (Saha & Alam Reference Saha and Alam2014).

2.3. Dimensionless forms of the balance equations and constitutive relations

For non-dimensionalization, we use the maximum number density ( $n_{R}=n_{c}$ ) as the reference scale for $n(\boldsymbol{x},t)$ , the granular temperature at the bottom plate as the reference temperature ( $T_{R}=T_{p}$ ), the particle diameter as the reference length ( $L_{R}=d$ ), $U_{R}=\sqrt{T_{p}/m}$ as the reference velocity and $t_{R}=L_{R}/U_{R}$ as the reference time:

and from here onwards all equations and symbols are written in dimensionless forms. The dimensionless balance equations for the number density, $x$ -momentum, $y$ -momentum and granular energy, respectively, are

The dimensionless forms of the transport coefficients are $p(n,T)=f_{1}(n)T$ , ${\it\mu}(n,T)=f_{2}(n)T^{1/2}$ , ${\it\lambda}(n,T)={\it\zeta}(n,T)-{\it\mu}(n,T)\equiv -{\it\mu}(n,T)$ , ${\it\kappa}(n,T)=f_{4}(n)T^{1/2}$ and $\mathscr{D}(n,T)=f_{5}(n)T^{3/2}$ , where $f_{i}(n)$ are dimensionless functions of the number density defined as

and $l(n)={\it\alpha}_{0}(n^{-1}-1)/(1-{\it\alpha}_{1}n)$ is the dimensionless mean free path, with ${\it\alpha}_{0}=\sqrt{3/32}$ and ${\it\alpha}_{1}=1-\sqrt{3/8}$ . In (2.18), $S$ is the scaled granular temperature at the bottom plate,

which is the same as the shaking strength as defined in (2.5). The last term on the right-hand side of (2.19),

represents viscous dissipation (or shear work), which was set to zero in the previous linear stability analysis (Eshuis et al. Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010, Reference Eshuis, van der Weele, Alam, van Gerner, van der Höf, Kuipers, Luding, van der Meer and Lohse2013) but is retained in the present nonlinear analysis. The importance of this term on nonlinear stability will be assessed in § 6.

2.4. Base state: granular Leidenfrost state

The quiescent base state is assumed to be steady and one-dimensional with no variation along the horizontal ( $x$ ) direction:

On substituting (2.23) into the hydrodynamic equations, (2.16)–(2.19), we obtain the following equations for the base state:

To complete the above system we need to specify three boundary conditions (Eshuis et al. Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010). The (dimensionless) granular temperature at the bottom wall ( $y=0$ ) is unity, and the heat flux, ${\it\kappa}\boldsymbol{{\rm\nabla}}T$ , is zero far away from the bottom wall ( $y\rightarrow \infty$ ). The third condition is associated with particle conservation, which can be written as an integral constraint

where $F=h_{0}/d$ is the number of particle layers at rest (2.3) and

is the effective height of the expanded fluidized bed (Jackson Reference Jackson2000; Rao & Nott Reference Rao and Nott2008). It should be noted that the bed height $H$ is not a constant, rather it depends on the shaking strength $S$ for specified values of $F=h_{0}/d=H_{0}$ (i.e. the bed height at rest) and the restitution coefficient ( $e$ ). Since $S$ is a measure of the injected energy, a larger value of $S$ would imply a larger expansion of the granular bed and hence a larger value of $H$ , satisfying the constraint (2.26) of particle number-density conservation. The base-state equations (2.24), (2.25) with boundary conditions are solved numerically using the Runge–Kutta method. For all numerical calculations, we varied the bed height ( $H$ , (2.27)), depending on $S$ , so as to compute the base state accurately such that the constraint (2.26) is satisfied within an error of $O(10^{-6})$ .

Figure 2. Profiles of the density (blue solid line) and temperature (red dashed line) of the Leidenfrost base state for $e=0.9$ and $(F,S)=$ (a) (6, 100) and (b) (6, 400). (c) Density profiles in semilog scale for different values of $S$ for an initial bed height of $F=6$ . (d) Effect of the initial bed height $F$ on the density profiles at $S=100$ ; see the text for a description of three regions (I–III) in the density profile for $F=16$ .

Figure 2(a,b) show the base-state density (blue solid line) and granular temperature (red dashed line) profiles for two sets of parameters $(F,S)=(6,100)$ and (6, 400). It is clear that the density profile has a maximum at some vertical location away from the bottom plate, but the corresponding granular temperature decreases monotonically with increasing height. It should be noted that the bed expands (i.e. $H$ , (2.27), increases) with increasing shaking strength $S$ , see figure 2(a,b) for $S=100$ and 400, respectively. A comparison of density profiles in the semilog scale for different values of $S$ in figure 2(c) clarifies the existence of a ‘ballistic’ dilute layer at the top of the bed in which the density decreases exponentially with height. The latter observation follows from the dilute limit of (2.24):

which holds for large $y\gg F$ , where $T_{b}$ is the average granular temperature in the ballistic layer. It should be noted that the density is very low at the top of the bed, with particles undergoing only ballistic motion there, but the corresponding temperature variation is much smaller (see figure 2 a), and can be approximated by its average value ( $T_{b}$ ) at the leading order.

For the same injected energy ( $S=100$ ) as in figure 2(a) but with increasing particle loading ( $F$ ), the effective height of the bed ( $H$ ) decreases, as is evident from a comparison of the density profiles in figure 2(d) for $F=6$ , 10 and 16. In all cases the density at the top of the bed decreases exponentially with height, but the thickness of this ‘ballistic’ layer decreases sharply with increasing $F$ . It should be noted that the maximum density in the ‘levitated’ dense layer and the density within the collisional dilute layer (beneath the dense layer) increase with increasing $F$ , leading to a more compact bed at the same energy input. Referring to the case of $F=16$ (the red curve in figure 2 d), we identify three distinct regions in the density profile marked by I, II and III, respectively: (I) a ‘collisional’ dilute layer near the vibrating plate, (II) a ‘levitated’ dense layer at some intermediate height and (III) a ‘ballistic’ dilute layer at the top of the bed. This is the so-called granular Leidenfrost state (Eshuis et al. Reference Eshuis, van der Weele, van der Meer and Lohse2005), or the floating cluster (Meerson et al. Reference Meerson, Pöschel and Bromberg2003). It is clear from figures 2(c,d) that the extents of the above three regions depend on both the particle loading ( $F$ ) and the shaking strength ( $S$ ). As stated before, the present work is focused on the nonlinear stability of these density-inverted base states, with the goal of ascertaining the resulting nonlinear patterns of granular convection.

3.1. Linear and adjoint eigenvalue problems

To analyse the stability of the Leidenfrost state, we superimpose a finite-amplitude perturbation on the static base state,

where the superscripts zero and prime denote the base-state and disturbance quantities, respectively. By substituting (3.1) into (2.16)–(2.19) and subtracting the base flow (2.24), (2.25), we obtain the nonlinear evolution equations for the disturbance vector $\boldsymbol{X}=(n^{\prime },u^{\prime },v^{\prime },T^{\prime })^{\text{T}}$ :

where $\mathscr{L}=[l_{ij}]$ is the linear operator (appendix A) and $\mathscr{N}_{2}$ and $\mathscr{N}_{3}$ are the quadratic and cubic nonlinear terms (appendix B), respectively.

The boundary conditions are taken as

where the boundary operators $\mathscr{B}_{0}$ and $\mathscr{B}_{H}$ are linear as defined by

It should be noted that we have approximated the semi-infinite domain $y\in (0,\infty )$ by a truncated finite domain $y\in (0,H)$ as in the calculation of our base state. It should be recalled from § 2.4 that the domain height $H$ (2.27) is chosen so as to satisfy the particle conservation constraint (2.26) to an accuracy of $O(10^{-6})$ . We would like to emphasize here that such a representation of an infinite domain by a truncated domain is quite common in a variety of free shear layers (Lakkaraju & Alam Reference Lakkaraju and Alam2007) as well as in granular flows (Forterre & Pouliquen Reference Forterre and Pouliquen2002) – this has advantages from the viewpoint of numerical implementation too (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988), see § 4.2. In (3.4), the temperature perturbation being zero at $y=0$ follows from our mean-field approximation of a constant temperature at the base, as discussed in § 2.1. It should be noted that these boundary conditions (3.4) were also used by Eshuis et al. (Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010, Reference Eshuis, van der Weele, Alam, van Gerner, van der Höf, Kuipers, Luding, van der Meer and Lohse2013) for the linear stability analysis of the same problem. We shall show in § 6.2 that changing the free-surface boundary conditions on (i) the perturbation temperature (e.g. a zero heat-flux condition, $\partial T^{\prime }/\partial y=0$ , at $y=H$ ) and (ii) the perturbation velocities does not noticeably affect our predictions on the nonlinear stability of the present problem.

The linear stability problem is obtained by neglecting the nonlinear terms in the disturbance equations (3.2):

Since the linear operator and the boundary conditions are translation invariant in $x$ and time, we can seek a normal-mode solution

where $\widehat{\boldsymbol{X}}$ and $c$ denote the complex eigenfunction and eigenvalue, respectively, and $k_{x}$ is the real streamwise wavenumber. By substituting (3.6) into (3.5) and (3.3) we obtain an eigenvalue problem

where $\unicode[STIX]{x1D647}=\mathscr{L}(\partial /\partial x\rightarrow \text{i}k_{x},\partial /\partial y\rightarrow \text{d}/\text{d}y,\dots )$ is the complex-valued linear differential operator, and the boundary operators transform as $B_{0}=\mathscr{B}_{0}(\partial /\partial y\rightarrow \text{d}/\text{d}y)$ and $B_{H}=\mathscr{B}_{H}$ . For the temporal development of perturbations (3.6), the flow is said to be stable, unstable or neutrally stable if the real part of any eigenvalue is negative, positive or zero ( $\text{Re}(c)<0,>0,=0$ ), respectively.

Since the linear stability operator is not self-adjoint, we need to formulate the adjoint eigenvalue problem whose solution is needed to compute the Landau coefficients (§ 3.2). The adjoint operator, $\unicode[STIX]{x1D647}^{\dagger }$ , for the linear stability problem is defined via

where the following definition of the inner product is adopted:

Here, $\boldsymbol{X}_{1}(y)$ and $\boldsymbol{X}_{2}(y)$ are two complex-valued vector functions, and the ‘tilde’ on a quantity denotes its complex conjugate. The explicit forms of $\unicode[STIX]{x1D647}^{\dagger }$ and the corresponding boundary conditions are determined from (3.8) via integration by parts. The adjoint eigenvalue problem can be written as

where $\widehat{\boldsymbol{X}}^{\dagger }=(\widehat{n}^{\dagger },\widehat{u}^{\dagger },\widehat{v}^{\dagger },\widehat{T}^{\dagger })$ is the adjoint eigenfunction, with associated boundary conditions

The elements of $\unicode[STIX]{x1D647}^{\dagger }$ are given in appendix A.

3.2. Amplitude expansion method and Landau equation

The normal-mode ansatz of linear stability (3.6) can be rewritten as

where $\text{e}^{at}$ accounts for the temporal variation of an infinitesimal disturbance. However, due to the finite size of the disturbances, a single mode interacts with itself and produces its harmonics which distort the base flow. Consequently, the perturbation amplitude can no longer be defined by the exponential term $\text{e}^{at}$ and the frequency of perturbation would also depend on its amplitude. Thus, the disturbance can be expressed as

where the second term is the complex conjugate (denoted by tilde) of the first term. Here, $A$ is the real amplitude of the perturbation, ${\it\theta}$ is the cycle of the fundamental wave and the summation is taken over all positive integers $k\geqslant 0$ . We require that the nonlinear problem at $O(A)$ must reduce to the linear stability problem with an exponential growth of amplitude, i.e. $A\rightarrow \text{e}^{at}$ for infinitesimal perturbations. Further, at $O(A^{2})$ we have only two choices: (i) the interaction of two linear modes giving the second harmonic and (ii) the interaction of the fundamental with its conjugate producing a distortion of the base flow. The above considerations suggest that the disturbance amplitude in (3.13) can be written in the form of a power series in amplitude,

where the superscripts $k$ and $m$ on $\boldsymbol{X}^{[k;m]}(y)$ are related to the Fourier mode and the order of the perturbation amplitude, respectively (Shukla & Alam Reference Shukla and Alam2011b ; Alam & Shukla Reference Alam and Shukla2013). For example, $[k;m]=[1;1]$ refers to the linear/fundamental mode which is of order $O(A)$ , $[2;2]$ is the second harmonic which is of order $O(A^{2})$ , $[0;2]$ is the distortion of the base state at order $O(A^{2})$ , $[1;3]$ is the distortion of the fundamental mode at order $O(A^{3})$ , and so on.

To deal with the time-dependent terms of the disturbance equations we assume that the amplitude of the perturbation evolves according to the Landau equation (Shukla & Alam Reference Shukla and Alam2011b ; Alam & Shukla Reference Alam and Shukla2013),

and the associated frequency ${\it\omega}={\it\omega}(A)$ evolves as

In (3.15), $a^{(0)}$ and $a^{(2)}$ are the linear growth rate and the real part of the first Landau coefficient, respectively. It should be noted that in the limiting case of ‘infinitesimal-amplitude’ perturbations ( $A\rightarrow 0$ ) the nonlinear terms in (3.15) and (3.16) can be neglected, resulting in

which corresponds to the well-known ansatz of exponential growth in the linear stability analysis. We note in passing that the general validity of the Landau equation (3.15), (3.16) to predict nonlinear instabilities has been demonstrated for granular shear flow by deriving the same equation using an altogether different method (centre manifold reduction; Shukla & Alam Reference Shukla and Alam2009).

By substituting (3.14)–(3.16) into the disturbance equations and separating the powers of the amplitude, we obtain the following system of differential equations:

where $\unicode[STIX]{x1D647}_{km}$ and $\unicode[STIX]{x1D647}_{k}$ are the linear operators and $\boldsymbol{G}_{km}$ is a vector representing the sum of linear and nonlinear terms, as given by

with $\boldsymbol{E}_{km}$ and $\boldsymbol{F}_{km}$ being the vectors representing the quadratic and higher-order nonlinear terms, respectively. In (3.18),

are the Landau coefficients, ${\it\delta}_{k1}$ is the Kronecker delta and $\boldsymbol{X}^{[1;1]}$ is the linear eigenfunction.

In the weakly nonlinear regime of small growth rates, $a^{(0)}\sim 0$ , the Landau coefficients (3.21) are determined from the solvability condition of (3.18) at odd orders ( $m=3,5,\dots$ ):

where $\widehat{\boldsymbol{X}}^{\dagger }$ is the adjoint eigenfunction as defined in (3.10). Knowing the Landau coefficients, the systems of equations (3.18) at different order $O(A^{m})$ can be solved sequentially for any Fourier mode ( $k=1,2,\dots$ ); see § 4.2 for the related numerical procedure.

4.1. First Landau coefficient and equilibrium solution

The nonlinear equilibrium/stationary solutions $A(t)=A_{e}$ of the Landau equation (3.15) are calculated by substituting $\text{d}A/\text{d}t=0$ . At cubic order $O(A^{3})$ , we obtain the trivial solution $A_{e}=0$ , representing the base flow, and two nontrivial solutions

where $a^{(0)}$ and $a^{(2)}$ are the growth rate and the real part of the first Landau coefficient, respectively. The non-zero solutions exist if $a^{(0)}$ and $a^{(2)}$ have opposite signs. Depending on the sign of $a^{(0)}$ and $a^{(2)}$ , two types of bifurcations can be classified: subcritical bifurcation ( $a^{(0)}<0$ and $a^{(2)}>0$ ) for which the base flow is linearly stable and supercritical bifurcation ( $a^{(0)}>0$ and $a^{(2)}<0$ ) for linearly unstable flows. The phase equation (3.16) provides information on whether the bifurcation is stationary ( $b^{(0)}={\it\omega}=0$ ) or oscillatory ( $b^{(0)}={\it\omega}\neq 0$ ), corresponding to pitchfork and Hopf bifurcations, respectively. It may be noted that the cubic Landau equation does not yield any stable equilibrium solution for the subcritical case.

The finite-amplitude equilibrium solution for perturbation flow fields (number density, velocity and granular temperature) at cubic order $O(A^{3})$ can be written as

where $\boldsymbol{X}^{[1;1]}$ , $\boldsymbol{X}^{[0;2]}$ , $\boldsymbol{X}^{[2;2]}$ , $\boldsymbol{X}^{[1;3]}$ and $\boldsymbol{X}^{[3;3]}$ are the fundamental mode, the distortion of mean flow, the second harmonic, the distortion of the fundamental and the third harmonic, respectively; the subscripts $r$ and $i$ denote the real and imaginary parts, respectively. For a stationary instability ( ${\it\omega}=b^{(0)}=0$ ), ${\it\theta}=2{\rm\pi}x/{\it\lambda}_{x}$ (note that $k_{x}=2{\rm\pi}/{\it\lambda}_{x}$ , where ${\it\lambda}_{x}$ is the wavelength along the $x$ -direction), the above equation simplifies to It should be noted that the nonlinear equilibrium solution has four different modal contributions ( $\boldsymbol{X}^{[0;2]}$ , $\boldsymbol{X}^{[2;2]}$ , $\boldsymbol{X}^{[1;3]}$ and $\boldsymbol{X}^{[3;3]}$ ), in addition to the fundamental mode, at cubic order. The numerical computation of these modal amplitudes is considered in § 4.2.

Addition of (4.2) to the base-state solution ( $\boldsymbol{X}^{0}$ ) yields the nonlinear solution for flow fields:

While presenting results in § 5.1 we will compare the above nonlinear solution with its linear counterpart

which involves only the fundamental mode scaled by its equilibrium amplitude.

4.2. Spectral-based numerical method

To determine the equilibrium amplitude (4.1) and the associated nonlinear pattern (4.3), the differential equations (3.18) for different Fourier modes ( $k=1,2,3,\dots$ ) at different orders in the perturbation amplitude $O(A^{m})$ must be solved, along with the first Landau coefficient (3.22). These equations are solved numerically by employing the spectral collocation scheme (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988), the details of which can be found in Shukla & Alam (Reference Shukla and Alam2011a ,Reference Shukla and Alam b ) and Alam & Shukla (Reference Alam and Shukla2013). For the sake of completeness of the present paper, we briefly recall the salient features of our numerical scheme.

At $O(A)$ we need to solve the linear stability equations in the form of the following generalized differential eigenvalue problem (put $k=m=1$ into (3.18)):

along with homogeneous boundary conditions:

The four ordinary differential equations in (4.6) are discretized along the $y$ -direction by implementing the staggered-grid spectral collocation scheme (Alam & Nott Reference Alam and Nott1998; Alam & Shukla Reference Alam and Shukla2013) which uses $M$ th order Chebyshev polynomials as basis functions. More specifically, while the mass balance equation is collocated at Gauss points (which are the zeros of the Chebyshev polynomials), the momentum and granular energy equations are collocated at Gauss–Lobatto points (which are the extrema of the Chebyshev polynomials). Since the Gauss points do not include the two boundary points, there is no need to impose any ‘artificial’ boundary condition for density (Alam & Nott Reference Alam and Nott1998); such ‘artificial’ boundary conditions are, however, needed when all equations are collocated at Gauss–Lobatto points, resulting in one spurious eigenvalue. We use a transformation

to transfer the variables from the physical grid $y\in (0,H)$ to the spectral grid ${\it\xi}\in (-1,1)$ . Upon implementing the spectral discretization of the derivative, the differential eigenvalue problem (4.6) is transformed into a generalized matrix eigenvalue problem:

where $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}=\unicode[STIX]{x1D644}$ are square matrices of order $(4M+3)$ , with $M$ being the degree of the Chebyshev polynomial. Equation (4.9) has been solved by the QZ-algorithm of MATLAB software, yielding all eigenvalues $c^{(0)}$ and associated ‘discrete’ linear eigenfunctions $\boldsymbol{X}^{[1;1]}=(n^{[1;1]},u^{[1;1]},v^{[1;1]},T^{[1;1]})$ . The same numerical procedure is followed to solve the adjoint eigenvalue problem (3.10), (3.11), yielding the adjoint eigenfunctions $\hat{\boldsymbol{X}}^{\dagger }=(\hat{n}^{\dagger },{\hat{u}}^{\dagger },\hat{v}^{\dagger },\hat{T}^{\dagger })$ and related eigenvalues.

The general form of equations for $k\neq 1$ and $m\geqslant k$ (namely (3.18)) can be written as $\unicode[STIX]{x1D647}_{km}\boldsymbol{X}^{[k;m]}=\boldsymbol{G}_{km}$ , which at $O(A^{2})$ leads to two equations:

for the distortion to mean flow ( $\boldsymbol{X}^{[0;2]}$ ) and the second harmonic ( $\boldsymbol{X}^{[2;2]}$ ), respectively. The inhomogeneous systems of ordinary differential equations (4.10), (4.11) are discretized using the same staggered-grid spectral collocation scheme as explained above. The resulting set of algebraic equations along with boundary conditions can be put into the form

where $\boldsymbol{X}$ is the discrete version of the respective modal amplitude function ( $\boldsymbol{X}^{[0;2]}$ or $\boldsymbol{X}^{[2;2]}$ ). To deal with the dense, nearly ill-conditioned coefficient matrix $\unicode[STIX]{x1D63C}$ , we have employed the singular-value decomposition (Golub & van Loan Reference Golub and van Loan1989) to solve (4.12).

All the above information is needed to calculate the first Landau coefficient from (3.22) which involves the ratio of two definite integrals:

To achieve spectral accuracy and superior convergence, these integrals are evaluated using a Gauss–Chebyshev quadrature formula (Shukla & Alam Reference Shukla and Alam2011a ; Alam & Shukla Reference Alam and Shukla2013) which uses a Chebyshev polynomial as an interpolating polynomial. Having calculated $c^{(2)}$ from (4.13), we need to solve the equations for (i) the distortion to the fundamental mode (put $k=1$ and $m=3$ into (3.18)) and (ii) the third harmonic (put $k=3=m$ into (3.18)):

respectively. Both involve second-order ordinary differential equations that can be put into the same form as (4.12) via spectral discretization of derivatives. The numerical solution of (4.14), (4.15) completes the solution of the nonlinear stability problem up to cubic order $O(A^{3})$ , which provides leading-order nonlinear corrections to the linear stability problem.

Before presenting results on nonlinear stability, we need to validate the numerical solution of the adjoint stability problem. It should be noted that the eigenvalues of the adjoint operator (3.10) are complex conjugates of those of the linear stability operator (3.7). To check this, we show the variation of the real part of the least stable eigenvalue ( $c^{(0)}=a^{(0)}+\text{i}b^{(0)}$ ) of linear (circles) and adjoint (stars) operators with wavenumber for two values of $(F,S)=(6,100)$ and (16, 400) in figures 3(a) and (b), respectively; it should be noted that the imaginary part of the eigenvalue is identically zero ( $b^{(0)}=0$ ) for both cases, implying that the instability is stationary. It is clear from these figures that the numerically obtained linear and adjoint eigenvalues overlap with each other, which, in turn, validates the adjoint stability problem (3.10), (3.11).

Figure 3. Variation of the real part of the least stable mode of linear (circles) and adjoint (stars) eigenvalues with wavenumber ( $k_{x}$ ) for ( $F,S$ ) (a) $=(6,100)$ and (b) $=(16,400)$ , with $e=0.9$ .

To check the convergence of the numerical method with respect to the degree of the Chebyshev polynomial (i.e. the number of collocation points $M$ ), the variation of the growth rate ( $a^{(0)}$ ) with $k_{x}$ , with other parameters being the same as in figure 3(a), is shown in figure 4; the stars, triangles and circles denote results for collocation points of $M=20$ , 50 and 100, respectively. It is seen that the growth rates are well matched for $M=50$ and 100; we have also verified that a similar agreement holds for the first Landau coefficient (not shown) too. Therefore, for all results presented in this paper the degree of the Chebyshev polynomial ( $M$ ) is fixed to 50.

Figure 4. The effect of the number of collocation points ( $M$ ) on the least stable mode. Parameter values are the same as in figure 3(a).

5.1. Numerical results: from ‘full’ convection to a state of ‘floating’ convection

The variation of the real part of the first Landau coefficient $a^{(2)}$ (denoted by stars) with wavenumber $k_{x}$ , along with the corresponding variation of the growth rate $a^{(0)}$ (circles) of the least stable mode, is shown in figure 6(a). The layer depth at rest is set to $F=6$ and the shaking strength to $S=100$ , with the restitution coefficient being $e=0.9$ (which is appropriate for glass balls). For this problem, we have verified that the instabilities are stationary, i.e. $b^{(0)}=b^{(2)}=0$ , leading to a stationary/pitchfork bifurcation from the underlying density-inverted Leidenfrost state. It is seen that $a^{(2)}<0$ for two ranges of $k_{x}$ over which the Leidenfrost state is linearly unstable ( $a^{(0)}>0$ ): (i) small wavenumbers of $k_{x}\in (0.0366,0.0535)$ and (ii) moderate wavenumbers of $k_{x}\in (0.0703,0.1748)$ . Therefore, the flow admits supercritical pitchfork bifurcations for the above two ranges of $k_{x}$ . The corresponding variation of the nonlinear equilibrium amplitude $A_{e}$ (see (4.1)) is displayed in figure 6(b), confirming the existence of $A_{e}$ over the above two ranges of $k_{x}$ . It is seen that $A_{e}$ increases from both the lower and upper bifurcation points, located at $k_{xl}\approx 0.0366$ and $k_{xu}\approx 0.1748$ , respectively; it should be noted further in figure 6(b) that the equilibrium solution emanating from $k_{xl}$ terminates at $k_{x}\approx 0.0535$ . Figure 6(a) confirms that there is a small window of wavenumbers $k_{x}\in (0.0535,0.07)$ over which $a^{(0)}>0$ but $a^{(2)}>0$ , and hence the equilibrium solutions (4.1) do not exist there; the vertical line in figure 6(b) corresponds to the location $k_{x}\approx 0.0703\equiv k_{d}$ at which $a^{(2)}=0$ (see figure 6 a), and hence the equilibrium amplitude diverges at $k_{x}=k_{d}$ .

Figure 6. Variations of (a) the growth rate of the least stable mode $a^{(0)}$ (circles) and the real part of the first Landau coefficient $a^{(2)}$ (stars) with wavenumber $k_{x}$ for a shaking strength of $S=100$ . (b) Variation of the equilibrium amplitude $A_{e}$ with $k_{x}$ ; see the text for details. The number of particle layers at rest is $F=h_{0}/d=6$ and the restitution coefficient is $e=0.9$ .

Figure 6(a) shows that the first Landau coefficient $a^{(2)}$ undergoes a sign reversal at a wavenumber of $k_{x}\approx 0.0535\equiv k_{r}$ at which $a^{(2)}\rightarrow \pm \infty$ (i.e. a jump discontinuity), and an equilibrium solution does not exist at this location. The former is, in fact, a signature of resonance, as we discuss below.

5.2. Qualitative comparison of patterns with experiments

A limited set of experiments has been conducted for a specific case for which a transition from the Leidenfrost state to convection is known (Eshuis et al. Reference Eshuis, van der Weele, van der Meer, Bos and Lohse2007) to occur. The primary purpose of these experiments is to make a qualitative comparison with the nonlinear predictions and to highlight the underlying differences, if any, between experiment and theory. The experimental set-up and method, similar to those of Eshuis et al. (Reference Eshuis, van der Weele, van der Meer, Bos and Lohse2007), are detailed in Ansari & Alam (Reference Ansari, Alam, Yu, Dong, Yang and Luding2013). The monodisperse spherical glass balls (diameter $d=1.0$  mm, with a standard deviation of less than 5 %; see the snapshot in figure 14(a) to have an idea about the size uniformity of these balls) are filled into a quasi-two-dimensional rectangular Plexiglas container of length $L_{x}/d=80$ and depth $D/d=5.5$ (see figure 1), with an open top. This glass-ball-filled container is mounted on an electromagnetic shaker and is shaken harmonically, (2.1), along the vertical direction. The initial layer depth is set to $F=6$ and the shaking amplitude to $a_{s}/d=3$ ; the shaking strength (2.4) is varied by varying the frequency of the shaker, $f_{s}$ , at a rate of $0.01\ \text{Hz s}^{-1}$ .

Figure 14. (a) A raw snapshot of particles undergoing convective motion in experiments with parameter values of $F=6$ , $S=120$ , $a_{s}/d=3$ and $L_{x}/d=80$ . Instantaneous (spatially) coarse-grained maps of (b) velocity vectors and (c) density.

When the shaking strength is increased from the Leidenfrost regime, a convection pattern with counter-rotating rolls emerges, as shown in the raw snapshot of particles in figure 14(a); the value of the shaking strength is $S=120$ . The corresponding velocity-vector map is displayed in figure 14(b), which shows four counter-rotating rolls. The pictures such as in panel (a) were captured with a high-speed camera at a framing rate of 1000 fps (frames per second) and were subsequently analysed with the aid of PIV (particle image velocimetry) software (Ansari & Alam Reference Ansari, Alam, Yu, Dong, Yang and Luding2013). It should be noted that panel (b) represents an instantaneous velocity map implying that we have not made any temporal averaging but only spatial averaging with a variable spatial box of $16\ \text{pixels}\times 16\ \text{pixels}$ to $64\ \text{pixels}\times 64\ \text{pixels}$ . (Approximately 10 pixels corresponds to one particle diameter $d=1\ \text{mm}$ .) The analogue of figure 14 for a higher shaking strength of $S=150$ is displayed in figure 15; it should be noted that the velocity vectors are superimposed on the raw image in figure 15(a).

Figure 15. The same as figure 14 but for $S=150$ . In panel (a) the raw image and the PIV velocity vectors are superimposed.

In addition to calculating the PIV velocity field, the coarse-grained density map has been obtained from the digitized version of the particle snapshots by calculating an ‘effective’ density field in terms of the average light intensity over a box of $10\ \text{pixels}\times 10$ pixels. This is displayed in figures 14(c) and 15(b) which represent the instantaneous spatially averaged density field at $S=120$ and 150, respectively. In each case, a Gaussian filter with a width of $2{\it\sigma}$ ( ${\it\sigma}$ is the standard deviation) was used to smoothen the density field. It is seen in figure 14(c) that there is a band of dense particle clusters around a vertical height of $y/d\approx 6.5$ , punctuated by dilute regions. On comparing figure 14(c) with its raw image in (a), we find that the particles shoot up through these relatively dilute regions. Once the Leidenfrost state becomes unstable at a critical shaking strength ( $S=S_{c}(F)$ ), the particles shoot up through the dense bed, creating relatively dilute channels at different horizontal locations, and rain down on two sides of each channel, thus creating two dense clusters at a lower elevation as well as two counter-rotating rolls. A comparison of figures 14(c) and 15(b) indicates that increasing the shaking strength results in a more modulated density field, with the row of dense clusters being located at a relatively higher elevation ( $y/d\approx 9$ in figure 15 c), as expected.

For parameter values as in the experiments of figure 14, the theoretical nonlinear patterns of velocity and number density are shown in figures 16(a) and (b), respectively. The length of the box has been chosen such that it fits four counter-rotating rolls (i.e. $L_{x}/d=2{\it\lambda}_{x}/d=4{\rm\pi}/k_{x}\approx 78.5$ ), as in the experiments of figure 14(b). The theoretical density pattern in (b) looks only qualitatively similar (with respect to large-scale features) to that in the experiments in figure 14(c); the latter pattern appears to be more tenuous, having complicated small-scale features. More importantly, the dense particle clusters are located at a much higher level in theory than in experiments. (Similar observations can be made from a comparison between figures 15 and 17 at $S=150$ .) The latter quantitative discrepancy between theory and experiment may be improved partially by relaxing the theoretical assumption of zero energy loss at the front and back walls of the quasi-2D container, which is assumed to be periodic in the lateral directions. This energy loss will reduce the overall energy input to the vibrating grains and can therefore be tied with a lower value of the ‘effective’ shaking strength ( $S$ ), which corresponds to the dense layer being located at a relatively lower height (see figure 2 c). A comparison between figures 16(b) and 17(b) for $S=120$ and 150, respectively, supports the above conclusion on the location of the dense layer. In any case, the quantitative discrepancies between theory and experiment provide enough impetus to relook at (i) the underlying assumptions, (ii) the missing ingredients and (iii) the phenomenological parameters in the constitutive relations of the hydrodynamic model. Additional experiments with an exhaustive analysis are recommended for a future work.

Figure 16. Theoretical predictions of nonlinear patterns for (a) the velocity and (b) the density fields in the ( $x,y$ )-plane for $F=6$ , $S=120$ and $k_{x}=0.16$ . The restitution coefficient for glass beads is taken to be $e=0.9$ .

Figure 17. The same as figure 16 but for $S=150$ .

6.1. Effects of bulk viscosity, Prandtl number and shear work

So far, all results have been presented by setting the bulk viscosity ( ${\it\zeta}=0$ ) to zero in the constitutive model (2.13). To understand the effects of ‘non-zero’ bulk viscosity on the nonlinear results presented in § 5, we consider two expressions for the dimensionless bulk viscosity:

where ${\it\phi}$ ( $=n{\rm\pi}d^{3}/6$ ) is the volume fraction of the particles. Equation (6.1) holds for a dense gas of elastically colliding spheres, and (6.2) for hard disks (Jenkins & Richman Reference Jenkins and Richman1985a ). The analogues of figure 6(a) with (6.1) and (6.2) are displayed in figures 18(a) and (b), respectively; other constitutive expressions, the boundary conditions as well as the base flow (the bulk viscosity does not appear in the base-flow equations (2.24), (2.25)) remain the same. It is clear from figures 18(a,b) that the bulk viscosity does not have a noticeable effect on the least stable mode $a^{(0)}$ and the first Landau coefficient $a^{(2)}$ ; only the equilibrium amplitude $A_{e}$ at smaller $k_{x}$ ( ${\sim}0.05$ ) appears to have been affected marginally.

Figure 18. The effect of the bulk viscosity on the variations of $a^{(0)}$ (circles), $a^{(2)}$ (stars) and $A_{e}$ (inset) with wavenumber $k_{x}$ at $S=100$ : (a) ${\it\zeta}\neq 0$ (6.1); (b) ${\it\zeta}\neq 0$ (6.2); other parameters are the same as in figure 6.

Figure 19. (a) The density dependence of the ‘effective’ Prandtl number for a dense gas of elastically colliding hard spheres; the circles represent the exact expression obtained from the kinetic theory (see text for details) and the solid line is a polynomial fit (6.3). (b) The effect of a density-dependent Prandtl number (6.3) on the variations of $a^{(0)}$ (circles), $a^{(2)}$ (stars) and $A_{e}$ (inset) with wavenumber $k_{x}$ . Other parameters are the same as in figure 6.

As elaborated in § 2.2, the shear viscosity ( ${\it\mu}$ ) was taken proportional to the thermal conductivity (2.13), with a constant value for the ‘effective’ Prandtl number ( $Pr=1.7$ ), which provided a good match (Eshuis et al. Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010) for the critical value of the shaking strength $S$ at the onset of convection (see figure 5). However, the Prandtl number should depend on the density as expected from the kinetic theory of dense gases. By employing the explicit expressions of the shear viscosity and the thermal conductivity for nearly elastic ( $e\sim 1$ ) spheres (Jenkins & Richman Reference Jenkins and Richman1985a ,Reference Jenkins and Richman b ), we have obtained an expression for the Prandtl number which is well fitted via a ninth-degree polynomial in ${\it\phi}$ (volume fraction of particles):

This is shown in figure 19(a) – the Prandtl number is almost constant near the dense limit, but decreases to a minimum at some intermediate density ( ${\it\phi}\sim 0.2$ ) and sharply increases thereafter, reaching its dilute limit. It should be noted that in figure 19(a) we have kept the average Prandtl number as 1.7, denoted by the dashed line, as in all results presented in § 5. The stability results with (6.3) are displayed in figure 19(b) and the results are marginally affected in comparison with figure 6 which was calculated with a constant $Pr=1.7$ . For example, the maximum growth rate is slightly larger for variable $\mathit{Pr}$ ; the upper and lower critical points are shifted to $k_{xu}\approx 0.179$ and $k_{xl}\approx 0.0351$ (in contrast to $k_{xu}\approx 0.175$ and $k_{xl}\approx 0.0365$ for constant $\mathit{Pr}$ ), respectively; the locations of 2:1 resonance and $a^{(2)}=0$ (i.e. diverging amplitude) are shifted to $k_{r}\approx 0.0531$ and $k_{d}\approx 0.0671$ , (in contrast to $k_{r}\approx 0.0535$ and $k_{d}\approx 0.0703$ for constant $\mathit{Pr}$ ), respectively. Another functional form for (6.3) was also tested (which holds for hard disks), but the results remain almost identical to figure 19(b). We have further verified that decreasing the average Prandtl number (say, to $Pr_{av}=1$ , with or without density variation) does not qualitatively change our results, but the growth rate of the least stable mode as well as the range of $k_{x}$ over which the instability occurs increases with decreasing $Pr_{av}$ . The latter observations (not shown) are expected since the effective shear viscosity ( ${\it\mu}\propto Pr{\it\kappa}$ ) decreases with decreasing $Pr_{av}$ , leading to more destabilization of the unstable mode. These results are omitted for the sake of brevity.

In most stability works on classical Rayleigh–Bénard convection as well as in other free convection flows (e.g. thermal plumes), it is often assumed that the shear work or the viscous dissipation term ( ${\it\Phi}_{vis}$ as in (2.22)) in the energy equation (2.19) is negligible (Gebhart et al. Reference Gebhart, Jaluria, Mahajan and Sammakia1988; Lakkaraju & Alam Reference Lakkaraju and Alam2007). This was also assumed to hold for granular convection in previous linear stability analysis (Eshuis et al. Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010). As mentioned in § 2, we have retained this term in our nonlinear analysis and the results presented so far refer to non-zero shear work. To check the ansatz of Eshuis et al. (Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010, Reference Eshuis, van der Weele, Alam, van Gerner, van der Höf, Kuipers, Luding, van der Meer and Lohse2013), we have repeated the calculations of figure 6 by setting the shear work to zero, and the results are shown in figure 20. A comparison of figures 6(a,b) and 20 confirms that there are no qualitative changes to either linear or nonlinear stability predictions irrespective of whether ${\it\Phi}_{vis}=0$ or ${\it\Phi}_{vis}\neq 0$ . There exist minor quantitative differences, however; for example, the Landau coefficient $a^{(2)}$ is slightly smaller (near the upper critical point $k_{x}\sim 0.17$ ) with ${\it\Phi}_{vis}=0$ .

Figure 20. Effects of zero shear work ( ${\it\Phi}_{vis}=0$ , (2.22)) on variations of $a^{(0)}$ (circles), $a^{(2)}$ (stars) and $A_{e}$ (inset) with wavenumber $k_{x}$ . Other parameters are the same as in figure 6, for which ${\it\Phi}_{vis}\neq 0$ .

6.2. Effects of free-surface boundary conditions

In all results presented so far, we have used Dirichlet conditions at the free surface, i.e. the perturbation temperature and velocities are zero ( $T^{\prime }=0=u^{\prime }=v^{\prime }$ ). This follows from the expectation that the perturbations decay exponentially with height such that

hold at $y=H$ . The above boundary conditions are equivalent to imposing Dirichlet conditions, $(u^{\prime },v^{\prime },T^{\prime })\sim \exp (-k_{x}H)$ , at $y=H$ . In fact, this decaying solution can be obtained by analysing the asymptotic behaviour of the linear stability equations at $y\rightarrow \infty$ as in the work of Forterre & Pouliquen (Reference Forterre and Pouliquen2002) on inclined chute flow – we omit the related details. Instead of Dirichlet conditions, one can also impose the Neumann conditions,

at the free surface $y=H$ . The linear stability calculations of figure 3(a) are repeated with the above two boundary conditions and the results are shown in figure 21(a). The triangles and circles represent (6.4) and (6.5), respectively, and the pluses represent results with $T^{\prime }=0=u^{\prime }=v^{\prime }$ at $y=H$ as in figure 3(a). As expected, the growth rates of the least stable modes for the boundary conditions (6.4) almost overlap with those for $T^{\prime }=0=u^{\prime }=v^{\prime }$ . However, the Neumann boundary conditions (6.5) yield slightly larger growth rates for the modes around $k_{x}=0.1$ . It is noteworthy that the overall features, including the range of $k_{x}$ over which the flow is unstable ( $a^{(0)}>0$ ), remain virtually unchanged irrespective of the choice among the above three free-surface boundary conditions.

Figure 21. (a) The effect of different free-surface boundary conditions on the growth rate of the least stable mode. (b) The effect of the zero heat-flux boundary condition on the variations of $a^{(0)}$ (circles), $a^{(2)}$ (stars) and $A_{e}$ (inset) with $k_{x}$ , see text for details. Other parameters are as in figure 6.

For the base-state calculations in § 2.4, the zero heat-flux boundary condition was used at the free surface; therefore, we checked our nonlinear stability calculations of figure 6(a) with $\text{d}T^{\prime }/\text{d}y(y=H)=0$ by keeping the remaining boundary conditions the same. For this case, the adjoint boundary conditions at the free surface (3.11) need to be changed to

The results are shown in figure 21(b), with other parameters as in figure 6. It is clear that there are minor changes in the wavenumber variations of (i) the first Landau coefficient $a^{(2)}$ , (ii) the location of the resonance point (marked by the sign reversal of $a^{(2)}$ ) and (iii) the equilibrium amplitude $A_{e}$ . These observations hold even if we use the boundary conditions (6.4). Therefore, our nonlinear results presented in § 5 are robust with respect to the choice of free-surface boundary conditions.

6.3. Range of validity of the present nonlinear analysis

In the above we have presented the nonlinear equilibrium amplitudes of the unstable modes for a large range of control parameters ( $k_{x}$ and $S$ ) based on the first Landau coefficient. As elaborated previously (Shukla & Alam Reference Shukla and Alam2011a ,Reference Shukla and Alam b , Reference Shukla and Alam2013; Alam & Shukla Reference Alam and Shukla2013), readers are again reminded that the present nonlinear results may not be valid far away from the bifurcation point since the radius of convergence of the Stuart–Landau series (3.15), (3.16) cannot be inferred from just one Landau coefficient. This requires the calculation of higher-order Landau coefficients, which is beyond the scope of the present paper.

Another related issue is the occurrences of nonlinear resonance phenomena, namely, the mean-flow and 2:1 resonances as discussed in § 5.1.1. This also puts restrictions on the range of validity of our nonlinear analysis: the present single-mode analysis is not valid at resonance points (Shukla & Alam Reference Shukla and Alam2011b ). The resonant mode interactions, in terms of coupled Landau equations, must be developed to obtain equilibrium amplitude near such singular points, and this too is beyond the scope of the present paper.

Ideally, one should choose a simple base flow (such as uniform shear flow) to address the convergence issues. The second Landau coefficient was calculated recently by Shukla & Alam (Reference Shukla and Alam2013) for the vorticity-banding modes in plane Couette flow. They found that the cubic- and quintic-order solutions for the equilibrium amplitude match with each other near the bifurcation point. They further showed that the quintic terms have a stabilizing effect far away from the supercritical bifurcation point in the sense that the equilibrium amplitude (far away from the critical point) decreases with the addition of the second Landau coefficient. The latter is good news since this implies that the radius of convergence is not small, at least for the specific case of plane Couette flow. For the present problem as well, we need to calculate a few higher-order Landau coefficients to check the convergence of the Stuart–Landau series in the future; however, it might take years of concerted effort to resolve this issue.