3.1. Density segregation theory for ternary mixtures

Consider a ternary granular mixture consisting of three different species. The particles of the three species have identical volumes ($V_p$) due to the same size of particles, and differ only in density (figure 2a). The masses of the heavy (green), medium (red) and light (blue) particles are $m_H$, $m_M$ and $m_L$, respectively, and $\rho _{i} = m_{i}/V_p$ is the density of the particles. Following Tripathi & Khakhar (Reference Tripathi and Khakhar2013), we use an effective medium approach and calculate the buoyant force on any particle by assuming it to be immersed in a continuum of density $\rho _m$ given as

(3.1)\begin{equation} \rho_{m}=\phi(\rho_{H} f_{H} + \rho_{M} f_{M} + \rho_{L} f_{L} ), \end{equation}

where $\phi$ is the total solids volume fraction, and $f_{i}=n_i/n$ is the number fraction of the $i{\rm th}$ species, with $n_{i}$ the number density of the $i{\rm th}$ species, and $n$ the total number density. The species volume fraction $\phi _i$ is equal to the product of the species number fraction and the particle volume. Since the sizes of the particles of all the species are the same, $\phi _i=n_iV_{p}$ and the total volume fraction $\phi =\phi _H+\phi _M+\phi _L$ equals the product of the total number density $n=n_H+n_M+n_L$ and particle volume $V_p$, i.e. $\phi =nV_p$. Hence the volume fraction of the $i$th species $\phi _i/\phi$ equals the number fraction $n_i/n=f_i$. Thus $f_i$ can also be interpreted as the concentration of the $i{\rm th}$ species per unit granular volume.

Figure 2. (a) Schematic of a ternary granular mixture containing particles of three different masses. High-density particles of mass $m_{H}$ are represented using green spheres, medium-density particles of mass $m_{M}$ are represented using red spheres, and low-density particles of mass $m_{L}$ are represented using blue spheres. (b) Gravity and buoyant force on the three types of particles in an effective medium of density $\rho _m$, where $g_y=g\cos \theta$ is the magnitude of the gravity component in the $y$-direction.

Consider a single heavy particle (green particle shown in figure 2b) immersed in this effective medium of density $\rho _m$. The net body force acting on this particle in the upward direction (i.e. in the positive $y$-direction) is

(3.2)\begin{equation} F_{H} = F_{B}-m_{H}g\cos\theta, \end{equation}

where $g \cos \theta$ is the magnitude of the $y$-component of gravitational acceleration, and $F_{B}$ is the buoyancy force exerted on this particle in the upward (i.e. positive $y$) direction by the granular medium (see figure 2b). To account for the presence of voids in the vicinity of this heavy particle, we hypothesise that the volume of the displaced granular fluid (shown using the dashed circle in figure 2b) is larger than the particle volume, so that $F_{B}={\rho }_{m}{V_{E}}{g\cos \theta }$, where ${V_{E}}$ is the volume of the displaced medium. The net force on the heavy particle in the upward $y$-direction is

(3.3)\begin{equation} F_{H}={\phi (\rho_{H} f_{H} + \rho_{M} f_{M} + \rho_{L} f_{L} )}{V_{E}}{g\cos\theta}-m_{H}g\cos\theta. \end{equation}

For the limiting case when the densities (and masses) of all three species are equal (i.e. $\rho _{H} = \rho _{M} = \rho _{L} = m_H/V_p$), there should be no net force on the particle, i.e. $F_{H} = 0$. Since $f_H+f_M+f_L=1$, we get the effective volume of the particles to be $V_{E}= V_p/\phi$. An identical expression for $V_E$ has been obtained in the case of the binary mixture by Tripathi & Khakhar (Reference Tripathi and Khakhar2013) and is consistent with the theoretical calculations of Kumar, Khakhar & Tripathi (Reference Kumar, Khakhar and Tripathi2019). This expression for $V_E$ indicates that the free void volume among the particles is shared equally by all the particles in the case of equal-size particles. Using this expression for $V_E$, we get the buoyancy force acting on the particles in the upward direction as

(3.4)\begin{equation} F_{B} = [f_{H} m_{H} + f_{M} m_{M} + f_{L} m_{L}]g\cos\theta, \end{equation}

and the net force on a heavy particle in the $y$-direction turns out to be

(3.5)\begin{equation} F_{H} ={-}m_{H} g\cos\theta + [f_{H} m_{H} + f_{M} m_{M} + f_{L} m_{L}]g\cos\theta. \end{equation}

Using $f_{H}=1- f_{M} - f_{L}$, the above expression simplifies to

(3.6)\begin{equation} F_{H} ={-}[(m_{H} - m_{M})f_{M} + (m_{H} - m_{L}) f_{L}]g\cos\theta. \end{equation}

The forces on the medium- and light-mass particles in the $y$-direction are obtained in a similar fashion, giving

(3.7)$$\begin{gather} F_{M} ={-}[(m_{M} - m_{H})f_{H} + (m_{M} - m_{L}) f_{L}]g\cos\theta, \end{gather}$$

(3.8)$$\begin{gather}F_{L} ={-}[(m_{L} - m_{M})f_{M} + (m_{L} - m_{H}) f_{H}]g\cos\theta. \end{gather}$$

Since $m_H>m_M>m_L$, (3.6) shows that the value of the net force on the heavy particle in the upward (i.e. positive $y$) direction is always negative. Similarly, the net upward force on the light particle always remains positive. The force direction on the medium-mass particle, however, depends upon the relative concentrations of heavy and light particles, and the mass differences between the particles. The net force acting on the heavy (light) particle leads to a downward (upward) motion of the particle that is balanced by the viscous drag offered by the granular layer.

Force balance equations in the flowing direction for the three types of particles can also be written in a manner similar to (3.6), (3.7) and (3.8). Such a force balance in the $x$-direction gives rise to a slip velocity of the heavy particle in the flow direction due to the resultant force of gravity and buoyancy acting in the $x$-direction. Due to this slip velocity, the heavy particles move with a slightly different velocity in the $x$-direction compared to other particles. Such slip velocities have been reported for different-size tracer particles van der Vaart et al. (Reference van der Vaart, van Schrojenstein Lantman, Weinhart, Luding, Ancey and Thornton2018). Our simulation results indeed show a very small but finite slip velocity between the heavy, medium and light particle velocities in the flow direction. A detailed calculation of the lift force for same-size different-density tracer particles is also possible using the approach of van der Vaart et al. (Reference van der Vaart, van Schrojenstein Lantman, Weinhart, Luding, Ancey and Thornton2018). However, as discussed below, we find that the effect of this lift force on the segregation can be captured adequately without knowing the detailed expression.

Following Tripathi & Khakhar (Reference Tripathi and Khakhar2011a, Reference Tripathi and Khakhar2013), we assume that the viscous drag force on a settling sphere in a granular fluid is proportional to the viscosity $\eta$ and settling velocity $v_i$, i.e. $F_D=-c{\rm \pi} \eta v_i d$, where $c=c(\phi )$ is the solids volume fraction dependent Stokes drag coefficient. For a heavy particle that is settling at a constant velocity through a surrounding mixture of particles, $F_H+F_D=0$, and the settling heavy particle achieves a terminal velocity

(3.9)\begin{equation} v_{H} = \dfrac{F_{H}}{c{\rm \pi}\eta d}. \end{equation}

It is evident that the above equation does not account for the presence of the lift force due to the slip velocity of the heavy particles. Let $\beta$ be the ratio of the magnitude of the lift force to the combined force of gravity and buoyancy $F_H$, i.e. $F_{Lift}=\beta F_H$. In general, the factor $\beta$ is not a constant and needs to account for the unknown dependence of the shear rate, granular viscosity, particle density and all other factors that influence the ratio $F_{Lift}/F_H$. By accounting for the lift force in the force balance equation in the $y$-direction, we obtain $F_H+F_{Lift}=c{\rm \pi} \eta v_H d$. Using $F_{Lift}=\beta F_H$, we get $F_H=c{\rm \pi} \eta v_H d/(1+\beta )$. Thus the explicit inclusion of the lift force in the flow direction effectively leads to a modified value of the Stokes drag coefficient $c(\phi )=c/\beta$. The DEM simulations in Tripathi & Khakhar (Reference Tripathi and Khakhar2011a, Reference Tripathi and Khakhar2013) use the total $y$-direction force on the heavy particles to obtain the values of the Stokes drag coefficient $c(\phi )$ and thus implicitly account for both the lift force due to slip velocity between the species as well as the net force due to the combined effect of gravity and buoyancy in the $y$-direction.

Settling of these heavy particles leads to the segregation flux $\dot {J}^{S}_{H}$ of heavy particles

(3.10)\begin{equation} \dot{J}^{S}_{H} = n_{H}(v_{H} - v_{m}), \end{equation}

where $v_{m} = f_{H}v_{H} + f_{M} v_{M} + f_{L}v_{L}$ is the local volume average mixture velocity in the downward direction. Using the expressions for segregation velocities for medium- and light-mass particles as $v_{M} = F_{M}/c{\rm \pi} \eta d$ and $v_{L} = F_{L}/c{\rm \pi} \eta d$, where $F_M$ and $F_L$ are given by (3.7) and 3.8, we get $v_{m} = 0$. The segregation flux of heavy particles is given as

(3.11)\begin{equation} \dot{J}^{S}_{H} ={-}n f_{H}\,\dfrac{[(m_{H} - m_{M})f_{M} + (m_{H} - m_{L}) f_{L}]g\cos\theta}{c{\rm \pi}\eta d}. \end{equation}

The diffusion flux $\dot {J}^{D}_{H}$ of heavy particles in the $y$-direction is

(3.12)\begin{equation} \dot{J}^{D}_{H} ={-} n D_H\,\dfrac{{\rm d}f_{H}}{{\rm d} y}, \end{equation}

where $D_H$ is the self-diffusion coefficient of heavy particles. The constitutive equations (3.11) and (3.12), along with a transport equation of the form

(3.13)\begin{equation} \dfrac{\partial (n f_H)}{\partial t} ={-} \dfrac{\partial}{\partial y}(\dot{J}_{H}^{S} + \dot{J}_{H}^{D}), \end{equation}

can be used to obtain the time-dependent concentration profiles of the heavy species. We consider the steady-state case for the flow and segregation in the layer to get

(3.14)\begin{equation} \dot{J}_{H}^{S} + \dot{J}_{H}^{D} = 0. \end{equation}

Substituting (3.11) and (3.12) in the above equation and rearranging gives the following equation for the steady-state concentration field of the heavy species:

(3.15)\begin{equation} \dfrac{{\rm d}f_{H}}{{\rm d} y} ={-} \dfrac{f_{H}[( m_{H} - m_{M}) f_{M} + (m_{H} - m_{L}) f_{L} ] g \cos \theta}{D c{\rm \pi}\eta d}. \end{equation}

Various studies (Savage Reference Savage1993; Utter & Behringer Reference Utter and Behringer2004; Tripathi & Khakhar Reference Tripathi and Khakhar2013; Fan & Hill Reference Fan and Hill2015) have shown that the diffusion coefficient for granular mixtures in case of shear flows varies linearly with the shear rate $\dot {\gamma }$, i.e.

(3.16)\begin{equation} D_H = b \dot{\gamma} d^2, \end{equation}

where $b$ is a constant. Using the expression for viscosity $\eta = |\tau _{yx}|/\dot {\gamma }$ along with (3.16), we get

(3.17)\begin{equation} \dfrac{{\rm d}f_{H}}{{\rm d} y} ={-} \dfrac{f_{H}[( m_{H} - m_{M}) f_{M} + (m_{H} - m_{L}) f_{L} ] g \cos \theta}{b c {\rm \pi}\,|\tau_{yx}|\,{d}^3}. \end{equation}

Our DEM simulation results show that the self-diffusion coefficient for the three species is identical, i.e. $D_H=D_M=D_L$, so the constant $b$ remains the same for all three species. Using this result, the equations governing the concentrations of medium- and light-mass particles are obtained as

(3.18)$$\begin{gather} \dfrac{{\rm d}f_{M}}{{\rm d} y} ={-} \dfrac{f_{M}[( m_{M} - m_{H}) f_{H} + (m_{M} - m_{L}) f_{L} ]g \cos \theta}{b c {\rm \pi}\,|\tau_{yx}|\,{d}^3}, \end{gather}$$

(3.19)$$\begin{gather}\dfrac{{\rm d}f_{L}}{{\rm d} y} ={-} \dfrac{f_{L}[( m_{L} - m_{M}) f_{M} + (m_{L} - m_{H}) f_{H} ]g \cos \theta}{b c {\rm \pi}\,|\tau_{yx}|\,{d}^3}. \end{gather}$$

Equations (3.17)–(3.19) show that the concentration of any one species depends upon the concentrations of the other two species and the mass (or density) difference between the species in addition to the gravity component and shear stress. As mentioned before, the Stokes drag coefficient $c$ depends upon the packing fraction $\phi$, and the values reported by Tripathi & Khakhar (Reference Tripathi and Khakhar2013) are used in the present study.

3.2. Extension to $N$-component mixtures

We now generalise the ternary density segregation theory to mixtures having $N$ different components. Following the approach mentioned in the previous subsection, the net force acting on a particle of the $i{\rm th}$ component in the $y$-direction is

(3.20)\begin{equation} F_{i} = \sum_{j=1}^{N}m_{j} f_{j} g\cos\theta - m_{i} g\cos\theta, \end{equation}

where $f_i$ is the concentration of the $i{\rm th}$ species per unit granular volume. Making use of the identity $\sum _{j=1}^{N}f_{j} = 1$, the above expression simplifies to

(3.21)\begin{equation} F_{i} ={-}\sum_{j=1}^{N} (m_{i} - m_{j})f_{j} g\cos\theta. \end{equation}

Equations (3.20) and (3.21) for the $N$-component mixtures are analogous to (3.5) and (3.6) of the ternary mixture. As before, the steady-state terminal velocity of the $i{\rm th}$ species particle is obtained by balancing the net body force with the drag force:

(3.22)\begin{equation} v_{i} = \dfrac{F_{i}}{c{\rm \pi}\eta d}. \end{equation}

Using these expressions for $v_{i}$, we obtain $v_{m} = 0$, and the segregation flux of the $i{\rm th}$ species is calculated as $\dot {J}^{S}_{i} = n f_{i} v_{i}$, giving

(3.23)\begin{equation} \dot{J}^{S}_{i} ={-}\dfrac{n f_{i}}{c{\rm \pi}\eta d}{\sum_{j=1}^{N} (m_{i} - m_{j}) f_{j} g\cos\theta}. \end{equation}

The diffusion flux $\dot {J}^{D}_{i}$ of the $i{\rm th}$ species is given as

(3.24)\begin{equation} \dot{J}^{D}_{i} ={-} n D\,\dfrac{{\rm d} f_{i}}{{\rm d} y}. \end{equation}

The time-dependent concentration $f_i$ of species $i$ can be obtained using the transport equation

(3.25)\begin{equation} \dfrac{\partial (n f_i)}{\partial t} ={-}\dfrac{\partial}{\partial y}(\dot{J}_{i}^{S} + \dot{J}_{i}^{D}). \end{equation}

For an incompressible flow with constant $\phi \ (=nV_p)$, we get an equation similar to the one given by Gray & Ancey (Reference Gray and Ancey2011) for multi-component size segregation:

(3.26)\begin{equation} \dfrac{\partial (f_i)}{\partial t} -\frac{\partial}{\partial y}\left(\dfrac{ f_{i}}{c{\rm \pi}\eta d}{\sum_{j=1}^{N} (m_{i} - m_{j}) f_{j} g\cos\theta}\right)= \frac{\partial}{\partial y}\left( D\,\dfrac{{\rm d} f_{i}}{{\rm d} y}\right). \end{equation}

Interestingly, this equation has the same species concentration dependence of the segregation flux as postulated in (2.28) of Gray & Ancey (Reference Gray and Ancey2011) for multi-component mixtures of different size particles. Equation (3.26) shows that the segregation coefficients (or segregation parameters) for multi-component mixtures of different density particles depend on the mass differences between the particles, size of the particles, drag coefficient, component of the gravitational acceleration in the segregation direction, and local viscosity. Using (3.26), the steady-state concentration profile $f_i$ of any species $i$ is given as

(3.27)\begin{equation} \dfrac{{\rm d} f_{i}}{{\rm d} y} ={-}\dfrac{f_{i}}{D c{\rm \pi}\eta d} \sum_{j=1}^{N} (m_{i} - m_{j}) f_{j} g\cos\theta. \end{equation}

Using the fact that the self-diffusion coefficients of different species are identical and are given by $D = b \dot {\gamma } d^2$ for shear flows, the species concentration is described using the equation

(3.28)\begin{equation} \dfrac{{\rm d} f_{i}}{{\rm d} y} ={-}\dfrac{ f_{i} g \cos{\theta} }{b c {\rm \pi}\, |\tau_{yx}|\,d^3} \sum_{j=1}^{N} (m_{i} - m_{j})f_{j}. \end{equation}

The above equation represents the generalised form of (3.17)–(3.19) of the ternary mixture.

3.3. Multi-component mixture rheology and momentum balance equations

Equation (3.27) clearly shows that the local species concentration depends upon the local viscosity of the granular medium. This presence of viscosity in the density segregation model leads to the coupling of segregation with the rheology. The dense flow rheology of granular systems is described by the empirical dependence of the shear stress to pressure ratio (referred to as the effective friction coefficient $\mu$) in terms of a non-dimensional inertial number $I$, i.e. ${|\tau _{yx}|}/{P}=\mu (I)$ (GDR MiDi 2004; Da Cruz et al. Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005). Jop, Forterre & Pouliquen (Reference Jop, Forterre and Pouliquen2006) suggested the following empirical form to capture the dependence of $\mu (I)$ on $I$:

(3.29)\begin{equation} \mu(I)=\mu_{s}+\frac{\mu_{m}-\mu_{s}}{1+I_{0}/I}, \end{equation}

where $\mu _{s}$, $\mu _{m}$ and $I_{0}$ are the rheological model parameters that need to be obtained from simulations/experiments, and the inertial number $I=\dot {\gamma } d/(P/\rho )^{1/2}$ depends upon the local shear rate $\dot {\gamma }$ and pressure $P$, in addition to the particle size $d$ and density $\rho$. Concerns about this rheology have been raised by various researchers (Barker et al. Reference Barker, Schaeffer, Bohórquez and Gray2015, Reference Barker, Schaeffer, Shearer and Gray2017; Heyman et al. Reference Heyman, Delannay, Tabuteau and Valance2017; Dsouza & Nott Reference Dsouza and Nott2020) since the $\mu \unicode{x2013}I$ rheology is found to be mathematically ill-posed at high and low inertial numbers. This ill-posedness of the rheological equation leads to grid-size-dependent results in the case of two-dimensional transient computations. However, these issues do not pose any problem for the steady-state flow property calculations in the present work. This rheological model is complemented with a dilatancy law describing the dependence of the solids volume fraction $\phi$ on the inertial number $I$. Most studies suggest a linear variation of $\phi$ with $I$:

(3.30)\begin{equation} \phi=\phi_{max,l}-(\phi_{max,l}-\phi_{min,l})I, \end{equation}

where $\phi _{max,l}$ and $\phi _{min,l}$ are model parameters.

Tripathi & Khakhar (Reference Tripathi and Khakhar2011b) have shown that the viscoplastic rheology of nearly monodisperse granular systems in the dense flow regime can be extended to binary granular mixtures in terms of a generalised inertial number $I_{mix}=\dot {\gamma } d_{mix}/(P/\rho _{mix})^{1/2}$, where $d_{mix}$ and $\rho _{mix}$ are the local volume-averaged size and density of the particles. In the case of binary mixtures having same-size and different-density particles, $d_{mix}=d$ and $\rho _{mix}=\rho _1 f_1 + \rho _2 f_2$. Using this generalised inertial number rheology, the authors were able to predict the coupled effect of rheology and segregation in steady fully developed chute flows of binary mixtures differing in density for a range of density ratios (Tripathi & Khakhar Reference Tripathi and Khakhar2013). The authors also showed the presence of normal stress differences and accounted for this effect by a normal stress difference parameter $a$ that is used to correct the pressure calculated from momentum balance equations.

In this study, we extend the binary mixture rheology to multi-component mixtures by accounting for the presence of all the species in the expression of local volume average density. The local average density of particles in the case of a multi-component mixture is calculated as

(3.31)\begin{equation} \rho_{mix}(y)=\sum_{i=1}^{N} \rho_i f_i, \end{equation}

which requires knowledge of the concentrations of all the species $f_i$. The calculation of the species concentration $f_i$ (using (3.28)) requires knowledge of the shear stress $\tau _{yx}$. This calculation of the shear stress (using (3.32)) requires knowledge of $\rho _{mix}$, which in turn depends upon the species concentration $f_i$ (see (3.31)). This interdependence of the concentration and the shear stress on each other is indicative of the inter-coupling of rheology and segregation with each other. In order to be able to predict the species concentrations due to the segregation of different density particles in steady fully developed chute flow, the momentum balance equations need to be solved along with the rheological model and density segregation model. The values of the rheological parameters are taken from Tripathi & Khakhar (Reference Tripathi and Khakhar2011b) and are reported in table 1.

Table 1. Values of the rheological model parameters (as in Tripathi & Khakhar Reference Tripathi and Khakhar2011b) used in this study.

Solving momentum balance equations for an incompressible, steady and fully developed flow over an inclined plane at inclination angle $\theta$ in the $x$- and $y$-directions leads to

(3.32)$$\begin{gather} \tau_{yx}(y) ={-}g \sin \theta \int_{y}^{H} \phi \rho_{mix}\, {{\rm d} y}, \end{gather}$$

(3.33)$$\begin{gather}P (y) = g \cos \theta (1 - a \tan \theta )\int_{y}^{H} \phi \rho_{mix} \,{{\rm d} y}, \end{gather}$$

where $H$ is the thickness of the flowing granular layer (known from the DEM simulations), and the parameter $a$ accounts for the anisotropy of the stress tensor due to non-zero normal stress differences (Tripathi & Khakhar Reference Tripathi and Khakhar2011b). Calculation of $\tau _{yx}$ and $P$ requires knowledge of $\phi$ and $\rho _{mix}$. Following Tripathi & Khakhar (Reference Tripathi and Khakhar2011b), the solids fraction $\phi$ is obtained by using the expression for the generalised inertial number $I_{mix}$ in the dilatancy law

(3.34)\begin{equation} \phi=\phi_{max,l}-(\phi_{max,l}-\phi_{min,l})I_{mix}. \end{equation}

The iterative procedure to solve the momentum balance equations along with the mixture rheological model and density segregation model is described in the next subsection.

3.4. Numerical solution of the equations

Using (3.32) and (3.33), the value of the effective friction coefficient at any inclination $\theta$ is calculated as

(3.35)\begin{equation} \frac{|\tau_{yx}|}{P}=\mu(I) = \dfrac{\tan\theta} { 1 - a \tan\theta}. \end{equation}

The inertial number $I$ is calculated by rearranging the empirical form of $\mu (I)$ as

(3.36)\begin{equation} I_{mix} =I_{o}\,\dfrac{\mu(I)-\mu_{s}}{\mu_{m}-\mu(I)}, \end{equation}

where we use the value of $\mu (I)$ obtained in terms of $\theta$ from the previous equation. The calculation of the shear rate $\dot {\gamma }$ is done from the previously calculated value of the inertial number $I_{mix}$ using the relation

(3.37)\begin{equation} \dot{\gamma}(y)=\dfrac{I_{mix}}{d}\sqrt{\dfrac{P}{\rho_{mix}}}. \end{equation}

Finally, the velocity profile $v_{x}(y)$ is obtained by integrating the shear rate as

(3.38)\begin{equation} v_x(y)= v_{slip} + \int_{0}^{y}\dot{\gamma} \,{{\rm d} y}, \end{equation}

where $v_{slip}$ is the slip velocity at the chute base. We use a rough and bumpy base so that $v_{slip} \approx 0$ for all the cases considered in this study. However, the prediction of the shear rate using (3.37) requires knowledge of the local mixture density $\rho _{mix}$ as well as the local pressure $P$. Since the calculation of these properties requires knowledge of species concentration $f_i$, we utilise the multi-component segregation model proposed before as follows.

Algorithm 1: Algorithm used for predicting properties of an N-component mixture

The segregation model presented before requires that the system of equations represented by (3.28) be solved simultaneously for all $N$ components of the mixture. These equations are solved iteratively for $N-1$ components in MATLAB using the ode45 solver, while the $N{\rm th}$ species concentration is calculated as $f_N = 1 - \sum _{i=1}^{N-1}f_i$. To solve these ordinary differential equations, the equations are expressed in symbolic form. While all other terms in (3.28) are in symbolic form, the shear stress $\tau _{yx}$ is calculated numerically using (3.31) and (3.32). Thus, given a set of values of $f_i$ at each $y$, we obtain values of $\tau _{yx}$ at each $y$. To represent the $\tau _{yx}(y)$ in symbolic form, we fit a polynomial using the poly-fit function in MATLAB, and obtain $\tau _{yx}(y)$ as a function of $y$. Note that this curve fitting is used only to represent the calculated shear stress data as a polynomial, hence any curve that is able to represent the data well suffices. Our simulation results show that while a linear fit seems to work well for density ratios up to 1.5, it fails for higher density ratios. Similarly, a quadratic fit also seems to work for most of the cases but fails in the case of mixtures with density ratios larger than 3. Given the monotonic variation of the shear stress, a cubic polynomial fit does not capture $\tau _{yx}$ variation with $y$. A quartic polynomial fit for $\tau _{yx}$ with $y$ works for all the cases, hence a fourth-order polynomial is fitted to $\tau _{yx}$.

The boundary conditions are required to solve the system of equations represented by (3.28). We use a guess value $f_{i0} = f_{i}(y=0)$, where $f_{i0}$ is the concentration of the $i{\rm th}$ species at the base $y=0$. Note that the species concentrations at the base are not known a priori and remain to be determined. To determine these $f_{i0}$ values, we use the species mass balance equation for each species. For the appropriate value of $f_{i0}$ to be used as a boundary condition, the total species concentration of the layer must be equal to $f_{i}^{T}$, i.e. $({1}/{H}) \int _{0}^{H} f_{i} \,{{\rm d} y} = f_{i}^{T}$, where $f_{i}^{T}$ is the total concentration of the $i{\rm th}$ species. This species balance condition is used as the convergence criterion for calculating the unknown $f_{i0}$ value at the base $y=0$. We use the following condition for convergence of $f_i$ profiles:

(3.39)\begin{equation} \sqrt{\sum_{i=1}^{N} \left|\dfrac{1}{H} \int_{0}^{H} f_{i}\, {{\rm d} y} - f_{i}^{T}\right|^2 } < 10^{{-}4}. \end{equation}

The guess values $f_{i0}$ are updated iteratively as follows:

(3.40)\begin{equation} f_{i0} = \begin{cases} f_{i0} + \Delta f_0, & \text{if} \ f_{i}^{T} > \dfrac{1}{H} \displaystyle\int_{0}^{H} f_{i}\, {{\rm d} y}, \\[10pt] f_{i0} - \Delta f_0, & \text{if} \ f_{i}^{T} < \dfrac{1}{H} \displaystyle\int_{0}^{H} f_{i}\, {{\rm d} y}, \end{cases} \end{equation}

until the convergence criterion given by (3.39) is satisfied. We choose $\Delta f_0$ as $10^{-5}$ in this study, and report the algorithm for species concentration and other mixture properties in Algorithm 1.

4.1. Effect of density ratios

Figure 6 shows the segregation of mixtures with different density ratios. The simulations are performed at inclination angle $\theta =25^\circ$ with equal compositions of each species, i.e. $f_L^T=f_M^T=f_H^T=1/3$. Figure 6(a) shows the snapshot for density ratio $\rho _H:\rho _M:\rho _L=1.50:1.25:1.00$; figure 6(b) shows the snapshot for density ratio $\rho _H:\rho _M:\rho _L=3:2:1$; figure 6(c) shows the snapshot for density ratio $\rho _H:\rho _M:\rho _L=5:3:1$; and figures 6(d)–6( f) show the concentration profiles for the respective cases below the snapshots.

Figure 6. Simulation snapshots and comparison of the theoretical predictions (lines) with simulation results (symbols) for ternary mixtures for different density ratios with equal composition (33.33 %) of each species in the ternary mixture at inclination angle $\theta = 25^\circ$.

For small differences in particle densities, figure 6(a) shows that the region near the free surface resembles a binary mixture with a few medium-density (red) particles dispersed in light (blue) particles. This is also evident from figure 6(d), where $f_L>f_M$ and $f_H\rightarrow 0$ near the free surface. The rest of the layer consists of a significant fraction of all three species; the concentration of heavy (green) particles increases, while that of light particles decreases towards the base. The concentration of medium-density particles increases from ${\approx }10\,\%$ near the free surface to ${\approx }40\,\%$ at $y=20$, and then shows little variation in the remaining part of the layer.

While the light particles do not undergo complete segregation for small density ratios, figures 6(b), 6(e) and 6(c), 6( f) show that a pure layer of light particles near the free surface is observable for higher-density ratios with $f_M=f_H=0$. A binary mixture of medium and heavy particles consisting of larger fractions of heavy and smaller fractions of medium-density particles is found in the lower part of the layer, with $f_L=0$ in both cases. This is in contrast to the results shown for the case of a low density ratio where the lower portion of the layer consisted of all three species. The middle portion of the layer consists of a ternary mixture rich in medium-density particles. The concentration of light particles keeps increasing, and that of heavy particles keeps decreasing with increasing distance from the base. The medium-density particle concentration shows non-monotonic behaviour with $y$ and shows the maximum concentration in the middle of the layer. The theoretical predictions of the concentrations of all the species, shown using solid lines, are found to be in excellent agreement with the theory for all the density ratios considered in this study. A comparison of figures 6(d), 6(e) and 6( f) indicates that the extent of segregation in the layer increases with the increasing density difference among the species.

Figures 6(g) and 6(h) show that the pressure and shear stress profiles vary linearly with $y$ near the free surface. This is attributed to the high concentration of light particles near the free surface, which leads to nearly constant bulk density of the layer near the free surface for all three cases. However, the nonlinear behaviour of the shear stress and pressure is exhibited in the lower part of the layer. As discussed before, this nonlinear behaviour is attributed to the increasing bulk density of the layer towards the base due to the increasing concentration of heavy particles near the base. With increasing density ratios, the magnitudes of shear stress and pressure increase due to increased bulk density of the layer. Figure 6(i) shows that an increase in the density ratios leads to a slight decrease in the layer velocity. These interesting variations of the shear stress, pressure and velocity are very well captured by the theory.

4.2. Influence of species total compositions

Figure 7 shows the results for three different mixture compositions of heavy, medium and light particles. Symbols represent the simulation data for density ratio $\rho _L:\rho _M:\rho _H = 1:2:3$ at inclination angle $\theta = 25^\circ$. The total concentration of one of the species in the mixture is $50\,\%$, while the concentration of the other two species is $25\,\%$. Figures 7(a), 7(b) and 7(c) show steady-state snapshots of the system corresponding to $50\,\%$ light (blue), $50\,\%$ medium-density (red) and $50\,\%$ heavy (green) particles, respectively. Figures 7(d)–7( f) show concentration profiles for the respective mixtures below the snapshots. In all the cases, a pure layer of light (blue) particles is observable near the free surface. The thickness of the pure light particles layer is around one-third of the layer thickness in case $f_L^T=50\,\%$ and gets reduced in proportion of the total composition in figures 7(e) and 7( f). Beneath this pure layer of light particles, the concentration of light particles keeps reducing, while that of heavy particles keeps increasing as one moves towards the base.

Figure 7. Simulation snapshots and comparison of the theoretical predictions (lines) with simulation results (symbols) for ternary mixtures for various compositions of each species with density ratios $\rho _{H} : \rho _{M}:\rho _{L} = 3:2:1$ at inclination angle $\theta = 25^\circ$.

The concentration profiles of heavy particles seem to be similar in mixtures having identical total heavy particle composition as shown in figures 7(d) and 7(e). Similarly, the light particle concentrations also appear to be very similar in mixtures with identical total light particle composition (figures 7e, f). The concentrations of medium-density particles, however, differ significantly from each other in figures 7(d) and 7( f), despite having identical total composition. This shows a non-monotonic variation with $y$, and displays a maximum concentration somewhere in the middle part of the layer. As before, the theoretical predictions shown using solid lines capture all these features very well, including the nonlinear variation of shear stress (figure 7g) and pressure (figure 7h). Figure 7(i) illustrates that the mixture velocity profiles show very little effect of the mixture composition, which is also reflected in the theoretical predictions.

Figures 8(a)–8(c) show the concentration variations for mixtures having high concentration of two species with very little concentration of the third species. Specifically, figure 8(a) shows the concentration variation for a mixture with $50\,\%$ composition of the light species, $40\,\%$ composition of the medium species, and $10\,\%$ composition of the heavy species. Figure 8(b) shows the concentration variation for a mixture with $10\,\%$ light, $50\,\%$ medium and $40\,\%$ heavy species. Similarly, figure 8(c) shows the concentration variation for a mixture with $40\,\%$ light, $10\,\%$ medium and $50\,\%$ heavy species. Figures 8(d)–8( f) show results for three different mixtures having species concentration combinations of $70\,\%$, $20\,\%$ and $10\,\%$. As before, the symbols represent the simulation data for density ratio ($\rho _L:\rho _M:\rho _H = 1:2:3)$ at inclination angle $\theta = 25^\circ$, and lines represent the theoretical predictions. The concentration profiles of heavy, medium and light particles for all these cases are captured very well by the theory. The variations of shear stress, pressure and velocity (not shown) are similar to those reported in figure 7. The theory captures accurately these mixture properties for all these cases as well. These results establish that the theory presented in this work is able to predict accurately various properties of interest for a wide range of compositions of ternary mixtures.

Figure 8. Comparison of the theoretical predictions (lines) with simulation results (symbols) for ternary mixture for various compositions of each species with the same density ratios $\rho _{H} : \rho _{M}:\rho _{L} = 3:2:1$ at inclination angle $\theta = 25^\circ$.