2.1. Bidisperse systems

We consider granular mixtures consisting of particles with two sizes – large grains with an average diameter of $d^{l}$ and small grains with an average diameter of $d^{s}$. We consider both two-dimensional systems of disks, as illustrated in figure 1, and three-dimensional systems of spheres. To eliminate the effect of density-based segregation (e.g. Tripathi & Khakhar Reference Tripathi and Khakhar2013) and isolate size-based segregation, all particles are made of the same material with mass density $\rho _{s}$, which represents the area density for disks and the volume density for spheres. Throughout, we utilize the notational convention in which we denote large-grain quantities using a superscript ${l}$ and small-grain quantities using a superscript ${s}$. The species-specific solid fractions – i.e. the areas occupied by each species per unit total area for disks and the volumes occupied by each species per unit total volume for spheres – are $\phi ^{l}$ and $\phi ^{s}$, respectively, and the total solid fraction is $\phi = \phi ^{l} + \phi ^{s}$. The concentration of each species then follows as $c^{l} = \phi ^{l}/\phi$ and $c^{s} = \phi ^{s}/\phi$, so that $c^{l} + c^{s} = 1$. The average mixture grain size is defined as the sizes of both species weighted by their concentrations, $\bar {d} = c^{l}d^{l} + c^{s}d^{s}$. We make the common idealization that the total area for dense systems of disks or total volume for dense systems of spheres does not change (Savage Reference Savage1998; Gray & Thornton Reference Gray and Thornton2005; Gray & Chugunov Reference Gray and Chugunov2006; Fan & Hill Reference Fan and Hill2011b), and therefore $\phi$ is idealized as constant at each point in space and at each instant in time during the segregation process. We have verified in our DEM simulations that area (or volume) dilatation at flow initiation occurs over a much shorter time scale than the process of segregation, and that, subsequently, the solid fraction field is approximately uniform both spatially and temporally. Moreover, we have verified that dilatancy-driven secondary flows (e.g. Dsouza & Nott Reference Dsouza and Nott2021) are not observed in the dense flows considered in this work. Therefore, the idealization of a constant solid fraction is reasonable. Throughout this study, we use $\phi =0.8$ for disks, and $\phi =0.6$ for spheres.

Figure 1. A representative schematic of a dense, bidisperse granular system consisting of two-dimensional disks.

Regarding the kinematics of flow, each species has an associated partial velocity, $v_i^{l}$ and $v_i^{s}$, and the mixture velocity is given by $v_i = c^{l}v_i^{l} + c^{s}v_i^{s}$. The mixture strain-rate tensor is then defined using the mixture velocity in the standard way: ${\mathsf{D}}_{ij} = (1/2)(\partial v_i/\partial x_j + \partial v_j/\partial x_i)$, where ${\mathsf{D}}_{kk}=0$ since we have assumed that the mixture area (or volume) does not change. The equivalent shear strain rate is defined as $\dot {\gamma } = (2{\mathsf{D}}_{ij}{\mathsf{D}}_{ij})^{1/2}$.

Then, the relative area (or volume) flux for each grain type, $\nu ={l}$ or ${s}$, is defined through the difference between its partial velocity and the mixture velocity as $w_i^{\nu } = c^{\nu } ( v_i^{\nu }-v_i )$, so that $w_i^{l} + w_i^{s}=0$. Conservation of mass for each species requires that ${\rm D}{c}^{\nu }/{\rm D}t + \partial w_i^{\nu }/\partial x_i = 0$, where ${\rm D}(\bullet )/{\rm D}t$ is the material time derivative. Due to the fact that $c^{l} + c^{s} = 1$, only one of $c^{l}$ and $c^{s}$ is independent. Therefore, we utilize $c^{l}$ as the field variable that describes the dynamics of size segregation in the following discussion, and the evolution of $c^{l}$ is governed by its conservation of mass equation

(2.1)\begin{equation} \dfrac{{\rm D}{c}^{l}}{{\rm D}t} + \dfrac{\partial w_i^{l}}{\partial x_i} = 0. \end{equation}

2.2. Stress and the equations of motion

We recognize that the symmetric Cauchy stress tensor $\sigma _{ij} = \sigma _{ji}$ represents the Cauchy stress of the mixture, rather than the partial stress of either species. Regarding stress-related quantities for the granular mixture, we define the pressure $P=-(1/3)\sigma _{kk}$, the stress deviator $\sigma _{ij}' = \sigma _{ij} + P\delta _{ij}$, the equivalent shear stress $\tau = (\sigma _{ij}'\sigma _{ij}'/2)^{1/2}$ and the stress ratio $\mu =\tau /P$. The Cauchy stress is then governed by the standard equations of motion

(2.2)\begin{equation} \phi\rho_{s}\frac{{\rm D}v_i}{{\rm D}t}= \frac{\partial \sigma_{ij}}{\partial x_j} + b_i, \end{equation}

where $\phi$ is the constant total solid fraction, and $b_i$ is the non-inertial body force per unit volume (typically gravitational). In order to close the system of equations, we require (i) rheological constitutive equations for the Cauchy stress $\sigma _{ij}$ and (ii) a constitutive equation for the flux $w_i^{l}$, each of which is discussed in the following subsections.

2.3. Rheological constitutive equations for bidisperse mixtures

In this section, we discuss the rheology of dense, bidisperse granular mixtures. Our strategy for formulating rheological constitutive equations for bidisperse mixtures is to relate mixture-related quantities, such as the Cauchy stress $\sigma _{ij}$ and the strain-rate tensor ${\mathsf{D}}_{ij}$, instead of specifying constitutive equations for species-specific partial stresses and then combining them to obtain the mixture stress.

The starting point of this discussion is the local inertial, or $\mu (I)$, rheology (MiDi Reference MiDi2004; da Cruz et al. Reference da Cruz, Emam, Prochnow, Roux and Chevoir2005; Jop et al. Reference Jop, Forterre and Pouliquen2005), which follows from dimensional arguments. For a dense, monodisperse system of dry, stiff grains with mean grain diameter $d$ subjected to homogeneous shearing, the local inertial rheology asserts that the stress ratio $\mu$ is given through the equivalent shear strain rate $\dot {\gamma }$ and the pressure $P$ through the dimensionless relationship $\mu =\mu _{loc}(I)$, where $I=\dot {\gamma }\sqrt {d^2\rho _{s}/P}$ is the inertial number, representing the ratio of the microscopic time scale associated with particle motion $\sqrt {d^2\rho _{s}/P}$ to the macroscopic time scale of applied deformation $1/\dot {\gamma }$. As shown by Rognon et al. (Reference Rognon, Roux, Naaïm and Chevoir2007) and Tripathi & Khakhar (Reference Tripathi and Khakhar2011), the inertial rheology function $\mu _{loc}(I)$ may be straightforwardly generalized from monodisperse to bidisperse systems by defining the inertial number for a bidisperse system as $I = \dot {\gamma }\sqrt {\bar {d}^2\rho _{s}/P}$, where the average mixture grain size for a bidisperse system $\bar {d}$ has been used in place of $d$ for a monodisperse system. Then, the same local rheology function $\mu _{loc}(I)$ utilized for the monodisperse system may be used for bidisperse systems without any changes to the parameters appearing in the fitting function. This approach neglects potential effects of new dimensionless quantities that arise in a bidisperse granular system, such as the grain-size ratio $d^{l}/d^{s}$, but has been shown to capture DEM data well (Rognon et al. Reference Rognon, Roux, Naaïm and Chevoir2007; Tripathi & Khakhar Reference Tripathi and Khakhar2011).

To demonstrate this point, consider DEM simulations of homogeneous, simple shearing of a dense, bidisperse system of disks, illustrated in figure 2(a) for the case of $d^{l}/d^{s}=1.5$ and a system-wide large-grain concentration of $c^{l}=0.5$. Details of the simulated granular systems, including grain interaction properties, for both two-dimensional disks and three-dimensional spheres are given in Appendix A.1. The large particles are dark grey, and the small particles are light grey. With the system-wide mean grain size denoted by $\bar {d}_0 = c^{l}d^{l} + (1-c^{l})d^{s}$, the rectangular domain has a length of $L=60\bar {d}_0$ in the $x$-direction and a height of $H=60\bar {d}_0$ in the $z$-direction, which is filled with $\sim$5000 flowing grains. Shearing is driven through the relative motion of two parallel, rough walls, which each consist of a thin layer of touching glued grains, denoted as black in figure 2(a). Walls consisting of glued grains are utilized to minimize slip between the walls and the adjacent granular medium and mitigate any potential boundary effects associated with slip. The bottom wall is fixed, and the top wall moves with a velocity $v_{w}$ along the $x$-direction. Following previous works in the literature (da Cruz et al. Reference da Cruz, Emam, Prochnow, Roux and Chevoir2005; Koval et al. Reference Koval, Roux, Corfdir and Chevoir2009; Kamrin & Koval Reference Kamrin and Koval2012; Zhang & Kamrin Reference Zhang and Kamrin2017; Liu & Henann Reference Liu and Henann2018; Kim & Kamrin Reference Kim and Kamrin2020), the $z$-position of the top wall is not fixed but continuously adjusted using a feedback scheme so that the normal stress applied by the top wall is maintained at a target value of $\sigma _{zz}(z=0)=-P_{w}$. Periodic boundary conditions are applied along the $x$-direction. For homogeneous simple shearing, no segregation will occur since the flow is homogeneous and no pressure or strain-rate gradients are present. We utilize the DEM procedures described in detail in Liu & Henann (Reference Liu and Henann2018) in order to extract the relationship between $\mu$ and $I$ for bidisperse mixtures with grain-size ratios of $d^{l}/d^{s} = 1.5, 2.0, 2.5$ and 3.0 and $c^{l}=0.5$. The simulated relationships are plotted in figure 2(b) using triangular symbols of different colours, along with the monodisperse data from Liu & Henann (Reference Liu and Henann2018) plotted as grey circles. The relationship between $\mu$ and $I$ for dense systems of disks is observed to be approximately independent of $d^{l}/d^{s}$. As for the monodisperse case, the DEM data for bidisperse mixtures of disks may be fitted by a linear, Bingham-like functional form

(2.3)\begin{equation} \mu_{loc}(I) = \mu_{s} + bI, \end{equation}

as shown by the solid line in figure 2(b), where $\mu _{s}=0.272$ and $b=1.168$ are the dimensionless material parameters for monodisperse disks (Liu & Henann Reference Liu and Henann2018).

Figure 2. (a) Configuration for two-dimensional DEM simulations of bidisperse simple shear flow. Upper and lower layers of black grains denote rough walls. Dark grey grains indicate large flowing grains, and light grey grains indicate small flowing grains. A $10\,\%$ polydispersity is utilized for each species to prevent crystallization. (b) The local inertial rheology ($\mu$ vs $I = \dot {\gamma }\sqrt {\bar {d}^2\rho _{s}/P}$) for monodisperse as well as bidisperse mixtures of disks for grain-size ratios of $d^{l}/d^{s} = 1.5$, 2.0, 2.5 and 3.0 and $c^{l}=0.5$. The solid black line represents the best fit to the monodisperse DEM data using (2.3) with $\mu _{s}=0.272$ and $b=1.168$. (c) The local inertial rheology for monodisperse and bidisperse mixtures of spheres for grain-size ratios of $d^{l}/d^{s} = 1.5$ and 2.0 and $c^{l} = 0.5$ along with the DEM data of Tripathi & Khakhar (Reference Tripathi and Khakhar2011). The solid black curve represents the best fit to the monodisperse DEM data using (2.4) with $\mu _{s}=0.37$, $\mu _2 = 0.95$, and $I_0=0.58$.

Similarly, we consider DEM simulations of homogeneous, simple shearing of dense, bidisperse systems of spheres. The simulation domain consists of a rectangular box of length $L=20\bar {d}_0$ in the $x$-direction (i.e. the shearing direction), width $W=10\bar {d}_0$ in the $y$-direction (i.e. the direction perpendicular to the plane of shearing) and height $H=40\bar {d}_0$ in the $z$-direction. The domain is filled with $\sim$10 000 flowing grains, and periodic boundary conditions are applied along both the $x$- and $y$-directions. The simulation domain is bounded along the $z$-direction by two parallel, rough walls, consisting of touching glued grains, and as for the case of disks, shearing along the $x$-direction and normal stress along the $z$-direction are applied by the walls. We perform DEM simulations of steady simple shearing for size ratios of $d^{l}/d^{s} = 1.5$ and $2.0$ for a system-wide large-grain concentration of $c^{l}=0.5$ as well as for the monodisperse case over a range of top wall velocities. The $\mu$ vs $I$ relationship extracted from DEM simulations for these cases along with data from the prior DEM study of Tripathi & Khakhar (Reference Tripathi and Khakhar2011) collapse quite well, as shown in figure 2(c), showing minimal dependence on $d^{l}/d^{s}$. This relationship for dense systems of spheres may be fitted using the nonlinear functional form of Jop et al. (Reference Jop, Forterre and Pouliquen2005) for $\mu _{loc}(I)$

(2.4)\begin{equation} \mu_{loc}(I) = \mu_{s} + \frac{\mu_2 - \mu_{s}}{I_0/I + 1}, \end{equation}

as shown by the solid curve in figure 2(c), where $\{\mu _{s} = 0.37, \mu _2 = 0.95, I_0 = 0.58\}$ are the dimensionless parameters for frictional spheres. (We note that these parameters are nearly the same as those determined by Zhang & Kamrin (Reference Zhang and Kamrin2017) for monodisperse frictional spheres.) In this way, one may capture the rheology of bidisperse mixtures of both disks and spheres in homogeneous simple shearing without introducing additional fitting functions or adjustable parameters beyond those used for the monodisperse case.

Despite the successes of the local inertial rheology in capturing steady, homogeneous shear flow, it has been well established in the literature that a local rheological modelling approach cannot be applied to a broad set of inhomogeneous flows that span the quasi-static and dense inertial flow regimes $(I\lesssim 10^{-1})$, such as annular shear flow (Koval et al. Reference Koval, Roux, Corfdir and Chevoir2009; Tang et al. Reference Tang, Brzinski, Shearer and Daniels2018), split-bottom flow (Fenistein & van Hecke Reference Fenistein and van Hecke2003) and gravity-driven heap flow (Komatsu et al. Reference Komatsu, Inagaki, Nakagawa and Nasuno2001). In these dense, inhomogeneous flows, significant deviation from a one-to-one constitutive relationship $\mu =\mu _{loc}(I)$ is observed (Koval et al. Reference Koval, Roux, Corfdir and Chevoir2009; Kamrin & Koval Reference Kamrin and Koval2012), stemming from the fact that local rheological modelling does not account for cooperative effects, which become dominant in the quasi-static regime. Therefore, to consider inhomogeneous flows of dense, bidisperse granular systems, it is necessary to generalize a non-local rheological modelling approach to the case of dense, bidisperse mixtures. In the present work, we focus attention on the NGF model of Kamrin & Koval (Reference Kamrin and Koval2012), which has been shown to robustly capture a variety of inhomogeneous, steady flows of monodisperse granular systems (Kamrin Reference Kamrin2019). As is standard in the NGF model, we introduce the granular fluidity $g$, which is a positive, scalar field quantity, and recognize that $g$ represents the fluidity of the mixture. (See Zhang & Kamrin (Reference Zhang and Kamrin2017) and Kim & Kamrin (Reference Kim and Kamrin2020) for further discussion of the kinematic description of the granular fluidity field for monodisperse granular systems.) Then, we utilize the steady-state form of the NGF model, which relates the stress state, the strain rate, and the granular fluidity through two constitutive equations: (i) the flow rule and (ii) the non-local rheology.

First, invoking the common idealization that the Cauchy stress deviator and the strain-rate tensor are co-directional (Rycroft, Kamrin & Bazant Reference Rycroft, Kamrin and Bazant2009), the flow rule relates the Cauchy stress tensor $\sigma _{ij}$, the strain-rate tensor ${\mathsf{D}}_{ij}$ and the granular fluidity through

(2.5)\begin{equation} \sigma_{ij} =-P\delta_{ij} + 2\dfrac{P}{g}{\mathsf{D}}_{ij}. \end{equation}

Taking the magnitude of the deviatoric part of (2.5) and rearranging leads to the following scalar form of the flow rule:

(2.6)\begin{equation} \dot{\gamma} = g\mu. \end{equation}

Second, the granular fluidity of the bidisperse mixture is governed by the following differential relation:

(2.7)\begin{equation} g = g_{loc}(\mu,P) + \xi^2(\mu) \frac{\partial^2 g}{\partial x_i \partial x_i}, \end{equation}

where $g_{loc}(\mu,P)$ is the local fluidity function and $\xi (\mu )$ is the stress-dependent cooperativity length. The grain size enters (2.7) through (i) the time scale associated with microscopic particle motion that appears in the local fluidity function $g_{loc}(\mu,P)$ and (ii) the length scale that scales the cooperativity length $\xi (\mu )$, and both of these roles must be considered when generalizing the NGF model from monodisperse to bidisperse systems.

The local fluidity function gives the granular fluidity during steady, homogeneous shear flow at a given state of stress and is related to the local inertial rheology function $\mu _{loc}(I)$. Denote the inverted form of the local inertial rheology function $\mu _{loc}(I)$ as

(2.8)\begin{equation} I_{loc}(\mu) = \left\{\begin{array}{@{}ll} \mu^{-1}_{loc}(\mu) & \text{if $\mu>\mu_{s}$,}\\ 0 & \text{if $\mu\le \mu_{s}$,} \end{array}\right. \end{equation}

which is a function of the stress ratio $\mu$. Then, following the generalization of the local inertial rheology discussed above, in which the time scale associated with microscopic particle motion is taken to be $\sqrt {\bar {d}^2\rho _{s}/P}$ for a bidisperse mixture, the local fluidity function is $g_{loc}(\mu,P) = \sqrt {P/\bar {d}^2\rho _{s}}\, I_{loc}(\mu )/\mu$. For the case of bidisperse disks, using (2.3), the local fluidity function is

(2.9)\begin{equation} g_{loc}(\mu,P) = \left\{\begin{array}{@{}ll} \sqrt{\dfrac{P}{\bar{d}^2\rho_{s}}} \, \dfrac{(\mu-\mu_{s})}{b\mu} & \text{if $\mu > \mu_{s}$,}\\ 0 & \text{if $\mu\le\mu_{s}$,}\end{array}\right. \end{equation}

with $\{\mu _{s}=0.272, b=1.168\}$, and for the case of bidisperse spheres, using (2.4), the local fluidity function is

(2.10)\begin{equation} g_{loc}(\mu,P) = \left\{\begin{array}{@{}ll} I_0\sqrt{\dfrac{P}{\bar{d}^2\rho_{s}}}\,\dfrac{(\mu-\mu_{s})}{\mu(\mu_2-\mu)} & \text{if $\mu>\mu_{s}$,}\\[10pt] 0 & \text{if $\mu\le\mu_{s}$,} \end{array} \right. \end{equation}

with $\{\mu _{s} = 0.37, \mu _2 = 0.95, I_0 = 0.58\}$. No additional adjustable parameters beyond those used to describe the local inertial rheology for the monodisperse case are introduced in the local fluidity function.

As discussed in several of our previous works (Henann & Kamrin Reference Henann and Kamrin2014; Kamrin & Henann Reference Kamrin and Henann2015; Liu & Henann Reference Liu and Henann2017), the manner in which the cooperativity length $\xi (\mu )$ depends on the stress ratio $\mu$ is also connected to the choice of the $\mu _{loc}(I)$ function. Without going into details here, the functional forms for the cooperativity length corresponding to (2.3) and (2.4) are

(2.11a,b)\begin{equation} \xi(\mu) = \frac{A\bar{d}}{\sqrt{|\mu - \mu_{s}|}}\quad\text{and}\quad \xi(\mu) = A\bar{d}\sqrt{\dfrac{(\mu_2-\mu)}{(\mu_2-\mu_{s})|\mu-\mu_{s}|}}, \end{equation}

respectively. The parameter $A$ is a dimensionless material constant, referred to as the non-local amplitude, which quantifies the spatial extent of cooperative effects. In the monodisperse case, the cooperativity length is directly proportional to the grain size $d$, and motivated by the success in generalizing the local inertial rheology, we follow an analogous approach to generalize the cooperativity length to the bidisperse case. We replace $d$ for monodisperse grains with $\bar {d}$ for bidisperse grains, resulting in the expressions for the cooperativity length (2.11a,b), in which $\xi (\mu )$ is proportional to $\bar {d}$. Regarding the non-local amplitude $A$, we also follow an approach that is analogous to the generalization of the local inertial rheology and utilize values of $A$ previously determined for monodisperse frictional disks and spheres – namely, $A=0.90$ as determined by Liu & Henann (Reference Liu and Henann2018) for monodisperse disks and $A=0.43$ as determined by Zhang & Kamrin (Reference Zhang and Kamrin2017) for monodisperse spheres. The generalization of the cooperativity length and the choices of $A$ for bisdisperse mixtures will be tested in later sections by comparing flow fields predicted by the NGF model with measured flow fields in DEM simulations of bidisperse, inhomogeneous flows.

2.4. Segregation model

The segregation model consists of the constitutive equation for the large-grain flux $w_i^{l}$. In the present work, we focus on dense flows in the absence of pressure gradients, and we take the large-grain flux $w_i^{l}$ to comprise two contributions: (i) a diffusion flux $w_i^{diff}$ and (ii) a shear-strain-rate-gradient-driven segregation flux $w_i^{seg}$, so that

(2.12)\begin{equation} w_i^{l} = w_i^{diff} + w_i^{seg}. \end{equation}

First, the diffusion flux acts counter to segregation to mix the species and is taken to be given in the standard form, in which the diffusion flux is driven by concentration gradients: $w_i^{diff} = -D ( \partial c^{l}/\partial x_i )$, where $D$ is the binary diffusion coefficient (Utter & Behringer Reference Utter and Behringer2004; Artoni et al. Reference Artoni, Larcher, Jenkins and Richard2021; Bancroft & Johnson Reference Bancroft and Johnson2021). Following prior works (e.g. Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021; Duan et al. Reference Duan, Umbanhowar, Ottino and Lueptow2021; Trewhela, Ancey & Gray Reference Trewhela, Ancey and Gray2021), based on dimensional arguments, we take the diffusion coefficient for a bidisperse system to be

(2.13)\begin{equation} D = C_{diff} \bar{d}^2 \dot{\gamma}, \end{equation}

where $C_{diff}$ is a dimensionless material parameter which remains to be calibrated. Therefore, we have that the diffusion flux is

(2.14)\begin{equation} w_i^{diff} =-C_{diff} \bar{d}^2 \dot{\gamma} \dfrac{\partial c^{l}}{\partial x_i}. \end{equation}

Second, regarding segregation, a major question is what field quantity drives the segregation flux in the absence of pressure gradients. Gradients of a number of kinematic quantities are possible – e.g. strain rate, velocity fluctuations or fluidity. For perspective, we note that recent works (Fan & Hill Reference Fan and Hill2011b; Hill & Tan Reference Hill and Tan2014; Tunuguntla, Thornton & Weinhart Reference Tunuguntla, Thornton and Weinhart2016; Tunuguntla et al. Reference Tunuguntla, Weinhart and Thornton2017) have shown that gradients in kinetic stress, which are defined through the velocity fluctuations, correlate well with segregation flux. In the present work, we adopt the simplest approach and hypothesize that the segregation flux is driven by gradients in the shear strain rate $\dot {\gamma }$ and take the segregation flux to be given in the following phenomenological form:

(2.15)\begin{equation} w_i^{seg} = C_{seg}^{S} \bar{d}^2 c^{l}(1-c^{l}) \dfrac{\partial \dot{\gamma}}{\partial x_i}. \end{equation}

The factor $c^{l}(1-c^{l})$ ensures that segregation ceases when the bidisperse mixture becomes either all large $(c^{l} = 1)$ or all small $(c^{l} = 0)$ grains, and the factor $\bar {d}^2$ is present for dimensional consistency. The quantity $C^{S}_{seg}$ is a dimensionless material property. While it is possible for $C^{S}_{seg}$ to depend on the size ratio $d^{l}/d^{s}$, we will demonstrate that this effect is negligible over the range of size ratios considered in the DEM simulations of § 4 and therefore treat $C^{S}_{seg}$ as a constant, dimensionless material parameter, which will be determined by fitting to DEM simulation results for disks and spheres, respectively.

Combining (2.14), (2.15) and (2.12) with conservation of mass (2.1), we obtain the following differential relation governing the dynamics of $c^{l}$:

(2.16)\begin{equation} \dfrac{{\rm D}{c}^{l}}{{\rm D}t} + \dfrac{\partial}{\partial x_i} \left(-C_{diff} \bar{d}^2 \dot{\gamma} \dfrac{\partial c^{l}}{\partial x_i} + C^{S}_{seg} \bar{d}^2 c^{l}(1-c^{l}) \dfrac{\partial \dot{\gamma}}{\partial x_i} \right) = 0 , \end{equation}

where $\{ C_{diff}, C^{S}_{seg} \}$ represent two constant dimensionless material parameters that remain to be determined.

We close this section by noting that the incompressibility constraint, the equations of motion (2.2), the non-local rheology (2.7) and the segregation dynamics equation (2.16) represent a closed system of equations for the velocity field $v_i$, the pressure field $P$, the fluidity field $g$ and the large-grain concentration field $c^{l}$, which may be used to simultaneously predict flow fields and the segregation dynamics in the absence of pressure gradients.

4.1. Vertical chute flow

Consider a dense, bidisperse granular mixture flowing down a long vertical chute with parallel, rough walls separated by a distance $W$ under the action of gravity $G$. This flow geometry has been utilized extensively in the literature to study dense flows of monodisperse, frictional disks (Kamrin & Koval Reference Kamrin and Koval2012; Liu & Henann Reference Liu and Henann2018) and spheres (Zhang & Kamrin Reference Zhang and Kamrin2017; Kim & Kamrin Reference Kim and Kamrin2020) as well as flows of bidisperse, frictional spheres (Fan & Hill Reference Fan and Hill2011a,Reference Fan and Hillb). Beginning with the case of bidisperse disks, the DEM set-up is shown in figure 4(a) for $W=60\bar {d}_0$, where $\bar {d}_0$ is the system-wide average grain size. In all cases for disks, we take the chute length to be $L=60\bar {d}_0$ and apply periodic boundary conditions along the $z$-direction. The parallel, rough walls consist of touching glued large grains, denoted as black in figure 4(a). The left vertical wall is fixed, and the right wall is fixed in the $z$-direction but can move slightly in the $x$-direction so as to maintain a constant compressive normal stress $P_{w}$ on the granular material, utilizing the same wall-position control method described in Liu & Henann (Reference Liu and Henann2018). We have verified that the chute length $L$ is sufficiently large, so that it does not affect the resulting flow and segregation fields and all fields are invariant along the $z$-direction. In the resulting flow fields, the only non-zero component of the velocity is $v_z$, which only depends on the cross-channel coordinate $x$. A typical steady velocity field is qualitatively sketched in figure 4(a), illustrating that the shear strain rate is greatest at the walls $(x = \pm W/2)$.

Figure 4. (a) Initial well-mixed configuration for two-dimensional DEM simulation of bidisperse vertical chute flow with $4327$ flowing grains. The chute width is $W=60\bar {d}_0$. As in figure 2, black grains on both sides represent rough walls. (Only large particles are used as wall grains here.) (b) Segregated configuration after flowing for a total simulation time of $\tilde {t} = t/( d^{s}\sqrt {\rho _{s}/P_{w}} )=4.3 \times 10^5$. (c) Spatio-temporal evolution of the large-grain concentration field. Spatial profiles of (d) the concentration field $c^{l}$ and (e) the normalized velocity field $(v_{cen}-v_z)\sqrt {\rho _{s}/P_{w}}$ at three times ($\tilde {t} =4 \times 10^3$, $4 \times 10^4$ and $4 \times 10^5$) as indicated by the horizontal lines in (c).

In all of our DEM simulations of vertical chute flow of bidisperse disks, we observe that the stress field quickly becomes independent of time, so that macroscopic inertia (i.e. the left-hand side of (2.2)) may be neglected. Moreover, we observe that the normal stresses are approximately equal, i.e. $\sigma _{zz}\approx \sigma _{xx}$. Therefore, due to the force balance along the $z$-direction, the equivalent shear stress field is $\tau (x) = |\sigma _{xz}(x)| = |\sigma _{zx}(x)| = \phi \rho _{s}G|x|$, where $x$ is measured from the centreline of the chute, and due to the force balance along the $x$-direction, the pressure field is $P(x) = -\sigma _{xx}(x) = P_{w}$. The stress ratio field then follows as

(4.1)\begin{equation} \mu(x) = \mu_{w}\left(\dfrac{|x|}{W/2}\right), \end{equation}

where $\mu _{w}=\phi \rho _{s}GW/2P_{w}$ is the maximum value of $\mu$, occurring at the walls $(x = \pm W/2)$. We note that while flow is driven by gravity, the pressure field is constant throughout the chute and no pressure gradients are present. Therefore, segregation occurs only due to shear-strain-rate gradients, enabling us to consider this effect in isolation.

Apart from the grain interaction properties that are held constant throughout this work (Appendix A.1), there are four important dimensionless parameters that fully describe each case of vertical chute flow of dense, bidisperse granular mixtures: (i) $W/\bar {d}_0$, the dimensionless chute width; (ii) $\mu _{w}$, the maximum stress ratio, which occurs at the walls and controls the total flow rate; (iii) $c^{l}_0(x)$, the initial large-grain concentration, which is not necessarily constant but can be a spatially varying field; and (iv) $d^{l}/d^{s}$, the bidisperse grain-size ratio. This list of system parameters $\{W/\bar {d}_0,$ $\mu _{w},$ $c^{l}_0,$ $d^{l}/d^{s}\}$ specifies the geometry, loads and initial conditions of a given case of vertical chute flow. As a representative base case, we consider the parameter group $\{W/\bar {d}_0 = 60,\mu _{w}=0.45,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$. We then run the corresponding DEM simulation starting from the well-mixed initial configuration shown in figure 4(a) and observe that, after a simulation time of $\tilde {t} = t/( d^{s}\sqrt {\rho _{s}/P_{w}} )= 4.3 \times 10^5$, the large, dark-grey grains segregate towards the regions near the walls where the shear strain rate is greatest, while the small, light-grey grains gather in bands just inside these regions, as shown in figure 4(b). A well-mixed core persists along the centre of the vertical chute where the shear strain rate is nearly zero. To obtain a more quantitative picture of the segregation process, we coarse grain the concentration field $c^{l}$ in both space and time and plot contours of the spatio-temporal evolution of $c^{l}$ in figure 4(c). The large-grain concentration field evolves rapidly in time during the initial stages of the segregation process. Then, over longer times, the evolution becomes slower. Spatial profiles of the concentration and velocity fields at three snapshots in time – specifically, $\tilde {t} = t/( d^{s}\sqrt {\rho _{s}/P_{w}} ) =4 \times 10^3$, $4 \times 10^4$ and $4 \times 10^5$ as indicated by the horizontal lines in figure 4(c) – are plotted in figures 4(d) and 4(e). These three snapshots correspond to early, medium, and late times with respect to the segregation process. The spatial $c^{l}$ profiles shown in figure 4(d) demonstrate the transition from a well-mixed state to a segregated state with large-grain-rich and small-grain-rich regions. In figure 4(e), the normalized velocity fields $(v_{cen} - v_z)\sqrt {\rho _{s}/P_{w}}$, relative to the velocity at the centre of the chute, $v_{cen} = v_z(x=0)$, show that the velocity field rapidly develops into a steady flow field, even while the segregation process is still ongoing, and the $c^{l}$ field continues to evolve.

At long times, near the end of the simulated time window ($\tilde {t} = t/( d^{s}\sqrt {\rho _{ s}/P_{w}} )\gtrsim 3 \times 10^5$), the concentration field evolves very slowly, so that ${\rm D}{c}^{l}/{\rm D}t \approx 0$, and the state of segregation may be regarded as quasi-steady. Therefore, according to (2.1) and (2.12) and the no-flux boundary condition at the walls, the total flux is approximately zero ($w_i^{l} = w_i^{diff} + w_i^{seg} \approx 0_i$) at each $x$-position, meaning that the segregation flux is approximately balanced by the diffusion flux in this quasi-steady flow regime. Then, using the expressions for the two fluxes, (2.14) and (2.15), this observation implies that

(4.2)\begin{equation} C_{diff}\bar{d}^2\dot{\gamma}\dfrac{\partial c^{l}}{\partial x} \approx C^{S}_{seg} \bar{d}^2c^{l}(1-c^{l})\dfrac{\partial \dot{\gamma}}{\partial x}. \end{equation}

The field quantities appearing in this expression may be obtained by coarse graining the DEM data in the quasi-steady flow regime. Therefore, since $C_{diff}$ has been previously determined, (4.2) may be used to determine the parameter $C^{S}_{seg}$ as follows. First, we acquire the field quantities $c^{l}$ (and hence $\bar {d}$) and $v_z$ by spatially coarse graining $152$ evenly distributed snapshots in time in the quasi-steady regime ($\tilde {t} > 3 \times 10^5$). (Results are not particularly sensitive to the quasi-steady regime criterion, and this value is only provided as a guideline.) Then, we arithmetically average these fields in time, yielding fields that only depend on the spatial coordinate $x$, and take spatial gradients, as described in Appendix A.2, to obtain ${\partial c^{l}}/{\partial x}$, $\dot {\gamma } = {\partial v_z}/{\partial x}$ and ${\partial \dot {\gamma }}/{\partial x} = {\partial ^2 v_z}/{\partial x^2}$. Next, as suggested by (4.2), we plot $C_{diff}\bar {d}^2\dot {\gamma }({\partial c^{l}}/{\partial x})$ vs $\bar {d}^2c^{l}(1-c^{l})({\partial \dot {\gamma }}/{\partial x})$ in figure 5(a), in which each symbol represents a unique $x$-position. A linear relation is observed, supporting our choice for the form of the constitutive equation for the segregation flux (2.15). In order to obtain further evidence for this choice, we consider four additional cases: (i) a lower flow rate case $\{W/\bar {d}_0 = 60,\mu _{w}=0.375,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$; (ii) a narrower channel case $\{W/\bar {d}_0 = 40,\mu _{w}=0.45,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$; (iii) a more large grains case $\{W/\bar {d}_0 = 60,\mu _{w}=0.45,c^{l}_0=0.75,d^{l}/d^{s}=1.5\}$; and (iv) a larger size ratio case $\{W/\bar {d}_0 = 60,\mu _{w}=0.45,c^{l}_0=0.5,d^{l}/d^{s}=3.0\}$. Coarse graining the quasi-steady fields for each case and including the field data in figure 5(a), we observe a strong linear collapse. Finally, the dimensionless material parameter $C_{seg}^{S}$ may be obtained from the slope of the linear relation in figure 5(a) (indicated by the solid line). We determine the numerical value for disks to be $C^{S}_{seg}=0.23$ and note that this value indeed appears to be independent of the grain-size ratio $d^{l}/d^{s}$ over the range of 1.5 to 3.0 for disks.

Figure 5. Collapse of $C_{diff}\bar {d}^2\dot {\gamma }({\partial c^{l}}/{\partial x})$ vs $\bar {d}^2c^{l}(1-c^{l})({\partial \dot {\gamma }}/{\partial x})$ for several cases of vertical chute flow of (a) bidisperse disks and (b) bidisperse spheres. Symbols represent coarse-grained, quasi-steady DEM field data, and the solid lines are the best linear fits using (a) $C^{S}_{seg}=0.23$ for disks and (b) $C^{S}_{seg}=0.08$ for spheres.

We carry out an analogous set of DEM simulations for dense, bidisperse mixtures of spheres to determine the value of $C^{S}_{seg}$ for spheres. The DEM set-up for spheres is similar to that shown in figure 4(a) for disks with a domain size of length $L=20\bar {d}_0$ in the $z$-direction, width $W$ in the $x$-direction, which is varied in our DEM simulations, and out-of-plane thickness $H=10\bar {d}_0$ in the $y$-direction. Periodic boundary conditions are applied along both the $y$- and $z$-directions for spheres, and a constant compressive normal stress $\sigma _{xx} = -P_w$ is applied using the same feedback scheme as utilized for disks. In DEM simulations of dense flows of spheres, we observe normal stress differences, in which the normal stresses $\sigma _{yy}$ and $\sigma _{zz}$ are slightly different from the prescribed value of $\sigma _{xx}= -P_w$, which is a widely reported feature of dense flows of spheres in the literature (e.g. Srivastava et al. Reference Srivastava, Silbert, Grest and Lechman2021). However, all normal stresses, and hence the pressure field, are spatially uniform, and segregation occurs only due to shear-strain-rate gradients. As for disks, the set of system parameters $\{W/\bar {d}_0,$ $\mu _{w},$ $c^{l}_0,$ $d^{l}/d^{s}\}$ specifies the geometry, loads and initial conditions for a given case of vertical chute flow. Here, the dimensionless parameter $\mu _{w} = \phi \rho _{s}GW/2P_{w}$ continues to control the total flow rate down the chute. We consider five different cases: (i) a base case $\{W/\bar {d}_0 = 60,\mu _{w}=0.51,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$; (ii) a lower flow rate case $\{W/\bar {d}_0 = 60,\mu _{w}=0.46,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$; (iii) a narrower channel and higher flow rate case $\{W/\bar {d}_0 = 50,\mu _{w}=0.59,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$; (iv) a more large grains and higher flow rate case $\{W/\bar {d}_0 = 60,\mu _{w}=0.58,c^{l}_0=0.75,d^{l}/d^{s}=1.5\}$; and (v) a larger size ratio and higher flow rate case $\{W/\bar {d}_0 = 60,\mu _{w}=0.58,c^{l}_0=0.5,d^{l}/d^{s}=2.0\}$. Each case involves $\sim$20 000 flowing grains. After a sufficiently long simulation time, a quasi-steady state is attained in each case, implying the flux balance (4.2). The quasi-steady field quantities appearing in (4.2) are extracted for 1000 snapshots in time using the coarse-graining techniques described in Appendix A.2 and then arithmetically averaged in time. The calculated quasi-steady diffusion flux $C_{diff}\bar {d}^2\dot {\gamma }({\partial c^{l}}/{\partial x})$ at discrete $x$-positions for each of the five cases is plotted vs the calculated quantity $\bar {d}^2c^{l}(1-c^{l})({\partial \dot {\gamma }}/{\partial x})$ in figure 5(b), and we observe a linear relation. Therefore, the form of the constitutive equation for the segregation flux (4.2) is also applicable to dense, bidisperse systems of spheres. The numerical value of the dimensionless material parameter $C^{S}_{seg}$ for spheres is determined from the slope of the linear relation to be $C^{S}_{seg} = 0.08$, which appears to be independent of the grain-size ratio $d^{l}/d^{s}$ over the range of 1.5 to 2.0 for spheres.

4.2. Annular shear flow

The constitutive equation for the segregation flux (2.15) and the fitted values of the material parameters should be general across different flow geometries. To test this for the case of disks, we apply the same process described in the preceding section for vertical chute flow to a different flow geometry – annular shear flow. In this flow geometry, flow is driven through motion of the boundary rather than by gravity, but as in vertical chute flow, the pressure field is spatially uniform, which eliminates hydrostatic pressure gradients so that only shear-strain-rate-gradient-driven size segregation occurs.

Our DEM simulations of annular shear flow of a dense, bidisperse granular mixture of disks follow the procedures utilized in prior works in the literature for monodisperse, frictional disks (Koval et al. Reference Koval, Roux, Corfdir and Chevoir2009; Kamrin & Koval Reference Kamrin and Koval2012, Reference Kamrin and Koval2014; Liu & Henann Reference Liu and Henann2018). Consider a dense, bidisperse granular mixture in a two-dimensional annular shear cell with rough circular walls of inner radius $R$ and outer radius $R_{o}$, as shown in figure 6(a) for the case of $R=60\bar {d}_0$. The inner and outer walls consist of rings of glued large grains, denoted as black in figure 6(a). The circumferential velocity of the inner wall is prescribed to be $v_{w}$, and its radial position is fixed. The outer wall does not rotate, and its radius $R_{o}$ fluctuates slightly so as to maintain a constant imposed compressive normal stress $P_{w}$, utilizing the wall-position control method used throughout this work (Koval et al. Reference Koval, Roux, Corfdir and Chevoir2009). As in Liu & Henann (Reference Liu and Henann2018), we simulate the full shear cell, as shown in figure 6(a), instead of applying periodic boundary conditions along the circumferential direction to a slice (Koval et al. Reference Koval, Roux, Corfdir and Chevoir2009; Kamrin & Koval Reference Kamrin and Koval2012, Reference Kamrin and Koval2014). In this flow geometry, all fields are axisymmetric (i.e. invariant along the $\theta$-direction), and the only non-zero component of the velocity field is the circumferential component $v_\theta$. Moreover, flow tends to localize near the inner wall with $v_\theta$ rapidly decaying with radial position, as qualitatively illustrated by the steady velocity field sketched in figure 6(a). We choose the outer radius $R_{o}$ to be sufficiently large so that it does not affect the resulting flow and segregation fields, and based on our experience, we take $R_{o}=2R$. Therefore, the role of the outer wall is simply to apply a far-field pressure, and otherwise, it does not affect the flow and segregation fields.

Figure 6. (a) Initial well-mixed configuration for two-dimensional DEM simulation of bidisperse annular shear flow with 40 108 flowing grains. The inner-wall radius is $R=60\bar {d}_0$, and the outer-wall radius is $R_{o}=2R$. The inner and outer walls consist of rings of glued large grains, denoted as black. (b) Segregated configuration after flowing for a total simulation time of $\tilde {t}=t/( R/v_{w} )=584$. (c) Spatio-temporal evolution of the large-grain concentration field. Spatial profiles of (d) the concentration field $c^{l}$ and (e) the normalized circumferential velocity field ${v}_\theta /v_{w}$ at three times ($\tilde {t}=5$, $50$ and $500$) as indicated by the horizontal lines in (c).

Regarding the stress field, in all of our DEM simulations of annular shear flow of bidisperse disks, we observe that the normal stresses are approximately equal, i.e. $\sigma _{rr}\approx \sigma _{\theta \theta }$, and spatially uniform, so that the force balance along the $r$-direction gives that the pressure field is $P(r) = -\sigma _{rr}(r) = P_{w}$. Since the pressure field is spatially uniform, all segregation in annular shear flow is due to shear-strain-rate gradients. The moment balance gives that the equivalent shear stress field is $\tau (r) = |\sigma _{r\theta }(r)| = |\sigma _{\theta r}(r)| = \tau _{w}(R/r)^2$, where $\tau _{w}$ is the inner-wall shear stress. It is important to note that $\tau _{w}$ is not directly prescribed in our DEM simulations. Instead, the inner-wall velocity $v_{w}$ is prescribed, and $\tau _{w}$ arises as a result. The stress ratio field is then

(4.3)\begin{equation} \mu(r) = \mu_{w}(R/r)^2, \end{equation}

where $\mu _{w}=\tau _{w}/P_{w}$ is the maximum value of $\mu$, occurring at the inner wall $(r = R)$.

As for vertical chute flow, there are four important dimensionless parameters that specify the geometry, loads and initial conditions for a given case of annular shear flow of dense, bidisperse granular mixtures: (i) $R/\bar {d}_0$, the dimensionless inner-wall radius; (ii) $\tilde {v}_{w} = (v_{w}/R)\sqrt {{\rm \pi} \bar {d}_0^2\rho _{s}/(4P_{w})}$, the dimensionless inner-wall velocity, which determines $\tau _{w}$ and hence $\mu _{w}$; (iii) $c^{l}_0(r)$, the initial large-grain concentration field; and (iv) $d^{l}/d^{s}$, the bidisperse grain-size ratio. We choose a representative base case of annular shear flow identified by the parameter set $\{R/\bar {d}_0 = 60,\tilde {v}_{w}=0.01,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$. The well-mixed initial configuration for the base-case DEM simulation is shown in figure 6(a), and the segregated configuration after driving flow for a total simulation time of $\tilde {t}=t/( R/v_{w} )=584$ is shown in figure 6(b). The large, dark-grey grains segregate into a ring near the inner wall, while the small, light-grey grains form a band just outside this region. Outside of these bands, where the shear strain rate is very small, the large and small grains remain well mixed. Contours of the spatio-temporal evolution of the coarse-grained concentration field $c^{l}$ are plotted in figure 6(c), illustrating the time evolution of the $c^{l}$ field. Spatial profiles of the concentration and velocity fields at three selected snapshots during the segregation dynamics ($\tilde {t}=t/( R/v_{w} )=5, 50$ and 500 as indicated by the dashed lines in figure 6c) are shown in figures 6(d) and 6(e). The radial $c^{l}$ profiles in figure 6(d) demonstrate the formation of large-grain-rich and small-grain-rich regions with a persistent well-mixed far field. The normalized velocity fields in figure 6(e) demonstrate that the flowing zone is localized near the inner wall with slow creeping flow observed far from the wall. As in the case of vertical chute flow, the velocity field quickly develops into a steady flow field, while the large-grain concentration field $c^{l}$ evolves over a longer time scale.

Again, at long times, near the end of the simulated time window ($\tilde {t}=t/( R/v_{w} )\gtrsim 500$), the concentration field evolves very slowly, which we identify as the quasi-steady regime. (We note that this is not a true steady state and that over even longer times, the concentration field will continue to evolve very slowly to more and more segregated states.) In this regime, the segregation and diffusion fluxes approximately balance at each $r$-position, implying that

(4.4)\begin{equation} C_{diff}\bar{d}^2\dot{\gamma}\dfrac{\partial c^{l}}{\partial r} \approx C^{S}_{seg} \bar{d}^2c^{l}(1-c^{l})\dfrac{\partial \dot{\gamma}}{\partial r} . \end{equation}

As for vertical chute flow, we spatially coarse grain the DEM data to obtain the $c^{l}$ and $v_\theta$ fields for $144$ evenly distributed snapshots in time in the quasi-steady regime ($\tilde {t} \gtrsim 500$), which are arithmetically averaged in time to obtain the quasi-steady $c^{l}(r)$ and $v_\theta (r)$ fields and then spatially differentiated to obtain the remaining field quantities in (4.4). Next, we plot $C_{diff}\bar {d}^2\dot {\gamma }({\partial c^{l}}/{\partial r})$ vs $\bar {d}^2c^{l}(1-c^{l})({\partial \dot {\gamma }}/{\partial r})$ in figure 7 with each symbol representing a unique $r$-position. Finally, this process is repeated for four additional cases of annular shear flow: (i) a lower inner-wall velocity $\{R/\bar {d}_0 = 60,\tilde {v}_{w}=0.001,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$; (ii) a smaller inner-wall radius $\{R/\bar {d}_0 = 40,\tilde {v}_{w}=0.01,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$; (iii) more large grains $\{R/\bar {d}_0 = 60,\tilde {v}_{w}=0.01,c^{l}_0=0.75,d^{l}/d^{s}=1.5\}$; and (iv) a larger size ratio $\{R/\bar {d}_0 = 60,\tilde {v}_{w}=0.01,c^{l}_0=0.5,d^{l}/d^{s}=3.0\}$, and the coarse-grained, quasi-steady fields are included in the data plotted in figure 7. Collectively, we observe a strong collapse to a linear relation. Crucially, the slope of the linear relation in figure 7 (indicated by the solid line) gives the same value for the dimensionless material parameter $C_{seg}^{S}$ obtained by fitting to vertical chute flow data, $C_{seg}^{S}=0.23$. This observation of agreement between the best fit values of $C_{seg}^{S}$ obtained using two different flow geometries provides support for our choice of the constitutive equation for the segregation flux (2.15) and the fitted value of $C^{S}_{seg}$ for disks. Having established that the parameter $C^{S}_{seg}$ is independent of the flow geometry, driving conditions and initial conditions, we henceforth regard $C^{S}_{seg}$ as a material parameter for a given dense granular system, analogous to how rheological material parameters such as $\mu _{s}$ are regarded. Of course, the values of $C^{S}_{seg}$ determined above for disks and spheres likely depend on grain interaction properties, such as the inter-particle friction coefficient, but elucidating this dependence is beyond the scope of the present work.

Figure 7. Collapse of $C_{diff}\bar {d}^2\dot {\gamma }({\partial c^{l}}/{\partial r})$ vs $\bar {d}^2c^{l}(1-c^{l})({\partial \dot {\gamma }}/{\partial r})$ for several cases of annular shear flow of bidisperse disks. Symbols represent coarse-grained, quasi-steady DEM field data, and the solid line is the best linear fit using $C^{S}_{seg}=0.23$.

5.1. Vertical chute flow

First, we describe in detail how the continuum model is solved to obtain predictions for the transient evolution of segregation and flow fields for the case of vertical chute flow of disks. In vertical chute flow, the stress field may be straightforwardly deduced from a static force balance, giving that the pressure field is uniform, $P(x) = P_{w}$, and that the stress ratio field $\mu (x)$ is given through (4.1) and, therefore, the balance of linear momentum (2.2) is satisfied and does not further enter the solution procedure. Continuum model predictions are obtained by numerically solving the remaining governing equations using finite differences. Summarizing the coupled boundary/initial-value problem for flow and segregation in the context of vertical chute flow, the unknown fields are the velocity field $v_z(x,t)$ and the accompanying strain-rate field $\dot {\gamma }(x,t) = \partial v_z/\partial x$, the granular fluidity field $g(x,t)$ and the large-grain concentration field $c^{l}(x,t)$. The governing equations are (i) the flow rule (2.6)

(5.3)\begin{equation} \dot{\gamma} = g\mu, \end{equation}

(ii) the non-local rheology (2.7)

(5.4)\begin{equation} g = g_{loc}(\mu,P_{w}) + \xi^2(\mu)\dfrac{\partial^2 g}{\partial x^2}, \end{equation}

with $g_{loc}$ and $\xi$ given through (2.9) and (2.11a), respectively, and (iii) the segregation dynamics equation (2.16)

(5.5)\begin{equation} \dfrac{\partial c^{l}}{\partial t} + \dfrac{\partial}{\partial x}\left(-C_{diff}\bar{d}^2\dot{\gamma}\dfrac{\partial c^{l}}{\partial x} + C^{S}_{seg}\bar{d}^2c^{l}(1 - c^{l})\dfrac{\partial \dot{\gamma}}{\partial x}\right) = 0, \end{equation}

where $\bar {d} = c^{l}d^{l} + (1-c^{l})d^{s}$.

The non-local rheology (5.4) requires non-standard boundary conditions for the granular fluidity field at the walls. Walls consisting of touching glued grains are employed in the DEM simulations throughout this work, and to select a reasonable fluidity boundary condition for this type of wall, we look to DEM simulation results for the simple shear flows of figure 2, in which we observe a spatially uniform strain-rate field with minimal deviation near the walls. Therefore, we elect to use a fluidity boundary condition that would predict a uniform strain-rate field when applied to these simple shear flows. There are two options: (i) a Dirichlet boundary condition wherein the fluidity at the boundary is a function of the stress state through the local fluidity function, i.e. $g = g_{loc}(\mu _{w},P_{w})$ at $x=\pm W/2$ for vertical chute flow, or (ii) a homogeneous Neumann fluidity boundary condition, i.e. $\partial g/\partial x$ at $x=\pm W/2$ for vertical chute flow. Other choices for the fluidity boundary condition would result in a non-uniform solution for the fluidity field (and hence the strain-rate field) when applied to simple shear flow, in which deviations from a uniform strain-rate field arise in boundary layers at the walls. We have tested both of these options for the fluidity boundary conditions and found that continuum model predictions are similar for all cases of vertical chute flow and annular shear flow considered in the present work. In all continuum model predictions that follow, we impose inhomogeneous Dirichlet fluidity boundary conditions at the walls. Regarding boundary conditions for (5.5), we impose no-flux boundary conditions at the walls, i.e. $w_x^{l} = -C_{diff} \bar {d}^2 \dot {\gamma } (\partial c^{l}/\partial x) + C^{S}_{seg} \bar {d}^2 c^{l}(1 - c^{l})(\partial \dot {\gamma }/\partial x) = 0$ at $x=\pm W/2$. Due to the time derivative in (5.5), an initial condition for the concentration field $c^{l}_0(x) = c^{l}(x,t=0)$ is required. In order to account for the concentration fluctuations inherent in the initial state, we obtain the coarse-grained $c^{l}$ field from the initial DEM configuration for each case and utilize this field as the initial condition field $c_0^{l}(x)$ in each of the respective continuum simulations.

Then, for a given case identified through a set of input parameters $\{W/\bar {d}_0,$ $\mu _{w},$ $c^{l}_0(x),$ $d^{l}/d^{s}\}$, we obtain numerical predictions of the continuum model utilizing finite differences as follows. First, at a given point in time, the concentration field $c^{l}(x)$ is known, allowing the average grain-size field to be calculated through $\bar {d}(x) = c^{l}(x)d^{l} + (1-c^{l}(x))d^{s}$. Using the stress ratio field $\mu (x)$ for vertical chute flow (4.1), the local fluidity $g_{loc}(\mu,P)$ and the cooperativity length $\xi (\mu )$ ((2.9) and (2.11a)) may be calculated at each spatial grid point. Then, the non-local rheology (5.4) may be used to solve for the fluidity field $g(x)$ at the current step, using central differences in space. The strain-rate field follows using (5.3), which may be integrated to obtain the velocity field $v_z(x)$. Next, (5.5) is used to determine the concentration field at the next time step utilizing the forward Euler method and central differences in space with one modification – the spatial derivatives of $c^{l}$ appearing in the diffusion flux term in (5.5) are treated implicitly in order to improve numerical stability. This completes one time step, and this process is repeated to step forward in time and calculate the transient evolution of the concentration and flow fields, $c^{l}(x,t)$ and $v_z(x,t)$. In our finite-difference calculations, we utilize a fine spatial resolution of $\Delta x \ll \bar {d}_0$, and we have verified that the time step is sufficiently small in order to ensure stable, accurate results.

We compare predictions of the continuum model against DEM data for all five cases of vertical chute flow considered in § 4.1 in order to test the generality of the model. Figures 8 and 9 summarize the comparisons for these five cases. The first columns of figures 8 and 9 show the spatio-temporal contours of the evolution of the $c^{l}$ field measured in the DEM simulations for each case, and the second columns show comparisons of the DEM simulations (solid black lines) and the continuum model predictions (dashed grey lines) for the $c^{l}$ field at four snapshots in time ($\tilde {t} = t/( d^{s}\sqrt {\rho _{s}/P_{w}} ) =4 \times 10^3$, $2 \times 10^4$, $1 \times 10^5$ and $4 \times 10^5$), indicated by the horizontal lines in the first columns of figures 8 and 9. Based on figures 8 and 9, the coupled model generally does a good job capturing the salient features of the evolution of the $c^{l}$ field across all cases. For instance, for the narrower chute case shown in figure 8(c), the segregation process nearly completes within the simulated time window with the mixed core along the centre of the chute nearly disappearing, and the continuum model prediction captures this observation well.

Figure 8. Comparisons of continuum model predictions with corresponding DEM simulation results for the transient evolution of the segregation dynamics for three cases of vertical chute flow of disks. (a) Base case $\{W/\bar {d}_0 = 60,\mu _{w}=0.45,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$. (b) Lower flow rate case $\{W/\bar {d}_0 = 60,\mu _{w}=0.375,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$. (c) Narrower chute width case $\{W/\bar {d}_0 = 40,\mu _{w}=0.45,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$. Additional cases are shown in figure 9. For each case, the first column shows spatio-temporal contours of the evolution of $c^{l}$ measured in the DEM simulations. The second column shows comparisons of the DEM simulations (solid black lines) and continuum model predictions (dashed grey lines) of the $c^{l}$ field at four time snapshots representing different stages of the segregation process: $\tilde {t}=4 \times 10^3$, $2 \times 10^4$, $1 \times 10^5$ and $4 \times 10^5$ in the sequence of top left, top right, bottom left, bottom right. The third column shows comparisons of the quasi-steady, normalized velocity profiles at $\tilde {t}=4 \times 10^5$ from DEM simulations and continuum model predictions.

Figure 9. Comparisons of continuum model predictions with corresponding DEM simulation results for the transient evolution of the segregation dynamics for two cases of vertical chute flow of disks. (a) More large grains case $\{W/\bar {d}_0 = 60,\mu _{w}=0.45,c^{l}_0=0.75,d^{l}/d^{ s}=1.5\}$. (b) Larger size ratio case $\{W/\bar {d}_0 = 60,\mu _{w}=0.45,c^{l}_0=0.5,d^{l}/d^{s}=3.0\}$. Additional cases are shown in figure 8. Results are organized as described in the caption of figure 8.

Regarding flow fields, comparisons of the quasi-steady, normalized velocity fields at $\tilde {t} = 4 \times 10^5$ from the DEM simulations and the continuum model predictions are shown in the third columns of figures 8 and 9. Since the velocity field evolves minimally during the segregation process, only the flow field at long time is shown. We note that the velocity field is well predicted in all cases, including the creeping regions far from the wall. This favourable comparison stems from the ability of the NGF model to transition seamlessly between the dense inertial flow regime, such as the region near the walls in vertical chute flow where $\mu >\mu _{s}$, and the quasi-static flow regime, such as the central region of the chute where $\mu <\mu _{s}$. Further, the good agreement provides support for the assumptions underlying our generalization of the NGF model to bidisperse granular systems discussed in § 2.3 – in particular, the choices to use the average grain size $\bar {d}$ in the expression for the cooperativity length (2.11a) and to continue to use the numerical value for the non-local amplitude determined for monodisperse systems, $A=0.90$, without refitting. We note that since the diffusion and segregation fluxes in our model depend strongly on the flow kinematics through the strain rate and its gradient, respectively, accurate predictions of flow kinematics obtained using the generalized NGF model are crucial to successfully capturing the dynamics of the concentration field using the segregation model (5.5). To illustrate this point, the results of figure 8(a) are compared with predictions of the segregation model when coupled with a regularized version of the $\mu (I)$ rheology rather than the NGF model in Appendix C.

Next, we make comparisons between continuum model predictions and DEM data for vertical chute flow of bidisperse spheres. The governing equations and boundary conditions are the same as those used above for bidisperse disks. That is, the fluidity field is governed by the non-local rheology (5.4) with Dirichlet fluidity boundary conditions at the walls ($g=g_{loc}(\mu _{w}, P_{w})$ at $x=\pm W/2$), and the concentration field is governed by the segregation dynamics equation (5.5) with no-flux boundary conditions at the walls ($w_x^{l}=0$ at $x=\pm W/2$). The only difference from the process described above for bidisperse disks is that the values $P_{w}$ and $\mu _{w}$ used in the continuum simulations are obtained from the coarse-grained stress fields in the DEM data for each case, rather than based on the nominal value of the compressive wall stress $P_{w}$ applied in the DEM simulation. This is done to account for the normal stress differences that arise for spheres, which slightly affect the predicted velocity fields and hence the consequent concentration fields. We consider three different cases: (i) the base case $\{W/\bar {d}_0 = 60,\mu _{w}=0.51,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$, (ii) the more large grains and higher flow rate case $\{W/\bar {d}_0 = 60,\mu _{w}=0.58,c^{l}_0=0.75,d^{l}/d^{s}=1.5\}$ and (iii) the larger size ratio and higher flow rate case $\{W/\bar {d}_0 = 60,\mu _{w}=0.58,c^l_0=0.5,d^l/d^s=2.0\}$, which are shown in figures 10(a), 10(b), and 10(c), respectively. The first column of figure 10 shows the spatio-temporal contours of the evolution of the $c^{l}$ field from the DEM data. The second column shows comparisons between the DEM simulations (solid black lines) and the continuum model predictions (dashed grey lines) for the $c^{l}$ field at four snapshots in time ($\tilde {t} = t/( d^{s}\sqrt {\rho _{s}/P_{w}} ) =5 \times 10^2$, $2 \times 10^3$, $1 \times 10^4$ and $4.5 \times 10^4$), indicated by the horizontal lines in the first column of figure 10. Lastly, comparisons of the quasi-steady, normalized velocity fields at $\tilde {t} = 4.5 \times 10^4$ from the DEM simulations and the continuum model predictions are shown in the third column of figure 10. The continuum model is able to capture the decaying velocity field quite well in all cases, and therefore, our choice to continue using the value of the non-local amplitude estimated for monodisperse systems, $A=0.43$, without readjustment works well for spheres. Overall, the coupled continuum model is capable of quantitatively predicting both the flow fields and the transient evolution of the segregation dynamics in vertical chute flow of bidisperse spheres.

Figure 10. Comparisons of continuum model predictions with corresponding DEM simulation results for the transient evolution of the segregation dynamics for three cases of vertical chute flow of spheres. (a) Base case $\{W/\bar {d}_0 = 60,\mu _{w}=0.51,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$. (b) More large grains and higher flow rate case $\{W/\bar {d}_0 = 60,\mu _{w}=0.58,c^{l}_0=0.75,d^{l}/d^{s}=1.5\}$. (c) Larger grain-size ratio and higher flow rate case $\{W/\bar {d}_0 = 60,\mu _{w}=0.58,c^{l}_0=0.5,d^{l}/d^{s}=2.0\}$. Results are organized as described in the caption of figure 8.

5.2. Annular shear flow

We utilize an analogous process to obtain continuum model predictions for the transient evolution of concentration and flow fields in annular shear flow of disks. In annular shear, based on the static force and moment balances, the pressure field is uniform, $P(r)=P_{w}$, and the stress ratio field is given by (4.3). The governing equations (5.4) and (5.5) are modified to appropriately account for the divergence and Laplacian operators in cylindrical coordinates. Also, since $v_{w}$, not $\mu _{w}$, is specified in our DEM simulations of annular shear, while $\mu _{w}$ is specified in our continuum simulations, we iteratively adjust the value of $\mu _{w}$ input into our continuum simulations in order to achieve the target value of $v_{w}$ in the predicted quasi-steady flow field. Otherwise, our process for obtaining numerical predictions from the continuum model is the same. Dirichlet fluidity boundary conditions and no-flux boundary conditions are imposed at the walls, and the initial concentration field is extracted from the initial DEM configuration for each case.

Then, we compare continuum model predictions against DEM data for the five cases of annular shear flow discussed in § 4.2 in order to further validate the model. Figures 11 and 12 summarize the comparisons for these fives cases and are organized in the same manner as figures 8 and 9. Again, the coupled, continuum model does a good job capturing the segregation dynamics and its dependence on the input parameters. Moreover, the quasi-steady velocity fields are well predicted by the NGF model in all cases, including the quasi-static, creeping region far from the inner wall, which is key to accurately capturing the evolution of the concentration field. We reiterate that all continuum model predictions are obtained using the same set of material parameters for disks (5.1). As for vertical chute flow, to illustrate the key role played by the NGF model in capturing the segregation dynamics, the results of figure 11(a) for the base case of annular shear flow are compared with predictions of the segregation model when coupled with a local rheological model in Appendix C.

Figure 11. Comparisons of the continuum model predictions with corresponding DEM simulation results for the transient evolution of the segregation dynamics for three cases of annular shear flow of disks. (a) Base case $\{R/\bar {d}_0 = 60,\tilde {v}_{w}=0.01,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$. (b) Lower inner-wall velocity case $\{R/\bar {d}_0 = 60,\tilde {v}_{w}=0.001,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$. (c) Smaller annular shear cell case $\{R/\bar {d}_0 = 40,\tilde {v}_{w}=0.01,c^{l}_0=0.5,d^{l}/d^{s}=1.5\}$. Additional cases are shown in figure 12. For each case, the first column shows spatio-temporal contours of the evolution of $c^{l}$ measured in the DEM simulations. The second column shows comparisons of the DEM simulations (solid black lines) and continuum model predictions (dashed grey lines) of the $c^{l}$ field at four time snapshots representing different stages of the segregation process: $\tilde {t}=5$, $50$, $200$ and $550$ in the sequence of top left, top right, bottom left, bottom right. The third column shows comparisons of the quasi-steady, normalized velocity profiles at $\tilde {t}=550$ from DEM simulations and continuum model predictions.

Figure 12. Comparisons of the continuum model predictions with corresponding DEM simulation results for the transient evolution of the segregation dynamics for two cases of annular shear flow of disks. (a) More large grains case $\{R/\bar {d}_0 = 60,\tilde {v}_{w}=0.01,c^{l}_0=0.75,d^{l}/d^{s}=1.5\}$. (b) Larger size ratio case $\{R/\bar {d}_0 = 60,\tilde {v}_{w}=0.01,c^{l}_0=0.5,d^{l}/d^{s}=3.0\}$. Additional cases are shown in figure 11. Results are organized as described in the caption of figure 11.