# Analysis of failure occurrence from direct simulations

Luc Sibille, François Nicot, Frédéric-Victor Donzé, Félix Darve

► To cite this version:

Luc Sibille, François Nicot, Frédéric-Victor Donzé, Félix Darve. Analysis of failure occurrence from direct simulations. European Journal of Environmental and Civil Engineering, 2009, 13 (2), pp.187-201. 10.1080/19648189.2009.9693099. hal-01005248

HAL Id: hal-01005248

https://hal.science/hal-01005248

Submitted on 21 Feb 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analysis of failure occurrence from direct simulations

Luc Sibille* — François Nicot** — Frédéric-Victor Donzé***
Félix Darve***

* Laboratoire GeM, Université de Nantes, ECN, CNRS
IUT de Saint-Nazaire, BP 420, F-44606 Saint-Nazaire cedex
luc.sibille@univ-nantes.fr

** Cemagref de Grenoble, Unité ETNA
BP 76, F-38402 Saint-Martin-d'Hères cedex
francois.nicot@cemagref.fr

*** Laboratoire 3S-R, INPG, UJF, CNRS
BP 53, F-38041 Grenoble cedex 9
{Frederic.Donze, Felix.Darve}@inpg.fr

**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.

**RÉSUMÉ.** Des recherches directionnelles en contrainte sont simulées avec la méthode des éléments discrets (MED). A partir de ces simulations, nous montrons que le modèle numérique discret présente une règle d'écoulement non associée. Par conséquent, le signe du travail du second ordre est calculé pour des états de contraintes inclus à l'intérieur de la condition limite de plasticité. Les conditions de développement d'une rupture lorsque le travail du second ordre s'annule sont discutées.

**KEYWORDS:** failure, second-order work, control parameter, flow rule, stress probes, discrete element method.

**MOTS-CLÉS :** rupture, travail du second ordre, paramètre de contrôle, règle d'écoulement, recherche directionnelle, méthode des éléments discrets.
# 1. 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:

$$d^2W = d\sigma : d\varepsilon \le 0 \quad [1]$$

As shown by Darve et al. (2004), $d^2W$ is essentially a directional quantity. In other words, a value of $d^2W$ does not correspond to a given mechanical state, but corresponds to a given loading direction (in stress or strain space) from a given stress-strain state. In addition, $d^2W$ can vanish inside the plastic limit condition only for non-associated materials, for associated materials $d^2W$ 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^2W \le 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 “stainability”, 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^2W \le 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^2W$ within the plastic limit condition.

# 2. 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
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_s^{max}/d_s^{min} = 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 $\varphi_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.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** a) numerical model with frictionless walls and frame axes; b) sketch of the intergranular interaction law</figcaption></figure>

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 ($\sigma_1, \sigma_2, \sigma_3$) and/or strains ($\varepsilon_1, \varepsilon_2, \varepsilon_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 axisymmetric 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 axisymmetric conditions ($\sigma_2 = \sigma_3$ and $\varepsilon_2 = \varepsilon_3$), thus stress states and strain states can be represented in the Rendulic (or axisymmetric) planes of stresses ($\sigma_1, \sqrt{2}\sigma_3$) or strains ($\varepsilon_1, \sqrt{2}\varepsilon_3$).
<table><caption>**Table 1.** *Characteristics of the numerical samples*</caption><thead><tr><th>Sample</th><th>k<sub>n</sub>/d<sub>s</sub><br>(MPa)</th><th>k<sub>t</sub>/k<sub>n</sub></th><th>φ<sub>c</sub><br>(deg)</th><th>Void ratio<br>e</th><th>Coordination number<br>z</th></tr></thead><tbody><tr><td>E1</td><td>356</td><td>0.42</td><td>35.0</td><td>0.618</td><td>4.54</td></tr><tr><td>E1</td><td>356</td><td>0.42</td><td>35.0</td><td>0.693</td><td>4.42</td></tr></tbody></table>

**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\vec{\sigma}|| = \sqrt{(d\sigma_1)^2 + (d\sigma_3)^2}$ and the orientation angle $\alpha$ (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 $\sigma_3$ and the shear stress ratio $\eta = q/p$ (with $q = \sigma_1 - \sigma_3$ and $p = (\sigma_1 + 2\sigma_3)/3$). The norm of the incremental stress loading $||d\vec{\sigma}||$ 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\vec{\sigma}$ is a tricky point. $||d\vec{\sigma}||$ 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\vec{\sigma}||$ 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 $\alpha$ is constant, and strain response is computed and characterised in Rendulic plane of strain increment by $||d\vec{\varepsilon}|| = \sqrt{(d\varepsilon_1)^2 + (d\varepsilon_3)^2}$ and the orientation $\beta$ (see Figure 2b). The envelope of the set of strain responses $d\vec{\varepsilon}$ computed for $0 \le \alpha \le 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\vec{\varepsilon}_e$ and a plastic part $d\vec{\varepsilon}_p$ such as:

$$d\vec{\varepsilon} = d\vec{\varepsilon}_e + d\vec{\varepsilon}_p \quad [2]$$

- the existence of an elastic limit surface $f$, of normal $\vec{n}$, separating an elastic tensorial zone and an elasto-plastic tensorial zone,

- the existence of a plastic potential $g$, of normal $\vec{m}$.

To compute the elastic responses to stress probes, the sliding at intergranular contacts is inhibited by imposing $\varphi_c \to 90$ deg. For the size of the stress increments imposed,
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 [2].

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Stress probes a) and strain responses b) defined in the Rendulic plane of stress and strain increments, respectively</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Decomposition of the strain response envelope into an elastic part and a plastic part; sample E1, $\sigma_3$ = 200 kPa and $\eta = 0.77$</figcaption></figure>

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 en-
velope with an elliptical shape centered with respect to the axis origin is typical of
a purely elastic behaviour (Gudehus, 1979). The plastic response envelope (repre-
sented with black points) forms a straight line meaning that the direction $\beta_p$ of plastic
strain increments $d\epsilon'_p$ is constant and independent of stress loading directions $\alpha$ (for
a given initial stress-strain state). Therefore, this plastic response envelope shows the
existence of a flow rule for the numerical samples used.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** a) norms of plastic, ||dε<sub>p</sub>||, and elastic, ||dε<sub>e</sub>||, strain responses; b) directions of plastic, β<sub>p</sub>, and elastic, β<sub>e</sub>, strains; sample E3, σ<sub>3</sub> = 200 kPa and η = 0.63</figcaption></figure>

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 βₚ is quasi constant (for
||dεₚ|| ≠ 0) confirming the existence of a flow rule, and with βₚ = 129 deg corre-
sponding 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 sur-
face *f*. Right and left tangents to the elastic limit surface can be deduced from Fig-
ure 4a: first plastic strains are found for α = 60 deg and last ones for α = 240 deg.
Right and left tangents are collinear and the normal $\vec{n}$ to the elastic limit surface is aligned along the direction $\alpha = 150$ deg. This latter direction corresponds to the maximum of $||d\vec{\epsilon}_p||$ as shown in Figure 4a. By this way we obtain a rough estimation of the direction of $\vec{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><img src="image_6.png" /><figcaption>**Figure 5.** Definition of stress directions leading to an elasto-plastic response in the framework of classical elasto-plasticity</figcaption></figure>

Figure 6 shows directions of $\vec{n}$ and $\vec{m}$ for different stress states and for samples E1 and E3. We remark for high value of $\eta$ that the elastic limit surface tends to align with the Mohr-Coulomb criterion (plotted with a thin continuous line).

## 4. Occurrence of diffuse failure

Since the second-order work $d^2W$ is essentially a directional quantity, it is useful to perform stress probes (as defined in Section 3) to check values of $d^2W$ 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^2W_{norm}$ is defined as (Darve et al., 2000):

$$d^2W_{norm} = \frac{d\vec{\sigma} \cdot d\vec{\epsilon}}{||d\vec{\sigma}|| ||d\vec{\epsilon'}||} \quad [3]$$

$d^2W_{norm}$ is equal to the cosine of the angle between $d\vec{\sigma}$ and $d\vec{\epsilon}$, consequently from Figure 5 it appears clearly that $d^2W_{norm}$ (and thus $d^2W$) can take negative values before the plastic limit condition only if the material is non-associated.
<figure><img src="image_7.png" /><figcaption>Figure 6. Comparison of directions of the plastic flow, $\vec{m}$, and of the normal $\vec{n}$ to the elastic limit surface</figcaption></figure>

Figure 7 shows a circular diagram of $d^2W_{norm}$ versus the stress direction $\alpha$ computed with samples E1 and E3 for a confinement of 100 kPa. In such diagram, an arbitrary constant $\rho$ is added to polar values of $d^2W_{norm}$ to verify:

$$ \forall \alpha, \quad d^2W_{norm}(\alpha) + \rho > 0 \qquad [4] $$

A dashed circle drawn in Figure 7 represents vanishing values of $d^2W_{norm}$. Inside the dashed circle $d^2W_{norm}$ is negative, outside it is positive. For both samples, we remark negative values of $d^2W_{norm}$ (or $d^2W$) for stress directions $\alpha$ 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 $\eta$ ($\eta = 0.82$ and 0.46 for E1 and E3 respectively, below these $\eta$ values $d^2W > 0$ for all stress-strain states and stress directions checked (Sibille et al., 2007; Sibille et al., 2008)). These shear stress levels (corresponding to mobilized friction angles $\phi_m = 21.0^\circ$ and $12.3^\circ$) are included within the Mohr-Coulomb criterion (defined by $\phi = 24.7^\circ$ and $21.2^\circ$ 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 $\alpha$ values are included in the cone than for the densest sample, (iii) first cone is found for a lower value of $\phi_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-
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).

<figure><img src="image_8.png" /><figcaption>Figure 7. Normalized second-order work, $d^2W_{norm}$, versus stress probe direction $\alpha$, for samples E1 and E3 at the same confining pressure $\sigma_3 = 100$ kPa</figcaption></figure>

For simulations of stress probes, numerical samples were fully stress controlled. The loading programme was defined for a given stress direction $\alpha$ by the control parameters $d\sigma_1 = cst_1$ and $d\sigma_3 = cst_3$ (plus $d\sigma_2 = d\sigma_3$ to verify the axisymmetric condition), where $cst_1$ and $cst_3$ are chosen such that the stress direction $\alpha$ is verified. The corresponding response parameters are $d\varepsilon_1$ and $d\varepsilon_3$. With these control parameters vanishing or negative values of $d^2W$ 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 $\sigma_1$ and $\sigma_3$ can be imposed for instance through the condition $d\sigma_1 - d\sigma_3/R = 0$, where $R = \cos \alpha / (\sqrt{2 \sin \alpha})$ for $\alpha \in ]180^\circ; 270^\circ]$. By restricting our analysis to conjugated parameters in the sense of energy (Nova, 1994), then:

$$ \sigma_1\varepsilon_1 + 2\sigma_3\varepsilon_3 = \varepsilon_1\left(\sigma_1 - \frac{\sigma_3}{R}\right) + (\varepsilon_1 + 2R\varepsilon_3)\frac{\sigma_3}{R} \quad [5] $$

Thus the loading programme can be defined for instance by:

$$ d\sigma_1 - \frac{d\sigma_3}{R} = 0 \quad \text{and} \quad d\varepsilon_1 + 2R d\varepsilon_3 < 0 \quad [6] $$

and the response parameters are $d\varepsilon_1$ and $d\sigma_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 en-
countered in very classical laboratory test, for instance a drained triaxial compression
corresponds to a mixed control mode ($\sigma_3 = 0$ and $d\epsilon_1 > 0$), and for an undrained
 triaxial compression a condition is imposed on a linear combination of principal strain
 components ($d\epsilon_v = d\epsilon_1 + 2 d\epsilon_3 = 0$).

Numerically, to impose the loading programme defined by [6] to the sample, sev-
eral 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:

$$
C_{\sigma}^{t} = \sigma_{1}^{t} - \frac{\sigma_{3}^{t}}{R} \quad \text{and} \quad C_{\epsilon}^{t} = \epsilon_{1}^{t} + 2R\epsilon_{3}^{t} \qquad [7]
$$

where $C_{\sigma}^{t}$ and $C_{\epsilon}^{t}$ are constants set at each step $t$ such that conditions defined in [6] are verified when the simulation ends; $\sigma_{1}^{t}$, $\sigma_{3}^{t}$, $\epsilon_{1}^{t}$ and $\epsilon_{3}^{t}$ 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><img src="image_9.png" /><figcaption>**Figure 8.** Principle of the loop to control linear combinations of principal stress or strain components</figcaption></figure>

In a first time we consider the stress direction characterised by $R = 1$ ($\alpha = 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\epsilon_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\epsilon_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 be visualized in the evolution of the axial $\sigma_1$ and radial $\sigma_3$ stresses; when $\sigma_3$ vanishes, $\sigma_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\epsilon_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\epsilon_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:

$$d\sigma_1 - \frac{d\sigma_3}{R} = 0 \quad \text{and} \quad d\epsilon_1 + 2R d\epsilon_3 = 0 \qquad [8]$$

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.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison of responses of sample E3 along a constant shear stress loading path control by $d\sigma_1 = d\sigma_3 < 0$ or by $dq = 0$ and $d\epsilon_v < 0$</figcaption></figure>

<table><caption>**Table 2.** Stress directions and belonging to cones of unstable stress directions ($\sigma_3 = 100$ kPa)</caption><thead><tr><th rowspan="2">α<br>(deg)</th><th rowspan="2">R</th><th colspan="2">∈ cone of unstable stress directions?</th></tr><tr><th>sample E1 (η = 0.82)</th><th>sample E3 (η = 0.46)</th></tr></thead><tbody><tr><td>200</td><td>1.94</td><td>no</td><td>no</td></tr><tr><td>220</td><td>0.843</td><td>no</td><td>yes</td></tr><tr><td>230</td><td>0.593</td><td>yes</td><td>yes (close to the limit of the cone)</td></tr><tr><td>250</td><td>0.257</td><td>no</td><td>no</td></tr></tbody></table>

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, $\sigma_3$, evo-
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²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.

α = 200 deg; R = 1.94

<figure><img src="image_11.png" /></figure>

α = 220 deg; R = 0.843

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Comparison of responses, to a perturbation in kinetic energy, computed with samples E1 and E3, for α = 200 and 220 deg</figcaption></figure>

**5. 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.

$$ \alpha = 230 \text{ deg}; R = 0.593 $$

<figure><img src="image_13.png" /></figure>

$$ \alpha = 250 \text{ deg}; R = 0.257 $$

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Comparison of responses, to a perturbation in kinetic energy, computed with samples E1 and E3, for $\alpha = 230$ and $250$ deg</figcaption></figure>

## 6. References

* Bardet J., "Numerical simulations of the incremental responses of idealized granular materials", *International Journal of Plasticity*, vol. 10, n° 8, p. 879-908, 1994.

* Calvetti F., Viggiani G., Tamagnini C., "A numerical investigation of the incremental behavior of granular soils", *Rivista Italiana di Geotecnica*, vol. 3, p. 11-29, 2003.

* Combe G., Roux J., "Strain versus stress in a model granular material: a devil's staircase", *Physical Review Letters*, vol. 85, n° 17, p. 3628-3631, 2000.

* Cundall P., Strack O., "A discrete numerical model for granular assemblies", *Geotechnique*, vol. 29, n° 1, p. 47-65, 1979.
Darve F., Laouafa F., “Instabilities in granular materials and application to landslides”, *Mech. Cohes. Frict. Mater.*, vol. 5, n° 8, p. 627-652, 2000.

Darve F., Servant G., Laouafa F., Khoa H., “Failure in geomaterials : continuous and discrete analyses”, *Computer Meth. Appl. Mech. Eng.*, vol. 193, n° 27-29, p. 3057-3085, 2004.

Donzé F., Magnier S., “Spherical discrete element code”, *Discrete Element Project Report 2*, Laboratory GEOTOP, Université du Québec à Montréal, 1997.

Gudehus G., “A comparison of some constitutive laws for soils under radially symmetric loading and unloading”, *Proc. 3rd Numer. Meth. in Geomechanics*, A. A. Balkema, Aachen, p. 1309-1323, 2-6 April, 1979.

Imposimato S., Nova R., “An investigation on the uniqueness of the incremental response of elastoplastic models for virgin sand”, *Mech. Cohes. Frict. Mater.*, vol. 3, p. 65-87, 1998.

Imposimato S., Nova R., “On the value of the second order work in homogeneous test on loose sand specimens”, in H.-B. M. et al. (ed.), *Bifurcation and Localization Theory in Geomechanics*, Swets and Zeitlinger, p. 209-216, 2001.

Nicot F., Darve F., Khoa H., “Bifurcation and second order-work in granular materials”, *Int. J. Numer. Anal. Meth. Geomech.*, vol. 31, p. 1007-1032, 2007.

Nova R., “Controllability of the incremental response of soil specimens subjected to arbitrary loading programmes”, *J. Mech. Behav. Mater.*, vol. 5, n° 2, p. 193-201, 1994.

Sibille L., Modélisations discrètes de la rupture dans les milieux granulaires, PhD thesis, Institut national polytechnique de Grenoble, Grenoble, 2006.

Sibille L., Donzé F., Nicot F., Chareyre B., Darve F., “From bifurcation to failure in a granular material, a DEM analysis”, *Acta Geotechnica*, vol. 3, n° 1, p. 15-24, 2008.

Sibille L., Nicot F., Donzé F., Darve F., “Material instability in granular assemblies from fundamentally different models”, *Int. J. Numer. Anal. Meth. Geomech.*, vol. 31, p. 457-481, 2007.