Segregation–rheology feedback in bidisperse granular flows: a coupled Stokes’ problem

Abstract

The feedback between particle-size segregation and rheology in bidisperse granular flows is studied using the Stokes’ problem configuration. A method of lines scheme is implemented to solve the coupled momentum and segregation equations for a normally graded particle size distributed bulk at constant solids volume fraction. The velocity profiles develop quickly into a transient state, decoupled from segregation yet determined by the particle size. From this transient state, the velocity profile changes due to the particles’ relative movement, which redistributes the frictional response, hence its rheology. Additionally, the particles’ relative friction is modified via a frictional coefficient ratio, by analogy with the particles’ size ratio. While positive values of this coefficient exacerbate the nonlinearity of the velocity profiles induced by size differences, negative values dampen this behaviour. The numerical solutions reproduce well the analytical solutions for the velocity profile, which can be obtained from the steady-state conditions of the momentum and segregation equations for the transient and steady states, respectively. Segregation–momentum balances and four characteristic time scales can be established to propose two non-dimensional quantities, including specific Schmidt and Péclet numbers that describe broadly the segregation–rheology feedback. The proposed scheme, theoretical solutions and non-dimensional numbers offer a combined approach to understand segregation and flow dynamics within a granular bulk, extensible across many flow configurations.

1. Introduction

Size segregation is responsible for important mobility feedback in granular flows (Barker et al. 2017; Gray 2018; Umbanhowar, Lueptow & Ottino 2019; Edwards et al. 2023). Since particle size $d$ is a defining parameter for most rheological models of dry granular materials (GDR-MiDi 2004; Jop, Forterre & Pouliquen 2006; Kamrin & Koval 2012; Henann & Kamrin 2013; Kim & Kamrin 2020), segregation-induced grain movement alters the very same mechanical response that generates segregation in the first place. Such feedback determines the flow's runout, the formation of levees, fingers, and blunt or bulbous fronts in debris flows, all tied to coarse grains segregating towards the flow's free surface (Woodhouse et al. 2012; Denissen et al. 2019; Edwards et al. 2023).

Numerical simulations and experiments have been pivotal to shed light on the rich segregation–rheology coupling in granular flows. Recent experiments in two- and three-dimensional shear boxes unravelled a scaling law for bidisperse particle-size segregation; asymmetric segregation fluxes depend on shear rate $\dot {\gamma }$, particle sizes, their size ratio $R_{d}$ and pressure $p$ (i.e. the inertial number $I$, implicitly), and local particle concentrations (Trewhela, Ancey & Gray 2021a; Trewhela, Gray & Ancey 2021b). These experimental results are consistent with previous observations on pressure-dependent segregation (Golick & Daniels 2009) and the intrinsic asymmetric nature of segregation fluxes (van der Vaart et al. 2015; Jones et al. 2018). Similar results, in terms of an inertial number $I$ dependence, were obtained in discrete element method simulations, by modelling a bedload transport flow (Chassagne et al. 2020) and a shear flow (Bancroft & Johnson 2021). Based on the empirical scaling law of Trewhela et al. (2021a), Barker et al. (2017) performed simulations of fully coupled granular flows down a slope and in a square rotating drum. They discussed the intricate segregation–rheology feedback, concluding that it is not clear how subtle changes in flow properties influenced the coupling mechanisms and their impact on the flows’ overall dynamics.

In 1850, G.G. Stokes proposed a problem concerning a semi-infinite layer of fluid in contact with a horizontal and infinitely long plate. The plate suddenly slides horizontally at $t=0$ with finite speed, therefore transferring shear to the fluid (Stokes 1850). This problem, known as Stokes’ first problem (or Rayleigh's problem), offers an idealized, yet practical configuration for the study of sheared flows (Rayleigh 1911). Besides the historical and theoretical interest that this problem has had for fluid mechanics (Stewartson 1951; Tokuda 1968; Preziosi & Joseph 1987; Jordan & Puri 2005; Jerome & Di Pierro 2018), its flow configuration has various applications for experimental and industrial facilities using granular materials, e.g. shear cells or conveyor belts (Golick & Daniels 2009; van der Vaart et al. 2018; Trewhela & Ancey 2021). Using the $\mu (I)$-rheology (GDR-MiDi 2004; Jop et al. 2006), Jerome & Di Pierro (2018) studied Stokes’ problem for granular flows, developing numerical and analytical solutions for momentum transfer and velocity profiles within a semi-infinite granular medium. Their analysis focused on a sheared layer $h$, which developed when applying shear $\tau _{w}$ at the upper boundary, and depended on a kinematic granular viscosity $\nu _{g}\sim d$, an observation also recovered with kinetic theory in gravitational granular currents (e.g. Larcher & Jenkins 2019) and discrete element method simulations (Jing et al. 2022). Intuitively, changes in particle-size distribution and subsequent segregation will significantly alter momentum diffusion within the bulk, so coupling this process with segregation is not straightforward.

This work investigates the segregation–rheology feedback, its coupling and impact on the overall dynamics of dense granular flows. Stokes’ first problem configuration, with the theoretical and numerical framework for the segregation–rheology feedback, is introduced in § 2. The numerical and steady-state solutions are presented in § 3, where emphasis is put on the crucial role of size ratio and frictional coefficient differences in the flow dynamics. Analytical solutions and characteristic time scales are sought, so that non-dimensional numbers are proposed to describe the segregation–rheology feedback. The conclusions in § 4 put focus on the found numerical and analytical solutions, the description provided by the non-dimensional balances and their applicability in general flow configurations.

2. Governing equations and numerical solutions

2.1. Stokes’ problem for a dense bidisperse granular flow

Consider a semi-infinite granular material formed by a dense bidisperse mixture of grains that share the same density $\rho _{*}$ but differ in size (figure 1). The particle sizes $d_{s}$ and $d_{l}>d_{s}$ define grain species $\nu =\{ s,l\}$ of small and large particles, respectively. The contribution of each species to the solids volume fraction $\varPhi$ is set by its partial volume concentration $\phi _{\nu }$ per unit volume of grains, so that $\phi _{s}+\phi _{l}=1$. A plate of infinite length lies initially at rest over the bulk and is set to slide horizontally with velocity $u_{w}$ at $t=0$, and a frictional coefficient $\mu _{w}$. Under these conditions and at constant $\varPhi$, such that $\rho =\varPhi \rho _{*}$, the mass and momentum conservation equations for the whole granular material (small and large particles) in the $x\unicode{x2013}z$ directions reduce to

(2.1a–c)\begin{equation} \frac{\partial u}{\partial x}=0,\quad \rho_{*}\varPhi\,\frac{\partial u}{\partial t}=-\frac{\partial\tau_{xz}}{\partial z}, \quad \frac{\partial p}{\partial z}=\rho_{*}\varPhi g,\end{equation}

where $g$ is the gravitational acceleration in the $z$-direction, $u$ is the component of the velocity field $\boldsymbol {u}$ in the $x$-direction, and $\tau _{xz}$ is the $x\unicode{x2013}z$ component of the deviatoric shear stress tensor $\boldsymbol {\tau }$. In dense granular flows, pressure $p$ and shear stress $\tau$ relate through a Coulombic frictional coefficient $\mu$ so that $\tau =\mu p$. From (2.1c), pressure $p=p_{w}+\rho _{*} gz$ can be considered lithostatic with $p_{w}$ imposed at the plate above. For various flow configurations, GDR-MiDi (2004) showed that $\mu$ depends on the dimensionless inertial number $I=\dot {\gamma }d/\sqrt {p/\rho _{*}}$, with $d$ the particle diameter and $\dot {\gamma }=2\,\|\boldsymbol {D}\|=|\text {d}u/\text {d}z|$ the shear rate. Barker & Gray (2017) introduced a partial regularization for the $\mu (I)$-rheology (Jop et al. 2006), so that its well-posedness range could be extended. This regularization can be introduced by considering the extended function

(2.2)\begin{equation} \mu(I)=\left\{\begin{array}{@{}ll} \sqrt{\alpha/\log (A/I)}, & \text{for } I\leq I_{1},\\[2pt] \dfrac{\mu_{1}I_{0}+\mu_{2}I+\mu_{\infty}I^{2}}{I+I_{0}}, & \text{for } I> I_{1}, \end{array}\right. \end{equation}

where $\mu _{1}$ is the static friction coefficient, $\mu _{2}$ is the dynamic friction coefficient, and $I_{0}$ are the original parameters employed by Jop et al. (2006), corresponding to particle-dependent or material constants. The remaining parameters $\mu _{\infty }$, $I_{1}$ and $\alpha$ are also material-dependent and were introduced by Barker & Gray (2017). While $\mu _\infty$ holds physical meaning in the context of rapid granular flows beyond $\mu _{2}$ (Holyoake & McElwaine 2012; Barker & Gray 2017), $\alpha$ and $I_{1}$ are just mathematical fitting constants. Finally, $A$ is a constant that guarantees continuity of the regularized function and is calculated using the latter parameters. All the frictional coefficients and material constants used for this work follow those in Barker et al. (2017) and are detailed in table 1.

Figure 1. Stokes’ problem for a dense bidisperse granular flow: an infinite plate shears a granular bulk made of small ($\nu =s$) and large ($\nu =l$) particles (each particle species is characterized by its diameter $d_{\nu }$ and partial volume concentration $\phi _{\nu }$). The plate suddenly moves with constant velocity $u_{w}$ at $t=0$, imposing pressure $p_{w}$ with a frictional coefficient $\mu _{w}$ at contact with the granular media. (a) An inversely graded initial condition (segregated, but unstable) for the granular bulk, with partial concentrations $\phi _{s}=1$ on top and $\phi _{s}=0$ ($\phi _{l}=1-\phi _{s}=1$) at the bottom. (b) Initially mixed granular bulk (unstable) condition, with uniform small particle concentration $\phi _{s}^{0}$ and an equivalent concentration-averaged diameter $\bar {d}=\phi _{s}d_{s}+\phi _{l}d_{l}$ through the bulk's sheared layer $h$.

Table 1. Intrinsic density $\rho _{*}$, solids volume fraction $\varPhi$, and frictional parameters $\mu _{1}$, $\mu _{2}$ and $I_{0}$ for the $\mu (I)$-rheology (Jop et al. 2006). Frictional coefficient $\mu _{\infty }$ with the parameters $\alpha$, $I_{1}$ and $A$ (computed) for the partial regularization of the $\mu (I)$-rheology were proposed and measured by Barker & Gray (2017) for $d_{l}=143\ \mathrm {\mu }$m glass beads.

The bidisperse nature of the considered granular material can be included via (i) a concentration-averaged diameter $\bar {d}$ in $I$ (Tripathi & Khakhar 2011; Chassagne et al. 2020; Bancroft & Johnson 2021), and (ii) a concentration-averaged frictional coefficient $\bar {\mu }$ in (2.2). Similarly to the concentration-averaged diameter $\bar {d}=\phi _{s}d_{s}+\phi _{l}d_{l}$ (Rognon et al. 2007; Tripathi & Khakhar 2011), a concentration-averaged $\bar {\mu }$ definition can be proposed in the form $\bar {\mu }=\mu _{s}\phi _{s}+\mu _{l}\phi _{l}$, where $\mu _{\nu }$ is the frictional coefficient for the $\nu$ species. Arranging the terms, this $\bar {\mu }$ definition can be expressed as a function of the ratio between the species’ frictional coefficients $R_{\mu }=\mu _{l}/\mu _{s}$, resulting in

(2.3)\begin{equation} \bar{\mu}=(1-\chi_{\mu}\phi_{s})R_{\mu}\mu_{s},\end{equation}

which prompts the definition $\chi _{\mu }=(R_{\mu }-1)/R_{\mu }\in (-\infty,1]$. This expression assumes that $I_{0}$ is invariant with size, and only the frictional coefficients $\mu _i$ ($i=\{1,2,\infty \}$) in table 1 are size-dependent.

As expected, the inclusion of the $\mu (I)$-rheology poses several difficulties for numerical schemes, particularly when the bulk is at rest or in steady state (Cawthorn 2011; Barker & Gray 2017; Jerome & Di Pierro 2018). To overcome this issue, the analysis can be focused on a sheared layer $h$, defined as

(2.4)\begin{equation} h=\left(\frac{\mu_{w}}{\mu_{1}}-1\right)\frac{p_{w}}{\rho g \varPhi}, \end{equation}

where $\mu _{w}$ and $p_{w}$ are the frictional coefficient and pressure applied at the upper plate, respectively. This shear layer $h$ definition was already employed in the solution of the monodisperse granular Stokes’ problem, but still poses difficulties for numerical schemes and solutions (Cawthorn 2011; Jerome & Di Pierro 2018).

2.2. Bidisperse particle-size segregation

Particle-size segregation in sheared dense granular flows is driven primarily by gravity and friction. This gravity-driven segregation is the result of two mechanisms: small particles fall through gaps via kinetic sieving (Middleton 1970), and large particles are squeezed by surrounding particles, expulsed from their position (Savage & Lun 1988). There are, however, other segregation mechanisms, many of which appear in different configurations and can be described employing kinetic theory and quantities such as granular temperature (e.g. Larcher & Jenkins 2019; Neveu et al. 2022). In this work however, segregation mechanisms are not discussed, and focus is put on the shear and gravity-driven particle-size segregation problem for dense granular flows under a frictional regime (Ancey & Evesque 2000; GDR-MiDi 2004). The segregation dynamics to be studied with the Stokes’ problem is well described by the mentioned shear-driven mechanisms and well documented in theoretical, numerical and experimental work (e.g. May et al. 2010; van der Vaart et al. 2015; Guillard, Forterre & Pouliquen 2016; Trewhela & Ancey 2021; Edwards et al. 2023), and has a determinant role for various geophysical granular flows (Delannay et al. 2017; Gray 2018).

For the Stokes’ problem in figure 1(a), the bidisperse segregation equations for $\phi _{s}$ and $\phi _{l}$ are

(2.5a,b) \begin{align} \frac{\partial \phi_{s}}{\partial t}+\boldsymbol{\nabla} \boldsymbol{\cdot} \left(f_{sl}\phi_{s}\phi_{l}\,\frac{\boldsymbol{g}}{|\boldsymbol{g}|}\right) = \boldsymbol{\nabla} \boldsymbol{\cdot} (\mathcal{D}_{sl}\,\boldsymbol{\nabla} \phi_{s}),\quad \frac{\partial \phi_{l}}{\partial t}-\boldsymbol{\nabla} \boldsymbol{\cdot} \left(f_{sl}\phi_{s}\phi_{l}\,\frac{\boldsymbol{g}}{|\boldsymbol{g}|}\right) = \boldsymbol{\nabla} \boldsymbol{\cdot} (\mathcal{D}_{sl}\,\boldsymbol{\nabla} \phi_{l}), \end{align}

where $f_{sl}$ is the segregation velocity magnitude, and $\mathcal {D}_{sl}=\mathcal {A}\dot {\gamma }\bar {d}^{2}$ is the particles’ diffusivity, with $\mathcal {A}=0.108$ and $\bar {d}$ a mean diameter (Utter & Behringer 2004). Note that the inclusion of the gravitational acceleration considers the fact that segregation in dense granular flows is a shear and gravity-driven process due mainly to the kinetic sieving mechanism (Middleton 1970). Recently, Trewhela et al. (2021a) presented an experimentally based scaling law for $f_{sl}$, obtaining

(2.6)\begin{equation} f_{sl}=\left(\frac{\hat{\rho}\hat{\mathcal{B}}\rho_{*}g \dot{\gamma}\bar{d}^{2}}{\mathcal{C}\rho_{*}g\bar{d}+p}\right) [(R_{d}-1)+\mathcal{E}(1-\phi_{s})(R_{d}-1)^{2}],\end{equation}

where $R_{d}=d_{l}/d_{s}$ is the particles’ size ratio, and $\hat {\rho }=(\rho _{*}-\rho _{*}^{f})/\rho _{*}$ is the relative density difference between grains and interstitial fluid of density $\rho _{*}^{f}$. The constants $\hat {\mathcal {B}}=\mathcal {B}/\hat {\rho }\approx 0.713$, $\mathcal {C}=0.271$, $\mathcal {E}=2.096$ were obtained from a least squares fit of their experimental data (Trewhela et al. 2021a). The concentration-averaged particle diameter $\bar {d}$ gives the segregation velocity $f_{sl}$ its asymmetric behaviour (van der Vaart et al. 2015), since after rearranging and introducing the size ratio $R_{d}$ in its definition,

(2.7)\begin{equation} \bar{d}=(1-\chi_d\phi_{s})d_{l},\end{equation}

where $\chi _d=(R_{d}-1)/R_{d}\in [0,1]$ is the asymmetry parameter that shifts the maximum segregation velocities for small and large particles (Gray 2018). The definition for $\bar {d}$ is replaced in (2.6) to obtain

(2.8)\begin{equation} f_{sl}=\left(\frac{\rho_{*}g\mathcal{B}\dot{\gamma} (R_{d}-1)}{p_{w}+\rho_{*}\varPhi g z}\right) [1+\mathcal{E}(1-\phi_{s})(R_{d}-1)](1-\chi_d\phi_{s})^{2}d_{l}^{2},\end{equation}

where $p=p_w+\rho _{*}\varPhi g z$ replaced the pressure $p$ in (2.6). The term $\mathcal {C}\rho _{*}g\bar {d}$ was replaced by $p_{w}$, since it was introduced by Trewhela et al. (2021a) to avoid numerical singularities when $p=0$ at the free surface, a purpose that $p_{w}$ serves as well by definition.

From these definitions, the theoretical appearance of the segregation–rheology feedback is clear. The shear rate $\dot {\gamma }$ and pressure $p$ are decisive for the segregation velocity $f_{sl}$ magnitude, and reciprocally, the particle concentrations $\phi _{\nu }$ determine the frictional coefficient $\mu$, explicitly and implicitly via the averaged frictional coefficient $\bar {\mu }$ and the inertial number $I$, respectively (Golick & Daniels 2009; Chassagne et al. 2020; Bancroft & Johnson 2021; Trewhela et al. 2021a).

2.3. Non-dimensional coupled system of equations

Equation (2.8) synthesizes the interdependence between flow and particle distribution, hence a rheology and particle-size segregation feedback. By replacing the definitions made in §§ 2.1 and 2.2, and introducing the non-dimensional variables $\hat {t}=u_{w}t/h$, $\hat {z}=z/h$, $\hat {u}=u/u_{w}$, $\hat {p}=p/p_{w}$, $\hat {f}_{sl}=f_{sl}/u_{w}$ and $\hat {\mathcal {D}}=\mathcal {D}_{sl}/(hu_{w})$ into (2.1a–c) and (2.5a,b), the coupled partial differential equation (PDE) system for $\hat {u}$ and $\phi _{s}$, hence for $\phi _{l}=1-\phi _{s}$, can be obtained:

(2.9a,b)\begin{equation} \frac{\partial \hat{u}}{\partial \hat{t}}+\frac{\partial}{\partial \hat{z}}\left((1-\chi_{\mu}\phi_{s})R_{\mu}\mu_{s} \hat{p}\right)=0, \quad \frac{\partial \phi_{s}}{\partial \hat{t}}+\frac{\partial}{\partial \hat{z}}\left(\hat{f}_{sl}\phi_{s}(1-\phi_{s})-\hat{\mathcal{D}}\, \frac{\partial\phi_{s}}{\partial \hat{z}}\right)=0,\end{equation}

where the coupling with the velocity $\hat {u}$ in (2.9a) comes from the non-dimensional segregation velocity $\hat {f}_{sl}$ (see (2.8)) and particle diffusivity $\hat {\mathcal {D}}$, which are both shear rate $\dot {\gamma }$-dependent (Utter & Behringer 2004; Trewhela et al. 2021a). A non-dimensional segregation equation for $\phi _{l}$ is not included in the PDE system due to $\phi _{l}$ being the reciprocal of $\phi _{s}=1-\phi _{l}$. Because of their nonlinear nature, (2.9a,b) do not admit analytical solutions, but they can be solved separately numerically for the uncoupled cases (Cawthorn 2011; Jerome & Di Pierro 2018; Trewhela et al. 2021a).

2.4. Numerical method

The coupled PDE system in (2.9a,b) is solved using the method of lines (MOL) (Schiesser & Griffiths 2009). Schaeffer et al. (2019) used a similar numerical approach as a proof of concept to study compressibility in granular flows. In this work, the robustness of MOL was first tested for the uncoupled segregation equation in (2.5a,b), as validated by Trewhela et al. (2021a, figure 13) against the experimental results of van der Vaart et al. (2015). Second, the uncoupled monodisperse granular flow equation in (2.9a) was solved and compared to the solution presented by Cawthorn (2011, figure 2.10) and Jerome & Di Pierro (2018) (figure 2a) for $p_{w}=100$ Pa, $u_{w}=1$ m s$^{-1}$ and $\mu _{w}=0.5$, with the parameters in table 1. The plotted profiles for $\hat {u}$ in figure 2(b) also resemble those of Guillard et al. (2016, figures 2c,e), which were obtained using discrete numerical simulations for $g=0$ and $g=9.81$ m s$^{-2}$ and for a similar flow configuration. This resemblance validates the presented theoretical framework and numerical scheme.

Figure 2. (a) Velocity profile $\hat {u}$ solutions (2.9a) at different times $\hat {t}$ for the monodisperse case ($\phi _{s}=0$) and lithostatic pressure distribution with $p_{w}=100$ Pa, $u_{w}=1$ m s$^{-1}$ and $\mu _{w}=0.5$, similarly to the solutions found by Cawthorn (2011) and Jerome & Di Pierro (2018). (b,c) Comparison of the velocity profiles $\hat {u}$ and small particle concentration $\phi _{s}$ solutions for different times $\hat {t}$ for linear (lithostatic) $p=p_{w}+\rho g\varPhi z$ (continuous line) and constant pressure $p=p_{w}$ (dashed line) distributions under the same values for the wall parameters $u_{w}$, $p_{w}$ and $\mu _{w}$ with those of (a).

In this work, a particular pressure $p=p_{w}$ condition is used for the numerical solutions. A linear pressure (lithostatic) distribution adds further complexity to the numerical scheme with no relevant insights into the segregation–rheology feedback. Velocity profiles $\hat {u}$ for both pressure distributions are compared in figures 2(b,c) to see the influence of these distributions in the computations. Despite the fact that there are slight differences when comparing the solutions, these are not important compared to the cost of compromising the method's robustness. Pressure is still relevant for segregation (as seen in (2.8) and discussed briefly in Appendix A), but overall segregation dynamics does not change substantially when a constant pressure distribution is imposed over the lithostatic distribution. The latter assumption is particularly justified when considering a highly confined granular bulk, in which case $p_{w}\gg \rho _{*}g\varPhi h$. In contrast, the condition where $p_{w}$ is neglected in comparison to $\rho _{*}g\varPhi z$ poses difficulties to the numerical schemes since segregation velocity $f_{sl}\rightarrow \infty$ as $z\rightarrow 0$.

Figures 2(b,c) give a first glance at the different time scales at which the segregation–rheology feedback operates. A first time scale is that of momentum diffusion within the granular bulk, related to the granular flow rheology, and visible in the fast development of the velocity profile in figure 2(b). The other two time scales are longer than the former and are related to particle segregation and diffusion processes. These processes are visible in the slow changes of particle concentration profiles $\phi _{s}$ in figure 2(c).

3. Results

3.1. Effect of size ratio in the segregation and velocity profiles

To demonstrate the crucial role of particle sizes and segregation in the shear flow dynamics, the vertical distributions of small particles $\phi _{s}$ and horizontal velocity $\hat {u}$ are computed as functions of time $\hat {t}$ for three particles size ratios $R_{d}=\{1.25,2,3.33\}$, keeping $p_{w}=100$ Pa, $\mu _{w}=0.5$ and $u_{w}=0.2$ m s$^{-1}$ constant (figure 3). The initially unstable arrangement of small particles over large ones influences the velocity $\hat {u}$ profiles; all reach a transient state at $\hat {t}\approx 10^{2}$, only to be then altered by segregation. From $\hat {t}>10^{2}$, the linear velocity profiles $\hat {u}$ for both species, and featuring a slope change at the small–large particle interface, changes into a nonlinear smooth distribution for the final steady-state condition. At initial times, segregation does not influence the velocity profile $\hat {u}$ or the time scale to attain the transient linear state. Therefore, momentum transfer is uncoupled from segregation, with the latter acting at longer time scales, an observation made in various experiments and flow configurations (Golick & Daniels 2009; Tripathi & Khakhar 2011; van der Vaart et al. 2015; Trewhela & Ancey 2021). Nonetheless, size ratio $R_{d}$ influences the transient linear state, determining the velocity vertical gradient, a result also observed by Tripathi & Khakhar (2011) and corroborated by Barker et al. (2017). The asymmetric nature of segregation is captured clearly in figures 3(a–c). As the size ratio $R_{d}$ increases, small particles segregate faster to the bottom compared to large particles arriving at the top. As a result, when asymmetry $\chi _{d}$ is increased, segregation fluxes are fast and dynamic, achieving the nonlinear steady-state swiftly (as observed for $R_{d}=3.33$, $\chi _{d}=0.7$ in figures 3c,f).

Table 2. Definition of the parameters used to set the presented numerical simulations in this work. The numerical values for the size asymmetry coefficient $\chi _{d}$ (hence particle-size ratio $R_{d}$) and the frictional asymmetry coefficient $\chi _{\mu }$ are used to set the particle parameters. The numerical values for pressure $p_{w}$, velocity $u_{w}$ and frictional coefficient $\mu _{w}$ at the upper plate or wall (see figure 1) resulted in a range of values for $h$ (see (2.4)) and set the granular rheology, together with the values of table 1. Finally, the empirical parameters $\mathcal {A}$, $\mathcal {B}$, $\mathcal {E}$ proposed by Trewhela et al. (2021a) and Utter & Behringer (2004) control the segregation dynamics in the numerical solutions.

Figure 3. (a–c) Vertical and temporal distributions of small particle concentrations $\phi _{s}$ for $u_{w}=0.2$ m s$^{-1}$, $p_{w}=100$ Pa, $\mu _{w}=0.5$ and $R_{d}=\{1.25,2,3.33\}$ values. (d–f) Small particle concentration profiles $\phi _{s}$ at different $\hat {t}$ times for the corresponding $R_{d}$ values. (g–i) Velocity profiles $\hat {u}$ at different $\hat {t}$ instants for the corresponding $R_{d}$ values, reflected in its asymmetrical coefficient $\chi _{d}=(R_{d}-1)/R_{d}$ value. The thicker line in each of (d–i) highlights the final solution of each profile. Animated solutions of (d–i) are in supplementary movie 1, available at https://doi.org/10.1017/jfm.2024.168.

Other striking differences appear when comparing the final small particle concentrations $\phi _{s}$. Larger $R_{d}$ values produce sharper profiles with strongly segregated $\phi _{s}$ profiles and a marked small–large interface. Subsequently, this segregated and nonlinear distribution is imprinted into the shape of the steady-state velocity profile $\hat {u}$, thus a direct result of particle-size segregation effectiveness. For size ratios $R_{d}\rightarrow 1$, apparent mixing (i.e. diffusive remixing) is observed. Controlled by diffusion, mild segregation produces quasi-linear $\hat {u}$ profiles, as shown in figure 3(d). In contrast to these quasi-linear profiles, $\hat {u}$ solutions for large $R_{d}$ show a sharp break in the velocity profile, drawing a rhomboid or lens geometry. The rhomboid is defined graphically by the linear responses expected for the normally graded and inversely graded conditions for particle-size distribution.

3.2. Initially mixed conditions and resulting profiles

The normally graded condition shown in figure 1(a) is a widely explored initial condition for segregation in experiments and numerical simulations (e.g. van der Vaart et al. 2015; Gray 2018; Chassagne et al. 2020). However, in the context of mixed bidisperse granular flows, the initially mixed condition illustrated in figure 1(b) also provides a valuable assessment on how segregation affects the rheology of granular mixtures. Figure 4 shows the resulting small particle concentration $\phi _{s}$ and velocity $\hat {u}$ profiles for three different initially mixed conditions of small and large particles at $\phi _{s}^{0}=25$, 50 and 75 % homogeneously mixed concentrations. To compare with the results shown previously in figures 3(c,f,i), an $R_{d}=3.33$ value is used. The steady state $\phi _{s}$ distribution and the $\hat {u}$ profile observed in figures 4(b,e,h) are identical to those observed in figures 3(c,f,i). This agreement proves that steady-state solutions are independent of the normally graded or mixed initial conditions. Such independence is just about the particles’ initial arrangement, and it is not extended to the average particle concentration $\phi _{s}^{0}$, a condition that results in differing profiles for $\phi _{s}$ and $\hat {u}$ (figure 4).

Figure 4. (a–c) Vertical and temporal distributions of small particle concentrations $\phi _{s}$ for $u_{w}=0.2$ m s$^{-1}$, $p_{w}=100$ Pa, $\mu _{w}=0.5$ and $R_{d}=3.33$ for initially-mixed conditions of small and large particles, $\phi ^{0}_{s}=\{0.25,0.5,0.75\}$. (d–f) Small particle concentration $\phi _{s}$ and (g–i) velocity $\hat {u}$ profiles at different $\hat {t}$ instants for the corresponding $\phi _{s}^{0}$ and the same wall parameters and $R_{d}$ value. The thicker line in each of (d–i) highlights the final solution of each profile. Animated solutions of (d–i) are available in supplementary movie 2.

As expected, there are important differences in the intermediate solutions for both $\phi _{s}$ and $\hat {u}$ due to the change in the initial conditions. Homogeneously mixed initial $\phi _{s}^{0}$ concentrations produce not only different velocity $\hat {u}$ and skewed – but similar – $\phi _{s}$ profiles, but also delay segregation. The delayed steady-state solutions observed in figures 4(a–c) for larger $\phi _{s}^{0}$ concentrations are tied inherently to the asymmetric nature of size segregation, i.e. the segregation flux $\boldsymbol {F}_{s}=f_{sl}\phi _{s}\phi _{l}\boldsymbol {g}/|\boldsymbol {g}|\sim (1-\chi _{d}\phi _{s})^{2}\phi _{s}(1-\phi _{s})$. Small particles at low $\phi _{s}$ concentrations segregate at a faster rate than large particles at low $\phi _{l}$ concentrations (Gajjar & Gray 2014; van der Vaart et al. 2015; Trewhela et al. 2021a). If segregation fluxes were to be symmetric $\boldsymbol {F}_{s}= \phi _{s}\phi _{l}\boldsymbol {g}/|\boldsymbol {g}|$, then segregation rate would be maximum at $\phi _{s}^{0}=0.5$, and at identical rates for $\phi _{s}^{0}=\{0.25,0.75\}$ (Gray 2018). Therefore, figure 4 reproduce the expected asymmetry. For $\phi _{s}^{0}=0.25$, segregation alters the velocity profile solution swiftly from the mixed condition (at $\hat {t}\approx 10^{2}$) into the steady-state condition at almost $\hat {t}\approx 10^{4}$. In the case of $\phi _{s}^{0}=0.75$, segregation also initiates at $\hat {t}\approx 10^{2}$, but steady-state solutions are reached at a slower pace, i.e. at $\hat {t}\approx 10^{5}$. These differences hint at a $\phi _{s}$-dependent time scale for segregation that controls the final duration of the segregation–rheology feedback.

3.3. Influence of differences in the particles’ frictional coefficients

An additional effect in the Stokes’ problem is explored by changing the small particles’ frictional coefficient $\mu _{s}$ while keeping the size ratio $R_{d}=2$ fixed, along with $p_{w}=300$ Pa, $\mu _{w}=0.5$ and $u_w=0.2$ m s$^{-1}$. Differences in $\mu _{s}$, relative to $\mu _{l}$ via $R_{\mu }$, introduce a non-zero $\chi _{\mu }$ value in (2.9a), which is null when small and large particles share the same $\mu (I)$-rheology parameters (table 1). Changes in $\chi _{\mu }$ lead to velocity differences at the small–large particles interface, compromising the numerical scheme due to $\hat {u}$ discontinuity. Numerical solutions are found only for slight changes in $\chi _{\mu }$, within the interval $\chi _{\mu }=[-0.033,0.013]$ for $R_{d}=2$. This interval depends on the size ratio as well; as the size ratio $R_{d}$ increases, the $\chi _{\mu }$ interval for the solutions narrows.

In figure 5, $\phi _{s}$ and $\hat {u}$ solutions are plotted for various $\chi _\mu$ values. While the small particle concentration $\phi _{s}$ solutions do not show significant differences between the different $\chi _\mu$ cases, the velocity profiles $\hat {u}$ display interesting behaviour. For negative $\chi _\mu$ values, $\mu _{s}$ is larger than $\mu _{l}$, resulting in the transient solution having inverted slopes compared to the case with positive $\chi _\mu$ values; see figures 5(f,j) for comparison. Intermediate negative values for $\chi _\mu$ seem to counteract the $\hat {u}$ profile nonlinearity generated by the size difference in $R_{d}$ (figures 5g,h). This finding is consistent with the proposed theoretical framework and could illustrate a way to obtain linear velocity profiles in granular flows with particles of varying sizes and frictional properties. Therefore, to control the particles’ frictional properties is to partially control the segregation–rheology feedback, despite having a segregated bulk with different grain sizes.

Figure 5. Numerical solutions for (a–e) small particle concentrations $\phi _{s}$, and (f–j) velocity profiles $\hat {u}$ with constant wall parameters $u_{w}=0.2$ m s$^{-1}$, $p_w=300$ Pa, $\mu _{w}=0.5$, $R_{d}=2$, and variable $\chi _{\mu }=\{-0.03,-0.02,-0.01,0,0.01\}$ for (a,f), (b,g), (c,h), (d,i) and (e,j), respectively. The thicker line in each plot highlights the final solution of each variable. Solutions plotted in (d,f–j) are shown as animations in supplementary movie 3.

3.4. Steady-state analytical solutions and comparison with numerical solutions

A solution for the transient velocity profiles $\hat {u}$ can be determined analytically from the steady-state condition in (2.9a), i.e. $\bar {\mu }=\mu _{w}$. This condition is replaced in the quadratic expression for $I$, obtained from the regularized $\mu (I)$-rheology in (2.2):

(3.1)\begin{equation} \mu_{\infty}I^{2}+(\mu_{2}-\mu_{w})I+(\mu_{1}-\mu_{w})I_{0}=0.\end{equation}

The positive solution for $I$ is then balanced with the inertial number $I$ definition to yield a differential equation for the velocity $\hat {u}$. This ODE is solved by separating variables and using the no-slip boundary condition $\hat {u}|_{\hat {z}=1}=0$ to yield

(3.2)\begin{equation} \hat{u}=\frac{h}{u_{w}\bar{d}}\,\sqrt{\frac{p_{w}}{\rho_{*}}}\, \mathcal{I}(\mu_{w})\,(1-\hat{z}),\end{equation}

where $\mathcal {I}(\mu _{w})=((\mu _{w}-\mu _{2})+\sqrt {(\mu _{w}- \mu _{2})^{2}+4\mu _{\infty }(\mu _{w}-\mu _{1})I_{0}})/2\mu _{\infty }$ is the aforementioned positive solution of (3.1), i.e. the inertial number $I$. Equation (3.2) results in the rhomboid formation, already described and discussed in § 3.1. If the concentration-averaged diameter $\bar {d}$ is evaluated with the normally-graded size distribution, $\phi _{s}\big |_{\hat {z}=0}^{\hat {z}=1/2}=1$ and $\phi _{s}\big |_{\hat {z}=1/2}^{\hat {z}=1}=0$, and the inversely-graded size distribution, $\phi _{s}\big |_{\hat {z}=0}^{\hat {z}=1/2}=0$ and $\phi _{s}\big |_{\hat {z}=1/2}^{\hat {z}=1}=1$, then the rhomboid can be diagrammed. The solution captured by (3.2) provides a good approximation of the transient velocity profile without considering segregation in its development. When segregation determines the velocity profile, a steady-state solution for the small particle concentration $\phi _{s}$ can be obtained. For that, the procedure used by Trewhela et al. (2021a) is followed, yet their solution is implicit for $\phi _{s}$. Here, a simpler differential equation is developed from $\mathcal {E}=0$ (a first-order approximation in $R_{d}-1$), resulting in the explicit solution

(3.3)\begin{equation} \phi_{s}=\frac{{\rm e}^{(\hat{z}+\mathcal{K})/\lambda}}{1+ {\rm e}^{(\hat{z}+\mathcal{K})/\lambda}},\end{equation}

where $\mathcal {K}=-1/2$ is the constant of integration, which can be calculated using the depth-averaged condition $\bar {\phi }_{s}=\int ^{1}_{0}\phi _{s}\,\text {d}\hat {z}=0.5$, and $\lambda =\mathcal {A}p_{w}/(\mathcal {B}\rho g h(R_{d}-1))$. This small particle concentration $\phi _{s}$ profile serves as input for the generalized ODE:

(3.4)\begin{equation} \frac{\text{d}\hat{u}}{\text{d}\hat{z}}=-\frac{h}{u_{w}\bar{d}}\,\sqrt{\frac{p_{w}}{\rho_{*}}}\,\mathcal{I}(\hat{\mu}_{w}), \end{equation}

where $\hat {\mu }_{w}=\mu _{w}(1-2.7(1-\phi _{s})\chi _{\mu })/(1-\chi _{\mu }\phi _{s})$ is a $\phi _{s}$- and $\chi _{\mu }$-dependent frictional coefficient, with $\bar {d}$ a function of the former as well. Due to this dependence, the ODE's solution is computed using a Runge–Kutta fourth-order method.

Figure 6 summarizes the core behaviour of the segregation–rheology feedback. The transient-state solution of (3.2) and its reciprocal stable state (red and blue toned rhomboids) are shown in figure 6(a). These rhomboids emphasize the slope change in $\hat {u}$ as $R_{d}$ increases, with the fully linear case $R_{d}=1$ as reference. In figure 6(b), the generalized solution in (3.4) is drawn for the $R_{d}=2$ case and various $\chi _{\mu }$ values, with the reference curve $\chi _{\mu }=0$ at the middle. Negative $\chi _{\mu }$ values tend to displace the steady-state velocity profile $\hat {u}$ towards the linear reference, acting as a counterpart for the deviation generated by the size ratio $R_{d}$. Conversely, positive $\chi _{\mu }$ values exacerbate the nonlinear behaviour of the velocity profiles $\hat {u}$, displacing the curves beyond the reciprocal stable-state curve (with large particles over small ones). Moreover, solutions attempting to pass beyond the stable state could explain the difficulties in obtaining solutions when the frictional asymmetric coefficient is $\chi _\mu >0.013$ for $R_{d}=2$. A detailed evolution of the velocity profile $\hat {u}$ solutions is shown in figure 6(c). The velocity profile $\hat {u}$ develops from the initial condition to the transient state $\hat {t}\sim 10^{2}$, adjusting well to the analytical solution of (3.2) (left rhomboid linear curves in figure 6(c), taken from $R_{d}=2$ in figures 6a,b). Then segregation alters this state, rearranging particles according to their size, achieving a steady state at $\hat {t}\approx 10^{5}$ that fits well the ODE solution in (3.4) (dot-dashed line in figure 6c). In figure 6(c) the segregation–rheology feedback ends up being described graphically by the changes in $\hat {u}$ from the transient to the steady state, hence analytically bounded within the $R_{d}$-dependent rhomboid region.

Figure 6. (a) Velocity profiles $\hat {u}$ for the steady-state condition in (3.2) for various $R_{d}=[1,4]$ values, represented in blue and red tones for $d_{s}$ and $d_{l}$, respectively. (b) Velocity profiles $\hat {u}$ of the numerically solved (3.4) for variable $\chi _{\mu }=[-0.05,0.05]$ and constant $R_{d}=2$, with its corresponding rhomboid from (a). (c) Temporal evolution of the velocity profile $\hat {u}$ for the $R_{d}=2$ and $\chi _{\mu }=0.01$ case compared with the corresponding steady-state solutions from (a) ((3.2), plotted as a rhomboid) and (b) ((3.4), plotted with a dot-dashed line). Note how the velocity profiles $\hat {u}$ first collapse around the rhomboid solution, and then converge to the steady-state solution (red dot-dashed line). Fully animated solutions (for $\chi _{\mu }=\{-0.03,-0.01,0.01\}$ and constant $R_{d}=2$), compared to the steady-state solutions of (3.2) and (3.4), are shown in supplementary movies 4–6.

3.5. Granular kinematic viscosity, time scales and segregation–momentum balances

From the ODE for the non-dimensional shear rate $\hat {\dot {\gamma }}=|\text {d}\hat {u}/\text {d}\hat {z}|$ in (3.4), a general definition for the granular kinematic viscosity $\nu _{g}$ can be formulated. Using the theoretical constitutive relation for a granular dynamic viscosity $\tau =\mu p=\eta _{g}\dot {\gamma }$ (Jing et al. 2022), the dimensional kinematic viscosity can be defined as

(3.5)\begin{equation} \nu_{g}=\frac{\hat{\mu}_{w}\bar{d}}{\mathcal{I}(\hat{\mu}_{w})}\, \sqrt{\frac{p_{w}}{\rho_{*}}}.\end{equation}

Compared to the granular viscosity definition of Jerome & Di Pierro (2018), in (3.5) this bidisperse granular viscosity is $R_{d}$-, $\chi _{\mu }$- and $\phi _{s}$-dependent (via $\bar {d}$ and $\hat {\mu }_{w}$), and is still determined by the wall conditions and the $\mu (I)$-rheology parameters (table 1). It is worth recalling that the $\mathcal {I}$ function is in fact the inertial number $I$, but for the Stokes’ problem analysed in this work, it is left as the function defined in (3.2).

The definition of $\nu _{g}$ allows the analysis of the different time scales at which the segregation–rheology feedback develops. Evidently, a first time scale appears from the non-dimensional time $\hat {t}$ definition, where $t_{\dot {\gamma }}\sim h/u_{w}$ can be defined as a convective or shear rate time scale, inferred from $1/\dot {\gamma }=(\text {d}u/\text {d}z)^{-1}\approx h/u_{w}$ (as in the definition of $I$; GDR-MiDi 2004). A second time scale can be associated with momentum diffusing through the granular bulk and derived directly from (3.5), yielding

(3.6)\begin{equation} t_{\nu_{g}}\sim \frac{h^{2}}{\nu_{g}}=\frac{h^{2}\,\mathcal{I}(\hat{\mu}_{w})}{\hat{\mu}_{w} \bar{d}}\,\sqrt{\frac{\rho_{*}}{p_{w}}}.\end{equation}

Then two more time scales appear from the particle-segregation-related processes. A segregation time scale is drawn by the time it takes particles of diameter $\bar {d}$ to segregate through the sheared layer $h$ with velocity $f_{sl}$:

(3.7)\begin{equation} t_{f_{sl}}\sim \frac{h}{f_{sl}}=n_{\bar{d}}P\,\mathcal{F}(R_{d},\phi_{s})^{-1}\,t_{\dot{\gamma}}, \end{equation}

where $n_{\bar {d}}=h/\bar {d}$ is the number of particle layers through the sheared layer depth, $\mathcal {F}(R_{d},\phi _{s})=\mathcal {B}(R_{d}-1)[1+\mathcal {E}(R_{d}-1)(1-\phi _{s})]$ is a function defined from the segregation velocity $f_{sl}$ function, and $P=p_{w}/(\rho _{*}g\bar {d})$ corresponds to a non-dimensional pressure (Trewhela et al. 2021a). The segregation time scale is proportional to $t_{\dot {\gamma }}$ and $p_{w}$, in agreement with the clear segregation dependence on shear rate and pressure (Golick & Daniels 2009; Trewhela et al. 2021b). Thus the segregation time scale $t_{f_{sl}}$ is the shear rate time scale $t_{\dot {\gamma }}$ modulated by pressure conditions at the wall, segregation, and the number of layers that particles must segregate through. Finally, the particle diffusion time scale can be defined as

(3.8)\begin{equation} t_{\mathcal{D}_{sl}}\sim \frac{h^{2}}{\mathcal{D}_{sl}}= \frac{n_{\bar{d}}^{2}}{\mathcal{A}}\,t_{\dot{\gamma}}, \end{equation}

which again is proportional to $t_{\dot {\gamma }}$. Although these time scales do not provide the actual time that the velocity $\hat {u}$ and small particle concentration $\phi _{s}$ profiles take to reach their corresponding steady-state solutions, which is out of the scope of this work, they serve to describe the progression of the segregation–rheology feedback.

For the numerical solutions computed in § 3.1, $t_{\dot {\gamma }}=0.0157$ s and the momentum diffusion time scale share similar values $t_{\nu _{g}}\approx 0.15$ s, yielding a non-dimensional time $\hat {t}=t_{\nu _{g}}/t_{\dot {\gamma }}\approx 10$ for the $\hat {u}$ profiles to develop into the transient solution. This value is close to what is shown in figure 3, with $\hat {t}\approx 10^2$ being the time at which the transient solution is achieved in the computations. In terms of segregation, there are larger differences in $t_{f_{sl}}\approx \{100,20,7\}$ s and $t_{\mathcal {D}_{sl}}=\{80,120,170\}$ s for $R_{d}=\{1.25,2,3.33\}$, respectively. These $t_{f_{sl}}$ values capture well the faster set in of segregation at $\hat {t}=t_{f_{sl}}/t_{\dot {\gamma }} \approx \{6\times 10^{3},10^{3},5\times 10^{2}\}$ for the same $R_{d}$ values. However, these $t_{f_{sl}}$ are balanced out by shorter or larger diffusive time scales $t_{\mathcal {D}_{sl}}$, which finally settle the concentration profiles into their steady-state solution. Consistent results are observed when analysing these time scales in the case of the initially mixed condition solutions at fixed $R_{d}=3.33$. These yield the same $t_{\dot {\gamma }}$ and $t_{\nu _{g}}$ values for all solutions, but $t_{f_{sl}}\approx \{3,7,20\}$ s and $t_{\mathcal {D}_{sl}}\approx \{100,170,310\}$ s increase with $\phi _{s}^{0}=\{0.25,0.5,0.75\}$, resulting in the monotonic increase of $\hat {t}=t_{\mathcal {D}_{sl}}/t_{\dot {\gamma }} \approx \{7\times 10^3,10^4,2\times 10^4\}$. These values are also similar to the times found in the numerical solutions (figure 4). In the case of the varying $\chi _\mu$ solutions, there are only slight differences in $t_{\nu _{g}}$. Values decrease from approximately 0.9 s to 0.7 s as $\chi _\mu$ goes from $-$0.03 to 0.01.

From these results and the analysed time scales, it is thus natural to compare the two diffusion processes acting in the bidisperse granular Stokes’ problem. With the definition of the granular momentum diffusivity time scale, or its diffusivity coefficient in $\nu _{g}$, the Schmidt number or the ratio of the momentum and particle diffusivity time scales is obtained directly as

(3.9)\begin{equation} Sc=\frac{t_{\nu_{g}}}{t_{\mathcal{D}_{sl}}}=\frac{\nu_{g}}{\mathcal{D}_{sl}}= \frac{\hat{\mu}_{w}}{\mathcal{A}I^{2}}, \end{equation}

where the inertial number is reintroduced for simplicity and generalization; for that matter, $I$ can be considered equal to $d\dot {\gamma }/\sqrt {p/\rho _*}$ or equal to the definition of $\mathcal {I}$ made after (3.2). It is worth noting that the $I^{2}$ term corresponds to the Savage or Coulomb number (Savage 1984; Ancey, Coussot & Evesque 2000). Low $Sc$ values indicate an important particle diffusivity (diffusive remixing), dominating over momentum diffusion, which at first glance can be achieved with low $\mu _{w}$, high $\mathcal {A}$ or high inertial numbers $I$. Conversely, a large $Sc$ points towards a more viscous response from the granular material, with a lesser role played by particle diffusion in the overall flow dynamics. Interestingly, this definition of the Schmidt number depends implicitly not only on the particles’ concentration-averaged diameter $\bar {d}$ (via $I$), hence their size ratio $R_{d}$ and size asymmetry coefficient, but also on the frictional asymmetry coefficient $\chi _{\mu }$ and the particles’ frictional parameters.

To illustrate these dependencies, figure 7(a) shows the Schmidt number $Sc$ as a function of both $\mu _{w}$ and an overall small particle concentration, which can be interpreted as the initial small particle concentration $\phi _{s}^{0}$ or the depth-averaged small particle concentration $\overline {\phi _{s}}=\int _{0}^{1}\phi _{s}\,\text {d}\hat {z}$. In figure 7(a), the contours’ slope is determined by $\chi _{\mu }$: positive $\chi _{\mu }$ yields a negative slope (dashed contours), and vice versa (filled contours). For the computed values in the domain $\mu _{w}=[0.4,1.2]$ and $\phi _{s}=[0,1]$, $Sc$ varies from $10^{3}$ to $10^{-1}$, decreasing linearly in two regions marked approximately by $\mu _{w}\lessgtr \mu _{2}=0.557$. This is a direct result of the partially regularized $\mu (I)$-rheology, which allows $\mu _{w}>\mu _{2}$, therefore extending the range for $Sc$. Overall, $Sc$ decreases with $\mu _{w}$, and increases slightly with higher $\phi _{s}^{0}$, when $\chi _{\mu }<0$. The former observation contradicts the linear dependence of $Sc$ on $\mu _{w}$, but is consistent with the fact that high $\mu _{w}$ values result in high $I$ values as well, since $Sc\sim I^{-2}$. (For more results shedding light on the role of $\mu _{w}$, see Appendix A.) While particle diffusivity scales with the shear rate $\mathcal {D}_{sl}\sim \dot {\gamma }$ (Utter & Behringer 2004), the granular viscosity scales inversely with the latter $\nu _{g}\sim \dot {\gamma }^{-1}$; so larger friction at the wall produce high shear rates and will result in low $Sc$. But as stated before, the segregation–rheology feedback is not a purely momentum- or particle-diffusion-dominated process; segregation also controls the flow dynamics.

Figure 7. (a) Filled contours for the Schmidt number $Sc=t_{\nu _{g}}/t_{\mathcal {D}_{sl}}=\nu _{g}/\mathcal {D}_{sl}$ as a function of both $\mu _{w}$ and the overall small particle concentration $\phi _{s}^{0}$, in the case $\chi _{\mu }=-0.03$. Dashed contours show the change of slope in the case $\chi _{\mu }=0.03$ for the same $Sc$ values of the filled contours. (b) Filled contours for the segregation–rheology Péclet number $Pe_{sr}=t_{f_{sl}}/t_{\nu _{g}}=hf_{sl}/\nu _{g}$ as a function of $\mu _{w}$ and the size ratio $R_{d}$, in the case $\chi _{\mu }=0.01$. (c) Variation of the segregation–rheology Péclet number $Pe_{sr}$ as a function of the wall friction $\mu _{w}$ and the frictional asymmetry parameter $\chi _{\mu }$. For comparison, the constant frictional parameter $\mu _{2}=0.557$ is plotted alongside the same $Pe_{sr}$, but using the non-regularized inertial number definition $\mathcal {I}_{nr}=I_{0}(\mu _{w}-\mu _{1})/(\mu _{2}-\mu _{w})$, which is independent of $\chi _{\mu }$ (dot-dashed thick line).

To understand better at which extent segregation controls the segregation–rheology feedback, a particle segregation–diffusion balance is borrowed from previous work. Following the definitions made by Gray (2018) and Trewhela et al. (2021a), a Péclet number for segregation $Pe$ can be defined as the ratio between the segregation velocity $f_{sl}$ and particles’ diffusivity $\mathcal {D}_{sl}$, i.e. $Pe=hf_{sl}/\mathcal {D}_{sl}$, or alternatively, $t_{f_{sl}}/t_{\mathcal {D}_{sl}}$. With this definition, a specific Péclet number for the segregation–rheology feedback can be introduced as

(3.10)\begin{equation} Pe_{sr}=\frac{Pe}{Sc}=\frac{hf_{sl}}{\nu_{g}}= \frac{t_{f_{sl}}}{t_{\nu_{g}}}=\frac{I^{2}}{\varPhi_{p} \hat{\mu}_{w}}\,\mathcal{F}(R_{d},\phi_{s}), \end{equation}

where $\varPhi _{p}=p_{w}/(\rho g h)$ is a non-dimensional pressure, which can be interpreted as an equivalent solids volume fraction. Again, this second non-dimensional quantity depends on $I^{2}$ (or the Savage number). Compared to the previous $Sc$ definition, the segregation–rheology Péclet number $Pe_{sr}$ depends on the size ratio $R_{d}$ and the empirical constants determined by Trewhela et al. (2021a). The proposed $Pe_{sr}$ somehow confirms the results discussed previously in §§ 3.1 and 3.3, where segregation profiles were sharper for large $R_{d}$ (figure 3) and positive $\chi _{\mu }$ values (figure 5).

Figures 7(b,c) show the variation of the segregation–rheology Péclet number $Pe_{sr}$ as a function of $\hat {\mu }_{w}$, $\mu _{w}$ and $R_{d}$. Broadly, the $Pe_{sr}$ dependence on $R_{d}$ is well captured in figure 7(b). As expected, $Pe_{sr}$ is minimum for $R_{d}\rightarrow 1$, and interestingly, larger size differences do not necessarily result in predominant segregation effects over momentum diffusion, since $Pe_{sr}=0.01$ for $R_{d}=4$ and $\mu _{w}=0.4$. Therefore, the results presented in the previous subsections are summarized conceptually in the $Pe_{sr}$ definition; the segregation–rheology feedback is relevant when a combination of both high friction and large size differences is in place. Nonetheless, high $\mu _{w}$ is still enough to promote segregation, even accompanied by low $R_{d}$. All this, as long as there is sufficient pressure $p_{w}$ or the bulk is compressed ($\varPhi _{p}$) enough at the wall to transfer shear without restraining both momentum and segregation. The latter opens the possibility of including compressibility effects in the proposed system and solutions. In figure 7(c), the $Pe_{sr}$ curves are plotted for various $\chi _{\mu }$ values, all of which show similar behaviour. However, when $Pe_{sr}$ is calculated using the non-regularized $\mu (I)$-rheology, i.e. $\mathcal {I}_{nr}=I_{0}(\mu _{w}-\mu _{1})/(\mu _{1}-\mu _{w})$, the $Pe_{sr}$ value diverges as $\mu _{w}\rightarrow \mu _{2}$ (thick dashed line in figure 7c). This shows the importance of having a well-posed extended formulation for the $\mu (I)$-rheology. When using the partially regularized $\mu (I)$-rheology, $Pe_{sr}$ continues to increase monotonically with $\mu _{w}>\mu _{2}$, but at a lower rate than for $\mu _{w}<\mu _{2}$, and without a clear cutoff value.

4. Conclusions

The segregation–rheology feedback in dense granular flows was explored numerically and theoretically via Stokes’ first problem for a bidisperse granular bulk at constant solids volume fraction. Using this paradigmatic problem for such a nonlinear and coupled system required the implementation of the method of lines, which proved to be a practical numerical scheme despite the issues inherited from the $\mu (I)$-rheology implementation. The solved PDE system incorporated the partially regularized $\mu (I)$-rheology of Barker & Gray (2017) and the scaling law for bidisperse size segregation of Trewhela et al. (2021a). The results focused on the role of size ratio $R_{d}$ in the development of velocity profiles, and included the effect of frictional differences or a frictional ratio $R_{\mu }$ between particle species via the introduced parameter $\chi _{\mu }=(R_{\mu }-1)/R_{\mu }$, in analogy with the segregation asymmetry coefficient $\chi _{d}=(R_{d}-1)/R_{d}$. Under a constant pressure $p_{w}$ condition, the numerical solutions showed a transient state, decoupled from the segregation process, achieving linear velocity profiles $\hat {u}$. Then segregation-induced feedback alters particle concentrations, thus changing the linear $\hat {u}$ state into a $\phi _{s}$- and $R_{d}$-dependent nonlinear solution for $\hat {u}$. The transient and steady-state solutions encompass a rhomboid region that expands as $R_{d}$ increases, and reduces to a single linear profile when $R_{d}\rightarrow 1$ (with no segregation). The steady-state solution is determined not only by $R_{d}$ but also by the particles’ frictional coefficient ratio $R_{\mu }$. While positive $\chi _\mu$ values exacerbate the solutions’ nonlinearity, negative values of $\chi _{d}$ dampen the effect of $R_{d}$ and favour a linear behaviour, counteracting the segregation-induced feedback. The numerical solutions converge well around the analytical expressions (3.2) and (3.4), which bound the feedback's solution range.

The kinematic granular viscosity $\nu _{g}$ was redefined based on the partially regularized $\mu (I)$-rheology and the concentration-averaged diameter $\bar {d}$, hence determined locally by particle concentrations $\phi _{\nu }$. Four characteristic time scales arising from shear rate $t_{\dot {\gamma }}$, granular momentum diffusion $t_{\nu _{g}}$, segregation $t_{f_{sl}}$ and particle diffusion $t_{\mathcal {D}_{sl}}$ were introduced to describe the progression of the segregation–rheology feedback. These time scales did not provide precise expressions for the time at which the variables attained the transient and steady states, but explained well the observed temporal differences in the presented numerical results. From $\nu _{g}$ and these four time scales, two non-dimensional numbers appeared naturally as a result of segregation–momentum balances: a Schmidt number $Sc=\nu _{g}/\mathcal {D}_{sl}$, and a specific segregation–rheology Péclet number $Pe_{sr}=hf_{sl}/\nu _{g}$. These two quantities not only showed consistency with the numerical and theoretical results, but also provided a broad description for the segregation–rheology feedback and its dependence on the plethora of parameters addressed here.

Although the coupled Stokes’ problem is an idealized case study, the developed analyses provide a robust numerical scheme for particle-size segregation in granular flows that is easy to compare with experimental results and applicable to the rheometry of dense granular flows (e.g. Golick & Daniels 2009; Trewhela & Ancey 2021). Particle-size segregation and local friction redistribution, as discussed in this work, may have consequences for the scales associated with shear bands (May et al. 2010), or other processes like jamming (Rietz & Stannarius 2008). Further research is required to extend this feedback analysis towards an extended range of size ratios and flow characteristics – compressibility, reverse segregation or other segregation mechanisms – that are still not considered but could complement or fit well in the presented approach (Schaeffer et al. 2019; Neveu et al. 2022).