Contact evolution in granular materials with inherently anisotropic fabric

This paper presents the results of two triaxial compression tests performed on approximately 9×103 oblate spheroids (lentils) with different initial orientations with respect to the loading axis (approximately 30° and 60°). Typical stress-strain information is complemented with x-ray tomography scans every 1 % strain. Starting from an initial labelling of particles, a new technique is used to obtain unprecedented detail of tracking of all the particles through time. This rich dataset is analysed from the perspective of inter-particle contacts (building on previous metrological work on this subject) revealing that the mean contact orientation in both samples rotates towards the direction of σ1 at a rate depending on the initial orientation. Different trends are observed for the alignment of contact orientation, with a significant evolution observed in the sample prepared at 30° which is not as pronounced as in the other sample.


Introduction
Layering in sands deposits is frequently observed in nature, resulting most often from deposition under gravity. This layering can be modelled as anisotropy in granular fabric (i.e., packing density, orientation of particles, orientation of contacts). It is important to distinguish this inherent (initial) anisotropy from induced anisotropy due to loading. Previous studies on 2D DEM simulations [1][2][3] show that inherent fabric anisotropy affects the material response. The results are in agreement with the predictions of the Anisotropic Critical State Theory (ACST) [4]: during shearing, the fabric increases its anisotropy and tends to align with the direction of the principal stress. However, the intrinsic mechanisms under which this alignment takes places have barely been explored experimentally.
The present paper focuses on the effect of inherent anisotropy on the evolution of the contact fabric, making use of x-ray tomography to image two specimens comprising approximately 9000 lentils during triaxial compression. These specimens are built with particles oriented at a specific bedding angle with respect to the loading axis; and the effect of inherent anisotropy is observed on the evolution of the contact network. Special attention is given to the mechanisms that allows the fabric alignment with the direction of the major principal stress.

Experiments
Two specimens of green lentils are prepared by depositing approximately 9000 lentils into membrane stretchers with diameter of 60 mm and height of 145 mm, tilted at 30°( LENGP30) and 60°(LENGP60) from the vertical axis. * e-mail: gustavo.pinzon@3sr-grenoble.fr θ X Y Z α Figure 1: Inclination (θ) and azimuth (α) of a 3D vector These specimens are then subjected to triaxial compression (allowing no rotation of the end platens, and imposing only vertical displacement of the lower platen) under a confining pressure of 50 kPa. Both specimens have an initial void ratio of ≈ 0.54. Lentils are close to being oblate ellipsoids (2.3 mm in diameter and 1.3 mm in height giving an aspect ratio of 1.77), and are characterised by a plane of symmetry. The vector normal to this plane is taken as the orientation vector of the particles, and is decomposed into an orthogonal coordinate system described by its inclination (θ) and azimuth (α) (see Fig. 1). The initial average orientation of the particles of the specimens, within a confidence interval of 99%, are θ 0 = 33.9 • (±0.8 • ) and α 0 = 169.3 • (±0.8 • ) for LENGP30, and θ 0 = 58.4 • (±0.8 • ) and α 0 = 94.2 • (±0.8 • ) for LENGP60. As observed, the sample preparation method allows the generation of specimens with a precise targeted initial orientation, and with similar variability between tests.
Each triaxial compression test is performed in-situ, using x-ray tomography to image the specimens at 1 % axial shortening intervals. Figure 2 (top) shows the macroscopic stress-strain response for both specimens (with q = σ 1 −σ 3 and p = 1 3 (σ 1 + 2σ 3 )), where the relaxation parts correspond to the stages at which loading is temporary halted and the x-ray scan is performed. The bottom graph shows the volumetric strain evolution with shortening, highlighting an initial contracting phase followed by the dilation of the specimens.

Particle tracking and contact measurement
Each 3D image resulting from a scan has a pixel size of 180 µm and dimensions 600 × 600 × 850 voxels. All the subsequent image-based analysis is performed using the open-source software spam [5]. With the objective of characterising contacts, the particles first need to be identified. In the first scan the solid phase is identified with a global normalised threshold (based on the peaks of the attenuation field histogram), and the individual lentils are labelled using a watershed algorithm, with manual postprocessing to ensure 100 % fidelity. Given that the scans are performed at small increments of axial shortening, the position of each particle is followed accurately, using a novel discrete digital image correlation technique based on [6,7], allowing a consistent labelling of the particles in all subsequent scans. The tracking in this case is total (i.e., always keeping the 0%-shortening image as the reference configuration), and subsequent labelled images are generated by applying the unique measured transformation to each particle in the initial labelled image, and used as markers for the re-segmentation of the deformed binary image. For illustrative purposes, the coloured insets in Fig.  2 show grains coloured by their total 3D rotation angle, revealing two friction cones with little deformation.
Establishing whether two particles are in contact presents unique challenges, which arise from the physical interpretation of a contact, since two particles can be infinitesimally close but not touching each other. Following the technique presented by Wiebicke et al. [8], a contact-detection algorithm is applied, making use of a global threshold to identify possible contacts and a local threshold to refine the detection process (since particles are mostly smooth, a local threshold of 0.7 is used). Additionally, a random walker algorithm is used to determine the plane of the contact zone, and the contact orientation is taken as the normal vector to the contacting plane [9].
Finally, a von Mises-Fisher (vMF) probability distribution [10,11] is used to describe the spatial orientation and concentration of the contacts. This probability distri- bution is analogous to a multivariate Gaussian distribution and is specially suited to normalized vectorial data. In the vMF distribution the vectors present a rotational symmetry around a mean direction (θ i , α i ), and their tightness around it is parametrised by the concentration parameter κ i -larger values of κ i imply a smaller spread of the data and a higher concentration towards the mean direction (θ i , α i ).

Contact evolution
As suggested in [12], the evolution of contacts can be tracked by studying: a) the preservation of contacts between the two states, b) the creation of new contacts due to particles moving closer to each other, and c) the loss of contacts while particles move away from each other. Since the particles preserve their labels throughout the test, each contact subset (i.e., preserved contacts C = , created contacts C + , and lost contacts C − ) is computed for each pair of subsequent scans, making use of all the contact from the whole contact network. Additionally, a von Mises-Fisher distribution is used to describe each contact subset using two parameters: i) the mean inclination (θ i , following the conventions of Figure 1), and ii) the spread parameter (κ i ). Analysis of the azimuth angle α i of each distribution shows very little changes. This (non trivial) finding is not discussed further here for the sake of brevity. Figure 3 (top) shows the evolution of the mean orientation θ i of each subset during the test, as well as the mean orientation of the particles for the two experiments -the error bars correspond to a confidence interval of 99% for all the measurements. The first remark is that the tilted sample preparation method indeed ensures that the mean orientation of particles is as imposed, and that the mean orientation of C = (i.e., most contacts) is very close to the orientation of the particles. As both specimens are deformed, the mean orientation of C = smoothly decreases, approximately halving from beginning to end of test (from , showing a clear tendency of alignment toward the direction of the major principal stress. The mean orientation of the particles follows the same trend, but with a slight delay. A pronounced decrease of the mean orientation of the distribution can also be observed for C − (starting from However, the mean orientation of C + varies much less, since only a small decrease is visible for LENGP60. Figure 3 (bottom) presents the concentration parameter κ i of each subset, which measures the anisotropy of the orientations of the subsets -the error bars correspond to a confidence interval of 99% for all the measurements. Considering again the preparation of the samples, the ini-tial states of both specimens (i.e., essentially κ = i ) have a remarkably similar initial spread of orientations.
Unlike the mean orientation, there is quite a dramatic difference in the evolution of the spread of contact orientations between the two experiments. In LENGP30 all three subsets become more concentrated, with the preserved being the most concentrated and the lost being the least concentrated. The difference (κ + i − κ − i ) appears to stay constant during the test. In LENGP60, concentrations remain roughly the same -with κ = i being higher than κ + i and κ − i (which are essentially equal throughout) -until about 22 % shortening, after which they all start to increase.

Discussion
Under deviatoric loading, it is expected that prolate spheroids particles will align with the direction of the major principal stress, and this is clearly confirmed by the presented results. Like particles, contacts align to σ 1 , and do so faster than particles; in both cases contacts on average are created closer to 0°and among the contacts that are lost the mean orientation is always higher than the contacts that are preserved. Remarkably, the initial orientation of the particles has a direct effect on the mean orientation of the contact loss process along the test: while the difference (θ = i -θ − i ) decreases in LENGP30, it increases in LENGP60 as the specimen shortens.
Furthermore, in the case closer to coaxiality with the major principal stress (LENGP30), the distribution of contact orientations tightens during compression as particles align. This is not observed in the other sample, perhaps due to the disorganisation created by the large rotations that the particles undergo -it is remarkable that rotations greater than 30°are possible without increasing the scatter in the orientation of contacts.

Conclusions
This paper explores how the contact orientation distribution evolves for inherently anisotropic granular materials. For this, two specimens of lentils are subjected to in-situ triaxial compression test, using x-ray tomography to image the specimens at several strain levels. Each specimen (LENGP30 and LENGP60) presents a unique initial orientation (θ 0 ≈ 30 • , θ 0 ≈ 60 • , respectively). A novel discrete digital image correlation technique facilitates the tracking of the particles throughout the test, enabling a consistency in labelling of the particles among the scans. For each pair of subsequent scans, the contact network is compared and divided into three subsets: preserved contacts, created contacts, and lost contacts.
The results sheds light on the mechanisms under which the contact orientation distribution evolves towards an alignment with the direction of the major principal stress. For each strain increment, new contacts are created ahead while misaligned contacts are erased behind. The initial orientation of the particles affects the concentration of the distribution, as well as the extent of the contacts creation/loss process during the test.
The presented experiments at 30 and 60°initial orientations have already been supplemented by experiments at 0°, 45°and 90°, whose analysis is ongoing, and which will complete the observations presented. The new particle tracking method adopted in this work is capable of tracking particles throughout the test with no particle loss. This is unprecedented (to the best of the authors' knowledge), and the obtained tracked dataset opens up a broad range of perspectives for quantification. Further paths for analysis currently planned include the correlation of contact half-life as a function of orientation relative to the direction of the major principal stress; correlating particle kinematics as a function of orientation misalignment; the coupling between local fabric changes and local strains, and strain localisation mechanisms as a function of initial orientation.