Linking macro-scale yielding and micro-scale response

The nature of the link between macro-scale yield and particle scale interactions is poorly understood. This contribution addresses this knowledge gap using discrete element method (DEM) triaxial test simulations performed with a modified version of the granular LAMMPS code, employing a simplified Hertz-Mindlin contact model. Each sample comprised monodisperse spheres and these samples were sheared from different initial density states. The positions of the yield surfaces were identified from the overall load:deformation behaviour, just as in a laboratory test. The overall sample responses observed in these DEM simulations confirmed that the size and positions of the sub-yield surfaces largely depend upon the sample density. The particle-scale parameters monitored included the coordination number, the second order fabric tensor and the proportion of sliding contacts (the sliding fraction). A correlation was established between the variation in coordination number with strain and stiffness degradation with strain. The sliding fraction data information showed that, in contrast to conventional understanding, the particles start to slide relative to each other even within the linear elastic zone


Introduction
Jardine [1] proposed a framework to describe the different types of yield associated with the load:deformation response of soil. Central to this framework is the existence of three yield surfaces, named 1 , 2 , and 3 which divide behaviour into three main zones, namely linear elastic behaviour (zone 1), non-linear elastic behaviour (zone 2) and elastoplastic behaviour (zone 3) respectively. Although hypothetical fundamental explanations for the mechanics underlying these yield surfaces have been proposed [2], the nature of the particle scale interactions that must determine these transition points is still not well-understood.
This con0tribution uses DEM (the Discrete Element Method) to addresses this knowledge gap focussing on the effect of specimen density. Triaxial compression tests are simulated and key particle-scale interactions are related to the overall material response. The particle-scale data obtained were related to the macro-scale yield points. The particlescale parameters monitored were the coordination number, the second order fabric tensor and the proportion of sliding contacts (the sliding fraction).

Simulation approach
Triaxial compression tests were carried out on samples with two different densities using a modified version of the granular LAMMPS code [3], with a simplified Hertz-Mindlin contact model. Each sample was created from an initial cloud of 34,986 noncontacting monodisperse spherical particles. The samples were enclosed in the vertical direction by two rigid, planar, horizontal boundaries; periodic boundaries were used in the lateral direction (Fig. 1a). To vary density, two values of inter-friction were used as the sample was isotropically compressed to an initial mean effective stress of ′ 0 = 1 ; the dense sample (RMD) was created using = 0, while = 0.15 was used to generate the loose sample (RML). These two samples were isotropically compressed to ′ 0 = 1000 using servo-control algorithm as described in [4]. The samples were then subject to strain controlled axial compression at a constant lateral stresses. A local damping coefficient of 0.1 was used and gravity was not applied during the triaxial compression. A summary of the DEM model parameters used in the simulations is outlined in Table  1. Both dynamic and static approaches were used to determine the shear modulus at specific strain levels. The deviator strain, was calculated as = 2/3( − ) where z and x are the major and minor normal strains. At a given stress state or strain level, the secant shear modulus, and tangent shear modulus, are given as: The shear stiffness was also dynamically measured using the shear wave velocity 'see [5]', , as follows: where is the sample density. The shear wave velocities were generated by exciting a rigid wall and are measured using the peak-peak method as illustrated in Fig. 1.

Identification of yield loci
The macro-scale yield points were determined using the procedure outlined in [1] and [6]. The yield point 1 was identified as the point at which a clear reduction was observed in the − relationship ( . . = / ; and = ′ z − 0.5 ( ′ y + ′ y )). Following Gasparre et al. [7], the point at which the relationship between deviator stress and axial strain become non-linear was located. The location of this point was then checked by unloading and ensuring that the resultant same stress-strain relationship coincided with the loading response (see Fig. 2).  Referring to Fig.3, where a plot of shear strain versus volumetric strain is shown for the RMD sample, the point at which the slope changes is taken as Y 2 . The yield point 3 , the point at which there is an abrupt change in the relationship between and ∆p′, (Fig  4a) was considered following documented laboratory work [6]. Referring to Fig. 4b, 3 coincides with a distinct kink in the − relationship which is widely used to map the large-scale yielding [8].    Fig. 6 illustrates the variation in both static and dynamic stiffness with shear strain; shear strain is equivalent to deviator strain. The three zones of behaviour outlined above are illustrated as zone 1, zone 2 and zone 3 on Fig. 6. It is obvious that each method gives a different stiffness value, and these differences are greater for sample RML. Exploring the origin of the gap between two methods is outside the scope of this study. It is clear from Fig. 6 that is close to the values of for the RMD sample; however, there is a big gap between and observed for the RML sample.

Fig. 6
Linking macro-scale yielding to trends in normalised stiffness measured using dynamic and static approach: (a) RMD sample; and (b) RML sample.

Linking the macro-scale yielding and coordination number
If the total number of contacts is and the total of particles is , the coordination number is defined as = 2 / . However, many grains exist in the systems that do not bear any force. Hence Thornton [9] introduced the mechanical coordination number, * which considers the number of particles with no contacts, 0 , and the number of particles with one contact, 1 : * = 2 − 1 − 0 − 1 Fig. 7 Linking macro-scale yield with coordination number data * : (a) RMD sample; and (b) RML sample.
Referring to Fig. 7, which illustrates the variation in both * and with shear strain, * is slightly higher than over the strain range considered, as would be expected. The values are close as the material is monodisperse, however the difference between * and is larger for the RML sample than for the RMD sample. Both * and correlate closely with . There is no clear correlation between the yield points ( 1 , 2 , 3 ) and changes in either * or . * remains constant throughout zone 1, and then starts to reduce in zone 3, before 3 is attained. experiences a slight rise and then starts to reduce when the max is observed. No change in the pattern of * / reduction is observed at point 3 .

Linking the macro-scale yielding and sliding fraction
The sliding fraction is defined as the ratio of the number of sliding contacts to the total number of contacts. Sliding occurs when the tangential component of the contact force, exceeds the product of the normal force component, and, (i. e. | | > ). Fig. 8 shows that sliding fractions strongly depend on sample density. In the RML sample, a larger percentage of contacts are initially sliding at small strains (within zone 1) when compared with RMD. However, the maximum sliding fraction is higher in RMD sample (35%) than in RML sample (25%). After the peak of the sliding fraction is attained RMD sample experiences a more rapid drop in sliding fraction than RML sample. At the largest strain level considered both samples attain a sliding fraction of about 20%. It has been postulated that particles do not slide within zone 1 and that particles start to slide relative each other in zone 2 (for example, [2], [6]). The data obtained here show that, in contrast to these earlier ideas, the particles may slide relative to each other even within zone 1. As indicated in Fig. 8 the sliding fractions vary from 2% to 8% in zone 1 for the RML sample. In both samples there is a more rapid increase in sliding fraction in zone 2 (i.e. sliding fraction increases up to almost 10 %). The sliding fraction reaches its peak at the stress state close to the largescale yielding and then drops when the stress state moves towards or beyond the large-scale yielding ( 3 ).

Linking macro-scale yielding and the second-order fabric tensor
The second-order fabric tensor (Φ ij ) is widely adopted to quantify fabric for DEM-generated samples; Φ ij , is given by Satake [10]: where is the contact normal orientation for contact and is the total number of contacts in the granular assembly.
The changes in contact orientations that result in measurable variations in fabric of samples were assessed. As the samples were sheared along the stress path in which lateral stress is kept constant and only vertical stress is increased, and so the fabric variation Φ /Φ is of interest. The fabric is nearly unchanged so long as behaviour of system is elastic within zone 1 and zone 2; the ratio Φ /Φ decreases as the stress state move towards the yield 3.

Conclusions
Compression tests and planar wave propagation have been performed, for which the correlation between the macro-scale yielding and particle-scale parameters is investigated. DEM data confirmed that macro-scale yielding largely depends on the sample density. There is no simple link between the yielding point defined from the macro-scale stress deformation response and either the sliding fraction and or the coordination number. However, the sliding fraction data information showed that, in contrast to conventional understanding, the particles starts to slide relative to each other even within the linear elastic zone. The amount of sliding at the contacts achieves a maximum value at large strains (i.e. when large-scale yielding occurs).