Strain dependent vorticity in sheared granular media

Displacement fields in sheared particle packings often display vortex-like structures that reveal essential features about the mechanical state of the collection of particles. There are several metrics to quantify these flow field features, yet extracting such quantitative metrics from flow field or particle tracking data involves making numerous choices on the time and length scales over which to average. Here we employ a much used experimental data set on sheared disk packings to explore how such arbitrary data mining choices affect the obtained results. We focus on calculating the strain dependent vorticity, as this metric is a differential method hence potentially sensitive to the way it is computed. We find that the total surface area with an absolute vorticity above a certain threshold approaches a plateau value as shear progresses. This plateau value exhibits a non-monotonic dependence on packing fraction. We also show which range of choices yields results that can support an analysis method independent, physical interpretation of the flow field data.


Introduction
Amorphous materials have long intrigued, among others due to the complex flow patterns these systems generate under shear. Substantial efforts have been dedicated to elucidate the origin and physical meaning of such nonaffine displacement fields. These fields have been studied via metrics, often non-local, like vortices [1][2][3][4][5][6], a two point velocity correlation function [7] and χ 4 , a four point displacement correlation function [8,9], the Okubo-Weiss parameter borrowed from turbulence [10], and a directional particle transport mechanism recently reported by the authors [11]. The benefit of probing the distribution of granular velocity and vorticity, is that they can be observed directly, as the vortex-like structures are embedded in the non-affine displacement fields of the granular materials. Unfortunately, extracting such quantitative metrics from flow field or particle tracking data involves making numerous choices on the time and length scales over which to average. It is not always obvious how such arbitrary choices affect the resulting data, which may impact the associated physical interpretation of the observations.
Here we employ a much used experimental data set on sheared granular media taken at Duke physics department to explore how such data analysis choices affect obtained results. We focus on calculating the strain dependent vorticity, as this metric is a differential method and hence most likely to be sensitive to the way it is computed. We show which range of choices yields results that can support physical interpretations of the data that are independent of the analysis approach. We use the curl of the displace- * e-mail: zhenghu@tongji.edu.cn ment field as the definition of the vorticity as this is the most common method to extract rotational displacement. We coarse grain the displacement field onto a regular grid before calculating the vorticity and evaluate the total vorticity above a given threshold for a set of different strain intervals. This allows us to verify the independence of this metric with respect to choices of the strain step used.

Methods
The data originates from numerous published works on a well studied setup that can be described briefly: we apply simple shear to bidisperse granular systems composed of photoelastic discs (Vishay PSM-4, diameters 0.64 cm and 0.80 cm) with inter-particle friction coefficient µ ≈ 0.7. Shear is performed in an apparatus that suppresses shear bands and related shear-induced density fluctuations entirely [11][12][13]. Particles were carried by the separate slats that form the shear box base. These slats move affinely in accordance with the applied shear (see Fig. 1). The optical properties of photoelastic discs reveals particle-scale contact forces when placed between a pair of crossed polarizers [14,15]. The system contained approximately 45 × 20 discs with a large to small number ratio 1 : 3.3 to prevent crystallization. Every run at given packing fraction φ, ranging from 0.69 to 0.816, was repeated five times; the initial stress-free state was prepared anew for each run. The precision on the packing fraction is about 0.001, as we can count the number of particles in the box, the accuracy of the packing fraction is affected by system and particle size measurement errors and boundary effects and is likely larger than 0.001.  Figure 1. Schematic of special apparatus for applying simple shear strain to a collection of photoelastic discs at the beginning (left) and the end (right) of shear. In the sketch, both the slats and particles are drawn much larger relative to the boundaries than in the real experiment. The x−y axes indicate the coordinate system in the lab frame, where simple shear is applied along the y axis.
Shear was applied quasi-statically in the y direction to a shear box ( Fig. 1). Starting from a stress-free state with a parallelogram box, the system was sheared by strain step of δγ = 0.0027. The total shear strain γ is defined as the total deformation of the box in the y direction divided by the box length l (see Fig. 1). Then the system was left to relax for six seconds, followed by taking two images, which reveal information on particle position and photoelastic response. Such a process -stepwise shearing, relaxing and imaging -was repeated until a total strain of 0.54 was achieved. From the images, we tracked particle positions and calculated particle displacements. Henceforth, we use the average diameter of particles in the systemD as the unit length in all results.
After tracking each particle for the whole experiment, we took out the cumulative affine displacements at each strain step by fitting particle displacements in y direction (∆y) as a linear function of particle x positions to obtain non-affine displacements, as shown by arrows in Fig. 2. We then determined non-affine particle velocities at a given strain by dividing their non-affine displacements from ∆ 2 strain steps backward to ∆ 2 strain steps forward by a chosen value ∆ (we set one shear strain step to be the unit time). Here ∆ = 2, 4, 6, 8, and 10 are studied. To obtain vorticity V, we then mesh the shear box with a mesh grid size to be approximately the small particle radius. The velocity field v at these mesh grids is smoothed by cubic interpolation. Finally the vorticity V simply follows the definition V(γ) = ∇ × v. The curl is then computed numerically on the grid; we thus measure vorticity per grid site, and the grid size is regular and interpolated based on irregularly spaced particle data.

Vorticity
To show examples of the non-affine displacement fields and associated vorticities, we plot spatially coarse-grained version of these quantities in Fig. 2 for a given packing with φ = 0.796 at γ = 0.0972 for step sizes ∆ = 2, 6, and 10. In general, we see that both positive and negative V exist in the non-affine displacement field in granular systems under linear shear, in accordance with literature [4]. Fig. 2 also shows that continuous regions with the same spin span over approximately 5 − 10 grain sizes, much bigger than the grid size we took to interpolate the nonaffine displacement field. Although we did not check for the dependency on the grid size, this result suggests that the latter does not play a role in the determination of the rotation. Though irregularly distributed spatially, regions with the same spin seem to occupy roughly half of the total area in the system. We can see that different values of ∆ have little affect on the obtained geometry of the vorticity field, e.g., values of V at any given location and the general size of regions with the same sign V. We see this directly from qualitatively comparing vorticity as determined by various ∆ shown in Fig. 2.

Spatial characterization of vorticity
We further quantitatively investigate the effect of choosing different ∆ on computing V by measuring the total surface area of regions with absolute value of vorticity, |V|, greater than |V| th . We choose this threshold to be the average |V| over all the grid sites from all the strain steps for the lowest packing fraction. We denote the thresholded surface area A V . Fig. 3(a) shows the evolution of A V as a function of γ at φ = 0.794 for various ∆. For any given ∆, A V first increases with γ and then reaches a plateau and fluctuates around that plateau. Strikingly, the onset γ at which A V reaches the plateau, γ onset , seems not to depend on ∆. For example, in Fig. 3a, γ onset is about 0.1. However, we can see a clear dependence on ∆ for the plateau value of A V . As shown in Fig. 3(a), we see an increase in the plateau value of A V as we increase the value of ∆. In addition, as ∆ gets larger, the increase of the plateau value seems to saturate.
We then look at the impact of the packing fraction on vorticity. Fig. 3(b) shows examples of A V as a function of γ for various φ with ∆ = 6. Generally, we see that γ onset for A V decreases as φ increases. In addition, we observe that the plateau value of A V varies with φ. To better characterize this A V (φ), we determine the plateau A V (φ) value, denoted asĀ(φ), by averaging A V for all γ greater than γ onset . We showĀ(φ) in Fig. 3(c). We see that for any given ∆,Ā(φ) shows a non-monotonic behavior as φ increases. For both lowest and highest φ,Ā(φ) is relatively low and peaks when φ ≈ 0.785. We attribute a lowĀ value for both low and high φ to different causes: At low φ, the average spacing between grains is relatively large and less force-bearing contacts are thus expected to occur during the process of shear. Therefore, there will be little nonaffine displacements, resulting in low vorticity in general, hence less regions with high |V|. At high φ, a lowĀ is more likely caused by little space allowed for large-scale particle rearrangement to occur. Similar to the observation in Fig. 3(a), Fig. 3(c) also shows thatĀ(φ) seems to increase but asymptotes with ∆. The physical interpretation of this result is that the observation of the existence of a particular volume fraction at which the vorticity is maximum is not a feature of the analysis settings, but an intrinsic property of the system. Interestingly, the φ for which the peak in A is observed is the same for which we observed a substantial decrease in the diffusive particle transport prop-erties [11]. Apparently the confinement effects above the peak φ suppress both non-affine displacement fields, suggesting a connection between the underlying mechanism that drives both vorticity and diffusion [6]. However, below the peak φ, we do observe thatĀ increases with φ, yet D 1 [11] is constant in this regime, suggesting a disconnect between vortices and non-affine diffusion effects.

Probability distribution function of vorticity
To better understand the behavior of A V andĀ and investigate more locally the effect of ∆ on V, we further measure the mean of absolute value of V, |V| , at different strain ranges, and we study the probability distribution function (PDF) of |V|. As mentioned in Sec. 3.2, we choose the low- est packing fraction (φ = 0.692) to obtain the threshold for determination of A V . Hence, we will focus on φ = 0.692 in this section. One would expect smaller |V| values for the beginning of shear as compared to later stages of shear as the onset of shear has not yet resulted in enough particle interactions and resultant non-affine displacements. This strain dependence may affect the calculation of |V|, so the question we ask here is: will the calculation of |V| threshold change when we change the range of strain used in the calculation of |V| th ? To answer this, we measure |V| for γ < 0.1 and γ > 0.1, as shown in Fig. 4(a). Indeed we see that |V| is smaller for γ < 0.1 than those calculated from γ > 0.1 and from all strain steps. In addition, we observe that the trend of |V| decreasing with ∆ increasing is not affected by the choice of strain range chosen for calculating |V| . Again we have observed that the differences in |V| for ∆ ≥ 6 is much smaller for ∆ < 6. Therefore, in addition to little effect on V from ∆ qualitatively as shown in Sec. 3.1 and Sec. 3.2, we expect the choice of ∆ will not change qualitative results for V when ∆ is relatively large (∆ ≥ 6). This is also observed in Fig. 3(c).
With the effect of ∆ on a global measurement |V| being the same regardless of γ range we test, we then turn our attention to the question on how different choices of ∆ change the local behavior of V? To answer this question, we calculate the PDF of |V|, P(|V|), from all the strain steps at φ = 0.692. Since values of V roughly decrease as ∆ increases, as evidenced in Fig. 4(a), we show that the PDF of |V| after rescaling P(|V|) and |V| by |V| in Fig. 4(b). We see that the PDF decays roughly exponentially with |V| for all ∆ tested. This general dependence of the PDF on |V| does not change for all ∆ tested here, hinting that there is no qualitative difference for various choices of ∆ when calculating vorticity. In addition, we see the tail of the PDF gets higher as ∆ increases. This indicates that there are more regions with |V| greater than V , which explains high A V andĀ for larger ∆, as shown in Sec. 3.2. Note that the seemingly fat tail of the PDF, especially for large ∆, is not caused by large contributions from |V| at larger strain steps since we see no qualitative difference in the PDF of |V| sampled from different γ range (results not shown).

Conclusions
In summary, we have studied the vorticity V observed in quasi-statically sheared disk packings. We varied the time step size ∆ used to determine vorticity to verify that extracting the vorticity is independent of this numerical setting. We observe that multiple vortices develop in the granular material under shear. Vortex size generally spans 5 to 10 particles. We see a clear dependence of vorticity V on shear strain γ. This is evidenced by the evolution of total surface area (A V ) of |V| above a chosen threshold |V| th as a function of γ, with a plateau value for A V at large γ.
This plateau value,Ā, shows a non-monotonic dependence on packing fraction φ and peaks at φ ≈ 0.785. In addition, we find that qualitative features of the vortices, such as their PDF and A V (φ) are not changed, although the quantitative values of A V are affected by the choice of ∆. For future work, it is of great interest to investigate underlying mechanisms that govern the evolution of vortices and the interaction between adjacent vortices and other mesoscopic features in sheared granular media, such as force chains or grain loops. Such understandings will help advance continuum descriptions of particulate systems.