Theory for shear-induced segregation of dense granular mixtures

It is well known that a mixture of different sized particles will segregate in a gravitational field. However, it has only recently been shown that a gradient of shear rate alone can drive segregation in dense sheared systems. In contrast with sparse energetic granular materials, in dense sheared systems, large particles segregate to regions with higher shear rates. In this paper, we develop a model for shear-induced segregation in dense mixtures of different sized particles. The model is comprised of two primary parts. The first involves the tendency of a gradient in kinetic stress—stress associated with velocity fluctuation correlations—to drive all particles toward regions of low shear rate. The second is essentially a kinetic sieving effect in which small particles are more likely than large particles to find voids into which they can travel. The two features together segregate small particles to regions of low shear rate and squeeze large particles in the opposite direction. We validate this model via three-dimensional discrete element method simulations in a vertical chute.


Introduction
Segregation of sheared dense granular mixtures is ubiquitous in both natural and industrial applications [1][2][3][4][5][6]. In some cases, segregation is undesirable, as in pharmaceutical industries, as the tendency for particles to segregate by property can hamper the mixing of pharmaceutical powders and consequently lead to inconsistent delivery of medication [1]. In those cases, the study of segregation is driven largely by a desire to inhibit it. In certain geophysical processes, understanding the specific mechanisms and rates of segregation would be helpful in predicting natural processes as in sediment transport in rivers [2][3][4] and erosion in upland regions associated with debris flows [5,6]. However, modeling segregation in dense sheared flows has proven difficult.
It is well known that a mixture of different sized particles tends to segregate when sheared either by the mechanical movement of solid boundaries (e.g. in an annular shear cell [7,8]) or by a component of gravity (e.g. in inclined plane flows [9][10][11] and in rotating drums [12][13][14][15][16][17][18][19][20][21][22]). Although we have cited several studies in this paper, over the years segregation in these sheared systems has presented such a range of interesting problems that the studies are simply too numerous to list here. The focus of modeling efforts regarding these segregation phenomena has largely been on the segregation associated with the driving force of gravity. In most cases, large particles in the granular mixtures will rise in the opposite direction from gravity and small particles sink, often referred to as 'reverse grading' in certain geophysical problems [23,24]. However, there are some exceptions in which large particles sink when we might expect them to rise (e.g. [19,21,25]). Although these tendencies appear to vary systematically in terms of their 3 dependence on particle properties such as size and density, there is no predictive physics-based model for this variation.
Several mechanisms have been proposed for gravity-driven segregation of high solid fractions of different sized particles (e.g. [9-11, 26, 27]). However the process is often associated with processes known as kinetic sieving and squeeze expulsion as in [9]. Essentially, gravity pushes all particles in one direction, but the structure associated with the high packing fraction prevents most particles from responding; the movement of particles in the direction of gravity is limited by the lack of sufficiently large 'holes' into which they may drop. Jostling associated with shearing between layers of particles opens up holes of variable sizes. Statistically, in a mixture of different sized particles it is more likely that small particles find holes sufficiently large for them to fall through these 'fluctuating sieves' via what is referred to as a kinetic sieving mechanism. Mass balance is achieved when the larger particles are pushed by these smaller particles migrating downward by what is sometimes referred to as a 'squeeze expulsion' mechanism. Savage and Lun [9] first developed a detailed statistical model to predict segregation trends based on this process. More recently, Gray and Thornton [10] developed a continuum framework for this process based in part on binary mixture theory. Essentially, the different ways in which different sized particles respond to a pressure gradient is represented in this model by the relative proportion of the pressure associated with gravity that each species bears compared to its local concentration. The framework extended to include diffusive mixing effects [11] has been shown to be effective in predicting segregation trends and other exciting dynamics [28] in granular flow.
In contrast with gravity-driven segregation in dense flows, relatively little attention has been given to the accompanying effect of a shear rate gradient, particularly in systems of higher solid fractions. The absolute magnitude of shear rate appears in some models for segregation in dense sheared granular flows (e.g. [9]). However, it appears only in the rate at which segregation occurs rather than in how its gradient might actually modify segregation trends. The effects of shear rate gradients on segregation have generally not been considered in segregation models for dense systems.
In related systems, a shear rate gradient has been shown to be responsible for segregation of different sized particles, particularly in sparser dry systems and in particle suspensions. For example, Conway et al [29] demonstrated that sparse mixtures of different sized particles sheared in an experimental Couette cell could be segregated purely by the dynamics associated with shear rate gradients. In their experiments the shear rate gradients gave rise to gradients in velocity fluctuations (commonly associated with a 'granular temperature'), and large particles segregated to the region of lowest shear rate and granular temperature. Leighton et al [30][31][32] and Abbott et al [33] showed that when concentrated particle suspensions were sheared in a circular Couette cell, all particles migrate to regions of lower shear rate, while larger particles segregate to the region of lower shear rate faster, similar to the case of sparse dry granular mixtures. These trends have been modeled effectively for relatively sparse mixtures of dry particle mixtures using kinetic theory (e.g. [34,35]). However, for denser granular systems, kinetic theory has been shown to be somewhat less effective [36,37] and, at times, to predict a trend opposite to that observed [38].
In dry dense granular materials, there is experimental and computational evidence that, in contrast with sparse granular materials and suspensions, large particles segregate toward the region of highest shear rate, also corresponding to regions of high granular temperature. In the 1970s, Bridgwater and colleagues performed some of the earliest experiments on this 4 phenomenon in an annular shear cell [7] and a reciprocating shear cell with flexible end walls [39]. They focused on the effect of the shear rate gradient on the segregation and found that large particles move to the region of highest shear rate in both cases. They did not report on variations in local solid concentrations and velocity fluctuations.
More recently, we performed simulations of dense sheared granular mixtures in a vertical chute to investigate how segregation according to gradients in shear rate varies with solid fraction in this system [38]. While the vertical chute does not eliminate gravity entirely, it provides a system where shear rate gradients are normal to gravity and advection does not appear to mix up the effects associated with the two. In a vertical chute with roughened walls, the shear rate and granular temperature are highest adjacent to vertical walls for a wide range of solid fractions (figure 1). We found that for relatively sparse granular mixtures, the large particles segregate away from the walls to the low-temperature, low-shear-rate region of the chute, similar to the results reported in [29,35]. For higher solid fractions, we found that the trend reverses and, similar to the results for other dense sheared granular mixtures [7,37,39,42], the large particles segregate towards the walls to the regions of high shear rate and high granular temperature (figure 2). We found that the transition in segregation trends coincides with the formation-at higher solid fractions-of a large-scale cluster of particles connected via interparticle contacts.
Based on these observations, we hypothesize that a combination of two mechanisms is responsible for the segregation of larger particles to regions where the shear rate is low. Firstly, the dynamics associated with the shear rate gradients pushes all particles toward regions of low shear rate and granular temperature, analogous to the way in which gravity pushes all particles downward. Secondly, the system-wide cluster of particles that forms in the higher solid fractions systems gives rise to dynamics associated with kinetic sieving. In its original form, the kinetic sieving mechanism segregates particles along pressure gradients associated with gravity. We hypothesize that in dense systems, where gravity plays no role in inducing a pressure gradient, the kinetic sieving mechanism segregates particles along pressure/stress gradients induced by the shear rate gradients forcing large particles to the walls and small particles to the middle of the chute.
In this paper, we consider our hypotheses in the context of a model for gravity-driven segregation of binary mixtures developed by Gray and colleagues [10,11] described briefly above. Compared with the model by Gray and colleagues involving a hydrostatic pressure gradient associated with gravity, we consider the effects of stress gradients arising from the shear rate gradients. In section 2, we show theoretically how the separate consideration of what have been called contact and kinetic stresses elucidates the mechanism of shear-induced segregation. In section 3, we present simulation results that we use to test the model. Specifically, in section 3.5 we suggest some simplifications that lead to a predictive form of the model and then compare the results with simulations. In section 4, we discuss some of the effects of the granularity of the system and the associated limitations on the continuum-style modeling we employ. We conclude and discuss future directions in section 5.

Theory
Our model concerns a binary mixture of different sized spherical particles with the same material density ρ m . We denote bulk Eulerian properties of each species with superscripts and those of the mixture of both species together as variables without superscript. For example, (b-f) Profiles of kinematic quantities in the y-direction for a 50/50 mixture of 2 and 3 mm particles at steady state (quantities averaged from t = 950-1000 s after the time at which the particles were initially released). The average solid fraction of the mixture f = 0.6. Data are calculated at spatial intervals in the y-direction of 2 mm, so the average trends are apparent. Higherresolution data analysis reveals structures related to layering of particles near the walls as will be discussed in section 4. (b) Solid fraction of the mixture f . (c) Streamwise velocity w. (d) Absolute value of shear rate |dw/dy|. (e) Granular temperature T = (u u + v v + w w )/3 and three components u u , v v , w w . (f) Averaged kinetic normal stress ρv v . In (e), the spatial resolution of the data is the same as that for the other plots (2 mm), although lines are used for clarity in distinguishing data sets. In (c), the line represents the linear least squares fit of w = a + b * (1 + cosh(y/c)) to the data, where (a, b, c) = (−0.3753 (m s −1 ), 0.0002207 (m s −1 ), 3.063 (mm)). This is similar to the profile found to fit the long-term average vertical motion of particles in vertically shaken containers in experiments performed by Knight et al [40] and in a stochastic model for void transport by Shinbrot et al [41]. we represent the local volume fraction filled with particles of species i by f i ; therefore, the local bulk mass density of species i, ρ i = ρ m f i . We denote the concentration of species i by The average local density of the mixture ρ = ρ i and the average local velocity u = u i φ i .

Conservation of mass and momentum for dense sheared mixtures
We first consider conservation of mass and momentum for the mixture: We denote the stress tensor by σ , using the relatively standard sign convention for stresses as, for example, noted in [43], and we denote body forces (per unit volume) F. We consider systems in which the mixture velocities and densities reach the steady state long before segregation. With this in mind, we set the temporal derivatives ∂ρ/∂t and ∂(ρu)/∂t to zero. We then consider the instantaneous value of each variable q at position r as a sum of the local temporal average q(r) and the difference between its instantaneous value and the average q (r, t) = q(r, t) − q(r) 7 (typically called 'Reynolds decomposition'; [44]). Considering this, we rewrite the momentum equation (2) for the jth component as As detailed in the appendix, we take the average of both sides of equation (3), use basic identities such as q = q, q = 0 and p + q = p + q and assume that temporal correlations between velocity fluctuations and densities are negligible, i.e. ∂ ∂ x k (u j ρ u k ), ∂ ∂ x k (u k ρ u j ) and ∂ ∂ x k (ρ u k u j ) ≈ 0. We consider the results in the context of pseudo-2D systems such as the vertical chute flow shown in figure 1 in which segregation occurs in a direction normal to the average flow (e.g. the y-direction), the flow exhibits uniformity in the other two directions (e.g. the xand z-directions) and the only body force (particle weight) is in the z-direction. Then equation (3) in the y-direction may be simply expressed by In this equation, −ρv v is often referred to as a Reynolds stress component (e.g. [44,45]) and ρv v is sometimes referred to as a streaming stress component (e.g. [45]) or a kinetic stress (e.g. [46]). We follow the terminology of [46] and denote σ k yy = ρv v as a kinetic stress. We also note that for dry macroscopic particles, the normal contact stress will be solely compressive, whereas in the standard sign convention for σ (e.g. [43]) normal stresses such as σ yy are positive only for tensile stress and negative for compressive stress. We therefore define a contact stress tensor σ c = −σ so that terms such as σ c yy are positive for our problem. Then we rewrite equation (4) as which indicates that the horizontal gradients in contact and kinetic stresses are exactly equal and opposite in sign to one another for the mixture. We note that ρv v scales roughly with T (figures 1(e) and (f)), so that equation (4) indicates that a gradient in T can be associated with a gradient in both σ k yy and σ c yy (of opposite signs).

Conservation of mass and momentum for species in a mixture and subsequent segregation
Next, we investigate the conservation laws for two species: Here, all of the terms used for species i are identical in form to those used for the mixture in equation (2), with the exception of the new term β i , which represents the interaction force exerted on species i by the other species. We note that σ i is the local stress borne by species i and that the total stress σ = σ i . By performing Reynolds decomposition and then averaging both sides of equation (7), we may use similar arguments to those we employed for the mixture, 8 as detailed in the appendix. With this, we rewrite equation (7) in the y-direction as We note that all terms in equations (5) and (8) are averaged and drop the overbar from this point forward, so that unless noted for each variable q refers to the average q.
The discussion above reflects much of the notation from mixture theory (e.g. [47][48][49]). However, there are a few exceptions, which we note here. In mixture theory the velocities of the constituents are assumed to be equal to the bulk velocity (u i = u), although this cannot be the case for segregating mixtures, so we allow u i to vary from u and set u = u i φ i as previously noted. Along these lines, as demonstrated in [20], the velocity fluctuation correlations that comprise the kinetic stress terms are different for different sized particles in a mixture, and typically |u s k u s j | > |u l k u l j |. (Here, and from this point forward, superscripts s and l refer to quantities associated with the small and large particles, respectively.) Additionally, mixture theory assumes that other fields are related by linear volume fraction scalings (e.g. for constituent i, mixture theory states that σ i = φ i σ ). In contrast, Gray and Thornton [10] proposed that for mixtures of different sized particles subjected to gravity, the associated lithostatic pressure for small and large constituents will not partition with linear volume fraction scalings. Considering these previous studies, we suggest that the partitioning of kinetic and contact stresses between the large and small particles is not equal to the associated solid fraction (σ c,i yy = φ i σ c yy and σ k,i yy = φ i σ k yy ). Instead, we propose that, similar to the model for partial pressures in a gravity-driven system proposed by Gray and Thornton [10], where ψ c,i and ψ k,i determine the proportion of normal contact and kinetic stresses carried by species i and are not necessarily equal to φ i . To ensure that i σ c,i yy = σ c yy and i σ k,i yy = σ k yy , two constraints on ψ c,i and ψ k,i must be satisfied: (i) ψ l + ψ s = 1 and (ii) if only one species is present, then it must support the entire local stress, so that when φ i = 1, ψ j = δ i j , where δ i j is the Kronecker delta function.
For the interaction term β i y , we propose a similar form to that proposed by Gray and Chugunov [11]: following [10,11]. The first two terms on the right-hand side of the equation ensure that, as in Darcy's law, the percolation process (in this case, the kinetic sieving process) is driven by intrinsic rather than partial stress gradients. The third term is a linear drag law similar to that provided by Morland [49] for the percolation of fluids-c D is the linear drag coefficient. The fourth term acts as a 'remixing force' (e.g. [11] and [49]) that drives grains of constituent i towards areas of lower concentration-d is the ordinary diffusion coefficient. Combining equations (9) and (10) with equation (5), we find that we can express a segregation flux of species i as where, to facilitate a physical interpretation of the governing features of this equation, we have introduced partitioning variables R c,i = ψ c,i /φ i and R k,i = ψ k,i /φ i to represent the ratio 9 between the fraction of the local stress (both contact and kinetic) carried by each species and its local concentration in the mixture. From this equation, we see that if R c,i , R k,i = 1, that is, the local stress of both types borne by each species is equal to its local concentration, there should be no segregation. If this is not true of both contact and kinetic stress, but R c,i = R k,i , still equation (11) indicates that there should be no segregation. However, if R c,i = R k,i , equation (11) indicates that, in the context of a nonzero stress gradient (∂σ k yy /∂ y = 0), segregation should occur. Further, the direction of segregation flux is indicated by the sign of (R c,i − R k,i )∂σ k yy /∂ y.

Segregation by shear compared with segregation by gravity
We now compare our model with that developed from the related theory for segregation associated with gravity by Gray and colleagues [10,11]. In their case, they considered shallow flows driven by gravity, down an inclined plane whose angle of inclination with the horizontal is ζ . For the discussion here, we letz be the direction perpendicular to the free surface. Since this is the direction of segregation, thez-direction in gravity-driven segregation is comparable to the y-direction in our discussion above; for both cases, the flow is assumed to be uniform in all other directions. In their development, Gray and colleagues assumed that the acceleration terms on the left-hand side of equations (2) and (7) are negligible and that F is the specific weight of the material, so that conservation of momentum for the mixture (equation (2)) in thez-direction may be written as This is comparable with our equation (4) where σ c zz is analogous to p in equation (12). They showed that their conservation of momentum equation in thez-direction for mixture species i (compared with equation (8)) could be written as where p i = ψ p,i p is the amount of local pressure borne by species i (ψ p,i = p i / p), is the term representing the interaction between the species. w i and w are the velocities in thẽ z-direction of species i and the mixture, respectively. As detailed in the appendix, these can be used to derive a segregation flux comparable with equation (11): where R p,i ≡ ψ p,i /φ i . We now compare the segregation associated with gravity predicted by equation (15) and the segregation associated with a shear rate gradient predicted by equation (11). In both cases, there is an apparent physical force that drives all particles in one direction. In the case of gravity alone (disregarding any shear rate gradients and associated kinetics), gravity pulls all particles downward, away from the free surface and in the absence of resistive forces they simply fall. In the case of a shear rate gradient alone, an associated gradient in kinetic stresses drives all particles along the kinetic stress gradient to lower values of σ k yy and T . In both cases, when the particles encounter one another, one constituent is driven more effectively along the relevant gradient than the other. The process by which this occurs for dense sheared systems may be seen as a kinetic sieving process allowing some particles (the smaller particles) to pass more freely through the matrix created by a dense arrangement of particles. In the case of segregation associated with gravity alone, equation (15) suggests that whichever constituent supports a greater fraction of the pressure than their representative concentration (i.e. R p,i > 1) will segregate upwards. In the case associated with a shear rate gradient alone, equation (11) suggests that whichever constituent bears more of the kinetic stress than the contact stress (R k,i > R c,i ) will segregate in the direction of decreasing values of σ k yy and T . Since large particles rise toward the free surface in gravity-driven flows, equation (15) suggests that R p,l > 1 and R p,s < 1, that is, that ψ p,l > φ l and ψ p,s < φ s , as detailed first in [10]. In dense sheared systems in situations when gravity can be neglected, larger particles segregate to high-T regions, essentially regions of high kinetic stress, so equation (11) suggests that R c,l > R k,l and R c,s < R k,s .

Summary of theory
In this section we demonstrated that in relatively simple sheared systems, in the direction of a shear gradient and no body force, a gradient in contact stresses exactly balances a gradient in kinetic stresses, as in equation (5): Whichever species bears a greater fraction of the kinetic stress than it does the contact stress (R c,i < R k,i ) moves toward lower kinetic stresses, as indicated by equation (11): Likewise, whichever species bears a greater fraction of the contact stress than it does the kinetic stress moves toward lower contact stresses. The segregation is limited by diffusion. This framework is in many ways analogous to a gravity-driven segregation model for free surface granular flow derived by Gray and colleagues. In that case, a local gradient in pressure in the direction normal to the free surface is associated with the specific weight of the material, as in equation (12): In that case, whichever species bears a greater fraction of the pressure than its representative fraction in the mixture (R p,i > 1) moves toward lower pressures, as indicated when equation (12) is substituted into (15): Again, this segregation is limited by diffusion. We test these predictions from equation (11) as well as the assumptions used to derive equation (11) using the data from vertical chute flow simulations, as detailed in the next section. Values of the numerical parameters used in the DEM simulations. The parameters were derived for particles having similar properties to glass (for details, please see, e.g., [52]), although to reduce the computational time we reduce the normal stiffness k n by a factor of O(10 3 ) and the material density of the particles ρ by a factor of 10. (Sample simulations with and without reducing k n and ρ result in similar segregation concentrations and other details.)

Discrete element method simulations
To investigate the assumptions and predictions associated with the theory in section 2, we investigate in detail the discrete element method (DEM) simulation results in a chute flow described briefly in the introduction (figures 1 and 2). As is typical in DEM simulations [50], we track each particle throughout each simulation, and therefore every force on each particle is obtainable at each time step. To model the forces between contacting particles, we use a nonlinear force model based on Hertzian and Mindlin contact theories [51], damping components based on experimental data as outlined in [52] and the Coulomb law of friction. The interparticle forces F = F n + F t each have components normal (F n ) and tangential (F t ) to the plane of contact: In these equations, we follow convention and model deformations resulting from interparticle contact as effective overlap in the directions normal and tangential to the plane of contact, δ n and δ t , respectively. In all cases in equation (20), n and t refer to directions normal and tangential to the plane of contact, respectively. V n = (dδ n /dt)n and V t = (dδ t /dt)t are relative velocities of contacting particles. n and t are unit vectors in each direction. k n , k t , η n , η t and µ are the interaction coefficients derived from material properties as described in [51] and [52]. Sample interaction coefficients for contacts in a mixture of 2 and 3 mm spheres can be found in table 1. We calculate these based partly on the material properties of glass, although to reduce computational time we reduce the stiffness and density as described in the caption of table 1. For the results described here, we simulate a 50/50 mixture by weight of 2 and 3 mm spheres with a 10% polydispersity to impede crystallization.
The boundary conditions for our simulations are those of a vertical chute of dimensions D = 40 mm, W = 50 mm and L = 50 mm in the x, y and z directions, respectively (see figure 1(a)). Our chute has one pair of vertical side walls (perpendicular to the y-direction), which are roughened using 2 mm spheres in a random close-packed arrangement. The boundaries are periodic in the z (vertical) and x directions. Although we perform simulations for several different solid fractions, for the results described here we focus on simulations for a system solid fraction f = 0.6, using a total of ≈ 8000 particles in the simulation. As in section 1, we denote the velocity u = ux + v y + wz according to the directions denoted in figure 1(a).
For each simulation, the particles are initially arranged randomly in the chute and then released with small random velocities. After their initial release, particles collide with one another and with the vertical walls as they initially accelerate downward. Dissipation of energy through interparticle and wall-particle interactions limits the velocity throughout the cell. As we noted in the introduction, while the vertical chute does not eliminate gravity entirely, it provides a system where shear rate gradients are normal to gravity and advection does not appear to mix up the effects associated with the two. We focus primarily on the dynamics as they vary in the y-direction, accordingly.

Calculations of kinematics and stresses
To calculate the relevant kinematics for this paper, we divide the chute into equal sized bins in the y-direction of width y. We calculate kinematic quantities by considering the contribution of each particle j within a bin of width y centered at y, as illustrated in figure 3. For example, for f n (y) and v n (y), the solid fraction and the velocity associated with species n are and v n (y) Here, i refers to the ith time step of which there are N , and j refers to the jth particle (of species n) that is partly or fully in this bin (at time step i). V n i j and v n i j are the volume and velocity of that particle, respectively. V bin is the total volume of the bin, V bin = DL y. A sketch illustrating parameters needed to calculate contact stresses associated with the contact between particles i and j. Note that l i j is the vector from the center of particle i to the center of particle j. F i j is the force of particle i exerted on particle j.
Similarly, we calculate σ k,n yy (y) = ρ n v n v n (y) (the y-component of the normal kinetic stress of species n) according to where v(y) = n ( f n (y)v n (y))/ n f n (y). As in section 2.2, σ k yy (y) = n σ k,n yy (y). To calculate the local contact stress at each position y, we consider each interparticle contact K in a bin of width y centered at y ( figure 4). Then, we sum the stresses associated with each interparticle contact in each region, as in [45,53]. Specifically, for the mixture we calculate Here, F i j K is the force of particle i on particle j associated with the K th contact in this bin, of which there are N c (y, τ ) at time step τ . There are N such time steps. l i j K is the vector from the center of particle i to the center of particle j ( figure 4). We calculate the component of interest of the contact stress tensor according to Compared with the calculation of partial kinetic stresses, which is relatively straightforward, calculating the partial contact stresses for each species takes a little more care, since a contact may involve particles of different species. We consider three types of contacts separately in calculating the species contact stresses. (i) Contacts between two small particles only contribute to the partial contact stress of the small particles, and we denote the stress associated with the K th such contact as σ c,ss K . (ii) Contacts between two large particles only 14 contribute to the partial contact stress of the large particles, and we denote the stress associated with the K th such contact as σ c,ll K . (iii) Contacts between one small particle and one large particles contribute to the contact stress of both species; we denote the stress associated with the K th such contact as σ c,sl K . There are many ways one could divide each instance of σ c,sl K between the partial stress of the species. We note that both the magnitudes of the interparticle forces acting on each of the two particles in contact are equal and the contact area over which the force is distributed is equal. Therefore, for each collision between a small and a large particle, we divided the contribution of stress to the partial stresses equally between the two species. Based on that, we calculate the partial contact stress at y for particles of species n at time step τ as In this equation, σ c,nn K denotes the contact stress associated with the K th contact between a particle of type n with another particle of the same species in a bin of width y centered at y, of which there are N i c (y). σ c,n j K denotes the contact stress associated with the K th contact between two particles of different species in a bin of width y centered at y, of which there are N j c (y). We calculate the average stress over N time steps: We note that this satisfies σ c (y) = σ c,l (y) + σ c,s (y), as specified in section 2.

Average kinematics
Upon evaluating the kinematics using a relatively coarse resolution (bin size y = 2 mm), we obtain the steady-state kinematic profiles shown in figure 1, which we summarized briefly in section 1. We discuss these results in a little more detail here. These steady-state kinematic profiles of the granular mixtures are similar to those in a mono-sized system measured in physical experiments and other simulations (e.g. [46,[54][55][56]). At the highest solid fractions with which this paper is concerned, f is roughly uniform, although near the walls the value of f decreases slightly ( figure 1(b)). The vertical velocity profile w(y) resembles that of a plug flow, and the shear rateγ = dw/dy is largest near the side walls (figures 1(c) and (d)). Both the granular temperature T ( figure 1(e)) and the component of the kinetic stress of concern σ k yy = ρv v (figure 1(f)) are highest adjacent to the walls, whereγ is the greatest (compare figures 1(d)-(f)).   In other words, the segregation patterns apparent in figures 2(b) and (c) are directly associated with these horizontal segregation fluxes rather than in other cases in which there is advection, which essentially 'rotates' the segregation pattern relative to the direction of the segregation flux (e.g. [37]).
To determine how the rate of segregation compares with the rate at which the kinematics of the mixture reaches the steady state, we consider two quantities. The first is the velocity averaged over the width of the chute, w ( figure 5(a)). The second is a measure of the segregation S as a function of time ( figure 5(b)). For the latter, we calculate the standard deviation of mean concentration S i of each species i at each time step (t): where N bin = 2500 is the number of bins in the y-direction, φ i j is the concentration of species i in bin j and φ i is mean concentration of this species in the system (0.5 for both species). Since φ l = φ s = 0.5, and φ l j (t) + φ s j (t) = 1, S l = S s , which we denote by S. In figure 5 we fit both w and S by exponentials described by f (t) = A − B exp(−t/τ ), where A, B and τ are fitting parameters given in the caption of figure 5. The results indicate that the mixture kinematics reaches near-steady-state conditions in less than 10 s, although segregation takes somewhat longer.

Average contact and kinetic stresses
In figures 6(a) and (b), we have plotted profiles of σ k yy (y) and σ c yy (y) averaged over the 50 s interval in the simulation after the mixture kinematics is essentially at steady state, but the segregation dynamics is still active. The profile of σ k yy (y) peaks near the rough walls and dips in the middle. As shown, it is well fitted by an exponential function σ k yy (y) = A exp(B|y|) whose fitting parameters A and B are given in the caption of figure 6. As one would expect from equation (5), σ c yy (y) follows the opposite trend: it is highest in the middle and dips near the walls, although it is significantly noisier than σ k yy (y). We believe that the noise is associated with the layering of the particles near the side walls as we discuss in section 4. We superpose the plot of σ c yy (y) with a line representing 7.29 − (A exp(B|y|)) where 7.29 is the value of σ c yy (y) in the center of the cell and A and B are the fitting parameters for σ k yy (y) = A exp(B|y|) given in the caption of figure 6. Considering the noise in this plot, we find reasonable agreement. We are currently considering a method to relate the variations in f to those in σ c yy (y) as we discuss in section 4. Nevertheless, for the theory developed in this paper we consider this to be evidence that the data support the theoretical prediction from equation (5).

Partial stresses
In figure 7, we plot the contact and kinetic stresses associated with the two species. In figure 7(a), we plot profiles of the contact stresses associated with collisions between two small particles (σ c,ss ), two large particles (σ c,ll ) and one large and one small particle (σ c,sl ). In figures 7(b) and (c), we plot profiles of the partial kinetic and contact stresses, respectively, for the two constituents. In figure 7(d), we plot the concentration profile of each constituent. These results indicate that, depending on the time and the region of the chute, either the small or the large particles may take up a higher fraction of the local stress. They also show that in locations in the chute where there is an equal percentage of both species, the partial stress of the smaller particles is slightly higher than that of the large particles. Now, we investigate how the contact and kinetic stresses partition for different sized particles compared with their representative volume fraction in the mixture. In figure 8, we plot R c,i = ψ c,i /φ i and R k,i = ψ k,i /φ i calculated for three different time intervals. These results are rather noisy, but they indicate that at nearly all times and all degrees of segregation, the small particles carry a higher fraction of the normal stresses than their concentration; the large particles carry a lower fraction of the normal stresses than their concentration. The results for the kinetic normal stresses are consistent with the results reported by Hill and Zhang [20], who showed that in dense sheared mixtures, small particles should have larger velocity fluctuations than large particles. By definition, this indicates that the small particles will carry more kinetic stress than the large particles. Our results for the contact stresses are inconsistent with the assumption presented by Gray and colleagues [10,11] that ψ i,L /φ i,L > 0 and ψ i,S /φ i,S < 0 for the gravity-induced segregation they studied. This discrepancy may be due to the different nature of shear-induced and gravity-induced segregation, an issue we are currently investigating. Now we consider whether or not, in light of these data, equation (11) predicts segregation trends consistent with the simulation results. In figure 9(a), we plot the profiles of R c,s − R k,s in the y-direction at three different time intervals: (i) near the beginning of the simulation when segregation is still active, (ii) near the end of the simulation when segregation is nearly complete and (iii) halfway between. We can see that for all cases R c,s − R k,s < 0 except for the central region where segregation is minimal. Here, R c,s − R k,s ≈ 0. Considering this and the observation that ∂σ k yy /∂ y > 0 for y < 0 and ∂σ k yy /∂ y < 0 for y > 0 (see figure 6(a)), equation (11) predicts that small particles segregate to the center. This is supported by the analogous consideration of R c,l − R k,l , primarily positive as indicated in figure 9(b). Considering this in equation (11) with the spatial dependence of ∂σ k yy /∂ y indicates that the large particles will segregate to regions close to side walls. These predictions agree with our simulation results as shown in figure 2. In figures 9(c), we plot the concentration of the small particles at the same three time intervals. Comparing these with figure 9(b) it appears that the magnitude of R c,l − R k,l increases with φ s . The dependence of R c,l − R k,l on φ s is plotted in figures 9(d) excluding the regions in the center and small regions close to boundaries. At early times, there is roughly a linear relation between R c,l − R k,l and φ s , although this is no longer true at later times. To develop a predictive form of the model, in the next section we approximate the relationship between R c,i − R k,i and φ j to be linear based on the early time results.

Theoretical prediction compared with simulation results
We now further evaluate the model outlined in section 2 by determining how well a prediction of the evolution of local species concentration φ i for species i matches the simulation results. We first consider the equation of conservation of mass for species i. With no gradients in the x and z directions, for constant material density ρ m , we can rewrite equation (6) as As we mentioned in the discussion above, we make the approximation that R c,l − R k,l = Bφ s , where B is a positive constant, based on early time results shown in figure 9(d). Then equation (11) may be rewritten for large particles as Upon substituting equation (30) into (29), we find that This continues to follow the formalism first presented by Gray and colleagues (e.g. [10,11]). They essentially suggested a form for the partitioning of the pressure between species that gave rise to the linear relationship R p,i − 1 = B φ j , where B is an undetermined parameter. Using this with equation (15) and the equation for the conservation of mass for species i with a spatial gradient only in thez-direction, 20 an equation analogous to equation (31) for the system studied by Gray and colleagues (e.g. [11]) could be written: The temporal-spatial profiles of the concentration of large particles can then be obtained by solving equation (31) numerically provided ρ(y), ∂σ k yy /∂ y, diffusivity D = d/c D and the segregation coefficient q = B/c D can be obtained. For this, we use the profiles of the mixture solid fraction ( figure 1(b)) and normal kinetic stresses in the y direction ( figure 6(a)) obtained from the simulations. We do not know D and q, so we choose constant values for these two parameters empirically by comparing the predictions obtained using different values of D and q to the steady-state simulation results.
For our numerical solution, we use the initial conditions consistent with a homogeneous mixture: φ s = φ l = 0.5 at t = 0 for all values of y (34) and no-flux conditions at the two walls We then discretize the problem and solve numerically by using a central difference scheme for spatial derivatives and the fourth-order Runge-Kutta method for time integration. Figures 10(a) and (b) show temporal-spatial profiles of the concentration of large particles from simulation data and theoretical predictions up to 1000 s, respectively. Based on trial and error we chose q = 2.5 × 10 −3 s and D = 0.05 mm 2 s −1 . In both the simulation results and theoretical predictions, the large particles segregate to the side walls and small particles segregate toward the center. In the middle slow creeping region, where the gradient of normal kinetic stresses is very small, the segregation process is much slower than in other regions. All of these indicate good qualitative agreement between theoretical predictions and simulation results, although they do not reach perfect quantitative agreement.
There are a number of steps that can be taken to improve the predictive power of the model. First, a more sophisticated relationship between R c − R k and the flow properties is needed for this model, requiring simulations or experiments with higher resolution in time and/or space, currently under way. The functional forms of the coefficients q and D should also be considered, including whether or not their values are dependent on local flow properties such as local shear rate. As mentioned by Wiederseiner et al [57], for flow down an inclined plane, D may weakly depend on shear rate. Natarajan et al [58] suggest that the diffusion coefficient D has a significant dependence on both the fluctuating velocities and shear rate in vertical chute flow.

Effects of granularity
We now briefly consider the effects of the granularity of the system, particularly the layering that occurs next to the walls.
To investigate this, we first plot some of the profiles of the kinematics and the stresses in figure 11. These are comparable to some of the plots in figure 1 and 6, although calculated at a much higher resolution (i.e. smaller bin size y). Regular oscillations in f (y) (figure 11(a)) suggest an ordering adjacent to the wall. The wavelength of the fluctuations is approximately 3 mm, the average diameter of the large particles, almost as if layers of beads organized in planes parallel to the walls were sliding vertically relative to one another. Similar phenomena have been previously reported for free surface sheared flows in drums (e.g. [59]) and are evidently apparent for f (y), dw/dy as well as the contact and kinetic stresses in chute flow. Since the manner in which the contact stress is calculated focuses the computed stress associated with each interparticle contact entirely at the point of contact, when particles are densely arranged in chains such as indicated in figure 11(a), σ c will be high where f is low and low where f is high. This does not, of course, account for the stress taken up within each solid particle, but only that associated directly with particle-particle contacts.
To obtain a measure of the fluctuations in f for comparison with fluctuations in σ c , we first perform two sets of fits for f (y). To simplify the discussion, we fit according to distance from one of the side wallsỹ = 23 mm − y. The first fit represents the overall apparent exponential decay in f (y) adjacent to the wall: f exp = f max − f o exp(−ỹ/y e ). The second fit represents the oscillations in f whose amplitude decays exponentially away from the wall: f osc = A exp (−ỹ/y o ) sin(2πỹ/λ + ϕ). The parameters are given in the caption of figure 12. In figure 12(a), we plot the result from the fit ( f exp + f osc versus y) on top of the simulation data, showing that the form of the fit represents the data fairly well. Then we consider the form of the fit from the kinetic stress that we used for the contact stress in   The oscillations are somewhat larger in σ c yy (y) than for f -by trial and error, we find the oscillations in σ c yy (y) to be approximately five times those in f . To test our hypothesis that oscillations in f correspond to oscillations in σ c yy (y), in figure 12(b) we plot σ c yy (y) exp + 5 f osc on top of the simulation data for σ c yy (y) and once again find reasonable agreement.
We conclude that the oscillations in σ c yy are directly associated with oscillations in f (y). This is an effect not captured in the continuum approach described in this paper, and its role in segregation is a topic currently under investigation.

Summary and discussion
In this paper, we have developed a hydrodynamic model for shear-induced segregation of binary mixtures of different sized particles in dense sheared flow. By considering two different types of stresses, the kinetic and contact stresses, and allowing the partitioning of stresses between constituents to be different from their concentration in the mixture, we show that segregation can be driven solely by the dynamics associated with a shear gradient. Specifically, we show how gradients in kinetic stresses that arise from a shear rate gradient in combination with a kinetic sieving mechanism will segregate large particles to regions of high velocity fluctuation correlations and small particles away from them.
Although the framework is reasonable for shear-induced segregation, and predictions appear to correlate reasonably well with observations, for a predictive model a deeper understanding of the rheology of dense sheared mixtures is needed. First, we need a relationship between R c,i − R k,i and flow properties such as particle concentrations and flow velocities to close the governing equations. For this, we have temporarily used a linear relationship between R c,i − R k,i and φ j (for species i and j, however, it is clear from the simulation results that the relationship changes as the segregation in the system evolves, so this needs to be improved). Additionally, it is not clear whether the drag force on the particles should be linear with relative velocity or, even if it is linear, what the coefficient of drag should be. Further, although we set the diffusivity D as a constant, as mentioned in [57] and [58] and demonstrated in [59], D may depend on the local shear rate as well as other factors. A more mechanistic way to obtain relationships for D and c D as they depend on the kinematics of the flow is necessary of a predictive model of shear-induced segregation.
We note that our model for shear-induced segregation is complementary to that of gravitydriven segregation developed by Gray and colleagues [10,11]. The models represent different phenomenologies that coexist in many systems, particularly those of free-surface gravity-driven flow commonly used to investigate segregation effects in granular mixtures. Therefore, it would be useful to consider our model of shear-induced segregation combined with the models for gravity-driven segregation [10,11] for a more complete model for predicting segregation in dense, gravity-driven flow in which shear rates are present.
Finally, it is interesting to consider these results in the context of the results from other systems such as shaken systems that were studied experimentally by Knight et al [40] and theoretically by Shinbrot et al [41]. In their studies, a container holding particles was tapped; the particles in the center moved upward with each tap, and the particles near the edges moved downward. The average flow profile derived from this movement was well fitted by a hyperbolic cosine similar to the profile we found in our chute flow ( figure 1(c)). Void transport has been used successfully as a basis for modeling other granular flows, such as granular flow past an obstacle by Caram and Hong [60]. Indeed, the physical motivation behind the framework we propose for segregation involves the statistical distribution of hole sizes (as in [9]). Considering that the stochastic model of hole transport proposed by Shinbrot et al [41] produces a velocity profile similar to ours, it could be that a similar model allowing for different hole sizes is capable of producing the segregation flux profiles we observe as well. We note that ρv v scales roughly with T (figures 1(e) and (f)), so that equation (A.5) indicates that a gradient in T can be associated with a gradient in both σ k yy and σ c yy (of opposite signs). In figure A.1, we plot the profiles of σ c yy and σ k yy with the profiles of the relevant neglected terms (ρ v v), (v ρ v ) and (ρ i v i v i ). In doing so, we find the neglected components to be several orders of magnitude below σ c yy and σ k yy and consider it reasonable that we neglected related terms here and in section 2.
Next, we investigate the conservation laws for two species: Here, all of the terms used for species i are identical in form to those used for the mixture in equation (A.2), with the exception of the new term β i , which represents the interaction force exerted on species i by the other species. We note that σ i is the local stress borne by species i and that the total stress σ = σ i . By performing Reynolds decomposition and then averaging both sides of equation (A.8) we may use similar arguments to those we employed for the mixture and obtain a somewhat simplified form of the conservation equation for each species. Again, we assume terms involving correlations between velocity and density fluctuations such as ∂ ∂ x k (u i j ρ i u i k ), ∂ ∂ x k (u i k ρ i u i j ) and ∂ ∂ x k (ρ i u i k u i j ) are negligible. We verify that this is reasonable for the relevant terms in figure A.1. Then, for the y-direction we obtain The first three terms in equation (A.9) are not so easily discounted for the individual species as are analogous terms for the mixture case, as we expect v i to be nonzero and changing during segregation, and the spatial derivative of ρ i to be nonzero after segregation. However, in dense sheared flows v i has been shown to be typically quite small except when ρ is uniform (e.g. [25] and [37]) and as may be seen in figure 2(d) ∂v i /∂ y is still quite small even during segregation. We therefore neglect these acceleration terms as well. With these assumptions we can rewrite equation We note that all terms in equations (A.6) and (A.10) are averaged and drop the overbar from this point forward, so that unless noted for each variable q alone refers to the average q. Similar to the model for partial pressures in a gravity-driven system proposed by Gray and Thornton [10], we propose that where ψ c,i and ψ k,i determine the proportion of normal contact and kinetic stresses carried by species i and are not necessarily equal to φ i . To ensure that i σ c,i yy = σ c yy and i σ k,i yy = σ k yy , two constraints on ψ c,i and ψ k,i must be satisfied: (i) that ψ l + ψ s = 1, and (ii) if only one species is present then it must support the entire local stress, so that when φ i = 1, ψ j = δ i j , where δ i j is the Kronecker delta function.
Considering equation (A.11), we rewrite equation (A.10) as For the interaction term β i y , we propose a similar form to that proposed by Gray and Chugunov [11]: following [10,11]. The first two terms on the right-hand side of the equation ensure that, as in Darcy's law, the percolation process (in this case, the kinetic sieving process) is driven by intrinsic rather than partial stress gradients. The third term is a linear drag law similar to that provided by Morland [49] for the percolation of fluids-c D is the linear drag coefficient. The fourth term acts as a 'remixing force' (e.g. [11] and [49]) that drives grains of constituent i towards areas of lower concentration-d is the ordinary diffusion coefficient.
We substitute equation (A.13) into (A.12) and obtain 0 = −ψ c,i ∂σ c (A.14) Combining this equation with equation (A.6), we find that we can express a segregation flux of species i as To facilitate a physical interpretation of equation (A.15), we introduce partitioning variables R c,i = ψ c,i /φ i and R k,i = ψ k,i /φ i , to represent the ratio between the fraction of the local stress (both contact and kinetic) carried by each species and its local concentration in the mixture. We then rewrite equation (A.15) as (A. 16)