Shear-induced diffusion in non-local granular flows

Shearing dense granular flows induces diffusion of grains. This mechanism underpins the rate of mixing [1], heat transfer [2] and segregation [3] in a variety of natural and industrial granular flows. It is usually modelled by a coefficient of self-diffusion, also called diffusivity D [m²/s].

In homogeneous shear flows, in which there is no gradient of shear rate, three relationships have been established from which diffusivity can be predicted:

$$\begin{eqnarray} D &\propto \delta v d, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{\delta v}{\dot \gamma d} &\propto I^{-\frac{1}{2}}; \quad I = \dot \gamma t_i, \quad t_i = d \sqrt{\rho/\sigma}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} bI &= \mu-\mu_s, \text{ for } \mu > \mu_s; \quad I = 0, \text{ otherwise.} \end{eqnarray} \tag{ 3 }$$

The first scaling relates the diffusivity to the velocity fluctuations $\delta v$ and grain size d [4–6]. It is consistent with a typical grain trajectory following a random walk of step d and frequency $\delta v/d$. The second scaling relates the velocity fluctuations to the inertial number I, which compares the shear rate $\dot \gamma$ to an inertial time t_i involving the normal stress σ in the flow, and the grain size d and density ρ. It is consistent with the development of clusters of jammed grains with size $\ell/d \propto I^{-\frac{1}{2}}$ [7–9]. At relatively high inertial numbers ($I\gtrsim 0.1$), this length scale reaches a minimum of $\ell =d$ and the velocity fluctuations are given by $\delta v \propto d \dot \gamma$. Accordingly, the diffusivity can be expressed as

$$\begin{equation} D \propto \begin{cases} d^2 \dot \gamma, & \text{for } I \gtrsim 0.1,\\ d^2 \dot \gamma \frac{1}{\sqrt{I}}, & \text{otherwise.} \end{cases} \end{equation} \tag{ 4 }$$

The scaling (3) is a local constitutive law that relates the inertial number to the level of stresses within the flow. μ is the ratio of shear-to-normal stress, $\mu_s$ a yield criterion and b a numerical constant [7,10,11]. Like Bingham fluids, this law indicates that there should be no flow ($\dot \gamma = 0$) if the shear stress is lower than a threshold, $\tau< \mu_s \sigma$. According to (1) and (2), there should then be no diffusion.

However, most granular flows are not homogeneously sheared, owing to some gradient of stresses. Little is known about the shear-induced diffusion in these heterogeneous flows.

Non-local effects arise in heterogeneous flows that cannot be captured by the local constitutive law alone. For instance, a flowing layer can induce some flow in a nearby layer where the stresses are below the yield criterion ($\mu< \mu_0$) [12–14]. We refer to such layers as sub-yield layers. A number of non-local models have been introduced to capture the profiles of inertial number in heterogeneous granular flows, including in sub-yield layers [15–22]. However, established non-local models such as the non-local fluidity model [18] and the gradient expansion model [19] are expressed in the form of a second-order PDE. Their predicted shear-rate and inertial number profile in a specific geometry therefore depends on the choice made for two boundary conditions, for instance the fluidity at the wall and its gradient, which are not usually known a priori.

The lack of an explicit expression for the shear rate profile in heterogeneous flows impedes predicting the shear-induced diffusion. Furthermore, even if the shear-rate were known, there is no evidence confirming or challenging the validity of the scalings (1) and (2) in sub-yield layers.

In this letter, we explore the properties of shear-induced diffusion in heterogeneous granular flows, with a special focus on sub-yield layers. To this aim, we have simulated a series of steady and heterogeneous granular flows in which non-local effects arise. Velocity fluctuations and diffusivity are systematically measured in different parts of the flow, in order to assess the validity of the scalings (1) and (2). We then seek to identify an explicit expression capturing the profile of shear rate in these flows.

We use a Discrete Element Method to simulate dense granular flows in three geometries: plane shear (PS), plane shear with gravity (PSG), and Poiseuille flows (PF). These geometries are illustrated in fig. 1 and detailed below.

Zoom In Zoom Out Reset image size

Table 1:.  Parameters of the simulated flows: Plane Shear (PS) from [9], the Plane Shear with Gravity (PSG) and Poiseuille flow (PF), and corresponding symbols used in figs. 1 and 3 to 5. Filled symbols correspond to layers above the yield criterion ($\mu(y)>\mu_0$) and open symbols correspond to sub-yield layers ($\mu(y)<\mu_0$).

Symbol	Geometry	$H/d$	$10^3\frac{P_w}{E}$	$v_w \sqrt{\frac{\rho}{P_w}}$	$\frac{g}{d}\sqrt{\frac{\rho}{P_w}}$
$\rlap{+}{\times}$	PS	120	1	–	0
	PSG	60	0.4	0.316	0.0095
	PSG	60	0.4	0.791	0.0126
	PSG	30	0.4	0.316	0.019
	PSG	30	0.4	0.791	0.019
	PF	80	1	–	0.01
	PF	80	1	–	0.0125
	PF	40	1	–	0.036

All tests involve grains that are 2D disks of average diameter d, mass m and density ρ. They interact by elastic, frictional and dissipative contacts characterised by a Young's modulus E, a coefficient of friction $\mu_g=0.5$ and a coefficient of restitution $e=0.5$ for normal impact. A uniform polydispersity of $\pm 20\%$ is introduced on the grain diameter to prevent crystallisation during shear. There is no contact adhesion and no fluid in the pore space. The value of these contact parameters only marginally influence the flow properties, as discussed in [7,8].

Plane shear with no gravity presents the advantage of producing homogeneous stresses and shear rate throughout the system. Biperiodic boundary conditions are used to avoid walls and the shear heterogeneities they induce [22]. A series of steady flows were performed by prescribing a constant normal stress $P=10^{-3}E$ and different values of shear rate $\dot \gamma$. These simulations enabled us to measure the local constitutive law $\mu(I)$ of the materials by averaging the shear-to-normal stress ratio and the inertial number spatially in the entire flow, and temporally over 25 shear deformations. The measured values of μ vs. I are shown in fig. 1(b). They are best fitted with the local constitutive law (3) using $\mu_0 = 0.26$ and $b=1.1$, which is consistent with previously reported values [7]. It is worth noting that flow occurs only if $\mu>\mu_0$ in this geometry.

Unlike homogeneous plane shear flows, PSG and PF involve some stress gradient that leads to a non-homogeneous shear and a spatial variation of the inertial number in one direction.

PSG is simulated in a periodic domain along the x-axis, while parallel walls bound the system in the y-direction (see fig. 1(c)). Walls are made of aligned contacting grains with average diameter 2d, which effectively prevents wall slip [18]. Wall grains do not rotate and move as a rigid body. The top wall can translate in both directions to prescribe a shear velocity v_w. It is also subjected to a vertical external normal stress P_w. The vertical wall's motion is governed by an inertial dynamic. Its acceleration $\ddot y$ is given at any time during the flow as $M \ddot y = Ld (P_i-P_w)$ where L is the length of the wall, M the total mass of the wall grains, and P_i the internal normal vertical stress due to contacts between wall and flowing grains. In steady states, the wall vertical position is nearly constant, with some fluctuations smaller than d. The bottom wall is immobile. Flowing grains are subjected to gravity g and the corresponding body force $f_b = \pi d^3 \rho g/4$ in the direction transverse to the shear. This produces a gradient of normal stress in the y-direction, while the shear stress is constant. The stress ratio μ is thus maximum at the top and minimum at the bottom. It is possible to tune the external normal stress P_w and shear velocity v_w in such a way that the flow is comprised of a top layer where ($\mu(y)>\mu_0$) and a bottom layer where ($\mu(y)<\mu_0$). The local constitutive law predicts that there should be no flow in this layer, and therefore no diffusion. However, fig. 1(d) shows that the inertial number is in fact not null in this layer, owing to non-local effects. This suggests that there may be some diffusion in this layer.

PF is also simulated in a periodic domain in the flow direction and between two parallel walls. Unlike PSG, walls do not move in the flow direction and produce no shear. Both walls are subjected to an external normal stress P_w and are free to move in the y-direction according to an inertial dynamic similar to that used in PSG. A body force f_b is applied on flowing grains, but this time in the flow direction x. This leads to a gradient of shear stress in the y-direction, while the normal stress is constant. As a result, the stress ratio is maximum near the walls and minimum at the centre. The body force and the applied pressure can be tuned in such a way that a middle layer develops that is below the yield criterion ($\mu(y) <\mu_0$), while the two layers near the walls are above the yield criterion ($\mu(y)>\mu_0$). Like PSG, fig. 1(f) shows that non-locallity induces some flow in the central zone while it is below the yield criterion, suggesting possible diffusion.

To highlight the diffusive behaviour in these flows, we measured the typical grain trajectory characterised by their mean square displacement. We took advantage of the time invariance of the steady flow and of their spatial homogeneity (at least in one direction), to measure an averaged mean square displacement defined as

$$\begin{equation} \Delta^2 (t) = \frac{1}{MN}\sum_{i=1}^{N} \sum_{j=1}^{M}(y_i(t_j) - y_i(t_j+t))^2, \end{equation} \tag{ 5 }$$

where t is a time interval, t_j a reference time and y_i the y position of grain i at a given time. Averaging is performed considering a series of $M = 100$ reference times selected at random during steady flows. It is also performed on N grains. In homogeneous plane shear flows, all grains can be included in this average, leading to a single MSD for one given flow. In contrast, it is expected that the MSD might depend on the initial position of the grains y(t_j) in heterogeneous flows. MSDs $\Delta^2(t,y)$ are then measured at different position y by averaging on grains initially located within strips of width d centred at y. In principle, grains may experience different values of shear rate as they diffuse through different position y, and could therefore experience a variable diffusivity. In the following, we will limit the analysis of the MSD to relatively small values ($\Delta \lesssim d$), so that these variations are negligible.

Figure 2 shows examples of MSD evolutions at different layers within flows in the PS, PSG and PF geometries. It appears that these MSDs first exhibit a power law $\Delta^2(y,t) \propto t^{2}$ at short time scales. This marks a super-diffusive behaviour, reflecting a ballistic (or constant speed) grain trajectory, as observed in [9,23]. In contrast, MSDs exhibit a normal-diffusive behaviour $\Delta^2(y,t) \propto t$ at long time scales. A coefficient of self-diffusion D can be measured in this regime using the Einstein formula [24]:

$$\begin{equation} \Delta^2(y,t) =2D t \end{equation} \tag{ 6 }$$

Figure 2 indicates that this normal-diffusive behaviour develops after a period of time proportional to the shear time $\dot \gamma^{-1}$. Seemingly, normal diffusive behaviour arises after approximately a tenth of shear deformation ($t\dot\gamma \gtrsim0.1$) in all layers and in all flow geometries. Then, the value of mean square displacement is also similar in all cases, approximately $\Delta^2\approx 10^{-2} d^2$, which corresponds to a typical grain displacement of $10^{-1}d$.

Zoom In Zoom Out Reset image size

Most importantly, MSDs exhibit two qualitatively different behaviours in layers above and below the yield criterion. Above yield, the super-diffusive regime is directly followed by the normal-diffusive regime. In contrast, a sub-diffusive regime develops after the super diffusive phase in sub-yield layers. This sub-diffusive phase is characterised by a plateauing of the MSD, which denotes a caging of the grains [25,26].

This caging only develops in sub-yield layers. In particular, it does not develop in homogeneous plane shear, even at low inertial numbers. As a consequence, this caging appears to be a distinguishing feature of the grain trajectories in sub-yield layers.

As a way to assess the validity of the scalings (1) and (2) in heterogeneous flows, we have measured the profiles of diffusivities D(y) and velocity fluctuations $\delta v(y)$ within flows in the PSG and PF geometries. Figure 3 shows how these quantities scale with one another, as well as with the local inertial time t_i(y) and shear rate $\dot\gamma(y)$. These results highlight the following three observations.

Zoom In Zoom Out Reset image size

The first observation is that the scaling (1) of the diffusivity with the velocity fluctuation is not valid everywhere in heterogeneous flows. This is evidenced in fig. 3(a). This scaling is valid for layers with the highest velocity fluctuations where results indicate $Dt_i d^2 \approx 0.1 v t_i/d$, which is equivalent to $D \approx 0.1 \delta v d$. These layers correspond to the zone of the flow above the yield criterion. In contrast, there is a neat breakdown of this scaling in sub-yield layers.

The second observation is that the velocity fluctuation scaling (2) is not valid everywhere in heterogeneous flows. This is evidenced in fig. 3(b), which suggests two limits. At high inertial numbers, data seemingly converge toward the scaling $\delta v t_i/d \propto I$, or $\delta v \propto d \dot \gamma$, which is similar to that measured in homogeneous shear flows in this range of inertial numbers. In layers with low inertial numbers, which are sub-yield layers, results suggest that the velocity fluctuations become independent of the shear rate and controlled by the inertial time:

$$\begin{equation} \delta v \propto d/t_i. \end{equation} \tag{ 7 }$$

In this range of inertial number, this scaling differs from the one measured in homogenous shear flow. It indicates that velocity fluctuations do not vanish in sub-yield layers even when the shear rate tends to zero. They would vanish in homogeneous plane shear, according to (2). To interpret this shear-rate independence of $\delta v$, it is important to consider that it corresponds to instantaneous velocity fluctuations. We will see below that velocity fluctuations averaged on a period of time of $\dot \gamma^{-1}$ do scale like $\dot \gamma d$, even in sub-yield layers.

The third observation is that there is a simple scaling between the diffusivity and the local shear rate in all layers of all tested PSG and PF flows. This scaling, evidenced in fig. 3(c) is

$$\begin{equation} D \approx 0.1 d^2 \dot \gamma. \end{equation} \tag{ 8 }$$

It differs from the diffusivity scaling measured in homogenous plane shear flows in this range of inertial numbers, which is $D \propto d^2 \dot \gamma/\sqrt{I}$.

Two conclusions can be drawn from these scalings and from the MSD evolutions. The first conclusion is practical: one can directly deduce the profile of diffusivity in a heterogenous granular flow from the shear rate profile, according to (8).

The second conclusion concerns the physical process controlling the diffusivity. In homogeneous shear, instantaneous velocity fluctuations control grain self-diffusion. The underlying mechanism is a random walk with a step size proportional to d during which grains experience a constant velocity $\delta v \propto \dot \gamma d$, which is sustained for a period of time $\dot \gamma^{-1}$. In sub-yield layers, the diffusivity is seemingly controlled by a different process. Grains still undergo a random walk of step size proportional to d, and frequency $\dot \gamma$, as evidenced by the MSD at long time scales. However, the elementary step of this walk is comprised of fast inertial displacements of typical velocity $\delta v = d/t_i$ lasting a period of time proportional to t_i, leading to a small displacement ($\Delta \ll d$). These displacements correspond to the super diffusive behaviour observed in figs. 2(b) and (c). However, these short-lived displacements seemingly produce caged trajectories at longer time scales, as evidenced by the plateauing of the MSD. This caging lasts approximately $\dot \gamma^{-1} - t_i$. As a result, the average velocity fluctuation of grains during one shear deformation defined as $\Delta v = \frac{1}{\dot \gamma^{-1}} \int_0^{\dot \gamma^{-1}} \delta v (t) \mathrm{d}t$ can be expressed as: $\Delta v \propto \delta v t_i/\dot \gamma^{-1} \propto \dot \gamma d$. The diffusion coefficient can then be expressed as $D \propto \Delta v d$. Consequently, the intensity of the instantaneous velocity fluctuations $\delta v$ no longer influences the diffusivity in sub-yield layers.

We now seek to identify a formula that explicitly defines the profiles of shear rate $\dot \gamma$ that are driving the diffusivity, as per eq. (8). The aim is to express the shear rate profile in terms of local stresses and position in the flow to ultimately infer the diffusivity profiles from these parameters.

A possible approach to predict the shear rate profiles is to combine the local constitutive law (3) with a non-local model [17–19]. However, existing non-local models are expressed in the form of a PDE and their predictions rely on the choice of boundary conditions, which are not always known a priori.

We follow here an alternative approach that has been recently proposed for amorphous solids, referred to as the scaling description [27–29]. These materials satisfy a Hershel-Buckley local constitutive law: they yield above a shear stress threshold ($\tau>\tau_0$), and then start to flow with a shear rate $\dot\gamma \propto (\tau-\tau_0)^\beta$, where $\tau_0$ and β are material-dependent parameters. Interestingly, amorphous solids also exhibit non-local effects that are similar to those in granular flows: some flow may exist in a zone below the yield criterion if an adjacent zone is flowing. The scaling description of such non-local effects consists in establishing a scaling for the local shear rate in terms of the distance to the yield, as follows [27,29]:

$$\begin{equation} \dot\gamma (y) = |\tau-\tau_0|^\beta \mathcal{F}_{\pm}(|\tau-\tau_0|^\nu|y-y_0|), \end{equation} \tag{ 9 }$$

where y is the position of the layer, y₀ is the position of the layer in the flow where $\tau(y_0)=\tau_0$, and ν is some exponent. $\mathcal{F}_+$ and $\mathcal{F}_-$ are referred to as scaling functions, which correspond to layers above ($\tau>\tau_0$) and below ($\tau<\tau_0$) the yield criterion, respectively.

We seek to adapt this approach to granular flows by considering the relevant frictional yield criteria and non-dimensional shear rate I, as

$$\begin{equation} I (y) = |\mu-\mu_0|^\beta \mathcal{F}_{\pm}(|\mu-\mu_0|^\nu|y-y_0| /d). \end{equation} \tag{ 10 }$$

Figure 4 shows that this scaling does capture the measurements in our simulated PSG and PF granular flows using $\beta = 1$ and $\nu=1$: data collapse on two curves, one for the layers below the yield criterion and one for the layers above the yield criterion. Above the yield criterion, $\mathcal{F}_+(x)$ becomes weakly dependent on x and of the order of 1, indicating that non-local effects may be neglected in these layers. In sub-yield layers, data suggests a transition from a power −0.5 to a power −2 for the function $\mathcal{F}_-(x)$. We introduce the following interpolation to capture these two regimes and their transition:

$$\begin{eqnarray} \mathcal{F}_- (x) &= \frac{A }{x^2+B\sqrt{x}}, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} x &= |\mu-\mu_0| |y-y_0|/d, \end{eqnarray} \tag{ 12 }$$

where A and B are the two fitting parameters. Figure 4 shows that this function satisfactorily captures the measurements with $A=2\times10^{-3}$ and $B=0.03$. Interestingly, x may be seen as a distance to the yield that includes a stress-wise distance $|\mu-\mu_0 |$ and a Euclidian distance $|y-y_0|/d$. Far from the yield, for $x\gg1$, the function $\mathcal{F}_- (x)$ tends to $\mathcal{F}_- (x) = Ax^{-2}$. Considering (10), the profile of inertial number is then given by

$$\begin{equation} I (y) \approx \frac{A}{|\mu-\mu_0|} \frac{d^2} {|y-y_0|^2}. \end{equation} \tag{ 13 }$$

Figure 5 shows how the profile of inertial number and the profile of diffusivity can be captured using the scaling prediction for the inertial number (10) and the scaling (8), in all layers of the heterogeneous flows.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that the scaling (13) is consistent with the self-activated mechanism underlying the non-local model introduced in [17]. This mechanism considers that plastic events may be triggered in sub-yield layers by stress fluctuations originating from remote flowing layers. In this model, it is thought that stress fluctuations decay as the distance to the power −2 from their origin. This is consistent with the scaling $I (y) \propto \frac{d^2} {|y-y_0|^2}$. Further still, it is thought that the magnitude of these stress fluctuations required to uncage a grain is proportional to $|\mu-\mu_0|$, which is consistent with the scaling $I (y) \propto \frac{1}{|\mu-\mu_0|}$.

Nonetheless, the scaling (10) and the choice of the exponent $\alpha=1$ and $\beta=1$ are empirical. While they lead to a satisfactory collapse of the data, it is possible that other functional forms or values of the power exponents could lead to a collapse of the data onto different functions $\mathcal{F}_{\pm}$.

This study points out that shear-induced diffusion is qualitatively different in granular layers flowing below and above the yield criterion.

Above the yield criterion, diffusivity is proportional to the instantaneous velocity fluctuations, which themselves are driven by the shear rate and possibly by the inertial number, as per (1) and (4). In contrast, diffusivity in sub-yield layers is not controlled by instantaneous velocity fluctuations, which become shear-rate independent and controlled by the inertial time. Grain trajectories are then found to exhibit fast inertial motions leading to some caging that lasts one shear deformation. Nonetheless, above and below the yield criterion, the shear-induced diffusivity is found to be directly proportional to the local shear rate, as per (8).

Predicting the diffusivity profiles in heterogeneous flows then reduces to the question of predicting their shear rate profile. This shear rate profile may be deduced from existing non-local continuum models, or by the scaling approach we introduced, which led to (13). This scalings can readily be used to resolve diffusion processes in non-local granular flows.

However, the scaling description we introduce here is purely empirical, and its physical basis remains to be established. Consequently, it is possible that different functional forms or different exponent values may equally well capture the shear rate profile in sub-yield layers.