High resolution incremental stress testing of crushable granular materials

The incremental behaviour of crushable sands is investigated by means of a Discrete Element Method (DEM) based analogue. The DEM sample is subjected to a comprehensive set of small stress perturbations in the triaxial plane in order to identify the basic mechanisms contributing to the corresponding strain response. Three contributions to incremental strains are distinguished: elastic, plastic-structural and plastic-crushing. The behaviour observed appears to be consistent with the classic tenets of elasto-plasticity. It is also shown that high resolution probing is required to identify significant features such as elastic anisotropy and irreversible effects on the tangent bulk and shear moduli. As a consequence, computational efficiency is therefore a must for numerical studies of incremental response.


Introduction
Experimental studies of the incremental behaviour of granular materials involve very significant difficulties. In a typical experimental program, a series of experiments are carried out, each starting from the same initial stress state. Each experiment then records the deformation response to a small increment of stress applied in a particular orientation in (usually triaxial) stress space.
Experimental challenges include the control of the tests, small strain measurements and the preparation of a series of identical samples (one for each probe). Consequently, the use of the discrete element method (DEM) numerical models is emerging as an alternative tool to advance the understanding of the incremental behaviour of granular materials. Prior studies include the 2D work of Bardet [1] and later 3D contributions such as the work of [2], [3]. While these contributions have provided useful information on the underlying constitutive behaviour of granular materials [4], they did not include the possibility of grain crushing. However, it is clear that, if the effect of grain crushing is also of interest, the experimental difficulties mentioned above are significantly increased, as well as the relative advantages of numerical alternatives.
In this contribution, which follows [5], [6], the incremental strain response of a crushable sand is investigated numerically using stress probes on a 3D DEM analogue. The small strain (elastic) stiffness matrix is obtained and the effect of crushing on the response is also analysed. Low resolution stress probes were sufficient to identify the stress dependency of pseudoisotropic elastic moduli G0 and K0 in results aligned with previous DEM studies [7] of wave propagation in granular soils. However, when attention is focused on more advanced constitutive features, such as (i) the appearance of stress-induced elastic anisotropy and (ii) the effects of particle breakage and particle sliding on the stiffness evolution of the irreversible response, higher resolution probes where required.

Contact model
The modelling approach adopted was outlined in [8]. Spherical particles whose rotation was inhibited were used to capture the rotational resistance that exists between non-spherical grains [9], [10]. The contact model follows uses a simplified Hertz-Mindlin formulation and Coulomb friction. The limit criterion at which breakage is activated for a given particle was formulated following [11]. A particle breaks if the force F, at any of its contacts is such that where lim is the limit strength of the material and AF is the contact area. Once the limit condition is reached, a particle, modelled as a sphere in the PFC DEM code used [12], will split into smaller inscribed tangent spheres. Using a simplified Hertz-Mindlin model to describe also the contact area eq. (1) results in    where r1 and r2 are the radii of the contacting spheres and Ei, νi are the Young's Moduli and Poisson's ratio respectively. Note that this breakage criterion does not involve exclusively the maximum force on the particle: there is a strong inbuilt dependency on the characteristics of the contacting particles. Full details of the model are reported in [8].

Numerical specimen and calibration
In this work, a cubic numerical specimen with a side length of 4mm is filled with about 10k randomly assembled rigid particles of different sizes. The particle size distribution (PSD) corresponding closely to that of Fontainebleau sand. Gravity is neglected in the simulations and the specimen boundaries were defined using smooth "wall" elements. Target stress values were attained by using a servo-control to adjust the wall positions. As the rigid boundaries were smooth, the principal axes of stresses and strains are coincident with the cube axes. The principal strains were calculated directly from the wall displacements, while the corresponding principal stresses were obtained from the boundary forces. The contact law parameters were calibrated by simulating drained triaxial compression of dense and loose specimens of Fontainebleau sand under 100 kPa confinement ( Figure 1a). The particle failure criterion parameters were calibrated by capturing the apparent yield point of the loading curve in a one dimensional high pressure compression test simulation ( Figure 1b) by [13]. The calibrated model parameters are summarized in Table 1.

Initial states for probing
15 initial states were selected for probing. 12 were reached by imposing four radial compression stress paths with stress ratios, = q/p′, of 0, 0.3 0.75 and 1 (3 points along each stress path in order to cover a wide domain in the triaxial (q/p′) plane. Additionally, 3 initial states were selected from a p′=0 triaxial compression. Figure 2 shows the stress path and the initial conditions in the compression plane of the radial compression tests while Figure 3 presents the stress path of the p′=0 triaxial compression. In the latter the iso-grading state index (IG) contours quantifies the amount of crushing that has taken place. The contours were obtained by tracking the PSD during the radial compression tests [10]. Table 2 and  Table 3 summarise the initial conditions of the radially compressed and constant p′ sheared numerical samples respectively. The meaning of the different stiffness parameters reported therein is clarified in section 3.1.

Fig. 2
Stress paths and initial states on the triaxial and compression plane for the radial compression tests.

Stress probing program
For axisymmetric (triaxial) loading, the number of independent stress and strain variables reduces to two, and a convenient graphical representation can be given in the Rendulic plane (Figure 4a). The stress probe magnitude is given by where 'z is the vertical stress, and the two horizontal stresses ('x and 'y) are equal for this axisymmetric case.
The stress probe magnitudes used here were less than 1.7% of the current principal stress. The stress probe direction is defined by the angle  between the horizontal axis and the stress increment vector as indicated in Figure 4a. The response is depicted using the incremental strain response envelope (RE) (Figure 4b).
The direction of strain increment is defined by the angle, between the horizontal axis and the strain increment vector. a) b) Fig. 4. a) Stress probes b) Response envelope to stress probes.
As described in [5], the incremental deformation observed,  can be decomposed into a reversible ("elastic") part  e , and two irreversible parts: pu ("plastic-uncrushable"), the irreversible strain increment in the case where no crushing occurs, and  pc , the crushing induced plastic strain increment. Therefore To achieve this decomposition, three simulations are required for each stress probe. The first simulation uses the crushable DEM model described above to measure . In the second simulation all the mechanisms responsible for plasticity (interparticle sliding, opening of contacts and particle crushing) are inhibited to give  e . In the final simulation only crushing is inhibited to measure  e + pu . The number of radial probes that are made from a given initial state controls the resolution with which the response envelope may be defined. In this work we first present results for low resolution probing and then for high resolution probing. For low resolution probing only 12 loading directions were considered. These directions are characterised by the  loading vectors reported in Figure 4a, which follow some classical loading programs (IC/E = isotropic compression/expansion; TC/E = 'triaxial' (axi-symmetric) compression/extension; DC/E = purely deviatoric compression/extension; RE/C = radial extension/compression) and are spaced -on average-at 45°. On the other hand, for high resolution probing 36 incremental directions, spaced at 10° are probed. Because 3 parallel probes are performed in each direction, a total of 108 probes are performed for each initial state thus explored.

Elastic (reversible) incremental response
The elastic response of soils is generally anisotropic, but the complexity with which this anisotropy can be addressed is constrained by each particular testing program, [14]. In this work all the initial states and probing directions are contained within the triaxial plane. In these circumstances the most general elastic stiffness response can be described by the formulation proposed by Graham & Houlsby [15] as and hence 1 1 In eq. (5) q='z-'x and p' is the mean effective stress. J is a parameter expressing the cross dependence of the shear strain on the mean pressure and the volumetric strain on the shear stress. If isotropy is assumed J is equal to 0 and the generalized bulk and shear moduli, K* and G* reduce to K0 and G0. In eq. (6) C1= 3G*/Det, C2= -J/Det, C3= K*/Det and Det=(K*3G*-J 2 ). a) d) Fig. 5. a) Stress probes b) Response envelope to stress probes represented in the triaxial and vol-dev planes respectively.
For triaxial conditions, it can be shown that the incremental elastic response envelope lies within an ellipse. The principal axes of this ellipse are rotated when the off-diagonal term J is different from 0. If isotropic behaviour is assumed, the principal axes of the ellipse are aligned with the volumetric and deviatoric strain axes. G0 and K0 can also be more simply calculated using the IC and the DC probes ( Figure 5) and the elastic component of the volumetric and deviatoric strains. Figure 6 illustrates the results of low-resolution probing for initial states B0, B03, B075 and B1 (refer to Table 2 for the initial conditions of these states). A theoretical isotropic elastic ellipse in the vol-dev plane has been fitted to the DEM results. Although there is some indication of anisotropy, the isotropic fit may be deemed acceptable in most cases.  Table 2 reports the corresponding G0 and K0 deduced from the isotropic fit, while Figure 7 represents  G0 and K0 as function of the mean effective stress for different stress obliquity values. Whilst G0 and K0 all increase with p' only the Poisson ratio appears to have a clear trend with respect to the stress obliquity ratio  Fig. 7.  G0 and K0 as a function of p'.
A close look to Figure 6 suggests that the REe of the states characterised by  appears to be slightly rotated, meaning that the elastic response is slightly anisotropic (hence J≠0). To investigate further this characteristic of the elastic material response, states A B and C (characterised by the same value of p') in Figure 3 were subject to high resolution probing. The elastic response obtained is reported in Figure 8. At this level of resolution, the gradual development of elastic anisotropy as increases can be tracked. Although not explored further here, it is noted that this elastic anisotropy evolution may be linked to the development of internal state variables, such as the deviatoric fabric dev=1 - which, as reported in [6], shows a similar trend during loading.    (7) and therefore, by comparing the slope of the REp and REc in the vol-dev plane, it is possible to infer that (i) both plastic mechanisms have a unique flow direction and (ii) the dilatancy induced by the plastic-uncrushable mechanism is smaller than the dilatancy component of the plastic-crushable one. This means that the plasticuncrushable mechanism is more deviatoric than the plastic-crushable one (i.e. the plastic-crushable is more volumetric). By plotting the incremental response in the q-dev and p'-vol planes ( Figure 10) it easier to observe how plastic mechanisms affect the stiffens of the material response. The stiffness decreases as you move from an elastic to an elasto-plastic and elasto-plastic-crush response. In particular, the crushable response is always more compliant than the elasto-plastic-uncrushable. Crushing appears to have a greater effect on the bulk modulus than on the shear one. In fact the volumetric plastic incremental deformations due to the plasticcrushable mechanism (c) are very similar in magnitude to those of the plastic-uncrushable mechanism (p). On the other hand, the plastic deviatoric incremental deformations are smaller for the crushing induced mechanism (c) than for the uncrushable mechanism (p). These observations are in line with the dilatancy effect described above and were possible to visualize only when using the high-resolution stress probes.

Conclusions
The incremental behaviour of crushable granular soils was investigated using DEM. Axisymmetric stress probes with varying direction in stress space were applied starting from various stress states. A procedure was proposed to separate the reversible and irreversible components of total strain increments when particle breakage occurs. A detailed analysis of the results shows that: i) The strain response envelop obtained by stress probing a DEM model where interparticle sliding is inhibited, can be used to obtain the full elastic stiffness matrix. ii) Low resolution elastic stress probes on the are sufficient to identify have a good estimate of G0 and K0 which are found to be dependent on the confinement pressure and follow the same trends obtained in other DEM studies using wave propagation techniques. iii) High resolution response envelopes composed by 36 incremental directions where required to characterise the level of elastic anisotropy. iv) The level of elastic anisotropy was found to be directly proportional to the deviatoric fabric. v) Stress probing when the material is allowed to deform irreversibly indicates the development of a component of strains causing the REs to change shape. vi) The direction of the plastic flow induced by either interparticle sliding or particle crushing is independent of the loading direction. vii) The dilatancy induced by the plastic-uncrushable mechanism is smaller than the dilatancy component of the plastic-crushable one (i.e. the plasticcrushable is more volumetric).