Cross-stream oscillations in the granular flow through a vertical channel

Abstract

The gravity flow of a granular material between two vertical walls separated by a width $2W$ is simulated using the discrete element method (DEM). Periodic boundary conditions are applied in the flow (vertical) and the other horizontal directions. The mass flow rate is controlled by specifying the average solids fraction $\bar {\phi }$, the ratio of the volume of the particles to the volume of the channel. A steady fully developed state can be achieved for a narrow range of $\bar {\phi }$, $\bar {\phi }_{max} \geq \bar {\phi } \geq \bar {\phi }_{cr}$, and the material is in free fall for $\bar {\phi } < \bar {\phi }_{min}$. For an intermediate range of $\bar {\phi }$ ($\bar {\phi }_{cr} > \bar {\phi } \geq \bar {\phi }_{min}$), there are oscillations in the horizontal coordinate of the centre of mass, velocity components and stress. As $\bar {\phi }$ decreases in the range $\bar {\phi }_{cr} > \bar {\phi } \geq \bar {\phi }_{min}$, the amplitude of the oscillations increases proportional to $\sqrt {\bar {\phi }_{cr} - \bar {\phi }}$ and the frequency appears to tend to a non-zero value as $\bar {\phi } \rightarrow \bar {\phi }_{cr}$, indicating a supercritical Hopf bifurcation. The relation between the dominant frequency and the higher harmonics of the position, velocity and stress fluctuations are explained using the momentum balance. It is found that dissipation in the inter-particle and particle–wall contacts is critical for the presence of an oscillatory state.

1. Introduction

The gravity flow of granular materials through standpipes, hoppers and bunkers is encountered in numerous unit operations such as reactors, mixers, separators and heat exchangers. Examining the pressure exerted by the flowing material on the walls has practical importance for developing a reliable design. Earlier studies reported fluctuations in the wall stresses (Rao & Venkateswarlu 1975; Tüzün & Nedderman 1985; Roberts & Wensrich 2002; Ramírez, Nielsen & Ayuga 2010; Wang et al. 2020). The ‘stick-slip’ motion of particles past the wall and ‘silo-quaking’ are some common sources of such fluctuations (Tüzün & Nedderman 1985; Tejchman & Gudehus 1993; Dhoriyani et al. 2006). Apart from this, instability in the flow can lead to a rich variety of patterns that affect the wall stresses significantly. Using the discrete element method (DEM) based simulations, the current work explores an instability-driven oscillatory flow and oscillatory wall stresses in gravity flow through a vertical channel.

Instability in granular flows is as fascinating as that in the hydrodynamic stability of simple fluids. Inelastic inter-particle collisions in these dry frictional materials, which are the major source of energy dissipation, drive an instability resulting in formation of anisotropic structures, clusters of particles, vortices and density waves (Hopkins & Louge 1991; Goldhirsch & Zanetti 1993; Raafat, Hulin & Herrmann 1996; Kudrolli, Wolpert & Gollub 1997; Kumaran 1998; Hua & Wang 1999; Luding & Herrmann 1999; Alam & Hrenya 2001; Liss, Conway & Glasser 2002; Conway & Glasser 2004; Mitrano et al. 2011; Fullmer & Hrenya 2017). Their presence in the flow has a crucial impact on stress states, which is analogous to turbulent stress in simple fluids, and mixing and segregation (Liss & Glasser 2001; Conway, Liu & Glasser 2006). The kinetic theory based models for inelastic hard-particles predict a strong connection between inelastic dissipation and instability (Babic 1993; Brey et al. 1998; Garzó 2005), which is in agreement with hard-particle simulations (Hopkins & Louge 1991; Goldhirsch & Zanetti 1993; Mitrano et al. 2011, 2013).

In a molecular gas, their fluctuating motions are due to their mean motions and true temperature, whereas the fluctuating motions in athermal granular materials, a measure of the ‘granular temperature’, are only because of the mean motions of the particles. In granular flow, the fluctuating velocities and the mean free path of the particles are strongly related to the nature of inter-particle collisions. If the collisions are significantly dissipative, there is a tendency to form clusters of particles where the frequency of collision increases resulting in a further decrease in the mean free path – this phenomenon is known as ‘clustering instability’. Luding & Herrmann (1999) have shown that the rate of cooling is higher at the beginning in freely cooling systems and subsequently reduces resulting in an inelastic collapse. If there is an energy source and the flow sustains, then rich patterns form.

Raafat et al. (1996) performed experiments with a narrow long vertical tube of circular cross-section. They observed three regimes: (i) at high mass flow rates, the material is very dilute and in free fall; (ii) the dense material flows down slowly and uniformly through the length of the pipe at low mass flow rates and (iii) the wave regime at intermediate mass flow rates where the stream breaks down into dense and dilute clusters in the flow direction, characterised as ‘density waves’. These waves appear after traversing a certain distance from the inlet of the pipe. Aider et al. (1999) observed an oscillatory velocity in the downward (flow) direction in the wave regime. They mentioned that the material is in free fall up to a distance from the inlet before density waves appear downstream inside the pipe, and the oscillatory velocity may resemble stick-slip motion.

To capture density waves, Hua & Wang (1999) measured the temporal variation of the density at a circular cross-section using electrical capacitance tomography, which is more accurate than the methods used earlier. Their results showed an oscillatory trend of the density with time at the central region. The earlier work reported the appearance of density waves in the flow direction whereas Hua & Wang (1999) observed these in a cross-section at a distance from the inlet. They proposed that the particle collisions are the only mechanism for the formation of density waves. The particles bounce back towards the centre after colliding with the wall forming a dense zone of high solids fraction $\phi$. Further collisions allow the particles to spread across the cross-section and hence a reduction in $\phi$ at the centre. The radial variation of the solids fraction seems asymmetric in some of their results. This may indicate that there is a possibility of cross-stream velocity and stresses oscillating with time, which has not been reported there. It must be noted that the studies mentioned above could not minimise the effect of air drag on the particle phase in their experiments. Further, the flow is quite dilute, with the maximum value of $\phi$ being only approximately 0.16.

Wang, Jackson & Sundaresan (1997) studied the stability of fully developed gravity flow bounded between two vertical walls at a distance $2W$ using the kinetic theory based model of Lun et al. (1984) and the boundary conditions of Johnson & Jackson (1987). The effect of air drag was neglected in the model. The stability criteria of the particle phase depends solely on the average solids fraction $\bar {\phi }$ (ratio of the volume of the particles to the volume of the channel) and the width of the channel $2W$. They reported that there are three dominant growing modes which are travelling waves. At an unstable base state, these waves propagate in the flow direction generating three types of density waves.

Later, Wang & Tong (2001) solved the unsteady equations based on kinetic theory numerically for a two-dimensional flow and analysed different types of density waves. Depending on the ratio of the height $H^*$ of the channel to its width $2W$, the distribution of $\phi$ with respect to the central plane is either symmetric or anti-symmetric. A cross-stream flow is evident, and the formation of dense and dilute clusters in the flow direction is shown by their results. These instabilities are shown to be of inertial nature and dependent on inelastic particle–particle collisions. However, they did not integrate the equations for a long period of time, and the values of $\bar {\phi } = 0.15$ and 0.4 used are very small. Intuitively, one can expect a very low resistance from the walls on the material for small values of $\bar {\phi }$, resulting in an accelerated flow when air drag is neglected. Surprisingly, this is not seen in their results. They stated that the inertial effects and significant dissipation due to inelastic collisions cause the instability; however, the dominant mechanism for instability is not clearly understood from their work.

Liss et al. (2002) performed simulations with rigid disks flowing inside a two-dimensional vertical channel and applied periodic boundary conditions in the flow direction. The degree of inelastic collision and roughness of particles are determined by the coefficients of normal and tangential restitution ($e_n$, $e_t$). For a small value of $(H^*/2W)$, the cross-stream flow and density waves are absent in their simulation. The reason is the following. The spatial variations of length scale higher than the dimension of a periodic domain $H^*$ are suppressed because of the use of periodic boundaries, and thus cannot be captured (Babic 1993; Allen & Tildesley 2017). The predictions of simulations by Liss et al. (2002) reflect a stronger dependence of the dimension of the periodic domain on the manifestation of different types of density waves, qualitatively in agreement with Wang & Tong (2001).

The work of Wang & Tong (2001) is based on the integration of the equations of the kinetic theory for smooth particles in a two-dimensional flow. The only source of dissipation is the normal coefficient of restitution $e_n$ for inter-particle collisions. Subsequently, Liss et al. (2002) used event-driven hard particle DEM simulations of disks, where the dissipation is modelled by the coefficients of restitution $e_n$ and $e_t$. The parameter $e_t$ accounts for the roughness of the particles. The present work is based on soft particle DEM simulations where dissipation can occur because of viscous damping in the normal and tangential directions, and the surface friction is also included. It is of interest how these different factors influence the dissipation and the instability.

The granular flow through a vertical channel was examined by Debnath, Kumaran & Rao (2019, 2022a); Debnath, Rao & Kumaran (2022b); they found a narrow range of the average solids fraction ($\bar {\phi }_{max}$, $\bar {\phi }_{cr}$) for which a steady flow can be obtained. The flow ceases when $\bar {\phi }$ is above $\bar {\phi }_{max}$, and there is an oscillatory flow followed by an accelerated flow when $\bar {\phi }$ is below $\bar {\phi }_{cr}$. Four different shear regimes were identified for a steady fully developed flow: a central plug region where there is no shear; a dense shearing zone adjoining the plug region where the solids fraction $\phi$ exceeds 0.587, which is the solids fraction for the arrested dynamics (Kumaran 2008); a loose shearing zone where $\phi$ is less the 0.587 and a wall shearing zone extending over a distance of 2–3 particle diameters from the wall where the shear rate is large. It was shown that the rheology is well explained using the hard-particle model in the loose shearing zone, but the particle stiffness affects the rheology in the dense shearing zone and the plug zone. An extended kinetic theory has been proposed by Berzi, Jenkins & Richard (2020) where the effect of the particle stiffness is incorporated, and applied for the present flow configuration by Debnath et al. (2022a) and Islam, Jenkins & Das (2022).

The present work demonstrates the presence of an oscillatory state when $\bar {\phi }$ is lower than that for the steady flow in a vertical channel. These oscillations are self-induced and the characteristics of these oscillations are examined. The flow configuration and simulation method are discussed in § 2, the analysis of the flow characteristics is the subject of §§ 3 and 4, and the important results are discussed in § 5.

2. Flow configuration and simulation method

The gravity flow of particles through a vertical channel of rectangular cross-section is simulated. The channel is of finite width $2W$ in the $x$ (cross-stream) direction, and periodic boundary conditions are used in the other two directions, the horizontal $z$ (spanwise) direction of width $B$ and the vertical $y$ (flow) direction of length $H$. A schematic of the flow configuration is shown in figure 1.

Figure 1. Flow in a vertical channel. The dimensions of its rectangular cross-section in the $z$ and $x$ directions are $B\times 2W$. (a) A static bed is filled up to height $H - \Delta H$ with a bottom. (b) In the flowing bed, the bottom is removed and periodic boundary conditions are applied in the $z$ and $y$ directions. The height of the flowing bed is $H$.

The initialisation protocol is the following. A flat bottom is placed below the channel initially and the particles are poured into the channel by raining them uniformly and allowing them to settle under gravity. The channel is filled up to a height $H - \Delta H$, and there is a vacant head space of height $\Delta H$ above the particles (figure 1a). The average solids fraction $\bar {\phi }$, the ratio of the volume of the particles to the volume of the channel, of the static bed is close to $\phi _{drp} = 0.64$, the solids fraction corresponding to the dense random packing. During the onset of flow, the flat bottom is removed and periodic boundary conditions are applied (figure 1b). In the flowing state, the value of $\bar {\phi }$ is set by changing $\Delta H$. The spherical particles are polydisperse to avoid crystallisation, with a number fraction 0.3 for 0.9 $d_p$, 0.4 for $d_p$ and 0.3 for 1.1 $d_p$. In the simulations, $(H-\Delta H)/d_p$ and $B/d_p$ are set to 60 and 40, respectively, and $2W/d_p$ is varied in the range 40–120.

The DEM is used to integrate the motion of particles and the flow is simulated using LAMMPS open source software (Plimpton 1995). The event driven simulation method for inelastic hard spheres is computationally more expensive for the present system. This is because the solids fraction $\phi$ is close to $\phi _{drp}$ in the centre of the channel and the time interval between successive collisions is very small. It results in longer simulation times to capture collisions of all particles. Apart from that, the former method cannot be used if there are enduring frictional contacts. Contact dynamics allows instantaneous collisions, but multi-body collisions are permitted within one time step. As we are interested in tracking particle overlaps and frictional contacts, the DEM has been used in the present work.

The linear spring-dashpot model is used, with spring constants ($k_n$, $k_t$) in the normal and tangential directions to the surfaces at contact, damping constants ($\xi _n$, $\xi _t$), and the coefficient of inter-particle friction $\mu _p$. Unless otherwise specified, the results correspond to $k_n = 10^6 \rho _p g d_p^2$, $k_t = 2/7 k_n$, $\xi _n = 180 \sqrt {g/d_p}$, $\xi _t = 1/2 \xi _n$ and $\mu _p = 0.5$, where $\rho _p$ is the intrinsic density of a particle and $g$ is the acceleration due to gravity. The values of the parameters for the particle–wall interactions are same as for the inter-particle interactions. Details of the DEM and methods to calculate properties are described by Debnath et al. (2022a,b).

The variation of the flow characteristics with $k_n$ in the plug region was examined by Debnath et al. (2022b) (see figure 10 and the accompanying discussion). It was shown that the particle stiffness has a significant effect on the temperature profile in the plug region. However, there was no change in the profiles of the solids fraction, velocity and shear rate for $k_n/(\rho _p g d_p^2) > 10^5$. As the interesting oscillatory dynamics is observed in the shearing zone close to the wall, the scaled spring stiffness has been chosen to $10^6$ in the present study. We have verified that an increase in the scaled spring stiffness to $10^7$ has no effect on the dynamics to within the simulation resolution.

All results are shown for the scaled quantities, where the length, time, velocity, granular temperature, force and stress are scaled by $d_p$, $\sqrt {d_p/g}$, $\sqrt {g d_p}$, $g d_p$, $\rho _p g d_p^3$ and $\rho _p g d_p$, respectively. Without loss of generality, $d_p$, $\rho _p$ and $g$ are set equal to 1.

3. Results

The flow dynamics is characterised by the variations in the coordinates and the velocities of the centre of mass, denoted by the subscript ${}_{cm}$,

(3.1)\begin{equation} \varpi_{cm} = \frac{\displaystyle\sum_{i = 1}^{N} m_i \varpi_i}{\displaystyle\sum_{i=1}^{N} m_i}, \end{equation}

where $N$ is the total number of particles in the simulation box and $m_i$ is the mass of particle $i$. The variation of the $x$ coordinate of the centre of mass $x_{cm}$ and the flow velocity of the centre of mass $-u_{y|cm}$ with time $t$ are shown in figure 2. The following states are observed.

1. (i) There is no flow for $\bar {\phi } > 0.62$. There is a very small initial velocity, but it decreases quickly to zero, and the material between the two walls reaches a jammed state.

2. (ii) For $\bar {\phi } = 0.60 < 0.62$, the material is in free fall initially, and $-u_{y|cm}$ increases linearly with time as expected. There are oscillations in $x_{cm}$ of very small amplitude till $t \equiv t^\ast = 100$. After $t^\ast = 100$, the material experiences equal resistance from both the walls where the walls bear the weight of the material, and the flow reaches a steady fully developed state.

3. (iii) For $\bar {\phi } = 0.57$, the initial transients are similar to those of $\bar {\phi } = 0.60$. The amplitude and frequency of $x_{cm}$ increase with time till $t^\ast \approx 1050$. There is a sharp transition in the profile of $-u_{y|cm}$ after $t^\ast \approx 1050$ where there is traction with the walls; the acceleration decreases nearly to zero. The flow is not steady, and $x_{cm}$ exhibits oscillations with frequency much higher than that of the initial transients. The amplitude and frequency of the oscillations attain a constant value. This is termed as an ‘oscillatory flow’. In this state, the upward forces exerted by the two walls on the material are not equal as $x_{cm}$ oscillates. However, the time average of the upward forces exerted by the walls are equal, and their sum is equal to the weight of the material.

4. (iv) If $\bar {\phi }$ is further reduced to 0.552, the material is in a state of vertical free fall till $t = 1500$, and this state is termed as an ‘accelerated flow’. The simulations have been performed till $t=4000$, but the flow continues to accelerate and there is no sign of a steady state. If a mass accelerates with gravitational acceleration $g\, (= 1)$ and zero initial velocity, its velocity after traversing a distance, say $H$, would be $\sqrt {2 g H}$. Hence, figure 2 implies that the flow reaches an oscillatory state after traversing a distance ${\sim } 7 \times 10^3 \, H$ for $H = 66$ and $\bar {\phi } = 0.57$. However, for $H = 68.15$ and $\bar {\phi } = 0.552$, the free fall persists even after traversing a distance of ${\sim } 110 \times 10^3 \, H$. Due to constraints on the computational time, the motion of the particles could not be simulated beyond $t = 4000$.

Figure 2. Variation of the $x$ coordinate of the centre of mass $x_{cm}$ (solid curves, left ordinate) and the flow velocity of the centre of mass $-u_{y|cm}$ (dashed (red) curves, right ordinate) for a channel of width $2W = 100$; $\diamond$, $\bar {\phi } = 0.625$ (no flow); $\triangle$, $\bar {\phi } = 0.60$ (steady fully developed flow); $\circ$, $\bar {\phi } = 0.57$ (oscillatory flow); $\square$, $\bar {\phi } = 0.552$ (accelerated flow). The time series data are shown for $f_{sf} = 1/(500\, \Delta t)$, where $f_{sf}$ is the sampling frequency (see Appendix A) and $\Delta t$ is the simulation time step.

When $\bar {\phi }$ is low, the resistance to flow from both the walls is not experienced significantly by the material in the transition period ($t < t^\ast$). The material drifts very slowly (in the cross-stream direction) from one wall to the other while flowing down. At some instances while drifting, the material is completely detached from both walls. This causes the free fall in the time period $t < t^\ast$. When both the walls exert equal resistance, the flow reaches a steady state, and if the effect is unequal, that drives the flow to a new oscillatory state.

The flow regime map for the width of the channel in the range $2W = 40\unicode{x2013}120$ and flat frictional walls is shown in figure 3. (A preliminary study can be found in the paper by Debnath et al. 2019.) As $2W$ increases, the range of $\bar {\phi }$ for the different regimes decreases. The boundary between no flow and a steady flow at $\bar {\phi } \equiv \bar {\phi }_{max} = 0.62$ is independent of $2W$. The boundaries between a steady and an oscillatory flow, and between an oscillatory and an accelerated flow depend on $2W$. The steady fully developed flow is observed in a very narrow range $\bar {\phi }_{cr} \leq \bar {\phi } \leq \bar {\phi }_{max}$. The features of the steady fully developed flow and the predictions of different classes of continuum models were discussed by Debnath et al. (2022a,b). An oscillatory flow is observed when $\bar {\phi }$ is decreased below $\bar {\phi }_{cr}$. An accelerated flow is observed when $\bar {\phi }$ is further decreased below a value $\bar {\phi }_{min}$, where the material is in free fall. The oscillatory regime in the range $\bar {\phi }_{min} \leq \bar {\phi } < \bar {\phi }_{cr}$ is the focus of the present work.

Figure 3. Flow regime map for a vertical channel of rectangular cross-section with flat frictional walls fixed at a distance $2W$.

3.1. Oscillatory regime

The centre of mass coordinate and velocities are shown with a greater magnification in figure 4. The $x$ component of the velocity of the centre of mass $u_{x|cm}$ oscillates with the same frequency as $x_{cm}$ and a phase difference of ${\rm \pi} /2$, because the velocity is a maximum when $x_{cm}$ passes through the centre of the channel. Here $-u_{y|cm}$ decreases when $|x_{cm}|$ is a maximum in either direction due to a larger traction with the wall, and increases when $x_{cm}$ passes through the centre of the channel. The frequency of $-u_{y|cm}$ is two times that of $x_{cm}$, and a possible reason will be provided shortly using the momentum balance equations.

Figure 4. In the oscillatory regime, the expanded view of $x_{cm}$, $u_{x|cm}$ and $-u_{y|cm}$ in a very small time interval for $2W = 100$ and $\bar {\phi } = 0.57$. The black and red curves correspond to the left and right ordinates, respectively; $x_{cm}$, black solid curve; $u_{x|cm}$, black dotted curve.

In the oscillatory state, a variable $\varpi$ is defined as the sum of its time-averaged ‘base’ state $\varpi ^b$, which is time-independent, and the oscillations about the base state $\varpi ^\prime$,

(3.2)\begin{equation} \varpi = \varpi^b + \varpi^\prime. \end{equation}

The time-averaging is done over $5\times 10^5$ time steps. Here, $x_{cm}^{b}$ and $u_{x|cm}^{b}$ are both zero, and $u_{y|cm}^{b} \neq 0$ is the time-averaged velocity in the vertical direction. In the time series profiles, the data are collected at a sampling frequency $f_{sf}$. The effect of the sampling frequency on the results is discussed in Appendix A.

The effect of the ratio of the height $H$ to width $2W$ has been examined in Appendix B to investigate whether the time series of $x_{cm}$ varies with the distance along the flow. The value of $H/2W$ is increased up to a maximum of $2.5$ for $2W = 80$ and $H = 200$. Even for this large height, no variation of $x_{cm}$ with time has been found at different locations in the flow direction, as shown in Appendix B.

3.2. Solids fraction

The variation of the solids fraction $\phi$ and the wall stresses are shown in figure 5. The time series of $\phi$ at the centre of the channel is nearly steady, and it has a high value of approximately 0.63–0.64 (the dotted curve in figure 5a). In contrast, $\phi$ at the the walls (the solid and dashed curves in figure 5a) shows oscillations in the range $\phi = 0.08\unicode{x2013}0.3$ with an average $\sim$0.16. (Note that the properties at the walls are obtained by averaging in a bin adjacent to the walls of thickness 1.5 $d_p$.) The wall shear and normal stresses (figure 5b) exhibit sharp maxima when $\phi$ passes through a maximum, and are close to zero when $\phi$ passes through a minimum. The maxima in the stresses are sharp, but the minima are shallow and close to zero. Therefore, the maximum is more than two times the average values shown by the dotted horizontal lines in figure 5(b).

Figure 5. A typical time series of (a) the solids fraction $\phi$ at the walls and (b) wall stresses $\sigma _{xx}$ (black, left ordinate) and $\sigma _{xy}$ (red, right ordinate). In panels (a) and (b), the solid curves represent the left wall at $x' \equiv x + W = 0$ and the dashed curves for the right wall at $x' = 2W$; the dotted curve in panel (a) is at the mid-plane $x' = W$. The solid horizontal line in panel (a) represents the time-averaged value of $\phi$; the dotted (black and red) horizontal lines in panel (b) represent the time-averaged values of $\sigma _{xx}$ and $|\sigma _{xy}|$. (c) Variation of $\phi$ with the distance from the left wall $x'$. In panel (c), $\diamond$ and $\square$ represent profiles for time $t_1 = 6001.2$ and $t_2 = 6003.6$ when the material is fully compressed to the left and right walls, respectively, $\circ$ represents the time-averaged base state $\phi ^b$, and the inset is a magnification of the profiles near the right wall. The parameter values are $2W = 100$ and $\bar {\phi } = 0.57$. The time series data are shown for $f_{sf} = 1/(500 \Delta t)$, where $f_{sf}$ is the sampling frequency (see Appendix A) and $\Delta t$ is the simulation time step.

The solids fraction $\phi$ and the normal and shear stresses at the two walls are anti-correlated (figure 5a,b), i.e. when $\phi$ is a minimum at one wall, it is a maximum at the other wall, implying that the material cyclically detaches from one wall and accumulates near the other wall. The same trend is reflected by the normal and shear stresses exerted by the material on the walls which is expected. The stresses are close to zero when the material detaches from a wall. There is a large traction at one wall for a short period of time when $\phi$ at that wall is close to 0.3; at the same instant, the other wall experiences nearly no traction when $\phi$ is approximately 0.08.

The profiles of $\phi$ corresponding to the time instants when the wall stresses are maximum and minimum are shown in figure 5(c). The time-averaged profile $\phi ^b$ is symmetric about the centre of the channel, while the other two are asymmetric and (almost) mirror images about the centreline. The distance between the time-averaged value of $\phi$ at the wall and the corresponding values of the profiles of $\phi$ for $t_1$ and $t_2$ (see caption of figure 5) are 4.75 and 3.25, respectively. The mean distance is 4.0, which is approximately twice $A_{x_{cm}} = 2.14 \pm 0.11$, where $A_{x_{cm}}$ is the amplitude of the dominant frequency of $x_{cm}$. When the material is fully compressed towards one wall, the thickness of the layer is found to be approximately $2(W-A_{x_{cm}})$, in accordance with $x_{cm}$ which traverses a maximum distance $2 A_{x_{cm}}$ in half a period of oscillation. The corresponding value of $\bar {\phi }$ in this layer of dimensions $H \times B \times 2(W-A_{x_{cm}})$ is $\approx \bar {\phi }_{cr}$. Hence, the effective width of the flowing layer turns out to be $2(W - A_{x_{cm}})$, two times the difference between the half-width of the channel and the amplitude of oscillations.

3.3. Flow velocity

As mentioned earlier, the time-averaged flow velocity $u^b_{y|cm} \neq 0$, and it exhibits a monotonic variation with $2W$ and $\bar {\phi }$ (figure 6). Here, $u^b_{y|cm}$ is shown both for the steady fully developed state ($\bar {\phi } \geq \bar {\phi }_{cr}$) and the oscillatory state ($\bar {\phi } < \bar {\phi }_{cr}$) – the two are separated by the pink filled circles in figure 6(a). The magnitude of $u^b_{y|cm}$ increases as $\bar {\phi }$ decreases and $2W$ increases. Figure 6(a) shows a clear discontinuity in the slope at the transition between the steady and the oscillatory states. This suggests a continuous transition between the two states. The scaled velocity, $-u^b_{y|cm}/\sqrt {2(W-A_{x_{cm}})}$, is plotted as a function of $(\bar {\phi }_{cr}-\bar {\phi })$ in figure 6(b), where $A_{x_{cm}}$ is the amplitude of the oscillation in $x_{cm}$. The data fall on a single curve for $(\bar {\phi }_{cr}-\bar {\phi }) \rightarrow 0^+$ in the oscillatory state, implying that the Froude number based on the average velocity and the effective width $2(W - A_{x_{cm}})$ is independent of the width $2W$. A similar collapse is not observed for the steady fully developed state for $\bar {\phi } \geq \bar {\phi }_{cr}$ as $\bar {\phi }_{cr}$ is a function of $2W$ (see figure 3). However, our earlier work (Debnath et al. 2022b) has shown that the scaling is valid if the scaled velocity $-u^b_{y|cm}/\sqrt {2W}$ is plotted as a function of $\bar {\phi }$ for the steady fully developed flow. This suggests a Beverloo-type correlation (Beverloo, Leniger & Van de Velde 1961) for the oscillating mass of particles with a characteristic length scale $(W-A_{x_{cm}})$. There is a departure from this scaling near the transition to the accelerated regime where $\bar {\phi }$ is close to $\bar {\phi }_{min}$.

Figure 6. Variation of (a) $u^b_{y|cm}$ with $\bar {\phi }$ and (b) its scaled values with $\bar {\phi }_{cr} - \bar {\phi }$. The magenta curve in panel (a) marks the boundary between the oscillatory and the steady fully developed states. In panel (b), the oscillatory and steady states correspond to the right and left of the vertical dashed line, respectively. Parameter values: $\diamond$, $2W = 60$; $\square$, $2W = 80$; $\circ$, $2W = 100$; $\triangle$, $2W = 120$.

3.4. Stress

The time-averaged stresses are shown in figure 7. From the momentum balance, the normal stress $\sigma ^b_{xx}$ is a constant across the channel, i.e. $\sigma _{xx} \equiv \sigma ^b_{xx} = {\rm constant}$, and the magnitude of $\sigma ^b_{xy}$ increases from the centre to the wall for the steady fully developed flow (Debnath et al. 2022a,b). The maximum value of $|\sigma ^b_{xy}|$ at the wall is equal to the weight of the material in one half of the channel, i.e. $|\sigma ^b_{xy}| (x = W) = W \bar {\phi }$ for both the steady and the oscillatory states. The stresses are continuous at the transition between the steady and the oscillatory states. There is a distinct change in the shape of the normal stress $\sigma ^b_{xx}$ at the transition; – $\sigma ^b_{xx}$ decreases as $\bar {\phi }$ is decreased in the steady fully developed flow. However, the normal stress shows only a small variation with $\bar {\phi }$ in the oscillatory state (figure 7a). As the variation of $\bar {\phi }$ is in a very small range (<6 %), the variation of $\sigma ^b_{xy}$ at the wall with $\bar {\phi }$ is small (figure 7b). The stresses vary with $2W$. In the oscillatory state, the ratio of the stresses and the width is independent of $2W$ and $(\bar {\phi }_{cr}-\bar {\phi })$, as shown by the red dashed curves in figure 7. Because $\bar {\phi }_{cr}$ varies with $2W$, the scaling is not suitable in the steady fully developed state for $(\bar {\phi }_{cr} - \bar {\phi }) \leq 0$. However, our earlier work (Debnath et al. 2022b) has shown that the scaling is suitable if it is a function of $\bar {\phi }$ for the steady fully developed flow.

Figure 7. Variation of the time-averaged stresses (a) $\sigma ^b_{xx}$ and (b) $\sigma ^b_{xy}$, and their scaled values with $\bar {\phi }_{cr} - \bar {\phi }$ on the left wall $x' = 0$; solid and dashed curves correspond to the left and right ordinates, respectively. The oscillatory and steady states correspond to the right and left of the vertical dashed line, respectively. Parameter values: $\diamond$, $2W = 60$; $\square$, $2W = 80$; $\circ$, $2W = 100$; $\triangle$, $2W = 120$.

It is interesting to note that $\sigma ^b_{xx}$ is approximately independent of $(\bar {\phi }_{cr} - \bar {\phi }) > 0$ (oscillatory state). A plausible reason is the following. The stress arises from inter-particle contacts and fluctuating velocities (see Debnath et al. 2022b, (2.4)). In the steady fully developed flow, $u_x$ is zero and the square of the velocity fluctuation $u'_x$ is negligible. Hence, the major contribution to $\sigma ^b_{xx}$ is due to the contact forces which decrease with $\bar {\phi }$ (intuitively) and perhaps saturate as $\bar {\phi } \rightarrow \bar {\phi }_{cr}^+$. The data for $2W = 60$ and 80 show this trend, whereas there is no clear trend for larger widths. Now, in the oscillatory state, the normal stress is due to the contact forces, the oscillatory velocity of the centre of mass in the $x$ direction and the fluctuating velocity of the particles about the latter. It is speculated that the oscillations in $u_x$ only affect the perturbed normal stress $\sigma ^\prime _{xx}$ and not the time-averaged base state $\sigma ^b_{xx}$. Because of this, $\sigma ^b_{xx}$ is independent of $\bar {\phi }$ in the oscillatory regime for a given width. As $2W$ increases, $\sigma ^b_{xx}$ increases for a given value of $\bar {\phi }_{cr} - \bar {\phi }$ and the former appears to be proportional to $2W$, consistent with the Janssen solution at large depths (Rao & Nott 2008).

As both the shear and normal stresses are proportional to $2W$, the wall friction coefficient $\mu ^b_w = |\sigma ^b_{xy}|(x = W)/\sigma ^b_{xx}$, appears to be a constant for the time-averaged base states in the oscillatory flow, which is independent of $\bar {\phi }$ and $2W$. In contrast, the coefficient of friction strongly depends on $\bar {\phi }$ and does not depend on $2W$ in the steady fully developed state (Debnath et al. 2022b). The wall stresses, $\sigma _{xx}$ and $\sigma _{xy}$, and their perturbed states, $\sigma _{xx}^\prime$ and $\sigma _{xy}^\prime$, are plotted against each other at many time instants (figure 8). The ratio, which is the coefficient of friction, is close to 0.45 for both $|\sigma _{xy}|/\sigma _{xx}$ and $|\sigma ^\prime _{xy}|/\sigma ^\prime _{xx}$. Though it appears that the friction law can be applied to the oscillatory component of the stress as well, a common friction coefficient may not be applicable to model the wall stresses in the oscillatory state.

Figure 8. Variations of $|\sigma _{xy}|$ with $\sigma _{xx}$ (black markers, left ordinate) and $|\sigma ^\prime _{xy}|$ with $\sigma ^\prime _{xx}$ (red markers, right ordinate) at the left ($\circ$) and the right ($\square$) walls for (a) $2W = 100$ and (b) $2W = 80$. The slope of the blue dashed lines in panels (a) and (b) is 0.45.

3.5. Effect of dissipation on oscillations

The effect of the dissipation on the oscillatory flow is examined. There are two sources of dissipation – frictional dissipation and damping dissipation. Two important observations are the following:

1. (i) It has been observed that if the coefficient of inter-particle friction is reduced from $\mu _p = 0.5$ to 0.25 (keeping inter-particle damping constant $\xi _n$ fixed), the material accelerates even if $\bar {\phi }$ is in the range $\bar {\phi }_{max} \geq \bar {\phi } \geq \bar {\phi }_{cr}$, and the flow does not reach a steady fully developed state. Hence, a decrease in the frictional dissipation has a destabilising influence on the flow and the choice of $\mu _p$ appears to be significant.

2. (ii) The effect of the inter-particle damping constants $\xi _{n,p-p}$ and the particle–wall damping constants $\xi _{n,p-w}$ are examined separately. The ratio of the tangential to the normal damping constants are set the same in all cases. For the cases studied here, it is found that the oscillations disappear and the flow is steady if the inter-particle damping is absent, i.e. $\xi _{n,p-p} = 0$ (keeping $\mu _p = 0.5$ fixed) in the range $\bar {\phi }_{cr} > \bar {\phi } \geq \bar {\phi }_{min}$. Hence, a reduction in the damping dissipation has a stabilising influence on the flow, in contrast to the role of the frictional dissipation. However, oscillatory states are present even if the particle–wall damping coefficient is set to zero, $\xi _{n,p-w} = 0$. Thus, damping in particle–wall collisions seems to have no effect on the oscillatory states. A detailed analysis of the friction and damping coefficients is provided in Appendix C.

4. Spectral analysis of the oscillations

4.1. $x_{cm}$ and $u^\prime _{y|cm}$

The amplitude–frequency spectra of $x_{cm}$ and $u^\prime _{y|cm}$ are shown in figure 9 for $2W=100$ and different values of $\bar {\phi }$. For $\bar {\phi } = 0.61$ (red curve), the frequency spectrum is noisy and exhibits no clear maxima. As $\bar {\phi }$ is decreased, one clear maximum appears for $\bar {\phi } = 0.59$ (blue), but the amplitude of oscillations is found to be approximately 5–7 times smaller than a particle diameter. As $\bar {\phi }$ is further decreased to $0.58$ and $0.57$ (black and magenta), there are multiple maxima, with the largest maximum (higher than a particle diameter) denoted by $A_{x_{cm}}$ and $A_{u^\prime _{y|cm}}$ at the dominant frequencies denoted by $f_{x_{cm}}$ and $f_{u^\prime _{y|cm}}$, respectively, and smaller maxima at the higher harmonics. The amplitudes of the higher harmonics of $x_{cm}$ and $u^\prime _{y|cm}$ are orders of magnitude smaller than that at the dominant one. For definiteness, $\bar {\phi }_{cr}$ (see figure 3) is chosen such that the largest amplitude of the dominant peak in the amplitude–frequency spectra of $x_{cm}$ is less than half a particle diameter. (Note that if the former criterion to obtain $\bar {\phi }_{cr}$ is slightly changed, i.e. the value is reduced from 0.5 to 0.2 particle diameter (60 % reduction), the change in $\bar {\phi }_{cr}$ is observed to be in the range 0.5–0.8 %).

Figure 9. Amplitude–frequency spectra for (a) $x_{cm}$ and (b) $u^\prime _{y|cm}$. Parameter values: $2W = 100$; red, $\bar {\phi } = 0.61$; blue ($\triangle$), $\bar {\phi } = 0.59$; black ($\square$), $\bar {\phi } = 0.58$; magenta ($\circ$), $\bar {\phi } = 0.57$; $f_{sf} = 1/(500 \Delta t)$.

For example, the dominant frequency of $x_{cm}$ is $f_{x_{cm}} = 0.18$ for $\bar {\phi } = 0.57$ and $2W = 100$, which corresponds to 18 Hz for $d_p = 1$ mm and $g = 10$ m s$^{-2}$, and the frequencies of the higher harmonics are $3, 5, 7, \ldots$ times the dominant frequency $f_{x_{cm}}$. As the variation of $x_{cm}$ with time is symmetric about its mean value (figure 4), the amplitude–frequency spectra of $x_{cm}$ contains maxima only at odd multiples of the dominant frequency $f_{x_{cm}}$. The dominant frequency of $f_{u^\prime _{y|cm}}$ is $2 f_{x_{cm}}$. This is because the vertical velocity has two minima when $x_{cm}$ passes through one maximum and one minimum, and the wall shear stress on either wall is the largest. The higher harmonics of $u^\prime _{y|cm}$ are $2, 3, \ldots$ times the dominant frequency $f_{u^\prime _{y|cm}}$, as the time series of $u^\prime _{y|cm}$ is not symmetric about its mean value.

4.2. Wall stress

The spectra for the perturbed normal and shear stresses, $\sigma ^\prime _{xx}$ and $\sigma ^\prime _{xy}$, at the walls are shown in figure 10, and the corresponding time series profiles are shown in the insets. The dominant frequency of the largest maximum in the spectrum of the normal and shear stress at both the walls is the same as that for $x_{cm}$, and other maxima are at $2, 3, 4, \ldots$ times the dominant frequency. The relation between the stress and the velocity spectra can be understood using the momentum balances integrated over the width of the channel. The time-averaged velocity, the time-averaged body force and the time-averaged stresses exerted by the walls have been subtracted to obtain the equations for the perturbations to the acceleration of the centre of mass.

Figure 10. Amplitude–frequency spectra for (a) $\sigma ^\prime _{xx}$ and (b) $\sigma ^\prime _{xy}$ at the walls. Parameter values: $2W = 100$ and $\bar {\phi } = 0.57$; $f_{sf} = 1/(5000 \Delta t)$; black curve ($\square$), left wall ($x = -W$); green curve ($\circ$), right wall ($x = W$); red curve ($\triangle$), the terms on the right-hand side in (4.1) and (4.2). The insets in panels (a) and (b) show the corresponding time series of the perturbed wall stresses.

1. (i) The momentum balance in the cross-stream $x$ direction (for the perturbed state) is

   (4.1)\begin{equation} 2W \bar{\phi} \frac{\text{d} u_{x|cm}}{\text{d} t} \equiv 2W \bar{\phi} \frac{\text{d}^2 x_{cm}}{\text{d} t^2} ={-} (\sigma^\prime_{xx}|_{x=W} - \sigma^\prime_{xx}|_{x={-}W}). \end{equation}
   The left-hand side of the above equations is the scaled mass times acceleration in the $x$ direction per unit area in the $x$–$z$ plane, and the right-hand side is the stresses (force per unit area) exerted by the walls. On the right-hand side of (4.1) is the difference in the normal stress at the left and the right walls, as the right wall exerts a normal force in the $-x$ direction and the force due to the left wall is in the $+x$ direction. The time series of the perturbed wall normal stresses and their difference are shown in the inset of figure 10(a). If one oscillation of the centre of mass $x_{cm}$ is considered, $\sigma ^\prime _{xx}$ on each wall passes through one maximum and one minimum. Thus, the difference of $\sigma ^\prime _{xx}$ between two walls passes through one maximum and one minimum, and hence the dominant frequency of the difference $(\sigma ^\prime _{xx}|_{x = W} - \sigma ^\prime _{xx}|_{x = -W})$ is the same as that of $x_{cm}$. As the time series of $\sigma ^\prime _{xx}$ on either wall is not symmetric about its time-averaged value with a sharp maximum and a shallow minimum, the spectrum contains all the frequencies $f, 2f, 3f, \ldots$. However, the difference in $\sigma ^\prime _{xx}$ between two walls is necessarily symmetric about its time-averaged value which is zero, and therefore it only contains the frequencies $f, 3f, 5f, \ldots$. Because of this, the spectra of $x_{cm }$ and $u_{x|cm}$ contain only odd frequencies.

2. (ii) The momentum balance in the $y$ direction (for the perturbed state), equivalent to (4.1), is

   (4.2)\begin{equation} 2W \bar{\phi} \frac{\text{d} u_{y|cm}^\prime}{\text{d} t} ={-} (|\sigma_{xy}^\prime |_{x=W} + |\sigma_{xy}^\prime |_{x={-}W}). \end{equation}
   On the right-hand side of (4.2) is the negative of the sum of the absolute values of the shear stresses on the two walls, because the shear stress acts opposite to the flow direction on both the walls. The perturbed shear stress $\sigma ^\prime _{xy}$ on each wall has one maximum and one minimum displaced by a phase ${\rm \pi}$ over one cycle of the centre of mass $x_{cm}$. Therefore, the sum of these two stresses contains two maxima and two minima over one cycle (figure 10b). This has frequency $2f$, if $f$ is the frequency of $x_{cm}$. Because of this, the dominant frequency of $u_{y|cm}^\prime$ is two times that of $x_{cm}$ and $u_{x|cm}$. As the sum of the perturbed shear stresses at the walls is not symmetric about its time-averaged value, the higher harmonics are at $4f, 6f, \ldots$.

4.3. Spectra near bifurcation

Let us consider the variation of the amplitude and frequency of the dominant peak with $\bar {\phi }$. As expected, the amplitude of $x_{cm}$, denoted by $A_{x_{cm}}$, increases as $\bar {\phi }$ decreases (the solid curves in figure 11a). The values are approximately independent of $2W$, and the dependence on $2W$ is accounted for by the variation of $\bar {\phi }_{cr}$ (figure 3). The amplitude of $u_{x|cm}$ is related to $A_{x_{cm}}$ by $2 {\rm \pi}f_{x_{cm}} A_{x_{cm}}$ if the variation of $x_{cm}$ with time $t$ is sinusoidal. The inset in figure 11(a) shows that this is indeed the case, even though the variation is not strictly sinusoidal. The variation of $A_{u^\prime _{y|cm}}$ with $(\bar {\phi }_{cr}-\bar {\phi })$ shows a similar trend as that of $A_{x_{cm}}$ (the dashed curves in figure 11a). The data can be fitted by the power law $A_{x_{cm}}, A_{u^\prime _{y|cm}} \propto (\bar {\phi }_{cr}-\bar {\phi })^{1/2}$, as shown by the pink lines in figure 11(a). The variation of the dominant frequency is shown as a function of $(\bar {\phi }_{cr} - \bar {\phi })$ in figure 11(b), where the frequency does not depend on $2W$ significantly. However, it decreases with $(\bar {\phi }_{cr}-\bar {\phi })$ which implies that the mass takes a longer time to travel a longer distance, as expected. The frequency appears to tend to a constant value for $\bar {\phi }_{cr} - \bar {\phi } \ll 1$. Both of these are characteristics of a supercritical Hopf bifurcation with a non-zero frequency at the bifurcation point.

Figure 11. Variation of (a) the amplitude and (b) the frequency of the dominant peak in the amplitude–frequency spectra of $x_{cm}$ and $u^\prime _{y|cm}$ with $(\bar {\phi }_{cr} - \bar {\phi })$. The solid and dashed curves are for $x_{cm}$ and $u^\prime _{y|cm}$, respectively; $\diamond$, $2W = 60$; $\square$, $2W = 80$; $\circ$, $2W = 100$; $\triangle$, $2W = 120$; $f_{sf} = 1/(500 \Delta t)$. The inset in panel (a) is for $2W = 100$. The magenta curves represent the fit $y = a x^{1/2}$ in panel (a), where $a$ is a constant; $a = 14$ for the solid curve and $1.5$ for the dashed curve.

4.4. Sustained oscillations

The oscillations cannot be modelled by considering the granular material as a block colliding sequentially with the two walls with inelastic or frictional collisions. The magnitude of the horizontal velocity will decrease after each successive collision with the wall, and the vertical velocity will increase due to the gravitational acceleration, resulting in an accelerated state. Even if the collisions are frictional, i.e. the vertical impulse at a wall collision depends on the horizontal momentum of the block, the system will eventually reach an accelerated state. For sustained oscillations, the exchange of the fluctuating energy between the vertical and the horizontal fluctuations is essential.

Consider a simple model for the change in the velocity fluctuations $v_x^\prime$ and $v_y^\prime$ of the material before and after collisions with the two walls,

(4.3)$$\begin{gather} (v_x^{\prime})^2|_{t + \tau} = e_x^2 (v_x^{\prime})^2|_{t} + k [(v_y^{\prime})^2 - (v_x^{\prime})^2 ] |_{t}, \end{gather}$$

(4.4)$$\begin{gather}(v_y^{\prime})^2|_{t + \tau} = e_y^2 (v_y^{\prime})^2|_{t} + k [(v_x^{\prime})^2 - (v_y^{\prime})^2 ] |_{t} + g \left.\frac{2 A v_y^\prime}{v_x^\prime} \right|_{t}, \end{gather}$$

where $\tau$ is the time period for successive collisions. Here, $e_x$ and $e_y$ are the effective coefficients of restitution (not to be confused with coefficients of restitution of a particle in a granular medium) for the decrease in the fluctuating velocities after successive collisions with the two walls. The term proportional to $k$ results in an exchange of energy from the vertical to horizontal direction for $(v_y')^2>(v_x')^2$ and vice versa. As the energy dissipation during collisions is expected, both $e_x$ and $e_y$ must be less than $1$, and they can be functions of the coefficient of friction and damping constants. Here, $k$ is positive so that there is a transfer of energy from the direction with higher energy to that with lower energy. The last term on the right in (4.4) is the total change in the energy per unit mass after collisions with the two walls due to the acceleration in the vertical direction. The change in momentum per unit mass is the product of the gravitational acceleration $g$ and the time $(2 A /v_x^\prime )$ required to travel the distance $2 A$ in the horizontal direction, where $A$ is the amplitude of the oscillations. The energy input per unit mass is the product of the change in momentum and the velocity $v_y^\prime$.

There is an oscillatory state if $v_x^\prime |_{t+\tau } = v_x^\prime |_{t}$ and $v_y^\prime |_{t+\tau } = v_y^\prime |_{t}$,

(4.5)$$\begin{gather} (v_x^\prime)^2 = m^2 (v_y^\prime)^2 = \frac{m 2 A g }{(1 - e_y^2 + k (1 - m^2))}, \end{gather}$$

(4.6)$$\begin{gather}(v_y^\prime)^2 = \frac{2 A g}{m (1 - e_y^2 + k (1 - m^2))}, \end{gather}$$

where $m^2 = (k/(1 - e_x^2 + k))$. As $e_x$ and $e_y$ are less than $1$ and $k$ is positive, solutions for the fluctuating velocities exist only if $k$ is non-zero. For $k=0$, there is an accelerating state where $v_x^\prime$ decreases to zero and $v_y^\prime$ diverges. Thus, this simple model demonstrates that the transfer of energy between the horizontal and vertical velocity fluctuations is essential for sustained oscillations; oscillations cannot be captured by a lumped parameter model where the granular material is treated as a block. The rotational degree of freedom of the particles and the inter-particle and particle–wall friction facilitate this transfer of energy, thereby stabilising the oscillations. Similarly, the inter-particle damping seems to disrupt the transfer of energy between the vertical and horizontal directions, thereby resulting in the absence of oscillations. The effect of the coefficient of friction and the damping constants on the flow was summarised in § 3.5, and is examined in Appendix C. To identify the parameter regimes and mechanisms for oscillations, a stability analysis that includes the nonlinear interactions is required. This is more complicated than the linear stability analysis and is not attempted here.

5. Discussion

The flow of a granular material through a vertical channel between two flat frictional walls is considered in the present analysis with the periodic boundary conditions in the flow and the spanwise directions. Different flow regimes are observed for a fixed value of $2W$, when the average solids fraction $\bar {\phi }$, the ratio of the volume of the solids to the volume of the channel, is varied. There is no flow when $\bar {\phi } > \bar {\phi }_{max}$. The steady fully developed flow ($\bar {\phi }_{max} \geq \bar {\phi } \geq \bar {\phi }_{cr}$) has been discussed by Debnath et al. (2022a,b). When $\bar {\phi }$ decreases below a minimum value $\bar {\phi }_{min}$, there is an accelerated flow where the material is in the state of vertical free fall. In between these two regimes, there is a new ‘oscillatory’ regime for $\bar {\phi }_{cr} > \bar {\phi } \geq \bar {\phi }_{min}$, where there are cross-stream oscillations of constant amplitude. The regime of the oscillatory flow is the focus of the present work. If the flow is initiated from a static configuration, there are initial transients where the material is in free fall and shows transverse oscillations of growing amplitude. Here the transient states and the accelerated flow ($\bar {\phi } < \bar {\phi }_{min}$) are not examined in detail.

When the material oscillates, the variation of the solids fraction in the cross-stream direction at different time instances is asymmetric. During half a period of an oscillation, it is compressed towards one wall and is very loosely packed near the other wall. That results in oscillations in the wall stresses. The resistance to the flow by the two walls are not equal instantaneously, even though they are equal in the time-averaged sense. The profiles of the wall normal and shear stresses with the time exhibit sharp maxima and shallow minima. The maximum stresses are 3–4 times the time-averaged values. This could be important for practical applications, as it is necessary to design equipment to withstand instantaneous stresses which are multiples of the time-averaged values.

The wall normal and shear stresses vary continuously as $\bar {\phi }$ is decreased. There is a transition from the steady fully developed regime to the oscillatory regime with a distinct change in the slope of the stresses with respect to $\bar {\phi }$. An interesting observation is that the time-averaged base state of the normal stress is independent of $\bar {\phi }_{cr} - \bar {\phi }$ in the oscillatory regime, in contrast to that in the steady fully developed regime. The analysis reveals a complex relationship between the wave forms and the spectra of the different quantities. If the dominant frequency of the centre of mass is $f$, some dynamical quantities oscillate with frequencies $f, 3f, 5f, \ldots$, some oscillate with frequencies, $2f, 4f, 6f, \ldots$, and some with frequencies $f, 2f, 3f, \ldots$. The reasons for these are explained on the basis of the momentum balances.

The amplitude of the $x$ coordinate of the centre of mass and the perturbed state of the vertical velocity of the centre of mass vary as $\sim (\bar {\phi }_{cr} - \bar {\phi })^{1/2}$. For lumped parameter systems, such a scaling suggests a supercritical Hopf bifurcation (Strogatz 1994). However, it is shown that the sustained oscillations cannot be explained on the basis of a lumped-parameter model, where the granular material is considered as a block colliding with the two walls. The exchange of the fluctuating energy from the vertical to horizontal velocity fluctuations is essential for the sustained oscillations. Particle rotation, inter-particle and particle–wall friction and dissipation seem to be essential for this transfer of energy.

Dissipation plays a key role in the flow to reach a steady state or an oscillatory state. There are two sources of dissipation – damping and friction. It is observed that the reduction in the damping dissipation has a stabilising influence, whereas the reduction in the frictional dissipation has a destabilising influence. For a fixed width of the channel, the stability of the flow appears to depend on $\bar {\phi }$ and the rate of dissipation, which depends on the inter-particle coefficient of friction $\mu _p$ and the inter-particle damping ($\xi _n$, $\xi _t$). An energy balance analysis was attempted by Debnath (2023) to examine the reasons for the oscillations. Using a pseudo thermal energy and mechanical energy balance, it was hypothesised that the different states are due to the relative magnitudes of the rate of energy input $\dot {E}_t$ due to gravity and the rate of dissipation $\dot {D}_t$ due to damping and friction. The material will be jammed for $\dot {D}_t > \dot {E}_t$ and accelerating for $\dot {D}_t < \dot {E}_t$. A steady flow is observed if $\dot {D}_t = \dot {E}_t$ instantaneously. An oscillating flow is expected if the time averages of $\dot {D}_t$ and $\dot {E}_t$ are equal, but they are not equal at each instant. The energy balance analysis of Debnath (2023) is preliminary and does not provide a detailed derivation for the mechanism of instability. A more rigorous approach such as the nonlinear stability analysis is suggested which has the potential to unravel the mechanism and to identify different parameter regimes.