2.1. The partially regularized $\mu (I)$ rheology

Consider a body of granular material containing particles of differing sizes, shapes and frictional properties, but with a constant intrinsic grain density $\rho _*$. The solids volume fraction $\varPhi$ is also assumed to be constant, so that the bulk density $\rho = \varPhi \rho _*$ is constant. Mass balance then implies that the bulk flow field $\boldsymbol {u}$ is incompressible,

(2.1)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{equation}

where $\boldsymbol {\nabla }$ is the gradient operator and $\boldsymbol {\cdot }$ is the dot product. Although shear and gravity-driven segregation relies on the formation of void spaces inside a flowing mixture, the bulk solids volume fraction in granular avalanches has been shown to remain approximately constant, and incompressibility is therefore a reasonable assumption (Tripathi & Khakhar Reference Tripathi and Khakhar2011). The momentum balance is

(2.2)\begin{equation} \rho\left(\frac{\partial\boldsymbol{u}}{\partial t} + \boldsymbol{u{\cdot}\nabla u}\right) ={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot}(2\eta\boldsymbol{D}) + \rho\boldsymbol{g}, \end{equation}

where $p$ is the pressure, $\boldsymbol{D}=\boldsymbol{\nabla u}+\boldsymbol{\nabla u}^T$ is the strain rate and $\boldsymbol {g}$ is the constant gravitational acceleration vector. The $\mu (I)$ rheology (GDR MiDi 2004; Jop et al. Reference Jop, Forterre and Pouliquen2006) relates the deviatoric stress to the strain-rate tensor $\boldsymbol {D}$ through a granular viscosity

(2.3)\begin{equation} \eta = \frac{\mu(I) p}{2\|\boldsymbol{D}\|}, \end{equation}

where the second invariant of the strain-rate tensor is

(2.4)\begin{equation} \|\boldsymbol{D}\| = \sqrt{\frac{1}{2} \mathrm{tr}(\boldsymbol{D}^2)}. \end{equation}

In terms of this second invariant, the inertial number becomes

(2.5)\begin{equation} I = \frac{2 d \|\boldsymbol{D}\|}{\sqrt{p/\rho_*}}, \end{equation}

where $d$ is the grain diameter. The $\mu (I)$ rheology is an empirical law that was formulated using dimensional analysis alongside evidence from DEM/DPM simulations and experiments across a range of steady-state monodisperse flow geometries (GDR MiDi 2004). Jop et al. (Reference Jop, Forterre and Pouliquen2005) derived the classical functional form

(2.6)\begin{equation} \mu(I) = \frac{\mu_{s} I_{0}+\mu_{d} I}{I_{0}+I},\end{equation}

from the basal friction measurements of Pouliquen & Forterre (Reference Pouliquen and Forterre2002). This starts at the static friction $\mu _s$ at $I=0$ and asymptotes to the dynamic friction $\mu _d$ for large inertial numbers. The range of inertial numbers over which this transition occurs is controlled by the parameter $I_0$. Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015) found that this formulation was well-posed for inertial numbers in the range $[I_1^N,I_2^N]$ provided that the difference $\mu _d-\mu _s$ was sufficiently large. The lower and upper neutral stability limits $I_1^N$ and $I_2^N$ are found by substituting (2.6) into equation (3.9) of Barker & Gray (Reference Barker and Gray2017) and determining the values of $I$ for which the condition $C=0$. When the inertial number is too low ($I< I_1^N$) or too high ($I>I_2^N$), the system leads to unbounded growth of small perturbations (ill-posedness) in the high wavenumber limit. This leads to grid-dependent results in numerical simulations (Barker et al. Reference Barker, Schaeffer, Bohorquez and Gray2015; Barker & Gray Reference Barker and Gray2017; Martin et al. Reference Martin, Ionescu, Mangeney, Bouchut and Farin2017), and indicates that important physics is missing in the model.

The issue of ill-posed behaviour can be resolved by the inclusion of higher gradients (Goddard & Lee Reference Goddard and Lee2017), non-locality (Kamrin & Koval Reference Kamrin and Koval2012; Kamrin Reference Kamrin2019) or compressibility (Barker et al. Reference Barker, Schaeffer, Shearer and Gray2017; Heyman et al. Reference Heyman, Delannay, Tabuteau and Valance2017; Schaeffer et al. Reference Schaeffer, Barker, Tsuji, Gremaud, Shearer and Gray2019), but these all add further complexity to the equations. A simple, but practical, approach that retains the structure of the incompressible Navier–Stokes equations (2.1)–(2.2), is to modify the shape of the $\mu (I)$ curve to increase the theory's range of well-posedness. In the small inertial number limit, Barker & Gray (Reference Barker and Gray2017) solved for the boundary between well- and ill-posed regions in parameter space to motivate the alternative $\mu (I)$ function,

(2.7)\begin{equation} \mu(I) = \begin{cases} \displaystyle \sqrt{\frac{\alpha}{\log\left(\dfrac{A}{I}\right)}} & \text{for } I \leq I_1^N, \\ \displaystyle \frac{\mu_{s} I_{0}+\mu_{d} I+\mu_{\infty} I^2}{I_{0}+I} & \text{for } I > I_1^N, \end{cases} \end{equation}

where $\alpha \lesssim 2$ ensures that the lower branch is well-posed for $I\in [0,I_1^N]$, $\mu _{\infty }$ is a new material constant and the constant

(2.8)\begin{equation} A=I_1^N\exp\left(\frac{\alpha(I_0+I_1^N)^2}{(\mu_s I_0 + \mu_d I_1^N + \mu_{\infty}(I_1^N)^2)^2}\right) \end{equation}

ensures that the function (2.7) is continuous at $I=I_1^N$. The function (2.7) implies that the friction $\mu =0$ at $I=0$. The partially regularized theory therefore does not have a static yield stress, but instead creeps for inertial numbers $I \leq I_1^N$. This makes little practical difference to numerical simulations, which already necessitate high viscosity regularization (see § 3).

Since $I_1^N\ll 1$, it follows that the upper branch of (2.7) closely approximates the original function (2.6) for $I_1^N< I\ll 1$. At large inertial numbers the upper branch asymptotes to $\mu =\mu _d+\mu _\infty I$, instead of tending to the constant value $\mu _d$. This is consistent with the high speed chute flow experiments of Holyoake & McElwaine (Reference Holyoake and McElwaine2012) and Barker & Gray (Reference Barker and Gray2017), which showed there was a linear dependence on $I$ at large inertial numbers. Such experiments therefore provide a means of determining the new material constant $\mu _\infty$. Barker et al. (Reference Barker, Rauter, Maguire, Johnson and Gray2021) showed that the creep state at small inertial numbers and the linear dependence of $\mu$ on $I$ at large inertial numbers significantly extends the range of well-posedness from $I=0$ to $I^N_2\simeq 17$. It is for this reason that the theory is known as the partially regularized $\mu (I)$ rheology, as it does not guarantee well-posedness, but instead makes the well-posed range much larger than that of the original theory (Jop et al. Reference Jop, Forterre and Pouliquen2006), where $I^N_2\simeq 0.28$ (Barker & Gray Reference Barker and Gray2017).

2.2. Generalized polydisperse segregation theory

The granular material is assumed to consist of $n$ discrete particle species $\nu$ that have differing particle diameters $d^\nu$, but the same intrinsic grain density $\rho _*$. Each species has a volume fraction per unit granular volume $\phi ^{\nu } \in [0,1]$, which satisfies the summation constraint

(2.9)\begin{equation} \sum_{\forall\nu} \phi^{\nu}=1. \end{equation}

All of the existing bidisperse and polydisperse theories for particle-size segregation (e.g. Bridgwater, Foo & Stephens Reference Bridgwater, Foo and Stephens1985; Savage & Lun Reference Savage and Lun1988; Dolgunin & Ukolov Reference Dolgunin and Ukolov1995; Gray & Thornton Reference Gray and Thornton2005; Gray & Chugunov Reference Gray and Chugunov2006; Fan & Hill Reference Fan and Hill2011; Gray & Ancey Reference Gray and Ancey2011; Schlick et al. Reference Schlick, Fan, Umbanhowar, Ottino and Lueptow2015; Gray Reference Gray2018; Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021) can be recast in the form of a general segregation–advection–diffusion equation for each species $\nu$,

(2.10)\begin{equation} \frac{\partial\phi^{\nu}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\phi^{\nu}\boldsymbol{u}) + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{F}^{\nu} = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathcal{D}}^{\nu}, \end{equation}

where $\boldsymbol {F^{\nu }}$ and $\boldsymbol {\mathcal {D}^{\nu }}$ are the segregation and diffusive flux vectors for species $\nu$, respectively. By summing (2.10) over each species and applying the summation constraint (2.9), the incompressibility condition (2.1) is satisfied provided

(2.11a,b)\begin{equation} \sum_{\forall\nu} \boldsymbol{F}^{\nu}=\boldsymbol{0} \quad \text{and} \quad \sum_{\forall\nu} \boldsymbol{\mathcal{D}}^{\nu}=\boldsymbol{0}. \end{equation}

For a bidisperse mixture, the segregation flux function should satisfy the constraint that no segregation occurs when the volume fraction of either species is zero (Bridgwater et al. Reference Bridgwater, Foo and Stephens1985). The simplest form therefore has a linear dependence on the species volume fractions, as proposed for a bidisperse mixture by Gray & Thornton (Reference Gray and Thornton2005). This was later generalized for a polydisperse mixture (Gray & Ancey Reference Gray and Ancey2011; Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021) by a summation of the segregation flux functions over each bidisperse sub-mixture, to give

(2.12)\begin{equation} \boldsymbol{F}^{\nu} = \sum_{\forall\lambda{\ne}\nu}f_{\nu\lambda}\phi^{\nu}\phi^{\lambda}\boldsymbol{e}_{\nu \lambda}, \end{equation}

where $f_{\nu \lambda }$ is the segregation velocity magnitude and $\boldsymbol {e}_{\nu \lambda }$ is the unit vector along the segregation direction. This satisfies the summation constraint (2.11) provided

(2.13a,b)\begin{equation} f_{\nu\lambda}=f_{\lambda\nu}, \quad \boldsymbol{e}_{\nu\lambda}={-}\boldsymbol{e}_{\lambda\nu}. \end{equation}

In Barker et al. (Reference Barker, Rauter, Maguire, Johnson and Gray2021) the diffusive flux vector is defined by analogy with the Maxwell–Stefan equations (Maxwell Reference Maxwell1867)

(2.14)\begin{equation} \boldsymbol{\mathcal{D}}^{\nu} = \sum_{\forall\lambda{\ne}\nu}\mathcal{D}_{\nu\lambda}(\phi^{\lambda}\boldsymbol{\nabla}\phi^{\nu} - \phi^{\nu}\boldsymbol{\nabla}\phi^{\lambda}), \end{equation}

so that, in general, the diffusion rate $\mathcal {D}_{\nu \lambda }$ may vary between sub-mixtures, while still satisfying (2.11), provided $\mathcal {D}_{\nu \lambda }=\mathcal {D}_{\lambda \nu }$. When the diffusion rate is the same across all sub-mixtures, (2.14) reduces to the standard equation for Fickian diffusion.

2.3. Segregation induced feedback on the bulk flow

The inertial number (2.5) is particle-size dependent. Hence, even if the particles have identical microscopic frictional properties, the local friction $\mu$ and, hence, the bulk flow behaviour will depend on the local particle size (GDR MiDi 2004). Particle-size segregation can therefore feedback on the local flow dynamics. By defining a volume fraction weighted mean particle size

(2.15)\begin{equation} \bar d=\sum_{\forall \nu} \phi^\nu d^\nu, \end{equation}

it is possible to generalize the inertial number definition (2.5) to polydisperse mixtures

(2.16)\begin{equation} I = \frac{2\bar d \|\boldsymbol{D}\|}{\sqrt{p/\rho_*}}. \end{equation}

This is sufficient to capture simple frictional feedback from the evolving particle-size distribution of otherwise identical particles (Rognon et al. Reference Rognon, Roux, Naaïm and Chevoir2007; Tripathi & Khakhar Reference Tripathi and Khakhar2011; Denissen et al. Reference Denissen, Weinhart, Te Voortwis, Luding, Gray and Thornton2019; Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021). This implies that, for the same flow depth, small particles will flow faster than large grains, so long as non-local wall effects do not come into play (Pouliquen & Forterre Reference Pouliquen and Forterre2002; Kamrin & Koval Reference Kamrin and Koval2014; Edwards & Gray Reference Edwards and Gray2015; Edwards et al. Reference Edwards, Russell, Johnson and Gray2019; Rocha, Johnson & Gray Reference Rocha, Johnson and Gray2019).

There are more complex effects in which the frictional properties of the particles differ due to their shape and/or microscopic surface properties. For each pure phase $\nu$, this can often be accounted for by modifying the parameters $\mu _{s}^\nu$, $\mu _{d}^\nu$, $\mu _{\infty }^\nu$ and $I_0^\nu$ in the friction law (2.7). The local effective friction of the mixture can then be defined as the volume fraction weighted average of the constituent frictions,

(2.17)\begin{equation} \bar\mu(I)=\sum_{\forall \nu} \phi^\nu \mu^\nu(I). \end{equation}

This leads naturally to a new definition for the granular viscosity (2.3), which replaces $\eta$ in the bulk momentum balance (2.2)

(2.18)\begin{equation} \bar\eta = \frac{\bar\mu(I) p}{2\|\boldsymbol{D}\|}. \end{equation}

The combination of (2.16) and (2.17) allows the frictional differences caused by variations in particle size and surface properties to feedback on the bulk motion. This is crucial for segregation-induced flow fingering, which frequently occurs in geophysical mass flows, and is responsible for longer run out (Pouliquen et al. Reference Pouliquen, Delour and Savage1997; Pouliquen & Vallance Reference Pouliquen and Vallance1999; Iverson & Vallance Reference Iverson and Vallance2001; Iverson Reference Iverson2003; Johnson et al. Reference Johnson, Kokelaar, Iverson, Logan, LaHusen and Gray2012; Woodhouse et al. Reference Woodhouse, Thornton, Johnson, Kokelaar and Gray2012; Kokelaar et al. Reference Kokelaar, Graham, Gray and Vallance2014; Baker et al. Reference Baker, Johnson and Gray2016b; Edwards et al. Reference Edwards, Rocha, Kokelaar, Johnson and Gray2023).

2.4. The feedback of the bulk flow on particle-size segregation and diffusion

Each particle species $\nu$ is transported by the bulk velocity field $\boldsymbol {u}$ in the segregation– advection–diffusion equation (2.10). The particle-size distribution is therefore directly affected by the local velocity $\boldsymbol {u}$. The local shear rate $\dot \gamma =2||\boldsymbol {D}||$ and the pressure $p$ also affect the segregation and diffusion through the functional dependence of the segregation velocity magnitude $f_{\nu \lambda }$ and diffusion rate $\mathcal {D}_{\nu \lambda }$ between each species $\nu$ and $\lambda$ pair.

Barker et al. (Reference Barker, Rauter, Maguire, Johnson and Gray2021) surveyed the existing literature for specific functional forms. In particular, dimensional analysis demands that the diffusion rate scales with the shear rate times the particle-size squared, times an arbitrary function of the inertial number. This scaling was implied by the analysis of Scott & Bridgwater (Reference Scott and Bridgwater1975) on the relationship between particle percolation velocity and diffusion, and has been extensively confirmed by experimental observations (Bridgwater Reference Bridgwater1980; Natarajan, Hunt & Taylor Reference Natarajan, Hunt and Taylor1995; Utter & Behringer Reference Utter and Behringer2004; Katsuragi, Abate & Durian Reference Katsuragi, Abate and Durian2010) and numerical simulations (Tripathi & Khakhar Reference Tripathi and Khakhar2013; Fan et al. Reference Fan, Schlick, Umbanhowar, Ottino and Lueptow2014; Cai et al. Reference Cai, Xiao, Zheng and Zhao2019). The simplest model for the diffusion rate between each species $\nu$ and $\lambda$ is therefore

(2.19)\begin{equation} \mathcal{D}_{\nu\lambda}= 2\mathcal{A}\|\boldsymbol{D}\| \bar d^2, \end{equation}

where $\mathcal {A}$ is a universal constant specified in table 1. Note that the diffusion is therefore assumed to be the same for each species pair, and hence, the diffusive flux (2.14) reduces to standard Fickian diffusion. Campbell (Reference Campbell1997) found that diffusion in granular shear flows was anisotropic, and Utter & Behringer (Reference Utter and Behringer2004) studied monodisperse Couette flow experiments using flat disc-shaped particles, from which they measured a radial diffusion coefficient of $\mathcal {A}=0.108$ and a tangential diffusion coefficient of $\mathcal {A}=0.223$. Evidence from bidisperse DEM/DPM shear-cell simulations with spherical particles implies that the anisotropy is only slight, and that instead $\mathcal {A}\in (0.03,0.05)$ (Tripathi & Khakhar Reference Tripathi and Khakhar2013; Cai et al. Reference Cai, Xiao, Zheng and Zhao2019; Artoni et al. Reference Artoni, Larcher, Jenkins and Richard2021; Bancroft & Johnson Reference Bancroft and Johnson2021). Therefore isotropic diffusion can be reasonably assumed, and $\mathcal {A}=0.04$ is used here.

Table 1. Parameter values used in the numerical simulations. The value of the universal constant $\mathcal {A}$ in the diffusivity (2.19) is based on the DEM/DPM simulations of (Tripathi & Khakhar Reference Tripathi and Khakhar2013; Cai et al. Reference Cai, Xiao, Zheng and Zhao2019; Artoni et al. Reference Artoni, Larcher, Jenkins and Richard2021; Bancroft & Johnson Reference Bancroft and Johnson2021). The values of the universal constants $\mathcal {B}$ and $\mathcal {C}$ in the segregation law (2.20) are based on those suggested by the single-intruder experiments of Trewhela et al. (Reference Trewhela, Ancey and Gray2021), although $\mathcal {B}$ is corrected for the absence of an interstitial fluid. The remaining non-dimensional constants $\mathcal {E}$, $a$ and $\phi ^s_c$ are also taken from Trewhela et al. (Reference Trewhela, Ancey and Gray2021).

Trewhela et al. (Reference Trewhela, Ancey and Gray2021) performed single-intruder shear-box experiments to derive a scaling law for segregation of large and small particles in a bidisperse mixture. This law is now generalized to polydisperse systems. For a sub-mixture composed of particles $\nu$ and $\lambda$, the segregation velocity magnitude is

(2.20)\begin{equation} f_{\nu\lambda}=\left(\frac{2\mathcal{B} \rho_* g \|\boldsymbol{D}\| \bar d^2 }{\mathcal{C}\rho_* g \bar d + p}\right) \mathcal{F}_{\nu\lambda}(\phi^\nu,\phi^\lambda,R_{\nu\lambda}), \end{equation}

where $\mathcal {B}$ and $\mathcal {C}$ are universal constants determined by Trewhela et al. (Reference Trewhela, Ancey and Gray2021) (given in table 1). The function $\mathcal {F}_{\nu \lambda }$ accounts for packing effects and the segregation-rate dependence on the generalized grain-size ratio

(2.21)\begin{equation} R_{\nu\lambda}=\hbox{max}\left(\frac{d^\nu}{d^\lambda},\frac{d^\lambda}{d^\nu}\right), \end{equation}

where $R_{\nu \lambda }=R_{\lambda \nu }>1$. The scaling law (2.20) is able to capture a variety of phenomena observed in slightly dilated, bidisperse, sheared granular flows, where kinetic sieving and squeeze expulsion are the dominant segregation mechanisms (Middleton Reference Middleton1970; Bridgwater et al. Reference Bridgwater, Foo and Stephens1985; Savage & Lun Reference Savage and Lun1988; Gray Reference Gray2018).

Kinetic sieving is a shear driven process that allows smaller particles to percolate down under the action of gravity, while squeeze expulsion describes the process in which particles are squeezed upwards to maintain bulk incompressibility. Dimensional analysis suggests that the segregation rate should have the dimensions of a velocity, and (2.20) assumes that $f_{\nu \lambda }$ behaves like $\dot \gamma \bar d=2||\boldsymbol {D}||\bar d$, which automatically implies that the segregation flux function $f_{\nu \lambda }\phi ^\nu \phi ^\lambda$ is asymmetric (Gajjar & Gray Reference Gajjar and Gray2014; van der Vaart et al. Reference van der Vaart, Gajjar, Epely-Chauvin, Andreini, Gray and Ancey2015). By performing experiments with a density matched interstitial fluid, Vallance & Savage (Reference Vallance and Savage2000) showed that segregation shuts off in the absence of gravity (see also Thornton, Gray & Hogg Reference Thornton, Gray and Hogg2006). The effect of gravity is included indirectly in the round bracketed term in (2.20) through a dependence on the reciprocal of the non-dimensionalized pressure $p/(\rho _* g\bar d)$. If the shear rate is constant, this form implies that the deeper one goes into the mixture, the lower the segregation rate becomes. This was observed directly in Trewhela et al. (Reference Trewhela, Ancey and Gray2021) experiments, i.e. the intruders did not rise linearly with time in response to a spatially uniform shear rate, but were either percolated down or squeezed up faster near the free surface. The constant $\mathcal {C}$ in (2.20) was included to avoid a singularity at zero pressure. Further evidence of the pressure dependence of segregation is provided in the papers of Golick & Daniels (Reference Golick and Daniels2009), Fry et al. (Reference Fry, Umbanhowar, Ottino and Lueptow2018) and Bancroft & Johnson (Reference Bancroft and Johnson2021).

Trewhela et al. (Reference Trewhela, Ancey and Gray2021) found that for moderate grain-size ratios of large and small particles (i.e. for $R_{sl}\in [1,4.17]$, where $s$ and $l$ denote small and large particles, respectively), the segregation rate of large intruders had a linear dependence on $(R_{sl}-1)$. However, the segregation rate of a small intruder obeyed a quadratic law of the form $(R_{sl}-1)+\mathcal {E}(R_{sl}-1)^2$, where $\mathcal {E}$ was a non-dimensional constant. The additional quadratic dependence arose because small intruders find it increasingly easy to percolate through a matrix of large grains as the particle-size ratio approaches the spontaneous percolation limit at around $R_{sl}\simeq 6$. Trewhela et al. (Reference Trewhela, Ancey and Gray2021) extended their model to intermediate concentrations by assuming that

(2.22)\begin{equation} \mathcal{F}_{sl}=(R_{sl}-1)+\mathcal{E}\varLambda(\phi^s)(R_{sl}-1)^2, \end{equation}

where the function

(2.23)\begin{equation} \varLambda(\phi^s)=\begin{cases} \displaystyle 1-\frac{\phi^s}{\phi^s_c} & \text{for } \phi^s \leq \phi^s_c, \\ \displaystyle 0 & \text{for } \phi^s > \phi^s_c, \end{cases} \end{equation}

gives the right $\mathcal {F}_{sl}$ dependence in the limit as $\phi ^s\rightarrow 0,1$ and shuts off the quadratic dependence when the small-particle concentration exceeds $\phi ^s_c$ (defined in table 1). Trewhela et al. (Reference Trewhela, Ancey and Gray2021) showed that this functional form was capable of quantitatively capturing the spatial and temporal evolution of the particle-size distribution measured in the refractive-index-matched shear-box experiments of van der Vaart et al. (Reference van der Vaart, Gajjar, Epely-Chauvin, Andreini, Gray and Ancey2015).

However, there is evidence from experiments and DEM/DPM simulations suggesting that for 50:50 mixtures of large and small particles in geometries for which the velocity and shear rate are functions of space and time, the segregation intensity is maximal near a grain-size ratio of $R_{sl}=2$ (Golick & Daniels Reference Golick and Daniels2009; Thornton et al. Reference Thornton, Weinhart, Luding and Bokhove2012), whereas the formulation (2.22) is monotonically increasing with the size ratio. Trewhela et al. (Reference Trewhela, Ancey and Gray2021) therefore suggested an alternative size-ratio dependency that captures this effect, within the framework of the incompressible segregation theory

(2.24)\begin{equation} \mathcal{F}_{sl}=\frac{(R_{sl}-1)+\mathcal{E}\varLambda(\phi^s)(R_{sl}-1)^2}{1 + a(R_{sl}-1)^2\phi^s\phi^l}, \end{equation}

where $a$ is a constant defined in table 1. The numerator in (2.24) is the same as (2.22), but the denominator now includes a reduction factor $1 + a(R_{sl}-1)^2\phi ^s\phi ^l$, which shuts off when $\phi ^s=0$ and $\phi ^l=0$, and reduces the segregation rate at intermediate configurations. The idea is that this factor parameterizes the enhanced packing of mixtures of particles of different sizes, which is thought to reduce the segregation rate. Trewhela et al. (Reference Trewhela, Ancey and Gray2021) showed that for a 50:50 mix, (2.24) predicts a peak segregation rate at $R_{sl}\simeq 1.66$, which is very close to the maximum at $R_{sl}=1.7$ observed in the DEM/DPM simulations of Thornton et al. (Reference Thornton, Weinhart, Luding and Bokhove2012). The bidisperse simulations in §§ 4–7 of this paper are performed using the functional form (2.24). Almost nothing is known about how the segregation rate and the packing behaves in polydisperse sheared granular mixtures. The tridisperse simulations presented in § 8 therefore use a trivial pairwise generalization of (2.24) as an illustration.

2.5. Sidewall friction

It has long been known that confining lateral sidewalls play an important role in granular flow experiments (Greve & Hutter Reference Greve and Hutter1993; Taberlet et al. Reference Taberlet, Richard, Valance, Losert, Pasini, Jenkins and Delannay2003; Jop et al. Reference Jop, Forterre and Pouliquen2005; Baker, Barker & Gray Reference Baker, Barker and Gray2016a). In rotating drum flows the additional friction introduced by sidewalls results in avalanches that are significantly thinner and faster than flows without sidewall friction (Hill et al. Reference Hill, Khakhar, Gilchrist, McCarthy and Ottino1999; Ottino & Khakhar Reference Ottino and Khakhar2000; Jop et al. Reference Jop, Forterre and Pouliquen2005; Mounty Reference Mounty2007). As a result, the two-dimensional square rotating drum simulations of Barker et al. (Reference Barker, Rauter, Maguire, Johnson and Gray2021) could not be compared directly to experiments. In principle, a Coulomb friction boundary condition could be imposed on the sidewalls and then the governing equations could be solved in three dimensions for any width of drum. However, in this paper the influence of sidewall friction is incorporated into the governing equations by width averaging the momentum balance equation (2.2). This reduces the problem to a two-dimensional one, which is much more computationally efficient.

Consider then a three-dimensional granular flow confined within a narrow channel between $y=0$ and $y=W$, where there can be slip at the sidewalls and only weak velocity gradients in the $y$ direction. Proceeding in a similar manner to Jop et al. (Reference Jop, Forterre and Pouliquen2005), assuming a Coulomb friction boundary condition at the lateral sidewalls implies

(2.25)\begin{equation} \boldsymbol{\tau n} ={-}\mu_W p\frac{\boldsymbol{u}}{|\boldsymbol{u}|}, \quad \text{at}\ y=0,W, \end{equation}

where $\boldsymbol {n}$ is an outward pointing normal to the sidewall, $\mu _W$ is a constant wall friction coefficient and the term $-\boldsymbol {u}/|\boldsymbol {u}|$ ensures that friction acts against the flow. The three-dimensional mass and momentum balance equations, defined analogously with (2.1) and (2.2), may then be integrated across the channel width in a similar fashion to the depth-averaged approach for a shallow system (Gray Reference Gray2001; Gray & Edwards Reference Gray and Edwards2014). Width-integrated variables are defined by

(2.26)\begin{equation} \tilde f = \frac{1}{W}\int^W_0 f\,\text{d}y, \end{equation}

for any variable $f$, and when variations in the $y$ direction are small, $\widetilde {fg}\simeq \tilde f\tilde g$ for any two variables $f$ and $g$. Integrating and applying the boundary conditions (2.25) results in a two-dimensional system for which the mass and momentum balances become

(2.27)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}=0, \end{gather}$$

(2.28)$$\begin{gather}\rho\left(\frac{\partial\tilde{\boldsymbol{u}}}{\partial t} + \tilde{\boldsymbol{u}}\,\boldsymbol{\cdot}\boldsymbol{\nabla}\tilde{\boldsymbol{u}}\right) ={-}\boldsymbol{\nabla} \tilde p + \boldsymbol{\nabla}\boldsymbol{\cdot}(2\tilde\eta\tilde{\boldsymbol{D}}) + \rho\boldsymbol{g} - \frac{2}{W}\mu_W p\frac{\tilde{\boldsymbol{u}}}{|\tilde{\boldsymbol{u}|}}, \end{gather}$$

where $\boldsymbol {\nabla }=(\partial /\partial x, \partial /\partial z)^T$, the two-dimensional gradient operator. Since there is no lateral flow and small velocity gradients in the $y$ direction, $\boldsymbol {u}/|\boldsymbol {u}|\simeq \tilde {\boldsymbol {u}}/|\tilde {\boldsymbol {u}}|$, and the tilde notation may be dropped on the understanding that the variables $(\boldsymbol {u}, p, \eta, \boldsymbol {D})$ from this point onwards refer to width-averaged quantities of the form (2.26). The final term in the momentum balance equation (2.28) models the influence of sidewall friction on the bulk flow. This alteration provides a straightforward method of capturing three-dimensional wall friction effects in a two-dimensional framework.

4.1. Experimental set-up

Experiments were conducted in a triangular rotating drum with equilateral sides of length $L=0.257$ m. This was formed from an aluminium outer frame that was confined between two transparent polymethyl methacrylate (PMMA) plates (separated by a gap of width $W=3\times 10^{-3}$ m) and screwed together. Before assembly, the PMMA walls were cleaned with anti-static spray and dried to prevent smaller particles sticking to the sidewalls. The glass beads used in the experiments were sourced from Sigmund Lindner GmbH. The larger green glass beads were sieved to a diameter $d^l = 600\unicode{x2013}800$ $\mathrm {\mu }$m, while the smaller red beads had a diameter $d^s = 300\unicode{x2013}400$ $\mathrm {\mu }$m. The desired volumes of large and small grains were measured separately, mixed together and then filled into the drum through a gap in the aluminium frame, which was then tightly closed with a stopper. The exact volume of grains in the drum can vary slightly from the individual combined large and small volumes due to enhanced packing of the mixture (Golick & Daniels Reference Golick and Daniels2009; Thornton et al. Reference Thornton, Weinhart, Luding and Bokhove2012). If necessary, additional grains were added to obtain the desired fill. For the particle-size ratios used in this paper these packing effects were small, and the theory, which is incompressible, appears to be able to capture the observed dynamics and the segregation patterns to a high degree of accuracy.

Particles of different sizes have a strong propensity to segregate. In order to obtain an approximately homogenous initial condition, the drum was oriented horizontally and shaken to induce segregation in the direction normal to the lateral transparent sidewalls. The large green beads therefore segregated to the front (observable) wall, with the small red beads concealed behind them against the rear wall. At any given cross-section of the drum an approximately uniform initial mixture of large and small particles was generated. The back of the drum was then attached to an ElectroCraft S642-1B/T servo motor and a modulation speed control unit to control the rotation rate. The axis of rotation coincided with the centre of the triangle. The drum was narrow enough that cross-slope particle diffusion, within the flowing avalanche, was sufficient to smooth out any lateral variation across it. The final pattern was therefore approximately uniform across the gap, except in the central core, which never entered the avalanche. This initial condition therefore has the additional advantage that the undisturbed core appears to the observer from the front as a uniform region of green particles, providing a strong contrast with the predominantly small red particles that are deposited adjacent to the core. The sidewall friction $\mu _W=\tan (15.5^\circ )$ is a material property of both the glass beads and the PMMA, and is measured empirically from the angle of failure for a bidisperse mixture of static grains on a gradually inclined PMMA surface (see table 2).

Figure 2 shows a series of photographs of the evolving pattern that develops when a 50:50 bidisperse mixture of large and small grains is rotated in a 70 % filled triangular drum rotating at 2.5 rpm. Supplementary movie 2 shows the full time-dependent evolution. The granular free surface, which was initially horizontal, inclines until a critical angle is exceeded, and a thin surface avalanche forms. The particles segregate very efficiently from one another in the surface avalanche, with the large particles rising to the surface and the small particle percolating down to the base, where they are first to deposit into the underlying solid rotating body of grains. As a result, the large green grains end up adjacent to the drum wall and the smaller red particles deposit closer to the central core. The complex flow field allows a quasi-periodic pattern to develop after just two rotations, which has small-particle rich arms that radiate outwards and point towards the corners of the triangle. A full discussion of the experimental results is deferred until the particle-size distribution of the numerical simulation is introduced in § 4.4.

Figure 2. Experimental photos of a 70 % filled drum with a 50:50 mixture of large green ($d^l = 600\unicode{x2013}800$ $\mathrm {\mu }$m) and small red ($d^s = 300\unicode{x2013}400$ $\mathrm {\mu }$m) glass beads rotating at 2.5 rpm. Solid-body rotation occurs in the first 2 seconds. The sequence starts at $t=2$ s and progresses in increments of 2 s (30$^\circ$ rotation), from left to right, top to bottom, up to $t=48$ s (two revolutions). Once the surface avalanche develops, the mixture gradually segregates into a structure of small-particle arms. The protruding stopper provides a reference point for the rotation. Supplementary movie 2 shows the full time-dependent evolution and is available in the online supplementary material.

4.2. Coordinate system, boundary and initial conditions for the simulations

In order to simulate the segregation of the grains in the experimental rotating drum, a fixed Cartesian coordinate system $Oxz$ is defined with the origin $O$ at the axis of rotation and the $z$ coordinate pointing in the opposite direction to gravity. The large- and small-particle diameters are assumed to take average values $d^l = 700$ $\mathrm {\mu }$m and $d^s = 350$ $\mathrm {\mu }$m. The large particles are therefore twice the size of the small grains. The drum is filled to 70 % of its volume with an initially homogeneous 50:50 mixture of large and small grains, so that $\varphi ^s=\varphi ^l=0.5$ everywhere within the granular material. In addition, at $t=0$ all the material is assumed to be in solid-body rotation with velocity

(4.1)\begin{equation} \boldsymbol{u} = \varOmega r\boldsymbol{\theta}, \end{equation}

where $\varOmega$ is the rotation rate, $r$ is the radial coordinate and $\boldsymbol {\theta }$ is the azimuthal unit vector (Gray Reference Gray2001). A rotation rate of 2.5 rpm ($\varOmega =-{\rm \pi} /12$ rad s$^{-1}$) is imposed, which places the drum in the rolling, or continuously avalanching, regime (Henein et al. Reference Henein, Brimacombe and Watkinson1983; Rajchenbach Reference Rajchenbach1990; Gray Reference Gray2001; Mellmann Reference Mellmann2001; Ding et al. Reference Ding, Forster, Seville and Parker2002; Yang et al. Reference Yang, Yu, McElroy and Bao2008). This implies that when the slope exceeds a critical angle, a liquid-like free-surface avalanche forms that continuously erodes and deposits grains with an underlying solid rotating granular body beneath.

It is useful to define the relative velocity to the drum as

(4.2)\begin{equation} \hat{\boldsymbol{u}}=\boldsymbol{u} - \varOmega r\boldsymbol{\theta}. \end{equation}

No-slip and no-penetration conditions are imposed on the triangular walls which, using the relative velocity, can be expressed as $\hat {\boldsymbol {u}}=\boldsymbol {0}$ on the drum walls. These conditions are applied on the rotating mesh using OpenFOAM's mesh-motion routines, with a structured triangular mesh containing $N^2=600^2$ cells. The parameters for segregation and diffusion are the universal constants summarized in table 1, while the frictional and remaining parameters for the glass beads and sidewall friction are specified in table 2. There are therefore no fitting parameters used in the simulations.

4.3. Bulk flow solutions and comparison to experiments

In the experiments and numerical simulations, the granular free surface is initially horizontal. As the drum rotates the granular material performs solid-body rotation until a critical angle is reached, and a liquid-like avalanche develops adjacent to the free surface. Due to sidewall friction, the critical free-surface inclination is significantly higher than the static friction angle $\zeta _s=\tan ^{-1}(\mu _s)$, which is usually the approximate angle of failure (Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021). The free-surface inclination remains relatively constant, with only subtle periodic variations as shown in figure 3 and supplementary movies 3–5. The entire flow then exhibits a quasi-periodic pulsating behaviour every 120$^\circ$ of rotation (8 s), due to the rotation of the drum and the continuously changing drum cross-section occupied by the grains.

Figure 3. Quasi-periodic motion repeating every 120$^\circ$ of revolution for a 70 % full triangular drum computed using the parameters from tables 1 and 2. For illustration, each field is shown at $t=42$, 44, 46 and 48 s, from left to right, corresponding to one third of a revolution. A black dot is plotted on one corner of the drum frame so that the changing orientation may be easily tracked. Panel (a) shows the modulus of the relative velocity field, (b) the pressure and (c) the base 10 logarithm of the inertial number. The white dashed line in (c) represents the position below which the high viscosity cutoff (3.9) becomes active. The full time-dependent evolution of each of the fields is shown in supplementary movies 3–5 (available online).

There can be subtle feedback effects from the evolving grain-size distribution, but they are not large at 70 % fill (Hill et al. Reference Hill, Khakhar, Gilchrist, McCarthy and Ottino1999; Mounty Reference Mounty2007). The simulation reproduces very closely the free-surface inclination observed experimentally (figure 2), including the slight S-shape, which becomes more pronounced when the free surface is longest, and the subtle inclination dip (at $t=42$ s in figure 3), which occurs just past the halfway point of the avalanche. Feedback between the bulk flow and the mixture composition plays a more pivotal role at fill levels near 50 %, when the strongly composition-dependent velocity field can lead to the formation of petal-like patterns in the particle distribution (Zuriguel et al. Reference Zuriguel, Gray, Peixinho and Mullin2006) for circular drums, and more complex patterns in triangular and square drums (Hill et al. Reference Hill, Khakhar, Gilchrist, McCarthy and Ottino1999; Khakhar et al. Reference Khakhar, McCarthy, Gilchrist and Ottino1999; Ottino & Khakhar Reference Ottino and Khakhar2000; Mounty Reference Mounty2007).

The length of the avalanche depends on the orientation of the drum. The orientations that generate the longest avalanches also produce the deepest and fastest flowing motion (see times $t=46$ and 48 s in figure 3a). The longer avalanche length means that more particles are entrained into the avalanche, which increases the flux and the peak velocities that occur about halfway down the avalanche. In the lower reach of the avalanche the particles are deposited back into the solid rotating body and the avalanche thins and decelerates. The super inclination of the slope, combined with acceleration down the increased avalanche length, may also contribute to the peak velocity.

The active avalanching layer is very thin, and is close to the depth observed in the experiments shown in figure 2 and supplementary movie 2. Below the avalanching layer the relative velocity $\hat {\boldsymbol {u}}$ rapidly decays into a quasi-static creep state. The pressure (figure 3b) has a lithostatic component, but with increasing depth the pressure gradient becomes increasingly gravity aligned. Since both the segregation velocity magnitude and the inertial number are strain-rate dependent, the inertial number (figure 3c) gives an accurate identification of the active layer in which segregation takes place. It appears slightly deeper than that implied by figure 3(a) since $\log _{10}(I)$ is plotted. It exhibits less variation in magnitude than $|\hat {\boldsymbol {u}}|$ and $p$ with changing orientation, in part due to the chosen logarithmic scale, but also because the strain-rate dependent inertial number is tied to velocity gradients rather than velocity. The white dashed line in figure 3(c) indicates the threshold for the high viscosity cutoff (3.9). Above the white dashed line the dynamic and creep regimes of the partially regularized rheology (2.7) are active, while below it the high Newtonian viscosity is applied. Since this quasi-static creeping region is predominantly in rigid body rotation, it has little influence on the overall flow dynamics.

4.4. Particle-size distribution

Figure 4 shows a sequence of images of the computed time-dependent evolution of the small-particle concentration in the 70 % filled drum with a 50:50 mix of large and small particles. Each image is directly comparable to the experimental photos in figure 2, and is separated by two seconds, or equivalently, 30$^\circ$ angular increments, beginning at $t=2$ s and ending after two revolutions. Supplementary movie 6 shows the full time-dependent evolution. In regions of solid-body-like rotation, the shear rate is very close to zero, so the scaling laws, (2.19) and (2.20) imply that both the diffusion and the segregation are negligibly small. The grains therefore stay at $\varphi ^s=0.5$ during their initial phase of rotation. The surface avalanche is, however, very efficient at segregating differently sized particles, because the shear rate is high and the pressure is low. This implies that the segregation velocity magnitude (2.20) is much higher than in the solid rotating body, and the grains can segregate. As the particles are deposited along the lower reach of the avalanche into the solid rotating body of grains beneath, the segregation rate once again shuts off, and the deposited particle-size distribution is rotated with the drum. The active region of segregation is therefore confined to a narrow boundary layer close to the free surface, i.e. the avalanche. The highest rate of segregation occurs adjacent to the free surface, as observed experimentally here (see § 4.5) and elsewhere (Gray & Hutter Reference Gray and Hutter1997; Khakhar, McCarthy & Ottino Reference Khakhar, McCarthy and Ottino1997). This effect is intensified by the fact that particles near the free surface have longer trajectories through the actively segregating layer before being deposited, and hence, are subject to more accumulated shear.

Figure 4. Small-particle concentration in a 70 % filled triangular rotating drum with a 50:50 mixture of large and small particles, computed using the parameters specified in tables 1 and 2. The simulation is shown over two full revolutions, beginning at $t=2$ s and progressing in increments of 2 s, or 30 $^\circ$ angular increments, from left to right, top to bottom. The excess air above the grains is coloured white. The images are directly comparable to the experimental photos in figure 2. A small black dot in one corner of the triangular frame provides a reference point for the clockwise rotation. Supplementary movie 6 shows the full time-dependent evolution, and this is compared side-by-side to the experiment in supplementary movie 7.

Since the fill level of the drum exceeds 50 %, a central core forms that is never entrained into the avalanche and remains at its initial concentration (Gray Reference Gray2001; Mounty Reference Mounty2007). The central core is a three-sided shape with curved sides, reminiscent of a Reuleaux triangle. The size of the core is in good agreement with the experiment (figure 2), which is a direct result of (i) the sidewall friction making the avalanche thinner and faster than equivalent flows without sidewall friction, and (ii) the segregation (2.20) and diffusion (2.19) laws confining the particle redistribution to the avalanching region. In contrast, in Barker et al. (Reference Barker, Rauter, Maguire, Johnson and Gray2021) simulations (which neglected sidewall friction) the avalanche was deeper and slower, and consequently, both the central core size and the segregation were underestimated compared with experimental observations (Hill et al. Reference Hill, Khakhar, Gilchrist, McCarthy and Ottino1999; Ottino & Khakhar Reference Ottino and Khakhar2000; Mounty Reference Mounty2007).

Over the first full rotation, the material separates out around the homogeneous core into distinct regions of predominantly large or small particles with a diffuse region separating them. The larger particles, that segregate to the surface of the avalanche, deposit close to the drum wall, while the smaller particles at the base of the avalanche, deposit adjacent to the central core. The combination of the complex time-dependent bulk flow (described in § 4.3) and the ordering mechanism of particle-size segregation then begins to generate a pattern that has small-particle rich arms. The first arm emerges after one complete rotation of the drum (figure 4), and the second arm then follows one third of a turn later. The time scale for the gradual formation of the pattern is captured by the numerical simulation, as shown by the excellent frame-by-frame match between figures 2 and 4, and supplementary movie 7. Over the second rotation of the drum, the arms pass through the avalanche again, the segregation becomes stronger and the pattern becomes more clearly defined. On each third of a revolution the pattern then approximately repeats. Simulations are therefore only performed for the first two drum revolutions.

The pattern exhibits a strong three-fold symmetry. To see this more clearly, it is useful to plot the $\varphi ^s=0.5$ contour, and then replot it twice more rotated by $\pm 2{\rm \pi} /3$ radians. This is done in figure 5 at $t=36$ and 44 s. In both cases this produces a triskelion-like pattern with three arms that start adjacent to the central core and radiate outwards in a clockwise sense ultimately pointing towards the corners of the drum. The outer length of the arm is roughly parallel to the adjacent wall of the triangle, culminating in a flattened top. In reality, the full triskelion structure cannot be seen, because one arm necessarily lies above the free surface, but it is interesting to observe the three-fold symmetry using this approach. At $t=36$ s, the symmetry is almost perfect, because each arm has passed through the avalanche an equal number of times. At $t=44$ s, one of the arms has passed through the avalanche an extra time, so the symmetry is not as good. The arm that has passed through the avalanche an extra time has a large rich crevice that penetrates closer to the central core, elongating the small-particle rich arm. This emphasizes the fact that the pattern is continuing to evolve, with the best symmetries obtained when all the arms have passed through the avalanche the same number of times. The same symmetry is also evident in the experiments in figure 2, despite it being harder to generate a spatially homogeneous initial mixture.

Figure 5. The rotational symmetry of the small-particle rich arms is shown by plotting the $\varphi ^s=0.5$ contour for the numerical simulations of a 70 % filled triangular rotating drum containing a 50:50 bidisperse mixture at $t=36$ s (a) and $t=44$ s (b). In both cases the original contour is plotted in black, the same contour rotated by $-2{\rm \pi} /3$ rad (i.e. rotation in the clockwise direction) in red and the contour rotated $2{\rm \pi} /3$ rad (anticlockwise) in blue. At $t=44$ s, one arm has been entrained through the avalanche an additional time, leading to a crevice in this arm near the central core.

4.5. Quantitative analysis and numerical convergence

To provide a quantitative comparison between the simulations and the experiments, this paper focuses on the small-particle concentration field. This can be approximately determined from experimental images using the following method. First, the drum is filled with a mixture of known small-particle concentration, shaken to achieve approximate uniformity and then photographed. The drum, containing this same mixture, is repeatedly shaken and photographed to reduce noise. This process is then repeated for a variety of mean small-particle concentrations across the range $\bar \varphi ^s\in [0,1]$. All the images are then cropped to remove the non-granular regions, and the colour profiles of the resulting raster images are processed in MATLAB to return a matrix giving the RGB (red-green-blue) intensities for each individual pixel. The ratio of the red intensity to the green intensity is determined for each pixel, and averaged over the entire domain. The red and green intensities exhibit a strong correlation with the mean small-particle concentration since they correspond to the colours of the small and large particles, respectively. This process returns a mean red to green pixel intensity ratio corresponding to each measured mean (over the granular domain) small-particle concentration $\bar \varphi ^s$, and the results for a given concentration can be reliably reproduced experimentally. The ratio can then be plotted against $\varphi ^s\in [0,1]$ and a best fit curve drawn to give a complete set of unique expected intensity ratios for every possible concentration. In fact, there is an approximately linear relationship between the intensity ratio and the particle concentration, and consequently, the mapping can be further improved.

Once the relationship between the intensity ratio and the concentration is assumed to be linear, the appropriate linear fit can be determined using data taken from only two concentrations. When the mean small-particle concentration $\bar \varphi ^s=0$ or 1, it follows that $\varphi ^s=\bar \varphi ^s$ at every point in the granular mixture. Therefore, the images with these two mean particle concentrations can be divided into $8\times 8$ pixel cells, where the cell size is chosen to ensure a reasonable number of particles per cell, and the concentration in every cell is known a priori. The mean colour intensity ratio in each individual cell is then calculated, and the mean cell intensities when $\varphi ^s=0$ and 1 can be used to produce the linear fit between intensity ratio and small-particle concentration. This method has three major advantages over using a greater number of points for $\bar \varphi ^s$ and taking the mean pixel intensity ratio over the entire drum. First, since it uses data from $8\times 8$ pixel cells, it can be used to derive the concentration in cells of the same size from other drum images, without any additional assumption of validity. Second, since the concentration in every cell is known when $\bar \varphi ^s=0$ or 1, the standard deviation of the cell intensity ratio can be trivially calculated and used to plot a set of worst fit lines that facilitate accurate error estimation. Error is due to particle shadows and uneven light reflections (Hill et al. Reference Hill, Khakhar, Gilchrist, McCarthy and Ottino1999), which may be reduced, but not eradicated, by an appropriately positioned lighting set-up. Finally, when the drums are filled with pure phases of small or large particles, attempts at attaining a relatively homogeneous mixture are not subject to inadvertent segregation and so the colour intensity data are more reliable than for intermediate particle concentrations.

Using the calculated mapping from colour intensity to concentration, images of the developing rotating drum patterns from figure 2 are divided into $8\times 8$ pixel cells, so that the mean red to green intensity ratio can be determined for the individual cell and a small-particle concentration value assigned using the set of expected intensity ratios. To account for the undisturbed core, which appears to have $\varphi ^s=0$ but in reality has $\varphi ^s=0.5$ due to the method of attaining the initial mixture composition (see § 4.1), this region is filled with black pixels for which an exception is made so that they are assigned the correct concentration value. Note that although some pixels elsewhere in the drum appear very dark due to gaps between particles adjacent to the transparent walls, over the $8\times 8$ pixel cells there is sufficient saturation of colour to avoid confusion with the central core region. This results in a projected small-particle concentration field for the entire drum. Using this method, the mean particle concentration of the experiments using an approximately 50:50 mixture of large and small particles is calculated to be $\bar \varphi ^s=0.49\pm 0.03$, matching very closely with the approximate true value of $\bar \varphi ^s=0.5$.

The concentration field for the 70 % filled rotating drum after $t=48$ s is shown in figure 6. The number of pixels per cell is chosen to give a detailed concentration field with strong image definition where each cell contains many particles. The subtle region of weaker small-particle concentration above the core in the crevice of the arm being entrained into the avalanche is more clearly visible in this image than in the pre-processed image, which is shown in the final panel of figure 2. As predicted by the theory, the predominantly large-particle regions adjacent to the drum walls are more strongly segregated than the small-particle regions adjacent to the central core, since segregation is strongest at the free surface, towards which the large particles are segregated. This was not evident in the pre-processed images, but has been verified by examining the field $(\varphi ^s-\bar \varphi ^s)^2$, which confirms that the predominantly large-particle regions deviate further from the mean particle concentration than the predominantly small regions.

Figure 6. Approximated small-particle concentration field at $t=48$ s for the experiment using a triangular drum 70 % filled with an approximately 50:50 mixture of large green and small red glass beads. The small-particle concentration is obtained using the projection method described in § 4.5.

To test the strength of the segregation and compare it directly to the simulation data, the segregation intensity is defined analogously to Danckwerts (Reference Danckwerts1952), as the standard deviation of $\varphi ^s$ normalised by the mean small- and large-particle concentrations, $\bar \varphi ^s(1-\bar \varphi ^s)$,

(4.3)\begin{equation} \mathcal{S}=\sqrt{\frac{\displaystyle\iint_\omega(\varphi^s-\bar\varphi^s)^2 \,\text{d}{\kern1pt}x\,\text{d}z}{\displaystyle\bar\varphi^s(1-\bar\varphi^s)\iint_\omega \,\text{d}{\kern1pt}x\,\text{d}z}}, \end{equation}

where $\omega$ is the part of the domain occupied by granular material. This implies that $\mathcal {S}=0$ for a homogeneously mixed material, and $\mathcal {S}=1$ when it is fully segregated. The computed and experimental segregation intensity are plotted as a function of time and at various grid resolutions in figure 7. Sidewall friction results in a free-surface avalanche that acts as a very thin boundary layer where all the segregation occurs. To correctly predict the segregation and, hence, the overall pattern formation, this boundary layer must be adequately resolved, which requires a very fine numerical mesh. For coarser meshes, the avalanche is under resolved and numerical diffusion therefore leads to an underestimation of the segregation intensity. Although the finest resolution, which uses 600$^2$ triangular cells, is not perfectly resolved, the two finest meshes yield similar results, and the time evolution of the segregation intensities demonstrates convergence of the solution with increasing refinement. Furthermore, the simulations appear to be converging on a solution within the error bounds of that derived from experiments. This comparison confirms the earlier observation that the segregation takes place over a similar time scale in simulations and experiments, even when the avalanche is severely under resolved. This means that the numerical simulations provide a strong qualitative and quantitative match with experimental data without the need for any fitting parameters.

Figure 7. Comparison of the segregation intensity $\mathcal {S}$, defined in (4.3), with time, in a 70 % filled drum with a 50:50 mix of large and small particles. The experimental data are computed every third of a revolution (i.e. every 8 s) and is shown using red circles. The error bars represent the $\pm$1 standard deviation in the RGB intensity mapping used to determine the small-particle concentration. The simulation data (coloured lines) is shown for grid resolutions $N=300$, 400, 500 and 600. This is equivalent to a total number $N^2$ mesh cells. In the inset, the segregation intensity of the simulations at $t=48$ s is plotted against the grid resolution, with the segregation intensity calculated for the experiment represented by the dashed blue line for comparison, while the dot-dashed red lines represent the experimental error of the intensity.

The simulation data shows periodic peaks and troughs in segregation intensity with a period of one third of a revolution. This can be expected given the variation in the avalanche dynamics with the shifting geometry (see § 4.3 and figure 3), which produces higher strain rates, longer particle trajectories through the avalanche and, hence, stronger segregation when the avalanche is at its longest. Both the experimental and computed segregation intensity begin to stabilize after the second revolution, with the final segregation intensity of $S=0.77$ at $t=48$ s at the highest resolution. After this time it will continue to plateau over many revolutions, but will never become fully periodic (Mounty Reference Mounty2007).