On the shape and size of granular roll waves

Abstract

This paper describes, from a theoretical point of view, the appearance and characteristics of granular roll waves in chute flow, and the maximal size these waves can attain for a given influx of material into the system. Granular roll waves are steady travelling wave solutions of the generalized Saint-Venant equations for flowing granular matter, appearing when the Froude number $Fr$ of the incoming flow exceeds a critical value, $Fr>Fr_{cr}$. We focus upon the phase space of the corresponding dynamical system, where the roll waves take the form of a stable limit cycle around an unstable fixed point; this limit cycle gives precise information on the size and periodicity of the roll wave. It is found that, for any given value of $Fr$, the limit cycle cannot become arbitrarily large because it is constrained by a homoclinic loop in phase space. Roll waves of larger amplitude can be generated by increasing the Froude number $Fr$.

1. Introduction

The study of roll waves has a long tradition in fluid mechanics (Cornish 1910, 1934; Dressler 1949; Whitham 1974; Needham & Merkin 1984; Merkin & Needham 1986; Kranenburg 1992; Balmforth & Mandre 2004; Fowler 2011). The granular counterpart of these waves, which we study here, has been the focus of many investigations during the past few decades, both experimentally and theoretically (Savage & Hutter 1989; Forterre & Pouliquen 2003; Forterre 2006; Di Cristo et al. 2009; Gray & Edwards 2014; Razis et al. 2014; Edwards & Gray 2015; Viroulet et al. 2018; Razis, Kanellopoulos & van der Weele 2019; Fei et al. 2021). The mainstream approach to describe roll waves is via a generalized version of the Saint-Venant equations for shallow water, adequately modified to account for the special characteristics of flowing granular matter.

Roll waves (see figure 1) are, together with the uniform flow and the granular monoclinal wave, among the few stable travelling waveforms that can occur in fully dynamic granular chute flow (Razis et al. 2019). The specification ‘fully dynamic’ is needed here, since in granular flows (in contrast to usual fluids) one may also encounter travelling waveforms which exhibit stopping regions, i.e. regions where the grains come to a temporary standstill and then have to be set into motion again via a mechanism that involves erosion and deposition. These intriguing waveforms have been studied in recent years by Edwards & Gray (2015), Edwards et al. (2017, 2021), Rocha, Johnson & Gray (2019) and Russell et al. (2019) but in the present paper we focus on fully dynamic travelling waves.

Figure 1. A periodic train of granular roll waves, travelling with wave speed $c$ on a chute with inclination angle $\zeta$. This is the principal travelling waveform for Froude numbers exceeding the critical value $Fr_{cr}$ given by (1.2). With $Fr_{b}$ we denote the local Froude number at the crest of the waves, just before the shock front, whereas $Fr_{a}$ is the Froude number in the trough immediately after the shock front ($Fr_{b} > Fr_{a} > Fr_{cr}$). The upper plot depicts the wave profile using identical scaling for the $x$ and $z$ axes, as it will appear in experiment, while the lower plot has been rescaled to demonstrate the features of the wave profile more clearly.

Another interesting class of waveforms that falls outside the scope of this paper is that of standing waves, the most iconic example of which is the so-called granular jump (Hákonardóttir & Hogg 2005; Boudet et al. 2007; Gray & Cui 2007; Johnson & Gray 2011; Faug 2015; Faug et al. 2015; Méjean, Faug & Einav 2017; Viroulet et al. 2017; Méjean et al. 2020; Kanellopoulos, Razis & van der Weele 2021). The spectrum of possible waveforms becomes even wider when one considers the fact that granular flows are often multi-disperse, with particles of varying sizes and shapes, introducing segregation effects and stratification. Indeed, granular flows can support an immense variety of waveforms, only matched by the multitude of wave phenomena encountered in the fields of Newtonian and non-Newtonian fluids combined.

Which travelling waveform will prevail in practice depends primarily on the value of the Froude number at the inlet of the chute. The Froude number is defined as

(1.1)\begin{equation} Fr(x,t)=\frac{\bar{u}(x,t)}{ \sqrt{ h(x,t) g \cos \zeta} } , \end{equation}

where $\bar {u}(x,t)$ denotes the depth-averaged velocity of the granular sheet, $h(x,t)$ stands for its height, $g$ is the gravitational acceleration and $\zeta$ is the inclination angle of the chute. As far as the waveforms studied in the present paper are concerned, i.e. travelling waves in the fully dynamic regime for monodisperse granular matter, the uniform flow and the monoclinal wave (being a travelling shock wave connecting two uniform flows) are stable for $Fr(x,t)< Fr_{cr}$ (Razis, Kanellopoulos & van der Weele 2018; Razis et al. 2019), whereas roll waves are found for $Fr(x,t)>Fr_{cr}$ (Forterre & Pouliquen 2003; Gray & Edwards 2014; Razis et al. 2014; Edwards & Gray 2015; Razis et al. 2019). The critical Froude number is given by the expression

(1.2)\begin{equation} Fr_{cr}=\tfrac{2}{3} \left( 1 - \varGamma \right) \end{equation}

and depends only on the Froude offset $\varGamma$, which plays a role in the friction between the chute and the granular sheet (Forterre & Pouliquen 2003; Edwards et al. 2017; Kanellopoulos 2021). This parameter $\varGamma$ is determined by the specifics of the experimental set-up and the granular material, see § 2.

Figure 1 shows the anatomy of a typical train of roll waves propagating downwards with constant shock speed $c$. Each member of the wave train consists of an elongated rising flank, gradually (and for the most part exponentially) building up to a peak, followed by a sudden drop in height, i.e. a shock front linking the peak to the lowest point of the wave. Immediately afterwards the wave begins to rise again. The positions ‘before’ and ‘after’ the shock front are denoted by the subscripts $b$ and $a$, respectively. The main objective of this paper is to describe, from a theoretical point of view, the reason for the appearance and the stability of the roll waves, together with their specific shape and periodicity, and the fact that they cannot grow beyond a certain maximal size depending on the influx of material into the system.

The paper is organized as follows. In § 2.1 we state the governing equations, namely, the mass and momentum balance – being two coupled partial differential equations (PDEs) for $h(x,t)$ and $\bar {u}(x,t)$ – together with appropriate constitutive relations that capture the special physical properties of flowing granular matter. In § 2.2 we non-dimensionalize the problem and identify the relevant dimensionless numbers. In § 2.3 it is shown how the governing equations, if one limits oneself to the case of travelling waves (which is the case we are interested in), reduce to a single ordinary differential equation (ODE) of second order for $h(\xi )$, where $\xi =x-ct$ is the so-called travelling wave variable. Subsequently, in § 2.4 we write this second-order ODE in the form of a dynamical system consisting of two coupled first-order ODEs. In the same subsection we determine the fixed points and the nullclines of this dynamical system.

In § 3 we then focus on the specific solution of the dynamical system that corresponds to the periodic train of roll waves. In phase space this solution takes the form of a stable limit cycle. In § 3.1 we determine the parameter range (especially with respect to the Froude number and the wave speed) where roll waves may expected to be found; § 3.2 deals with the birth of the limit cycle, which takes place via a Hopf bifurcation; in § 3.3 we discuss the particular shape of the limit cycle, which is intimately connected to the shape of the nullclines; and in § 3.4 we explain why the limit cycle (for a given Froude number of the incoming flow) cannot grow beyond a certain maximum size. It turns out to be an immediate consequence of a homoclinic bifurcation taking place in the vicinity of the limit cycle. In § 3.5 we show that these results remain qualitatively the same when we choose different values of $\zeta$ (the inclination angle of the chute) and $\varGamma$ (the offset of the Froude number), but that the quantitative changes may be significant. Throughout the paper, we continuously indicate how our findings in phase space are linked to the birth, shape and maximum amplitude of the roll waves in physical space. This is not only of interest from a fundamental point of view but also of practical relevance. Indeed, the present study may be anticipated to contribute to our understanding of roll wave formation in large scale geophysical flows (such as debris flows or rock avalanches) and to answer the critical question concerning the maximal amplitude that these waves can attain. The spontaneous formation of waves can seriously amplify the destructive potential of such flows.

In § 4 we study how the roll wave pattern reacts to an increase of the Froude number of the incoming flow, which in experiment can be regulated by enhancing the inclination angle $\zeta$ or by further opening the inlet at the top of the chute. It is found that the roll waves increase in amplitude with growing Froude number. Finally, in § 5 we summarize our main results and discuss the prospects for testing them experimentally.

2. Governing equations

2.1. The granular Saint-Venant equations and the associated constitutive relations

In analogy with the shallow water approximation for normal fluids, and justified by the fact that the wavelength of the waves under consideration is orders of magnitude larger than the thickness of the sheet, which in turn is markedly larger than the wave amplitude (cf. figure 1), we describe the dynamics of the flowing granular sheet in terms of its height $h(x,t)$ and depth-averaged velocity $\bar {u}(x,t)$. The latter is defined as

(2.1)\begin{equation} \bar{u}(x,t) = \frac{1}{h(x,t)}\int_{0}^{h(x,t)} u(x,z,t) \,{\rm d}z, \end{equation}

where $u(x,z,t)$ is the detailed velocity profile, which also depends on the ordinate $z$. In the description, any velocity variations in the $z$-direction are averaged out, while flow variations in the cross-wise direction are ignored altogether, hence both $h(x,t)$ and $\bar {u}(x,t)$ depend only on $x$ and $t$, cf. figure 1. Furthermore, we treat the granular material as an incompressible fluid (Savage & Hutter 1989), assuming constant density throughout the flowing sheet.

The height $h(x,t)$ and the depth-averaged velocity $\bar {u}(x,t)$ are governed by a system of two coupled, nonlinear PDEs (Gray & Edwards 2014; Razis et al. 2014; Edwards & Gray 2015; Razis et al. 2018; Edwards et al. 2019; Razis et al. 2019; Edwards et al. 2021; Kanellopoulos 2021), namely, the mass conservation

(2.2)\begin{equation} \frac{\partial h}{\partial t} + \frac{\partial}{\partial x} (h \bar{u}) = 0, \end{equation}

and the momentum balance

(2.3)\begin{equation} \frac{\partial}{\partial t} (h \bar{u}) + \frac{\partial}{\partial x} (h \bar{u}^{2}) = gh\sin\zeta - \frac{\partial}{\partial x}\left(\frac{1}{2}gh^{2}\cos\zeta\right) - \mu(h, \bar{u})gh \cos\zeta + \frac{\partial}{\partial x}\left(\nu h^{3/2}\frac{\partial \bar{u}}{\partial x}\right). \end{equation}

The left-hand side of (2.3) is the Stokesian derivative of the depth-averaged momentum per unit volume, divided by the constant density. The terms on the right-hand side of (2.3) represent the various depth-integrated forces (per unit volume and divided by the density) acting on the flowing material: (i) the gravity component along the $x$-direction (the direction of the flow); (ii) the force arising from variations in $h(x,t)$, i.e. the negative gradient of the depth-averaged pressure; (iii) the frictional force exerted by the chute bed, modelled as a (height- and velocity-dependent) friction coefficient $\mu (h,\bar {u})$ multiplied by the normal reaction force acting on the granular sheet; and (iv) the viscous-like diffusive term arising from depth averaging the in-plane stresses in the sheet.

The last two terms of (2.3) contain constitutive relations that are specific for granular matter, and it is here that the above equations deviate from their counterparts for Newtonian fluids such as water. For the friction coefficient $\mu (h,\bar {u})$ we adopt the modified expression introduced by Edwards et al. (2017) that is valid for granular matter in the fully dynamic regime, i.e. without stopping regions

(2.4)\begin{equation} \mu(h,\bar{u}) = \tan \zeta_1 + (\tan \zeta_2 - \tan \zeta_1)\left(1 + \frac{\beta h}{\mathcal{L} (Fr+\varGamma)}\right)^{{-}1}. \end{equation}

Here, $\zeta _1$, $\zeta _2$, $\mathcal {L}$, $\beta$ and $\varGamma$ are experimental parameters obtained from independent experiments (Edwards et al. 2021). The angles $\zeta _{1}$ and $\zeta _{2}$ denote the two inclinations delimiting the interval in which uniform granular flows are possible: for $\zeta < \zeta _{1}$ the granular sheet remains at rest, whereas for $\zeta > \zeta _{2}$ it flows downwards in an accelerated fashion. The length scale $\mathcal {L}$ has its origins in the functional form of the experimental fit to the critical angle curves that determine the friction force, as described by Pouliquen & Forterre (2002). It represents a characteristic flow thickness (usually ranging from $1$ to $2$ particle diameters) over which the transition between $\zeta _{1}$ and $\zeta _{2}$ takes place in the above friction law, (2.4). Its exact value depends on the properties of the grains as well as on the bed roughness (Pouliquen & Forterre 2002; Edwards et al. 2017). The parameter $\beta$ and the Froude number offset $\varGamma$ are closely related to each other and are given by (Forterre & Pouliquen 2003)

(2.5)\begin{equation} Fr = \beta \frac{h}{h_{{stop}}(\zeta)}-\varGamma, \end{equation}

where $h_{{stop}}(\zeta )$ denotes the thickness of the deposition layer that is left on the chute (inclined at an angle $\zeta$) after a uniform flow of height $h$ has passed over it. A first version of (2.5) (without the $\varGamma$ offset) was derived by Pouliquen (1999) by experimentally measuring the Froude number as a function of the ratio $h/h_{{stop}}$ for four different systems of glass beads, revealing that this ratio is independent of the bead size and the roughness of the bed. The Froude offset $\varGamma$ was introduced by Forterre & Pouliquen (2003) in order to fit the experimental results concerning the mean velocity for different materials (sand and glass) into a single law, thus formulating (2.5). As stated in GDR-MiDi (2004), one must be cautious in using this relation for $h$-values close to $h_{stop}$, since the fit involved may lose its accuracy in this limit. In the present study this poses no problem, however, since for travelling waves without stopping regions one has to ensure anyhow that $h$ remains everywhere well above $h_{stop}$; so any possible ambiguities in the value of $\varGamma$ close to $h_{stop}$ can hardly affect our results. For an accurate comparison with the literature, we note that the factor $1/\sqrt {\cos \zeta }$ was absent from the definition of the Froude number $Fr$ in many of the pioneering studies.

In Edwards et al. (2017) a parameter $\beta _{*}$ was also introduced, indicating the value of the Froude number corresponding to a sheet thickness $h=h_{*}(\zeta )$ with $h_{{stop}}(\zeta )< h_{*}(\zeta )< h_{{start}}(\zeta )$, where $h_{{start}}(\zeta )$ denotes the height for which a static layer is mobilized when the inclination is increased to the angle $\zeta$. Physically, $h_{*}(\zeta )$ is the height of the thinnest (and therefore slowest) possible steady uniform flow which leaves a deposit of smaller thickness $h_{{stop}}(\zeta )$. The parameter $\beta _{*}$ is also affected by the Froude offset $\varGamma$, since it must satisfy $\beta _{*}>\beta -\varGamma$ (Edwards et al. 2019). The granular flow is guaranteed to remain fully dynamic as long as $Fr(x,t) \geqslant \beta _{*}$ everywhere along the chute (Edwards et al. 2017).

Another parameter that is characteristic for granular flow is the diffusive coefficient $\nu$ in the fourth term of the right-hand side of the momentum balance (2.3). We adopt the analytical expression derived by Gray & Edwards (2014)

(2.6)\begin{equation} \nu = \nu(\zeta) = \frac{2 \mathcal{L} \sqrt{g} \sin\zeta }{9 \beta \sqrt{\cos\zeta}}\gamma(\zeta),\quad \textrm{where}\ \gamma(\zeta) = \frac{\tan \zeta_2 - \tan\zeta}{\tan\zeta - \tan \zeta_1}. \end{equation}

This expression is physically relevant in the interval $\zeta _1 < \zeta < \zeta _2$, decreasing from infinity at $\zeta = \zeta _1$ to zero at $\zeta = \zeta _2$. For inclination angles outside this interval the expression (2.6) becomes negative and cannot be used.

For completeness, we note that the advective term $\partial _x (h \bar {u}^{2})$ on the left-hand side of the momentum balance (2.3) is often multiplied by a shape factor $\alpha$ (of order $1$) to obtain an optimal correspondence with experiments (Saingier, Deboeuf & Lagrée 2016; Lagrée et al. 2017). However, for shallow granular flows with relatively small Froude number, as in the current paper, the solutions are quite insensitive to the value of this shape factor (Saingier et al. 2016; Viroulet et al. 2017) and therefore we choose to suppress it. Moreover, for shallow granular flows a certain amount of slipping motion at the bottom is almost unavoidable, even on a rough chute; this means that the difference in velocity between the lower and the upper layers is less pronounced, or in other words, that the velocity profile is not too far from plug flow. In this context, we note that in a study by Forterre (2006) (see in particular figure 1 in that paper) it was found that a simple plug flow with shape factor $\alpha =1$ fitted the experimental data better than the Bagnold velocity profile, for which $\alpha$ would be $5/4$ (GDR-MiDi 2004; Börzsönyi, Hasley & Ecke 2005; Gray & Edwards 2014).

2.2. Non-dimensional form of the equations

At this point, it is opportune to non-dimensionalize the above set of equations. We will take the height $h_0$ of the ‘uniform base flow’ as the characteristic length scale. This is the thickness of the sheet if the material released from the inlet is left to distribute itself evenly along the entire span of the chute. Practically, in the absence of obstacles, this height will be attained at some point not far from the inlet. The corresponding depth-averaged velocity $\bar {u}_0$ defines the characteristic time scale $h_0/\bar {u}_0$. So our non-dimensional variables are $\tilde {x}= x /h_{0}$ and $\tilde {t} = t / (h_{0}/\bar {u}_{0})$, and the height and depth-averaged velocity are rescaled as follows (we use tildes throughout to denote the non-dimensional quantities)

(2.7a,b)\begin{equation} \tilde{h}= \frac{h}{h_{0}},\quad \tilde{\bar{u}}=\frac{\bar{u}}{\bar{u}_{0}}. \end{equation}

The continuity equation (2.2) then becomes

(2.8)\begin{equation} \frac{\partial \tilde{h}}{\partial \tilde{t}} + \frac{\partial}{\partial \tilde{x}} (\tilde{h} \tilde{\bar{u}}) = 0 ,\end{equation}

and the momentum balance (2.3) takes the form

(2.9)\begin{equation} Fr_{0}^{2} \left(\frac{\partial}{\partial \tilde{t}} (\tilde{h} \tilde{\bar{u}}) + \frac{\partial}{\partial \tilde{x}} (\tilde{h} \tilde{\bar{u}}^{2})\right) = \tilde{h}\left(\tan\zeta-\mu(\tilde{h},\tilde{u})\right)- \tilde{h}\frac{\partial \tilde{h}}{\partial \tilde{x}} + \frac{Fr_{0}^{2}}{Re}\frac{\partial}{\partial \tilde{x}}\left(\tilde{h}^{3/2}\frac{\partial \tilde{\bar{u}}}{\partial \tilde{x}}\right), \end{equation}

where $Fr_0$ is the Froude number of the uniform base flow that is supposed to exist far upstream

(2.10)\begin{equation} Fr_{0}=\frac{\bar{u}_{0}}{\sqrt{h_{0}g\cos\zeta}}, \end{equation}

and $Re$ is given by the following dimensionless expression, being the granular analogue of the Reynolds number (introduced by Gray & Edwards (2014) on the basis of the so-called $\mu (I)$-rheology), again defined in terms of the base flow conditions

(2.11)\begin{equation} Re=\frac{\bar{u}_{0}\sqrt{h_{0}}}{\nu(\zeta)}. \end{equation}

For consistency, also the friction coefficient $\mu$ is now expressed as a function of the new dimensionless parameters

(2.12)\begin{equation} \mu(\tilde{h},\tilde{\bar{u}}) = \tan \zeta_1 + (\tan \zeta_2 - \tan \zeta_1)\left(1 + \frac{\beta h_{0}\tilde{h}}{\mathcal{L} \left( \varGamma+Fr_{0} ( \tilde{\bar{u}}/\sqrt{\tilde{h}} ) \right)}\right)^{{-}1}.\end{equation}

Here, we note that the height $h_0$ and the velocity $\bar {u}_0$ appearing in the above expressions, are in fact related. For the uniform base flow all derivatives with respect to $x$ and $t$ vanish, and the momentum equation (2.3) reduces to a balance between gravity and friction, $\tan \zeta = \mu (h,\bar {u})$, or equivalently (in dimensionless parameters) $\tan \zeta = \mu (\tilde {h},\tilde {\bar {u}})$. Inserting $\tilde {h}=1$ and $\tilde {\bar {u}}=1$ in the latter equality, using the expression (2.12) for $\mu (\tilde {h},\tilde {\bar {u}})$ (and (2.10) for the Froude number $Fr_0$ appearing there), we get the following relation between $\bar {u}_0$ and $h_0$:

(2.13)\begin{equation} \bar{u}_0 = \frac{\beta \sqrt{g \cos \zeta}}{\mathcal{L} \gamma(\zeta)}\,h_0^{3/2} - \varGamma \sqrt{g \cos \zeta}\,h_0^{1/2}. \end{equation}

With this, the Froude number $Fr_0$ and the granular Reynolds number $Re$ can be expressed as functions of $h_0$ alone, and are in fact seen to be related to each other (Kanellopoulos 2021)

(2.14a,b)\begin{equation} Fr_0 = \frac{\beta h_0}{\mathcal{L} \gamma(\zeta)}-\varGamma,\quad Re =\frac{\mathcal{L} \gamma(\zeta) \sqrt{g \cos \zeta}}{\beta \nu(\zeta)}\,Fr_0(Fr_0+\varGamma). \end{equation}

The above expression for $Fr_0$, which is fully compatible with (2.5) for $h=h_0$ and $h_{stop}(\zeta ) = \mathcal {L} \gamma (\zeta )$, shows that the Froude number of the base flow is an increasing function of $h_0$. That is, for fixed $\zeta$, thicker base flows correspond to higher Froude numbers.

2.3. Travelling waves: reducing the problem to an ODE

The periodic train of granular roll waves is a travelling waveform, which propagates without change of shape and at a constant shock speed $\tilde {c}$, and therefore its height and depth-averaged velocity will have the general form $\tilde {h}(\tilde {x},\tilde {t}) = \tilde {h}(\tilde {x}-\tilde {c}\tilde {t}) = \tilde {h}(\tilde {\xi })$ and $\tilde {\bar {u}}(\tilde {x},\tilde {t}) = \tilde {\bar {u}}(\tilde {x}-\tilde {c}\tilde {t}) = \tilde {\bar {u}}(\tilde {\xi })$. Inserting this form in the governing equations (2.8)–(2.9), these become ODEs. Specifically, the dimensionless mass balance equation (2.8) takes the form

(2.15)\begin{equation} -\tilde{c}\tilde{h} + \tilde{h}\tilde{\bar{u}} ={-}K, \end{equation}

where the non-negative integration constant $K=\tilde {h}(\tilde {c}-\tilde {\bar {u}})$ corresponds to the dimensionless flux of material in the co-moving frame per unit width of the channel. Given the fact that for the base flow, which will be established at some distance from the inlet (depending among other things on the inclination angle $\zeta$), we have $\tilde {h}=1$ and $\tilde {\bar {u}}=1$, this constant is equal to $K=\tilde {c}-1$ along the whole chute.

Now, by solving (2.15) with respect to $\tilde {\bar {u}}$, we find the depth-averaged velocity $\tilde {\bar {u}} = \tilde {c}-K\tilde {h}^{-1} = \tilde {c}-(\tilde {c}-1)\tilde {h}^{-1}$ and its derivatives $\tilde {\bar {u}}_{\tilde {\xi }},\tilde {\bar {u}}_{\tilde {\xi }\tilde {\xi }}$ in terms of $\tilde {h}$. Substituting these expressions into the non-dimensional momentum balance (2.9), the velocity is eliminated and we obtain the following second-order ODE for $\tilde {h}(\tilde {\xi })$:

(2.16)\begin{equation} \frac{{\rm d}^{2} \tilde{h}}{{\rm d} \tilde{\xi}^{2}} - \frac{1}{2 \tilde{h}} \left(\frac{{\rm d} \tilde{h}}{{\rm d} \tilde{\xi}} \right)^{2} + \frac{ Re \tilde{h}^{3/2} }{ Fr_0^{2} (\tilde{c}-1) } \left[ \left( \frac{ Fr_0^{2} (\tilde{c}-1)^{2} }{\tilde{h}^{3}} -1 \right) \frac{{\rm d} \tilde{h}}{{\rm d} \tilde{\xi}} + \tan \zeta - \mu(\tilde{h}) \right]=0, \end{equation}

where also the friction coefficient $\mu$ is now expressed exclusively in terms of $\tilde {h}$

(2.17)\begin{equation} \mu(\tilde{h}) = \tan \zeta_1 + \frac{\tan \zeta_2 - \tan \zeta_1}{1 + \displaystyle\frac{\gamma(\zeta)\tilde{h}(Fr_0+\varGamma)}{\varGamma+Fr_0(\tilde{c} \tilde{h}-\tilde{c}+1) \tilde{h}^{{-}3/2}}}. \end{equation}

Equation (2.16) first appeared in a study by Gray & Edwards (2014) on granular roll waves and also featured in Razis et al. (2019), where the transition from granular monoclinal waves to roll waves was discussed. The distinguishing difference here is that the friction coefficient (2.17) is enriched by the inclusion of the Froude offset parameter $\varGamma$ (Edwards et al. 2017; Kanellopoulos 2021).

Roll waves are nonlinear, finite-amplitude travelling waveforms that emerge as an instability to the uniform flow. Close to their birth, when their amplitude is still small, these waves are well described by a linear approximation of the governing equations (2.8)–(2.9) around the uniform solution $(\tilde {h}_0,\tilde {\bar {u}}_0) = (1,\tilde {\bar {u}}_0)$. This linear analysis yields the wave speed $\tilde {c} = 1 \pm 1/Fr_0 + \textrm {h.o.t.}$, (where h.o.t. stands for higher order terms), or in dimensional form $c = u_0 \pm \sqrt {h_0 g \cos \zeta } + \textrm {h.o.t.}$, in which one may recognize the familiar expression for the velocity (in the laboratory frame) of ‘gravity waves’ propagating on the surface of a fluid moving at speed $u_0$. As the roll waves gain height, as a result of the aforementioned instability, the nonlinear terms become increasingly important, in such a way as to mitigate their growth and eventually arrest it altogether, producing a fully ripened roll wave pattern. This sequence of events is in fact at the core of the present study, and will be treated in detail in the next sections.

2.4. The dynamical system formulation

Equation (2.16) can equivalently be written as two coupled first-order ODEs, forming a dynamical system

(2.18a)\begin{gather} \frac{{\rm d}\tilde{h}}{{\rm d} \tilde{\xi}} = \tilde{s} = f(\tilde{h},\tilde{s}), \end{gather}

(2.18b)\begin{gather}\frac{{\rm d} \tilde{s}}{{\rm d} \tilde{\xi}} = {\frac{\tilde{s}^{2}}{2\tilde{h}}} - \frac{Re \tilde{h}^{3/2}}{Fr_0^{2}(\tilde{c}-1) } \left[ \left( \frac{Fr_0^{2} (\tilde{c}-1)^{2}}{\tilde{h}^{3}} -1 \right) \tilde{s} + \tan \zeta - \mu(\tilde{h}) \right] = g(\tilde{h},\tilde{s}). \end{gather}

The fixed points of the above dynamical system lie at the intersections of the nullclines in the $(\tilde {h},\tilde {s})$ phase space. The nullclines are defined as the curves where the direction field ${\rm d}\tilde {s}/{\rm d}\tilde {h} = g(\tilde {h},\tilde {s})/f(\tilde {h},\tilde {s})$ attains infinite and zero slope, respectively,

(2.19a)\begin{gather} f(\tilde{h},\tilde{s}) = \tilde{s} = 0, \end{gather}

(2.19b)\begin{gather}g(\tilde{h},\tilde{s}) = {\frac{\tilde{s}^{2}}{2\tilde{h}}} - \frac{Re \tilde{h}^{3/2}}{Fr_0^{2}(\tilde{c}-1) } \left[ \left( \frac{Fr_0^{2} (\tilde{c}-1)^{2}}{\tilde{h}^{3}} -1 \right) \tilde{s} + \tan \zeta - \mu(\tilde{h}) \right] = 0. \end{gather}

In order to clearly distinguish between the two, we will call $f(\tilde {h},\tilde {s})=0$ the ‘vertical-motion nullcline’ and $g(\tilde {h},\tilde {s})=0$ the ‘horizontal-motion nullcline’.

Wherever the nullclines intersect, we have a fixed point. From (2.19a) we immediately see that any fixed point necessarily corresponds to a flat profile of the flowing sheet, with $\tilde {s}=0$ (zero slope). Substituting this (2.19a) into (2.19b), we recover the condition for uniform flow $\tan \zeta = \mu (\tilde {h})$, which we already discussed in connection with (2.13). This condition yields two different values of $\tilde {h}$, namely $\tilde {h}=1$ and $\tilde {h} =\tilde {h}_{+} < 1$; the latter is not a constant, but rather a function of the parameters $\tilde {c}$, $Fr_0$ and $\varGamma$, as can be seen from the following expression that relates these quantities:

(2.20)\begin{equation} \tilde{c} = \frac{1}{1-\tilde{h}_{+}}\,\frac{1}{Fr_0} \left((Fr_0+\varGamma)(1-\tilde{h}_{+}^{5/2})-\varGamma(1-\tilde{h}_{+}^{3/2}) \right). \end{equation}

The associated fixed points $(\tilde {h},\tilde {s})=(1,0)$ and $(\tilde {h}_{+},0)$ (see figure 2) correspond to a uniformly flowing sheet with constant thickness and velocity. The fixed point $(1,0)$ represents the familiar base flow near the inlet, whereas the second fixed point $(\tilde {h}_{+},0)$ represents a thinner uniform flow. Both these uniform flows can occur simultaneously on the chute, linked by a travelling shock front, forming a ‘monoclinal’ flood wave or an undular variation of this waveform (Razis et al. 2018, 2019; Kanellopoulos 2021).

Figure 2. The fixed points $(\tilde {h}_{+},0)$ and $(1,0)$ (indicated by open dots) in the $(\tilde {h},\tilde {s})$ phase plane for $Fr_0=0.7$ and $\tilde {c}=2.42457 (=\tilde {c}_{H}-0.004)$, together with the vertical-motion nullcline $\tilde {s}=0$ (blue) and the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$ (magenta, given by (2.26)). The black dotted line depicts the almost vertical contour of zero divergence of the vector field, which is discussed in § 3.3. The grey background arrows denote the direction field of the dynamical system and the large black arrows signify this field evaluated on the nullclines. The system parameters are as in table 1 and the inclination angle is chosen to be $\zeta =36^{\circ }$.

Table 1. Values of the system parameters as measured by Edwards et al. (2017) for a flow of carborundum particles over a bottom layer of glass beads glued on a chute. The inclination angle of the chute $\zeta$ must be chosen between $\zeta _1$ and $\zeta _2$ to get a steady, non-accelerated flow; in the present paper we take $\zeta = 36^{\circ }$ as a typical case, but for completeness we will also investigate other values on the interval ($\zeta _1,\zeta _2$), see § 3.5. For the given parameter values, the critical Froude number where the uniform flow gives way to roll waves is $Fr_{cr} = (2/3) (1- \varGamma ) = 0.40$, cf. (1.2).

To determine the local stability properties of either fixed point we consider the Jacobian matrix

(2.21)\begin{equation} J(\tilde{h},\tilde{s}) =\left.\begin{pmatrix} \dfrac{\partial f}{\partial \tilde{h}} & \dfrac{\partial f}{\partial \tilde{s}} \\ \dfrac{\partial g}{\partial \tilde{h}} & \dfrac{\partial g}{\partial \tilde{s}} \end{pmatrix}\right\vert_{(\tilde{h},\tilde{s})} = \left.\begin{pmatrix} 0 & 1 \\ \dfrac{\partial g}{\partial \tilde{h}} & \dfrac{\partial g}{\partial \tilde{s}} \end{pmatrix}\right\vert_{(\tilde{h},\tilde{s})}, \end{equation}

where the elements of the matrix must be evaluated at the fixed point in question. The corresponding eigenvalues are

(2.22)\begin{equation} \lambda_{a,b}(\tilde{h},\tilde{s}) = \frac{1}{2} \left( \textrm{Tr}\, J \pm \sqrt{(\textrm{Tr}\,J)^{2}-4\,\textrm{det} \,J} \right) =\frac{1}{2} \left( \frac{\partial g}{\partial \tilde{s}} \pm \sqrt{ \left( \frac{\partial g}{\partial \tilde{s}} \right)^{2} + 4 \frac{\partial g}{\partial \tilde{h}} }\,\right). \end{equation}

The above eigenvalues can be evaluated analytically, yielding four expressions (two for each fixed point) depending on $\tilde {c}$, $Fr_0$, $\varGamma$ and the angles $\zeta _1$, $\zeta$ and $\zeta _2$. The corresponding formulas are

(2.23)\begin{align} \lambda_{a,b}(\tilde{h},\tilde{s}) &= \frac{1}{2Fr_0^{2}(\tilde{c}-1)\tilde{h}^{3/2}} \left[ Re \tilde{h}^{3} + Fr_0^{2} (\tilde{c}-1) \left( \tilde{h}^{1/2}\tilde{s}-Re(\tilde{c}-1) \right)\right.\nonumber\\ &\quad \pm \left\{Re^{2} \tilde{h}^{6} + Fr_0^{4} (\tilde{c}-1)^{2} \left( Re^{2} (\tilde{c}-1)^{2}+4Re(\tilde{c}-1) \tilde{h}^{1/2}\tilde{s} -\tilde{h}\tilde{s}^{2}\right)\right.\nonumber\\ &\quad +2Fr_0^{2} Re (\tilde{c}-1) \tilde{h}^{3} \left(\left(4\tilde{h}^{1/2}\tilde{s}-Re(\tilde{c}-1)\right)\right.\nonumber\\ &\quad \left.\left.\left.+\tilde{h}^{1/2}\left(3\mu(\tilde{h})-3\tan\zeta+2\tilde{h}\mu'(\tilde{h})\right) \right) \right\}^{1/2} \right], \end{align}

where $\mu (\tilde {h})$ is given by (2.17). Since these expressions are rather long, we will simply evaluate their numerical values whenever we address the issue of the stability, with specific system parameters ($\varGamma, \zeta _1, \zeta _2$). For definiteness, in the rest of this study we adopt the parameter values measured experimentally by Edwards et al. (2017), concerning the flow of carborundum particles on a bed of glass beads glued to the chute. The values are given in table 1. This constitutes a realistic setting to demonstrate our results. Other sets of parameters are also possible (see e.g. table for the flowing sand studied by Rocha et al. 2019), yet the changes this produces are only quantitative, while the qualitative trend is preserved.

Before we proceed, let us provide the explicit expressions of the trace and the determinant of the Jacobian matrix at the fixed point $(\tilde {h},\tilde {s})=(1,0)$, which we will use in the analysis to follow

(2.24)\begin{gather} \textrm{Tr} \,J(1,0) \equiv \left.\frac{\partial g}{\partial \tilde{s}}\right\vert_{(1,0)} = \frac{Re \left( 1 - Fr_0^{2} (\tilde{c}-1)^{2} \right)}{Fr_0^{2} (\tilde{c}-1)} , \end{gather}

(2.25)\begin{gather}\displaystyle \textrm{det} \,J(1,0) \equiv\left. -\frac{\partial g}{\partial \tilde{h}}\right\vert_{(1,0)} ={-}\frac{(\tan\zeta -\tan\zeta_1)(\tan\zeta_2-\tan\zeta) Re \left( \left( \tilde{c}-\frac{5}{2} \right) Fr_0-\varGamma \right)}{(\tan\zeta_2-\tan\zeta_1) (Fr_0+\varGamma) Fr_0^{2} (\tilde{c}-1) } . \end{gather}

For the same reason, we give the analytical expression for the horizontal-motion nullcline associated with (2.19b). This is a quadratic equation for $\tilde {s}$, with the following two roots:

(2.26)\begin{equation} \tilde{s}_{N} = \frac{Re\tilde{h}^{5/2}}{Fr_0^{2} (\tilde{c}-1)} \left\{\frac{Fr_0^{2}(\tilde{c}-1)^{2}}{\tilde{h}^{3}}-1 \pm \sqrt{ \left( \frac{Fr_0^{2}(\tilde{c}-1)^{2}}{\tilde{h}^{3}} -1 \right)^{2} + \frac{2Fr_0^{2} (\tilde{c}-1)}{Re \tilde{h}^{5/2}} \left( \tan \zeta - \mu(\tilde{h}) \right) } \right\}. \end{equation}

The nullcline in question consists of two branches, as can be seen in figure 2, where $\tilde {s}_{N}$ as given by (2.26) is represented by the magenta curves. These branches will play a central role in our analysis.

In closing, we note that the condition (2.19b) for the horizontal-motion nullcline corresponds to (2.16) without the second-order derivative ${\rm d}^{2}\tilde {h}/{\rm d}\tilde {\xi }^{2}$. This means that on the nullcline $\tilde {s}_{N}$ the second-order ODE governing the travelling waveforms reduces locally to a first-order ODE. Or equivalently, since the term ${\rm d}^{2}\tilde {h}/{\rm d}\tilde {\xi }^{2}$ stems from the diffusive term in the momentum balance (2.9), on the horizontal-motion nullcline the full problem locally reduces to its inviscid counterpart.

3. Roll waves as a limit cycle

3.1. On the parameter range where roll waves may be expected

The local stability of both fixed points with $\varGamma =0$ has been addressed in Razis et al. (2019). As we will see, the inclusion of the $\varGamma$ offset makes a marked quantitative difference yet it does not give rise to any qualitative changes. We will show this in detail for $\varGamma = 0.4$ (the value of table 1). In § 3.5 we will also discuss other values of $\varGamma$, corresponding to different experimental set-ups, which lead to the same conclusion.

Regarding the fixed point $(\tilde {h},\tilde {s})=(\tilde {h}_{+},0)$ given by (2.22) we find that the eigenvalues $\lambda _{a,b}(\tilde {h}_{+},0)$ are real and of opposite sign for all relevant parameter values, showing that $(\tilde {h}_{+},0)$ is always a saddle point.

The stability properties of the fixed point $(\tilde {h},\tilde {s})=(1,0)$ are more intricate. In the interval $\beta _{*}\leqslant Fr_0 < Fr_{cr}$ (where the uniform flow is stable, and where also the monoclinal flood wave is found) it is an unstable node or, in the immediate vicinity of $Fr_{cr}$, an unstable spiral (Razis et al. 2019; Kanellopoulos 2021). The monoclinal flood wave manifests itself as a heteroclinic connection in phase space between the two fixed points $(1,0)$ and $(\tilde {h}_{+},0)$. For $Fr_0 \geqslant Fr_{cr}$ (i.e. the regime where roll waves exist) the fixed point $(1,0)$ can be either an unstable spiral, a centre, a stable spiral or a stable node, depending on the parameter values. Figure 3, for the typical value $Fr_0=0.7$, shows the changing characteristics of the fixed point $(1,0)$ for varying values of the non-dimensional wave speed $\tilde {c}$.

Figure 3. The trajectories in the complex plane, for decreasing $\tilde {c}$ and at a fixed Froude number $Fr_0=0.7$, of the two eigenvalues $\lambda _{a,b}(1,0)$ of the fixed point $(\tilde {h},\tilde {s})=(1,0)$. In the insets above the main plot we sketch how the local stability characteristics of the fixed point $(1,0)$ change as the eigenvalues vary. The five schematic pictures below the main plot depict how, simultaneously, the orbits surrounding the fixed point $(1,0)$ interact with those around the neighbouring saddle point $(\tilde {h}_{+},0)$ denoted by ‘S’. The stable and unstable manifolds of the latter are indicated in red and blue, respectively. The letter ‘H’ stands for Hopf bifurcation (where $(1,0)$ is momentarily a centre point and gives birth to a limit cycle) and ‘hom’ denotes the critical event where the stable and unstable manifolds of the saddle connect to form a closed homoclinic loop. In the region between ‘H’ and ‘hom’ (dashed) the fixed point $(1,0)$ is surrounded by a stable limit cycle, see the middle picture below the main plot; this limit cycle corresponds to a periodic train of roll waves.

From the discussion below (2.15) we know that $\tilde {c}=1+K>1$. Physically, this means that the travelling waveforms of the granular sheet are always propagating faster than the carrying base flow. For the specific case of roll waves, we can also give an upper bound for $\tilde {c}$. Since roll waves manifest themselves as limit cycles in phase space, and since a limit cycle cannot be formed around a saddle point, it follows that the determinant of the Jacobian matrix $J(1,0)$, given by (2.25), must be positive in the roll wave regime. This means that the factor $(\tilde {c}-5/2)Fr_0-\varGamma$ in the numerator of this determinant must be negative, since all other factors in the expression are positive. This gives us the upper bound $\tilde {c}<5/2+\varGamma /Fr_0$. So roll waves must be sought for $\tilde {c}$ values within a sub-interval of

(3.1)\begin{equation} 1<\tilde{c}<\frac{5}{2}+\frac{\varGamma}{Fr_0}. \end{equation}

Now, the positive determinant alone does not ensure the existence of a stable limit cycle. Also needed is the presence of a Hopf bifurcation, where the fixed point $(1,0)$ turns from a stable spiral into an unstable one (momentarily and locally taking the form of a centre point). The Hopf bifurcation can be found by setting the trace of $J(1,0)$ equal to zero, see (2.24): $1-Fr_0^{2} (\tilde {c}-1)^{2} =0$. Thus one finds that the Hopf bifurcation takes place at the following value of $\tilde {c}$, cf. the analyses by Gray & Edwards (2014) and Razis et al. (2019):

(3.2)\begin{equation} \tilde{c}=\tilde{c}_{H}=1+\frac{1}{Fr_0}. \end{equation}

It may be noted that the equation $1-Fr_0^{2} (\tilde {c}-1)^{2} =0$ also gives a second solution $\tilde {c}=1-1/Fr_0$, yet this one is discarded since it lies outside the interval (3.1) of admissible $\tilde {c}$ values. The local stability properties of the fixed point $(1,0)$ in the neighbourhood of the Hopf bifurcation at $\tilde {c}_{H}$, while keeping $Fr_0$ fixed, show that $(1,0)$ is an unstable spiral for $\tilde {c}<\tilde {c}_{H}$ and a stable one for $\tilde {c}>\tilde {c}_{H}$. This means that a stable limit cycle (i.e. a roll wave solution) exists only for $\tilde {c}<\tilde {c}_{H}$. In fact, as we will see, they appear only in a very narrow interval of $\tilde {c}$-values directly below $\tilde {c}_H=1+1/Fr_0$ (cf. figure 3).

Before we proceed, let us comment upon a subtle point concerning the dimensionless wave speed $\tilde {c}_{H}$ in (3.2). This relation may give the erroneous impression that roll waves originating from a base flow with large $Fr_0$ have small wave speeds (and vice versa) since $\tilde {c}_{H}$ is a decreasing function of $Fr_0$. In reality, however, the opposite is true, as was shown by Razis et al. (2014).

To get the correct picture, one has to consider the dimensional wave speed $c_{H}$, which reads as follows (with $\bar {u}_0$ given by (2.13) and using the expression $h_0 = (\mathcal {L} \gamma (\zeta )/\beta ) (Fr_0+\varGamma )$, cf. (2.14a,b)):

(3.3)\begin{equation} c_H=\tilde{c}_{H}\bar{u}_0=(1+Fr_0)\sqrt{(Fr_0+\varGamma) \frac{\mathcal{L} \gamma(\zeta) }{\beta}g \cos\zeta}, \end{equation}

which is an increasing function of $Fr_0$. This tells us that roll waves stemming from a base flow with higher $Fr_0$ (and hence higher $h_0$, cf. (2.14a,b)) have a larger amplitude, and therefore a larger shock speed. This is also observed in practice.

In the non-dimensional context, however, if we fix the value of $Fr_0$ (and hence $h_0$) one has to follow a path through $(Fr_0, \tilde {c})$ parameter space with decreasing $\tilde {c}$ to witness the birth of roll waves, and to see their amplitude gradually grow. This will be the type of path that we follow in our analysis.

3.2. The birth and evolution of the limit cycle

Figure 3, based on the analytical expressions of (2.22), depicts the trajectories of the eigenvalues – also known as the root locus – of the fixed point $(1,0)$ for decreasing values of $\tilde {c}$, while the Froude number $Fr_0$ is kept constant.

The sketches above the central plot of figure 3 show the changing character of the fixed point in question. We also depict, in the five figures below the main plot, the simultaneously changing interplay between $(1,0)$ and the saddle $(\tilde {h}_{+},0)$, here denoted by the letter ‘S’. The roll waves occur in the narrow dashed interval to the right of the imaginary axis, where $(1,0)$ has just become an unstable spiral. In this interval, which is only a small part of the section where $(1,0)$ is an unstable spiral, the fixed point $(1,0)$ happens to be surrounded by a stable limit cycle, corresponding to a periodic train of roll waves. Figure 4 shows that, in the $(Fr_0,\tilde {c})$ parameter space, the region of stable limit cycles (grey zone) is indeed remarkably narrow: the curves along which the Hopf bifurcation and the homoclinic event take place are seen to lie very close to each other.

Figure 4. The $(Fr_0,\tilde {c})$ parameter plane, with the region in which roll waves occur shaded in grey. The upper bound of this region (denoted by ‘H’) is given by (3.2) where the Hopf bifurcation takes place, cf. figure 3, and the lower bound corresponds to the homoclinic event (‘hom’). For Froude numbers below $Fr_{cr} = (2/3) (1-\varGamma ) = 0.4$ (hatched region) no roll waves exist; the corresponding critical value of $\tilde {c}$ follows from (3.1) and equals $3.5$, so the curves ‘H’ and ‘hom’ both emanate from the critical point $(Fr_0,\tilde {c})=(Fr_{cr}, 5/2 + \varGamma /Fr_{cr})=(0.4, 3.5)$. On the scale of the main figure the two curves ‘H’ and ‘hom’ seemingly coincide, therefore we also show an inset in which we magnify the section around $Fr_0=0.7$ (which is one of the typical values used in the text). The system parameters are as in table 1 and the inclination angle of the chute is $\zeta =36^{\circ }$.

In figure 5 we present the limit cycle for three different values of $\tilde {c}$ and for fixed Froude number $Fr_0=0.7$. In figure 5(a,b), with $\tilde {c}=\tilde {c}_{H}-10^{-4}$, the limit cycle – the blue closed curve surrounding $(1,0)$ – has just been born and hence the corresponding roll wave train (plot to the right) exhibits an almost sinusoidal behaviour with relatively small amplitude. Somewhat further away from the Hopf bifurcation, at $\tilde {c}=\tilde {c}_{H}-10^{-3}$, figure 5(c,d), the limit cycle is seen to have increased in size. The amplitude and the wavelength of the corresponding roll waves have increased accordingly, and the wave pattern starts to visibly deviate from a sinusoidal one. In figure 5(e,f), at $\tilde {c}=\tilde {c}_{H}-4\times 10^{-3}$ (again further away the Hopf bifurcation), the amplitude and wavelength have grown further and the asymmetry of the limit cycle has become very pronounced, corresponding to the characteristic roll wave pattern observed in practice.

Figure 5. Phase space portraits and corresponding wave profiles for three successive values of $\tilde {c}$ at fixed $Fr_0=0.7$, as shown in the accompanying section of the $(Fr_0, \tilde {c})$ parameter plane: (a,b,g(i)) at $\tilde {c}=\tilde {c}_{H}-10^{-4}=2.42847$, the limit cycle has just been born from $(1,0)$ and is still quite symmetric. The waves are nearly sinusoidal and of small amplitude. (c,d,g(ii)) At $\tilde {c}=\tilde {c}_{H}-10^{-3}=2.42757$, the limit cycle has become bigger and starts to develop an asymmetry. The roll waves are therefore of larger amplitude and wavelength, and their shape deviates from the harmonic. (e,f,g(iii)) At $\tilde {c}=\tilde {c}_{H}-4\times 10^{-3}=2.42457$, the limit cycle has grown further and now exhibits a marked asymmetry. A practically straight segment appears on its topmost part, which is responsible for the nearly exponential shape of the rising flanks of the roll waves. The latter are now seen to have an even larger amplitude and wavelength. The background grey arrows denote the direction field of the dynamical system. The system parameters are as in table 1 and the inclination angle is $\zeta =36^{\circ }$.

A striking feature of the limit cycle is the almost straight segment on its topmost part, which (as we shall see in the next subsection, cf. figure 6) can be attributed to the horizontal-motion nullcline $\tilde {s}_N$. This almost straight segment of the limit cycle is responsible for the nearly exponentially rising flank of the roll waves. This can be seen by writing down the mathematical expression for the observed straight segment above the fixed point $(1,0)$

(3.4)\begin{equation} \tilde{s} = \frac{{\rm d}\tilde{h}}{{\rm d}\tilde{\xi}} = \tilde{s}_1 + a (\tilde{h}-1), \end{equation}

where $a$ is the small positive slope of the segment, and $\tilde {s}_1$ denotes the $\tilde {s}$-value of the segment at $\tilde {h}=1$. The above equation comprises a first-order linear ODE for $\tilde {h}(\tilde {\xi })$, which is readily solved to give

(3.5)\begin{equation} \tilde{h}(\tilde{\xi}) = 1 - \frac{\tilde{s}_1}{a} + C {\rm e}^{a \tilde{\xi}}, \end{equation}

where $C$ is the constant of integration. This expression (3.5) is exponentially rising with $\tilde {\xi }$, demonstrating that the straight segment in the $(\tilde {h},\tilde {s})$ plane translates to an exponential portion of the wave profile $\tilde {h}(\tilde {\xi })$.

Figure 6. The path traced out by the limit cycle (solid blue curve) in the upper half-plane is seen to closely follow a path in the immediate vicinity of the two branches of the horizontal-motion nullcline $\tilde {s}=\tilde {s}_N$ (magenta dashed curves, given by (2.26)). In fact, the vector field forces all trajectories in this region through a narrow channel around the left branch until they cross the contour of zero divergence (dotted), after which they start to drift apart. Subsequently, the limit cycle crosses the right branch of the horizontal-motion nullcline and curves downwards. The background grey arrows denote the direction field of the dynamical system; the parameters are as in table 1 and the inclination angle is $\zeta =36^{\circ }$.

One of the advantages of the dynamical systems approach presented here, is that it allows for the analytical determination of the wavelength of roll waves at the moment of their birth. Close to the Hopf bifurcation (i.e. when $\tilde {c} \simeq \tilde {c}_{H}$) the periodic wave train is almost perfectly sinusoidal, as seen in figure 5(b). Its wavelength is

(3.6)\begin{equation} \varLambda \simeq \varLambda_H = \frac{2{\rm \pi}}{|\lambda_{H}(1,0)|}, \end{equation}

where $|\lambda _{H}(1,0)|$ stands for the magnitude of the eigenvalues $\lambda _{a,b}(1,0)$ at the Hopf bifurcation. These eigenvalues are purely imaginary. For the parameters used here we have $\lambda _H(1,0) = \pm 0.506 i$. At the Hopf bifurcation itself, i.e. for $\tilde {c} = \tilde {c}_{H}$, the selected wavelength is precisely $\varLambda _H$, but at that moment the size of the limit cycle – and hence the amplitude of the roll waves – is still nil. In the case of figure 5(a,b), slightly beyond the Hopf bifurcation, we find numerically that the wavelength is $\varLambda = 12.45$; this is indeed close to the value at the exact Hopf bifurcation, $\varLambda _H = 12.40$, given by (3.6).

3.3. The influence of the nullclines on the shape of the limit cycle

Having examined the birth and evolution of the roll waves, we now focus on their characteristic shape, showing that this is dictated to a large extent by the form of the nullclines. Figure 5(e,f) shows that the roll waves are comprised of an elongated rising flank followed by a relatively abrupt shock front. The former corresponds to the upper part of the limit cycle, lying in the positive $\tilde {s}$ half-plane, while the latter is associated with the parabolic-like section in the negative half-plane. Along the upper part (the straight segment) the dynamics is slow, whereas the lower part is traversed very quickly.

Here, we will show that the overall shape of the limit cycle is heavily influenced by the nullclines of the dynamical system. The nullclines are the curves in phase space where the slope ${\rm d}\tilde {s}/{\rm d}\tilde {h}$ is either zero or infinite. As we have shown in § 2.4, the vertical-motion nullcline (where ${\rm d}\tilde {s}/{\rm d}\tilde {h}$ diverges) coincides with the horizontal axis $\tilde {s} =0$, see (2.19a). The horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$ (where ${\rm d}\tilde {s}/{\rm d}\tilde {h}=0$), which is given by (2.26), has a rather more intricate form and consists of two branches. In figure 6 the nullclines are shown for the specific parameter values $Fr_0=0.7$ and $\tilde {c}=\tilde {c}_{H}-0.004406$. The dotted line in this figure denotes the contour of zero divergence, i.e. the curve along which the expression

(3.7)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \begin{pmatrix} f(\tilde{h},\tilde{s}) \\ g(\tilde{h},\tilde{s}) \end{pmatrix} = \frac{\partial f(\tilde{h},\tilde{s})}{\partial \tilde{h}}+\frac{\partial g(\tilde{h},\tilde{s})}{\partial \tilde{s}} \end{equation}

becomes zero. To the left of this line the divergence (3.7) is negative, whereas to its right it is positive.

Let us first concentrate on the aforementioned nearly straight segment in the upper part of the limit cycle. This is seen to trace out a path adjacent to the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$. The limit cycle enters the upper half of the phase plane vertically, when it crosses the $\tilde {h}$-axis, but soon it has to re-adjust its orientation – always following the arrows of the vector field – when it approaches the left branch of the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$.

Along this left branch, the divergence of the vector field is negative, so all nearby trajectories are squeezed together. As we see in figure 6, they converge into a narrow channel in the proximity of the horizontal-motion nullcline $\tilde {s}_{N}$, following its course until they are catapulted into the gap between the two branches and pass through the contour of zero divergence (the dotted line in figure 6). After this point, they start to drift away from each other, with some trajectories heading upwards and others going downwards. The limit cycle belongs to the latter category. The reason for this is that, for the current parameter values, the unstable manifold of the saddle in $(\tilde {h}_{+},0)$ (which embraces the limit cycle, cf. the central sketch in the lower part of figure 3) curls downwards. This will be discussed further in the next subsection.

When it starts to curve downwards, the trajectory of the limit cycle quickly intersects the right branch of the horizontal-motion nullcline. Here, the second derivative ${\rm d}^{2}\tilde {h}/{\rm d}\xi ^{2}$ momentarily vanishes (which corresponds to an inflection point of the associated wave profile) and immediately afterwards, the trajectory bends towards the South and makes its way to the negative half-plane. It crosses the $\tilde {h}$-axis vertically, so at this point $\tilde {h}(\tilde {\xi })$ attains its maximum value, corresponding to the crest of the roll wave.

The above description, in which the limit cycle has a pronounced straight segment, holds when the limit cycle approaches its maximal size and is restrained in its growth by the horizontal-motion nullcline (and by the unstable manifold emanating from the saddle point). We illustrate this in figure 7 for three successive values of the Froude number: $Fr_0=0.7$, $Fr_0=0.9$ and $Fr_0=1.1$, and corresponding values of $\tilde {c}$ close to $\tilde {c}_{hom}(Fr_0)$, where the subscript ‘hom’ refers to the homoclinic event that takes place at this value, as we will see shortly. The three limit cycles in the left column are very close to their maximal size and the role of the nullclines is seen to be profound.

Figure 7. (a,d,g) The maximal size limit cycles in phase space (blue solid curves) along with the two branches of the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$ (magenta curves) for Froude numbers $Fr_0=0.7$, $Fr_0=0.9$ and $Fr_0=1.1$, respectively. The corresponding wave speeds are $\tilde {c}_{hom}(0.7)=\tilde {c}_{H}(0.7)-0.004406 \simeq 2.424165$, $\tilde {c}_{hom}(0.9)=\tilde {c}_{H}(0.9)-0.002679 \simeq 2.108432$ and $\tilde {c}_{hom} (1.1)=\tilde {c}_{H}(1.1)-0.00179804 \simeq 1.907293$. (b,e,h) The associated fully ripened wave profiles. (c,f,i) ‘Escape’ of the limit cycles (through the gap between the two branches of the horizontal-motion nullcline) for $\tilde {c}$-values that are just below the aforementioned critical ones, namely: $\tilde {c}= 2.424161$ for $Fr_0=0.7$, $\tilde {c}= 2.108431$ for $Fr_0=0.9$ and $\tilde {c}= 1.907292$ for $Fr_0=1.1$. The above choices of $Fr_0$ and $\tilde {c}$ are indicated in the $(Fr_0,\tilde {c})$ diagram in (j–m); the regions marked by the arrows in the first plot are given in more detail in the three magnifications to the right. The other parameters are invariably taken as in table 1 and the inclination angle is $\zeta =36^{\circ }$.

In the positive $\tilde {s}$ half-plane, the two branches of the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$ (magenta dashed curves) together constitute the nearly straight upper border of the limit cycle. The orbit is forced to follow the left branch, where the divergence of the vector field is negative, and starts to drift away on the right hand side, where the divergence of the vector field is positive.

In the negative $\tilde {s}$ half-plane, where $\tilde {s}=\tilde {s}_{N}$ runs almost vertically, this nullcline is responsible for redirecting the trajectory of the limit cycle: at its lowest point it turns from South-West to North-West. This specific point corresponds to the spot of steepest descent in the shock front of the roll wave. (In fact, it is one of the two inflection points of the wave profile, where ${\rm d}^{2}\tilde {h}/{\rm d}\xi ^{2}=0$. These inflection points are located at the intersections of the right branch of the horizontal-motion nullcline with the limit cycle. One of these is the present point of steepest descent, while the other one is the aforementioned point close to the wave crest.)

On the other hand, the intersections of the limit cycle with the vertical-motion nullcline $\tilde {s}=0$ (where the slope of the direction field diverges) correspond to the highest and lowest point of the roll wave pattern. The nullcline in question is merely the horizontal axis. That is to say, the limit cycle (as well as any other trajectory) always crosses the $\tilde {h}$-axis vertically, towards the South for $\tilde {h}>1$ and towards the North for $\tilde {h}_{+}<\tilde {h}<1$, cf. figure 6.

Among the three cases shown in figure 7, the roll wave depicted in the top row (where $Fr_0=0.7$) has the smallest amplitude and the largest wavelength. When the Froude number increases, as in the second and the third row, where $Fr_0=0.9$ and $Fr_0=1.1$, respectively, the two fixed points of the dynamical system are seen to move further away from each other. As a result, the limit cycle follows the horizontal-motion nullcline $\tilde {s}=\tilde {s}_{N}$ to larger positive values of $\tilde {s}$, i.e. it reaches larger slopes, which means that the rising flank of the roll wave has more space to develop itself. This is the reason for the continuing increase in amplitude observed in the second and third rows.

The depicted limit cycles are of ‘maximal size’ for the given Froude numbers: for values of $\tilde {c}$ below the critical ones of figure 7, and with the same numerical precision, the limit cycles break down. In the rightmost column of figure 7, the phase space trajectories are now seen to escape to infinity in the North-East direction, so the dynamical system ceases to give physically meaningful solutions for the Froude numbers in question. Physically this means that the profiles shown in figure 7 are the largest possible roll waves for the given $Fr_0$ values. This will be discussed in more detail in § 3.4.

As a side note, we observe (from the values of $\tilde {c}$ given in the caption of figure 7) that the difference between $\tilde {c}_{H}(Fr_0)$ and $\tilde {c}_{hom}(Fr_0)$ gradually diminishes for growing $Fr_0$. In the limit of large Froude numbers $Fr_0$, the interval of $\tilde {c}$-values for which roll waves exist becomes truly narrow, in agreement with the findings of Gray & Edwards (2014).

3.4. On the maximal size of the roll waves

In this subsection we explain, from the dynamical systems perspective, why roll waves cannot become arbitrarily large for fixed parameter values and a given base flow. More specifically, we will address the question: Why does the limit cycle for a given Froude number $Fr_0$ exhibit a maximal size? To answer this, we will study what happens when the parameter $\tilde {c}$ is pushed below its threshold value. In practice, the system will never spontaneously push itself beyond the threshold, and the present discussion will illustrate why this is so: the system would then enter a regime where no bounded travelling wave solutions exist.

In the rightmost column of figure 7 we present the breakdown of the limit cycle when the value of $\tilde {c}$ is decreased below the critical value $\tilde {c}_{hom}(Fr_0)$, for each of the three cases $Fr_0=0.7$, $0.9$ and $1.1$. Responding to the change in $\tilde {c}$, the limit cycle tends to expand, but in the process it runs into the ‘escape channel’ in the gap region between the two branches of the horizontal-motion nullcline $\tilde {s}_N$. In these plots, we also see that the branches of $\tilde {s}_{N}$ reach higher positive values of $\tilde {s}$ for growing Froude number, as mentioned earlier, allowing the spiralling trajectories to expand further in that direction.

The reason behind the limit cycle's escape lies in the interplay between the stable and unstable manifolds of the saddle point $(\tilde {h}_{+},0)$. In figure 8 we illustrate this for the case of $Fr_0=0.7$, where we follow a vertical path through a small part of the $(Fr_0,\tilde {c})$ parameter space.

Figure 8. The changing interplay between the stable (red) and unstable (blue) manifolds of the saddle point $(\tilde {h}_{+},0)$ for $Fr_0=0.7$, at three values of $\tilde {c}$ as shown in the depicted section of the $(Fr_0, \tilde {c})$ parameter plane. In (b,d,f), a close-up of the area of interest is given, showing how the interplay determines the fate of the limit cycle. (a,b,g(i)) For $\tilde {c}_{hom}(0.7) < \tilde {c} < \tilde {c}_{H}(0.7)$, the unstable manifold (blue curve) is seen to dive below the stable manifold (red curve) and ends on the stable limit cycle surrounding $(1,0)$. (c,d,g(ii)) For $\tilde {c} = \tilde {c}_{hom}(0.7)$, the manifolds join precisely together and form a homoclinic orbit. The maximal size limit cycle has been indicated as a dotted curve. (e,f,g(iii)) For $\tilde {c} < \tilde {c}_{hom}(0.7)$, the unstable manifold sweeps over the stable one and flies off, whereas the stable manifold now forms a heteroclinic connection between $(1,0)$ and the saddle $(\tilde {h}_{+},0)$. The background grey arrows denote the direction field in phase space; the system parameters are as in table 1 and the inclination angle of the chute is $\zeta =36^{\circ }$.

In figure 8(a,b) (top row), when the value of the wave speed is $\tilde {c}_{hom}(0.7) < \tilde {c} < \tilde {c}_{H}(0.7)$, the unstable manifold of the saddle (blue line) curls itself below the saddle's stable manifold (red line) and closes in upon the stable limit cycle around $(1,0)$. Meanwhile, the stable manifold comes down from the positive $\tilde {s}$ half-plane; it describes a wide curve around the region containing the limit cycle and finally arrives in the saddle $(\tilde {h}_{+},0)$.

By contrast, when $\tilde {c} < \tilde {c}_{hom}(0.7)$, as depicted in figure 8(e,f) (third row), the stable manifold (red) begins at the unstable spiral point $(1,0)$. It spirals outwards and ends, like before, in the saddle $(\tilde {h}_{+},0)$ (establishing a heteroclinic connection). The unstable manifold (blue) now sweeps over its stable counterpart, flying off in the North-East direction, and – importantly – opening the aforementioned ‘escape channel’. To see this more clearly, in the right column of figure 8 we have expanded the relevant strip of the phase plane around $\tilde {s}=0$. Thus the limit cycle ceases to exist. All the orbits in the region of the former limit cycle (with the exception of the saddle's stable manifold itself) eventually escape to the North-East.

In the intermediate threshold case, $\tilde {c} = \tilde {c}_{hom}(0.7)$, the unstable and stable manifolds join precisely onto each other, thus forming a homoclinic loop, see figure 8(c,d) (second row). This homoclinic loop corresponds to a granular solitary wave with an unusually high amplitude and infinite wavelength (Razis et al. 2019). One cannot expect to witness this solitary wave in practice, however, since it exists only at this exact value $\tilde {c}_{hom}(Fr_0)$, which requires an unrealistic amount of fine tuning. Moreover, the homoclinic orbit is unstable to any small perturbations.

At this critical state, the limit cycle has attained its maximal size, and its upper segment coincides with the corresponding part of the homoclinic loop. In figure 8(c,d) the limit cycle is represented by the dotted blue line, and in figure 9(a) we depict it as a solid blue line. Evidently, even at this maximum size, the limit cycle is still distinctly smaller than the homoclinic loop. This can also be seen from the associated wave profiles (the roll wave train and the solitary wave, respectively) depicted in figure 9(b). Interestingly, despite their difference in amplitude, both of these waveforms have the same wave speed $\tilde {c}_{hom}(Fr_0)$. In fact, the rising flanks of the roll wave pattern (solid blue) have exactly the same form as the corresponding part of the large solitary wave (black dashed curve). The instability of the homoclinic loop implies that also the roll wave pattern becomes unstable at this critical event.

Figure 9. (a) The maximal size limit cycle (solid blue line) and the unstable homoclinic loop (black dashed line) in phase space, cf. figure 8(c,d). (b) The corresponding roll wave train of maximal amplitude (solid blue curve) and the unstable solitary waveform (black dashed curve). The rising flanks of the roll wave train have precisely the same shape as the associated part of the solitary wave. Both these waveforms travel at dimensionless wave speed $\tilde {c}_{hom}$. The solitary waveform only exists at this very specific value of $\tilde {c}_{hom}$. For the roll wave train, this same value of $\tilde {c}_{hom}$ is associated with the maximum amplitude the roll waves can attain for the given Froude number $Fr_0$. The system parameters are as in table 1 and the inclination angle of the chute is $\zeta =36^{\circ }$; in view of the extreme fine tuning of $\tilde {c}$ needed to create the homoclinic loop, and since the associated wave profiles are very unstable, the waves presented here have been generated using a value of $\tilde {c}$ marginally larger than $\tilde {c}_{hom}(Fr_0)$.

From the physical point of view, the breakdown of the limit cycle implies that for a specific value of $Fr_0$, the amplitude of the granular roll wave cannot grow beyond a certain size, with a well-defined wave speed value $\tilde {c}_{hom}(Fr_0)$. This is in agreement with the findings of Razis et al. (2014), who determined the maximum amplitude of roll waves in the absence of viscous-like effects. To witness larger roll waves, one should open the inlet at the top of the chute further or increase the inclination angle $\zeta$, thereby raising the value of $Fr_0$. The point $(1,0)$ then corresponds to a higher uniform flow and, with the new set of parameter values, a higher roll wave will be formed. This is described below in § 4.

3.5. The influence of varying the system parameters $\zeta$ and $\varGamma$

Up to now, we have varied the parameters $\tilde {c}$ and $Fr_0$. It is natural to ask to what extent the above scenario (concerning the birth of the roll waves and the subsequent ripening to their maximal size) is affected by a change in the other system parameters. So here we discuss the inclination angle $\zeta$ and the Froude offset $\varGamma$. It will turn out that changes in these two parameters cause only quantitative changes; the scenario as a whole remains qualitatively the same.

As far as $\zeta$ is concerned, figure 10 shows how the thickness of the region between ‘H’ (the birth of the roll waves) and ‘hom’ (where they attain their maximal size) depends on the inclination angle. At the boundaries of the inclination interval, i.e. for $\zeta =\zeta _1$ and $\zeta =\zeta _2$, the thickness becomes zero. This stands to reason, because at $\zeta =\zeta _1$ the sheet is still stagnant and hence the waves have no opportunity to establish themselves. For $\zeta =\zeta _2$, on the other hand, the uniform base flow breaks down (since the flow is on the verge of turning into an uncontrollable avalanche), so here the waves cannot develop because they lose the substrate on which to do so. In the intermediate interval of $\zeta$ values, however, it is seen that the roll waves always follow the same sequence from birth (zero amplitude) to maximal size. The three detailed plots in figure 10 show that the thickness of their region of existence is largest for $\zeta \approx 42^{\circ }$. Physically, this means that (for an experimental set-up with the parameters of table 1) roll waves with the largest possible amplitude may be anticipated to be found around this inclination angle.

Figure 10. The three-dimensional $(Fr_0,\tilde {c},\zeta )$ parameter space, being an extension of the $(Fr_0,\tilde {c})$ plane shown in figure 4, with the region in which roll waves occur shaded in grey. The upper bound of this region (the green surface denoted by ‘H’, where the roll waves are born via a Hopf bifurcation) is given by (3.2); the lower bound ‘hom’ is where the roll waves attain their maximal size, which is determined by the homoclinic event as shown in figure 7. The three plots in the lower row depict the region in question – with the numerically calculated values of ‘hom’ – for three cuts through the three-dimensional plot at $Fr_{0}=0.7$, $0.9$, and $1.1$; note the different scaling along the $\tilde {c}$-axis in these plots. The system parameters are as in table 1.

To illustrate the influence of the parameter $\varGamma$, in figure 11 we show three stages of the sequence between ‘H’ and ‘hom’ for an entirely different set of system parameters (see table 2) corresponding to the experimental set-up used by Rocha et al. (2019), who studied the flow of sand on a chute with ballotini beads glued to the bottom. The inclination angle of the chute, which must be chosen between $\zeta _1$ and $\zeta _2$, is taken to be $\zeta = 34^{\circ }$, just as in the actual experiment. For the given parameter values, the critical Froude number that must be exceeded to get roll waves is $Fr_{cr} = (2/3) (1- \varGamma ) = 0.107$, cf. (1.2). The wave profiles in figure 11 have been generated for a base flow with $Fr_0=0.7$, which amply fulfils the condition $Fr_0>Fr_{cr}$. The values of the dimensionless wave speed in the three plots are (a) $\tilde {c} = \tilde {c}_H - 0.5 \times 10^{-4}$, (b) $\tilde {c} = \tilde {c}_H -7.0 \times 10^{-4}$ and (c) $\tilde {c} = \tilde {c}_H -24.28 \times 10^{-4}$, where $\tilde {c}_H=2.42857$ marks the birth of the roll waves. It is seen that the ripening of the roll waves from birth to maximal size follows the same trend as described in the previous subsection (see especially figure 5) for the system associated with the parameter values of table 1.

Figure 11. Roll wave profiles for the experimental set-up employed by Rocha et al. (2019), concerning the flow of sand on a chute. The system parameters of this set-up are given in table 2. The inclination angle of the chute is $\zeta = 34^{\circ }$ and the Froude number of the incoming flow is $Fr_0=0.7$, well above the critical value $Fr_{cr} = 0.107$ for this set-up. The depicted wave profiles correspond to three successive values of the wave speed: (a) $\tilde {c} = \tilde {c}_H - 0.5 \times 10^{-4}$, (b) $\tilde {c} = \tilde {c}_H -7.0 \times 10^{-4}$ and (c) $\tilde {c} = \tilde {c}_H -24.28 \times 10^{-4}$, where $\tilde {c}_H=2.42857$ marks the birth of the roll waves.

Table 2. The system parameters of the experimental set-up used by Rocha et al. (2019) for the flow of red sand over a bottom layer of ballottini beads glued on a chute.

For completeness, we have also checked other sets of system parameters, such as those corresponding to the flow of glass beads from the aforementioned study by Rocha et al. (2019) (for which $\varGamma =0$). In all cases studied, we found that the roll waves undergo the same qualitative sequence of events. So all the evidence suggests that the scenario is generic and quite insensitive to the precise values of $\zeta$ or $\varGamma$.

4. The response of roll waves to changes in the influx

In this section we study what happens when, starting from a situation in which one has a fully ripened roll wave of maximum size, we change the Froude number $Fr_0$. In experiment, this can be done either by changing the flux of granular material at the inlet (thereby modifying the level $h_0$ of the base flow) or by adjusting the inclination angle $\zeta$ of the chute. In order to see how the flow reacts to such a change, we have to abandon the dynamical systems approach and return to the original set of PDEs (2.2)–(2.3). After all, although the dynamical system is a powerful tool for studying the system in its equilibrium state, it provides little or no information about how this equilibrium is reached.

The fact that the only stable waveform for $Fr_0 > Fr_{cr} > \beta _{*}$ (i.e. in the fully dynamic regime) is the roll wave (Razis et al. 2019), means that – as long as we stay in this regime – any change in the inflow can only result in another roll wave. To verify this, we conduct the following numerical experiment. Starting from a fully ripened roll wave for $Fr_0=0.7$ and dimensional wave speed $c=0.228\ {\rm m}\ {\rm s}^{-1}$ (corresponding to $\tilde {c}=2.42505 \simeq \tilde {c}_{H}(0.7)-0.003525$ in non-dimensional units, in the grey zone of figure 4), we force the system into a new state – with a higher or a lower Froude number – by raising or lowering the thickness of the base flow by $0.3$ mm. We call these the positively and negatively perturbed profiles, respectively. These perturbed profiles, one at a time, are used as the initial condition $h(x,0)$ for the PDEs (2.2)–(2.3). As far as the initial velocity condition $\bar {u}(x,0)$ is concerned, we use the dimensional analogue of the mass balance (2.15) with the aforementioned value $c=0.228\ {\rm m}\ {\rm s}^{-1}$.

So we run the PDE system (2.2)–(2.3) for three different settings: (a) the reference system (without any subtractions or additions) at its proper Froude number $Fr_0=0.7$, (b) the positively perturbed system with an initial condition in which we have added $0.3$ mm to the reference profile, corresponding to $Fr_0=0.767252$ and (c) the negatively perturbed system, where we start from an initial condition that lies $0.3$ mm below the reference profile, corresponding to $Fr_0=0.588516$. The numerical scheme we employ is the method of lines (Schiesser 1991), described in detail in Razis et al. (2019).

The numerical experiment is conducted on a domain equal to one wavelength of the reference roll wave, with periodic boundary conditions. This means that also the positively and negatively perturbed profiles are confined to this same domain, which differs from the natural wavelengths they would attain on a chute of infinite length. This is admittedly a rigid aspect of the current numerical experiment, yet it is a necessary measure to keep the computer time within reasonable bounds. On a longer chute, the perturbed profiles would be prone to a slow coarsening process (Razis et al. 2014), in which the waves would tend to cluster, redistributing the granular matter in the evolving wave train and thereby altering the values of $h_0$ and $Fr_0$. It would take too much computer time to reach the steady state in that case.

Figure 12 shows our results. The dashed curves represent the initial profiles $h(x,0)$ (being the same shape for all three experiments, but superimposed on base flows of different height) while the solid curves are the numerical solutions of the PDEs after $t=10$ s. It is seen that the roll wave of the reference system (a) remains unchanged. This is as expected and – importantly – shows the stability of the solution obtained by the dynamical systems approach. The perturbed initial profiles on the other hand (b,c), after $10$ s are seen to have developed into roll waves with increased and reduced amplitude, respectively, around their given base flow levels. Meanwhile, the wave speeds of these perturbed profiles have also adjusted themselves accordingly.

Figure 12. Evolution of the roll wave profiles for different values of the influx, i.e. different values of the Froude number $Fr_0$. Depicted are the profiles after $10$ s (solid lines) together with the initial profiles (dashed lines) for: the reference Froude number $Fr_0=0.70$ (a), increased Froude number $Fr_0=0.77$ (b) and reduced Froude number $Fr_0=0.59$ (c). We have imposed periodic boundary conditions, with domain length equal to the wavelength of the reference wave. The system parameters are as in table 1 and the inclination angle of the chute is $\zeta =36^{\circ }$.

In short, within the limitations of the numerical experiment (which does not allow for a change in wavelength), the system has aptly reacted to the modified influx by forming roll waves of the correct amplitude and wave speed to comply with the new value of $Fr_0$.

5. Concluding remarks

Here, we briefly recapitulate our main findings and comment upon the feasibility of testing the results of this paper in experiment.

The central point of the paper has been to describe, via a dynamical systems view, the appearance and the ripening of roll waves in granular chute flow. We find that roll waves – which in phase space take the form of a limit cycle – are born by means of a Hopf bifurcation of the fixed point $(1,0)$. This fixed point represents the incoming base flow, and the Hopf bifurcation corresponds to its destabilization through the development of a periodic wave train. This is the onset of the so-called roll wave instability; see also the seminal paper by Gray & Edwards (2014). One of the new results gained from the dynamical systems approach has to do with the wavelength selection at the moment of birth of the roll waves: the wavelength $\varLambda$ immediately after the Hopf bifurcation can be expressed analytically in terms of the magnitude of the eigenvalue of the fixed point $(1,0)$, see (3.6).

After birth, the roll waves grow. The present approach allows us to relate the characteristic shape of the growing roll waves (or equivalently, the limit cycle) to the form of the nullclines of the dynamical system, in particular to the horizontal-motion nullcline $\tilde {s}_N$, which is given analytically by (2.19).

Finally, the roll waves attain their maximal size when they meet a homoclinic loop, formed by the stable and unstable manifolds of the second fixed point of the dynamical system, $(\tilde {h}_{+},0)$. This second fixed point – a saddle – corresponds to uniform flow conditions with a smaller thickness than the incoming flow. The homoclinic loop is associated with an unstable solitary wave with a considerably larger amplitude than the roll waves; see figures 8(c,d) and 9. When the limit cycle meets this homoclinic loop (which happens along the upper part of the limit cycle, which is associated with the exponentially rising flank of the roll wave), the limit cycle becomes unstable. The orbit that previously traced out the limit cycle now escapes to infinity, as shown in figure 7. This escape demarcates the maximum size the limit cycle (and hence the roll wave) can acquire given a certain Froude number $Fr_0$ of the base flow.

The whole evolution of the roll waves, from their birth until they reach their maximal size, takes place in a surprisingly narrow region in parameter space. This region, between the Hopf bifurcation (H) and the homoclinic event (‘hom’), is shaded in the $(Fr_0, \tilde {c})$ parameter plane of figure 4. The same region is also sketched in the three-dimensional parameter space $(Fr_0, \tilde {c}, \zeta )$ of figure 10.

Throughout our analysis, we adopt the frictional law (2.4) – or in its non-dimensional form, (2.12) – that involves the offset $\varGamma$ of the Froude number. To the best of our knowledge this is the first time that this has been done for roll waves. The offset is found to cause no qualitative changes to the sequence of events from the birth of the roll waves to the stage when they attain their maximal size, although its quantitative effect is considerable. This qualitative robustness of the evolution of granular roll waves under changes in the offset $\varGamma$ has been further examined in § 3.5, where we also found that the evolution is equally robust under variations of the inclination angle $\zeta$ (as long as one stays within the interval $\zeta _1 < \zeta <\zeta _2$).

The aforementioned ‘escape’ of the limit cycle is not the end of the roll waves. It only means that for a given value of $Fr_0$ there is a maximum to the amplitude which the roll waves can attain. We performed a numerical experiment to show that, when the influx of material on the chute is enhanced (i.e. when we increase the value of $Fr_0$), the roll waves on the chute spontaneously adapt their amplitude and wave speed to a correspondingly larger value; see figure 12. Inversely, when we decrease the influx (thereby lowering the value of $Fr_0$), the amplitude of the roll waves self-regulates itself by becoming smaller.

In order to put the results of the present paper to the test, we recommend to use an experimental set-up with $Fr_{cr}$ as small as possible, such that the condition for the appearance of roll waves ($Fr_0>Fr_{cr}$) will be easy to satisfy. This means, since $Fr_{cr} = (2/3) (1-\varGamma )$, that it will be advantageous to choose a material with a large value of $\varGamma$.

Such a set-up was used in the experiments of Rocha et al. (2019), who – among other things – studied the flow of red sand with $\varGamma =0.84$, implying that $Fr_{cr}=0.107$. Indeed, in ‘movie 1’ from the online supplementary material to that paper, snapshots of which are shown on the cover of the journal, a wavy pattern is clearly discernible. Given the fact that the flow in this video has reached a steady state, and that the Froude number certainly exceeds the critical value $0.107$ (since $\beta _{*}=0.11$ for the set-up in question and the authors explicitly state that the flow is fully dynamic, i.e. $Fr_0>\beta _{*}$), we strongly believe that this wavy pattern is a periodic train of roll waves. Therefore, in our view, the set-up used by Rocha et al. (2019) is perfectly suited for testing the results of the present study. For instance, for $Fr_0$ considerably larger than $Fr_{cr}$, our theory predicts that the dimensionless wave speed $\tilde {c}$ will be very close to the Hopf value $\tilde {c}_H = 1+1/Fr_0$ and we think that this might be checked in a relatively straightforward way on this specific chute. That is to say, for any $Fr_0 \gg Fr_{cr}$ we predict that the observed (dimensional) speed of the roll wave train should be close to the value given by (3.3).

As a final note, we want to point out that the methodology of the present paper is not restricted to granular fluids. It can also be applied to ordinary Newtonian fluids, where one may anticipate to find qualitatively similar results. Pioneering work on roll waves in water using the dynamical systems approach has been performed by Needham & Merkin (1984) and Merkin & Needham (1986). One may expand on this by incorporating the realistic diffusive term proposed by Kranenburg (1992) and taking into account also more recent results on the dynamics of roll waves in normal fluids, such as those by Balmforth & Mandre (2004).