## Abstract

Preventing segregation in flowing granular mixtures is an ongoing challenge for industrial processes that involve the handling of bulk solids. A recent continuum-based modelling approach accurately predicts spatial concentration fields in a variety of flow geometries for mixtures varying in particle size. This approach captures the interplay between advection, diffusion and segregation using kinematic information obtained from experiments and/or discrete element method (DEM) simulations combined with an empirically determined relation for the segregation velocity. Here, we extend the model to include density-driven segregation, thereby validating the approach for the two important cases of practical interest. DEM simulations of density bidisperse flows of mono-sized particles in a quasi-two-dimensional-bounded heap were performed to determine the dependence of the density-driven segregation velocity on local shear rate and particle concentration. The model yields theoretical predictions of segregation patterns that quantitatively match the DEM simulations over a range of density ratios and flow rates. Matching experiments reproduce the segregation patterns and quantitative segregation profiles obtained in both the simulations and the model, thereby demonstrating that the modelling approach captures the essential physics of density-driven segregation in granular heap flow.

## 1. Introduction

Granular materials with different particle properties tend to segregate spontaneously when they are flowing [1–4] or vibrating [5–7]. Such segregation is frequently encountered in industrial processes that involve handling bulk solids [8,9] as well as in geophysical transport such as debris flows [10], pyroclastic flows [11] and mineral transport [12]. Thus, modelling and predicting segregation is important, yet accurate models that can be broadly applied are only just now being developed.

Among different particle properties that can drive segregation, particle size [13] and density [14] are critical factors. The focus of this work is the segregation due to differences in particle density, which can occur in vibrated granular mixtures [15–19], free surface flows [20–24] and vertical chute flows [25]. In gravity-driven free surface flows, particles with lower density are more likely to rise to the free surface while particles with higher density are more likely to segregate to the bottom of the flowing layer, resulting in segregation patterns such as a segregated core or streaks of heavier particles in rotating tumblers [14,21–23]. While particle-based simulation methods can reproduce density-driven segregation phenomena on a small scale, an accurate continuum-based model would be of clear practical and theoretical value.

Various continuum models have been proposed for predicting segregation in granular flows. Bridgwater [26] developed a continuum-based model that uses a segregation velocity based on percolation due to particle size differences. Savage & Lun [13] applied a statistical mechanics approach based on kinetic sieving and squeeze expulsion mechanisms and derived a size segregation velocity related to the particle properties, the local concentration and the shear rate, which was assumed constant. Variations of these continuum models have been applied to other geometries such as chute [24,27–36] and annular shear [37] flows, and some have achieved qualitative agreement with simulations and experiments. The continuum models are also evolving with more details being considered. Gray & Chugunov [28] successfully included the effect of diffusion into the continuum modelling framework. Marks *et al.* [35] included the effect of local shear rate on the segregation flux and achieved qualitative agreement with simulations. Recently, Tripathi & Khakhar [20], Tunuguntla *et al.* [24] and Gray & Ancey [36] incorporated bidisperse density segregation into the continuum model. These studies applied segregation velocities proportional to the normalized density difference and other parameters. Although the models showed a degree of agreement with simulations, validation of the segregation velocity or segregation patterns was not considered in detail. Here, we propose a continuum model for density bidisperse segregation with a different approach to defining the segregation velocity, and we validate it by comparing predictions of this approach to both discrete element method (DEM) simulations and experiments.

In our recent work, we developed a continuum-based approach for predicting segregation of granular materials that achieves quantitative agreement with simulations and experiments of size bidisperse, multidisperse and polydisperse granular materials in different geometries [38–41]. Similar to the previous models discussed above, this model is also based on the transport equation
*L* and depth *δ*), where segregation occurs in most gravity-driven flows. We define *x* as the streamwise direction (0<*x*<*L*) and *z* as the normal direction (−*δ*<*z*<0), with *z*=0 at the surface of the flowing layer. Terms with subscripts *i* represent the properties for species *i* in a bidispere mixture (*h* for heavy and *l* for light in the case of density segregation), while terms without subscripts represent the average flow properties of both species. The concentration of species *i* is defined as *c*_{i}=*f*_{i}/*f*, where *f*_{i} is the volume fraction for species *i*, and *f* is the total volume fraction of both species. The mean two-dimensional (2D) velocity field is *D* is the diffusion coefficient. The segregation velocity, *w*_{s,i}, is defined as the relative normal velocity component of species *i* with respect to the total normal velocity component of both species: *w*_{s,i}=*w*_{i}−*w*. In the segregation term (∂/∂*z*)(*w*_{s,i}*c*_{i}), only flows normal to the free surface are considered, as segregation occurs primarily in this direction, and the gradient of concentration in the streamwise direction is small. Unlike previous approaches where unknown free constants are retained in the models [27,28], or closure relations are proposed by coupling the transport equation with constitutive theories [20,35], we inform our model with physical control parameters and kinematic parameters acquired from DEM simulations [40] or experiments in a similar way as explored in Thornton *et al.* [33] for segregation in a chute flow. Thus, no fitting parameters are needed. While we consider the quasi-2D case here, the model can be extended to fully three-dimensional systems.

Here, this continuum model is used to predict the segregation of density bidisperse granular materials in a quasi-2D-bounded heap, which typically exhibits complicated kinematics [42–45] with different segregation patterns including stratified layers of the two types of particles [46], fully segregated particles [47,48] and mixed particles with no segregation [49]. Here we study the continuous flow regime for which stratification does not occur [49]. DEM simulations are performed to determine the kinematics of density bidisperse flows, and experiments are performed to verify the results of simulations and theoretical predictions. In §2, the quasi-2D-bounded heap geometry, the DEM simulation methods and the experiments are described. In §3, we show that the simulations quantitatively reproduce the experimental results and discuss the flow kinematics. In §4, the continuum model (equation (1.1)) is non-dimensionalized and solved numerically. The results are compared with simulation and experimental results for different cases, and the influence of physical control parameters on segregation is discussed. Conclusions are presented in §5.

## 2. Simulation and experimental methods

### (a) Discrete element method simulations

In DEM simulations, the translational and rotational momenta of each particle are tracked using integration of Newton's Second Law. As in our previous work [38–40,45], the normal force model used in this research is the linear-spring dashpot model [50–52], in which the normal contact force between two particles is *ϵ* and *V*_{ij} represent the overlap and relative velocity between two contacting particles *i* and *j*, respectively. The unit normal vector between two particles is *m*_{eff}=*m*_{i}*m*_{j}/(*m*_{i}+*m*_{j}) denotes the effective mass. The normal stiffness *k*_{n} and damping *γ*_{n} are determined from the restitution coefficient *e* and binary collision time *t*_{c}: *γ*_{n}=−ln(*e*)/*t*_{c}, where ln is the natural logarithm. The tangential force model is given by the linear spring model with Coulomb friction [51], which can be expressed as *t*_{s} is the initial contact time and *μ*, and the unit vector in the tangential direction is *ρ* often have different surface and elastic properties. Here, to ensure that density is the only driving mechanism for segregation in simulations, identical material properties are applied except for the density. The binary collision time is *t*_{c}=10^{−3} s, which is small enough to accurately describe the flow of hard spheres [45]. The restitution coefficient is *e*=0.9 and the friction coefficient is *μ*=0.2. These values are selected so that the dynamic repose angle *α* in simulations matches the dynamic repose angle in experiments. Negligible differences in segregation occur over a range of values for *t*_{c}, *e* and *μ*, indicating that the simulations are relatively insensitive to the specific values used within the range of feed rates and density ratios tested in this specific geometry. For particle-wall contact, the same force models and material properties are used with the walls, which have infinite mass and radius. The time step for the simulations is set to *t*_{c}/40=2.5×10^{−5} s, which ensures numerical stability [40].

A schematic of the simulation geometry is shown in figure 1. The quasi-2D-bounded heap consists of two parallel plates with two bounding endwalls. The width between the bounding endwalls is *W*=0.5 m and gap thickness between the parallel plates is *T*=0.016 m. To save computation time, the bottom wall is inclined at an angle *α*_{0}=16°, which is smaller than the dynamic repose angle *α* (which ranges from 18° to 21° in different cases). Initially, the bottom wall is covered with a layer of immobilized particles. After the particles flowing into the system form a 10–15 particle diameters thick layer, the velocity profiles and concentration profiles in the flowing layer become steady, indicating that the effect of the bottom wall can be neglected. In simulations, density bidisperse particles enter the system at a volume feed rate of *Q* and volume ratio of 1:1. The particle diameter *d* is uniformly distributed with a variance of ±0.1*d* to reduce crystallization. Particles of 2, 3 and 4 mm diameters are simulated. The flow of mixed particles onto the heap is controlled by a rising gate. The rising gate eliminates bouncing particles caused by free fall of the particles [49], which can influence density segregation dramatically in small systems like this one. The gate, located at *W*_{F}=0.06 m, controls the vertical rise velocity *v*_{r}=*Q*/*WT* of the heap surface. For data analysis, we neglect flow in the feed zone and the area affected by the feed zone which extends to *x* is in the streamwise direction, *y* is in the thickness direction and *z* is normal to the surface of the flowing layer. *u*, *v* and *w* are the velocities in the *x*, *y* and *z* directions, respectively.

Simulations with different feed rates, density ratios and particle sizes were performed on an Nvidia GTX 780 graphic card (Graphics Processing Unit) with a parallelized DEM algorithm, as in our previous work [40]. Details of the kinematics are discussed in §3.

### (b) Experiments

To validate the DEM simulations and theoretical predictions, experiments were performed with equal diameter particles with different densities, as indicated in table 1. For each particle type, 100 sample particles were randomly selected and the diameter of each and the total weight were measured. Table 1 lists the average and standard deviation of the diameter and the density calculated by dividing the sum of the volume of each individual particle by the total weight measured for each particle type.

The geometry of the experimental system is the same as the simulation geometry (*W*=0.5 m and *T*=0.016 m). Particle mixtures were held in a hopper and fed into the system by an auger feeder (Acrison, Inc., NJ, USA) at the desired volume feed rate. The rising gate was implemented as a vertical metal bar lifted by a linear actuator (Firgelli Automations Inc., WA, USA) with a control board (Firgelli Technologies Inc, Canada). The experiments were recorded using a high-speed camera (Point Grey Research Inc., Canada) with frame rates up to 400 frames s^{−1}. Video images were obtained during steady filling of the heap at the downstream end of the flowing layer in contact with the vertical bounding wall and were analysed to provide concentration profiles of the segregation pattern in the fixed bed and velocity profiles in the flowing layer.

The average image intensity was used to calculate the concentration profile in the streamwise direction for particles in the fixed bed below the flowing layer. The region in figure 2*a* outlined by the white box was analysed to characterize the final segregation pattern achieved during steady filling. The boxed image was rotated by the repose angle and sheared into a rectangular domain so that each column of the image has the same streamwise coordinate [49], as shown in figure 2*c*. The average image intensity at each streamwise location was calculated from the image. Reference image intensities of pure heavy particles and pure light particles were used to calibrate the grey scale.

Particle tracking velocimetry (PTV) was used to determine the velocity profiles in the flowing layer. In this case, the portion of the system to be analysed extended to the surface of the flow (shifting the box in figure 2*a* upward so its top edge coincided with the surface of the flowing layer). In close-up images of steady heap flow and with proper lighting conditions, the steel particles can be identified as dark regions with small specular highlights (bright spots) on them, and the ceramic particles can be identified as white or grey spheres. This allows us to apply a Matlab-based PTV code [54] to filter noise and identify the centre positions of all the particles (figure 2*b*). Using a series of images, we computed the velocity of every particle and obtained the streamwise and normal velocity profiles at various locations along the length of the flowing layer for steady filling of a steel and ceramic particle mixture. The resulting velocity field was used to validate the simulation results. Small errors could potentially result from particles that are not at the wall but visible in gaps between particles at the wall and are therefore not well lighted, as shown in figure 2*b*. However, the misidentified particles are fewer than 5% of the total particles identified, and this error is minimized by spatial and temporal averaging, so it only causes slight fluctuation in the resulting velocity field, but not any systematic error. Thus, any possible error does not influence the comparison with DEM simulation results.

### (c) Validation of the simulations

Results from a typical DEM simulation and experiment with 3 mm steel and ceramic particles are compared in figure 3. In both cases, the particles are mixed in the inlet region on the left and become more segregated downstream. More ceramic particles flow to the end of the heap forming a region with high ceramic particle concentration. The angle of repose for the simulation (21.4°) is also similar to that for the experiment (22.1°). A quantitative comparison of light particle concentration (*c*_{l}) profiles at the bottom of the flowing layer (deposited on the heap) versus position (figure 3*c*) shows good agreement between simulation and experiment, demonstrating that the DEM simulation is able to capture the physics of bidisperse density segregation. Validation of the kinematics of the flow in DEM is described in the next section.

## 3. Kinematics of density bidisperse flow

### (a) Streamwise velocity

The streamwise velocity for the steel and ceramic particle example described in §2c in both simulations and experiments calculated using the volume average binning method [45] is shown in figure 4. For the binning, we use equal and non-overlapping bins of Δ*x*×Δ*y*×Δ*z*. In simulations, we use Δ*x*=0.01 m, Δ*y*=0.016 m and Δ*z*=0.001 m. In experiments, we use Δ*x*=0.02 m and Δ*z*=0.004 m (only particles adjacent to the wall can be observed). In the binning process, we bin the particles based on their partial volumes inside each bin in the *z*-direction, similar to the method used by Freireich *et al.* [55]. By doing this, along with time averaging over about 5 s, we are able to use small bin sizes in the *z*-direction to capture more details while still maintaining accuracy [45]. Figure 4*a* shows the free surface streamwise velocity, *u*_{s}, along the length of the flowing layer. Results from the experiment and the simulation agree well, exhibiting a nearly linear decrease along the streamwise direction, which is again consistent with a uniform deposition of particles on the heap with an approximately constant flowing layer thickness [45]. Figure 4*b* shows the flowing layer thickness *δ*(*x*) along the streamwise direction based on the streamwise velocity profile, calculated using the criteria *u*(*x*,−*δ*)=0.1*u*(*x*,0) [45]. The flowing layer thickness remains almost constant at 7–8.5 particle diameters for most of the length of the flowing layer, except near the downstream end, again consistent with previous results [45]. For simplicity in the theoretical model, a constant flowing layer depth *c* for various positions along the length of the flowing layer. The velocity profiles at different streamwise positions in the simulation agree with PTV results and collapse to a single curve, indicating a self-similar exponential velocity profile. The streamwise velocity profiles measured here are consistent with previous results for monodisperse and size bidisperse experiments and simulations [45], so the same exponential expression for the streamwise velocity is used here:
*a*, and an exponential dependence on the normal direction, consistent with the self-similar velocity profiles in figure 4*c*. Here, *k* is a scaling constant set to *u*(*x*,−*δ*)=0.1*u*(*x*,0) [38,45]. To verify the general applicability of equation (3.1), the velocity profiles at *x*/*L*=0.5 are plotted in figure 4*d* for nine simulation cases with different feed rates, density ratios *R*_{D}=*ρ*_{h}/*ρ*_{l} and particle diameters (table 2), along with the exponential fit *u*/*u*_{s}=*e*^{kz/δ}. The collapse of the data to the exponential fit demonstrates that equation (3.1) describes the self-similar streamwise velocity profiles in the density bidisperse quasi-2D-bounded heap flows studied here.

### (b) Normal velocity

Using the same method, normal velocity profiles were extracted from the simulations and the experiment. Figure 5*a* shows a comparison of normal velocities between simulation and experiment for the same case shown in figure 4*c*. Because normal velocity is typically an order of magnitude smaller than streamwise velocity, the data are more scattered. Yet there is reasonable agreement between the simulation and the experiment. In the coordinate system moving upward with the rise velocity *z*/*δ*=0) and decreases to *z*/*δ*=−1). Based on equation (3.1) and the continuity equation, the normal velocity is [45]
*x*/*L*=0.5 for the nine different simulation cases in table 2 are plotted along with *b*. The results from the simulations collapse and are quite similar to the theoretical profile, confirming that equation (3.2) is a reasonable approximation of the normal velocity profiles of density bidisperse flows in quasi-2D-bounded heaps.

### (c) Segregation velocity

Previous studies of size segregation indicate that kinetic sieving and squeeze expulsion are the dominant segregation mechanisms in gravity-driven free surface flows [13,32,33]. Here, analogous phenomena are observed in density bidisperse flows. When voids are generated due to shear, particles with higher density are more likely to fall into voids below them, while particles with lower density are more likely to be squeezed up to voids above them. Two typical examples of these processes from experiment are shown in figure 6. Figure 6*a* shows a sequence of images in which a steel particle falls into a void generated below it. Figure 6*b* shows a sequence of images in which a ceramic particle is pushed into a void above it while its original place is taken by a steel particle. An explanation of these phenomena invokes a force imbalance between the gravitational force and contact forces from neighbouring particles, such that a heavier particle on average experiences a net force in the gravitational direction and a lighter particle a net force in the direction opposite gravity. This has been referred to as ‘buoyancy’ in previous studies [20,56].

Though the segregation mechanism at work here results from density differences instead of size differences, the essence of kinetic sieving and squeeze expulsion appears to be similar, and the resulting segregation patterns for density segregation are similar to those for size segregation in quasi-2D-bounded heaps [38,40,45,49]. This suggests that the factors that drive density segregation are the same as for size segregation. These factors include the local shear rate, *c*_{i}, which determines the nature of the contact forces a particle might encounter. To examine the relation between segregation velocity and these factors, the segregation velocity *w*_{s,i} is plotted as a function of *a*. For both light particles and heavy particles, the data suggest an approximately linear relation between *w*_{s,i} and *S*_{D}, defined as the segregation length scale, is the slope of the fitted line for the dependence of *w*_{s,i} on *a*. Analogous to the segregation length scale in size bidisperse systems [38], *S*_{D} is positive for light particles and negative for heavy particles. For the data in figure 7*a*, the characteristic segregation length scales for light (*S*_{D,l}=0.150 mm) and heavy (*S*_{D,h}=−0.151 mm) particles are nearly identical in magnitude but have opposite signs due to volume conservation [38]. Although equation (3.3) has the same functional form as in size segregation [38–40], it represents the segregation velocity resulting from differences in particle densities instead of particle sizes. Thus the segregation length scale, *S*_{D}, should vary with density related properties. Values for *S*_{D} were found for 20 simulation cases using 50 : 50 mixtures with density ratio *R*_{D} ranging from 0.1 to 10 and particle diameters of 2, 3 and 4 mm at *q*=0.0022 m^{2} s^{−1}. Simulations with different feed rates were also performed, but no significant difference in *S*_{D} was found, as expected. As shown in figure 7*b*, *S*_{D} varies with the density ratio *R*_{D} and the particle diameter *d*. The relation between *S*_{D} and *R*_{D} can be approximated by
*C*_{D} is a constant with the value 0.081. 40 data points are shown in figure 7*b*, since each simulation produces two data points: one for *S*_{D,h} (corresponding to *R*_{D}>1) and one corresponding to *S*_{D,l} (corresponding to *R*_{D}<1). This functional relation differs from the assumed linear dependence of the segregation flux on the density difference *ρ*_{h}−*ρ*_{l} used in previous studies [20,24]. Equation (3.4) expresses the relation of *S*_{D} to particle differences in a manner analogous to size segregation [40], but the constant *C*_{D} is roughly three times smaller than the corresponding constant, *C*_{S}=0.26, measured in size segregation [40], which, in equation (3.4), indicates that the density ratio needs to be raised to approximately the power of three to produce equivalent segregation velocity for the same size ratio. This agrees well with the model by Tunuguntla *et al.* [24] where the segregation velocity scales as the third power of the size ratio and the first power of the density ratio, though it is inconsistent with the model by Marks *et al.* [35] where the size ratio and the density ratio are interchangeable. In an experimental study of combined size and density segregation [21], size segregation dominates density segregation unless the size ratio is much smaller than the density ratio, confirming the difference in the contribution to the segregation velocity of the size and the density differences.

For DEM simulations with a density ratio *R*_{D}≥5, there is a gradual change in the overall kinematics of the flow in the downstream portion of the heap as *R*_{D} increases. At *R*_{D}=10 (figure 8), the segregation in the upstream portion of the flow still results from local buoyancy, generating a segregation pattern with strongly segregated regions of light and heavy particles and a clear interface between. However, in the downstream portion a global flow occurs. The resistance of the light particles to the motion of the heavy particles is not significant at this high-density ratio, so the heavy particles undercut the light particles towards the end of the heap, pushing the bulk of light particles upward, such that the free surface is no longer flat. Such penetration is analogous to group intruder penetration [57], where the depth of penetration is related to the intruder's speed and density. Since density ratios this large are uncommon in industrial granular flows, and this phenomenon introduces more complexity into the kinematics of the flow, we focus on cases with density ratio *R*_{D}≤5 in this study. Note, however, that we are still able to extract local data for density segregation in the flow for *R*_{D}>5, thus accounting for these data points in figure 7*b*.

### (d) Diffusion

The diffusion coefficient of the mixture in the normal direction, *D*, was determined by tracking the non-affine portion of particle trajectories using the mean square displacement as a function of time, 〈Δ*Z*(Δ*t*)^{2}〉[38]. The diffusion coefficient was then calculated based on 〈Δ*Z*(Δ*t*)^{2}〉=2*D*Δ*t* [38,58]. An example simulation result, shown in figure 9, demonstrates that in density bidisperse flows, the diffusion coefficient is shear rate-dependent, which is consistent with previous studies in dense granular flows [20,58]. In this study, we use the spatial average of the diffusion coefficient over the entire flowing layer domain for the continuum model. We previously showed that for size-disperse granular materials, using the average value of *D* provides sufficient accuracy to successfully apply the theory, though it is possible to use a locally varying value for *D* in the theory [38].

## 4. Predictions of the theoretical model

### (a) Non-dimensionalization and boundary conditions

When applying the transport equation (1.1) to modelling density bidisperse segregation in a quasi-2D-bounded heap, it is convenient to non-dimensionalize the equation using the non-dimensionalized parameters [38,40]
*i*:
*x* are neglected as assumed previously [20,24,38,39], since these terms are small in comparison with other terms in the equation so long as *δ*/*L*≪1. The sign of the segregation term is positive for light particles and negative for heavy particles. The non-dimensional shear rate is *Pe*=2*qδ*/*DL*, which represents the ratio of a diffusion time scale (*t*_{d}=*δ*^{2}/*D*) to an advection time scale (*t*_{a}=*Lδ*/2*q*). The other non-dimensional parameter *Λ*=*S*_{D}*L*/*δ*^{2}, where *S*_{D}=|*S*_{D,l}|/2+|*S*_{D,h}|/2, represents the ratio of an advection time scale (*t*_{a}) to a segregation time scale (*t*_{s}=*δ*^{3}/2*qS*_{D}). These two non-dimensional parameters depend on control parameters (feed rate *q* and flowing layer length *L*) and kinematic parameters (flowing layer thickness *δ*, diffusion coefficient *D* and segregation length scale *S*_{D}), and they represent the interplay of advection, diffusion and segregation. Boundary conditions are also identical to previous studies for size bidisperse flow [38,40]. At the inlet, the particles are well mixed, so *w*=−*v*_{r} cos *α*, and no particles leave the flowing layer at the top surface. At the downstream boundary, advection, diffusion and segregation included in equation (4.2) are in the normal direction. Thus, no boundary condition is needed. Here, we note that with our method of informing the model with kinematics measured from simulations, the model itself does not include the feedback between velocity field and segregation as was considered by Marks *et al.* [35]. However, we are solving for steady state solutions where the kinematics do not change, and the kinematics are measured directly from the steady state bidisperse flow simulations. As a result, the kinematics included in our model represent the final result of the interaction between segregation and kinematics. With the velocity profiles (equations (3.1), (3.3)), equation (4.2) can be solved numerically for steady state flow using an operator splitting method with a mapping method for the advection step and the implicit Crank–Nicolson method for the diffusion and segregation step [38,59,60].

### (b) Validation of the theoretical model

To demonstrate that the theoretical model quantitatively predicts bidisperse density segregation in the quasi-2D-bounded heap, theoretical, experimental and simulation results for three example cases with different particle density ratios *R*_{D} are shown in figure 10. Each contour subplot (left and centre columns) represents the concentration of light particles *c*_{l} in the flowing layer extending horizontally from *c*_{l}, deposited onto the heap at the bottom of the flowing layer [38], see right column in figure 10. The root-mean-square deviation (RMSD) is applied to quantify the quality of the theoretical predictions and is calculated as *c*_{j,t} denotes the data points from the theoretical predictions and *c*_{j,i} denotes corresponding data points from the reference (*i*=*e* for experiments and *i*=*s* for simulations), and *n* is the total number of data points in each case. In all three cases, the RMSD_{t,e} and RMSD_{t,s} are small, indicating that the theoretical predictions match the experimental and simulation results.

The theoretical predictions are determined completely by the two dimensionless parameters *Pe* and *Λ*, in the same manner as in previous work [28,33]. *Pe* describes the interplay between advection and diffusion: as *Pe* becomes larger, advection dominates diffusion, causing the interface between segregated heavy and light particles to become sharper and more easily distinguishable. We note that the *Pe* we use here is different from the Péclet number defined in Thornton *et al.* [33] for chute flows, where it represents the interplay between segregation and diffusion. *Λ* describes the interplay between segregation and advection. For larger *Λ*, segregation is stronger so particles tend to segregate before they flow very far down the heap. The influence of these two parameters on segregation has been investigated in detail in the context of size bidisperse systems [38]. Since the form of the theoretical model here is identical to that in previous work for size segregation [38,40] (except that the percolation length scale *S* is replaced by the density segregation length scale *S*_{D}), the discussion is not repeated here.

### (c) Predictions of segregation under different physical control parameters

Since *Pe* and *Λ* depend on physical control parameters, it is interesting to explore how theoretical predictions of segregation change when the physical control parameters are varied. Among the parameters, density ratio *R*_{D} and feed rate *q* have the greatest influence on the segregation [40]. Figure 11 shows a series of theoretical predictions and simulation results for the concentration of the light particles at the bottom of the flowing layer (which are deposited on the heap) for different *R*_{D} and *q*. In all cases, the diffusion coefficient, *D*, is the average diffusion coefficient based on the simulation for each set of conditions and *S*_{D} is from the relation shown in figure 7*b*. In figure 11*a*–*d*, it is clear that increasing *q*, which results in *Pe* increasing and *Λ* decreasing, results in less segregation (the transition from mixed particles to pure light particles occurs further downstream), and a sharper transition between segregated heavy and light particles, which is made more readily apparent in the figures with the aid of a vertical dashed benchmark line. In figure 11*e*–*h*, increasing *R*_{D}, which results in *Λ* increasing while *Pe* varies only a small amount, leads to an obvious increase in segregation; there is almost no segregation with *R*_{D}=1.11 but strong segregation with *R*_{D}=3.33. In all cases, RMSD_{t,s} is small, indicating that the theoretical predictions match the simulation results well, again demonstrating that the theory is capable of accurately predicting segregation when the physical control parameters are varied.

To further demonstrate the generality of the theory and the form of the segregation velocity, cases with different inlet particle concentrations were simulated and predicted using the theory, as shown in figure 12. In both cases, *D* is from simulations with *S*_{D} is from figure 7*b*, rather than using these values from the simulations for the *S*_{D}=0.27 mm), the theoretical predictions of the concentration distribution of the light particles in the flowing layer (row 2 in figure 12) agree reasonably well with the simulation results (row 3). The comparison of the light particle concentration profiles at the bottom of the flowing layer from simulations and theory (row 4) demonstrates that the theory is also accurate for particle mixtures with different volume ratios, even when using parameters obtained from a 50 : 50 mixture.

While examining the effect of different volume ratios on the segregation, we found a difference in the segregation length scale *S*_{D} calculated in the simulations with *S*_{D} calculated in each simulation, we recalculated *Λ* (*Pe* is independent of *S*_{D}): *S*_{D}=0.34 mm, *Λ*=0.28 for *S*_{D}=0.25 mm, *Λ*=0.21 for *S*_{D}=0.27 mm, *Λ*=0.23 for *Pe* and *Λ* recalculated in the two cases are also shown in figure 12*g*,*h* as dashed curves. For *S*_{D} from equation (3.4). For *S*_{D} determined for

The difference in *S*_{D} for different inlet concentrations is intriguing, because it indicates that the segregation for a few heavy particles in many light particles is stronger than the segregation for a few light particles in many heavy particles. This is analogous to recent work which shows that small particles segregate faster when surrounded by large particles than vice versa [61,62]. This asymmetry can possibly be explained in terms of the way that a heavy (or small) particle is able to continually push its way downward in the gravitational direction while waiting for a void below it to open when it is surrounded by light (or large) particles. By contrast, a light (or large) particle can only wait for the combination of a void opening above it at the same time as surrounding particles are pushing it upward against gravity when it is surrounded by heavy (or small) particles.

This asymmetry suggests that *S*_{D} depends on local particle concentrations. To further explore this, nine simulations with *R*_{D}=3.04 and *q*=0.0027 m^{2} s^{−1}, and *d*=2 mm were performed with the inlet light particle concentration varying from 0.1 to 0.9. The segregation velocity shown in figure 13 includes data from all nine simulations. To reduce the noise, we averaged the data into 100 equal-sized bins along the horizontal axis, and these averaged results are also shown in figure 13 as two yellow curves. The data for different inlet light particle concentrations are consistent with each other, forming two continuous curves (one for light particles and one for heavy particles). Upon closer examination, the curve for heavy particles, when compared with the curve for light particles, has a slightly smaller slope when 1−*c*_{i} is small, corresponding to a high local concentration of heavy particles (case A in the figure), and a slightly larger slope when 1−*c*_{i} is large, corresponding to a lower local concentration of heavy particles (case B in the figure). This asymmetry in density segregation reveals the reason for the difference in *S*_{D} measured in cases with different inlet light particle concentrations: *S*_{D} for *S*_{D} measured for *et al.* [35] and Tunuguntla *et al.* [24] should automatically include this nonlinearity, but its effect was not discussed in detail. Recent studies on asymmetry in size segregation used a cubic function to represent the asymmetric flux [61,63] resulting in better agreement with experiments [61]. For density segregation, this asymmetry could also be taken into account in the theory by using an *S*_{D} that is a function of *c*_{i}, which we are currently investigating.

## 5. Conclusion

In this study, we have demonstrated that our recently developed continuum model for size bidisperse systems [38–40] accurately predicts granular segregation for density bidisperse systems, specifically for bounded heap flow, though it is likely applicable to other flow geometries and multi- or polydisperse particle distributions, as we have already shown for size segregation [39]. Using experimental techniques and DEM simulations to investigate the kinematics, we developed an approximation for the segregation velocity that depends on local shear rate *c*_{i}, and a dimensional parameter defined as the density segregation length *S*_{D}, which depends primarily on the density ratio *R*_{D} and, to a lesser extent, on the local particle concentration. The model is based on the transport equation and includes the interplay of advection, diffusion and segregation. In the model, no arbitrary fitting parameters are needed as the system configurations are determined by two dimensionless parameters *Pe*=2*qδ*/*DL* and *Λ*=*S*_{D}*L*/*δ*^{2}, which depend only on physical control parameters and kinematic parameters measured from simulations (or experiments, if available). The theoretical predictions quantitatively agree with results from both simulations and experiments under different physical control parameters.

Compared with our model for size segregation [38,40], the primary difference in this study is the segregation length scale in the equation for the segregation velocity, which is related to the density ratio here, as opposed to the size ratio for size segregation. However, apart from this constant, the two models are identical and the resulting segregation patterns for size and density segregation are very similar. This suggests that although the driving force for granular segregation is different in the two cases, the shear-generated segregation mechanism for gravity-driven free surface flows is similar. This also suggests that the model has potential to predict combined size and density segregation [24,36]. The asymmetry observed in density segregation for different inlet concentrations indicates that more accurate predictions will require that variations of *S*_{D} with concentration be included in the theoretical model. Moreover, the model is not limited to quasi-2D-bounded heaps. For other typical geometries with gravity-driven free surface flows, the continuum model combined with the segregation velocity defined in equation (3.3) should also be applicable. To adapt the continuum model to segregation in other geometries, similar to studies of segregation in rotating tumblers [14,64], the corresponding boundary conditions need to be included and kinematic information for the specific geometry needs be used. This can be readily accomplished using constitutive theories [20,35] or direct measurements from simulations or experiments [38,39]. This approach has already been used to successfully model size segregation in circular tumblers [39] and chutes [41]. Given the similarities between size and density segregation, it is quite likely that the model will accurately predict density segregation in tumblers, chutes, unbounded heaps and even three-dimensional geometries.

## Data accessibility

Data supporting this article are available upon request within 5 years of publication of the article from Dr Richard M. Lueptow (r-lueptow{at}northwestern.edu).

## Authors' contributions

H.X. carried out the experiments, DEM simulations, and the theoretical calculations and drafted the manuscript. P.B.U., J.M.O. and R.M.L. designed the study, offered suggestions in discussions and edited the manuscript. R.M.L. supervised the project.

## Competing interests

We declare we have no competing interests.

## Funding

We received no funding for this study.

## Acknowledgements

The authors thank Austin B. Isner for the use of his GPU-based DEM simulation code along with Ben J. Freireich, Yi Fan and Karl Jacob from the Dow Chemical Company and John Hecht and Vidya Vidyapati from the Procter and Gamble Company for helpful discussions.

- Received December 15, 2015.
- Accepted June 10, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.