Failure in granular materials based on acoustic tensor: a numerical analysis

We investigate localization in granular material with the support of numerical simulations based upon DEM (Distinct Element Method). Localization is associated with a discontinuity in a component of the incremental strain over a plane surface through the condition of the determinant of the acoustic tensor to be zero. DEM simulations are carried out on an aggregate of elastic frictional spheres, initially isotropically compressed and then sheared at constant pressure p0. The components of the stiffness tensor are evaluated numerically in stressed states along the triaxial test and employed to evaluate the acoustic tensor in order to predict localization. This occurs in the pre-peak region, where the aggregate hardens under the circumstance to be incrementally frictionless: it is a regime in which the tangential force does not change as the deformation proceedes and, consequently, the deviatoric stress varies only with the normal component of the contact force.


Introduction
In a recent contribution, La Ragione et al. [1] analyze and predict localization in an ideal granular material based upon micro-mechanical considerations. They point out that failure is associated with a regime of deformation in which the incremental tangential forces become negligible. It is the possibility to describe the phenomena at particle level, including a more sophisticated kinematics between particles and their equilibrium [2,3], that allows to predict when localization occurs and under which conditions. This is done in the context of the model proposed by Rudnicki and Rice [4], Vardoulakis [5,6], in which localization is associated with the vanishing of the acoustic tensor. A different approach is proposed by Nicot and coauthors [7,8] as they focus on the vanishing of the secondorder work and the possibility to predict diffuse bifurcation. This phenomena may precede localization [9] but both second-order work and determinant of the acoustic tensor are zero when localization occurs. Here we focus on the acoustic tensor and evaluate the effective moduli via numerical simulations, following Recchia et al. [10]. We show that localization is plausible in a regime of deformation where the incremental material behavior is governed by normal forces.
In the present contribution, we briefly review the theory and the essential condition for localization. Next we employ numerical simulations to mimic a triaxial compression test for an ideal granular material and we evaluate the components of the stiffness tensor in stressed states along the main path. Finally, we use such components to verify under which conditions localization may occur.

Theory
We focus on an ideal elastic frictional aggregate of particles that has been isotropically compressed and next sheared. We label with y 1 the vertical axis of compression while in the orthogonal plane we consider the other two directions y 2 and y 3 such that y 1 , y 2 and y 3 is a proper orthogonal, rectangular Cartesian coordinate system. The loading is axially symmetric so we restrict our attention to the y 3 -y 1 plane.
Following [4], localization occurs when the determinant of the acoustic tensor vanishes: where l l l = (0, -sin(Ψ), cos(Ψ)) is the unit vector normal to the plane of discontinuity, see Fig. 1. The stiffness tensor, A i jqp , relates the increments in stress to the increments in total average strain as and assumes the form Because of the loading condition, the aggregate is transversely isotropic and depends on the six coefficients η i [11]. Eq. (1) is equivalent to (T 1 cos 4 Ψ + T 2 cos 2 Ψ + T 3 )(η 6 cos 2 Ψ + η 3 ) = 0, which leads to the solutions [1]: with and Localization occurs when the discriminant in Eq. (5) is zero. Then, the angle of localization is In the next section we carry out numerical simulations to evaluate the components A i jqp of the stiffness tensor of the aggregate.

DEM simulation
DEM represents granular materials as packings of solid particles (often spheres) that are allowed to overlap [12]. The particles in contact interact with their neighbors via repulsive springs, resulting in momentum exchange between particles. If the contact forces, acting on a particle, are known the problem is reduced to the integration of Newton's equations of motion for the translational and rotational degrees of freedom of that particle.
The system considered here is a random assembly of identical, frictional, elastic spheres in which gravity is neglected. Particles interact through a non-central force with Normalized deviatoric stress, q/p 0 (black symbols), versus normalized deviatoric strain, γ/p 0 , and its components q N /p 0 (magenta symbols) and q T /p 0 (blue symbols). a normal component that follows the non-linear Hertz law and a tangential component that includes a bilinear relation characterized by an elastic resistance followed by a Coulomb sliding. During the simulation, the average stress is calculated following the Cauchy relation for molecular systems (e.g. [10]).

Preparation protocol
We refer to material characteristics typical of glass spheres, with shear modulus G s = 29GPa, Poisson's ratio, ν = 0.2, and radius R = 0.1mm. At the beginning of the simulation, particles are randomly generated in a periodic cubic cell and then isotropically compressed to achieve a consolidated state. Compression happens in two stages. Initially the sample is compressed with friction coefficient µ = 0, until a solid volume fraction φ 0.64 is reached; then, after relaxation, a second compression follows with µ = 0.5, that brings the aggregate to the target isotropic pressure p 0 = 200kPa, with φ 0.64 (see [13] for details). In this reference isotropic configuration, the volumetric strain associated with the pressure p 0 is ∆ 0 = 1 × 10 −3 .

Axial-symmetric compression
After compression, the aggregate is sheared by applying strain along the axial direction y 1 , while the pressure p 0 = 200kPa is kept constant by means of a servomechanism [13]. The test is carried out with µ = 0.5. To ensure quasi-static conditions, the compression is performed with a sequence of small strain steps of the order of ∆ 11 = 4 × 10 −6 followed by relaxation steps. In the latter, particles are allowed to dissipate kinetic energy and reach intermediate equilibrium states.
Along the compression path, we can extract the deviatoric strain as γ = ( 22 + 33 )/4 − 11 /2 and measure the deviatoric stress q = (σ 22 + σ 33 )/2 − σ 11 , as well as its normal and tangential parts, q N and q T , associated, with the normal and tangential contact forces, respectively [14]. In Figure 2 we plot the normalized deviatoric stress q/p 0 against the normalized deviatoric strain γ/∆ 0 . In the same figure, the partitions q N /p 0 and q T /p 0 are shown as well. After an initial stage, it is clear that the deviatoric stress is almost identifiable with its normal component q N , with the increments in q T almost negligible. In particular, for γ/∆ 0 0.3 the tangential part of the deviatoric stress remains constant, i.e.q T = 0. This regime is characterized by a constant tangential force so that the resulting incremental force is only associated with its normal component. La Ragione et al. [1] define this regime incrementally frictionless. As well known in elasto-plasticity theory material yields both before and after the peak stress. We are in the pre-peak region in which plasticity is associated with q T constant while q N still increases indicating an hardening behavior of the aggregate untilq N = 0 [14], [15].

Stiffness tensor
Along the compression path, we consider different stress states and evaluate the non zero components of the stiffness tensor. As shown in Recchia et al. [10], for each point, we apply strain perturbations ∆ pq in the forward direction, i.e. consistent with loading along the triaxial stress path, so that if particles were sliding or rolling during the axial loading, they continue to do so during the probe. Then we measure the change in stress ∆σ i j between the final state (after sufficient relaxation) and the stress before the perturbation. Finally, the stiffness components are given as For example, the moduli A 1111 and A 3311 are obtained by applying an infinitesimal perturbation along y 1 , ∆ 11 0, with ∆ 22 = ∆ 33 = 0, that leads to stress change along y 1 and y 3 , thus A 1111 = ∆σ 11 /∆ 11 and A 3311 = ∆σ 33 /∆ 11 . (11) Similarly, the moduli A 1133 and A 3333 are calculated with a perturbation along y 3 , with ∆ 11 = ∆ 22 = 0, as A 1133 = ∆σ 11 /∆ 33 and A 3333 = ∆σ 33 /∆ 33 . (12) Finally, we verify that due to axial-symmetry.
As an example, in Figure 3 we show the evolution of the stiffness component A 1111 with the strain amplitude. As reported in [10,16,17], the first plateau refers to the elastic response in which the perturbation is so small to prevent sliding and rolling among particles. The second plateau is associated with an inelastic, yet linear, incremental response, that preserves the micro-and macromechanisms happening during axial loading. We choose the latter as reference to build the acoustic tensor ( i.e., we employ values of A i jkh obtained with strain perturbations ∆ i j 5 × 10 −5 ) because we refer to probes that are similar to those considered during the unixial loading, when localiaztion may occur.

Results: acoustic tensor
With the knowledge of the A i jkh components at several stress states along the triaxial deformation, we can evaluate T 1 , T 2 and T 3 and calculate the discriminant of the acoustic tensor in Eq. (5), T 2 2 − 4T 1 T 3 , and its evolution with the strain.
In Figure 4 we plot the discriminant normalized by the particle shear modulus G s versus the normalized shear strain γ/∆ 0 , evaluated in five stresses states, as indicated also in Figure 2. The discriminant grows monotonically and changes sign for γ/∆ 0 > 0.35 suggesting that localization may occur. It is interesting to note that the possibility of localization, thus failure, is here predicted in a region of deformation that precedes the stress peak. This is the regime in which q T is constant while q N keeps growing. Here, Recchia et al. [10] find, at contact level, that particles slide and roll confirming no increment in the tangential forces and supporting the idea that the increment in the contact forces occurs only along the normal component. It is expected that in a later stage shear bands will occur with major fabric changes and dilatancy, with particles failing into (larger) gaps steadily. However such analysis is beyond the scope of the present study.
In the context of micro-mechanics, the recent contribution by Karapiperis et al. [18] provides a useful tool to have a detailed analysis at particle level and test theory like this. It may be also relevant to test how the shape of the particles in the aggregate, here idealized as identical elastic spheres, influences the onset of localization [19].

Conclusion
Along a triaxial test on an aggregate of identical, elastic, particles we have studied the evolution of the acoustic tensor and the possibility that the condition for localization is met. The analysis has been conducted by means of DEM numerical simulations. The components of the stiffness tensor have been evaluated by probing stressed, anisotropic, states via strain increments. We show that the determinant of the acoustic tensor vanishes before the peak of deviatoric stress is reached [20], in a regime characterized by the incrementally frictionless behavior of the aggregate. A sequel is in progress to visualize per-particle incremental shear strain, that should indicate the presence of a localization band in the same stress-strain region where failure is predicted.