Evidence of a non-local φ(I) response

Granular dilatancy has been previously characterised through a simple linear relationship between the packing fraction and dimensionless shear rate. However, this relationship was developed for granular flows in a simple shear cell geometry. Here we examine inertial volume changes in a shear cell with gravity, a vertical chute, and a pseudo-2D hopper. In so doing, we show that the packing fraction displays both a local and non-local response, analogous to what is typically observed for the stress ratio μ.


Introduction
One challenge when describing granular flows is their inherent dilatancy. Granular flows will expand or contract in response to changes in shear or pressure. However, many continuum descriptions of granular flow rely on an assumption of a uniform, or practically uniform, packing fraction (φ) [1][2][3][4]. Perhaps because of this assumption, examination of changes in granular packing has been highly limited.
Perhaps the most well known approach to describing granular dilatancy is the linear φ(I) relationship [5] which relates the packing fraction to the inertial number (I) through a simple linear equation: where φ c and A represent fitted constants and I is defined as |γ| d ρ s /P where d is the local particle diameter, ρ s is the particle density and |γ| and P are the equivalent shear rate and pressure, respectively (both of which were defined as in [6]). This relationship is not the only φ(I) relationship and later works have suggested that φ ∝ I a or φ ∝ 1/(1 + AI) [7,8]. These are not discussed in depth as they represent only a slight improvement over the linear φ(I).
The linear φ(I) relationship defined in [5] was based on discrete element simulations of simple 2D shear cells. This result has been validated using similar 2D and 3D geometries [9,10]. Evaluation of changes in packing outside of simple shear systems has been relatively limited, only a few other systems with fairly narrow ranges of values for I have been considered [7,11]. The purpose of this work is to examine granular dilatancy using discrete element simulations (DEM) for various geometries over a broad I range, in order to establish whether the linear φ(I) relationship holds true for these systems. * e-mail: daniel.holland@canterbury.ac.nz

Methods
For this work four different geometries have been considered, shown in Figure 1. These are, a shear cell, a shear cell with gravity, a vertical chute, and a pseudo-2D Hopper. DEM simulations of these geometries were run using the open source software LIGGGHTS [12]. Simulations were run using a soft-sphere Hookean contact model, details of which can be found in the LIGGGHTS documentation [13]. The particle properties used for these simulations are given in Table 1. The timestep used was 1 µs. A size distribution of 0.8 − 1.2d m was also implemented to prevent crystallisation. The characteristic velocity, Young's modulus and Poisson's ratio (which together in LIGGGHTS define the spring constant) were chosen to ensure the simulations met the criteria for rigid particles, as defined in [14]. Additionally, a simulation (not shown) of the shear cell with gravity (see Figure 1b) was run with a Hertzian contact model and the same properties. For this simulation, the rigid criteria was defined as in [10], which is one of the few works that have replicated Equation 1 in 3D. This simulation produced functionally identical results to the equivalent Hookean simulation.
The two shear cell geometries and the vertical chute are periodic along the x and y axes with the top and bottom walls made of dense layers of particles with properties as in Table 1. All three geometries were run with a total of  (Figure 1a) was run with a series of different wall velocities (V w ) and an applied pressure (P w ) of 300 or 100 Pa. The shear cell with gravity ( Figure 1b) was run with P w =200 Pa, V w =4 m s −1 and G=9.81 m s −2 . The vertical chute (Figure 1c) had the wall particles fixed so that the chute's width along the z axis was constant. It was run with G set to 5 m s −2 . Finally, the hopper ( Figure 1d) had smooth frictional walls with properties identical to the particles and was periodic along the y and z axes. Particles leaving through the outlet would cycle back to the top. The hopper was filled with 53125 particles to a height of approximately 160 mm. The simulated geometry had space to allow for particles to free-fall before landing on the bed. It was run with G=9.81 m s −2 .
To extract data from the simulations coarse graining was used. This approach has been covered in depth by previous works [15][16][17] and so will not be detailed here. Notable to our particular approach is the usage of a Lucy function as the weighting function: where r is a position vector and c is the cutoff, here set to 3d.
When evaluating the shear cell, shear cell with gravity or chute geometries, values were calculated at different positions along the z axis (a total of 30 z positions were considered). At each z position, properties were calculated at 16 points spaced evenly across the xy plane and averaged to obtain the layer value. For the hopper geometry, the origin was placed at the outlet center. The solid fraction was evaluated at 30 positions on the z axis (with x=0) starting at 1.2 mm above base of the hopper up to a height of 15 mm. At each position, 4 points, evenly spaced along the y axis, were averaged to get the value at that z position. For all systems, the values at each position were subsequently averaged over 200 time steps, which were output at 0.1 s (100,000 DEM time steps) intervals. When coarse graining, the edge of the coarse graining function's cutoff was kept 1.5 times the maximum particle radius from the walls.

Results
In order to first validate our results, Figure 2 shows the φ and I values extracted from a series of shear cell simulations. Least squares regression was used to fit Equation 1 to this data using the 'polyfit' function in Matlab [18]. The fitting gave φ c =0.592 and A=0.127. Figure 2 shows that, over the evaluated I range, the packing fraction from each simulation collapses onto a straight line. There is some slight curvature, and the fit could be improved through the use of one of the alternative φ(I) relationships, such as a power law. However, the difference is fairly minimal. Thus, in a simple shear cell the linear relationship between φ and I is recovered, as was seen in prior work [10]. At low I values the packing fraction approaches 0.6, far from the random close packing limit of 0.64 [19], suggesting that even slow moving shear cells are still far from the maximum possible density. Now that it has been established that we can successfully recover the linear φ(I) relationship, next we look at   Figure 3b, which shows the behaviour on a semi-log scale. As I decreases the limiting φ value approached, and the nature of the approach, varies across the three geometries. Thus, the φ(I) response is geometry dependent.
Given the response that has been observed in Figure 3, an analogy can be drawn to the non-local rheology of gran- The local µ(I) rheology states that the stress ratio µ (defined as |σ| /P where |σ| is the equivalent shear stress, defined as in [6]) is a function of the inertial number [4]. Notably, this relationship holds only if µ > µ s where µ s is a yield value. Flow in regions where µ < µ s breaks from the µ(I) rheology and can be considered non-local [20]. Figure 4 shows µ and I extracted from our simulations. For the shear cell simulations, these were fit to a simple form of the µ(I) rheology [4], using Matlab's 'lsqnonlin' function [18]. This gave µ s =0.374 and the two fitting parameters, I 0 and µ 2 , equal to 0.767 and 1.075 respectively. Figure 4 shows that the shear cell data displays solely a local response, smoothly approaching a constant value as I approaches 0 (hence why Equation 3 was fit to this data). However, the shear cell with gravity, vertical chute and the hopper all display both a local and non-local response. The local response is seen when µ > µ s and all of the values collapse onto the plotted curve. The non-local response is seen for low I values where µ drops steeply. Figure 4b shows that, when plotted on a semi-log scale the non-local response differs between the geometries, preventing it from being easily characterised using a fitted model. The non-local response shown in Figure 4 mirrors what is observed in Figure 3. Indeed, in Figure 3 the values where µ > µ s have been plotted with filled in markers showing that the non-local region of the flow corresponds to the region where where φ deviates from the linear φ(I) rheology. Thus it can be inferred that the response of µ and φ are analogous and the deviation in φ represents a nonlocal response in the packing fractions. This analogy is not perfect as the pseudo-2D hopper shows a notable spike in φ while µ > µ s . This geometry is also transient in the Lagrangian sense, with a non-uniform flow profile along the z axis, and it's possible transient affects could lead to the observed discrepancy. As such rather than characterising the φ response based on µ s it might be more reasonable to use φ c . Thus, the packing fraction is divided into a local response, wherein φ < φ c and a potentially non-local response when φ > φ c . It is this possibly non-local response which displays the dependence on geometry, similar to µ, which makes describing it difficult. Recent work has shown that, the non-local µ(I) response may be potentially collapsed through incorporating an additional dependence on granular temperature [21]. A similar dependence may be necessary to capture the non-local φ(I) response.
When dealing with dense granular flows, the extent to which treating φ as constant affects the accuracy of any continuum model is unclear. However, given the inherently interconnected nature of flow properties, it is difficult to claim that there is no effect. Thus, there is benefit to being able to accurately describe granular packing. Our results show that trying to describe granular dilatancy using a simple linear φ(I) model can prove effective, but only when dealing with flows that display purely local behaviour. For more complex flows, we have shown that there is a non-local response that is not captured by the simple description of granular dilatancy. This necessitates the development of more complex equations of state in order to properly capture the behaviour of dense granular systems.