1 Introduction
Granular materials are inherently difficult to mix. When agitated, it is common for particles of differing size, density or shape to separate autonomously (Brown Reference Brown1939). This effect is termed segregation, and is a common feature of many natural and industrial granular flows such as landslides (Zhang et al. Reference Zhang, Xu, Huang and Chang2011), debris flows (Iverson Reference Iverson2003), snow avalanches (Bartelt & McArdell Reference Bartelt and McArdell2009) and industrial operations (Williams Reference Williams1976). When granular material is tapped or vibrated, this same observation is often referred to as the ‘Brazil nut effect’, where large particles (Brazil nuts being the largest nut in a bag of mixed nuts) rise to the top of a system of smaller particles (Möbius et al. Reference Möbius, Lauderdale, Nagel and Jaeger2001). A proposed mechanism responsible for this segregation, for which there are many competing propositions (Knight, Jaeger & Nagel Reference Knight, Jaeger and Nagel1993; Hong, Quinn & Luding Reference Hong, Quinn and Luding2001; Huerta & Ruiz-Suárez Reference Huerta and Ruiz-Suárez2004), is the formation of new voids during movement of the grains (Savage & Lun Reference Savage and Lun1988). If the system is being compressed, particles will attempt to move into these new void spaces. As it is easier for a small particle to fill a void than a larger one, there is a preference for small particles to accumulate in areas of relatively high void formation (Fan & Hill Reference Fan and Hill2011).
The Brazil nut effect, originally explained in terms of the local rearrangement of particles, was subsequently found to also be induced by the convection of the bulk particles (Knight et al. Reference Knight, Jaeger and Nagel1993). In these cases, bulk convection eventually carried all of the particles towards the top of the system, but large intruder particles were somehow trapped there. In other geometries, a coupling of advection and segregation has been observed (Khakhar, McCarthy & Ottino Reference Khakhar, McCarthy and Ottino1997; Hill et al. Reference Hill, Gilchrist, Ottino, Khakhar and McCarthy1999a ,Reference Hill, Khakhar, Gilchrist, McCarthy and Ottino b ). In the context of gravity currents, we often observe a collection of the largest particles at the snout of the flow, i.e. at the furthest point downslope, with medium sized particles at the free surface, and fines accumulated at the base of the flow (Gray & Ancey Reference Gray and Ancey2009). Laterally, large grains are deposited to form levees, due to an interaction between vertical segregation and lateral advection, which strongly affect the dynamics of the flow (Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012; Woodhouse et al. Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012; Baker, Johnson & Gray Reference Baker, Johnson and Gray2016). We know therefore that complex segregation patterns generally evolve in all of the available spatial dimensions of a given flow, and that we cannot simply decouple the segregation dynamics from the bulk flow. Direct measurement of the motion of each grain during oscillatory shear of a bidisperse mixture using refractive index matched scanning has shown that individual particles segregate at a velocity that depends on their size, with small particles moving in large steps and larger particles at a more constant rate (van der Vaart et al. Reference van der Vaart, Gajjar, Epely-Chauvin, Andreini, Gray and Ancey2015).
The first theoretical models of segregation by size in granular flows described the system using statistical mechanics or by applying kinetic theory to sparse granular flows (Jenkins & Mancini Reference Jenkins and Mancini1987; Savage & Lun Reference Savage and Lun1988; Jenkins & Mancini Reference Jenkins and Mancini1989). More recently, the kinetic theory description of segregation has been extended to dense flows (Jenkins & Yoon Reference Jenkins and Yoon2002). Other models have diverged in methodology, describing the segregation phenomenon using continuum mechanics (Dolgunin & Ukolov Reference Dolgunin and Ukolov1995; Gray & Thornton Reference Gray and Thornton2005; Gray & Ancey Reference Gray and Ancey2011; Marks, Rognon & Einav Reference Marks, Rognon and Einav2012; Hill & Tan Reference Hill and Tan2014; Tunuguntla, Bokhove & Thornton Reference Tunuguntla, Bokhove and Thornton2014). These continuum models in general have assumed that the direction of segregation is known a priori, where large particles go ‘up’ and small particles go ‘down’ (or vice versa). Efforts have been made to describe the direction of segregation in more complex systems, but only insofar as segregation is allowed in a single direction (Gray & Kokelaar Reference Gray and Kokelaar2010; Fan & Hill Reference Fan and Hill2011; Hill & Tan Reference Hill and Tan2014). Additionally, the coupling of orthogonal advection and segregation has been studied (Thornton & Gray Reference Thornton and Gray2008; Gajjar et al. Reference Gajjar, van der Vaart, Thornton, Johnson, Ancey and Gray2016).
A typical assumption in such models is that the medium is a mixture, where each constituent phase of the mixture is composed of particles with a uniform size. In Gray & Thornton (Reference Gray and Thornton2005), a keen insight was made that if the stress of the mixture was not carried proportionately by each phase, this stress imbalance would cause segregation. Recently, measurement of the micro-mechanical behaviour of these systems has shown that there are two subcategories of this stress, the contact stress and the kinetic stress, and that they are distributed between grains differently. The intrinsic contact stress, which comes from sustained contacts between particles, has been shown to be equal among the phases, and so cannot be driving segregation (Fan & Hill Reference Fan and Hill2011; Voivret Reference Voivret2013; Weinhart et al. Reference Weinhart, Luding, Thornton, Yu, Dong, Yang and Luding2013; Staron & Phillips Reference Staron and Phillips2015). The other component of stress in granular systems, which comes from the collision of particles, is known as the kinetic stress in direct analogy to turbulence in fluid flows. The same micro-mechanical investigations have found that this component of the stress does indeed vary with particle size, although the magnitude of this stress is in general much smaller than the contact stress. This implies that while stationary, the system is at equilibrium and segregation will not occur, but during motion, when the kinetic stress is non-zero, segregation is promoted. The specification of kinetic stress-induced segregation was first included in a continuum model in Fan & Hill (Reference Fan and Hill2011), and later modified in Weinhart et al. (Reference Weinhart, Luding, Thornton, Yu, Dong, Yang and Luding2013), Hill & Tan (Reference Hill and Tan2014), where the effects of stress partitioning were examined in detail.
Here, we show using a continuum model that if the relative magnitude of the kinetic stress is a function of the grainsize, then we can describe both the time-averaged and the fluctuating components of the segregation velocity purely using kinematic variables – namely the velocity and kinetic stress of an equivalent monodisperse system. We go on to present a simple two-dimensional experiment where such kinematic observations can be made, and then compare predictions and observations of the segregation behaviour of individual intruder particles in two spatial dimensions.
In the following, the kinetic stress arguments presented in Hill & Tan (Reference Hill and Tan2014) will be extended to a polydisperse (many sized) material. This serves several important purposes. Firstly, by including the physical size of the constituent particles in the analytic description, we can create a framework that genuinely predicts the behaviour of arbitrary mixtures of particles, rather than having to measure a new fitting parameter for each pair of constituents. Secondly, this method removes the binning effects associated with the assumption that each constituent is represented by a single particle size, rather than a narrow distribution, as shown in Marks et al. (Reference Marks, Rognon and Einav2012). Finally, a polydisperse material can be extended to additionally treat other mechanisms, such as grain crushing, as shown in Marks & Einav (Reference Marks and Einav2015), and agglomeration and melting.
2 Continuum model
As has been shown previously (Marks et al.
Reference Marks, Rognon and Einav2012), the analytic description of a mixture of arbitrarily sized particles can be concisely stated in the context of population balance equations (Ramkrishna Reference Ramkrishna2000). Such a theoretical description is constructed from a five-dimensional space, comprised of external space
$\boldsymbol{x}=\{x,y,z\}$
, time
$t$
and an internal coordinate
$s$
, which describes the grainsize of the particles. A polydisperse granular material can then be described as a single material, with an additional property
$\unicode[STIX]{x1D719}(\boldsymbol{x},\boldsymbol{s},\boldsymbol{t})$
that is a probability density function which characterises the grainsize distribution at any point in space. For example, if the material segregates perfectly by size, the grainsize distribution will approach a delta function. As in mixture theory (Morland Reference Morland1992), we maintain a distinction between the partial value of a property, which is the amount of that property within a representative volume element, and the intrinsic value, which is the density of that property per unit concentration.
For the polydisperse granular material described here, the partial density of the material,
$\unicode[STIX]{x1D70C}(\boldsymbol{x},s,t)$
, is defined as the mass of particles in the grainsize interval
$[s,s+\text{d}s]$
per unit volume of the mixture. The volume fraction of each grainsize is defined by the probability density function
$\unicode[STIX]{x1D719}=\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}^{\ast }$
, where
$\int \unicode[STIX]{x1D719}~\text{d}s=1$
and
$\unicode[STIX]{x1D70C}^{\ast }$
is the intrinsic density, defined as
$\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}^{M}\unicode[STIX]{x1D702}$
, where
$\unicode[STIX]{x1D70C}^{M}$
is the material density and
$\unicode[STIX]{x1D702}$
is the total solid fraction of all grains. As we will here only consider the case of a material with varying size, but equal in all other properties, we set
$\unicode[STIX]{x1D70C}^{M}$
to be independent of
$s$
. Finally, the bulk density is defined as
$\bar{\unicode[STIX]{x1D70C}}=\int \unicode[STIX]{x1D70C}~\text{d}s=\unicode[STIX]{x1D70C}^{\ast }$
. Similarly, partial, intrinsic and bulk contact stress are defined by
$\unicode[STIX]{x1D748}=\unicode[STIX]{x1D719}\unicode[STIX]{x1D748}^{\ast }$
and
$\bar{\unicode[STIX]{x1D748}}=\int \unicode[STIX]{x1D748}~\text{d}s$
, respectively. Using this notation, we can formulate conservation of momentum for this system, following arguments in Bedford & Drumheller (Reference Bedford and Drumheller1983), and further detailed in appendix A, as
$$\begin{eqnarray}\frac{\text{D}(\unicode[STIX]{x1D70C}\boldsymbol{u})}{\text{D}t}=\unicode[STIX]{x1D70C}\boldsymbol{g}-\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{\ast }+\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D737},\end{eqnarray}$$
where
$\boldsymbol{u}(\boldsymbol{x},s,t)=\{u_{x},u_{y},u_{z}\}$
is the velocity field,
$\boldsymbol{g}$
the acceleration due to gravity,
$\text{D}/\text{D}t$
the material derivative,
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }$
the divergence operator in external space
$\boldsymbol{x}$
and
$\unicode[STIX]{x1D737}(\boldsymbol{x},s,t)$
is an interaction term which describes how momentum is exchanged between different grainsizes, where by definition
$\bar{\unicode[STIX]{x1D737}}=\int \unicode[STIX]{x1D737}~\text{d}s=\mathbf{0}$
, such that there is no net momentum flux into or out of the flowing material due to particle interactions. Following Hill & Tan (Reference Hill and Tan2014), we can then decompose the velocity field into a mean and fluctuating part as
$\boldsymbol{u}=\langle \boldsymbol{u}\rangle +\boldsymbol{u}^{\prime }$
, respectively, where
$\langle \boldsymbol{u}\rangle =1/T\int _{t-T/2}^{t+T/2}\boldsymbol{u}\,\text{d}t$
, for some appropriate time window
$T$
, as
$$\begin{eqnarray}\frac{\text{D}}{\text{D}t}\left(\unicode[STIX]{x1D70C}\langle \boldsymbol{u}\rangle +\unicode[STIX]{x1D70C}\boldsymbol{u}^{\prime }\right)=\unicode[STIX]{x1D70C}\boldsymbol{g}-\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{\ast }+\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D737}.\end{eqnarray}$$
We now take a time average of (2.2), again over the same time window, and allow for the fact that in general time averages of products of fluctuating quantities do not vanish (and dropping the
$\langle \rangle$
notation for terms other than
$\boldsymbol{u}$
). We additionally assume that the flow is incompressible (isochoric), which after some rearrangement yields
$$\begin{eqnarray}\frac{\text{D}(\unicode[STIX]{x1D70C}\langle \boldsymbol{u}\rangle )}{\text{D}t}=\unicode[STIX]{x1D70C}\boldsymbol{g}-\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{\ast }-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{k}+\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D737},\end{eqnarray}$$
where
$\otimes$
is the dyadic product,
$\unicode[STIX]{x1D748}_{k}=\unicode[STIX]{x1D70C}\langle \boldsymbol{u}^{\prime }\otimes \boldsymbol{u}^{\prime }\rangle$
is the partial kinetic stress, which represents the additional stress induced by particle collisions and we have assumed that for all variables other than the velocity field, their fluctuations are uncorrelated. We note at this stage some recent measurements of the stress scaling between different sizes of grains in a polydisperse mixture. Several authors have found that the intrinsic contact stress in both a stationary (Voivret Reference Voivret2013) and flowing (Hill & Tan Reference Hill and Tan2014; Staron & Phillips Reference Staron and Phillips2015) medium is independent of grainsize, such that
$\unicode[STIX]{x1D748}^{\ast }(\boldsymbol{x},s,t)\equiv \bar{\unicode[STIX]{x1D748}}(\boldsymbol{x},t)$
. Conversely, it appears that the kinetic stress is a function of the grainsize, such that there is some scaling law, which is as yet unknown, which characterises how these stresses are carried (Hill & Tan Reference Hill and Tan2014). In direct analogy to Marks et al. (Reference Marks, Rognon and Einav2012), we here close the system of equations by assuming that the intrinsic kinetic stress scales with the bulk kinetic stress as
$\unicode[STIX]{x1D748}_{k}^{\ast }=f(\boldsymbol{x},s,t)\bar{\unicode[STIX]{x1D748}}_{k}$
, where
$f$
defines the scaling law for kinetic stress, and we require that
$\int \unicode[STIX]{x1D719}f~\text{d}s=1$
to conserve bulk momentum. Using these definitions, we integrate (2.3) over
$s$
, retrieving the bulk behaviour of the mixture, and find that
$$\begin{eqnarray}\frac{\text{D}(\bar{\unicode[STIX]{x1D70C}}\bar{\boldsymbol{u}})}{\text{D}t}=\bar{\unicode[STIX]{x1D70C}}\boldsymbol{g}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D748}}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D748}}_{k},\end{eqnarray}$$
where
$\bar{\boldsymbol{u}}=\int \unicode[STIX]{x1D719}\langle \boldsymbol{u}\rangle ~\text{d}s$
. This recovers the usual statement of conservation of bulk momentum, such as shown in Hill & Tan (Reference Hill and Tan2014). We further define the time-averaged segregation velocity as
$\hat{\boldsymbol{u}}=\langle \boldsymbol{u}\rangle -\bar{\boldsymbol{u}}$
, and using this definition we substitute (2.4) into (2.3), while also setting
$\unicode[STIX]{x1D735}(\unicode[STIX]{x1D719}f)\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D748}}_{k}=0$
, following the same logic as in appendix A leading to (2.2) where a similar term created spurious additional diffusion, giving
$$\begin{eqnarray}\frac{\text{D}(\unicode[STIX]{x1D70C}\hat{\boldsymbol{u}})}{\text{D}t}+\bar{\unicode[STIX]{x1D70C}}\bar{\boldsymbol{u}}\frac{\text{D}\unicode[STIX]{x1D719}}{\text{D}t}=\unicode[STIX]{x1D719}(1-f)\unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D748}}_{k}+\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D737}.\end{eqnarray}$$
A further constitutive assumption is required to close the system, which is the selection of an appropriate interaction term,
$\unicode[STIX]{x1D737}$
. As in Gray & Thornton (Reference Gray and Thornton2005), we assume a linear drag between different species, such that species which flow faster or slower than the bulk velocity will feel a drag force, so that
$$\begin{eqnarray}\unicode[STIX]{x1D737}=-\unicode[STIX]{x1D719}c\hat{\boldsymbol{u}},\end{eqnarray}$$
where
$c$
is a parameter which controls the magnitude of the drag force, with units of inverse time. The direct measurement of
$c$
has not been attempted here, but could be done using the method proposed in Guillard, Forterre & Pouliquen (Reference Guillard, Forterre and Pouliquen2016). At steady state (when
$\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}t=0$
and
$\unicode[STIX]{x2202}\hat{\boldsymbol{u}}/\unicode[STIX]{x2202}t=0$
), under the assumption that the time-averaged segregation velocity is small compared to the bulk velocity (
$\hat{\boldsymbol{u}}\ll \bar{\boldsymbol{u}}$
), it can be shown from conservation of mass that both terms on the left-hand side of (2.5) vanish, as shown in appendix B. Under these conditions,
$\hat{\boldsymbol{u}}$
reduces to
$$\begin{eqnarray}\hat{\boldsymbol{u}}=\frac{1-f}{\bar{\unicode[STIX]{x1D70C}}c}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D748}}_{k}.\end{eqnarray}$$
This relation allows for prediction of the segregation velocity for particles of any size if the bulk kinetic stress field is known. This further implies that any particles which are subject to less than the average kinetic stress, i.e.
$f<1$
, will segregate in the direction of the bulk kinetic stress gradient, while those that fluctuate more than the average will segregate against it.
Qualitatively, equation (2.7) predicts the segregation observed during inclined chute flow, where there exists a kinetic stress gradient perpendicular to the slope, and we observe large particles accumulating at the free surface. Additionally, as shown in Hill & Tan (Reference Hill and Tan2014), this term predicts well the behaviour in vertical chute flow.
In the following, we describe an experimental apparatus wherein we measure time-averaged and fluctuation velocities of a closed system of particles subject to continuous shear in a perpetual avalanche geometry. We firstly measure the velocity fields of a monodisperse set of particles, and then using the relation (2.7), predict the velocity fields of arbitrarily sized intruder particles, and compare these with experimental measurements of the behaviour of intruder particles. Because complex feedback mechanisms exist between the grainsize distribution and the bulk kinetic stress (Jenkins & Mancini Reference Jenkins and Mancini1987), we probe the system at
$\lim _{\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D6FF}(s-\bar{s})}$
, i.e. the limit where a single intruder particle is introduced to a field of monodisperse particles of grainsize
$\bar{s}$
, where
$\bar{s}=\int \unicode[STIX]{x1D719}s~\text{d}s$
.
3 Experiment
The experimental apparatus consists of two glass plates, separated
$3.20\pm 0.10~\text{mm}$
apart, sitting
$0.2~\text{mm}$
above a GT2 timing belt,
$6~\text{mm}$
wide, with teeth
$2~\text{mm}$
deep, as shown in figure 1. A single layer of nylon particles, with bulk density
$0.92~\text{g}~\text{cm}^{-3}$
, laser cut into regular pentagonal prisms, with height
$3.0\pm 0.1~\text{mm}$
, is placed within the void between the two glass plates, pentagonal sides facing the glass. Pentagons are used so that a very narrow range of sizes (measured by micrometer to have side lengths
$3.90\pm 0.10~\text{mm}$
) does not cause crystallisation. Rigid walls are placed at each end of the cell, perpendicular to the timing belt, spaced
$L=496.0\pm 1.0~\text{mm}$
apart, such that the particles cannot leave the cell. A DC motor is used to turn the timing belt, and the subsequent motion of the belt causes the particles to be sheared. At low shear rates, intermittent avalanches are observed at the free surface. With increasing belt velocity, there is a transition from intermittent to continual avalanching. Further increase of the belt velocity creates a gaseous state, where particles are rarely in contact, near both the timing belt and the free surface. A similar experimental geometry was used in Perng, Capart & Chou (Reference Perng, Capart and Chou2006), where the kinematics of the flow in such a geometry were studied in detail. This is a two-dimensional analogue of the experimental apparatus described in Gajjar et al. (Reference Gajjar, van der Vaart, Thornton, Johnson, Ancey and Gray2016), where index matching was used to observe similar breaking size-segregation waves in three dimensions.
Figure 1. Schematic of the front and side views of the perpetual avalanche experimental apparatus. A single layer of pentagonal prism-shaped particles are held between two glass plates and sheared from below by a conveyor belt. The system is tilted with respect to gravity such that the base of the system is parallel to the free surface.
Here we report on experiments at a single belt velocity of
$u_{b}=0.190\pm 0.010~\text{m}~\text{s}^{-1}$
, slope angle of
$29.8\pm 0.2^{\circ }$
(as close as possible to the dynamic angle of repose at this velocity, adjusting the slope manually until this condition was reached) and filling depth of
$H\approx 100~\text{mm}$
, or approximately 30 pentagonal circumradii. At this state the inertial number is
$I=2u_{b}\bar{s}/(H\sqrt{gH})\approx 0.016$
(assuming no slip between the conveyor belt and the flow, as observed experimentally), in the dense fluid regime, and we observe a continuous avalanche with minimal saltation of particles at the free surface and relatively uniform packing fraction throughout the flow, such that the assumption of isochoric flow is reasonable. We record with a Photron SA5 high-speed camera at 1000 frames per second, and a resolution of
$1024\times 256$
pixels, and use PIVlab (Thielicke & Stamhuis Reference Thielicke and Stamhuis2014) to measure coarse-grained velocity fields from 210 000 pairs of recorded images, with pairwise spacing of 0.001 s, covering a time span of 3.5 min. These measured velocity fields have a spatial resolution of 1.9 mm. By averaging in time, we are able to retrieve both
$\langle \boldsymbol{u}\rangle$
and
$\boldsymbol{u}^{\prime }$
. Due to the steady nature of the flow in this geometry, we are able to average over the full duration of recording to recover time-averaged properties. The relevant continuum fields are shown in figure 2, which summarises the kinematics of the system.
Figure 2. Steady state distributions of a flow of uniformly sized particles. (a) Time-averaged downslope velocity,
$\bar{u}_{x}$
$(\text{m}~\text{s}^{-1})$
. (b) Time-averaged cross-slope velocity,
$\bar{u}_{z}$
$(\text{m}~\text{s}^{-1})$
. (c) Magnitude of the shear strain rate,
$|\dot{\unicode[STIX]{x1D6FE}}|$
$(1~\text{s}^{-1})$
. (d) Kinetic pressure,
$\text{tr}(\bar{\unicode[STIX]{x1D748}}_{k})/2$
(
$\text{m}^{2}~\text{s}^{-2}$
). The black pentagon indicates the physical size of the particles in the experiment.
3.1 Segregation patterns
We now consider the case of single nylon pentagonal intruder particles flowing within a mass of monodisperse particles. By filming at 30 fps, for 100 000 frames (approximately 600 belt rotations) we collect a probability density map for the location of the centroid of a given intruder particle. The maps do not depend on the initial positions of the intruder particle, and have been averaged over several initial positions. Such maps are shown for a range of intruder sizes
$s$
, with
$s/\bar{s}=4.9$
,
$3.4$
,
$2.1$
,
$1.6$
and
$1.2$
, in figure 3. Intruder particles close to the size of the bulk particles do not segregate significantly, and as the intruder size increases, its location localises onto a fixed point on the right-hand side of the cell. Phenomenologically, we expect that a large particle will ‘rise’ when sheared, and so will flow towards the top right corner, where it will be subducted towards the left-going flow at the base. Subsequently, it will escape the rapid flow near the timing belt, doing so faster for larger intruders, forming a breaking size-segregation wave (Gajjar et al.
Reference Gajjar, van der Vaart, Thornton, Johnson, Ancey and Gray2016). Larger particles are then predicted to do smaller and smaller cycles at the right-hand side, which we observe in figure 3. The white lines on figure 3(a) depict streamlines of the time-averaged velocity of the intruder particle,
$\langle \boldsymbol{u}\rangle$
, as measured directly from particle tracking. For a video of this behaviour, please see the supplementary material available at https://doi.org/10.1017/jfm.2017.419.
We note that the large particles are attracted to a region of low shearing, whilst the smaller particles are repelled from it. Whether this behaviour is coincidental, or the signature of a competing mechanism for causing segregation, is at this stage unknown. Further work is required to investigate this possible alternative mechanism for segregation.
Figure 3. Locations and velocities of single intruder particles. (a) Normalised two-dimensional histogram of centroid of intruder particle. Top to bottom, the cases of an intruder particle with size
$s/\bar{s}=4.9$
,
$3.4$
,
$2.1$
,
$1.6$
and
$1.2$
, as shown to scale by the black pentagons. White lines are streamlines of corresponding time-averaged velocity of intruder particles as measured from experiment. (b) Normalised two-dimensional histogram of location of numerical tracer particles after 100 s of flow. White lines are streamlines of time-averaged velocity profile predicted using (2.7) from data recorded without an intruder particle, as shown in figure 2. For all cases
$c=22$
Hz.
4 Theoretical prediction
Given that we have measurements of the mean and fluctuating velocities, we can predict the velocity of an intruder particle using (2.7) as
$$\begin{eqnarray}\boldsymbol{u}=\bar{\boldsymbol{u}}+\frac{1-f}{\bar{\unicode[STIX]{x1D70C}}c}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D748}}_{k}+\boldsymbol{u}^{\prime },\end{eqnarray}$$
once we have defined
$f$
. Let us consider the collision of two grains
$a$
and
$b$
of radius
$1$
and
$R$
, respectively, and with the same density. Conservation of linear momentum in the centre-of-momentum reference frame implies that the velocities of each particle after impact are related by
$v_{a}=-R^{D}v_{b}$
, where
$D$
is the dimension. A small particle (
$R<1$
) will therefore have a higher fluctuating velocity than a larger one, and we can assume that
$f$
is a decreasing function of
$s$
, normalised such that
$\int \unicode[STIX]{x1D719}f~\text{d}s=1$
. This effect has been observed both numerically and experimentally in bidisperse systems varying in size, for both sparse and dense regions of flow (Hill & Zhang Reference Hill and Zhang2008; Staron & Phillips Reference Staron and Phillips2015). Under these conditions, we expect that the kinetic stress field scales such that small particles (
$f>1$
) advect against the kinetic stress gradient and large particles (
$f<1$
) advect with it. We therefore choose to define
$f$
as
$f=\bar{s}_{H}/s$
, where
$\bar{s}_{H}=1/\int (\unicode[STIX]{x1D719}/s)~\text{d}s$
is the harmonic mean grainsize. This simplistic assumption satisfies the two requirements outlined above. Additionally it qualitatively reproduces the asymmetry found from direct numerical measurement, see appendix C. We can then predict the flow field of an arbitrarily sized intruder particle by assigning the value of
$c$
, which here we take to be
$c=22$
Hz, which produces a reasonably good visual agreement with the experimental observations.
Figure 3(b) shows predictions for the velocity field and histogram of position of the intruder using (4.1) and the fields shown in figure 2. To make this prediction, we require additionally a value for the fluctuating component of the velocity,
$\boldsymbol{u}^{\prime }$
. Because of the definition of
$f$
, we can additionally relate
$\boldsymbol{u}^{\prime }$
to experimental measurements of the monodisperse particles, by randomly drawing from a Maxwell–Boltzmann distribution with root-mean-squared velocity equal to
$|\boldsymbol{u}^{\prime }|=\sqrt{f}|\bar{\boldsymbol{u}}^{\prime }|$
, being applied in a random direction at each time increment. All fields are described using third-order spatial interpolation and fourth-order Runge–Kutta temporal integration. Initially, 200 000 markers are spread throughout the system, and a histogram of their locations after 100 s of flow is shown on figure 3(b). It is evident that the continuum theory successfully predicts that segregation drives a single large intruder particle towards a fixed point attractor.
A further prediction of the velocity field for an intruder particle of size
$s=3\bar{s}$
is shown in figure 4. This velocity field has near zero divergence everywhere inside the domain of the experimental apparatus, and it is therefore interesting to analyse how particles are found to accumulate in a single region. For this reason, the finite-time Lyapunov exponent (calculated as
$\text{FTLE}=\ln (\sqrt{\unicode[STIX]{x1D706}})/|T|$
, where
$\unicode[STIX]{x1D706}$
is the maximum eigenvalue of the right Cauchy–Green deformation tensor and
$T$
is the time interval, see equation (12) from Shadden, Lekien & Marsden (Reference Shadden, Lekien and Marsden2005)) is shown in the inset of figure 4, computed using third-order spatial interpolation of the predicted time-averaged velocity field
$\langle \boldsymbol{u}\rangle$
and fourth-order Runge–Kutta temporal integration, integrated forwards in time over 45 s. Positive and negative finite-time Lyapunov exponents represent sources and sinks of particles, respectively. The negative value covering most of the domain indicates that large particles are predicted to converge towards the fixed point attractor at the centre of the inset. The larger the intruder particle, the more it is attracted towards this attractor, as shown in figure 3. Observations of a similar spiral pattern were made in three-dimensional flows in Johnson et al. (Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012). This behaviour may explain the rather perplexing results of Thomas (Reference Thomas2000), where systematic forensic excavation of segregating chute flows found moderately large particles (up to a size ratio of 3.5) at the surface, with even larger particles trapped below the surface, possibly due to a breaking size-segregation wave. This inference must be treated with caution, however, as the emplacement process of such particles is hard to discern by examining the deposit (Branney & Kokelaar Reference Branney and Kokelaar1992).
Figure 4. Predicting segregation using information from the monodisperse case shown in figure 2. Prediction of the intruder velocity field for a large particle with
$s=3\bar{s}$
, and the physical size indicated by the black pentagon. Inset: maximum finite-time Lyapunov exponent integrated forwards in time over 45 s.
5 Conclusions
We have here shown a new analytic description of grainsize segregation in flowing granular systems. Under the assumption that segregation is driven by kinetic stress gradients, predictions of segregation patterns were made using purely kinematic descriptors. We observed the formation of a fixed point attractor for large intruder particles in the theoretical predictions in the studied geometry. Additionally, experimental evidence of such segregation has been shown and compared with the theoretical predictions. The theory presented here is relevant to a large number of particulate flows, which for the first time is applicable to flows in which the direction of segregation is a function of the local flow conditions, and may follow complex paths in three spatial dimensions. This will be applicable to many natural and industrial flows, such as avalanche flow, levee formation and silo discharge. This new description of the segregation mechanism is also applicable to vibrated systems, such as the Brazil nut effect, where it is understood that there exist complex feedback mechanisms between segregation and bulk convection.
Acknowledgements
This work was partly supported by the Research Council of Norway through its Centres of Excellence funding scheme, project number 262644 and grant 213462/F20.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2017.419.
Appendix A
Following Ramkrishna (Reference Ramkrishna2000), we consider the particle state domain that contains the space domain
$\unicode[STIX]{x1D6EC}_{\boldsymbol{r}}$
whose elements are the three-dimensional vectors,
$\boldsymbol{x}=(x,y,z)$
, as well as the grainsize domain
$\unicode[STIX]{x1D6EC}_{s}$
whose elements are scalars
$s$
. Assuming that particles do not change size, mass and momentum conservation in the entire domain can be written as
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\int _{\unicode[STIX]{x1D6EC}_{s}(t)}\int _{\unicode[STIX]{x1D6EC}_{\boldsymbol{r}}(t)}\,\unicode[STIX]{x1D70C}\,\text{d}\boldsymbol{x}\,\text{d}s=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\int _{\unicode[STIX]{x1D6EC}_{s}(t)}\int _{\unicode[STIX]{x1D6EC}_{\boldsymbol{r}}(t)}\,\unicode[STIX]{x1D70C}\boldsymbol{u}\,\text{d}\boldsymbol{x}\,\text{d}s=\int _{\unicode[STIX]{x1D6EC}_{s}(t)}\int _{\unicode[STIX]{x1D6EC}_{\boldsymbol{r}}(t)}\,\boldsymbol{F}\,\text{d}\boldsymbol{x}\,\text{d}s, & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{F}(\boldsymbol{x})$
is the total force per unit volume, which will be discussed below. Using a generalisation of Reynold’s transport theorem to general vector spaces, and because the domain of these integrals is arbitrary and the integral is continuous, we can express the mass and momentum conservation in their local form as
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{u})=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}\boldsymbol{u})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{u}\otimes \boldsymbol{u})=\boldsymbol{F}, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D735}=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x+\unicode[STIX]{x2202}/\unicode[STIX]{x2202}y+\unicode[STIX]{x2202}/\unicode[STIX]{x2202}z$
is the differential operator and
$\otimes$
is the outer product. These two equations constitute a general framework that can describe the segregation dynamics of polydisperse granular material with a continuous grainsize distribution. Typically, three forces are considered to act on the granular assembly – a gravitational body force, a stress gradient and an interaction term, as
$$\begin{eqnarray}\boldsymbol{F}=\unicode[STIX]{x1D70C}\boldsymbol{g}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}+\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D737}.\end{eqnarray}$$
The nature of the stress gradient used above deserves further investigation. We decompose this partial stress gradient as
$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}=\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{\ast }+(\unicode[STIX]{x1D735}\unicode[STIX]{x1D719})\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{\ast }.\end{eqnarray}$$
Each of the terms contributes towards particle motion. The first term on the right-hand side causes particles to move if there is a stress gradient acting on that material, whilst the second represents stress-induced diffusion of particles. Under isotropic loading (
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{\ast }=\mathbf{0}$
), it is clearly unphysical for particles to diffuse due to this term, and so we choose to set this second term to zero, as done for immiscible fluid mixtures in Bedford & Drumheller (Reference Bedford and Drumheller1983). This yields the statement of conservation of momentum (2.1) used in the main text.
Appendix B
Conservation of mass, equation (A 3), can be time averaged and integrated over all sizes to give
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\bar{\unicode[STIX]{x1D70C}}\bar{\boldsymbol{u}})=0.\end{eqnarray}$$
Substituting this back into (A 3) and assuming
$\boldsymbol{u}\approx \bar{\boldsymbol{u}}$
, we recover
$\text{D}\unicode[STIX]{x1D719}/\text{D}t=0$
. As we have already assumed an incompressible (isochoric) flow field, we can state that
$\text{D}\unicode[STIX]{x1D70C}/\text{D}t=0$
, which can be time averaged and then integrated over all sizes to yield
$\text{D}\bar{\unicode[STIX]{x1D70C}}/\text{D}t=0$
. Together, these two relations imply that
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\boldsymbol{u}}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}}=0$
. By further assuming that
$\unicode[STIX]{x2202}\hat{\boldsymbol{u}}/\unicode[STIX]{x2202}t=0$
, we can state that
$\text{D}\hat{\boldsymbol{u}}/\text{D}t=0$
. Under these conditions, the left-hand side of (2.5) can be expressed as
$$\begin{eqnarray}\frac{\text{D}(\unicode[STIX]{x1D70C}\hat{\boldsymbol{u}})}{\text{D}t}+\bar{\unicode[STIX]{x1D70C}}\bar{\boldsymbol{u}}\frac{\text{D}\unicode[STIX]{x1D719}}{\text{D}t}=\bar{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D719}\cancel{\frac{\text{D}\hat{\boldsymbol{u}}}{\text{D}t}}+\bar{\unicode[STIX]{x1D70C}}\hat{\boldsymbol{u}}\cancel{\frac{\text{D}\unicode[STIX]{x1D719}}{\text{D}t}}+\unicode[STIX]{x1D719}\hat{\boldsymbol{u}}\cancel{\frac{\text{D}\bar{\unicode[STIX]{x1D70C}}}{\text{D}t}}+\bar{\unicode[STIX]{x1D70C}}\bar{\boldsymbol{u}}\cancel{\frac{\text{D}\unicode[STIX]{x1D719}}{\text{D}t}}=\mathbf{0}.\end{eqnarray}$$
Appendix C
The kinetic stress scaling function assumed here,
$f=\bar{s}_{H}/s$
, reproduces the behaviour noted in Tunuguntla, Weinhart & Thornton (Reference Tunuguntla, Weinhart and Thornton2016), where the relative partial kinetic stress fractions scale in an asymmetric manner with particle size. To make such a comparison, we can investigate the case of a bidisperse grainsize distribution, with two constituents
$\unicode[STIX]{x1D708}=a$
and
$b$
, such that
$\unicode[STIX]{x1D719}(s)=\unicode[STIX]{x1D6F7}_{a}\unicode[STIX]{x1D6FF}(s-s_{a})+\unicode[STIX]{x1D6F7}_{b}\unicode[STIX]{x1D6FF}(s-s_{b})$
, where
$\unicode[STIX]{x1D6F7}^{\unicode[STIX]{x1D708}}$
is the volume fraction of constituent
$\unicode[STIX]{x1D708}$
,
$\unicode[STIX]{x1D6F7}_{b}=1-\unicode[STIX]{x1D6F7}_{a}$
and
$\unicode[STIX]{x1D6FF}$
is the delta function. We can then evaluate both
$f^{a}=\unicode[STIX]{x1D6F7}^{a}f(s=s_{a})$
and
$f^{b}=\unicode[STIX]{x1D6F7}^{b}f(s=s_{b})$
, and plot them against their respective volume fraction, as shown in figure 5. The values of kinetic stress imbalance agree well with measured values, as obtained via discrete element simulations, see figure 8 from Tunuguntla et al. (Reference Tunuguntla, Weinhart and Thornton2016). Values are shown for size ratios up to 1000 for illustrative purposes only. The validity of the scaling law for size ratios larger than those explored here is purely conjecture.
Figure 5. Relative partial stress fractions as a function of volume fraction for bidisperse mixtures. Blue lines represent the smaller phase and red lines the larger phase, with the size contrast between the two being
$s_{b}/s_{a}=1.5,3,10,100$
and
$1000$
.
