## Abstract

Quasi-two-dimensional bounded heap flow is a useful model for many granular flows in industry and nature. It belongs to a family of free surface flows—inclined chute flow, rotating tumbler flow and unbounded heap flow—but differs from the others in that uniform deposition of particles onto the static bed results in the uniform rise of the heap. The kinematics, however, are only partially understood. We performed discrete element method simulations to study granular flows in quasi-two-dimensional bounded heaps. The experimentally validated computational results show a universal functional form for the streamwise velocity profile for both monodisperse and bidisperse systems when velocities and coordinates are scaled by the local surface velocity and the local flowing layer thickness. This holds true regardless of streamwise location, feed rate, particle size distribution and, most surprisingly, the local particle concentration for bidisperse flows. The local surface velocity decreases linearly in the streamwise direction, while the flowing layer thickness remains nearly constant; both quantities depending only on local flow rate and local mean particle diameter. Additionally, the velocity profile normal to the overall flow, which is important in understanding segregation, can be predicted analytically from the streamwise velocity and matches the simulation results.

## 1. Introduction

Understanding the flow of non-cohesive particles, common in many industrial and geophysical situations, is challenging. No universal governing equations, similar to the Navier–Stokes equations for Newtonian fluid flow, exist for these granular flows. It is therefore unsurprising that several canonical model flows, including heap flow, rotating tumbler flow and chute flow are frequently studied to better grasp the physics of granular flows. Here, we study bounded heap flow in which, unlike unbounded heap flow, an endwall forces the free surface height to increase with time. This flow geometry is shown in figure 1*a*. Typically, to study the flow a quasi-two-dimensional geometry is used in which the spanwise extent is small (O(10) particle diameters). In addition, most experiments consider a half heap where grains are fed from one side of the apparatus as shown in figure 1*a*.

Bounded heap flow differs in several key aspects from other canonical quasi-two-dimensional free surface granular flows, including rotating tumblers, unbounded heaps and inclined chutes [1] (figure 1*b*–*d*). For sufficiently large volumetric feed rate, *Q*, the surface flow is continuous (non-avalanching) and, because flowing grains are stopped by the endwall at the downstream end of the heap, the free surface rises steadily and uniformly along the length of the heap at a rise velocity, *v*_{r}. The local flow rate decreases linearly along the streamwise direction, resulting in a streamwise velocity gradient in the streamwise direction. By contrast, in unbounded heap flow and inclined chute flow, there is no endwall to stop the flow, so the local kinematics are invariant in the streamwise direction (i.e. d*u*/d*x*=0, where *u* is the streamwise velocity and *x* is the streamwise direction; figure 1*c*–*d*). While flow in rotating tumblers is similar to bounded heap flow in that *u* varies in the streamwise direction [2–5], the location of the free surface in the rotating tumbler remains fixed, as in unbounded heap and inclined chute flows (figure 1*b*–*d*). Thus, bounded heap flow is unique in comparison with these other flows in that it is never fully-developed and its free surface continually rises.

There have been extensive studies of kinematics in the canonical free surface flows mentioned earlier, for both monodisperse systems (e.g. [1–11] and references therein) and bidisperse systems [12–16], but only a few studies [17–19] have investigated the kinematics of bounded heap flow. By assuming a constant depth-averaged streamwise velocity along the flowing layer at each feed rate, Boutreux *et al.* [17,18] and Khakhar *et al.* [19] concluded that the local thickness of the flowing layer in bounded heap flows decreases along the streamwise direction. However, details of the kinematics such as profiles of streamwise and normal velocity components and solids volume fraction at different streamwise locations were not measured.

In addition to monodisperse granular systems, the kinematics of bidisperse systems during heap flow is of fundamental importance in driving segregation. Polydisperse granular materials tend to segregate during heap formation, often resulting in inhomogeneous final particle distributions. Heap segregation occurs in many contexts, particularly in industrial applications, and in most cases, the segregation is unwanted. Therefore, understanding the mechanisms of and developing a predictive model for segregation in heap flows is desirable. Previous studies on segregation of bidisperse mixtures of different-sized particles in quasi-two-dimensional bounded heap flow identified three different final particle configurations: stratified, in which there are layers of large and small particles [20–25]; segregated, in which small particles form the central portion of the heap and large particles form the outer portion [26–32] and mixed, in which the particles do not segregate [20,33]. Our recent experiments [34] showed that the transition between stratified and segregated states was controlled by the two-dimensional feed rate (*q*_{0}=*Q*/*T*, where *T* is the gap thickness between the two bounding side walls of the quasi-two-dimensional heap), while the transition between segregated and mixed states was determined by the heap rise velocity (*v*_{r}=*q*/*W*, where *W* is the horizontal width of the heap from the feed location to the outer bounding wall).

These different final particle configurations and the transitions between them are closely associated with the unique characteristics of the kinematics in bounded heap flow. For example, in a bounded heap, small particles percolate to the bottom of the flowing layer and deposit into the static bed, whereas large particles accumulate in the downstream region. The degree of segregation is influenced by the heap rise velocity [34], presumably due to the velocity of small particle percolation through the depth of the flowing layer compared with the rise velocity. In contrast, in unbounded heaps or inclined chutes, small particles segregate to the bottom of the flowing layer but continue to flow until they reach the end of the flowing layer. In rotating tumblers, particles exit the flowing layer and remain in the same streamlines before re-entering the flowing layer. As a result, in these flows, only the segregation rate is influenced by the flow rate but not the degree of segregation of the final state [35,36] as is the case for bounded heap flow.

Local particle distributions for segregating bidisperse mixtures in bounded heap flow have been modelled by Shinohara *et al.* [28] and Boutreux & de Gennes [17] by incorporating the percolation/segregation model proposed by Williams [26,27] and Bridgwater and co-workers [29,37,38] into a continuum framework. Although these two models predict some general features of heap segregation/stratification [21,31,32], the lack of details concerning the kinematics in the bounded heap flow necessitates unverified assumptions (e.g. the collision model in the study of Boutreux & de Gennes [17] or the velocity ratio of different sub-layers in the studies of Shinohara *et al*. [28] and Shinohara & Enstad [39]) with fitting parameters determined from experiments or simulations. Characterizing the kinematics of segregating bidisperse particles is necessary to enable better segregation models for bounded heap flow.

In this research, we perform a computational study of quasi-two-dimensional bounded heap flow using the discrete element method (DEM) to investigate the kinematics of both monodisperse and bidisperse systems. We examine the continuous flow regime in which the segregated state occurs [34]. The kinematics of intermittent avalanches, which can result in stratified states, is not examined here, although it has been explored in other studies [7,40,41]. In §2, the simulation technique and geometry are introduced, and comparisons with experiments are made to validate the simulations in terms of local particle distributions at steady-state and velocity profiles in the flowing layers. In §3, kinematics for both monodisperse and bidisperse systems are presented, including velocity profiles in both the streamwise and normal directions, and the local thickness of the flowing layer at particle size ratios, *R*, from 1 to 3 (where *R* is the ratio of large particle diameter, *d*_{l}, to small particle diameter, *d*_{s}) and a range of feed rates, *Q*. We find that at fixed *R* and *Q* the local surface streamwise velocity, *u*_{s}, and local mean shear rate, , decrease approximately linearly along the streamwise direction, while the thickness of the flowing layer, *δ*, remains roughly constant. We also find a scaling law that collapses all streamwise velocity profiles at different feed rates and particle size distributions onto a single curve. Section 4 presents our conclusions.

## 2. Simulation approach and experimental validation

### (a) Simulation method and geometry

The DEM is used here to simulate quasi-two-dimensional bounded granular heap flow. In DEM simulations, the translational and rotational motion of each particle are calculated by integrating Newton's Second Law. The forces between particles are repulsive and are non-zero only when particles are in contact. We used a linear-spring dashpot force model [42–45] to calculate the normal force between two contacting particles. This consists of two parts: a normal elastic force and a normal viscous damping force such that . Here, *ϵ* and *V*_{ij}=*v*_{i}−*v*_{j} denote the overlap and relative velocity of two contacting particles *i* and *j*, respectively. represents the unit vector in the normal direction between particles *i* and *j*, and *m*_{eff}=*m*_{i}*m*_{j}/(*m*_{i}+*m*_{j}) is the reduced mass. *k*_{n} and *γ*_{n} characterize the stiffness and damping of the granular materials, respectively, and are related to the collision time Δ*t* and restitution coefficient *e* by *γ*_{n}=−ln*e*/Δ*t* and [43,44]. For the tangential force, a linear spring at the contact point between two particles provides a restoring force. If this restoring force is larger than the Coulomb friction force, the spring is ‘cut’ and the force is sliding friction based on Coulomb's law. The tangential force can therefore be expressed as . Here, the tangential displacement *β* is given by [46], where *t*_{s} is the initial contact time between two particles. is the relative tangential velocity of two particles and is the unit vector in the tangential direction. The tangential stiffness is [43]. The velocity-Verlet algorithm [44] is used to update the positions and velocities of the particles.

The quasi-two-dimensional bounded heap simulated here is sketched in figure 2. To save computational cost, we simulate only the steady filling stage, where the heap contacts the bounding endwall and rises steadily and uniformly, which is similar to the experimental set-up of Drahun & Bridgwater [29]. To accomplish this, the bottom wall of the domain is inclined at an angle *θ*=24^{°} to horizontal, which is slightly less than the dynamic angle of repose *α* in our previous experiments [34]. Upon filling, particles in contact with the inclined bottom wall are immobilized to increase the effective wall friction similar to the physical situation where particles deposit on the heap and stop flowing. When the heap is sufficiently deep, approximately 10 particle diameter which corresponds to time *t*_{0} (≈10*d*_{l}/*v*_{r}), the effect of the bottom wall on the flowing layer is negligible, and the flow reaches a steady state comparable to our experiments. The width (*W*=0.46 m) and gap thickness (*T*=0.013 m) of the simulated domain are the same as in our previous experiments [34]. Particles are fed into the left end of the domain 0.1 m above the bottom wall at a volumetric feed rate *Q*. The width of the feed zone, the initial velocity and the packing density of the particles in the feed stream are varied to achieve different *Q*. We performed several trial computational runs and found the kinematics and segregation insensitive to these feed parameters at constant *Q*.

The material density of the simulated particles, *ρ*, is 2500 kg m^{−3} and restitution coefficient, *e*, is 0.8. Particle–particle and particle–wall friction coefficients are set to *μ*=0.4. Note that the side walls are frictional and flat to allow direct comparisons to experiments. To decrease computation time, the binary collision time is set to Δ*t*=10^{−3} s, consistent with previous simulations [47] and sufficient for modelling hard spheres [48] based on comparison with results for Δ*t*=10^{−4} s. The integration time step is Δ*t*/100=1.0×10^{−5} s to assure numerical stability. The particles are millimetre-sized, with each particle species *i* having a mean particle diameter between 1 and 3 mm. To prevent crystallization, the particles have a uniform size distribution with a variance of (0.1*d*_{i})^{2}. There are up to one million particles in our simulations depending on *Q* and particle size. Simulations typically run for 15–100 s of physical time, and the time to reach steady-state *t*_{0} is approximately one-third of the total simulation time.

As shown in figure 2, two different coordinate systems are used. The first, the fixed laboratory coordinate system (*x*_{lab},*y*_{lab},*z*_{lab}) with origin at the left bottom corner of the domain, is used to measure particle distributions of the final states. The second coordinate system (*x*,*y*,*z*) is rotated by the dynamic angle of repose *α* of the flowing particles and moves upward so that the *x*-axis is always at the free surface. It is used for the time-averaged kinematics in the flowing layer. The *x*-axis is along the streamwise direction; the *y*-axis is perpendicular to the side walls and the *z*-axis is normal to the free surface. The origin in this moving coordinate system is located at the intersection between the front wall, the end of the feed zone and the free surface. The (*x*,*y*,*z*) velocities in the flowing layer are (*u*,*v*,*w*), respectively.

### (b) Validation of discrete element method simulation

The DEM simulations in this study are validated by comparing the depth dependence of the streamwise velocity and the spatial variation of species concentration to the same quantities obtained from experiments. These variables are average values calculated over the duration of the steady state of the simulation (*t*>*t*_{0}) which corresponds to a heap rise of approximately 50 particle diameters. Agreement between streamwise velocity profiles obtained in experiments using Particle Tracking Velocimetry [5] and simulations is excellent and will be discussed in §3*c* after we present results from simulations.

Figure 3 shows a comparison of final states between DEM simulations and our previous experiments [34] for identical conditions with a bidisperse mixture of different-sized particles. We plot the profiles of volume concentration of small particles *c*_{s}=*f*_{s}/(*f*_{s}+*f*_{l}) in the steady filling stage along the *x*-direction excluding the flowing layer, where *f*_{s} and *f*_{l} are the solids volume fraction of small and large particles, respectively. The volume concentration profile for the experiment is obtained using image processing techniques, where the local particle concentration of each species is calculated by first measuring the intensity of each pixel and then calibrating by the intensity of mono-sized particles for each species. Note that the concentration profile measured from the experiment represents the particle distributions at the front wall. To measure the concentration profiles in the simulation, the domain is divided into Δ*x*=10 mm wide slices in the *x*-direction and averaged in the *y*-direction over particles close to the side wall (*y*<2*d*_{l}) to calculate the volume concentration of each species in each bin. Figure 3 shows that excellent agreement is obtained between the simulations and experiments for a mixture of 1 and 2 mm particles at *Q*=22 and 80 cm^{3} s^{−1}. We also find excellent agreement between simulations and experiments for other flow rates and particle size distributions. Finally, concentration and kinematics vary slightly in the *y*-direction due to wall effects (see appendix A), which is similar to what occurs in unbounded heap flow [49,50], though the scaling described later in this paper is unaffected.

## 3. Results

### (a) Free surface, dynamic repose angle and rise velocity

Unlike other free surface granular flows, in bounded heap flow the free surface rises as particles accumulate on the heap. Therefore, to measure time-averaged kinematics such as the velocity profile and thickness of the flowing layer at different locations, the free surface location needs to be determined at each instant of time. To do so, the computational domain is divided into equal, non-overlapping bins of Δ*x*×Δ*y*×Δ*z*, where Δ*x*=10 mm, Δ*y*=*T* and Δ*z*=0.5 mm. The local solids volume fraction in each bin was calculated as *f*_{i}=*V* _{i}/*V* _{bin}, where *V* _{i} is the fractional volume of all particles located in bin *i*, and *V* _{bin}=Δ*xΔyΔz* is the bin volume. Figure 4*a* shows profiles of solids volume fraction in the *z*-direction at different streamwise locations for a mixture of equal volumes of 1 and 2 mm particles for steady flow (*t*_{0}=10 s) at *Q*=22 cm^{3} s^{−1}. At each streamwise location, the solids volume fraction is *f*≈ 0.65 in the static portion of the heap for small *z*_{lab} (for example, in the region *z*_{lab}<0.17 m for *x*_{lab}=0.1 m). Moving upward in the heap, the solids volume fraction decreases slightly in the creeping region of the heap due to slow re-arrangement of particles (0.17 m <*z*_{lab}<0.185 m for *x*_{lab}=0.1 m). Moving upward through the flowing layer (0.185 m <*z*_{lab}<0.2 m for *x*_{lab}=0.1 m) the solids volume fraction continues to decrease slowly until it decreases rapidly to zero close to the free surface. Above the free surface in the upstream region of the heap (e.g. *z*_{lab}>0.205 m for *x*_{lab}=0.1 m), the solids volume fraction is non-zero due to the bouncing of some particles after impact on the heap.

Based on the profiles of solids volume fraction, the location of the free surface, *z*_{s}, at each streamwise position can be estimated based on a cutoff value of solids fraction *f*_{c} (similar to [51]). *f*_{c} is typically selected at an intermediate value, since too small a value includes bouncing particles in the flowing layer and too large a value includes only static particles. In this paper, we use *f*_{c}=0.35 because the measured free surface only varies by one particle diameter when *f*_{c} is changed between 0.1 and 0.5. The location of the entire free surface at different times is plotted in figure 4*b*, indicating a uniformly sloped surface. The dynamic angle of repose of the heap, *α*, at different times was determined by calculating the slope of the free surface using a linear fit. As shown in the inset of figure 4*b*, *α* remains nearly constant (25.69^{°}±0.07^{°} for the mixture of 1 and 2 mm particles at *Q*=22 cm^{3} s^{−1}) during the entire course of the steady-state portion of the simulation.

When flow is continuous (as opposed to intermittently avalanching), the heap rises steadily and uniformly. Figure 4*c* shows the free surface position as function of time at three different streamwise locations. The linear increase of *z*_{s} with increasing time confirms that the free surface rises steadily, as the slope of *z*_{s} versus *t* is the rise velocity. To facilitate locating the free surface at different times and streamwise locations, we determine the free surface *z*_{s}(*x*_{lab},*t*_{0}) at the initial time *t*_{0} for the steady state using the linear fit as in figure 4*b*. Then subsequent free surface locations are calculated as . Absolute scaled differences between the calculated free surface and the measured free surface from simulations, , is always less than 1, indicating that this approach is accurate.

With the free surface location and the dynamic angle of repose known as functions of time, we can transform from the laboratory (*x*_{lab},*y*_{lab},*z*_{lab}) to the instantaneous moving reference frame at the free surface (*x*,*y*,*z*). From here on, the kinematics at each time are measured in the (*x*,*y*,*z*) coordinates, with corresponding velocity field (*u*,*v*,*w*). The *z*-component of the rise velocity in this coordinate system is *v*_{r}′=*v*_{r}cos*α*.

### (b) Kinematics of monodisperse systems

#### (i) Streamwise velocity

We measure the local, time-averaged, streamwise velocity *u* during steady state for each simulation run (and the normal velocity *w*) by dividing the domain into equal, non-overlapping bins of size Δ*x*×Δ*y*×Δ*z*, where Δ*x*=10 mm, Δ*y*=*T*, and Δ*z*=1 mm. The streamwise velocity in bin *i* averaged over *δt* is given by , where *V* _{ij} is the fractional volume of particle *j* located in bin *i*, *u*_{j} is the velocity component in the streamwise direction of particle *j*, and *N*_{i} is the total number of particles that are partially or entirely located in bin *i* during *δt*. Figure 5*a* shows profiles of *u* at steady state (*t*>*t*_{0}) in the *z*-direction at four different streamwise locations averaged over *δt*=10 s for 1.5 mm monodisperse particles with *Q*=22 cm^{3} s^{−1}. The velocity profiles exhibit two regimes: a rapid decrease of *u* from a maximum value at the free surface in most of the flowing layer and a slow decay to the creeping region near the bottom of the flowing layer. In the quasi-static region below the flowing layer, particles move at small non-zero velocities O(1 mm s^{−1}). The ‘error’ bars in the plots represent the standard deviations of time-averaged *u* (associated with granular temperature), and show larger standard deviations for larger *u*, which is consistent with previous results (e.g. Jain *et al.* [5] in rotating tumblers). The streamwise velocity profiles resemble those in unbounded heap flow [1,8,49,50,52] or rotating tumbler flow [4,5,53]. However, in bounded heap flow, the streamwise velocity at constant depth decreases in the streamwise direction.

Figure 5*b* shows profiles of shear rate in the depth direction calculated from the velocity profiles in figure 5*a*. The shear rate at constant depth decreases in the streamwise direction. In a region that is only a few particles thick near the free surface, the shear rate increases slightly to a maximum value as *z* decreases. Below this region, the shear rate decreases smoothly to zero as the static portion of the heap is approached. The shear rate profile in the depth direction is neither constant nor exponential.

#### (ii) Kinematics along the streamwise direction

Since the local flow rate in bounded heap flow decreases linearly along the streamwise direction, kinematic properties may also change along the streamwise direction. In figure 6*a*, the local, time-averaged, two-dimensional flow rate is plotted as a function of *x* and shows a linear decrease along the streamwise direction as expected due to the uniform deposition of material on the heap. The flow rate reaches zero at the downstream endwall at *x*=0.42 m. Near *x*=0, the small deviation from a purely linear decrease is due to a loss of flux as a result of excluding the bouncing particles in the calculation of *q*. Figure 6*b* shows that the velocity at the free surface, *u*_{s}, also decreases approximately linearly in the streamwise direction. The depth-averaged shear rate, , is calculated by averaging the local shear rate ∂*u*/∂*z* over the flowing layer at each *x*. decreases from a maximum near *x*=0 to nearly zero at the downstream endwall as shown in figure 6*c*.

To determine the thickness of the flowing layer, it is necessary to locate the bottom of the flowing layer *z*_{bottom} (or equivalently, the boundary between the flowing layer and the quasi-static region). To do this, we tested three different methods based on the streamwise velocity profiles. In the first method, similar to the earlier studies [8,54,55], the velocity profile at each streamwise location is fit to *u*(*x*,*z*)=*u*_{o}(*x*)exp(*z*/*z*_{o}), where *u*_{o} is the nominal surface velocity at each *x* and *z*_{o} is a characteristic depth, to which *z*_{bottom} is proportional. In the second method, the bottom of the flowing layer is determined by extrapolating the approximately linear part of the velocity profile to zero [1]. In the third method, a cutoff value *u*_{c}, proportional to *u*_{s} at each *x* (e.g. 0.1*u*_{s}) determines *z*_{bottom}. Figure 6*d* shows *δ* as a function of *x* determined by these three methods. The values calculated for *z*_{bottom} and, consequently, *δ* are nearly identical, provided that the appropriate scale factor (*δ*=2.3*z*_{o}) for the first method and cutoff value (0.1*u*_{s}) for the third method are used. Similar to other free surface flows, *δ* spans a few particle diameters (6*d* to 7*d* for most of the length of the flowing layer shown in figure 6*d*) [3,56,57]. The flowing layer thickness decreases slightly along most of the length of the flowing layer, but decreases more at the end of the heap.

Khakhar *et al.* [19] proposed a relation between *δ* and in bounded heap flow along the streamwise direction
3.1where *L* is the length of the flowing layer and *δ*_{L} is the thickness of the flowing layer at *x*=*L*. Using *v*_{r}, and *δ*_{L} from simulation, we found that *δ* calculated from equation (3.1) matches well with the values measured directly from the simulation, as shown in figure 6*d*. However, our results show a decrease of along the streamwise direction, contrary to Khakhar *et al.* [19] where was predicted to be constant in the streamwise direction. This discrepancy is likely due to the measurement technique used by Khakhar *et al.* [19] to determine *δ* in which they manually located the bottom of the flowing layer using a constant streamwise velocity cutoff instead of the relative value used here.

#### (iii) Streamwise and normal velocity scaling

The streamwise velocity profiles at all streamwise locations collapse onto a single curve when *u* is normalized by *u*_{s} and *z* by *δ*, as shown in figure 7*a*. Throughout most of the flowing layer (|*z*/*δ*|<1), the scaled velocity *u*/*u*_{s} decreases rapidly from 1 at the free surface. Deeper in the flowing layer, *u*/*u*_{s} decreases approximately exponentially, similar to unbounded heap flow or rotating tumbler flow [1,8,55,58]. This exponential tail can be seen more clearly in figure 7*b*, where *u*/*u*_{s} is plotted on a log scale. In deeper regions (|*z*/*δ*|<2), velocities are small (|*u*/*u*_{s}|<0.01) but non-zero, though there is substantial scatter. Note that in unbounded heap flow or rotating tumbler flow, the velocity can be scaled by at different flow rates [1], where *d* is the particle diameter for monodisperse systems. However, since *u*_{s} changes with *x* in bounded heap flow (figure 5*a*), this constant scaling for *u* is insufficient.

A difference between bounded heap flow and other free surface flows is that particles in the flowing layer have a positive average velocity normal to the free surface due to the rise of the heap. For bidisperse flows, this positive normal velocity contributes to the separation of small and large particles in the streamwise direction which in turn controls the final segregation state of the heap. This does not occur in other free surface flows due to the unchanging location of the free surface (figure 1). Because of the important role of the vertical movement of the flowing layer in the final segregation configuration in the bounded heap, an analytical expression for the normal velocity of the heap is useful. To this end, we start with the conservation of mass
3.2and assume a linear decrease of the streamwise velocity *u* in the *x*-direction based on the results in §3*b*(i). Then, to first approximation,
3.3where *u*(0,*z*) gives the depth dependence at *x*=0. Substituting equation (3.3) into equation (3.2), integrating with the boundary condition *w*=0 at *z*=0 in the moving reference frame, and noting that *w* is a function of *z* only (uniform rise of the heap), an expression for the normal velocity *w*(*z*) in the flowing layer is obtained
3.4

The normal velocity profile can be calculated based on equation (3.4), if *u*(0,*z*) is known. Even though the often-used linear or exponential relations do not capture the exact functional form of the streamwise velocity as shown in figure 5, we try both forms for *u*(0,*z*), a linear expression *u*(0,*z*)=*u*(0,0)(1+*z*/*δ*) and an exponential expression *u*(0,*z*)=*u*(0,0)e^{kz/δ}, to obtain analytical solutions for the normal velocity. Here, *u*(0,0) is the surface velocity at *x*=0, and *k* is the ratio of the thickness of the flowing layer to the characteristic length scale of the exponential fit, where *k*=2.3 in this study. *u*(0,0) is associated with the two-dimensional flow rate *q* at *x*=0 by . For the linear streamwise velocity profile, *u*(0,0)=2*q*_{0}/*δ*, and for the exponential streamwise velocity profile, *u*(0,0)=*kq*_{0}/*δ*(1−e^{−k}). Substituting *u*(0,0), *q*_{0}=*v*_{r}*L*, and *u*(0,*z*) into equation (3.4), the normal velocity for the linear streamwise velocity profile is
3.5

The exponential streamwise velocity profile yields
3.6where and . Note that both equations (3.5) and (3.6) automatically satisfy the boundary condition *w*=−*v*_{r}′ at *z*=−*δ*.

Figure 8 shows theoretical predictions from equations (3.5) and (3.6) (curves), and time-averaged normal velocity measured from simulations in the flowing layer for 1.5 mm particles at *Q*=22 cm^{3} s^{−1} at different streamwise locations. Both simulation data and theoretical curves show that normal velocity in the moving reference frame decreases from zero at the free surface to −*v*_{r}′ at the bottom of the flowing layer, though the data are somewhat scattered at different streamwise locations due to the stochastic nature of the flow. The theoretical predictions from both equations (3.5) and (3.6) agree with the simulation data, even though some deviations exist. For example, in the upper region of the flowing layer (|*z*/*δ*|<0.4), both analytical solutions slightly underpredict normal velocity. Nonetheless, both analytical solutions provide reasonable predictions for the normal velocity.

### (c) Kinematics at different feed rates and particle size distributions

The preceding section showed that a universal streamwise velocity profile exists for a monodisperse system at *Q*=22 cm^{3} s^{−1}. However, it is not clear if this scaling remains applicable when the feed rate and particle size distributions change. The particle size distribution is characterized by size ratio *R* (*R*=1 for monodisperse systems and *R*>1 for bidisperse systems) and local particle mean diameter: , where *n*_{s}=*c*_{s}*R*^{3}/(*c*_{s}*R*^{3}+*c*_{l}) and *n*_{l}=*c*_{l}/(*c*_{s}*R*^{3}+*c*_{l}) are local number fractions, and *c*_{s} and *c*_{l} are local volume concentrations of small and large particles, respectively.

To test whether the streamwise velocity profile could be collapsed onto a single curve for other feed rates and bidisperse mixtures, a total of 11 simulations were performed at four different feed rates including four monodisperse systems and seven bidisperse systems with equal volumes of each species as listed in table 1. The streamwise velocity profiles at different streamwise locations for all simulations are plotted in figure 9. Once again, all simulation data collapse onto a single curve identical to that in figure 7 and independent of feed rate and particle size distribution, when *u* is normalized by *u*_{s} and *z* is normalized by *δ*. Clearly, the scaling is valid over a broad range of flow rates and particle size distributions. Furthermore, streamwise velocity profiles from experiments obtained using Particle Tracking Velocimetry^{1} show excellent quantitative agreement with the scaled simulation results (open circles in figure 9 are experimental results for 3 mm monodispere glass particles at *x*/*L*=0.5 and *Q*=14 cm^{3} s^{−1}). In addition, we have considered the streamwise velocity profiles for different gap thicknesses and different locations across the gap and found that they too collapse onto the single curve in figure 9 using the same scaling relation (not shown).

We further investigate the feed rate and particle size distribution dependence of surface velocity, shear rate, and flowing layer thickness along the streamwise direction in figure 10. The surface velocity non-dimensionalized by the surface velocity at *x*=0, *u*(0,0)=2*q*_{0}/*δ* from §3*b*(iii), is *u*_{s}*δ*/(2*q*_{0}), which can be plotted along the streamwise direction at all feed rates and particle size distributions. As shown in figure 10*a*, plotting *u*_{s}*δ*/(2*q*_{0}) along the streamwise direction collapses the data onto a single curve that decreases linearly along the *x*-direction to the end of the flowing layer, though deviations occur close to the feed zone (*x*=0) which are probably due to particles bouncing as they fall onto the heap. The dimensionless shear rate is plotted in figure 10*b*. Again, the data for the different simulation runs collapse onto a single curve that decreases linearly with *x*. Since the scaling of *u*_{s} and depend on the local thickness of the flowing layer *δ*, its variation with flow rate and particle size distribution needs further investigation. Figure 10*c* shows *δ* as a function of *x*/*L* at different flow rates and particle size distributions. Except for close to the feed zone (*x*<0.2) where particle bouncing effects occur, the flowing layer thickness generally increases as the feed rate increases for the same particle size distribution (namely, the same size ratio and particle diameter of each component). However, the profiles of the scaled thickness of the flowing layer in the streamwise direction (figure 10*d*) show no clear trends when feed rates and particle size distributions are changed.

The dependence of *δ* on flow rate and particle size has been studied in other free surface flows such as rotating tumbler flow or unbounded heap flow for monodisperse systems (e.g. [1,57,59] and references therein). In all cases, , where is the dimensionless flow rate. In figure 11, we plot as a function of local *q** for *x*/*L*≥0.2 on a log–log scale and find that data from all simulations appear to follow a power-law relationship between and *q** with *a*≈0.22. This value for *a* differs from that reported by GDR Midi [1] and Renouf *et al.* [59], where *a*≈0.5, but is close to what Pignatel *et al.* [57] found at comparable *q** in rotating tumblers, where *a* ranged from 0.1 to 0.45. Perhaps, what is most important here is that figure 11 can be used to predict the flowing layer thickness, at least approximately, over a wide range of flow rates and particle sizes. Note that this scaling relation is also valid for different gap thicknesses and different locations across the gap.

## 4. Discussion and conclusions

We have shown that streamwise velocity profiles in continuous, bounded heap flow at different feed rates and particle size distributions collapse onto a single curve when depth is scaled by the local flowing layer thickness *δ* and *u* is scaled by the local surface velocity *u*_{s} (figure 9), where *u*_{s} depends on the feed rate *q*_{0} and *δ* (figure 10*a*). This result demonstrates a universal functional form for the streamwise velocity in bounded heap flow for both monodisperse and bidisperse systems with the local flowing layer thickness being the only scaling parameter at a given feed rate. The fact that streamwise velocity profiles are independent of the local particle size distributions associated with segregation of bidisperse particles along the streamwise direction is different from other free surface flows such as inclined chute flow [13,16] or rotating tumbler flow [14]. Specifically, Rognon *et al.* [13] found that in bidisperse systems the segregating larger particles in the upper portion of the flowing layer have a relative smaller velocity compared with the monodisperse systems as they slide over the smaller particles in the lower portion of the flowing layer, which results in a different velocity profile from a monodisperse system. Tripathi & Khakhar [16] further showed that varying the local concentration of large particles can result in a non-monotonic change of streamwise velocity, where a minimum streamwise velocity occurs at a mass fraction of 70 per cent large particles. They also showed that velocity profiles depend on the size ratio of the two species. Although not explicitly stated, Hill & Zhang [14] showed that in a rotating tumbler with a bidisperse mixture, the streamwise velocities at the upper portion of the flowing layer decrease as more large particles segregate to this region. The underlying mechanism for the different effects of particle size distribution on the streamwise velocity between bounded heap flow and other free surface flows is not clear, but might be associated with the linear decrease in local flow rate due to uniform particle deposition into the static bed. For certain parameter values, this effect seems to dominate over the different mobilities of large and small particles in the flow when they segregate.

As we have discussed, the flowing layer thickness can be determined from the velocity profile using any of several criteria as shown in §3*b*(ii). However, our results demonstrate that the flowing layer thickness as a scaling parameter is insensitive to the specific functional form of the streamwise velocity profile and the specific measurement approach. The local flowing layer thickness depends only on the local flow rate and the local average particle diameter (figure 11). For a monodisperse system, the velocity field is well defined by these two variables. For a bidisperse system, however, local mean particle diameter changes as segregation occurs, so determining the streamwise velocity for bidisperse system requires a prediction of the local concentration of each species.

Although we have found a scaling for the streamwise velocity in bounded quasi-two-dimensional heap flow, unresolved issues remain. First, while experimental and simulation results provide a suitable empirical velocity profile, to theoretically determine the functional form of the streamwise velocity profile, a constitutive law for the stress is necessary. With recent progress on the rheology of dense granular flow, a local rheological model similar to that proposed by Jop *et al.* [60] or a non-local model similar to that proposed by Kamrin & Koval [61] might be applied to bounded heap flow. Second, quasi-two-dimensional bounded heap flow can be thought of as a section of three-dimensional bounded heap flow [34]. The fact that the scaling we found does not depend on gap thickness or the location across the gap in quasi-two-dimensional flow (§3*c*) suggests that this scaling could be directly extended to three-dimensional flow where the sidewall does not exist.

The mean flow kinematics of the bidisperse granular systems presented here provide a necessary precursor for modelling segregation in bounded heap flow. The motion of each species in the segregating mixture can be determined by superposing the relative motion of the species—the segregation velocity—onto the mean flow. Thus, the local concentration of each species can be determined using the transport equation for each species and considering the combined effects of advective fluxes related to the mean flow, the segregation fluxes and the diffusive fluxes due to particle collisions. This approach has been used to model size segregation in other free surface flows, such as inclined chute flow [62–65]. However, in these studies, the kinematics are relatively simple: the flow is fully developed (∂*u*/∂*x*=0), mean flow has no normal component (*w*=0), and the streamwise velocity profile is linear. In contrast, bounded heap flow is more complicated in that the mean flow decelerates in the streamwise direction (∂*u*/∂*x*<0) and the normal velocity profile is non-zero and nonlinear. As a result, the advection of each species due to mean flow may become important, and the shear rate could influence both the segregation velocity of each species [66–68] and the diffusion coefficients [69]. These effects probably need to be included when modelling segregation in bounded heap flow.

## Funding statement

We gratefully acknowledge financial support from The Dow Chemical Company.

## Acknowledgements

We are grateful for helpful discussions with Karl Jacob and Ben Freireich.

## Appendix A. Wall effects

In quasi-two-dimensional unbounded heap flow or inclined chute flow, side walls influence flow kinematics and rheology [49,50]. Particles are slowed by wall friction resulting in a blunt velocity profile with only a slight decrease in velocity near the side walls. In bounded heap flow, the side walls have a similar effect on the streamwise velocity. Figure 12*a* shows the streamwise surface velocity *u*_{s} averaged over different portions of the gap between side walls. Although surface velocities in different portions of the gap exhibit a similar linear decrease in the streamwise direction, surface velocities close to the wall (data averaged over 4 mm thick slices in the *y*-direction adjacent to both side walls) are roughly 10 per cent smaller than mean surface velocities, while those in the central region of the heap (4<*y*<8.7 mm) are 10 per cent larger. The strength of segregation is also affected by the side walls. As shown in figure 12*b*, the profiles of volume concentration of small particles *c*_{s} in the streamwise direction show enhanced segregation near the side walls. Although *c*_{s} decreases from the side walls to the centre in the upstream region of the heap (*x*_{lab}/*W*<0.6), *c*_{s} does not vary in the spanwise direction in the downstream region. This indicates that a small degree of horizontal segregation occurs in the lateral direction, presumably because small particles segregate towards the side walls due to wall exclusion of large particles or shear-induced segregation [70]. However, as mentioned in §3*c*, the scalings for the streamwise velocity and the flowing layer thickness are valid at different gap thicknesses, *T*, and different locations across the gap. Therefore, in this work, we report kinematic quantities averaged over the entire gap (0<*y*<*T*).

## Footnotes

↵1 Videos of particle flow at the wall in a 36

*d*wide by 26*d*high region were acquired at 300 fps using a Casio EX-F1 camera. Before processing, images were shifted vertically so that the free surface was at a fixed location. Then velocity fields were calculated using the Particle Tracking Velocimetry technique [5].

- Received April 15, 2013.
- Accepted June 18, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.