Analysis of failure occurrence from direct simulations

ABSTRACT Stress probes are simulated with the discrete element method (DEM). From these simulations we show that the numerical discrete model presents a non-associative flow rule. As the material is non-associate, sign of the second-order work is checked for stress states included within the plastic limit condition. Conditions of occurrence of failure when the second-order work vanishes are discussed.


Introduction
Classically in geomechanics, failure phenomenon is associated with the existence of limit stress states. Experimental observations show that these limit stress states cannot be exceeded. If an operator tries to impose an additional stress loading increment from a limit stress state, the strain response will be very "large" (or even "unbounded"), the material fails. Notions of limit state and failure have been generalized by Darve et al. (2004) for stress states strictly included within the plastic limit condition. In the framework of rate-independent materials, generalized limit states are interpreted as homogeneous bifurcation states: according to the control parameters chosen, different responses are possible from the bifurcation state, some of the responses could correspond to the failure of the material. These bifurcation states (or equivalently, mixed limit states) are detected when the second-order work is zero or negative: As shown by Darve et al. (2004), d 2 W is essentially a directional quantity. In other words, a value of d 2 W does not correspond to a given mechanical state, but corresponds to a given loading direction (in stress or strain space) from a given stressstrain state. In addition, d 2 W can vanish inside the plastic limit condition only for non-associated materials, for associated materials d 2 W vanishes on the plastic limit condition. Darve et al. (2004) have shown that a whole bifurcation domain in stress space (including stress states for which d 2 W ≤ 0 for one or several stress directions) exist for the Hostun sand.
From different basis and approaches, Nicot et al. (2007), with the notion of "sustainability", and Nova (1994) (or Imposimato et al. (1998), with the notion of "controllability", also discussed the occurrence of failure for stress states inside the plastic limit condition. All these approaches are not detailed in this paper and we encourage readers to refer to the quoted papers. However these approaches converge toward identical conditions to be fulfilled together, for the occurrence of failure: (i) stress state belongs to the bifurcation domain; (ii) loading direction is characterised by d 2 W ≤ 0; (iii) control parameters are the ones allowing to trigger the failure. These conditions have been deduced from analytical analyses, and their necessity in failure occurrence has not been systematically verified. Therefore, we present in this paper a numerical verification of these conditions, from direct simulations based on the discrete element method (Cundall et al., 1979). Non-associativeness of the numerical discrete model is first discussed, since it is a necessity for vanishing of d 2 W within the plastic limit condition.

Numerical discrete model
Direct simulations were performed on a numerical discrete model with the code SDEC developed by Donzé et al. (1997), and based on the discrete element method Analysis of failure occurrence 189 Cundall et al. (1979). The 3D numerical model has a cubical shape, it is composed of about 10,000 polydisperse spheres (Figure 1a). Let's denote d s the sphere diameter, then d max s /d min s = 4.75. Spheres are rigid, they can slightly overlap according to an interaction contact law presented in Figure 1b. In the direction normal to the tangent contact plane, the relation is purely elastic and no tensile force is allowed. In the direction included in the tangent contact plane, the relation is elastic perfectly plastic. Only three mechanical parameters are introduced in the model at the contact scale; the normal k n and tangential k t stiffnesses, and the friction angle ϕ c . From sphere positions and the interaction contact law, contact forces are computed at each time step. Then, new sphere positions and orientations are determined by applying the Newton's second law from resultant forces and torques acting on each sphere. An explicit scheme is used to integrate the Newton's second law. The macroscopic stress-strain state of the numerical model (or numerical sample) is imposed through six frictionless walls whose positions are controlled at each time step to follow the prescribed loading programme. Since there is no tangential forces in wall-grain contacts, stress and strain principal directions coincide with the normals to the walls (Figure 1a). Each principal value of stresses or strains can be controlled; either directly for strains by adjusting the wall displacements, either indirectly for stresses thanks to a closed-loop control (since only wall positions or displacements are controlled). Consequently, all loading programmes defined with principal values of stresses (σ 1 , σ 2 , σ 3 ) and/or strains (ε 1 , ε 2 , ε 3 ) can be applied. The strain state is determined from wall positions and the stress state from wall-grain contact forces.
In this paper, analyses are carried out from simulations performed on two numerical samples E1 and E3. The samples differ only in their initial density, their characteristics are given in Table 1. During an axisymetric triaxial compression (characterised by a constant radial stress), the densest sample E1 is dilatant, whereas the sample E3 is essentially contractant. All simulations presented hereafter are performed in axisymetric conditions (σ 2 = σ 3 and ε 2 = ε 3 ), thus stress states and strain states can be represented in the Rendulic (or axisymetric) planes of stresses (σ 1 , √ 2 σ 3 ) or strains (ε 1 , √ 2 ε 3 ).

Non associative flow rule
To investigate the incremental constitutive behaviour of the numerical samples, it is useful to plot the response envelopes to strain or stress probes as suggested by Gudehus (1979). Stress probes, defined in the Rendulic plane of stress increments by the norm d σ = (dσ 1 ) 2 + (dσ 3 ) 2 and the orientation angle α (see Figure 2a), are performed from initial stress-strain states reached after drained triaxial compressions. The initial stress states are characterised by the confining pressure σ 3 and the shear stress ratio η = q/p (with q = σ 1 − σ 3 and p = (σ 1 + 2σ 3 )/3). The norm of the incremental stress loading d σ is set equal to 1 kPa (i.e. of the order of 10 −2 times the mean pressure p imposed for the simulations). Since the response of the numerical model to a given loading is non-linear, the choice of the norm of d σ is a tricky point. d σ should be chosen as small as possible to limit the influence of the non-linearity of the strain response on the shape of the response envelope. Nevertheless the response of a discrete element model to a given loading is marked by successive quasi-static phases followed by inertial phases (related to a modification of the contact network) (Combe et al., 2000). Hence d σ should be large enough to include some of these inertial phases (Sibille, 2006) directly involved in the macroscopic constitutive behaviour of the numerical model. For each stress probe, orientation α is constant, and strain response is computed and characterised in Rendulic plane of strain increment by d ε = (dε 1 ) 2 + (dε 3 ) 2 and the orientation β (see Figure 2b). The envelope of the set of strain responses d ε computed for 0 ≤ α ≤ 360 deg (by step of 10 deg) from a given initial stress-strain state is called the response envelope. Such a response envelope is shown with crosses in Figure 3.
In the following, the computed response envelopes are discussed by making reference to the classical elasto-plasticity framework. Thus, we assume: -the classical decomposition of strains into an elastic part d ε e and a plastic part d ε p such as: -the existence of an elastic limit surface f , of normal n, separating an elastic tensorial zone and an elasto-plastic tensorial zone, -the existence of a plastic potential g, of normal m.
To compute the elastic responses to stress probes, the sliding at intergranular contacts is inhibited by imposing ϕ c → 90 deg. For the size of the stress increments imposed, Analysis of failure occurrence 191 this later condition is confirmed to be sufficient, by verifying the vanishing of residual strains after a cycle of loading/unloading, and by checking the lack of sliding and opening contacts. Once the elastic strain response is determined, the plastic strain response is computed from Equation In Figure 3 are presented a total response envelope, computed with sample E1, and its decomposition into elastic and plastic response envelopes. The elastic response envelope with an elliptical shape centered with respect to the axis origin is typical of a purely elastic behaviour (Gudehus, 1979). The plastic response envelope (represented with black points) forms a straight line meaning that the direction β p of plastic strain increments d ε p is constant and independent of stress loading directions α (for a given initial stress-strain state). Therefore, this plastic response envelope shows the existence of a flow rule for the numerical samples used. Another way to represent a strain response envelope consists in plotting the norm and the direction of strain responses versus the stress loading direction. Figures 4a and 4b display such representations for elastic and plastic strain response envelopes computed with sample E3. Figure 4b shows clearly that β p is quasi constant (for d ε p = 0) confirming the existence of a flow rule, and with β p = 129 deg corresponding to the direction of the normal m to the plastic potential.
In Figure 5 are presented the stress directions corresponding to an elasto-plastic response in the sense of classical elasto-plasticity. By following these directions counter-clockwise, the first and the last directions are tangent to the elastic limit surface f . Right and left tangents to the elastic limit surface can be deduced from Figure 4a: first plastic strains are found for α = 60 deg and last ones for α = 240 deg.
Right and left tangents are collinear and the normal n to the elastic limit surface is aligned along the direction α = 150 deg. This latter direction corresponds to the maximum of d ε p as shown in Figure 4a. By this way we obtain a rough estimation of the direction of n, but this estimation is enough to conclude about the non-associativeness of the numerical model (see also Bardet (1994) and Calvetti et al. (2003)). Therefore the second-order work should vanish for stress states included in the plastic limit surface.  Figure 6 shows directions of n and m for different stress states and for samples E1 and E3. We remark for high value of η that the elastic limit surface tends to align with the Mohr-Coulomb criterion (plotted with a thin continuous line).

Occurrence of diffuse failure
Since the second-order work d 2 W is essentially a directional quantity, it is useful to perform stress probes (as defined in Section 3) to check values of d 2 W with respect to stress loading directions. Values of second-order work are normalized to allow comparisons between different stress-strain states and different sample densities. The normalized second-order work d 2 W norm is defined as (Darve et al., 2000): d 2 W norm is equal to the cosine of the angle between d σ and d ε, consequently from Figure 5 it appears clearly that d 2 W norm (and thus d 2 W ) can take negative values before the plastic limit condition only if the material is non-associated.  Figure 7 shows a circular diagram of d 2 W norm versus the stress direction α computed with samples E1 and E3 for a confinement of 100 kPa. In a such diagram, an arbitrary constant ρ is added to polar values of d 2 W norm to verify: A dashed circle drawn in Figure 7 represents vanishing values of d 2 W norm . Inside the dashed circle d 2 W norm is negative, outside it is positive. For both samples, we remark negative values of d 2 W norm (or d 2 W ) for stress directions α grouped in a cone of unstable stress directions (Imposimato et al., 2001;Darve et al., 2004). Cones shown here are the first ones found with respect to the shear stress level η (η = 0.82 and 0.46 for E1 and E3 respectively, below these η values d 2 W > 0 for all stressstrain states and stress directions checked (Sibille et al., 2007;Sibille et al., 2008)). These shear stress levels (corresponding to mobilized friction angles φ m = 21.0˚and 12.3˚) are included within the Mohr-Coulomb criterion (defined by φ = 24.7˚and 21.2˚for E1 and E3 respectively).
Through Figure 7 are summarized the three main characteristics of the influence of density on cones of unstable stress directions. (i) the cone is more opened for the loosest sample, (ii) for the loosest sample lower α values are included in the cone than for the densest sample, (iii) first cone is found for a lower value of φ m with the loosest sample than with the densest one. An identical qualitative influence of the density on cones has been highlighted by computations made with a phenomenological constitutive relation fitted on dense and loose Hostun sands (Darve et al., 2004;Sibille et al., 2007). This simple comparison of cones of unstable stress directions with respect to the density, corroborates the fact that, experimentally failure of a loose sand sam-Analysis of failure occurrence 195 ple from a stress state included within the Mohr-Coulomb criterion is much easier to observe than with a dense sample (the bifurcation domain and the range of unstable stress directions are wider). For simulations of stress probes, numerical samples were fully stress controlled. The loading programme was defined for a given stress direction α by the control parameters dσ 1 = cst 1 and dσ 3 = cst 3 (plus dσ 2 = dσ 3 to verify the axisymetric condition), where cst 1 and cst 3 are chosen such that the stress direction α is verified. The corresponding response parameters are dε 1 and dε 3 . With these control parameters vanishing or negative values of d 2 W are found but no failure was observed (after a stress probe the numerical sample reach a new equilibrium state).
We verify hereafter that the control mode plays a fundamental role in failure occurrence. The objective is to impose a given stress direction to numerical samples through different control parameters than previous one. The ratio between σ 1 and σ 3 can be imposed for instance through the condition dσ 1 − dσ 3 /R = 0, where R = cos α / √ 2 sin α for α ∈ ]180 deg; 270 deg]. By restricting our analysis to conjugated parameters in the sense of energy (Nova, 1994), then: Thus the loading programme can be defined for instance by: and the response parameters are dε 1 and dσ 3 /R. This loading programme ensures a mixed control mode (in stress and strain) with control parameters defined as linear combinations of principal stress or strain components. Such control modes are encountered in very classical laboratory test, for instance a drained triaxial compression corresponds to a mixed control mode (σ 3 = 0 and dε 1 > 0), and for an undrained triaxial compression a condition is imposed on a linear combination of principal strain components (dε v = dε 1 + 2 dε 3 = 0).
Numerically, to impose the loading programme defined by [6] to the sample, several time steps are necessary (even if increments of control parameters are "small"). At each time step t of the discrete element method (DEM) the following conditions are verified: where C t σ and C t ε are constants set at each step t such that conditions defined in [6] are verified when the simulation ends; σ t 1 , σ t 3 , ε t 1 and ε t 3 correspond to the stress-strain state computed at step t. Figure 8 presents a simplified diagram of the numerical loop followed, involving conditions defined in [7]. Full details can be found in Sibille (2006).

Figure 8. Principle of the loop to control linear combinations of principal stress or strain components
In a first time we consider the stress direction characterised by R = 1 (α = 215.3 deg). This direction is included in the cone of unstable stress directions for sample E3 at the stress state considered (Figure 7). For R = 1 control parameters defined in [6] write: dq = 0 and dε v < 0 (dilatancy imposed). Thus the sample follows a constant shear stress path (for instance, such a stress path can be followed by a material point in a slope when pore-water pressure increases: dq = 0 and dp < 0). Figure 9 shows a comparison of responses of E3 between a full stress control and a control defined by [6] in terms of kinetic energy and stresses. The kinetic energy of the sample is computed as the sum of the kinetic energy of each grain. When the sample is fully stress controlled the kinetic energy stays very low and finally vanishes, the sample reaches a new equilibrium at a stress state close to the initial one. When the sample is controlled by dq = 0 and dε v < 0 a burst of kinetic energy is computed characterised by a monotonic increase until values of two order of magnitude larger than maximum kinetic energy computed during fully stress controlled probes. We observe a sudden vanishing of stresses (as a sudden liquefaction). Due to the dynamic response, the importance of inertial terms can by visualized in the evolution of the axial σ 1 and radial σ 3 stresses; when σ 3 vanishes, σ 1 is about 40 kPa.
This result verifies that failure occurrence from a bifurcation point, detected with the second-order work criterion, depends on the choice of the control parameters. For instance, when the sample is fully stress controlled no failure can occur before the plastic limit condition as shown by Nova (1994). Darve et al. (2004) have shown that generalized limit states exist strictly inside the plastic limit condition. For the particular case of constant shear stress paths on loose sand, the authors have shown that a maximum dilatancy cannot be exceeded (as a stress state corresponding to the plastic limit condition cannot be exceeded). The choice of the loading condition dε v < 0 has been motivated by this demonstration. We verify, for this very particular loading programme, that effectively a generalized limit state cannot be exceeded (note that if dε v > 0 was imposed no failure would be observed, in the same way as a "stress unloading" applied from a plastic limit state does not lead to the failure).
In the following, we pay attention to the importance of the stress direction, together with the influence of density, on failure occurrence. Four stress directions are considered and their belonging to cones of unstable stress directions for sample E1 and E3 (see Figure 7) are summarized in Table 2. The control parameters are those defined in expression [6], but no variation of loading parameters is imposed. Thus the loading programme is defined by: and apply to numerical samples through the algorithm presented in Figure 8.
If simulations are run in these conditions there is no evolution of the samples, they stay at their initial mechanical state governed by the prescribed control parameters. A perturbation is necessary to conclude about the importance of the stress direction. Simulations are performed without gravity. Therefore, at a given mechanical state, some spheres float within the pores. Samples are perturbed by imposing, at a given time step t, an instantaneous velocity in a random direction to eight floating grains. Samples are virtually split into eight sub-parallelepipeds, each perturbed floating grain is chosen randomly in each sub-parallelepiped respectively. The velocity imposed to each grain is computed such that the value of kinetic energy provided is equal for each grain. The total value of external input of kinetic energy is 10 −5 J, which is small with respect to the maximum value of kinetic energy computed for fully stress controlled probes: 10 −4 J.  Responses of samples E1 and E3 to the perturbation with respect to R value are presented in Figures 10 and 11 in terms of kinetic energy and radial stress, σ 3 , evo-Analysis of failure occurrence 199 lutions. For R = 1.94 and 0.257 (α = 200 and 250 deg respectively) samples E1 and E3 reach a new equilibrium at a stress state close to the initial one. These directions are not included in cones for both samples. For R = 0.843 (α = 220 deg) failure is observed only for sample E3, this direction is not included in the cone for sample E1. For R = 0.593 (α = 230 deg), failure is observed for both samples, this direction corresponds to negative values of d 2 W for both samples. However the burst of kinetic energy is not obvious for sample E3, for which direction α = 230 deg is close to the limit of the cone. These results verify that failure can occur only along loading directions included in the cone of unstable stress directions. These unstable stress directions are closely related to the density of the granular material.

Conclusion
Since the discrete element method involves few hypotheses, it constitutes a good tool to verify prediction deduced from analytical developments. Moreover, for very particular loading programmes, as those defined in this paper, simulations are easier to realize, than real experiments. The presented results confirm that from bifurcation points detected with the second-order work criterion, for stress directions included in cones of unstable stress directions, and for particular control parameters, failure can develop suddenly, characterised by a dynamic response of the material and the vanishing of stresses. These failures occur from stress states inside the Mohr-Coulomb criterion and could explain, for instance, landslides not predicted by a classical approach.