2.1. Force balance on the large intruder

The configuration is sketched in figure 1. The large particle is of diameter $d_l$, volume $V_l$ and density $\rho ^p$ in a bed of height $h$ made of small particles. Below this layer the grains are in the quasi-static regime. Applying Newton's second law, the vertical Lagrangian equation of the intruder can be expressed as

(2.1)\begin{equation} \rho^p V_l\dfrac{\textrm{d} w^l}{\textrm{d}t} = P - \varPi_f + f^f_{d} + f^p_{d} - \varPi_p - f_l. \end{equation}

In (2.1), the large intruder is submitted to six forces (see figure 1): its weight $P = - \rho ^p V_l g \cos \theta$, the fluid buoyancy force $\varPi _f = -\rho ^f V_l g \cos \theta$, the granular buoyancy $\varPi_p$ due to the reaction of the small particles, the drag force exerted by the fluid $f_d^f$, the drag force exerted by the small particles $f_d^p$ and the lift force $f_l$ responsible for the ascent of the large intruder (Guillard et al. Reference Guillard, Forterre and Pouliquen2014, Reference Guillard, Forterre and Pouliquen2016). In the present configuration, the slope angle is low ($\tan \theta = 0.1$) and the streamwise gravity component is negligible, making the contribution of the shear stress gradient to the granular buoyancy and lift forces negligible (see (1.2) from Guillard et al. Reference Guillard, Forterre and Pouliquen2016). Thus, the pressure gradient contribution of (1.2) is dominant and the granular buoyancy force is defined as

(2.2)\begin{equation} \varPi_p = \dfrac{V_l}{\varPhi_{max}}\displaystyle\frac{\partial p^s}{\partial z}, \end{equation}

where $p^s$ is the overburden pressure of the surrounding small particles. The lift force is

(2.3)\begin{equation} f_l = V_l \mathcal{F}_{l}(\mu)\displaystyle\frac{\partial p^s}{\partial z}, \end{equation}

where $\mathcal {F}_{l}(\mu ) = (1-\exp ({-70(\mu -\mu _c)}))$ will be called the empirical segregation function and is proposed based on the empirical function $\mathcal {F}(\mu )$ (see (1.3a,b)) in order to only account for the lift force part of (1.2) (see Appendix A).

Figure 1.

Vertical component of the forces acting on a large intruder; $\boldsymbol {\varPi }_{\boldsymbol {f}}$ (blue solid line) is the buoyancy due to the fluid, $\boldsymbol {\varPi }_{\boldsymbol {p}}$ (——, black solid line) is the granular buoyancy and $\boldsymbol {f}_{\boldsymbol {l}}$ (——, red solid line) is the lift force identified by Guillard et al. (Reference Guillard, Forterre and Pouliquen2014, Reference Guillard, Forterre and Pouliquen2016). The particle is also submitted to the drag forces $\boldsymbol {f}^{\boldsymbol {p}}_{\boldsymbol {d}}$ (——, light green solid line) and $\boldsymbol {f}^{\boldsymbol {f}}_{\boldsymbol {d}}$ (blue solid line), respectively due to the interaction with small particles (Tripathi & Khakhar Reference Tripathi and Khakhar2013) and the fluid.

While segregating at a velocity $w^l$, the large particle is submitted to a fluid drag force. Because the particulate Reynolds number based on the vertical velocity $Re_p = d_l \rho ^f w^l/ \eta ^f$ is very small in the bed, fluid inertial effects are negligible at the particle length scale and the vertical fluid drag force $f^{d}_f$ may be approximated by the Stokes law (Stokes Reference Stokes1851)

(2.4)\begin{equation} f^f_{d} = 3 {\rm \pi}\eta^f d_l (w^f - w^l), \end{equation}

where $\eta ^f$ is the fluid dynamic viscosity and $w^f$ is the vertical velocity of the fluid.

During its segregation motion, the intruder is also submitted to frictional forces from the surrounding small particles. This results in a particle drag force $f^p_{d}$ modelled as proposed by Tripathi & Khakhar (Reference Tripathi and Khakhar2013) (see (1.4)) as

(2.5)\begin{equation} f^{p}_{d} = c {\rm \pi}\eta^{p} d_l \left(w^{s}-w^l\right), \end{equation}

where $w^s$ is the vertical velocity of small particles and the drag coefficient $c$ is first approximated as a constant equal to 3 (Tripathi & Khakhar Reference Tripathi and Khakhar2013).

2.2. Dimensionless equation for the large intruder

In order to identify the dominant terms in (2.1), it is made dimensionless using classical scalings for granular flows

(2.6a–d)\begin{equation} w^k = \sqrt{d_l g} \tilde{w}^k, \quad z = d_l \tilde{z}, \quad t = \sqrt{d_l / g} \tilde{t} \quad \text{and} \quad p^k = \rho^p d_l g \widetilde{p^k}, \end{equation}

where $k=s$, $l$ or $f$ respectively for the surrounding small particles, the large intruder and the fluid. Introducing these variables and (2.3) in (2.1) and taking into account that $\cos \theta \sim 1$, the dimensionless form of the large intruder Lagrangian equation of motion can be written as

(2.7)\begin{equation} \dfrac{\textrm{d} \widetilde{w^l}}{\textrm{d}\tilde{t}} = \dfrac{\widetilde{w^f}-\widetilde{w^l}}{St^f} + \dfrac{\widetilde{w^s}-\widetilde{w^l}}{St^p} - \mathcal{F}_{l}(\mu) \displaystyle\frac{\partial\widetilde{p^s}}{\partial\tilde{z}}. \end{equation}

Equation (2.7) contains two dimensionless numbers. The first one is the fluid Stokes number

(2.8)\begin{equation} St^f = \dfrac{\rho^p d_l W}{6 c \eta^f}, \end{equation}

in which $W=\sqrt {d_l g}$ is the characteristic velocity of the large particle. This Stokes number compares the inertia of the large intruder with the viscous friction exerted by the fluid. Similarly, the granular Stokes number

(2.9)\begin{equation} St^p = \dfrac{\rho^p d_l W}{6 c \eta^p}, \end{equation}

compares the inertia of the large intruder with the contact friction exerted by the small particles in the vicinity of the intruder.

Assuming a classical bedload configuration, the fluid flows inside the porous matrix of the granular bed. Only the first layer of particles at the top is in a dense flow regime. Below this layer, $\mu <\mu _c$, meaning that the grains are in the quasi-static regime. Typical values of the granular viscosity for dense granular flows are very high compared with the water viscosity (typical ranges span from $10^3\ \textrm {Pa}\ \textrm {s}$ at the bed surface to $10^6\ \textrm {Pa}\ \textrm {s}$ at the bed bottom). This results in a fluid Stokes number $St^f$ much larger than the granular Stokes number $St^p$ whatever the height into the bed. Therefore, the first term on the right-hand side of (2.7), representing the fluid drag force, can be neglected. In addition, while the large intruder is rising, it only modifies the small particle bed structure locally. Therefore, it is assumed that the vertical velocity of the intruder does not disturb the vertical bulk velocity, implying that $w^s \sim 0$ (Tripathi & Khakhar (Reference Tripathi and Khakhar2011) showed that the bulk flow around the intruder was disturbed until a distance of 1 diameter of the large intruder for higher shear stress). Focusing on the position of the intruder in the quasi-static part of the bed, it can be deduced from (2.7) that the total solid volume fraction is constant, $\varPhi = cste$. Therefore, for this configuration, a simple equation for the vertical velocity of the large intruder can be written as

(2.10)\begin{equation} \dfrac{\textrm{d} \widetilde{w^l}}{\textrm{d}\tilde{t}} + \dfrac{1}{St^p} \widetilde{w^l} = \dfrac{\rho^p - \rho^f}{\rho^p} \varPhi \mathcal{F}_{l}(\mu). \end{equation}

In this dimensionless equation the fluid only acts through the fluid density coming from the granular pressure gradient. Therefore, the fluid can be considered as inert and (2.10) should be valid to model the vertical velocity of an intruder segregating in a dry granular flow (in this case the granular pressure gradient does not include the fluid density).

Equation (2.10) allows one to identify the main size segregation mechanisms and shows that the segregation of a large intruder can be seen as a simple relaxation process with characteristic time $St^p$.

3.1. Three-dimensional general governing equations

Following Jackson (Reference Jackson1997, Reference Jackson2000), the mass and momentum balance equations for each class are given by

(3.1)\begin{gather} \displaystyle\frac{\partial\epsilon \rho^f \boldsymbol{u}^f}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\epsilon \rho^f \boldsymbol{u}^f\right), \end{gather}

(3.2)\begin{gather}\displaystyle\frac{\partial\varPhi^i \rho^p \boldsymbol{u}^i}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \varPhi^i \rho^p \boldsymbol{u}^i\right), \end{gather}

(3.3)\begin{gather}\displaystyle\frac{\partial\epsilon \rho^f \boldsymbol{u}^f}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left( \epsilon \rho^f \boldsymbol{u}^f \otimes \boldsymbol{u}^f \right) = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{S}^f - \epsilon \rho^f \boldsymbol{g} - n_l \boldsymbol{f}_{f\rightarrow l} - n_s\boldsymbol{f}_{f\rightarrow s} \end{gather}

(3.4)\begin{gather}\displaystyle\frac{\partial\varPhi^i \rho^p \boldsymbol{u}^i}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left( \varPhi^i \rho^p \boldsymbol{u}^i \otimes \boldsymbol{u}^i \right) = \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{S}^i - \varPhi^i \rho^p \boldsymbol{g} + n_i \boldsymbol{f}_{f\rightarrow i} + n_i \boldsymbol{f}_{\delta\rightarrow i}, \end{gather}

where $f$ is the fluid and indices $i=l,s$ denote the large particle phase and the small particle phase, respectively ($\delta = l$ if $i = s$ and $\delta = s$ if $i = l$). Here, $\varPhi ^l$ and $\varPhi ^s$ are the volume fractions for the large and small grains and verify $\varPhi ^s + \varPhi ^l = \varPhi$ where $\varPhi$ is the volume fraction of the mixture, i.e. the total solid volume fraction. Consequently, the fluid volume fraction is $\epsilon = 1- \varPhi ^l - \varPhi ^s$ and $\boldsymbol {S}^k$ is the stress tensor associated with phase $k$ with $k = l, s \text { or } f$. They can be separated into pressure and shear stress contribution

(3.5)\begin{equation} \boldsymbol{S}^k ={-}p^k \boldsymbol{I} + \boldsymbol{\tau}^k, \end{equation}

where $\boldsymbol {\tau }^k$ is the shear stress tensor and $p^k$ is the pressure of phase $k$. It should be noted that, for a solid phase, the static pressure arises from the enduring contacts between the particles. Thanks to the mixture model approaches (Morland Reference Morland1992), it can be assumed that each particle phase carries the total overburden pressure $p^m$ according to their local volume fraction as

(3.6)\begin{equation} p^i = \dfrac{\varPhi^i}{\varPhi} p^m, \end{equation}

where $m$ denotes the mixture made of small and large particles. The total overburden pressure $p^m$ is computed using the formulation proposed by Johnson & Jackson (Reference Johnson and Jackson1987). For further information, the reader is referred to Chauchat et al. (Reference Chauchat, Cheng, Nagel, Bonamy and Hsu2017) and Chauchat (Reference Chauchat2018).

The momentum equations (3.3) and (3.4) contain two terms coming from the momentum exchange between the different phases: $n_i \boldsymbol {f}_{f\rightarrow i}$ and $n_i \boldsymbol {f}_{\delta \rightarrow i}$. The term $n_i \boldsymbol {f}_{f\rightarrow i}$ is the averaged value of the resultant forces exerted by the fluid on the particles of phase $i$. Jackson (Reference Jackson2000) showed that, for a collection of immersed particles, this interaction force can be written as

(3.7)\begin{equation} n_i \boldsymbol{f}_{f\rightarrow i} ={-} \varPhi^i \boldsymbol{\nabla} p^f + n_i \boldsymbol{f}^{f\rightarrow i}_{d}, \end{equation}

where $\varPhi ^i \boldsymbol {\nabla } p^f$ is the buoyancy force exerted by the fluid phase on the particles and $n_i \boldsymbol {f}^{f\rightarrow i}_{d}$ is the particle-averaged viscous drag force between the particles and the fluid phase. The term $n_i \boldsymbol {f}_{\delta \rightarrow i}$ is the averaged value of all interacting forces between large and small particle phases. It can be directly expressed in three dimensions from the local segregation force of Tripathi & Khakhar (Reference Tripathi and Khakhar2013) and Guillard et al. (Reference Guillard, Forterre and Pouliquen2016).

Therefore, the developed model is general and can be applied to 3-D configurations. For simplicity and for the purposes of the present study, the model will be only developed for a 1-D uniform flow.

3.2. Simplified 1-D vertical multi-phase flow model

The multi-phase flow model (3.1) to (3.4) is simplified by considering a uniform flow in the streamwise direction. From now, all the variables only depend on the vertical position $z$. Therefore, the spatially averaged velocity of the phase $k$ can be written as $\boldsymbol {u}^k = u^k(z) \boldsymbol {e}_x + w^k(z) \boldsymbol {e}_z$. The mass conservation equations simplify to

(3.8a,b)\begin{equation} \displaystyle\frac{\partial\epsilon}{\partial t} + \displaystyle\frac{\partial\epsilon w^f}{\partial z} = 0 \quad \text{and} \quad \dfrac{\partial {\varPhi^i}}{\partial {t}} + \dfrac{\partial{{\varPhi^i}w^i}}{\partial{z}}=0, \end{equation}

and the momentum balance equations in the vertical direction are

(3.9)\begin{gather} \rho^f \left[\displaystyle\frac{\partial\epsilon w^f}{\partial t} + \displaystyle\frac{\partial\epsilon w^f w^f }{\partial z}\right] ={-}\epsilon \displaystyle\frac{\partial p^f}{\partial z} - \epsilon \rho^f g \cos \theta - n_l\langle f^{f\rightarrow l}_{d}\rangle - n_s\langle f^{f\rightarrow s}_{d}\rangle, \end{gather}

(3.10)\begin{gather}\rho^p \left[\displaystyle\frac{\partial\varPhi^l w^l}{\partial t} + \displaystyle\frac{\partial\varPhi^l w^l w^l}{\partial z}\right] ={-}\displaystyle\frac{\partial p^l}{\partial z} - \varPhi^l \displaystyle\frac{\partial p^f}{\partial z} - \varPhi^l \rho^p g \cos \theta + n_l\langle f^{f\rightarrow l}_{d}\rangle + n_l \langle\, f_{s\rightarrow l}\rangle , \end{gather}

(3.11)\begin{gather}\rho^p \left[\displaystyle\frac{\partial\varPhi^s w^s}{\partial t} + \displaystyle\frac{\partial\varPhi^s w^s w^s}{\partial z}\right] ={-}\displaystyle\frac{\partial p^s}{\partial z} - \varPhi^s \displaystyle\frac{\partial p^f}{\partial z} - \varPhi^s \rho^p g \cos \theta + n_s\langle f^{f\rightarrow s}_{d}\rangle + n_s \langle\, f_{l\rightarrow s}\rangle. \end{gather}

In the two last equations, the solid pressures $p^l$ and $p^s$ are given by (3.6). To solve these equations, it is necessary to prescribe closures for the spatially averaged fluid/grain interaction and grain–grain interactions, and for the granular and fluid pressures.

Considering the fluid–grain interaction, both small and large granular phases interact with the fluid phase through $\varPhi ^i \partial p^f/ \partial z$ and the drag force $n_i\langle f^{f\rightarrow i}_{d}\rangle$. For an assembly of particles, the spatial averaging of the vertical total drag force applied by the fluid gives

(3.12)\begin{equation} n_i \langle f^{f\rightarrow i}_{d}\rangle = \dfrac{\varPhi^i \rho^p }{t_i}\left(w^f-w^i\right), \end{equation}

where $t_i = \rho ^p d_i^2 (1-\varPhi )^{3} / 18 \eta ^f$ is the particle response time and $d_i$ is the particle diameter of phase $i$. The factor $(1-\varPhi )^{3}$ is a correction proposed by Richardson & Zaki (Reference Richardson and Zaki1954) to take into account hindrance effects. Since the drag is linear, the spatial averaging is simply the drag force applied on one particle (given in (2.4)) multiplied by the number of particles per unit volume $n_i = \varPhi ^i/ V^i$ (Jackson Reference Jackson2000), where $V^i$ is the volume of a single particle of phase $i$.

The granular phases also interact with each other and the grain–grain interaction closure should be prescribed in the model. For a single large grain in a bath of small particles, it has been shown in § 2 that small particles exert three forces on a large intruder,

(3.13)\begin{equation} f_{s\rightarrow l} = f_{d}^p + \varPi_p + f_{l}. \end{equation}

In the case of spatially averaged equations, the force balancing the weight $\varPi _p$ already appears in the term $-{\partial p^l}/{\partial z}$. Hence, the granular interaction forces of (3.13) reduce to the granular drag and the granular lift force.

To extend these forces to a collection of large particles, the interaction force $f_{s\rightarrow l}$ is spatially averaged. Since this force is linear, it amounts to multiplying $f_{s\rightarrow l}$ by the number of large particles per unit volume $n_l = 6 \varPhi ^l / {\rm \pi}d_l^3$. Therefore, the total solid interaction force exerted by the small particles on the large ones is given by

(3.14)\begin{equation} n_l \langle\, f_{s\rightarrow l}\rangle = \dfrac{\varPhi^l \rho^p}{t_{ls}} \left(w^s-w^l\right) + \varPhi^l \mathcal{F}_{l}(\mu) \displaystyle\frac{\partial p^m}{\partial z}, \end{equation}

where $t_{ls} = \rho ^p d_l^2 / 6 c \eta ^p$ is the particle response time for the drag force between small and large particles. According to Newton's third law, the force exerted by the large particles on the small ones (3.11) is

(3.15)\begin{equation} n_s \langle\, f_{l\rightarrow s}\rangle ={-} n_l \langle\, f_{s\rightarrow l}\rangle . \end{equation}

The solid mixture phase is made of both particle phases and is noted with $i=m$. Its momentum balance is obtained by summing (3.10) and (3.11). Since the mixture phase does not distinguish between small and large particles, the solid interaction forces should not appear in this equation. Equation (3.15) ensures that these forces vanish when developing the mixture momentum equation.

The proposed volume-averaged multi-phase flow model describes size segregation of a bi-disperse mixture immersed in a fluid. This represents an improvement upon the model of Thornton et al. (Reference Thornton, Gray and Hogg2006), which was based on semi-empirical parametrisation of the interparticle forces between small and large particles. The present model provides closures based on forces applied on a single particle, and bridges the gap between granular-scale processes and continuum modelling in size segregation. This important result will be used in the following to derive an advection–diffusion model for size segregation.

5.1. DEM investigation of Chassagne et al. (Reference Chassagne, Maurin, Chauchat, Gray and Frey2020b)

In this section, the configuration explored and the main results of Chassagne et al. (Reference Chassagne, Maurin, Chauchat, Gray and Frey2020b) are briefly presented. The authors investigated grain-size segregation in turbulent bedload transport using a coupled fluid–DEM originally developed by Maurin et al. (Reference Maurin, Chauchat, Chareyre and Frey2015) and used to study bedload rheology (Maurin et al. Reference Maurin, Chauchat and Frey2016) and the slope influence (Maurin, Chauchat & Frey Reference Maurin, Chauchat and Frey2018). For further details, the interested reader is referred to Chassagne et al. (Reference Chassagne, Maurin, Chauchat, Gray and Frey2020b). The 3-D bi-periodic DEM set-up consisted in depositing a layer of small particles over large ones, and letting the particles be entrained by the fluid flow at a fixed Shields number. The latter is the dimensionless fluid bed shear stress $Sh = \tau ^f /[(\rho ^p-\rho ^f)g d_l]$ and was taken equal to $0.1$. The bed slope was fixed to $10\,\%$, which is representative of mountain streams. The size ratio was taken as $r=1.5$ with small particles of diameter $d_s = 4\ \textrm {mm}$ and large particles of diameter $d_l = 6\ \textrm {mm}$. The large and small particles were assimilated to a number of layers, $N_l$ and $N_s$. The number of layers of a given class represents the height, in terms of particle diameter of this class, occupied by particles if the concentration was equal to the random close packing ($\varPhi _{max} = 0.61$). In this way, the bed height at rest was defined as $h = N_l d_l + N_s d_s$ and was fixed to $h=10 d_l$ with a random close packing volume fraction $\varPhi _{max} = 0.61$ (profile in figure 2a). In the study of Chassagne et al. (Reference Chassagne, Maurin, Chauchat, Gray and Frey2020b), different simulations have been performed with $N_s$ varying from $0.01$ (a few isolated particles) to $N_s=2$. In this section it was decided that the comparison would be made with $N_s = 1.5$. The bulk response of the granular mixture to this fluid forcing is represented by the dimensionless mixture streamwise velocity profile in figure 2(a). The inset is a semi-log plot of the dimensionless velocity profile and shows that it is exponential. As shown in figure 2(b), the linearity of the curve in the semi-log plot confirms that the shear rate is exponentially decreasing in the quasi-static part of the bed (delimited by the two horizontal black dashed lines). As expected for a uniform flow, the mixture shear stress $\tilde {\tau }_{xz}^m$ shown in figure 2(c) is linear with depth. For both quantities, the following fits were proposed and plotted as a red dotted line in figure 2:

(5.1a,b)\begin{equation} \tilde{\dot{\gamma}}^m = \gamma_0 \,\textrm{e}^{\tilde{z}/s_0} \quad \text{and} \quad \tilde{\tau}^m_{xz} = a_0 \tilde{z} + \tau_0 .\end{equation}

The simulations performed by Chassagne et al. (Reference Chassagne, Maurin, Chauchat, Gray and Frey2020b) in this configuration were focused on the downward segregation of small particles. It was observed that the layer of small particles percolates rapidly for $\tilde {z}> 8.5$ (flowing layer) and then slows down below (this limit is marked in figures 2(b) and 2(c)). Chassagne et al. (Reference Chassagne, Maurin, Chauchat, Gray and Frey2020b) showed that the small particles are advected downward like a travelling wave into the bed made of large particle with a layer of constant thickness. As figure 2(e) shows, the small particle concentration has a Gaussian-like shape and remains self-similar in time while segregating. The centre of mass of small particles, $\tilde {z}_c$, is therefore representative of the dynamics of the entire layer. Chassagne et al. (Reference Chassagne, Maurin, Chauchat, Gray and Frey2020b) observed that the small particle layer travels down as a logarithmic function of time

(5.2)\begin{equation} \tilde{z}_c(t) ={-}a_1 \ln \tilde{t} + b, \end{equation}

where $a_1$ is a constant characterising the segregation velocity ($\textrm {d}\tilde {z}_c(t)/\textrm {d} \tilde {t} = - a_1/\tilde {t}$). The authors demonstrated that this logarithmic descent of small particles is a consequence of the dependency of the segregation velocity on the inertial number as

(5.3)\begin{equation} \dfrac{\textrm{d} \tilde{z}_c}{\textrm{d} t} \propto I^{0.85}(\tilde{z}_c). \end{equation}

Using (5.3) in the framework of the advection–diffusion model of Thornton et al. (Reference Thornton, Gray and Hogg2006) and Gray & Chugunov (Reference Gray and Chugunov2006) (see (1.5)) it was shown that the advection coefficient could be written as

(5.4)\begin{equation} S_{r} = S_{r0} I^{0.85}, \end{equation}

where $S_{r0} = 0.049$, extending the inertial number dependency found by Fry et al. (Reference Fry, Umbanhowar, Ottino and Lueptow2018) to bi-disperse size segregation in the quasi-static regime. Then, with the help of a travelling wave method, they evidenced that the small particles percolate as a layer and with a self-similar concentration profile because the ratio between the advection coefficient and the diffusion coefficient is constant. The Péclet number reads

(5.5)\begin{equation} Pe = \dfrac{S_{r}}{D}, \end{equation}

and is therefore constant with $\tilde {z}$, so that the diffusion coefficient has to have the same dependency on the inertial number as the segregation coefficient

(5.6)\begin{equation} D=D_{0} I^{0.85}, \end{equation}

where $D_0$ is taken as $D_0 = 0.01$ and should be pressure independent (Fry et al. Reference Fry, Umbanhowar, Ottino and Lueptow2019).

Figure 2.

Profiles and configuration from the DEM simulations. (a) Streamwise mixture velocity profile in the bed (——, blue solid line) and mixture volume fraction (——, green solid line). The inset is the semi-log plot of the velocity profile. (b) Solid shear rate (——, blue solid line) and the corresponding fit $\tilde {\dot {\gamma }}^m = \gamma _0 \,\textrm {e}^{\tilde {z}/s_0}$ (- - - -, red dashed line) with $\gamma _0 = 1.64 \times 10^{-7}$ and $s_0 = 0.74$. (c) Solid shear stress and the corresponding fit $\tilde {\tau }_{xz}^m = a_0 \tilde {z} + \tau _0$ (- - - -, red dashed line) with $a_0=-0.078$ and $\tau _0=0.91$. The top and lower boundaries of the quasi-static bed are represented by (- - - -, black dashed lines). (d) Sketch of the numerical experiment with the input profiles for the rheology. (e) Concentration profile of small particles at the initial state for the multi-phase flow model (——, green solid line), the DEM (- - - -, black dashed line) and the mixture concentration profile $\phi$ (——, black solid line).

This work also demonstrated that the dynamics of the fine particle layer is controlled by its bottom position, which acts as a lower bound for the segregation velocity. In this way, the particles in the layer cannot segregate faster.

The numerical resolution of the 1-D multi-phase model a priori requires us to solve the granular rheology in order to estimate the granular viscosity required to evaluate the granular drag force contribution. The goal of the present study is to focus on grain-size segregation modelling. In addition, the results of Chassagne et al. (Reference Chassagne, Maurin, Chauchat, Gray and Frey2020b) were obtained in the quasi-static regime, for which there is still no consensus regarding granular rheology. For these reasons, the granular viscosity is directly determined by the DEM results here. This makes it possible to focus on the effect of the segregation model and to put aside potential discrepancies linked to a non-accurate description of the granular rheology. The friction coefficient is therefore computed from the DEM results using the definition

(5.7)\begin{equation} \mu = \dfrac{|\tau^m_{xz}|}{p^m}. \end{equation}

Similarly, the granular viscosity is computed using the definition

(5.8)\begin{equation} \eta^p = \dfrac{|\tau^m_{xz}|}{\left\vert \dot{\gamma^m_{xz}}\right\vert}. \end{equation}

As shown by the red dashed lines in figures 2(b) and 2(c), the fits presented in (5.1a,b) match the DEM results in the region $8.5>\tilde {z}>3$. Thus, using these expressions in (5.8) provides an accurate estimate of the granular viscosity for the particle–particle drag closure. For this reason, the validation of the multi-phase flow model will only be carried out in this part of the bed. Therefore, the initial state consists in placing the small particles in the upper limit of the quasi-static part with the centre of mass $\tilde {z}^0_c = 8.5$. The sketch of this configuration is shown in figure 2(d). Figure 2(e) shows the small particle concentration profile of this numerical set-up at the initial state. The initial concentration is taken with a Gaussian fit on the DEM initial concentration (figure 2e), and ensures that the mass of particles is the same in the DEM and in the continuum simulations.

5.2. Comparison with the multi-phase flow model

The system of partial differential equations (3.8a,b)–(3.11) is solved numerically for the configuration shown in figure 2(d) with the initial concentration profile of figure 2(e). In these equations, the empirical segregation function is $\mathcal {F}_l(\mu )$ and the drag coefficient $c$ is equal to 3, as suggested by Tripathi & Khakhar (Reference Tripathi and Khakhar2013). Because the fluid is incompressible, there is no equation of state for the fluid pressure. Nevertheless, remembering that $\epsilon +\varPhi =1$ and defining the volume-averaged velocities, $w = \epsilon w^f+ \varPhi w^m$, it can be demonstrated that the particle–fluid mixture is incompressible. A PISO (pressure implicit with splitting of operators) algorithm classically developed for thw incompressible Navier–Stokes equations is used to solve the pressure–velocity coupling. As the fluid pressure $p^f$ is the sum of the hydrostatic pressure and of the excess pore pressure $p^f = \overline {p^f} + \rho ^f g z$, the PISO algorithm consists in solving the momentum balance equations without $\overline {p^f}$ in a predictor step. Then, using the predicted velocity fields, a Poisson equation is solved to find $\overline {p^f}$. Once the pressure is found, the velocity fields are corrected. This kind of algorithm has already been used to model sediment transport in Chauchat et al. (Reference Chauchat, Cheng, Nagel, Bonamy and Hsu2017) and Chauchat (Reference Chauchat2018).

Figure 3 shows the results of the spatio-temporal evolution of the small particle concentration for the multi-phase flow model (figure 3a) and for the DEM (figure 3b). First, it can be seen that the dynamics predicted by the multi-phase flow model is similar to the DEM. The position of the bottom of the layer is approximately the same in both cases. More quantitatively, the centre of mass $\tilde {z}_c$ of the small particle layer as a function of time is compared with the DEM in figure 4(a). After a first transient phase ($\tilde {t}>1\times 10^3$), the centre of mass position is linear in the semi-log plot, indicating that the logarithmic descent observed in the DEM simulation is well reproduced by the multi-phase flow model. The slope of the curve, representing coefficient $a_1$ of (5.2) is $0.68$ in the DEM simulation and $0.49$ for the multi-phase flow model, corresponding to an error of $28\,\%$. In addition, figure 4(b) shows that, in both models, the bottom of the layer is positioned at the same depth indicating that the multi-phase model reproduces well the bottom controlled behaviour observed by Chassagne et al. (Reference Chassagne, Maurin, Chauchat, Gray and Frey2020b) with DEM simulations. However, the Gaussian-like profile is not reproduced by the multi-phase flow model and a wider profile is obtained. In figure 3(a) the maximum concentration $\max (\phi ^s)$ (indicated by $\blacksquare$ in figure 4b) is almost two times smaller than the one predicted by the DEM simulation, while the extent of the small particle layer is much larger. These results indicate that, with the current parametrisation the multi-phase flow model is relevant to qualitatively predict the segregation dynamics. However, the error on the segregation velocity and the discrepancies on the concentration profile clearly show that the model needs to be improved.

Figure 3.

Spatio-temporal plot of the small particle concentration in the bed obtained with (a) the multi-phase model, (b) the DEM simulation.

Figure 4.

Results of the simulation for the multi-phase flow model (-$\cdot$-$\cdot$-, red dashed dot), the advection–diffusion model of (4.8) (——, green solid line) and the DEM (- - - -, black solid line). (a) Temporal evolution of the centre of mass. (b) Concentration profiles of small particles at the end of the simulation ($\tilde {t} \simeq 60\,000$); $\blacksquare$, represents the maximum concentration of the profiles.

5.3. Evaluation of the advection–diffusion model

In order to determine the ability of the advection–diffusion model to reproduce the same results as the multi-phase flow model, (4.8) is solved numerically. The resolution strategy is based on a conservative Godunov schemes where a no flux condition is applied at the bottom and on top of the vertical domain. A vertical discretisation of $\textrm {d}\tilde {z} = h / 80$ is taken and the time step is computed in order to satisfy the Courant–Friedrichs–Lewy (CFL) condition. The initial solution is the same than in the multi-phase flow (figure 2e).

The numerical solution at time $\tilde {t} = 60\,000$ is plotted in figure 4. Both the centre of mass (figure 4a) and the concentration profile (figure 4b) are almost superimposed with the multi-phase flow model solution, only slightly differing by numerical diffusion. Therefore, both models can be considered as strictly equivalent, meaning that the physical behaviour of the segregation forces added to the momentum equation of the small particles is accurately predicted by the single advection and diffusion coefficients of equation (4.8). Moreover, when deriving the advection–diffusion equation, it was assumed that the mixture volume fraction $\varPhi = cste$, the vertical acceleration of the small particles, the vertical mixture velocity and the fluid drag were negligible. The strong agreement between the models demonstrate that these assumptions are valid.

Equation (4.8) represents an important step in the upscaling since the behaviour of small particles can be predicted without solving the entire multi-phase flow model, providing a speed-up of one thousand for the numerical resolution of the equations. In the light of this result, the bi-disperse segregation problem can be simply viewed as a competition between the advection coefficient $S_r$ and the diffusion coefficient $D$ of small particles. The interplay between both coefficients has been explored in more details by Fan et al. (Reference Fan, Schlick, Umbanhowar, Ottino and Lueptow2014b).