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 Reference Rao and Venkateswarlu1975; Tüzün & Nedderman Reference Tüzün and Nedderman1985; Roberts & Wensrich Reference Roberts and Wensrich2002; Ramírez, Nielsen & Ayuga Reference Ramírez, Nielsen and Ayuga2010; Wang et al. Reference Wang, Liang, Guo, Chen, Liu, Ma, Chen and An2020). The ‘stick-slip’ motion of particles past the wall and ‘silo-quaking’ are some common sources of such fluctuations (Tüzün & Nedderman Reference Tüzün and Nedderman1985; Tejchman & Gudehus Reference Tejchman and Gudehus1993; Dhoriyani et al. Reference Dhoriyani, Jonnalagadda, Kandikatla and Rao2006). 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 Reference Hopkins and Louge1991; Goldhirsch & Zanetti Reference Goldhirsch and Zanetti1993; Raafat, Hulin & Herrmann Reference Raafat, Hulin and Herrmann1996; Kudrolli, Wolpert & Gollub Reference Kudrolli, Wolpert and Gollub1997; Kumaran Reference Kumaran1998; Hua & Wang Reference Hua and Wang1999; Luding & Herrmann Reference Luding and Herrmann1999; Alam & Hrenya Reference Alam and Hrenya2001; Liss, Conway & Glasser Reference Liss, Conway and Glasser2002; Conway & Glasser Reference Conway and Glasser2004; Mitrano et al. Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011; Fullmer & Hrenya Reference Fullmer and Hrenya2017). 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 Reference Liss and Glasser2001; Conway, Liu & Glasser Reference Conway, Liu and Glasser2006). The kinetic theory based models for inelastic hard-particles predict a strong connection between inelastic dissipation and instability (Babic Reference Babic1993; Brey et al. Reference Brey, Dufty, Kim and Santos1998; Garzó Reference Garzó2005), which is in agreement with hard-particle simulations (Hopkins & Louge Reference Hopkins and Louge1991; Goldhirsch & Zanetti Reference Goldhirsch and Zanetti1993; Mitrano et al. Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011, Reference Mitrano, Dahl, Hilger, Ewasko and Hrenya2013).
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 (Reference Luding and Herrmann1999) 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. (Reference Raafat, Hulin and Herrmann1996) 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. (Reference Aider, Sommier, Raafat and Hulin1999) 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 (Reference Hua and Wang1999) 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 (Reference Hua and Wang1999) 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 (Reference Wang, Jackson and Sundaresan1997) 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. (Reference Lun, Savage, Jeffrey and Chepurniy1984) and the boundary conditions of Johnson & Jackson (Reference Johnson and Jackson1987). 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 (Reference Wang and Tong2001) 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. (Reference Liss, Conway and Glasser2002) 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 Reference Babic1993; Allen & Tildesley Reference Allen and Tildesley2017). The predictions of simulations by Liss et al. (Reference Liss, Conway and Glasser2002) 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 (Reference Wang and Tong2001).
The work of Wang & Tong (Reference Wang and Tong2001) 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. (Reference Liss, Conway and Glasser2002) 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 (Reference Debnath, Kumaran and Rao2019, Reference Debnath, Kumaran and Rao2022a); Debnath, Rao & Kumaran (Reference Debnath, Rao and Kumaran2022b); 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 Reference Kumaran2008); 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 (Reference Berzi, Jenkins and Richard2020) where the effect of the particle stiffness is incorporated, and applied for the present flow configuration by Debnath et al. (Reference Debnath, Kumaran and Rao2022a) and Islam, Jenkins & Das (Reference Islam, Jenkins and Das2022).
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.
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 Reference Plimpton1995). 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. (Reference Debnath, Kumaran and Rao2022a,Reference Debnath, Rao and Kumaranb).
The variation of the flow characteristics with 
$k_n$ in the plug region was examined by Debnath et al. (Reference Debnath, Rao and Kumaran2022b) (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}$,
\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.
(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.(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.(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.(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$.
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. Reference Debnath, Kumaran and Rao2019.) 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. (Reference Debnath, Kumaran and Rao2022a,Reference Debnath, Rao and Kumaranb). 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.
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.
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$,
\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).
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. Reference Debnath, Rao and Kumaran2022b) 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 Reference Beverloo, Leniger and Van de Velde1961) 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}$.
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. Reference Debnath, Kumaran and Rao2022a,Reference Debnath, Rao and Kumaranb). 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. Reference Debnath, Rao and Kumaran2022b) has shown that the scaling is suitable if it is a function of 
$\bar {\phi }$ for the steady fully developed flow.
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. Reference Debnath, Rao and Kumaran2022b, (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 Reference Rao and Nott2008).
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. Reference Debnath, Rao and Kumaran2022b). 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.
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:
(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.(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 %).
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.
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.
(i) The momentum balance in the cross-stream
$x$ direction (for the perturbed state) is
(4.1)The left-hand side of the above equations is the scaled mass times acceleration in the\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}$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.(ii) The momentum balance in the
$y$ direction (for the perturbed state), equivalent to (4.1), is
(4.2)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\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}$\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.
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,
$$\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}$$
$$\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}$,
$$\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}$$
$$\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. (Reference Debnath, Kumaran and Rao2022a,Reference Debnath, Rao and Kumaranb). 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 Reference Strogatz1994). 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 (Reference Debnath2023) 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 (Reference Debnath2023) 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.
Acknowledgements
We are grateful to the Supercomputer Education and Research Centre (SERC), IISc for providing computational support.
Funding
This work was supported by funding from the MHRD and the Science and Engineering Research Board, Government of India (grant no. CRG/2018/002673 and SR/S2/JCB-31/2006) and Synopsys.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Effect of the sampling frequency
The error bars in the results are calculated as follows. When the flow oscillates with a frequency after a period of time, samples are collected over a time interval, marked by the blue box in figure 12(a). This time interval is further divided into ten sub-intervals (red boxes in figure 12b). For each sub-interval, the amplitude and frequency are calculated (figure 12c). The error bars represent 95 % confidence limits and are obtained from these ten sub-intervals. In many results, the error bars are smaller than the size of the symbols and are not shown.
Error bar calculation for the amplitude and frequency of the properties of the centre of mass.
To check the effect of sampling frequency on the results, the data are collected for four different sampling frequencies, 
$f_{sf} = 1/(500 \Delta t)$, 
$1/(1000 \Delta t)$, 
$1/(2500 \Delta t)$ and 
$1/(5000 \Delta t)$ in a time interval, where 
$\Delta t$ is the simulation time step 
$1.2\times 10^{-4}$. The variations of the amplitude and frequency corresponding to the dominant peak in the amplitude-frequency spectrum of 
$x_{cm}$ with 
$\bar {\phi }$ is independent of 
$f_{sf}$ (figure 13).
For 
$2W = 100$, variation of (a) 
$A_{x_{cm}}$ and (b) 
$f_{x_{cm}}$ with 
$\bar {\phi }$ for different sampling frequency 
$f_{sf}$: 
$f_{sf}$ = 
$1/(500 \Delta t)$ (
$\circ$), 
$1/(1000 \Delta t)$ (
$\square$), 
$1/(2500 \Delta t)$ (
$\triangle$), 
$1/(5000 \Delta t)$ (
$\diamond$).
Appendix B. Effect of the height to width ratio
The effect of system size on the dynamics of the oscillatory flow is examined in this section. Earlier studies have reported that the ratio of the height to the width of the channel, 
$(H/2W)$, affects the pattern of the density wave in the flow and the cross-stream directions (Wang & Tong Reference Wang and Tong2001; Liss et al. Reference Liss, Conway and Glasser2002). For example, Wang & Tong (Reference Wang and Tong2001) predicted a symmetric density wave type 1 for 
$(H/2W) = 71$, a symmetric density wave type 2 for 
$(H/2W) = 24$ and an antisymmetric density wave for 
$(H/2W) = 2.7$. In all cases, the variations in the flow direction are evident. The simulations of disks flowing through a two-dimensional channel by Liss et al. (Reference Liss, Conway and Glasser2002) captured similar types of waves for 
$(H/2W) \geq 2$; they are absent in 
$(H/2W) \leq 1$, where 
$2W = 33.3$ and 
$66$.
The present study considers a three-dimensional system which has an additional degree of freedom. The values of 
$\bar {\phi }$ are significantly higher and the results in the earlier sections are shown for 
$(H/2W)$ in the range 0.5–1.1. The results discussed earlier imply that the confinement effect may not become significant for 
$2W \geq 80$. Therefore, a few simulations were carried out with larger heights for 
$2W = 80\unicode{x2013}120$ to investigate the effect of 
$(H/2W)$ on the results. The system with the largest height considered is 
$(H/2W) = 2.5$ for 
$2W = 80$ and 
$\bar {\phi } = 0.565$.
The wall normal stress at five different heights for larger values of 
$(H/2W)$ in the range 2–2.5 is shown in figure 14. Simulations for these different heights show identical variations in time, and there is no change in the normal stress with stream-wise positions. This indicates that there is no stream-wise variation for a system where the ratio 
$(H/2W)$ is 
$\lesssim$2.5, and the material oscillates only in the cross-stream direction. Because of computational constraints, we could not simulate systems having 
$(H/2W) > 2.5$ and 
$2W \geq 80$. It is speculated that the variations in both the flow and cross-stream directions may appear in systems if the value of 
$(H/2W)$ is much higher. Simulating such large systems is beyond the scope in the present work.
Variation of the wall normal stress 
$\sigma _{xx}|_{x^\prime = 0}$ with the time 
$t$ at different heights 
$y/d_p =$ 20 (
$\circ$), 60 (
$\bigtriangledown$), 100 (
$\triangle$), 140 (
$\square$), 180 (
$\diamond$). Parameter values: black curves in panel (a), 
$(H/2W) = 2.5$, 
$2W = 80$, 
$\bar {\phi } = 0.565$; red curves in panel (b), 
$(H/2W) = 2$, 
$2W = 100$, 
$\bar {\phi } = 0.565$; 
$f_{sf} = 1/(5000 \Delta t)$.
Appendix C. Effect of dissipation and friction
In the soft-sphere DEM simulations, there are two sources of dissipation, namely viscous damping 
$\xi _n$ and Coulomb friction 
$\mu _p$. If the friction coefficients, 
$\mu _p$ and 
$\mu _{wall}$, vanish, a state of vertical free fall is attained, at least up to the time the equations are integrated. A few simulations are performed by setting 
$\mu _p = \mu _{wall} = 0.25$ without changing 
$\xi _n$, and only an accelerated state has been observed even if 
$\bar {\phi }$ is in the range 
$\bar {\phi }_{cr} \leq \bar {\phi } \leq \bar {\phi }_{max}$. Therefore, a stable flow becomes unstable if the frictional dissipation is reduced. Henceforth, we set 
$\mu _p = \mu _{wall} = 0.5$ and examine the roles of the inter-particle and particle–wall damping on the oscillations.
The black solid curve in figure 15 corresponds to an inter-particle damping coefficient 
$\xi _{n,p-p} = 180 > 0$ and a particle–wall damping coefficient 
$\xi _{n, p-w} = 180 > 0$ (case A). Here sustained oscillations are observed after some time, as discussed earlier. In figure 15, the black dotted, blue and red curves (referenced to the left ordinate) correspond to (
$\xi _{n,p-p} > 0$, 
$\xi _{n, p-w} = 0$) (case B), (
$\xi _{n,p-p} = 0$, 
$\xi _{n, p-w} > 0$) (case C), and (
$\xi _{n,p-p} = 0$, 
$\xi _{n, p-w} = 0$) (case D), respectively (see table 1). (Note that the tangential damping constant, 
$\xi _t$, is set to 
$2/7 \; \xi _n$; hence, if 
$\xi _n$ is set to zero, 
$\xi _t = 0$.) We expect the extent of dissipation to be in the order case A > case B > case C > case D.
Variation of 
$x_{cm}$ and 
$u_{y|cm}$ with the time 
$t$. Left ordinate: black solid curve, case A; black dotted curve, case B; blue solid curve, case C; red solid curve, case D (see table 1). Right ordinate: magenta dashed curve, case A; magenta dotted curve, case B; blue dashed curve, case C; red dashed curve, case D. Parameter values: 
$2W = 100$, 
$\bar {\phi } = 0.57$.
Damping coefficients for different cases for inter-particle and particle–wall interactions.
Case B is qualitatively similar to case A, except that the onset of oscillation is earlier and 
$-u_{y|cm}$ is lower in the former (see the dashed magenta curve (case A) and the dotted magenta curve (case B) referenced to the right ordinate in figure 15).
For the cases examined here, the oscillations disappear and a steady fully developed state is attained when the particle–particle damping is zero irrespective of the nature of particle–wall damping (cases C and D), and 
$-u_{y|cm}$ is higher for case C compared to case D. In the transition period, the material is always in free fall before the flow reaches either a steady fully developed state or an oscillatory state. The flow reaches a steady state quickly in case D where the dissipation due to damping is the least. Hence, figure 15 implies that a reduction in the damping dissipation has a stabilising influence, opposite to the role of the frictional dissipation.
The 
$y$ component of the velocity of the centre of mass is the highest for case A and lowest for case D (see the right ordinate in figure 15). The centreline velocity as a function of 
$\bar {\phi }_{cr} - \bar {\phi }$ shows the same trend in figure 16(a), irrespective of the width of the channel. This is counter-intuitive as the dissipation due to damping is higher for case A, and hence a lower velocity would be expected. The higher velocity for case A is mainly due to the higher slip velocity, as shown in figure 16(b). The variation of the velocity across the channel, which is an indication of the mean shear rate within the flow, is approximately similar for cases A and D (figure 16c). However, the slip velocity averaged over cycles for the oscillatory state is approximately two times higher for case A (oscillatory state) than that for case D (steady fully developed state).
Variation of (a) the time-averaged centreline velocity 
$(-u^b_y)|_{x = 0}$, (b) the slip velocity 
$(-u^b_y)|_{x = W}$, (c) their differences, (d) the normal stress (black curves, left ordinate) and its scaled value (red curves, right ordinate) with 
$(\bar {\phi }_{cr} - \bar {\phi })$: solid and dashed curves are for cases A and D, respectively; 
$\square$, 
$2W = 80$; 
$\circ$, 
$2W = 100$; 
$\triangle$, 
$2W = 120$.
The normal stress and its scaled value is measurably lower for case A in comparison to case D (figure 16d). This may be because of the contact stress which is expected to be higher for case D as there is no viscous damping. As the shear stress required to balance the weight of the material does not vary for a fixed value of 
$\bar {\phi }$, this indicates that the effective coefficient of friction at the wall is higher for case A than that for case D.
A comparison of the time-averaged base state of the solids fraction 
$\phi ^b$, velocity 
$u^b_y$, rate of deformation 
$\frac {1}{2} |\text {d}u_y/\text {d}\kern0.7pt x|^b$ and granular temperature 
$T^b$ among cases A to D is shown as a function of the distance 
$x^\prime = x + W$ from the left wall in figure 17. Here, the domain has been separated into zones following the work of Debnath et al. (Reference Debnath, Rao and Kumaran2022b) – the plug zone at the centre where the solids fraction is a constant and 
$\geq$0.63, the adjoining dense shearing zone where the solids fraction is in the range 0.587–0.63 and the loose shearing zone where the solids fraction is less than 
$0.587$.
Variation of the time-averaged quantities with the distance 
$x'$ from the left wall for 
$2W = 100$: (a) 
$\phi ^b$; (b) average velocity; (c) rate of deformation 
$\frac {1}{2} |\text {d}u_y/\text {d}\kern0.7pt x|^b$ and (d) granular temperature 
$T^b$. Parameter values: 
$\triangle$, case A; 
$\times$, case B; 
$+$, case C; 
$\circ$, case D; black curves, 
$\bar {\phi } = 0.57$; red curves, 
$\bar {\phi } = 0.565$.
Figure 17(a) shows that 
$\phi ^b$ is approximately similar for all the cases. However, 
$\phi ^b$ in the plug zone is slightly higher for cases A and B, where it is close to the dense random packing 
$\phi _{drp} = 0.64$. In contrast, 
$\phi ^b$ in the plug for cases C and D is close to 
$0.63$. We have verified that there is no crystalline structure in either case, so the higher packing for the oscillatory state is not due to crystallisation. Close to the wall, 
$\phi ^b$ is lower for cases A and B where there are oscillations. Because of the oscillation, 
$\phi$ at the left wall is very low when the centre of mass is at its right extreme, and vice versa. This results in a lower 
$\phi ^b$ near the wall. As 
$\bar {\phi }$ is the same for cases A to D, 
$\phi ^b$ in the plug region for cases A and B is expected to be higher than that for cases C and D. The significantly higher velocity for cases A and B than cases C and D is shown in figure 17(b). Figure 17(c) shows that the rate of deformation within the plug region decreases to zero within the simulation resolution, as indicated by the large error bars. The temperature at the plug is clearly non-zero because the error bars are smaller than the symbols (figure 17d). The granular temperature for cases A and B is higher than that for cases C and D in the plug and dense shearing zones, but lower near the walls in the loose shearing zone. The higher temperature in the dense and plug zones could be due to the higher rate of deformation in the dense shearing zone adjoining the plug for cases A and B. In the loose shearing zone, both the rate of deformation and the granular temperature are lower for cases A and B in comparison to cases C and D.
Thus, the different shear regimes for different cases are qualitatively similar and the locations of the boundaries are also the same to within simulation accuracy. However, the solids fraction, rate of deformation and granular temperature are higher near the centre, and lower at the wall for cases A and B where there are oscillations in comparison to cases C and D where there is no oscillation.
Open access
