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 Reference Fenistein and van Hecke2003). A self-similar velocity profile has not been established in previous studies for the granular flow in a vertical channel, though Pouliquen & Gutfraind (Reference Pouliquen and Gutfraind1996) 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. Reference Beverloo, Leniger and Van de Velde1961), 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 (Reference Cawthorn2011) 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. (Reference Mohan, Nott and Rao1999) 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. (Reference Barker, Zhu and Sun2022) 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. (Reference Barker, Zhu and Sun2022); 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 (Reference Cawthorn2011).

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 (Reference Cawthorn2011) and Barker et al. (Reference Barker, Zhu and Sun2022). 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. Reference Beverloo, Leniger and Van de Velde1961) 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. (Reference Barker, Zhu and Sun2022), 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 Reference Kumaran and Bharathraj2013; Bharathraj & Kumaran Reference Bharathraj and Kumaran2017). 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. (Reference Mohan, Rao and Nott2002), where the slip length was considered a constant.

The slip velocity defined in Shojaaee et al. (Reference Shojaaee, Brendel, Török and Wolf2012) 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. (Reference Barker, Zhu and Sun2022). 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 Reference Kumaran2009c,Reference Kumarand) 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. (Reference Chialvo, Sun and Sundaresan2012) 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. (Reference Barker, Zhu and Sun2022).

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 (Reference Jop2015). 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. (Reference Barker, Zhu and Sun2022) 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. (Reference Barker, Zhu and Sun2022) 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. (Reference Barker, Zhu and Sun2022) 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. (Reference Barker, Zhu and Sun2022); 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. (Reference Debnath, Kumaran and Rao2022), the nonlinear expression (B3) for $\mu (I)$ was used to solve the models of Barker et al. (Reference Barker, Schaeffer, Shearer and Gray2017) and Schaeffer et al. (Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019) 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 Reference Kumaran2009c). 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. (Reference Berzi, Jenkins and Richard2020) 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. (Reference Berzi, Jenkins and Richard2020) were fitted for plane shear simulation and an inclined plane flow in Berzi & Jenkins (Reference Berzi and Jenkins2015), 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. Reference Debnath, Kumaran and Rao2022). More work is required to resolve these issues.

Figure 14.

A comparison of the model of Berzi et al. (Reference Berzi, Jenkins and Richard2020) (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$.