Different shear regimes in the dense granular flow in a vertical channel

Abstract

The steady dense granular flow in a vertical channel bounded by flat frictional walls in one horizontal direction and with periodic boundary conditions in the other horizontal and vertical directions is studied using the discrete element method. The shape of the scaled velocity profile is characterized quantitatively by a universal exponential function, and the ratio of the maximum and slip velocities is independent of the average volume fraction $\bar {\phi }$ and the channel width $W$. For sufficiently wide channels, the maximum and slip velocities increase proportional to $\sqrt {W}$, and the thickness of the shearing zones increases proportional to $W$. There are four zones in the flow, each with distinct dynamical properties. There is no shear in the plug zone at the centre, but there is particle agitation, and the volume fraction $\phi$ is lower than the random close packing volume fraction $\phi _{rcp}$. In the adjoining dense shearing zone, $\phi$ is greater than the volume fraction for arrested dynamics $\phi _{ad} = 0.587$, and the granular temperature and shear rate depend on the particle stiffness. Adjacent to the dense shearing zone is the loose shearing zone where $\phi < \phi _{ad}$. Here, the properties do not depend on the particle stiffness, and the constitutive relations are well described by hard-particle models. The rheology in the loose shearing zone is similar to the dense flow down an inclined plane. There is high shear and a sharp decrease in $\phi$ in the wall shearing zone of thickness about two particle diameters, where the particle angular velocity is different from the material rotation rate due to the presence of the wall.

1. Introduction

The flow of particles in a vertical channel of infinite length bounded by parallel walls under the influence of gravity is one of the simplest examples of a granular flow. This is an approximation for the flow in the central region of a long vertical bin away from the free surface and the bottom exit. The flow is expected to be unidirectional and fully developed, and the velocity in the vertical direction is a function of only one cross-stream direction, if we assume that the system is infinite in the second cross-stream direction. This flow has been studied in experiments (Savage 1979; Nedderman & Laohakul 1980; Natarajan, Hunt & Taylor 1995; Pouliquen & Gutfraind 1996; Ananda, Moka & Nott 2008) and using theoretical models (Goodman & Cowin 1971; Savage 1979; Gutfraind & Pouliquen 1996; Mohan, Nott & Rao 1997, 1999; Barker, Zhu & Sun 2022; Debnath, Kumaran & Rao 2022). There is consensus that the flow consists of a central plug zone, where the material flows without shearing, and relatively thin layers close to the wall where there is shearing. These flows are usually very dense, and the volume fraction is close to random close packing, though there is a decrease in the volume fraction in the shearing zones close to the wall. From the momentum balance, the normal stress in the cross-stream direction is a constant. If the decrease in the volume fraction near the walls is neglected, then the shear stress increases proportional to the distance from the centre of the channel. The precise volume fraction is difficult to measure in experiments as it is near close packing, but the flow velocities have been measured. Despite the simplicity of the configuration, there is as yet no established model for the dependence of the velocity profile, volume fraction and shear layer thickness on the system parameters such as the channel width, particle diameter and average volume fraction.

In experimental configurations, the flow rate in the fully developed unidirectional flow far from the exit depends on the width of the exit slot and the exit conditions. The discharge rate at the exit is described by the Beverloo correlation (Beverloo, Leniger & Van de Velde 1961), where the velocity is proportional to the square root of the width of the exit slot. The discharge rate at the exit slot determines the average flow velocity far upstream of the exit slot, and the average volume fraction at that location is set by the flow velocity, overburden, and the particle–particle and particle–wall interactions. If the width of the channel is much smaller than the height, then the weight of the overburden at a given depth from the free surface is balanced by the frictional force exerted by the side walls (the Janssen effect; Janssen 1895), and the pressure is independent of depth; this is in contrast to the linear increase of pressure with depth for normal fluids. For a steady unidirectional flow, the total height of the channel does not affect the flow. Therefore, the flow profile in the central region of the channel should be amenable to a description based on the average volume fraction, the wall friction and the dimensions in that region, instead of the total height or the discharge rate fixed at the exit.

To examine a steady fully developed unidirectional flow, we consider a simpler configuration in simulations, where the channel is bounded by two flat frictional walls in the cross-stream direction, and periodic boundary conditions are imposed in the flow and spanwise directions. The average volume fraction of particles in the channel is fixed by the initial loading before the onset of flow. If the average volume fraction is small, then it is expected that the particles will accelerate continuously due to gravity, and no steady state will be reached. If the average volume fraction is near random close packing, then there will be a jammed state in which there is no flow. Here, we examine an intermediate range of the average volume fractions where there is a steady fully developed flow.

The salient feature of the flow in a channel is the presence of a plug zone at the centre where there is no shearing, and shear layers at the walls. Two possible states are predicted by frictional models, a static state and a plug flow with an indeterminate velocity. This is because frictional models do not have any rate dependence, and they do not contain any intrinsic length scale. A microscopic length scale has been included in different ways in models for dense granular flows. In the Cosserat models (Mohan et al. 1999; Mohan, Rao & Nott 2002), the particle spin is incorporated as an additional field; the length scale in this model is the distance from the walls over which there is a difference between the particle angular velocity and the material rotation rate due to the wall effect. Granular fluidity models (Aranson & Tsimring 2001; Henann & Kamrin 2013) define an additional field that is a measure of the extent of ‘fluidity’ of the material at a location. The fluidity is determined from a conservation equation similar to a diffusion equation (Henann & Kamrin 2013) or an order parameter equation describing the cross-over between flowing and static states (Aranson & Tsimring 2001). A different approach was followed by Dsouza & Nott (2020), where the flow rule was integrated over a representative volume with an effective mesoscopic size, to obtain constitutive relations containing the Laplacian of the rate of deformation tensor. However, these models require boundary conditions to be specified, and these boundary conditions are not prescribed. There is ambiguity in how the boundary conditions are to be imposed in different configurations. Some of these models have third and higher derivatives of the velocity, hence additional boundary conditions are required.

The incompressible $\mu (I)$ models (GDR-MiDi 2004; Jop, Forterre & Pouliquen 2006) attempt to incorporate rate dependence in the constitutive relations through the inertia parameter $I$. The friction coefficient $\mu$, the ratio of the shear and normal stresses, is expressed as a function of $I$, the shear rate non-dimensionalized by the square root of the normal stress along with suitable powers of the particle diameter and density. Universal relations, presented later in (B3) and (B4), have been proposed for the relation between the stress, volume fraction and inertia parameter. While the $\mu (I)$ rheology has been used widely for diverse applications, it is ill-posed in the sense that it is unstable to small wavelength perturbations. Well-posed models have been proposed by incorporating compressibility and inertia parameter in the flow rule and the yield function (Barker et al. 2017; Schaeffer et al. 2019).

Recently, Debnath et al. (2022) compared the predictions of some of these models with simulation results. None of the models predict all the profiles well, and some of the models could not be solved for some values of the parameters used. For example, the model of Barker et al. (2017) predicts that the volume fraction is almost a constant across the channel. The occurrence of negative fluidity in the model of Henann & Kamrin (2013) in some cases is an unrealistic feature. The model of Dsouza & Nott (2020) predicts thick shear layers and higher normal stress compared to results obtained using the discrete element method (DEM).

Another issue in applying the non-local model of Kamrin & Koval (2012) is the following. The vertical flow considered by Kamrin & Koval (2012) and Kim & Kamrin (2020) is unusual, because the pressure exerted on the walls is specified, and not the loading. Therefore, the width is permitted to vary, while the pressure is maintained as constant. This is rather unrealistic, because in a practical realization, one of the walls is not fixed, but is connected to springs that maintain a constant normal stress. This configuration was used specifically because the $\mu (I)$ parameter can be calculated directly from the known normal stress, and the shear stress can be determined from the body force per unit volume. In contrast, we are examining a flow where the average volume fraction is specified, but the normal stress is not.

The flow down an inclined plane is a configuration that has been studied widely using particle-based simulations (Silbert et al. 2001; Mitarai & Nakanishi 2005; Brewster et al. 2008; Reddy & Kumaran 2007, 2010). In three-dimensional simulations of particle assemblies, there is no flow when the angle of inclination is below a minimum value of about $20^\circ$. There is steady flow over a relatively small range of angles, ${\sim }20^\circ \unicode{x2013} 24^\circ$, and the flow becomes unstable when the angle of inclination exceeds about $25^\circ$. There are several intriguing features that are observed for a steady flow. From momentum balance, the ratio of the normal and shear stresses is equal to the tangent of the angle of inclination. The volume fraction is found to be independent of height, in contrast to the intuitive expectation of greater compaction with increasing depth due to the higher overburden. The ‘Bagnold law’ (Bagnold 1954, 1956) describes the relation between the stress and the velocity profile in the bulk, and the ‘granular temperature’ (defined in (2.3)) increases linearly with depth. A novel transition between an ordered and a disordered flow due to an increase in the base roughness has also been reported (Kumaran & Maheshwari 2012; Kumaran & Bharathraj 2013; Bharathraj & Kumaran 2017; Silbert et al. 2002). The relation between the stress and shear rate is found to obey the Bagnold law in both the ordered and disordered flows, though the Bagnold coefficients are very different.

The Bagnold law, which states that the stress is proportional to the square of the shear rate, follows from dimensional analysis if the inverse of the shear rate is the only time scale in the problem. The assumption is that the particle interactions are instantaneous, and this assumption is valid if the period of the interactions is much smaller than the inverse of the shear rate. The instantaneous collision assumption is made in kinetic theory of granular flows, where the particle velocity distribution function is defined in a manner similar to the molecular distribution function in the kinetic theory of gases (Lun et al. 1984; Jenkins & Richman 1985; Kumaran 1998). The mass, momentum and energy conservation equations are derived from the Boltzmann equation for the velocity distribution function. The energy conservation equation for the granular temperature contains an additional dissipation term due to particle interactions which is not present for a molecular fluid. In kinetic theory, it is also assumed that the system is dilute, and correlations between colliding particles are neglected due to the ‘molecular chaos’ assumption. Constitutive relations for dense granular flows have been derived using the Enskog procedure (Sela, Goldhirsch & Noskowicz 1996; Sela & Goldhirsch 1998; Kumaran 2004, 2006), and attempts have been made to go beyond the molecular chaos assumption and incorporate correlations in the pre-collisional velocities of the particles (Goldhirsch & van Noije 2000; Kumaran 2009a,b). With these refinements, it has been shown that the salient features of a dense granular flow in an inclined channel are captured adequately by models based on the hard-particle approximation (Kumaran 2008, 2014).

For an inclined plane flow, as the inclination angle is decreased, the maximum volume fraction in the flowing state is about 0.585–0.59, which is denoted the volume fraction for arrested dynamics, $\phi _{ad}$. The physical meaning of $\phi _{ad}$ is that an assembly of perfectly hard spheres (in which the repulsive force is zero when two particles are not in contact, and infinite when they are in contact) cannot be sheared if the volume fraction $\phi$ exceeds $\phi _{ad}$. For $\phi \geq \phi _{ad}$, shearing is possible only if the particles have finite stiffness. For $\phi < \phi _{ad}$, the Bagnold law based on the hard-particle model appears to apply even for dense granular flows. Kinetic models have been modified by incorporating a mesoscopic length scale representing particle chains to predict $\phi _{ad}$ in a sheared dense granular flow (Jenkins 2006, 2007; Berzi & Jenkins 2015).

The change in the rheology at $\phi = \phi _{ad}$ is consistent with the simulation studies of Kumaran (2009c,d). These showed that $\phi _{ad}$ for a dense sheared inelastic hard-particle fluid is a function of the coefficients of restitution ($e_n,e_t$) in the directions perpendicular and parallel to the surfaces at contact. The value of $\phi _{ad}$ is in the range 0.585–0.59 for smooth ($e_t = -1$) inelastic particles with $e_n$ less than 0.6, and 0.581–0.612 for perfectly rough ($e_t = 1$) inelastic particles with $e_n \geq 0.6$. It is significantly lower than the random close packing volume fraction $\phi _{rcp} = 0.64$ for an elastic hard-sphere fluid in the absence of shear. In contrast, the volume fraction $\phi$ in the plug zone could be significantly higher than $\phi _{ad}$ in the vertical channel. This raises the question of whether the nature of the flow in a vertical channel is qualitatively different from the flow down an inclined plane, or whether $\phi$ in the shearing zones is below $\phi _{ad}$ and the hard-particle model can be applied.

Berzi, Jenkins & Richard (2019, 2020) have separated the flow through an inclined chute with bumpy base and frictional side walls into two regimes: (i) a dense erodible layer near the bumpy base where $\phi \geq \phi _{c}$, a critical volume fraction; and (ii) a loose collisional layer above the erodible zone where $\phi < \phi _c$. The concept of $\phi _c$ can be linked to $\phi _{ad}$, which will be discussed shortly. The material in the erodible layer moves very slowly with negligible particle fluctuations compared to that in the collisional layer. The stress and dissipation rate in the collisional flow are adapted from the kinetic theory of granular flows with coefficients modified to incorporate the finite duration of contacts (Berzi & Jenkins 2015). In the erodible region, the rate-independent stresses arise because of sustained contacts, hence a rate-independent component is added to the pressure (Berzi et al. 2020). The form of the rate-independent term in pressure is motivated from the plane shear studies done by Chialvo, Sun & Sundaresan (2012). The latter have shown that the rate-independent effects are dominant for $\phi \geq \phi _c$ in the limit of vanishing shear rate, and $\phi _c$ does not depend on $e_n$. It varies with the coefficient of interparticle friction $\mu _p$; for example, $\phi _c = 0.587$ for $\mu _p = 0.5$, which Berzi et al. (2019, 2020) have adopted to mark the boundary between the erodible and collisional layers.

It is interesting to note that the range of $\phi _c$ for $\mu _p$ varying between $0$ and $1$ in the plane shear simulations of soft particles having finite stiffness (Chialvo et al. 2012; Berzi & Vescovi 2015) is similar approximately to $\phi _{ad}$ in plane shear simulations of hard inelastic rough particles having infinite stiffness (Kumaran 2009c, d). Chialvo et al. (2012) state that $\phi _c$ is similar to $\phi _{ad}$, and this volume fraction is an important parameter to explain arrested dynamics or jamming. In the vertical channel flow, we examine whether the distinction into different zones – collisional flow and erodible bed – in an inclined plane flow is similar to that between the plug and shear layers for a vertical channel flow. The volume fraction for arrested dynamics is considered to be $\phi _{ad} = 0.587$ in the channel flow for definiteness, as we performed soft-sphere DEM simulations with $\mu _p = 0.5$.

In the vertical channel flow, the simulations are carried out to examine the flow mechanics at the particle scale, to examine which of the different approximations apply in a vertical channel flow. An important issue is whether there are multiple zones with distinct flow regimes, the nature of the rheology in these zones, and whether properties like the volume fraction and the shear rate vary continuously across these zones. Of interest is whether the relation between the friction coefficient and the shear rate varies in different zones, and the nature of the variation. Also of interest is the effect of particle stiffness on rheology in different zones. Another relevant quantity in the kinetic-theory-based models is the granular temperature, which is the mean square of the fluctuating velocity of the particles. We examine whether the variation of the granular temperature is continuous across the channel, and whether the temperature is non-zero in the plug zone. The nature of the velocity profile in the shearing zone is also examined here, to determine the universality and scaling with respect to the channel width. The relation between the slip velocity and the shear rate at the wall, and their variation with channel width, is also determined.

The simulation methodology is summarized briefly in Appendix A, and the flow configuration and preparation protocols are discussed in § 2. Appendix B summarizes the theoretical formulations that are compared with simulations in § 3. The important conclusions are summarized in § 4.

2. Flow configuration

Studies are carried out on the flow through a vertical channel of rectangular cross-section with dimensions $H \times W \times B$, as shown in figure 1. The channel is bounded by two flat frictional vertical walls in the $x$ direction at $x = 0$ and $x = W$, and gravitational acceleration acts in the vertical $y$ direction. Periodic boundary conditions are applied in the vertical $y$ direction and in the $z$ direction in the horizontal plane. Attention is confined to steady fully developed flows that are symmetric about the central plane at $x = W/2$, therefore results are shown only for $0 \leq x \leq W/2$. The width of the channel is varied in the range $60 \, d_p \leq W \leq 150 \, d_p$, and the other two dimensions are $B = 40 \, d_p$ and $H = 60 \, d_p + \Delta H$, where $d_p$ is the nominal particle diameter, and the choice of $\Delta H$ is discussed shortly. The total number of particles is $1.7 \times 10^5$ for a system of width $60 \, d_p$, and $4.1 \times 10^5$ for a system of width $150 \, d_p$. A polydisperse mixture of particles is used, where 30 % of the particles have diameter $0.9 \, d_p$, 40 % have diameter $d_p$, and 30 % have diameter $1.1 \, d_p$.

Figure 1. The configuration and coordinate system for analysing the flow in a vertical channel bounded by two flat frictional walls. The hatched surfaces are the flat frictional walls in the $x$ direction, and periodic boundary conditions are applied in the other two directions.

The preparation protocol is shown in figure 2. Initially, a flat surface is placed at the bottom of the channel. The channel is filled by raining the particles uniformly under the effect of gravity up to a height $H - \Delta H$. After the particles settle, the bottom is removed and the periodic boundary conditions are imposed in the vertical $y$ and spanwise $z$ directions. Here, the value of $\Delta H$ is chosen such that a desired value of the average volume fraction $\bar {\phi }$, the ratio of the total volume of the particles to the volume of the channel, can be obtained. For $W = 60\unicode{x2013}150 \,d_p$, a final steady state is achieved for $\bar {\phi }$ in the range $\bar {\phi }_{cr}$ to $\bar {\phi }_{max}$. When $\bar {\phi }$ is greater than $\bar {\phi }_{max} = 0.62$, the channel jams and there is no flow. When $\bar {\phi }$ is decreased below $\bar {\phi }_{cr}$, there is an oscillatory and then an accelerating flow (Debnath, Kumaran & Rao 2019). The value of $\bar {\phi }_{cr}$ is found to be 0.59 for $W = 100\unicode{x2013}150 \,d_p$. The present study is restricted to the steady fully developed flows.

Figure 2. Preparation protocol for the channel flow. An impenetrable surface is placed at the bottom (a), and particles are filled in up to a height $(H - \Delta H)$. The bottom is removed (b), periodic boundary conditions are applied in the $y$ and $z$ directions, and the flow is evolved to reach steady state.

Raafat, Hulin & Herrmann (1996) studied the flow of sand grains in a capillary tube, where the ratio of the particle size to the tube diameter was relatively small, in the range 6–30. The average volume fraction $\bar {\phi }$ was in between the dilute regime where there is a free fall of grains, and the dense slow flow regime. At these intermediate average volume fractions, density waves are observed in the channel, in the form of bubbles and clogs. For the relative small size of the tube used, finite-size effects are likely to be important. This is different from the present analysis, where the size ratio is in the range 60–150, and $\bar {\phi }$ is in the dense flow regime. Due to this, density waves are not observed here.

The time required for the flow to attain a steady fully developed state is $2 \times 10^3 \sqrt {d_p/g}$, and the averaging time for calculating the properties is $6 \times 10^2 \sqrt {d_p/g}$. To calculate the variation of the properties in the cross-stream $x$ direction, the channel is divided into bins of width $d_p$, except for the bin near the wall, which is set to $1.5 \, d_p$. The velocity, angular velocity and granular temperature fields are calculated as

(2.1) $$\begin{gather} \boldsymbol{v} = \frac{\sum_{i=1}^N m_i \boldsymbol{v}_i}{\sum_{i=1}^N m_i}, \end{gather}$$

(2.2) $$\begin{gather}\boldsymbol{\varOmega} = \frac{\sum_{i=1}^N I_i \boldsymbol{\varOmega}_i}{\sum_{i=1}^N I_i}, \end{gather}$$

(2.3)$$\begin{gather}T = \frac{1}{6N} \sum_{i=1}^N \left(m_i (\boldsymbol{v}_i-\boldsymbol{v})^2 + I_i (\boldsymbol{\varOmega}_i - \boldsymbol{\varOmega})^2\right), \end{gather}$$

where N is the number of particles whose centres lie in the bin, $m_i$ and $I_i = \frac {1}{10} m_i d_i^2$ are the mass and moment of inertia of particle $i$, and $\boldsymbol {v}_i$ and $\boldsymbol {\varOmega }_i$ are the linear and angular velocities of that particle. The stress tensor $\boldsymbol {\sigma }$ is

(2.4)\begin{equation} \boldsymbol{\sigma} = \frac{1}{V} \sum_{i=1}^N \left[ \left( \sum_{j \neq i}^{n_{contacts}} \frac{\boldsymbol{F}_{ij} \,\boldsymbol{x}_{ij}}{2} \right) + m_i (\boldsymbol{v}_i - \boldsymbol{v}) (\boldsymbol{v}_i - \boldsymbol{v}) \right]. \end{equation}

Here, $V$ is the volume of the bin, and the summation $j$ is carried out over all contacts of particle $i$. The first term in (2.4) is the contribution to the stress due to interparticle interactions, and the second term is the kinetic contribution due to the particle fluctuating velocities. The rate of dissipation of energy per unit volume $D$ is calculated as

(2.5)\begin{equation} D = \frac{1}{V \tau} \sum_{i=1}^N \sum_{j \neq i}^{n_{contacts}} \int_0^\tau(\boldsymbol{F}_{ij} \boldsymbol{\cdot} \boldsymbol{v}_{ij}) \,{\rm d}t, \end{equation}

where $\tau$ is the period of the simulation. The rate of dissipation of energy is also a local quantity, evaluated for all the particles with centres within a bin.

The DEM and the simulation parameters are described in Appendix A. In the results, the particle density is considered to be $1$, so that the mass dimensions in all quantities are scaled appropriately.

3. Results

3.1. Velocity and stress profiles

There is a steady flow for a small range of the average volume fraction $\bar {\phi } = 0.59\unicode{x2013} 0.61$. The average velocity fluctuates in time for $\bar {\phi } = 0.62$, and there is no flow when $\bar {\phi }$ exceeds $0.62$. The typical time series for the velocity of the centre of mass $u_{cm}$, scaled by its time-averaged value $\bar {u}_{cm}$, is shown in figure 3. It is observed that there are fluctuations of 1–2 % for $\bar {\phi } = 0.62$. In contrast, for $0.59 \leq \bar {\phi } \leq 0.61$, the variation in the velocity of the centre of mass is less than 1 %, indicating that there is uniform steady flow.

Figure 3. A typical time series of the centre of mass velocity $u_{cm}$ divided by its time-averaged value $\bar {u}_{cm}$ for a channel of width 100 particle diameters and for average volume fraction $\bar {\phi } = 0.62$ ($\circ$), $\bar {\phi } = 0.61$ ($\triangle$), $\bar {\phi } = 0.6$ ($\nabla$) and $\bar {\phi } = 0.59$ ($\diamond$).

The velocity profiles, shown as functions of cross-stream distance in figure 4(a), exhibit a plug-like behaviour in the central region, a shear zone close to the wall and significant slip at the wall. Both the slip velocity $u_{slip}$ and the maximum velocity $u_{max}$ increase as $\bar {\phi }$ decreases, and as the channel width $W$ increases. Both $u_{slip}$ and $u_{max}$ increase by a factor of $5$ or more when $\bar {\phi }$ is decreased from $0.62$ to $0.59$. There is also a significant increase in $u_{max}$ and $u_{slip}$, by a factor of $3$ when $W$ is increased from $60 \, d_p$ to $150 \, d_p$. From figure 4, it is also evident that the length scale for the shearing region at the wall is also a function of $\bar {\phi }$ and $W$.

Figure 4. The velocity $u_y$ as a function of $x/d_p$ (a) and the scaled velocity $(u_y - u_{slip})/(u_{max}-u_{slip})$ as a function of $x/\delta _{0.95}$. Here, $u_{slip}$ is the slip velocity at the wall, $u_{max}$ is the maximum velocity at the centre and $\delta _{0.95}$ is the distance at which $(u_y - u_{slip})/(u_{max}-u_{slip}) = 0.95$. The average volume fractions are $\bar {\phi } = 0.62$ ($\circ$), $\bar {\phi } = 0.61$ ($\triangle$), $\bar {\phi } = 0.6$ ($\nabla$) and $\bar {\phi } = 0.59$ ($\diamond$). The profiles are shown in black for $W = 60 \, d_p$, blue for $W = 100 \, d_p$, green for $W = 120 \, d_p$ and brown for $W = 150 \, d_p$. The red dashed line in (b) is the exponential fit equation (3.1).

In figure 4(b), the scaled velocity is defined as $(u_y - u_{slip})/(u_{max} - u_{slip})$. This is plotted as a function of the scaled distance from the wall, $x/\delta _{0.95}$, where $\delta _{0.95}$ is the distance at which $(u_y - u_{slip})/(u_{max}-u_{slip}) = 0.95$. Plotted in these scaled coordinates, the scaled velocity follows a universal profile that is independent of $\bar {\phi }$ and $W$. Thus the scaled velocity profile has a universal form, though $u_{slip}$, $u_{max}$ and $\delta _{0.95}$ vary significantly with $\bar {\phi }$ and $W$. The red dashed line in figure 4(b) is the fit

(3.1)\begin{equation} \frac{u_y-u_{slip}}{u_{max}-u_{slip}} = 1 - \exp(- 3 x/\delta_{0.95}). \end{equation}

This exponential function provides an excellent fit for the scaled velocity profile for the entire range of $\bar {\phi }$ and $W$ examined here. The average velocity $\bar {u}$ is then expressed in terms of $u_{max}$, $u_{slip}$ and $\delta _{0.95}$ as

(3.2)\begin{equation} \bar{u} = u_{max} - \frac{2 \delta_{0.95} (u_{max}-u_{slip})(1 - \exp(- 3 W /2 \delta_{0.95}))}{3 W}. \end{equation}

There are some granular flows for which the scaled velocity profiles are self-similar, such as the Bagnold profile for the flow down an inclined plane, or the error function profiles in shearing zones in split-bottom geometries (Fenistein & van Hecke 2003). A self-similar velocity profile has not been established in previous studies for the granular flow in a vertical channel, though Pouliquen & Gutfraind (1996) fitted an exponential profile for the velocity in the shearing zone. The exponential scaling in the latter resulted from the assumption that the momentum transport is an activated process, where the probability of yielding depends on the difference between the stress and a yield stress. The shear stress and the velocity are then given by an exponential function.

In figures 5(a) and 5(b), $u_{slip}$ and $u_{max}$ are shown as functions of $\bar {\phi }$ for different values of $W$. Both $u_{slip}$ and $u_{max}$ decrease by about one order of magnitude when $\bar {\phi }$ is increased from $0.59$ to $0.62$. In these figures, $u_{slip}$ and $u_{max}$ scaled by $\sqrt {W}$ are shown by the brown lines referenced to the right ordinate. It is observed that there is an excellent collapse of the data for $u_{slip}$ and $u_{max}$ when scaled by $\sqrt {W}$ for $W \geq 100 \, d_p$, though the collapse is not as good for $W = 60 \, d_p$. This is expected if the Froude number for $u_{slip}$ and $u_{max}$ based on $W$ is a constant. The ratio $u_{max}/u_{slip}$ is shown in figure 5(c). It is remarkable that this ratio is nearly a constant, varying by about 2 % in a very small range 1.37–1.40 for all values of $W$ and $\bar {\phi }$.

Figure 5. Plotted as a function of $\bar {\phi }$: (a) $u_{slip}$ (black lines referenced to left ordinate) and $(u_{slip}/\sqrt {W})$ (brown lines referenced to right ordinate); (b) $u_{max}$ (black lines referenced to left ordinate) and $(u_{max}/\sqrt {W})$ (brown lines referenced to right ordinate); (c) the ratio of $u_{max}$ and $u_{slip}$; and (d) the shear layer thickness $\delta _{0.95}$ (black lines referenced to left ordinate ) and $\delta _{0.95}/W$ (brown lines referenced to right ordinate). The channel widths are $60 \,d_p$ ($\circ$), $100 \,d_p$ ($\triangle$), $120 \,d_p$ ($\nabla$) and $150\, d_p$ ($\diamond$).

The thickness of the shear zone $\delta _{0.95}$ is shown as a function of $\bar {\phi }$ in figure 5(d). This thickness decreases as $\bar {\phi }$ increases and as $W$ decreases, indicating that $\delta _{0.95}$ is dependent on $\bar {\phi }$ and $W$. Shown in brown lines in figure 5(d), referenced to the right ordinate, is the ratio $\delta _{0.95}/W$. This ratio is independent of the channel width for $W \geq 100\,d_p$, but not for $W = 60\,d_p$. This indicates that $\delta _{0.95}$ increases approximately proportional to $W$ for the flow in the vertical channel. A consequence of these scaling relations is that the average flow velocity in (3.2) is proportional to $\sqrt {W}$, provided that $W$ is greater than about $100 \, d_p$. Thus the Froude numbers based on the average velocity, maximum velocity and slip velocity are independent of the channel width, and they depend only on the average volume fraction. This scaling law for the velocity is similar to the Beverloo correlation (Beverloo et al. 1961), which postulates that the average velocity through an orifice of width $D$ is independent of $d_p$ for $d_p \ll D$. From dimensional analysis, the flow velocity has to be scaled by $\sqrt {gD}$.

Cawthorn (2011) has obtained an analytical solution of an incompressible $\mu (I)$ model for the vertical channel flow. The assumptions used are linear variation of $\mu$ with $I$, no variation in $\phi$, and a no-slip boundary condition at the rough walls. This solution predicts that the velocity in the plug is proportional to $W^{3/2}$, and the scaled thickness of the shear layer $\delta _{0.95}/W$ is a constant for high flow rate. In contrast, Mohan et al. (1999) have used an asymptotic analysis for purely rough walls, and their Cosserat model predicts that the shearing zone thickness is proportional to $W^{1/3}$ in the limit $d_p \ll W$. Recently, Barker et al. (2022) have obtained an analytical solution using a compressible $\mu (I)$ based model, and compared it with their DEM results. They have performed simulations in a periodic cell, where the gravitational acceleration is in opposite directions in the two halves of the cell separated by the vertical mid-plane. This results in a zero-velocity condition at the mid-plane and at the periodic boundaries, which is equivalent to a no-slip boundary condition at the wall. Linear approximations for the dependence of $\mu$ and $\phi$ on $I$ are used in the model of Barker et al. (2022); this results in an exponential variation of $\phi$ with distance from the wall. There is a difference between the exponential profiles predicted by the model and simulation results. However, the model prediction for the velocity profile, which is the sum of linear and exponential functions of the distance from the wall, is in reasonable agreement with simulation results. The average velocity and shear layer thickness are proportional to $W^{3/2}$ and $W$, respectively, which are similar to the values in Cawthorn (2011).

The average velocity is proportional to $W^{3/2}$ for the following reason. The maximum shear stress at the wall $\sigma _{xy}$ is proportional to $W$ from the momentum balance equation, if the volume fraction variation in the shearing zones is neglected. If the wall friction coefficient is a constant, then the normal stress $\sigma _{xx}$ is also proportional to $W$. If a constant friction constant is substituted in the linear $\mu (I)$ model, then the shear rate at the wall scales as $\sqrt {\sigma _{xx}/\rho d_p^2} \sim g^{1/2}W^{1/2}/d_p$ for large $I$, where $d_p$ is the particle diameter. The shear layer thickness increases proportional to $W$ in the models of Cawthorn (2011) and Barker et al. (2022). Consequently, the maximum and average velocities are proportional to $g^{1/2} W^{3/2}/d_p$, where $g$ is the acceleration due to gravity. It should be noted that this $W^{3/2}$ scaling can be derived in the limit $I \gg 1$ only if a linear relation is postulated between $\mu$ and $I$. If the nonlinear relationship is used where the friction coefficient asymptotes to constants for large $I$, then the scaling is more complicated.

For fixed channel width, the velocity scaling $g^{1/2} W^{3/2}/d_p$ implies that the maximum velocity increases as the particle diameter decreases, and it diverges as the inverse of the particle diameter for $W/d_p \gg 1$. In contrast, the Beverloo correlation (Beverloo et al. 1961) is based on the assumption that the particle diameter $d_p$ is not a relevant length scale when it is small compared to the orifice diameter. If the channel width is substituted for the orifice diameter in the Beverloo correlation, then the velocity scaling $(g W)^{1/2}$ is obtained in the limit $W \gg d_p$. The Beverloo scaling is observed in the present study, where smooth frictional walls are used and there is significant slip at the wall. If our results are found to be robust, then they imply that the nature of the wall may not change the scaling for $\delta _{0.95}$, but can change qualitatively the scaling of the velocity and the flow rate.

Another issue that needs to be examined is the ratio of the wall roughness scale and the particle diameter. A standard procedure for generating rough walls is to use frozen spherical particles at the wall, and the wall particle diameter is usually considered equal to that for the flowing particles. In Barker et al. (2022), periodic boundary conditions and a reverse flow were used for generating rough walls. One future avenue of research is to transition from rough to smooth walls by decreasing successively the wall particle diameter relative to the flowing particle diameter, and examining the scaling of the average velocity separately with the flowing and wall particle diameters (Kumaran & Bharathraj 2013; Bharathraj & Kumaran 2017). This would provide some insight into whether the wall roughness or the flowing particle diameter is the relevant length scale for the shear rate at the wall.

The slip velocity $u_{slip}$ and $u_{max}-u_{slip}$ are proportional to $\sqrt {W}$, and $\delta _{0.95}$ is proportional to $W$, and the ratio $u_{slip}/({\rm d} u_y/{{\rm d}\kern0.7pt x})$ at the wall is proportional to $W$. Therefore, if the slip boundary condition ${\rm d} u_y/{{\rm d}\kern0.7pt x} = u_{slip}/l_s$ is used at the wall, and the velocity gradient is calculated from (3.1), then the slip length $l_s$ increases proportional to $W$. This is different from earlier studies of Mohan et al. (2002), where the slip length was considered a constant.

The slip velocity defined in Shojaaee et al. (2012) is the relative velocity at the point of contact between a particle and the wall. The contact velocity of the particle surface includes the translational velocity of the particle centre and the cross-product of the rotation velocity and the vector distance between the point of contact and particle centre. However, in the current study, the slip velocity is assumed to be the velocity of the interval closest to the wall, not the surface velocity of the particle at contact. The particle angular velocity is also reported separately, and it is found that the angular velocity deviates from one half of the vorticity close to the wall.

The normal and shear stress profiles are shown in figure 6. As required by the momentum balance, the normal stress in the cross-stream direction, $\sigma _{xx}$, is a constant across the channel. This constant value increases significantly as $\bar {\phi }$ increases. The shear stress $\sigma _{xy}$ is close to a linear function of distance from the centre; there is departure from linearity in the shear zone due to the variation in $\phi$. There is very little change in the slope of $\sigma _{xy}$ as $\bar {\phi }$ is changed, because the variation in $\bar {\phi }$ is less than 5 %.

Figure 6. The normal stress $\sigma _{xx}$ (black, left vertical axis) and shear stress $\sigma _{xy}$ (brown, right vertical axis) in the flow of a granular material in a vertical chute of width $100 \, d_p$. The average volume fractions are $\bar {\phi } = 0.62$ ($\circ$), $\bar {\phi } = 0.61$ ($\triangle$), $\bar {\phi } = 0.6$ ($\nabla$) and $\bar {\phi } = 0.59$ ($\diamond$). The black vertical lines are boundaries between the plug and the dense shearing zones, the blue vertical lines are the boundaries between dense and loose shearing zones, and the red vertical line is the boundary between the loose shearing zone and the wall shearing zone.

The cross-stream normal stress $\sigma _{xx}$ is shown as a function of $\bar {\phi }$ for different channel widths $W$ in figure 7(a). There is a significant increase in $\sigma _{xx}$ as $W$ is increased, and as $\bar {\phi }$ is increased. However, the ratio of the normal stress to the channel width, $(\sigma _{xx}/W)$, shown by the brown lines referenced to the right ordinate in figure 7(a), varies very little as $W$ is varied for $W > 60 \, d_p$. This implies that $\sigma _{xx}$ increases approximately proportional to $W$. To put this in perspective, recall that $\sigma _{xy}$ at the wall increases proportional to $W$ if the decrease in $\phi$ close to the wall is neglected from the momentum balance. This is in agreement with the finding of Barker et al. (2022). When both $\sigma _{xx}$ and $\sigma _{xy}$ at the wall are proportional to $W$, the friction coefficient $\mu _{wall}$ at the wall should be independent of $W$. Figure 7(b) shows that this is indeed the case. The wall friction coefficient is independent of $W$, but it decreases as $\bar {\phi }$ increases. The latter trend may be understood by noting that small changes in $\bar {\phi }$ affect $\sigma _{xy}$ less strongly than $\sigma _{xx}$. Thus the value of $\sigma _{xx}$ across the channel is determined from the wall shear stress required to balance the weight of the material per unit area of the wall and the wall friction coefficient that is a function of $\bar {\phi }$ for specific wall and particle properties.

Figure 7. (a) The cross-stream normal stress $\sigma _{xx}$ as a function of $\bar {\phi }$ is shown by the black lines referenced to the left ordinate, and the ratio $\sigma _{xx}/W$ is shown by the brown lines referenced to the right ordinate. (b) The wall friction coefficient $\mu _{wall}$ as a function of $\bar {\phi }$. The channel widths are $60 \, d_p$ ($\circ$), $100 \, d_p$ ($\triangle$), $120 \, d_p$ ($\nabla$) and $150 \, d_p$ ($\diamond$).

3.2. Comparison with kinetic models

The volume fraction $\phi$ as a function of cross-stream distance $x/d_p$ is shown in figure 8(a) for channel width $W = 100 \, d_p$. The region at the centre of the channel is the ‘plug’ zone where $\phi$ is a constant. The boundary of the plug zone is indicated by the black vertical lines for the four different values of $\bar {\phi }$. It is interesting to note that $\phi$ in the plug is ${\sim }0.63$, which is discernibly lower than the random close packing volume fraction $\phi _{rcp} = 0.64$ for a monodisperse system. This indicates that the plug zone is not a jammed assembly of particles, but there is particle agitation even in this plug zone, and the granular temperature $T$ in this region is measurable, as shown in figure 8(c). When $\bar {\phi }$ is decreased, the volume fraction in the plug zone is nearly a constant, but the thickness of the plug zone decreases. At the lowest average volume fraction $\bar {\phi } = 0.59$, the plug zone occupies less than half of the channel width $W$. It should be noted that there is no layering in this region even at $\bar {\phi } = 0.62$. It is also verified that there is no crystallization; the $Q_6$ order parameter (Kumaran 2009c,d) in this region is below 0.4.

Figure 8. (a) The volume fraction $\phi$, (b) the shear rate ${\rm d}u_y/{\rm d}\kern0.7pt x$, (c) the granular temperature $T$, and (d) the rates of shear production of energy (black) and dissipation (red) due to particle interactions for $W = 100 \, d_p$. The average volume fractions are $\bar {\phi } = 0.62$ ($\circ$), $\bar {\phi } = 0.61$ ($\triangle$), $\bar {\phi } = 0.6$ ($\nabla$) and $\bar {\phi } = 0.59$ ($\diamond$). The black vertical lines are boundaries between the plug and dense shearing zones; the blue vertical lines are boundaries between the dense and loose shearing zones; and the red line is a boundary between the loose and wall shearing zones.

There is a transition from the plug zone to a ‘dense shearing’ zone when $\phi$ decreases below 0.63 at the black vertical lines in figure 8(a). A change in the slope of $\phi$ is evident at the boundary of the plug and dense shearing zones, and $\phi$ decreases continuously in the dense shearing zone. The blue vertical lines are drawn in figure 8(a) where $\phi = \phi _{ad} = 0.587$, where $\phi _{ad}$ is the volume fraction for arrested dynamics. This is the boundary between the dense and ‘loose shearing’ zones. Very close to the wall, within a distance of about two particle diameters, there is a wall shearing layer with a steep decrease in $\phi$, marked by a red vertical line. The distinctions between the plug and dense shearing zones (marked by black vertical lines) and the loose and wall shearing zones (marked by a red vertical line) are clearer in figure 9, where clear breaks are seen in the profiles of the ratios of the angular velocity to shear rate, and the shear rate to square root of the granular temperature $\sqrt {T}$. For the highest average volume fraction $\bar {\phi } = 0.62$, there appears to be no loose shearing zone as $\phi$ decreases to $\phi _{ad} = 0.587$ in the wall shearing zone. For lower values of $\bar {\phi }$, there is a clear distinction between the dense and loose shearing zones at $\phi = \phi _{ad}$, and the loose and wall shearing zones.

Figure 9. The ratios of (a) $\varOmega _z$ to $\omega _z$, and (b) $({\rm d} u_y/{{\rm d}\kern0.7pt x})/\sqrt {T}$ versus $x/d_p$, for $W = 100\,d_p$ and $\bar {\phi } = 0.62$ ($\circ$), $\bar {\phi } = 0.61$ ($\triangle$), $\bar {\phi } = 0.6$ ($\nabla$) and $\bar {\phi } = 0.59$ ($\diamond$).

The velocity profiles shown in figure 4 do not exhibit any perceptible change between the plug and shearing zones. However, the shear rate shown in figure 8(b) exhibits a clear distinction between the dense shearing zone (between the black and blue vertical lines) where there is a monotonic decrease in the shear rate, and the plug zone to the right of the black vertical lines. In the plug zone, the standard deviations in the measurements, shown by the error bars, are comparable to the mean values, indicating that the shear rate is zero to within the simulation accuracy. In contrast, in the dense shearing zone, the error bars are smaller than the mean values, and the profiles are continuous, even when the shear rate is 4–5 orders of magnitude smaller than that at the wall. There is also a small but distinct change in the curvature of the shear rate at the boundary (blue vertical lines) between the dense and loose shearing zones. Whereas the shear rate profiles have negative curvature in the dense shearing zone, they have positive curvature in the loose shearing zone. This change in curvature is significant as it implies a change in the rheology of the flow as $\phi \rightarrow \phi _{ad}$. From the momentum balance, the shear stress is approximately a linear function of distance from the centre, if we neglect the variations in $\phi$. If the shear rate profile is concave upwards, then the fluid is shear thickening with power-law exponent greater than $1$, whereas if it is concave downwards, then it is shear thinning with power-law exponent less than $1$. Therefore, an inflection point in the shear rate curve at $\phi = \phi _{ad}$ in figure 8(b) implies a transition from shear thickening to shear thinning behaviour. Similarly, the plane shear simulations in Chialvo et al. (2012) have captured a significant change in rheology, where the flow is quasi-static in the limit of very small shear rate for $\phi > \phi _{ad}$, and inertial for $\phi < \phi _{ad}$, and an intermediate regime is observed for $\phi = \phi _{ad}$.

The granular temperature $T$ profiles shown in figure 8(c) are continuous all the way from the wall to the centre of the channel. The standard deviation is always much smaller than the mean value, even at the centre where the temperature is more than five orders of magnitude smaller than that at the wall. The non-zero granular temperature indicates that the particles are not jammed, but they do have fluctuating motion. This is consistent with the $\phi$ profile in the plug zone in figure 8(a), where $\phi$ is lower than $\phi _{rcp}$. In the shearing zone, the granular temperature decreases as $\bar {\phi }$ increases. In contrast, in the plug zone, the granular temperature appears independent of $\bar {\phi }$, with the exception of $\bar {\phi } = 0.62$ where the granular temperature is noticeably higher. This could be attributed to the significantly higher normal stress $\sigma _{xx}$ shown in figure 6(a) – this higher normal stress seems to be a consequence of higher agitation of the particles in the plug zone.

Figure 10. The variation of (a) $\phi$, (b) $u_y$, (c) ${\rm d}u_y/{{\rm d}\kern0.7pt x}$ and (d) $T$, for $k_n = 10^5$ ($\circ$), $k_n = 10^6$ ($\triangle$), $k_n = 10^7$ ($\nabla$) and $k_n = 10^8$ ($\diamond$), and $\bar {\phi } = 0.61$ (black lines) and $\bar {\phi } = 0.59$ (brown lines). The red line is the boundary between the wall and loose shearing zones; the black and brown vertical lines from left to right are the boundaries between the loose and dense shearing zones, and the dense shearing zone and plug zone, respectively.

Figure 11. (a) The friction coefficient $\mu = |\sigma _{xy}/\sigma _{xx}|$ as a function of the inertial number $I$, and (b) $0.64 - \phi$ as a function of $I$, for $W = 100 \, d_p$ and $\bar {\phi } = 0.62$ ($\circ$), $\bar {\phi } = 0.61$ ($\triangle$), $\bar {\phi } = 0.6$ ($\nabla$) and $\bar {\phi } = 0.59$ ($\diamond$). The data in the dense shearing zone are shown in brown lines, the loose shearing zone by black lines and the wall shearing zone by red lines. The dotted blue lines are (B3) and (B4), respectively, and the inclined black dotted line in (b) shows a slope of $1$.

The rates of shear energy production (black) and dissipation due to particle interactions (red) are shown in figure 8(d). The energy production rate is significantly higher than the energy dissipation rate in the dense and loose shearing zones. The rates of production and dissipation of energy decrease as $\bar {\phi }$ is increased, exhibiting the same trend as the granular temperature in the shearing zones. In the plug zone, there is virtually no shear production of energy because the shear rate is zero to within the simulation resolution. However, the rate of dissipation of energy is clearly measurable and approximately a constant in this region. Since there is no shearing in the plug zone, there is no shear production of energy. The temperature in the plug zone is determined from a balance between the rate of conduction of energy and the rate of dissipation due to particle interactions. The dashed lines in figure 8(c) are fits using the conduction equation (B13), with specified values of the temperature and zero flux at the centre. The best fit for the conduction length in these curves is $l_c \approx 11.5\,d_p$. These fits provide good predictions for the temperature in the plug zone, and the temperature profile departs from these fits in the dense shearing zone where there is energy production due to shear.

Two other kinematic measures are shown in figure 9. The particle angular velocity is a non-conserved variable in a fixed coordinate reference frame for the flow of frictional particles. The angular velocity $\varOmega _z$ scaled by ${\rm d}u_y/{{\rm d}\kern0.7pt x}$ is shown as a function of $x/d_p$ in figure 9(a). For the present unidirectional flow, the vorticity is equal to $\omega _z = {\rm d} u_y/{{\rm d}\kern0.7pt x}$. When the particle rotation rate and the local material rotation rate are the same, the ratio $\varOmega _z/\omega _z$ equals $\frac {1}{2}$. However, $\varOmega _z$ could be very different from half of the vorticity in regions of thickness comparable to the particle diameter at the boundaries of the flow and in a region where the material rotation rate is negligible. In figure 9(a), there is a lot of scatter in the data in the plug zone; this is because the magnitude of the angular velocity and the shear rate are close to the simulation resolution in this region, and small errors in either of these quantities result in large variations in the ratio. However, this ratio is close to $\frac {1}{2}$ in both the dense and loose shearing zones, indicating that the particles are rotating with the same angular velocity as the material in these zones, as expected from the continuum approximation for a material with no microscopic torques. The ratio differs from $\frac {1}{2}$ in the wall shearing zone as well, due to the effect of wall collisions.

The ratio ${d}_p({\rm d} u_y/{{\rm d}\kern0.7pt x})/\sqrt {T}$ is shown in figure 9(b). This ratio has a behaviour similar to that for the shear rate, but it seems to have a positive slope at the wall, in contrast to the shear rate, which has a negative slope at the wall. This indicates that the rheology in the wall shearing zone is different from that in the loose shearing zone, and a continuum approximation may be difficult to formulate for this region as it extends over a distance equal to only two particle diameters. This ratio is continuous across the loose and dense shearing zones, though it is shown later, in figure 13, that the dependence on $\phi$ is very different.

The effect of particle stiffness $k_n$ on the volume fraction $\phi$, velocity $u_y$, shear rate ${\rm d}u_y/{{\rm d}\kern0.7pt x}$ and temperature $T$ are shown in figure 10 for channel width $W = 100 \, d_p$ and for two different average volume fractions $\bar {\phi }=0.61$ and $0.59$. In all cases, $\phi$ and $u_y$ in figures 10(a) and 10(b) are independent of $k_n$ for $k_n \geq 10^5$, except for one case corresponding to variation of $\phi$ with $x/d_p$ for $k_n = 10^5$ and $\bar {\phi } = 0.59$. The profiles of ${\rm d}u_y/{{\rm d}\kern0.7pt x}$ and $T$ in figures 10(c) and 10(d) suggest a more complicated picture. These are independent of $k_n$ in the wall and loose shearing zones for $\phi < \phi _{ad} = 0.587$. However, there is some dependence of ${\rm d}u_y/{{\rm d}\kern0.7pt x}$ on $k_n$, and an even stronger dependence of $T$ on $k_n$ in the dense shearing and plug zones. Here, $T$ in the plug zone decreases as $k_n$ increases, and it appears to depend only on $k_n$ and not on $\bar {\phi }$ in this region. This is expected because $\phi$ in this region is approximately 0.63 for all $\bar {\phi }$, and there is no shearing in this region. There are enduring contacts between particles, and the divergence of the shear stress is balanced by the gravitational force density due to the contacts. The monotonic decrease of $T$ with $k_n$ suggests that $T$ will decrease to zero in the hard-particle limit, and the plug zone will be in a jammed state. When $k_n$ is finite, there is agitation within the plug zone resulting in a non-zero temperature. As $T$ in the dense shearing zone also depends on $k_n$, the dense shearing zone cannot be described using hard-particle models. However, in the loose shearing zone and the wall shearing zone, the shear rate and granular temperature are independent of the particle stiffness, hence these zones can be described using the hard-particle models.

The friction coefficient $\mu = \sigma _{xy}/\sigma _{xx}$ is shown as a function of the dimensionless ‘inertial number’ $I = d_p |{\rm d} u_y/{{\rm d}\kern0.7pt x}|/\sqrt {\sigma _{xx}/\rho _p}$ in figure 11(a). In this figure, $\mu$ in the dense shearing zone is shown by the brown lines, the loose shearing zone by the black lines and the wall shearing zone by the red lines. A log scale is used on the horizontal axis, in order to highlight features that are not visible readily in a linear scale. For $\bar {\phi } = 0.62$, $\mu$ is lower than for $0.59 \leq \bar {\phi } \leq 0.61$. For the latter, $\mu$ approaches a finite value for the lowest values of $I$ discernible in the dense shearing zone. It scales approximately as a linear function of $\log (I)$ as $I$ becomes small, and is independent of $\bar {\phi }$ for $I \lesssim 10^{-3}$ in the dense shearing zone. However, it exhibits a definite dependence on $\bar {\phi }$ for $I \gtrsim 10^{-3}$ in the loose shearing zone and the wall shearing zone. The qualitative behaviour of the $\mu$–$I$ curves is the same for other channel widths, and is not shown. The inertial number varies over a large range ($O(10^{-5})\unicode{x2013}O(10^0)$) from the centre to the wall in the flow through a vertical channel, and the linear relation between $\mu$ and $I$ does not apply to the whole domain, in contrast to the assumption of Barker et al. (2022).

The scaling of $\phi$ with $I$, shown in figure 11(b), also reveals interesting features. Instead of plotting $\phi$ itself, which results in a rather featureless graph, we have shown the difference between $\phi _{rcp} = 0.64$ and $\phi$. It depends only on $I$ in the loose shearing zone for $\phi < \phi _{ad} = 0.587$, and depends on I and $\bar{\phi}$ in the dense shearing zone for $\phi \geq \phi _{ad}$. A clear break in the slope of the graph is visible at $\phi = \phi _{ad}$, which is the boundary between the dense and loose shearing zones. The slope of the $\phi$–$I$ curves is larger in the loose shearing zone in comparison to the dense shearing zone.

Figure 12. The variation of $-\log (\phi )$ with $x/d_p$ for $W = 100 \, d_p$ (solid lines) and $W = 150 \, d_p$ (dashed lines), and $\bar {\phi } = 0.62$ ($\circ$), $\bar {\phi } = 0.61$ ($\triangle$), $\bar {\phi } = 0.6$ ($\nabla$) and $\bar {\phi } = 0.59$ ($\diamond$).

The dotted blue lines in figures 11(a) and 11(b) are the predictions of (B3) and (B4), with the parameters $\mu _s = 0.384$, $\mu _\infty = 0.625$, $I_0 = 0.3$, $\phi _{min}=0.4$ and $\phi _{max} = 0.6$ reported in Jop (2015). It is evident that there is a quantitative difference between the simulation results and the parameters commonly used in the $\mu$–$I$ models. More importantly, there is a difference in the shape of the curves, which is not captured adequately when drawn on a linear axis. The simulation results for $\mu$ exhibit a slow logarithmic dependence on $I$ for the lowest values accessible in the dense shearing zone, in contrast to the constant value predicted by the $\mu$–$I$ model. Similarly, the volume fraction $\phi$ in the dense shearing zone increases to 0.63 in the limit of low $I$, but (B4) predicts a constant value of 0.6 in the limit $I \rightarrow 0$. The decrease in the volume fraction below $\phi _{ad}$ in the loose shearing zone is captured reasonably by (B4).

Figure 13. Plots of (a) $\sigma _{xx}/T$, (b) $\sigma _{xy}/(|{\rm d} u_y/{{\rm d}\kern0.7pt x}|\,T^{1/2})$ and (c) $D/T^{3/2}$ as functions of $(0.64 - \phi )$ for $\bar {\phi } = 0.62$ ($\circ$), $\bar {\phi } = 0.61$ ($\triangle$), $\bar {\phi } = 0.6$ ($\nabla$) and $\bar {\phi } = 0.59$ ($\diamond$), and $W = 100 \, d_p$. The data in the dense shearing zone are shown in black, in the loose shearing zone in brown and in the wall shearing zone in red. The blue circles, from left to right, are the results for the flow down an inclined plane for angles of inclination $20^\circ \unicode{x2013}24^\circ$. The blue dotted lines fit to the data diverging at $\phi = \phi _{ad} = 0.59 \pm 0.002$.

The comparison of our data with that of Barker et al. (2022) should be made with care, because the results of the latter are equivalent to the flow through a vertical channel with purely rough walls that enforce a no-slip boundary condition. The wall shearing zone observed in the current study for the flat frictional walls may not be present if the walls are rough. Additionally, the linear approximations used in Barker et al. (2022) to obtain an analytical solution may not be valid in the whole domain. For example, the linear relation between $\mu$ and $I$ is not valid in the loose and wall shearing zones where $\phi < \phi _{ad}$, and the linear relation between $\phi$ and $I$ does not capture the dense shearing zone where $\phi \geq \phi _{ad}$. Barker et al. (2022) obtained an exponential variation of $\phi$ with the distance from the wall in their analytical solution. However, a clear linear variation is not observed in the shearing zone when $\log (\phi )$ is plotted as a function of $x/d_p$, even if the first two or three points in the wall shearing zone are excluded, as shown in figure 12. The velocity profile obtained in the analytical solution is the combination of a linear function and an exponential function in Barker et al. (2022); it does not fit our data as well as the exponential function (3.1) (see figure 4b). Hence the fits of the functional forms for the volume fraction and velocity profiles are not shown here. In Debnath et al. (2022), the nonlinear expression (B3) for $\mu (I)$ was used to solve the models of Barker et al. (2017) and Schaeffer et al. (2019) numerically. However, this does not result in better agreement between the DEM results and the model predictions for $\mu (I)$ and $\phi (I)$.

Figure 13 shows different quantities of interest in the kinetic-theory-based models as a function of the volume fraction $\phi$. Here, the focus is on $\phi$ near $\phi _{rcp} = 0.64$, so the data are plotted as a function of $0.64 - \phi$ on a log-log scale. In this figure, the values in the dense shearing zones are shown by the black lines, the loose shearing zones by the brown lines, and the wall shearing zones by the red lines. The constitutive relations (B8), (B9) and (B10) are examined in figures 13(a), 13(b) and 13(c), respectively. If these constitutive relations are valid, then the functions $p_\phi = \sigma _{xx}/T$ in figure 13(a), $\mu _\phi = \sigma _{xy}/T^{1/2} ({\rm d} u_y/{{\rm d}\kern0.7pt x})$ in figure 13(b), and $D_\phi = D / T^{3/2}$ in figure 13(c) should depend only on the volume fraction $\phi$. In all cases, it is observed that the dependence of these quantities is different in the dense shearing zone in comparison to the loose shearing zone. There is significant scatter in the dense shearing zone; however, there is less scatter in the loose shearing zone, and the quantities are approximately functions of $\phi$ only in this zone. This indicates that the kinetic-theory-based models, which assume that these quantities depend only on $\phi$, do not apply in the dense shearing zone for $\phi \geq \phi _{ad}$, but they do seem to apply in the loose shearing zone for $\phi < \phi _{ad}$. This supports the conclusion that the flow in the loose shearing zone with $\phi$ less than $\phi _{ad} = 0.587$ can be described by a hard-particle model, and it cannot be used for the dense shearing zone.

The blue circles in figure 13 are the results for the dense granular flow down an inclined plane with a free surface. As explained in the Introduction, the highest volume fraction that can be attained in this unconfined flow is around $\phi = \phi _{ad}$. The volume fraction $\phi$ is a constant in the bulk due to the momentum conservation condition, and it has been found that the quantities in figure 13 depend only on $\phi$. The flow is stable for a relatively small range of angles of inclination, from about $20^\circ$ to about $24^\circ$. The results for this range of angles, separated by $1^\circ$ intervals, are shown by the five blue dots in each panel of figure 13. The results for the inclined plane flows are in agreement with the results for the loose shearing zone for $\phi < \phi _{ad}$, indicating that the dynamics in the loose shearing zone is the same as that in a dense granular flow down an inclined plane.

The coefficients $p_{\phi }$, $\mu _{\phi }$ and $D_{\phi }$ for an inclined plane flow shown by the blue circles in figure 13 diverge at $\phi = \phi _{ad}$ of about $0.585 \pm 0.005$. The value $\phi _{ad} = 0.587$ used in the current study to mark the boundary between dense and loose shearing zones is slightly lower, by less than 1 %, than the value $0.592$ obtained by an extrapolation of the hard-particle simulations (Kumaran 2009c). In figure 13, the divergence, shown by blue dotted curves proportional to $\chi = (\phi _{ad}-\phi )^{-1}$, provides the best fit for the loose shearing zone at $\phi _{ad} = 0.59 \pm 0.002$; the agreement is quite good for $p_\phi$ and $\mu _\phi$, but not as good for $D_\phi$. However, it is clear that the divergences in the stresses and the dissipation rate occur at nearly the same volume fraction $\phi = \phi _{ad}$ for the inclined plane flow and the loose shearing zone. For the inclined plane flow that is unconfined, the divergence in these functions, $p_\phi, \mu _\phi$ and $D_\phi$, leads to flow cessation. For the flow in the vertical channel, as $\phi$ increases beyond $\phi _{ad}$, the divergence in these functions is cut off by a cross-over to the dense shearing zone, where the hard-particle model does not apply.

The results for $p_\phi$, $\mu _\phi$ and $D_\phi$ are compared with the model of Berzi et al. (2020) in figure 14. The qualitative trends of these functions are captured in the loose shearing zone and the dense shearing zone; however, the quantitative differences are evident. The order of magnitude differs in some cases. Their model underestimates the coefficients in the loose shearing zone for $\phi < \phi _{ad} = 0.587$, though the slopes of the curves are approximately the same. There is less disagreement in the dense shearing zone except for $D_{\phi }$. The quantitative difference could be because the functions of Berzi et al. (2020) were fitted for plane shear simulation and an inclined plane flow in Berzi & Jenkins (2015), and not for a vertical chute. Another issue related to the conduction flux in their model is the following. In the channel flow, if the conduction flux at the interface where $\phi _{ad} = 0.587$ has to be matched, then the slope of the $\phi$ profile appears to be discontinuous, which does not agree with the trend predicted by the DEM results (Debnath et al. 2022). More work is required to resolve these issues.

Figure 14. A comparison of the model of Berzi et al. (2020) (red lines) and the results of the present simulations (black lines) for the functions (a) $\sigma _{xx}/T$, (b) $\sigma _{xy}/(|{\rm d} u_y/{{\rm d}\kern0.7pt x}|\,T^{1/2})$ and (c) $(D/T^{3/2})$, as functions of $0.64 - \phi$ for $\bar {\phi } = 0.61$ ($\triangle$) and $\bar {\phi } = 0.59$ ($\diamond$), and $W = 100 \, d_p$.

4. Conclusions

One of the simplest granular flows that could be envisaged is the flow in a vertical rectangular channel bounded by walls in one direction and with periodic boundary conditions in the other two directions. Even this flow exhibits complex structure–rheology relationships. Previous studies have classified the flow broadly into a plug zone at the centre and a shearing zone near the wall.

The important conclusions of the present study are as follows.

1. (i) The flow is studied for different channel widths $W=60\unicode{x2013}150\,d_p$. There is a steady flow for a relatively small range of the average volume fractions $\bar {\phi }$, from about $0.62$ to $0.59$ for $W \geq 100 \,d_p$. The material is static for $\bar {\phi } > 0.62$, and there is no steady state for $\bar {\phi } < 0.59$.

2. (ii) The velocity profile exhibits significant slip at the wall. Both the slip velocity $u_{slip}$ and the maximum velocity $u_{max}$ are proportional to $\sqrt {W}$ for $W \geq 100 \,d_p$. The ratio of the maximum and slip velocities, $u_{max}/u_{slip}$, is found to be nearly a constant in the range 1.38–1.4 for the entire range of channel widths and average volume fractions.

3. (iii) The width of the shear layer, $\delta _{0.95}$, is defined as the distance from the wall at which $u_y - u_{slip}$ is 95 % of its maximum value. A specific exponential form, (3.1), is found to provide an excellent fit for all of the profiles for different values of $W$ and $\bar {\phi }$. The width $\delta _{0.95}$ varies with $\bar {\phi }$, and it increases proportional to $W$ for $W \geq 100 \,d_p$.

4. (iv) The scaling of the velocity and boundary layer thickness with $W$ does not seem to apply for the smallest width $W = 60 \,d_p$. The breakdown of universality at $W = 60\, d_p$ is likely due to the finite-size effects. If the ratio of the width and particle diameter is not sufficiently large, there are two length scales, the particle diameter and the channel width, hence properties may not depend only on the channel width.

5. (v) A detailed analysis of the flow dynamics reveals the presence of multiple zones. In the ‘plug’ zone at the centre of the channel, the volume fraction $\phi$ has a constant value 0.63. The velocity is a constant and the shear rate is zero in this region, to within the simulation accuracy. However, the granular temperature is non-zero, though it could be 4–5 orders of magnitude smaller than that at the wall.

6. (vi) In the plug zone, the granular temperature decreases as the particle stiffness increases. This suggests that the system could be jammed in the limit of infinite stiffness, and a hard-particle approximation cannot be used to describe this region, subject to the following caveats.

7. (vii) Adjacent to the plug zone is a ‘dense shearing’ zone, where the particle volume fraction exceeds $\phi _{ad}$. Here, the shear rate and granular temperature are 3–4 orders of magnitude smaller than the values at the wall. However, the shear rate and the granular temperature are both finite and measurable, and the rate of production of energy due to shear exceeds the rate of dissipation due to particle interactions. The particle angular velocity in this region is equal to half of the vorticity, indicating that the particle and material rotation rates are equal. The shear rate and granular temperature in this region depend on the stiffness of the particles.

8. (viii) In the loose shearing zone between the dense and wall shearing zones, the volume fraction $\phi$ is below $\phi _{ad}$. The shear rate and the granular temperature in this region are independent of the particle stiffness, indicating that the hard-particle models can be used to describe the flow dynamics in this region. It is found that (B8) and (B9) for the stresses, and (B10) for the energy dissipation rate, formulated from the kinetic-theory-based models, are applicable in the loose shearing zone. The values of the coefficients in these constitutive relations are in agreement with those obtained from studies on the granular flow down an inclined plane.

9. (ix) The coefficients $p_\phi$, $\mu _\phi$ and $D_\phi$ in (B8)–(B10) initially diverge as the volume fraction $\phi$ is increased beyond $\phi _{ad}$ in a manner similar to the values for the flow down an inclined plane. However, a transition to the dense shearing zone cuts off this divergence in the vertical channel flow. These functions are in qualitative agreement with those proposed by Berzi et al. (2020), though there are quantitative differences, possibly because the functions of Berzi et al. (2020) were formulated for an inclined plane flow.

10. (x) The simulation results are compared with the predictions of the $\mu$–$I$ model, which relates the friction coefficient $\mu = |\sigma _{xy}/\sigma _{xx}|$ to the inertial number $I$ by (B3), and the volume fraction to $I$ by (B4). The model predictions are different qualitatively from the simulation results in the limit of low $I$. Whereas (B3) predicts that $\mu$ approaches a constant value for $I \rightarrow 0$, the simulation results in figure 11 show a logarithmic decrease down to the lowest value of $I$ accessible in the dense shearing zone. The friction coefficient can be approximated as a function of $I$ in the dense shearing zone, but it varies with the average volume fraction also in the loose shearing zone. The volume fraction $\phi$ is quite well predicted by (B4) in the loose shearing zone, but there is a difference between the model prediction and the simulation results in the dense shearing zone.