2.1. Bidisperse systems

We consider bidisperse granular systems made up of grains of two different sizes: large grains with average diameter $d^{l}$ and small grains with average diameter $d^{s}$. As in Liu et al. (Reference Liu, Singh and Henann2023), our focus is size-based segregation, and to eliminate density-based segregation, all grains are taken to be made of the same material with mass density $\rho _{s}$. Both three-dimensional systems of bidisperse spheres and two-dimensional systems of bidisperse disks are considered in this work but, for brevity, in what follows, we describe the continuum framework for three-dimensional systems of spheres. In the mixture theory approach utilized here, several quantities are introduced for each species: the large grains and the small grains. Large-grain quantities are denoted with a superscript ${l}$ and small-grain quantities with a superscript ${s}$. The solid fractions, i.e. the volumes occupied by each species per unit total volume, are denoted as $\phi ^{l}$ and $\phi ^{s}$ for large and small grains, respectively, and the total solid fraction is $\phi = \phi ^{l} + \phi ^{s}$. The concentrations of each species are defined as $c^{l} = \phi ^{l}/\phi$ and $c^{s} = \phi ^{s}/\phi ^{s}$, so that $c^{l} + c^{s} =1$. The average grain size is defined as $\bar {d} = c^{l} d^{l} + c^{s} d^{s}$.

2.2. Kinematics of flow

The velocity fields for each species are denoted as $v^{l}_i$ and $v^{s}_i$, and the mixture velocity field is defined as $v_i = c^{l} v^{l}_i + c^{s} v^{s}_i$. The strain-rate tensor for the mixture is defined as the symmetric part of the gradient of the mixture velocity: ${\mathsf{D}}_{ij} = (1/2)(\partial v_i / \partial x_j + \partial v_j/ \partial x_i )$. In this work, we focus on dense flows spanning the quasi-static and dense inertial flow regimes $(I\lesssim 10^{-1})$. Across our DEM simulations, we observe that the total solid fraction remains approximately uniform both spatially and temporally, and any volume change at flow initiation occurs over a much shorter time scale than the process of segregation. Therefore, we make the common idealization that dense flow of the mixture proceeds at constant volume (e.g. Gray & Thornton Reference Gray and Thornton2005; Gray & Chugunov Reference Gray and Chugunov2006; Gray & Ancey Reference Gray and Ancey2011; Gajjar & Gray Reference Gajjar and Gray2014; Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021; Liu et al. Reference Liu, Singh and Henann2023). (We note that this idealization is not appropriate for flows in the dilute, or collisional, flow regime ($I\gtrsim 10^{-1}$), in which the total solid fraction can evolve both spatially and temporally.) Under the constant-volume idealization, the mixture velocity field is divergence free, $\partial v_i/\partial x_i=0$; the strain-rate tensor is deviatoric, ${\mathsf{D}}_{kk}=0$; and the total solid fraction $\phi$ remains constant. Based on our DEM simulations, throughout this study, we take $\phi =0.6$ for spheres and $\phi =0.8$ for disks. The equivalent shear-strain rate is defined as $\dot {\gamma } = (2 {\mathsf{D}}_{ij} {\mathsf{D}}_{ij})^{1/2}$.

2.3. Balance of mass

Using the mixture velocity field and the species-specific velocity fields, the relative volume fluxes for large and small grains are defined as $w^{l}_i = c^{l} (v^{l}_i - v_i)$ and $w^{s}_i = c^{s} (v^{s}_i - v_i)$, respectively, so that $w^{l}_i + w^{s}_i = 0_i$. Conservation of mass applied to the large grains requires that

(2.1)\begin{equation} \frac{{\rm D} c^{l}}{{\rm D} t} + \frac{\partial w^{l}_i}{\partial x_i} = 0, \end{equation}

where ${\rm D}(\bullet )/{\rm D}t$ is the material time derivative. In what follows, we utilize $c^{l}$ as the field variable that describes the dynamics of size segregation. A constitutive equation is then required for the relative volume flux $w^{l}_i$. In our previous work (Liu et al. Reference Liu, Singh and Henann2023), we proposed a constitutive equation for $w^{l}_i$ that accounts for diffusion and shear-strain-rate-gradient-driven segregation. In § 2.5 we review this model and extend it to account for pressure-gradient-driven segregation.

2.4. Stress and the equations of motion

The stress-related fields for the mixture are defined as follows: the symmetric Cauchy stress tensor $\sigma _{ij} = \sigma _{ji}$, 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 equations of motion are

(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). Rheological constitutive equations are then required for the Cauchy stress $\sigma _{ij}$. In our previous work (Liu et al. Reference Liu, Singh and Henann2023), we generalized the NGF model for dense granular flow from monodisperse to bidisperse systems. In this work, we continue to utilize the generalized NGF model to describe the rheological behaviour of dense, bidisperse granular flows, which, for brevity, is summarized in Appendix A.

2.5. Segregation model

In this section we discuss the constitutive equation for the relative volume flux $w^{l}_i$. In dense granular flows a bidisperse mixture tends to segregate due to both shear-strain-rate gradients and pressure gradients. Moreover, diffusion acts counter to segregation and tends to mix the two species. Accordingly, we take the relative volume flux to be comprised of three contributions as follows:

(2.3)\begin{equation} w^{l}_i = w^{diff}_i + w^{S}_i + w^{P}_i. \end{equation}

Here $w^{diff}_i$ is the diffusion flux, $w^{S}_i$ is the shear-strain-rate-gradient-driven segregation flux and $w^{P}_i$ is the pressure-gradient-driven segregation flux.

The first two contributions have been examined in detail in our prior work (Liu et al. Reference Liu, Singh and Henann2023) and are briefly recapped here. First, we take the diffusion flux to be driven by concentration gradients as $w^{diff}_i = -D (\partial c^{l} / \partial x_i )$, where $D$ is the binary diffusion coefficient. In a dense, bidisperse granular mixture, the binary diffusion coefficient follows a well-established scaling (e.g. Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021; Duan et al. Reference Duan, Umbanhowar, Ottino and Lueptow2021; Trewhela et al. Reference Trewhela, Ancey and Gray2021; Liu et al. Reference Liu, Singh and Henann2023) with the average grain size $\bar {d}$ and the shear-strain rate $\dot {\gamma }$ as $D = C_{diff} \bar {d}^2 \dot {\gamma }$, where $C_{diff}$ is a dimensionless material parameter. Therefore, the diffusion flux is

(2.4)\begin{equation} w^{diff}_i ={-}C_{diff} \bar{d}^2 \dot{\gamma} \frac{\partial c^l}{\partial x_i}. \end{equation}

Second, in Liu et al. (Reference Liu, Singh and Henann2023) we isolated shear-strain-rate-gradient-driven segregation in two flow configurations with uniform pressure and proposed the following phenomenological constitutive equation for the flux $w^{S}_i$:

(2.5)\begin{equation} w^{S}_i = C^{S}_{seg} \bar{d}^2 c^{l}(1-c^{l}) \frac{\partial \dot{\gamma}}{\partial x_i}. \end{equation}

Here $C^{S}_{seg}$ is another dimensionless material parameter. In (2.5) the pre-factor depends on the average grain size through $\bar {d}^2$ for dimensional consistency and involves a symmetric dependence on $c^{l}$ through $c^{l}(1-c^{l})$, which ensures that segregation ceases when the bidisperse mixture becomes either all large ($c^{l}=1$) or all small ($c^{l}=0$) grains. While a symmetric dependence of the pre-factor on $c^{l}$ was sufficient to capture the DEM data of Liu et al. (Reference Liu, Singh and Henann2023) in the absence of pressure gradients, it contrasts with the asymmetric dependence that has been invoked to model pressure-gradient-driven size segregation in the literature (e.g. Gajjar & Gray Reference Gajjar and Gray2014; Tunuguntla, Weinhart & Thornton Reference Tunuguntla, Weinhart and Thornton2017; Jones et al. Reference Jones, Isner, Xiao, Ottino, Umbanhowar and Lueptow2018; Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021; Trewhela et al. Reference Trewhela, Ancey and Gray2021) and in the next paragraph of the present work.

Next, we consider the segregation flux associated with pressure gradients, $w^{P}_i$. We hypothesize that the segregation flux $w^{P}_i$ is driven by gradients in the pressure $P$ and adopt the following phenomenological form for the constitutive equation for $w^{P}_i$:

(2.6)\begin{equation} w^{P}_i ={-}C^{P}_{seg} \frac{\bar{d}^2 \dot{\gamma}}{P} c^{l}(1-c^{l}) (1- \alpha + \alpha c^{l}) \frac{\partial P}{\partial x_i}. \end{equation}

Here $C^{P}_{seg}$ and $\alpha$ are dimensionless parameters dependent upon mixture properties (e.g. the grain-size ratio). The dependence of the pre-factor on $\dot {\gamma }c^{l}(1-c^{l})$ ensures that (2.6) satisfies several minimal requirements – namely, that the pressure-gradient-driven segregation flux is zero when there is no flow ($\dot \gamma =0$) or when the mixture becomes either all large ($c^{l}=1$) or all small ($c^{l}=0$) grains. Next, the dependence of the pre-factor on $\bar {d}^2/P$ ensures dimensional consistency. Finally, in contrast to (2.5), the dependence of (2.6) on the factor $(1 - \alpha + \alpha c^{l})$ introduces an asymmetric dependence on $c^{l}$ into the flux constitutive equation.

Several comments on the flux constitutive equation (2.6) are warranted.

1. (i) The dependence of (2.6) on $c^{l}$ through the asymmetric flux function $f(c^{l}) = c^{l}(1-c^{l}) (1 - \alpha + \alpha c^{l})$ follows directly from the work of Gajjar & Gray (Reference Gajjar and Gray2014), where $\alpha \in [0,1]$ is a parameter that controls the amount of asymmetry (denoted as $\gamma$ in Gajjar & Gray Reference Gajjar and Gray2014). As discussed in Gajjar & Gray (Reference Gajjar and Gray2014), the flux function $f(c^{l})$ has the following properties: (1) for $\alpha =0$, the symmetric flux function is recovered; (2) for $\alpha \in (0,1]$, the function's maximum is skewed from $c^{l}=0.5$ towards $c^{l}=1$; (3) for $\alpha \in (0, 0.5]$, it is convex; and (4) for $\alpha \in (0.5, 1]$, it is non-convex with a single inflection point. Previously, a value of $\alpha =0.46$ was obtained in Gajjar & Gray (Reference Gajjar and Gray2014) by fitting to the experiments of Wiederseiner et al. (Reference Wiederseiner, Andreini, Épely-Chauvin, Moser, Monnereau, Gray and Ancey2011), and a value of $\alpha =0.89$ was determined by fitting to experiments in van der Vaart et al. (Reference van der Vaart, Gajjar, Epely-Chauvin, Andreini, Gray and Ancey2015). In what follows, we estimate $\alpha$ along with $C^{P}_{seg}$ for both simulated, frictional spheres and disks by fitting to DEM simulations.

2. (ii) We note that the scaling of the pressure-gradient-driven segregation flux (2.6) with $\bar {d}^2\dot {\gamma }/P$ is quite similar to the scalings for the segregation velocity reported in Trewhela et al. (Reference Trewhela, Ancey and Gray2021) and Jing et al. (Reference Jing, Ottino, Umbanhowar and Lueptow2022) based on experiments and discrete simulations, respectively, of the dynamics of a single intruder in a flowing dense granular medium, and the scaling of Trewhela et al. (Reference Trewhela, Ancey and Gray2021) was utilized in the continuum simulations of Barker et al. (Reference Barker, Rauter, Maguire, Johnson and Gray2021). As pointed out in these works, combining a segregation flux that is inversely proportional to the pressure with a pressure-independent diffusion flux enables a model to capture the attenuation of segregation with increasing pressure, which was observed in the discrete simulations of Fry et al. (Reference Fry, Umbanhowar, Ottino and Lueptow2018, Reference Fry, Umbanhowar, Ottino and Lueptow2019) in flows involving weak strain-rate gradients. However, one difference is that Trewhela et al. (Reference Trewhela, Ancey and Gray2021) introduces a second term in the denominator to prevent a singularity when the pressure is equal to zero, such as at a free surface. This additional term in the denominator depends on the product of the magnitude of the pressure gradient and the average grain size $\bar {d}$ and is multiplied by an additional (small) dimensionless constant. Here, we follow a similar approach but directly add a small constant to the pressure field when dealing with free-surface flows in § 4.1.

3. (iii) Finally, it has been well established that the pressure-gradient-driven segregation flux should depend on the grain-size ratio $d^{l}/d^{s}$ (e.g. Schlick et al. Reference Schlick, Fan, Isner, Umbanhowar, Ottino and Lueptow2015; Tunuguntla et al. Reference Tunuguntla, Weinhart and Thornton2017; Jones et al. Reference Jones, Isner, Xiao, Ottino, Umbanhowar and Lueptow2018; Barker et al. Reference Barker, Rauter, Maguire, Johnson and Gray2021; Duan et al. Reference Duan, Umbanhowar, Ottino and Lueptow2021; Trewhela et al. Reference Trewhela, Ancey and Gray2021). However, the aims of this paper are to combine pressure-gradient-driven and shear-strain-rate-gradient-driven segregation fluxes within a coupled model for size segregation and flow and then to use this model to examine the interplay between the mechanisms. Therefore, we consider a single grain-size ratio of $d^{l}/d^{s}=1.5$ throughout, and considering the dependence of (2.6) on $d^{l}/d^{s}$ is beyond the scope of this paper. We expect that the substantial progress in the literature to characterize the role of the grain-size ratio (e.g. Trewhela et al. Reference Trewhela, Ancey and Gray2021) may be incorporated into (2.6).

Combining (2.4), (2.5) and (2.6) with the balance of mass equation (2.1), we obtain the following differential relation governing the dynamics of $c^{l}$:

(2.7) $$\begin{gather} \dfrac{{\rm D}{c}^{l}}{{\rm D}t} + \dfrac{\partial}{\partial x_i} \left({-}C_{diff} \bar{d}^2 \dot{\gamma} \frac{\partial c^l}{\partial x_i} + C^{S}_{seg} \bar{d}^2 c^{l}(1-c^{l}) \frac{\partial \dot{\gamma}}{\partial x_i} \right.\nonumber\\ \left.- C^{P}_{seg} \frac{\bar{d}^2 \dot{\gamma}}{P} c^{l}(1-c^{l}) (1- \alpha + \alpha c^{l}) \frac{\partial P}{\partial x_i} \right) = 0. \end{gather}$$

The material parameters associated with the segregation model are the set $\{C_{diff}, C^{S}_{seg}, C^{P}_{seg}, \alpha \}$. As discussed in Liu et al. (Reference Liu, Singh and Henann2023), the parameter $C_{diff}$ may be determined from measurements of the mean square displacement during homogeneous simple shearing, and the parameter $C^{S}_{seg}$ may be determined by examining size segregation in flows, in which the pressure field is spatially uniform. Based on the DEM simulations of Liu et al. (Reference Liu, Singh and Henann2023), these parameters are $\{C_{diff}=0.045, C^{S}_{seg}=0.08\}$ for frictional spheres and $\{C_{diff}=0.20, C^{S}_{seg}=0.23\}$ for frictional disks over a range of modest grain-size ratios including $d^{l}/d^{s}=1.5$. The remaining material parameters $\{C^{P}_{seg}, \alpha \}$, which are associated with the pressure-gradient-driven segregation flux, will be determined for frictional spheres and disks in the following section.

4.1. Inclined plane flow

We first consider inclined plane flow to test predictions of the coupled continuum model against corresponding DEM results. As discussed above, when there are no normal stress differences, the stress field may be obtained from a static force balance, which implies that the stress ratio field is uniform and given by $\mu (z) = \tan \theta$ and that the pressure field is linear in $z$ and given by $P(z) = \phi \rho _{s} G z \cos \theta$. However, the normal stress differences that arise in dense flows of spheres induce a slightly higher, uniform value of $\mu$ and a slightly lower slope in the pressure field $P(z)$. Accordingly, to control for this effect when working with dense flows of spheres, in our continuum simulations, we utilize the values of the uniform stress ratio and the slope of the pressure field obtained from the coarse-grained stress fields in the DEM data for each case, rather than the nominal values of $\tan \theta$ and $\phi \rho _{s} G \cos \theta$, respectively. Moreover, to avoid a singularity in the pressure-gradient-driven segregation flux (2.6) at the free surface where $z=0$, we add a small constant to the pressure field corresponding to the weight of a layer of $(1/4)\bar {d}_0$ thickness, i.e. $(1/4)\phi \rho _{s} G \bar {d}_0 \cos \theta$. This approach is quite similar to that of Trewhela et al. (Reference Trewhela, Ancey and Gray2021), who directly incorporate this small constant into the denominator of their pressure-gradient-driven flux equation. We have verified that the subsequently presented results are insensitive to the exact choice of this constant, so long as it is sufficiently small. With the relevant stress-related fields determined in this way, 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. The remaining unknown fields in inclined plane flow are the velocity field $v_x(z,t)$ with the associated strain-rate field $\dot {\gamma }(z,t) = \partial v_x/\partial z$, the granular fluidity field $g(z,t)$ and the concentration field $c^{l}(z,t)$. The governing equations are (1) the flow rule (A2),

(4.2)\begin{equation} \dot{\gamma} = g \mu; \end{equation}

(2) the non-local rheology (A3),

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

where $g_{loc} (\mu, P)$ and $\xi (\mu )$ are stress-dependent functions given through (A4) and (A5)$_1$, respectively; and finally (3) the segregation dynamics equation (2.7),

(4.4)$$\begin{gather} \dfrac{\partial c^{l}}{\partial t} + \frac{\partial}{\partial z} \left(- C_{diff} \bar{d}^2\dot{\gamma} \frac{\partial c^{l}}{\partial z} + C^{S}_{seg}\bar{d}^2c^{l}(1-c^{l})\frac{\partial \dot{\gamma}}{\partial z} \right.\nonumber\\ \left. - C^{P}_{seg} \frac{\bar{d}^2 \dot{\gamma}}{P} c^{l}(1-c^{l})(1 - \alpha + \alpha c^{l}) \frac{\partial P}{\partial z} \right) = 0, \end{gather}$$

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

For the fluidity boundary conditions, we impose a homogeneous Neumann boundary condition at the free surface, i.e. $\partial g/ \partial z =0$ at $z=0$, due to the similarity of a flat free surface to a symmetry plane. The bottom surface is comprised of touching, glued grains in our DEM simulations, and the selection of a fluidity boundary condition for this type of boundary is discussed in our prior work (Liu et al. Reference Liu, Singh and Henann2023). Following this work, we impose a Dirichlet fluidity boundary condition at the bottom surface wherein the fluidity is a function of the stress state through the local fluidity function, i.e. $g=g_{loc}(\mu (z=H), P(z=H))$ at $z=H$. For the segregation dynamics equation (4.4), we apply no flux boundary conditions at both the free surface and the bottom surface, i.e. $w^{l}_z = - C_{diff} \bar {d}^2 \dot {\gamma } (\partial c^{l}/ \partial z ) + C^{S}_{seg} \bar {d}^2 c^{l} (1- c^{l})( \partial \dot {\gamma } / \partial z ) - C^{P}_{seg} ({\bar {d}^2 \dot {\gamma }}/{P})c^{l} (1- c^{l}) (1 - \alpha + \alpha c^{l}) (\partial P/ \partial z) =0$ at $z=0$ and $H$. Finally, an initial condition for the segregation dynamics equation (4.4) – i.e. $c^{l}_0(z) = c^{l}(z, t=0)$ – is necessary. To account for spatial fluctuations in $c^{l}_0(z)$, we coarse grain the concentration field in the DEM configuration at $t=0$ for each case and use this field as the initial condition in the corresponding continuum simulation.

Then, we obtain numerical predictions of the continuum model for a given case of inclined plane flow using finite differences as follows. At a given point in time, the concentration field $c^{l}(z)$ is known, and thus, the average grain size may be calculated as $\bar {d}(z) = c^{l}(z)d^{l} + (1-c^{l}(z))d^{s}$. Since the pressure distribution $P(z)$ and stress ratio distribution $\mu (z)$ are known, the local fluidity $g_{loc}(\mu, P)$ and the cooperativity length $\xi (\mu )$, given through (A4) and (A5)$_1$, respectively, can also be calculated at the spatial grid points. The non-local rheology equation (4.3) can then be solved by discretizing the Laplacian term $\partial ^2 g /\partial z^2$ using central differences in space. Once the fluidity field has been calculated at all grid points, the strain-rate field can be calculated using (4.2), and the velocity field $v_x(z)$ can be obtained by integrating the strain-rate field. Next, the segregation dynamics equation (4.4) is used to update the concentration field by (1) using the Euler method for temporal discretization, (2) treating both of the segregation fluxes explicitly, and (3) treating the diffusion term implicitly. We have verified that both the spatial and temporal resolutions employed are sufficiently refined, so as to ensure the numerical accuracy and stability of our finite-difference scheme. With this, one step of numerical integration is completed and the unknown fields $c^{l}$, hence $\bar {d}$, are obtained at the next time step. The same procedure is repeated over multiple time steps to calculate the transient evolution of the concentration and flow fields, $c^{l}(z,t)$ and $v_x(z,t)$, over the desired simulation time window.

We compare model predictions of the transient evolution of the concentration field and the steady-state flow field against corresponding DEM simulation results for all four cases of inclined plane flow of dense, bidisperse spheres considered in § 3. The comparisons are summarized in figures 3 and 4. For each case, the spatiotemporal evolution of the DEM-measured $c^{l}$ field is shown in the first column of figure 3 or figure 4, and the second columns show comparisons of snapshots of the $c^{l}$ field measured in DEM simulations against corresponding continuum model predictions at four different time instants – $\tilde {t} = t/\sqrt {\bar {d}_0/G} = 200, 2000, 10\,000$, and $15\,000$ – as indicated by the vertical lines on the contour plots in the first columns. The model does a good job predicting the segregation dynamics across different cases of inclined plane flow and can capture the variations in the evolution of the $c^{l}$ field as the input parameters are changed. The third columns show comparisons of the steady-state velocity field predicted by the continuum model with the corresponding DEM-measured velocity fields. The Bagnold-like profile is captured well in all cases. We note that in inclined plane flow, the local inertial, or $\mu (I)$, rheology described in Appendix A.1 would be sufficient to predict the flow fields, and a non-local rheological approach is not necessary for these flows. However, as shown in our prior work (Liu et al. Reference Liu, Singh and Henann2023) and in the subsequent section on planar shear flow with gravity, non-local rheological modelling is crucial in many other flow configurations. One benefit of the NGF model is that it automatically reduces to the local inertial rheology when appropriate, so for consistency, we have continued to use the NGF model when considering inclined plane flow. Finally, we have also tested predictions of the coupled continuum model against corresponding DEM results for dense, bidisperse flows of frictional disks, and the model does an equally good job simultaneously capturing the segregation dynamics and the steady-state flow fields for disks as for spheres.

Figure 3. Comparisons of continuum model predictions with corresponding DEM simulation results for the transient evolution of the segregation dynamics for two cases of inclined plane flow of bidisperse spheres: (a) base case $\{H/\bar {d}_0=50, \theta =26^\circ, c^{l}_0 = 0.50, d^{l}/d^{s} = 1.5\}$ and (b) lower inclination angle case $\{H/\bar {d}_0=50, \theta =24^\circ, c^{l}_0 = 0.50, d^{l}/d^{s} = 1.5\}$. For each case, plots (a i,b i) show spatiotemporal contours of the evolution of $c^{l}$ measured in DEM simulations. Plots (a ii–v,b ii–v) show comparisons of the DEM simulations (solid black lines) and continuum model predictions (dashed grey lines) of the $c^{l}$ field at four different time instants during the segregation process: $\tilde {t} = 200$, $2000$, $10\,000$ and $15\,000$ in the sequence of (a ii,b ii), (a iii,b iii), (a iv,b iv), (a v,b v). Plots (a vi,b vi) show comparisons of the steady-state, normalized velocity profiles at $\tilde {t} = 15\,000$ from DEM simulations and continuum model predictions.

Figure 4. Comparisons of continuum model predictions with corresponding DEM simulation results for the transient evolution of the segregation dynamics for two cases of inclined plane flow of bidisperse spheres: (a) more large grains case $\{H/\bar {d}_0=50, \theta =26^\circ, c^{l}_0 = 0.75, d^{l}/d^{s} = 1.5\}$ and (b) thicker layer case $\{H/\bar {d}_0=60, \theta =26^\circ, c^{l}_0 = 0.50, d^{l}/d^{s} = 1.5\}$. Results are organized as described in the caption of figure 3.

4.2. Planar shear flow with gravity

We have hypothesized that the proposed constitutive equation (2.6) for the pressure-gradient-driven segregation flux and the dimensionless material parameters $\{C^{P}_{seg},\alpha \}$ are general across different flow geometries. To test this hypothesis, we consider a new flow configuration: planar shear flow with gravity acting orthogonal to the shearing direction, as illustrated in figure 5(a). We did not use DEM data from flows in this configuration to estimate the parameters $\{C^{P}_{seg},\alpha \}$ and do no further parameter adjustment in this section, so comparisons of continuum model predictions of the transient evolution of the large-grain concentration field and the steady-state flow field against DEM simulation results may be regarded as independent validation tests of the model.

Figure 5. Representative base case for planar shear flow with gravity. (a) Initial well-mixed configuration for three-dimensional DEM simulation of a bidisperse mixture of spheres with $\sim$15 000 flowing grains. The thickness of the flowing layer is $H=50 \bar {d}_0$. Black grains represent the rough top and bottom walls. Dark grey grains indicate large flowing grains and light grey grains indicate small flowing grains. (b) Segregated state after flowing for a simulation time of $\tilde {t} = t/(\ell /v_{w}) = 600$. (c) Spatiotemporal evolution of the large-grain concentration field. Spatial profiles of (d) the concentration field $c^{l}$ and (e) the normalized velocity field $v_x/v_{w}$ at three times ($\tilde {t} = 10$, $100$ and $600$) as indicated by the vertical lines in (c).

We first describe this flow configuration in detail and discuss its important characteristics. Consider a semi-infinite layer of a dense, bidisperse granular mixture between two rough parallel walls, separated by a distance $H$ along the $z$ direction. For the case of bidisperse spheres, the DEM set-up is shown in figure 5(a) for $H=50 \bar {d}_0$. The top and bottom walls, indicated as black grains in figure 5(a), are made of layers of touching, glued grains. The top wall imposes a compressive normal stress $P_{w}$ on the flowing grains along the $z$ direction and is specified to move with a constant velocity $v_{w}$ along the $x$ direction. The compressive normal stress applied by the top wall is maintained at its target value $P_{w}$ using the feedback scheme described elsewhere in the literature (e.g. 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). The bottom wall is specified to remain fixed and imposes a no-slip boundary condition at the bottom. The gravitational body force with acceleration of gravity $G$ acts along the $z$ direction. Furthermore, we employ periodic boundary conditions along the flow ($x$) and lateral ($y$) directions, eliminating any lateral wall/boundary effects and leading to a simple one-dimensional flow in which the only non-zero velocity field component is $v_x$ that depends only on the $z$ coordinate. In all DEM simulations of spheres, we take the length of the simulation domain along the $x$ direction to be $L = 20 \bar {d}_0$ and the length along the $y$ direction to be $W = 10 \bar {d}_0$, as shown in figure 5(a).

Regarding the stress field, the force balance along the $x$ direction gives that the shear stress magnitude is spatially uniform and equal to the shear stress imparted by the moving wall $\tau _{w}$, i.e. $|\sigma _{xz}(z)| = |\sigma _{zx}(z)| = \tau _{w}$. We note that instead of directly prescribing $\tau _{w}$ in our DEM simulations, the top-wall velocity $v_{w}$ is prescribed, and the shear stress $\tau _{w}$ arises as an output. The force balance along the $z$ direction gives that $\sigma _{zz}(z) = -(P_{w} + \phi \rho _{s} G z)$. As discussed in § 3, we observe that $\sigma _{xx}(z) \approx \sigma _{zz}(z)$ in DEM simulations of two-dimensional systems of disks, so that $\tau (z) = |\sigma _{xz}(z)| = \tau _{w}$, $P(z) = -\sigma _{zz} = P_{w} + \phi \rho _{s} G z$ and $\mu (z) = \tau (z)/P(z) = {\mu _{w}}/{(1 + z/\ell )}$, where $\mu _{w} = \tau _{w}/P_{w}$ is the maximum value of the stress ratio occurring at the top wall ($z=0$) and $\ell =P_{w}/\phi \rho _{s} G$ is the loading length scale. The loading length scale is defined as the ratio of the top-wall pressure $P_{w}$ to the gravitational body force $\phi \rho _{s} G$, which may be interpreted as the height of a granular layer that applies a pressure due to its weight equivalent to $P_{w}$, and is a distinct length scale from the dimensions $H$, $L$ and $W$. Again, while the coarse-grained stress fields in DEM simulations of disks are quite close to these analytical fields determined from force balances, in DEM simulations of spheres we observe normal stress differences, and while this leads to slight differences in the coarse-grained $\mu (z)$ and $P(z)$ fields for spheres, the $z$ dependence of these fields is consistent with the expressions derived above. Since the stress ratio is greatest at the top wall ($z=0$) and decays with $z$, we expect there to be a flowing region just beneath the top wall with the velocity field $v_x(z)$ decaying as $z$ increases, as qualitatively sketched in figure 5(a). Importantly, the loading length scale $\ell$ is the only length scale appearing in the analytical expression for the stress ratio field $\mu (z) = {\mu _{w}}/{(1 + z/\ell )}$. Therefore, $\ell$ and not the layer thickness $H$ is the relevant length scale in planar shear flow with gravity and affects the characteristic size of the flowing region beneath the top wall. We have verified that all dimensions of the simulation domain $\{H=50 \bar {d}_0, L = 20\bar {d}_0, W = 10\bar {d}_0\}$ are sufficiently large, so that they do not affect the flow fields or the segregation dynamics.

The important dimensionless parameters that describe planar shear flow with gravity are (1) ${\ell }/{\bar {d}_0}=P_{w}/\phi \rho _{s}G\bar {d}_0$, the dimensionless loading length scale; (2) $\tilde {v}_{w}=({v_{w}}/{\ell }) \sqrt {{\rho _{s} \bar {d}_0^2}/{P_{w}}}$, the dimensionless top-wall velocity that determines shear stress $\tau _{w}$ and the stress ratio $\mu _{w}$ at the top wall; (3) $c^{l}_0(z)$, the initial large-grain concentration field; and (4) $d^{l}/d^{s}$, the bidisperse grain-size ratio. Thus, the parameter set $\{\ell /\bar {d}_0, \tilde {v}_{w}, c^{l}_0, d^{l}/d^{s}\}$ fully specifies a given case of planar shear flow with gravity. We consider a representative base case for spheres corresponding to the parameter set $\{\ell /\bar {d}_0=18,\tilde {v}_{w}=0.02, c^{l}_0=0.50, d^{l}/d^{s} = 1.5\}$. In this flow configuration the strain rate is greatest just under the top wall, where the pressure is lowest. Therefore, the shear-strain-rate-gradient-driven flux and the pressure-gradient-driven flux cooperate and drive large grains towards the top of the layer. The segregated state after running the base-case DEM simulation to a time of $\tilde {t} = t/(\ell /v_{w}) = 600$ is shown in figure 5(b), illustrating that the large grains (dark grey) indeed segregate towards the top wall, leaving a small-grain-rich region (light grey) underneath. Further away from the wall, the evolution of segregation is very slow due to the low strain rates in that region, and the bidisperse granular material remains well mixed. For a more quantitative picture, we coarse grain the concentration field $c^{l}$ in both space and time, and the spatiotemporal evolution of $c^{l}$ is shown in the contour plot in figure 5(c). Lastly, snapshots of the $c^{l}$ field and the non-dimensionalized velocity field $v_x/v_{w}$ at three different times – $\tilde {t} = 10, 100$ and $600$ – are shown in figures 5(d) and 5(e), respectively. Again, we observe that the velocity field quickly reaches a steady state, while the concentration field evolves over the full simulation time window.

Next, we solve the coupled continuum model for the transient evolution of the concentration and velocity fields in this geometry, following a process analogous to that described in § 4.1 for inclined plane flow. To control for the effects of normal stress differences in dense flows of spheres, the pressure at the top wall $P(z=0)$ and the slope of the pressure field are obtained from the coarse-grained pressure field in the DEM data for each case. These values are slightly different than the nominal values of $P_{w}$ and $\phi \rho _{s}G$ and give rise to a slightly adjusted value of the loading length scale, and we utilize these adjusted values in the stress fields, $\mu (z)$ and $P(z)$, in our continuum simulations. The governing equations are the flow rule (4.2), the non-local rheology (4.3) and the segregation dynamics equation (4.4), wherein $\mu (z)$ and $P(z)$ are given as inputs to the model. For the boundary conditions, we impose Dirichlet fluidity boundary conditions at both the top and bottom walls, i.e. $g=g_{loc}(\mu (z=0),P(z=0))$ at $z=0$ and $g=g_{loc}(\mu (z=H),P(z=H))$ at $z=H$, and no-flux boundary conditions at the top and bottom walls, i.e. $w^{l}_z=0$ at $z=0$ and $z=H$. The initial concentration field $c^{l}_0(z)$ is extracted from the initial DEM configuration for each case, and we use the same finite-difference-based numerical approach described above for inclined plane flow. Finally, since $v_{w}$ is prescribed in the DEM simulations of planar shear flow with gravity, while $\mu _{w}$ is specified in the corresponding continuum simulations, we iteratively adjust the value of $\mu _{w}$ input into the continuum simulations to match the target value of $v_{w}$ in the predicted steady-state velocity field.

We compare continuum model predictions against corresponding DEM simulation results for bidisperse spheres using the parameter set given in (4.1). To broadly exercise the model, we consider four different cases of planar shear flow with gravity: (1) the base case $\{\ell /\bar {d}_0=18,\,\tilde {v}_{w}=0.02, \, c^{l}_0=0.50, \, d^{l}/d^{s} = 1.5\}$; (2) a more large grains case $\{\ell /\bar {d}_0=18, \tilde {v}_{w}=0.02, \, c^{l}_0=0.75, \, d^{l}/d^{s} = 1.5\}$; (3) a larger loading length scale and lower top-wall velocity case $\{\ell /\bar {d}_0=36,\,\tilde {v}_{w}=0.01 c^{l}_0=0.50, \, d^{l}/d^{s} = 1.5\}$; and (4) a lower top-wall velocity case $\{\ell /\bar {d}_0=18,\,\tilde {v}_{w}=0.01 \, c^{l}_0=0.50, \, d^{l}/d^{s} = 1.5\}$, and results for these cases are shown in figures 6(a,b) and 7(a,b), respectively. The first columns show the spatiotemporal evolution of the concentration field $c^{l}$ measured from the DEM simulations. In the second columns, we compare snapshots of the concentration field $c^{l}$ predicted by the continuum model (dashed grey lines) with the corresponding coarse-grained DEM fields (solid black lines) at four different time instants – $\tilde {t} = t/(\ell /v_{w}) = 20, 100, 300$ and $600$ in figure 6 and $\tilde {t} = 10, 50, 150$ and $300$ in figure 7 – as indicated by the vertical lines on the contour plots in the first columns for each case. Finally, the steady-state normalized velocity fields (corresponding to the last time instant listed above for each case) obtained from DEM simulations and the corresponding continuum model predictions are shown in the last columns of figures 6 and 7. The continuum model is able to capture the decaying velocity fields in all cases as well as predict the transient evolution of the segregation process. A similar quality of comparison was observed in analogous tests for dense, bidisperse disks in several cases of planar shear flow with gravity.

Figure 6. Comparisons of continuum model predictions with corresponding DEM simulation results for the transient evolution of the segregation dynamics for two cases of planar shear flow with gravity of bidisperse spheres: (a) base case $\{\ell /\bar {d}_0=18, \tilde {v}_{w}=0.02, c^{l}_0=0.50, d^{l}/d^{s} = 1.5\}$ and (b) more large grains case $\{\ell /\bar {d}_0=18, \tilde {v}_{w}=0.02, c^{l}_0=0.75, d^{l}/d^{s} = 1.5\}$. For each case, plots (a i,b i) show spatiotemporal contours of the evolution of $c^{l}$ measured in DEM simulations. Plots (a ii–v,b ii–v) show comparisons of the DEM simulations (solid black lines) and continuum model predictions (dashed grey lines) of the $c^{l}$ field at four different time instants during the segregation process: $\tilde {t} = t/(\ell /v_{w}) = 20, 100, 300$ and $600$ in the sequence of (a ii,b ii), (a iii,b iii), (a iv,b iv), (a v,b v). Plots (a vi,b vi) show comparisons of the steady-state, normalized velocity profiles at $\tilde {t}=600$ from DEM simulations and continuum model predictions.

Figure 7. Comparisons of continuum model predictions with corresponding DEM simulation results for the transient evolution of the segregation dynamics for two cases of planar shear flow with gravity of bidisperse spheres: (a) larger loading length scale and lower top-wall velocity case $\{\ell /\bar {d}_0=36,\tilde {v}_{w}=0.01, c^{l}_0=0.50, d^{l}/d^{s} = 1.5\}$ and (b) lower top-wall velocity case $\{\ell /\bar {d}_0=18, \tilde {v}_{w}=0.01, c^{l}_0=0.50, d^{l}/d^{s} = 1.5\}$. Results are organized as described in the caption of figure 6. The time instants shown in (a ii–v,b ii–v) correspond to $\tilde {t} = 10, 50, 150$ and $300$, and the time instant shown in (a vi,b vi) corresponds to $\tilde {t}=300$.

At this stage, it is instructive to contrast features of planar shear flow with gravity with inclined plane flow to emphasize why this flow configuration represents a non-trivial validation test for a coupled model for size segregation and flow. At a high level, planar shear flow with gravity is boundary driven, while inclined plane flow is gravity driven. This difference in the manner in which flow is driven manifests in stark differences in the consequent flow fields. Namely, the Bagnold-like profile of inclined plane flow involves only modest strain-rate gradients, leading to segregation phenomenology that is dominated by pressure gradients. Moreover, the rheological behaviour of a dense granular medium in inclined plane flow is dominated by local grain-inertia effects that are captured by the local inertial, or $\mu (I)$, rheology. This is in contrast to planar shear flow with gravity, in which the flow fields decay rapidly with the distance from the moving, top wall. While local grain-inertia effects contribute in the flowing region under the top wall, non-local, cooperative effects also contribute and become increasingly dominant with the distance from the top wall, as the flow field transitions into the quasi-static, creeping region. Therefore, utilizing a non-local rheology, such as the generalized NGF model for bidisperse flows, is crucial to accurately capture the velocity fields in planar shear flow with gravity. It is also notable that the magnitudes of the strain-rate gradients are comparatively greater in planar shear flow with gravity than in inclined plane flow, and the comparative contributions of the two mechanisms of segregation, i.e. pressure gradients and shear-strain-rate gradients, are unclear. We return to this point in § 5. Due to the distinctly different flow fields and the relative directions of the two segregation fluxes, i.e. competitive versus cooperative, inclined plane flow and planar shear flow with gravity represent markedly different flow configurations, so the observation that the coupled continuum model is capable of capturing segregation dynamics and flow fields in both configurations with a single set of parameters is encouraging.