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.


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.
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 φ.
Further collisions allow the particles to spread across the cross-section and hence a reduction in φ 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 φ 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 φ (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 φ 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 φ = 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 φ, 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 ( φmax , φcr ) for which a steady flow can be obtained.The flow ceases when φ is above φmax , and there is an oscillatory flow followed by an accelerated flow when φ is below φ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 φ exceeds 0.587, which is the solids fraction for the arrested dynamics (Kumaran 2008); a loose shearing zone where φ 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 φ 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.

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.
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 − H, and there is a vacant head space of height H above the particles (figure 1a).The average solids fraction φ, the ratio of the volume of the particles to the volume of the channel, of the static bed is close to φ 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 φ is set by changing 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 − 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 φ is close to φ 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 (ξ n , ξ t ), and the coefficient of inter-particle friction μ p .Unless otherwise specified, the results correspond to k n = 10 6 ρ p gd 2 p , k t = 2/7k n , ξ n = 180 g/d p , ξ t = 1/2ξ n and μ p = 0.5, where ρ 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 /(ρ p gd 2 p ) > 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 , d p /g, gd p , gd p , ρ p gd 3 p and ρ p gd p , respectively.Without loss of generality, d p , ρ p and g are set equal to 1.

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 , 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 (i) There is no flow for φ > 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 φ = 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 ≡ t * = 100.After t * = 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 φ = 0.57, the initial transients are similar to those of φ = 0.60.The amplitude and frequency of x cm increase with time till t * ≈ 1050.There is a sharp transition in the profile of −u y|cm after t * ≈ 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 φ 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 √ 2gH.Hence, figure 2 implies that the flow reaches an oscillatory state after traversing a distance ∼7 × 10 3 H for H = 66 and φ = 0.57.However, for H = 68.15 and φ = 0.552, the free fall persists even after traversing a distance of ∼110 × 10 3 H.Due to constraints on the computational time, the motion of the particles could not be simulated beyond t = 4000.
When φ is low, the resistance to flow from both the walls is not experienced significantly by the material in the transition period (t < t * ).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 * .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-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 φ for the different regimes decreases.The boundary between no flow and a steady flow at φ ≡ φ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 φcr ≤ φ ≤ φ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 φ is decreased below φcr .An accelerated flow is observed when φ is further decreased below a value φmin , where the material is in free fall.The oscillatory regime in the range φmin ≤ φ < φcr is the focus of the present work.

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 π/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 state, a variable is defined as the sum of its time-averaged 'base' state b , which is time-independent, and the oscillations about the base state , The time-averaging is done over 5 × 10 5 time steps.Here, x b cm and u b x|cm are both zero, and u b y|cm / = 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.

Solids fraction
The variation of the solids fraction φ and the wall stresses are shown in figure 5.The time series of φ 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, φ at the the walls (the solid and dashed curves in figure 5a) shows oscillations in the range φ = 0.08-0.3with an average ∼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 φ passes through a maximum, and are close to zero when φ 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).
The solids fraction φ and the normal and shear stresses at the two walls are anti-correlated (figure 5a,b), i.e. when φ 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 φ at that wall is close to 0.3; at the same instant, the other wall experiences nearly no traction when φ is approximately 0.08.The profiles of φ corresponding to the time instants when the wall stresses are maximum and minimum are shown in figure 5 about the centreline.The distance between the time-averaged value of φ at the wall and the corresponding values of the profiles of φ 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 ± 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 2A x cm in half a period of oscillation.The corresponding value of φ in this layer of dimensions H × B × 2(W − A x cm ) is ≈ φ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.

Flow velocity
As mentioned earlier, the time-averaged flow velocity u b y|cm / = 0, and it exhibits a monotonic variation with 2W and φ (figure 6).Here, u b y|cm is shown both for the steady fully developed state ( φ ≥ φcr ) and the oscillatory state ( φ < φcr ) -the two are separated by the pink filled circles in figure 6(a).The magnitude of u b y|cm increases as φ 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 / 2(W − A x cm ), is plotted as a function of ( φcr − φ) 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 ( φcr − φ) → 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 φ ≥ φcr as φcr is a function of 2W (see figure 3).However, our earlier work (Debnath et al. 2022b)  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 φ is close to φmin .xx decreases as φ is decreased in the steady fully developed flow.However, the normal stress shows only a small variation with φ in the oscillatory state (figure 7a).As the variation of φ is in a very small range (<6 %), the variation of σ b xy at the wall with φ 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 ( φcr − φ), as shown by the red dashed curves in figure 7.Because φcr varies with 2W, the scaling is not suitable in the steady fully developed state for ( φcr − φ) ≤ 0. However, our earlier work (Debnath et al. 2022b) has shown that the scaling is suitable if it is a function of φ for the steady fully developed flow.
It is interesting to note that σ b xx is approximately independent of ( φcr − φ) > 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 σ b xx is due to the contact forces which decrease with φ (intuitively) and perhaps saturate as φ → φ+ 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 σ xx and not the time-averaged base state σ b xx .Because of this, σ b xx is independent of φ in the oscillatory regime for a given width.As 2W increases, σ b xx increases for a given value of φcr − φ 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 μ b w = |σ b xy |(x = W)/σ b xx , appears to be a constant for the time-averaged base states in the oscillatory flow, which is independent of φ and 2W.In contrast, the coefficient of friction strongly depends on φ and does not depend on 2W in the steady fully developed state (Debnath et al. 2022b).The wall stresses, σ xx and σ xy , and their perturbed states, σ xx and σ xy , 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 |σ xy |/σ xx and |σ xy |/σ 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.

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 μ p = 0.5 to 0.25 (keeping inter-particle damping constant ξ n fixed), the material accelerates even if φ is in the range φmax ≥ φ ≥ φ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 μ p appears to be significant.(ii) The effect of the inter-particle damping constants ξ n,p−p and the particle-wall damping constants ξ 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. ξ n,p−p = 0 (keeping μ p = 0.5 fixed) in the range φcr > φ ≥ φ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, ξ 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.

x cm and u y|cm
The amplitude-frequency spectra of x cm and u y|cm are shown in figure 9 for 2W = 100 and different values of φ.For φ = 0.61 (red curve), the frequency spectrum is noisy and exhibits no clear maxima.As φ is decreased, one clear maximum appears for φ = 0.59 (blue), but the amplitude of oscillations is found to be approximately 5-7 times smaller than a particle diameter.As φ 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 y|cm at the dominant frequencies denoted by f x cm and f u y|cm , respectively, and smaller maxima at the higher harmonics.The amplitudes of the higher harmonics of x cm and u y|cm are orders of magnitude smaller than that at the dominant one.For definiteness, φ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 φcr is slightly changed, i.e. the value is reduced from 0.5 to 0.2 particle diameter (60 % reduction), the change in φcr is observed to be in the range 0.5-0.8%).
For example, the dominant frequency of x cm is f x cm = 0.18 for φ = 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, . . .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 y|cm is 2f 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 y|cm are 2, 3, . . .times the dominant frequency f u y|cm , as the time series of u y|cm is not symmetric about its mean value.

Wall stress
The spectra for the perturbed normal and shear stresses, σ xx and σ 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, . . .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.
(i) The momentum balance in the cross-stream x direction (for the perturbed state) is 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, σ xx on each wall passes through one maximum and one minimum.Thus, the difference of σ xx between two walls passes through one maximum and one minimum, and hence the dominant frequency of the difference ( is the same as that of x cm .As the time series of σ 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 , . ... However, the difference in σ 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 , . ... 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 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 σ xy on each wall has one maximum and one minimum displaced by a phase π 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 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 , . ...

Spectra near bifurcation
Let us consider the variation of the amplitude and frequency of the dominant peak with φ.As expected, the amplitude of x cm , denoted by A x cm , increases as φ 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 φcr (figure 3).The amplitude of u x|cm is related to A x cm by 2π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 y|cm with ( φcr − φ) 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 y|cm ∝ ( φcr − φ) 1/2 , as shown by the pink lines in figure 11(a).The variation of the dominant frequency is shown as a function of ( φcr − φ) in figure 11(b), where the frequency does not depend on 2W significantly.However, it decreases with ( φcr − φ) 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 φcr − φ 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 y|cm with ( φcr − φ).The solid and dashed curves are for x cm and u y|cm , respectively; , 2W = 60; , 2W = 80; •, 2W = 100; , 2W = 120; f sf = 1/(500 t).The inset in panel (a) is for 2W = 100.The magenta curves represent the fit y = ax 1/2 in panel (a), where a is a constant; a = 14 for the solid curve and 1.5 for the dashed curve.

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 and v y of the material before and after collisions with the two walls, where τ 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 (2A/v x ) required to travel the distance 2A 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 .
There is an oscillatory state if where . 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 decreases to zero and v y 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.

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 φ, the ratio of the volume of the solids to the volume of the channel, is varied.There is no flow when φ > φmax .The steady fully developed flow ( φmax ≥ φ ≥ φcr ) has been discussed by Debnath et al. (2022a,b).When φ decreases below a minimum value φ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 φcr > φ ≥ φ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 ( φ < φ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 φ 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 φ.An interesting observation is that the time-averaged base state of the normal stress is independent of φcr − φ 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 , . .., some oscillate with frequencies, 2f , 4f , 6f , . .., and some with frequencies f , 2f , 3f , . ... 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 ∼ ( φcr − φ) 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 φ and the rate of dissipation, which depends on the inter-particle coefficient of friction μ p and the inter-particle damping (ξ n , ξ 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 Ėt due to gravity and the rate of dissipation Ḋt due to damping and friction.The material will be jammed for Ḋt > Ėt and accelerating for Ḋt < Ėt .A steady flow is observed if Ḋt = Ėt instantaneously.An oscillating flow is expected if the time averages of Ḋt and Ė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.were carried out with larger heights for 2W = 80-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 φ = 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 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 ≥ 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.

Appendix C. Effect of dissipation and friction
In the soft-sphere DEM simulations, there are two sources of dissipation, namely viscous damping ξ n and Coulomb friction μ p .If the friction coefficients, μ p and μ 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 μ p = μ wall = 0.25 without changing ξ n , and only an accelerated state has been observed even if φ is in the range φcr ≤ φ ≤ φmax .Therefore, a stable flow becomes unstable if the frictional dissipation is reduced.Henceforth, we set μ p = μ 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 ξ n,p−p = 180 > 0 and a particle-wall damping coefficient ξ 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 (ξ n,p−p > 0, ξ n,p−w = 0) (case B), (ξ n,p−p = 0, ξ n,p−w > 0) (case C), and (ξ n,p−p = 0, ξ n,p−w = 0) (case D), respectively (see table 1).(Note that the tangential damping constant, ξ t , is set to 2/7 ξ n ; hence, if ξ n is set to zero, ξ t = 0.) We expect the extent of dissipation to be in the order case A > case B > case C > case D.  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 φcr − φ 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).
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 φ, 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 φ b , velocity u b y , rate of deformation 1 2 |du y /dx| b and granular temperature T b among cases A to D is shown as a function of the distance x = x + W from the left wall in figure 17.Here, the domain has been separated into zones following the work of Debnath et al. (2022b) -the plug zone at the centre where the solids fraction is a constant and ≥0.63, the adjoining dense shearing zone where the solids fraction is in the range 0.587-0.63and the loose shearing zone where the solids fraction is less than 0.587.
Figure 17(a) shows that φ b is approximately similar for all the cases.However, φ b in the plug zone is slightly higher for cases A and B, where it is close to the dense random packing φ drp = 0.64.In contrast, φ 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, φ b is lower for cases A and B where there are oscillations.Because of the oscillation, φ 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 φ b near the wall.As φ is the same for cases A to D, φ 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.

Figure 1 .
Figure 1.Flow in a vertical channel.The dimensions of its rectangular cross-section in the z and x directions are B × 2W.(a) A static bed is filled up to height H − 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.

Figure 2 .
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; , φ = 0.625 (no flow); , φ = 0.60 (steady fully developed flow); •, φ = 0.57 (oscillatory flow); , φ = 0.552 (accelerated flow).The time series data are shown for f sf = 1/(500 t), where f sf is the sampling frequency (see Appendix A) and t is the simulation time step.

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

Figure 4 .
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 φ = 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.

Figure 5
Figure 5.A typical time series of (a) the solids fraction φ at the walls and (b) wall stresses σ xx (black, left ordinate) and σ xy (red, right ordinate).In panels (a) and (b), the solid curves represent the left wall at x ≡ 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 φ; the dotted (black and red) horizontal lines in panel (b) represent the time-averaged values of σ xx and |σ xy |.(c) Variation of φ with the distance from the left wall x .In panel (c), and represent profiles for time t 1 = 6001.2and t 2 = 6003.6when the material is fully compressed to the left and right walls, respectively, • represents the time-averaged base state φ b , and the inset is a magnification of the profiles near the right wall.The parameter values are 2W = 100 and φ = 0.57.The time series data are shown for f sf = 1/(500 t), where f sf is the sampling frequency (see Appendix A) and t is the simulation time step.

Figure 6 .
Figure 6.Variation of (a) u b y|cm with φ and (b) its scaled values with φcr − φ.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: , 2W = 60; , 2W = 80; •, 2W = 100; , 2W = 120.

Figure 7 .
Figure 7. Variation of the time-averaged stresses (a) σ b xx and (b) σ b xy , and their scaled values with φcr − φ 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: , 2W = 60; , 2W = 80; •, 2W = 100; , 2W = 120.

3. 4 .
Stress The time-averaged stresses are shown in figure 7. From the momentum balance, the normal stress σ b xx is a constant across the channel, i.e. σ xx ≡ σ b xx = constant, and the magnitude of σ b xy increases from the centre to the wall for the steady fully developed flow (Debnath et al. 2022a,b).The maximum value of |σ b xy | at the wall is equal to the weight of the material in one half of the channel, i.e. |σ b xy |(x = W) = W φ 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 σ b xx at the transition; -σ b

Table 1 .
Damping coefficients for different cases for inter-particle and particle-wall interactions.
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, φ = 0.57.