Micro-scale anisotropy of contacts and pores in granular media

This paper considers the anisotropy associated with the interparticle contact network and interconnected pore space, which is quantified through a scalar quantity based on the fabric tensor. The study focused on the anisotropy associated with subsets of contacts and pores. It is shown that the stress ratio can be directly related to the geometric anisotropy of strong and non-sliding contacts, and that the anisotropy of pores is correlated with the anisotropy of contacts.


Introduction
Understanding the micro-scale origins of the complex macroscopic behaviour of granular media is essential to develop better predictive models.This requires the characterisation of micro-scale features including the network of interparticle contacts and the interconnected pore space.Numerous studies have suggested that anisotropy of micro-scale features plays a crucial role in the behaviour of granular materials [1][2][3].While most prior studies have focused on the contact network, this paper explores the evolving anisotropy of both contacts and pores.Interparticle contacts contribute to the stress transmission process, while pores influence deformation characteristics and hence, both play distinct roles in the behaviour of granular materials.Micro-scale features are explored through numerical simulations of two-way cyclic axisymmetric tests on monodisperse assemblies of spherical particles.A novel feature of this study is that anisotropy is quantified for subsets of contacts and pores, thereby providing an indication of which micro-scale features have dominant contribution on the macroscopic behaviour.

Numerical Simulations
Discrete element simulations were conducted using the open-source software LIGGGHTS [4].Details of these simulations are presented in [5].Axisymmetric simulations mimicking a drained triaxial test were carried out for loose (CP-L) and dense (CP-D) samples.Samples were generated by randomly inserting particles within a cubical domain (with no contact between particles), and then compressing the assembly by applying isotropic pressure on the rigid frictionless walls of the domain.The current simulations were conducted using a monodisperse assembly e-mail: a.sufian@unsw.edu.auThe interparticle friction, μ, is varied during isotropic compression to achieve samples with different pre-shear densities, but a constant value, μ = 0.25, was used for all shear simulations.The numerical analyses simulate quasistatic cyclic shearing (with an inertia number of I ≈ 10 −3 ) to a target axial strain of |ε a | = 20% through loading, unloading and reloading phases (Fig. 1).The macroscopic shear stress (Fig. 1a) and volumetric strain (void ratio, e, Fig. 1b) show behaviour that is qualitatively consistent with the measured behaviour of dense and loose sands.A critical state condition is approached in the loading phase, where the loose and dense sample exhibit similar stress state and void ratio at large axial strain.It should be noted that the void ratio was computed over a subset of the sample to avoid boundary irregularities associated with the rigid walls.

Anisotropy of Interparticle Contacts
The anisotropy of interparticle contacts can be quantified through a fabric tensor [6,7].This study is particularly concerned with the anisotropy associated with subsets of the interparticle contact network as forces are predominantly transmitted through a sparse subset of contacts.
In this study, a quad-partition is considered comprising four subnetworks: (i) strong and sliding contacts; (ii) strong and non-sliding contacts; (iii) weak and sliding contacts; and (iv) weak and non-sliding contacts, where the subnetworks are denoted {sμ, sη, wμ, wη} respectively.Strong contacts are defined as those contacts with an interparticle normal force ( f n ) greater than the mean, f n ≥ f n , while weak contacts are those with f n < f n .Sliding con-tacts are at the Coulomb friction limit, that is, The anisotropy of these subnetworks is quantified by considering the directional distribution of contact normals.This follows the stress-force-fabric relationship proposed in [5] for partitioned subnetworks, where further details can be found.If n is a unit vector along the contact normal direction, then the probability density function of contact normals with orientation n in the k th subnetwork (where k can have the values {sμ, sη, wμ, wη}) can be expressed by the second-order approximation: where the subscripts {i, j} are indices of the vector and tensorial quantities.D k i j is termed the deviatoric directional probability density tensor in the k th subnetwork which is given by: where is termed the moment tensor (or fabric tensor) for interparticle contacts.δ i j is the Kronecker delta and N k c is the number of contacts in the k th subnetwork.A scalar measure of anisotropy can be calculated by considering the secondinvariant of the directional probability density tensor: a k c provides a measure of geometric anisotropy, as it is purely concerned with the arrangement of contacts rather than the forces transmitted through the contacts.Figures 2a-2d illustrate the evolution of a k c for each of the four subnetworks in the CP-D simulation (similar trends occur for the CP-L sample).Note that a signed measure of anisotropy is considered, where positive values imply preferential orientation in the axial direction, while negative values indicate preferential orientation in the radial direction.
Figures 2a and 2b show that contacts in the {sμ, sη} subnetworks are preferentially orientated in the direction of loading, while contacts in the {wμ, wη} subnetworks are orientated perpendicular to the direction of loading (Fig. 2c-2d).This reflects the force chain description of the contact network [8], where contacts with relatively large force magnitudes form columnar structures, which are braced laterally by contact with relatively small force magnitudes.
The subnetwork of strong and non-sliding contacts, {sη}, is clearly the dominant contributor to the geometric anisotropy of the contact network (Fig. 2b).What is particularly interesting about the behaviour of this subnetwork is that it virtually mimics the stress-strain response in Figure 1a.In fact, the relationship between anisotropy of the strong and non-sliding contacts and the stress ratio can  be expressed as q p = 2 5 a sη c .This is demonstrated in Figure 3, which indicates that stress ratio can be directly related to the geometric configuration of strong and non-sliding contacts throughout the two-way cyclic test.The potential implications of this relationship is discussed further in [5].

Anisotropy of Pores
The anisotropy of pores can be quantified in a similar manner to Section 3 but requires a clear and unambiguous identification of individual pore bodies.In this study, individual pores are identified using the modified Delaunay tessellation detailed in [9].Initially, a classical Delaunay tessellation is performed to identify tetrahedral Delaunay cells, where the vertices of the tetrahedron are the particle centres.These Delaunay cells are merged together according to a connectivity criterion based on the inscribed sphere associated with each Delaunay cell.If the inscribed spheres of two adjacent Delaunay cells overlap, then the Delaunay cells are merged together (Fig. 4).This results in pores that are described as polyhedral cells, providing a physically intuitive delineation of the pore space that bet- ter captures local density fluctuations.This approach has been shown to describe water retention characteristics in monodisperse assemblies [10].
In order to quantify anisotropy of pores, the orientation of individual pores must be calculated.This is achieved using the pore orientation tensor [9]: where t α and a α are the unit normal vector and surface area of the triangular faces forming the polyhedral pore (Fig. 4) and N s is the number of faces.The major principal orientation can be determined by an eigen-decomposition.If p is the unit vector along the major principal direction, then the probability density function of pores with an orientation p is given by: This is analogous to the expression in Equation 1, and in a similar manner to Section 3: where is the moment tensor for the pore space.N k p is the number of pores in the k th subset of pores.As per the contact subnetworks considered in Section 3, the anisotropy associated with distinct subsets of pores will be determined.Two dual-partitions are considered: (i) merged and unmerged pores, denoted {m, u} respectively, and (ii) contracting and dilating pores, denoted {c, d}.Merged and unmerged pores are delineated through the modified Delaunay tessellation mentioned above.Unmerged pores have tetrahedral geometry, while merged pores are polyhedral (Fig. 4).This partition also reflects the contribution of pore volume, as unmerged pores have relatively small pore volume (and vice versa).
Contracting and dilating pores are delineated by considering the incremental change in pore volume.This is calculated by performing the classical Delaunay tessellation at time step t + Δt, and applying this tessellation to the previous time step t.If the incremental time step Δt is small, then the error in applying the same Delaunay tessellation to two time steps is minimal.The incremental volume change for a given pore is simply the summation of volume change associated with the corresponding Delaunay cells.
A scalar measure of anisotropy in pore orientation is defined by: The evolution of a k p for both dual-partitions of the pore space is shown in Figure 5 for the CP-D simulation.Unmerged pores have a preferential orientation perpendicular to the loading direction (Fig. 5a), while merged pores have a slight degree of anisotropy in the direction of loading (Fig. 5b).This observation can be explained by considering the particle contacts along the edges of the polyhedral pore cell.Merged pores (which have large pore volumes) can only be sustained by particle contacts with substantial interparticle normal forces.It has been observed that the average interparticle normal force increases with pore volume.Section 3 demonstrated that contacts with relatively large normal force magnitudes tend to align parallel with the loading direction and hence merged pores will align in this direction.Further, Figure 5c shows that contracting pores are approximately isotropically orientated, while dilating pores tend to align perpendicular to the loading di-rection (Fig. 5d).This agrees with the intuitive notion of lateral expansion under compression in a drained triaxial test.
These observations suggests that anisotropy of pores and contacts are somehow related.This is confirmed in Figure 6, which shows the correlation between anisotropy parameters of contacts and pores, along with their respective Pearson correlation coefficient.Certain parameters (for example, a wη c vs. a u p ) exhibit a strong linear correlation throughout the entire two-way cycle.Other sets of parameters have strong correlation for the individual loading, unloading and reloading branches (for example, a wμ c vs. a u p ). Figure 6 suggests that links can be established between the contact and pore space networks.The anisotropy of unmerged pores, a u p , appears to be highly correlated with the anisotropy of weak contacts (a wμ c and a wη c ).This suggests that unmerged pores predominantly comprise weak contacts, which agrees with the above discussion.Of particular interest is the observations that the anisotropy of dilating pores is correlated with strong contacts, while the anisotropy of contracting pores is somewhat correlated with weak contacts.This provides a connection between deformation characteristics (which may be controlled by the behaviour of pores) and stress transmission (which is governed by contacts).This presents an interesting angle to approach the constitutive behaviour of granular materials and is the subject of ongoing research.

Conclusion
This study quantifies the anisotropy of contacts and pores in a numerical simulation of a drained triaxial test.The anisotropy of subsets of contacts and pores are considered.This revealed that the stress ratio can be described uniquely as a function of geometric anisotropy of strong and non-sliding contacts and that the anisotropy of pores can be related to the anisotropy of contacts.

Figure 1 :
Figure1: Stress-strain and volumetric response for numerical simulations.Note that the colour scheme for loading, unloading, and reloading, along with the line type distinguishing CP-L and CP-D samples is maintained for the rest of the paper, unless stated otherwise.

Figure 2 :
Figure 2: Evolution of geometric anisotropy for the interparticle contact network in the CP-D simulation.

Figure 3 :
Figure 3: Unique relationship between stress ratio and the anisotropy associated with strong and non-sliding contacts.

Figure 4 :
Figure 4: Example of an individual pore comprising two Delaunay cells.Note that the particle sectors have been removed leaving just the pore space.

Figure 5 :
Figure 5: Evolution of anisotropy in pore orientation for the CP-D simulation.

Figure 6 :
Figure 6: Correlation between anisotropy parameters for contacts and pores.The Pearson correlation coefficient is also presented.