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

HAL
open science

# 3D bifurcation analysis in geomaterials. Investigation of the second order work criterion

Florent Prunier, Farid Laouafa, Félix Darve

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

## To cite this version:

Florent Prunier, Farid Laouafa, Félix Darve. 3D bifurcation analysis in geomaterials. Investigation of the second order work criterion. European Journal of Environmental and Civil Engineering, 2009, 13 (2), pp.135-147. 10.1080/19648189.2009.9693096. hal-01714887

HAL Id: hal-01714887

https://hal.science/hal-01714887

Submitted on 22 Feb 2018

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 3D bifurcation analysis in geomaterials

## Investigation of the second order work criterion

Florent Prunier* — Farid Laouafa** — Félix Darve*

* Laboratoire 3S-R, INPG-UJF.CNRS, BP 53 38041 Grenoble cedex 9 - France.
prunier@geo.hmg.inpg.fr

** INERIS, Unité Risques Naturels, BP 2, 60550 Verneuil en Halatte - France.
farid.laouafa@ineris.fr

**ABSTRACT.** In this paper, a study of failure in geomaterials is proposed through the analysis of the second order work criterion. The analysis is restricted to the material point scale. Then a phenomenological approach based on “macroscopic” constitutive models is adopted. In a first part, an analytical investigation of this criterion is proposed. General 3D equation of instability cones as well as of the 3D bifurcation domain limit are given for every incrementally piece-wise linear constitutive model. In the second part, these instability cones and bifurcation domain are displayed for constitutive models of Darve. Then a physical interpretation of these results is proposed. Noticeably, it is proved that, when the second order work vanishes along a given loading path, an extension of the notions of failure limit condition as well as of flow rule can be defined for proper conjugated variables. It is also shown, that the Lode’s angle (and not only the mean pressure associated with a deviatoric level) is also important according to bifurcation. Finally we conclude that this approach constitutes a good starting point when investigating bifurcation in geomaterials, and opens new horizons in experimental testing.

**RÉSUMÉ.** Dans cet article, une étude de la rupture dans les géomatériaux est proposée par l’analyse du critère du travail du second ordre. Cette analyse se limite à l’échelle du point matériel. Une approche phénoménologique basée sur des lois de comportement “macroscopic” est adoptée. Dans une première partie, ce critère est étudié analytiquement. L’équation générale à 3 dimensions dans l’espace des contraintes principales des cônes d’instabilités ainsi que celle de la limite du domaine de bifurcation sont données pour toute loi incrémentalement linéaire par morceaux. Dans la seconde partie, les cônes d’instabilités ainsi que le domaine de bifurcation pour les lois de comportement de Darve sont présentés dans cet espace. Puis une interprétation physique de ces résultats est proposée. Il est notemment prouvé, que lorsque le travail du second ordre s’annule le long d’un chemin de chargement donné, une extension des notions de la condition limite de rupture ainsi que celle de règle d’écoulement peuvent être définis pour un certain jeu de variables conjuguées. Il est aussi montré que l’angle de Lode (et
pas uniquement la pression hydrostatique associé au déviateur) est aussi important concernant
la bifurcation. Finalement, nous concluons que cette approche est bien appropriée à l'étude de
la bifurcation dans les géomatériaux, et ouvre de nouveaux horizons pour les expérimentateurs.

KEYWORDS: granular materials; bifurcation; second order work criterion; constitutive equations

MOTS-CLÉS : matériaux granulaires, bifurcation, travail du second ordre, loi de comportement
# 1. Introduction

The aim of this paper is to study failures in rate independent geomaterials. It is now widely recognized, that different failure modes can develop strictly inside the plasticity limit condition for such materials (soils, concrete, granular media) (Lade, 1992), (Nova, 1991), (Darve et al., 1987), (Darve et al., 2004b). One of the most important feature of these materials according to the failure, is their non associativeness. Normality rule is not fulfilled, and as a consequence, the constitutive tensor loses its symmetry (Hill, 1967), (Mandel, 1966), (Mróz, 1963), (Mróz, 1966). Because of this feature, strain localisation in shear bands (Desrues et al., 2004), (Vardoulakis et al., 1995), as well as diffuse mode of failure (Darve et al., 2004a),(Khoa et al., 2006), (Darve et al., 2007), can develop before reaching Mohr-Coulomb criterion.

Methods which allow to analyse such kind of failures, are generally based on loss of uniqueness of particular constitutive equation solutions. In that sense, criteria issue from these analyses are particular type of bifurcation criteria by opposition to stability criteria. Motions equations are not investigated, and as a consequence, stability of an equilibrium state (or even of a trajectory of a solid) is not studied. Nevertheless, loss of uniqueness of such equation solutions characterize the occurrence of a discontinuity in the response mode by keeping continuous loading parameters. This is the basic feature of a bifurcation.

Limit of plasticity and localization criterion of Rice (Rice, 1976), (Rudnicki et al., 1975), are the most common bifurcation criteria used by engineers. Nevertheless, as already written above, strain localization can appear before the plasticity limit condition, but in addition, diffuse failure modes which are neither described by the plasticity condition, nor by the Rice criterion can occur. The Hill's sufficient condition of stability seems describing accurately these different failure modes. Furthermore, Bigoni and Hueckel (Bigoni et al., 1991) have proved that this condition is reached before Rice condition and limit plasticity one.

This condition, stated by Hill (Hill, 1958), reads the following expression:

$$ \int_V \delta s_{ij} d\left(\frac{\partial u_j}{\partial x_i}\right) dV > 0 \quad [1] $$

with V the volume of the body at stage t, $s_{ij}$ the transpose of the first Piola-Kirchoff tensor, and $d(\frac{\partial u_j}{\partial x_i})$ the deformation velocity. In the following of this paper, the analysis is performed only at the material point scale. Moreover, small strain assumption is done and the geometrical effects are neglected. Thus the following condition is derived:

$$ d^2 W = d\sigma_{ij} d\epsilon_{ij} > 0 \quad \forall \|d\epsilon\| \neq 0 \quad [2] $$

This is the so called second order work criterion. $\sigma_{ij}$ and $\epsilon_{ij}$ denotes respectively, the components of the Cauchy stress tensor, and of the small strain tensor. It is worth noting that those quantities are linked by their constitutive relation. In the following of the paper, bifurcation in geomaterials is investigated through this criterion.

In a first part, some analytical developments and results about the second order
work are proposed for the class of incrementally piece-wise linear constitutive models. In fact, this class of models covers the majority of those proposed in the literature for rate independent materials. Secondly, a physical interpretation of these results is given, and numerical illustrations are presented with the constitutive models of Darve (Darve et al., 1995).

## 2. Some analytical results

In this section, an analytical analysis of relation [2] is proposed. To do so, following notations are used. $\underline{M}$ is the constitutive matrix which links $d\sigma$ to $d\varepsilon$ ($d\sigma = \underline{M}d\varepsilon$). $\underline{N}$ is the one which links $d\varepsilon$ to $d\sigma$. If $\underline{M}$ is invertible, $\underline{N} = \underline{M}^{-1}$. Finally, we define $\underline{du}$, $\underline{dv}$ the vectors of mixed linear combinations of $d\sigma$ and $d\varepsilon$ and $\underline{S}$ the constitutive matrix (built with the help of $\underline{M}$, or $\underline{N}$) which links $\underline{du}$ and $\underline{dv}$. This formulation allows the description of any loading path on an homogeneous sample. For example, if we want to describe an undrained triaxial compression test<fnref n="1" /> driven with the axial force, the following problem can be built:

$$ \begin{bmatrix} dq \\ d\varepsilon_v \end{bmatrix} = \underline{S} \begin{bmatrix} d\varepsilon_1 \\ d\sigma_3 \end{bmatrix} \qquad [3] $$

with $dq = d\sigma_1 - d\sigma_3$, and $d\varepsilon_v = (d\varepsilon_1 + 2d\varepsilon_3)$. In that way, $\underline{du} = {}^t[dq, d\varepsilon_v]$, and $\underline{dv} = {}^t[d\varepsilon_1, d\sigma_3]$ are conjugated variables according the second order work. By splitting matrices $\underline{M}$ and $\underline{N}$ in a symmetrical part $\underline{M}_s$ and $\underline{N}_s$ and a skew symmetrical part, the second order work takes the following expression :

$$ d^2 W = {}^t \underline{du} \cdot \underline{dv} = {}^t \underline{d\sigma} \cdot \underline{d\varepsilon} = {}^t \underline{d\sigma} \underline{N}_s \underline{d\sigma} = {}^t \underline{d\varepsilon} \underline{M}_s \underline{d\varepsilon} \qquad [4] $$

Thus, the sign of the second order work is strongly linked to the positiveness of $det(\underline{N}_s)$ or $det(\underline{M}_s)$. Positiveness of the quadratic form ${}^t \underline{d\sigma} \underline{N}_s \underline{d\sigma}$ gives a geometrical representation of the positiveness of $d^2 W$ in the stress rate space, while ${}^t \underline{d\varepsilon} \underline{M}_s \underline{d\varepsilon}$ in the strain rate space. Let us now analyse the condition :

$$ {}^t \underline{d\sigma} \underline{N}_s \underline{d\sigma} = 0 \qquad [5] $$

For any incrementally piece-wise linear constitutive relation, $\underline{N}$ can take the following form in a given tensorial zone<fnref n="2" /> and in principal axes:

$$ \underline{N} = \begin{bmatrix} \frac{1}{E_1} & -\frac{\nu_{12}}{E_2} & -\frac{\nu_{13}}{E_3} \\ -\frac{\nu_{12}}{E_1} & \frac{1}{E_2} & -\frac{\nu_{32}}{E_3} \\ -\frac{\nu_{13}}{E_1} & -\frac{\nu_{23}}{E_2} & \frac{1}{E_3} \end{bmatrix} \qquad [6] $$

<footnote n="1">1. assuming the fact that the pore water and each grain of the skeleton are incompressible, which is the case for the considered pressures</footnote>

<footnote n="2">2. A tensorial zone is a domain of the loading space in which the incremental constitutive relation is linear. A classical elasto-plastic model has two tensorial zones: one for plastic loading and the second one for elastic unloading.</footnote>
Hence, equation [5] reads :

$$ \frac{d\sigma_1^2}{E_1} + \frac{d\sigma_2^2}{E_2} + \frac{d\sigma_3^2}{E_3} - \left(\frac{\nu_{12}}{E_1} + \frac{\nu_{21}}{E_2}\right) d\sigma_1 d\sigma_2 - \left(\frac{\nu_{32}}{E_3} + \frac{\nu_{23}}{E_2}\right) d\sigma_3 d\sigma_2 - \left(\frac{\nu_{13}}{E_1} + \frac{\nu_{31}}{E_3}\right) d\sigma_1 d\sigma_3 = 0 \quad [7] $$

This equation, is the one of an elliptical cone in the stress rate space. Nevertheless, some degenerated form can exist depending on the positiveness of the $N_s$ eigenvalues. If we assume that all eigenvalues are positive at the virgin state, and are evolving continuously with the loading parameter, four cases have to be taken into account. Let us denote $(\lambda_1, \lambda_2, \lambda_3)$ the three eigenvalues of $N_s$, with the additional condition: $\lambda_1 < \lambda_2 < \lambda_3$. The four possibilities are the following :

1) $\lambda_1 > 0 : \det(N_s) > 0$. As a result, excepted the trivial solution ($d\sigma_1 = d\sigma_2 = d\sigma_3 = 0$), equation [7] is not fulfilled. The stress-strain state is “stable”.

2) $\lambda_1 = 0, \lambda_2 > 0 : \det(N_s) = 0$. Then, equation [7] admits a straight line (in the direction of the eigenvector of $\lambda_1$) as solution. Second order work vanishes according to this loading direction.

3) $\lambda_1 < 0, \lambda_2 > 0 : \det(N_s) < 0$. Equation [7] admits an elliptical cone as solution. Second order work vanishes for loading directions upon the cone, is already negative for loading directions included inside the cone and positive elsewhere.

4) $\lambda_1 < 0, \lambda_2 = 0, \lambda_3 > 0 : \det(N_s) = 0$. Equation [7] admits two secant planes as solution. Second order work vanishes on these planes, is positive in two of the subdivisions of the stress rate space, and negative in the two others. However, this solution has never been reached with constitutive models used by co-authors.

The **Figure 1** illustrates the previous discussion. Thus, the directional feature of this criterion observed by Darve and Laouafa (Darve et al., 2000), (Laouafa et al., 2002) is proved here. This means that, as $\det(N_s) \le 0$ or conversely $\det(M_s) \le 0$, it exists some loading directions for which bifurcation occurs. As a remark, this directional feature is stronger for non elastic behaviour. In fact, real solutions are the ones (or part of ones) of conditions (2) to (4) which are fulfilled in the proper tensorial zone. If we call the limit of the bifurcation domain, the set of stress-strain states for which a unique unstable loading direction exists, this limit can be describe by the following condition :

$$ \min_{i=1,...,n} \left( \det \left( \underline{N_s} \right)_i \right) = 0 \quad \text{with} \quad \underline{u}_i \subset \underline{Z}_i \quad [8] $$

with $n$, the number of tensorial zones (of the constitutive model), $u_i$ the eigenvector
<figure><img src="image_4.png" /><figcaption>Figure 1. Solutions of equation : $\lambda_1 X^2 + \lambda_2 Y^2 + \lambda_3 Z^2 = 0$.</figcaption></figure>

corresponding to the vanishing eigenvalue, and $Z_i$ the tensorial zone considered.
To conclude with this section, we will prove that :

$$ \det(\underline{M_s}) = 0 \Leftrightarrow \det(\underline{N_s}) = 0 \qquad [9] $$

as equation [4] impose it. In fact, we have stated that :

$$ \det(\underline{M_s}) = \frac{\det(\underline{N_s})}{(\det(\underline{N}))^2}, \qquad [10] $$

which is equivalent to :

$$ \det(\underline{N} \cdot \underline{M_s}) = \det(\underline{N}^{-1} \cdot \underline{N_s}) \qquad [11] $$

Nicot (Nicot, 2008) proposed the following proof of relation [11]. Considering the left hand side of equation [11], we have:

$$ \underline{N} \cdot \underline{M_s} = \frac{1}{2} (\underline{I} + \underline{N} \cdot \underline{t} \underline{N}^{-1}) \qquad [12] $$
and by considering the right hand side of equation [11], it follows:

$$ \underline{N}^{-1} \cdot \underline{N}_s = \frac{1}{2} (\underline{I} + \underline{N}^{-1} \cdot \underline{t} \underline{N}) \quad [13] $$

thus, as $det(\underline{A}) = det(\underline{t} \underline{A})$ for any matrix $\underline{A}$, it follows :

$$ det(\underline{N} \cdot \underline{M}_s) = det(\frac{1}{2}(\underline{I} + \underline{N}^{-1} \cdot \underline{t} \underline{N})) = det(\underline{N}^{-1} \cdot \underline{N}_s) \quad [14] $$

This result means that bifurcation analysis can be performed without any restriction on the use of the constitutive relation, whereas it is not the case for the plasticity limit condition. $det(\underline{M}) = 0$ is not equivalent to $det(\underline{N}) = 0$. In the case of associated material, relation [10] is still valid, but condition [9] not any more, because in this case the bifurcation limit coincides with the plasticity limit condition.

### 3. Physical interpretation

In this section, numerical results, as well as a physical interpretation of analytical developments presented in the previous part are given.

Numerical results, are displayed with the constitutive models of Darve (Darve et al., 1995). Without going into details of these models, we just recall that they are not based on the classical concepts of the elasto-plasticity. Decomposition of the strain in an elastic and plastic part is not assumed, and no plastic potentials are defined. According to the first model, the non linear relation which links strain rate to stress rate is directly described by an incrementally non-linear relation. The second model is a simplification of the first one, and becomes incrementally piece-wise linear with eight tensorial zones. In principal axes, these models are written as follows :

$$ \begin{bmatrix} d\epsilon_1 \\ d\epsilon_2 \\ d\epsilon_3 \end{bmatrix} = \frac{1}{2} [N^+ + N^-] \begin{bmatrix} d\sigma_1 \\ d\sigma_2 \\ d\sigma_3 \end{bmatrix} + \frac{1}{2 \|\underline{d\sigma}\|} [N^+ - N^-] \begin{bmatrix} d\sigma_1^2 \\ d\sigma_2^2 \\ d\sigma_3^2 \end{bmatrix} \quad [15] $$

and

$$ \begin{bmatrix} d\epsilon_1 \\ d\epsilon_2 \\ d\epsilon_3 \end{bmatrix} = \frac{1}{2} [N^+ + N^-] \begin{bmatrix} d\sigma_1 \\ d\sigma_2 \\ d\sigma_3 \end{bmatrix} + \frac{1}{2} [N^+ - N^-] \begin{bmatrix} |d\sigma_1| \\ |d\sigma_2| \\ |d\sigma_3| \end{bmatrix} \quad [16] $$

with

$$ N^{\pm} = \begin{bmatrix} \frac{1}{E_1^{\pm}} & -\frac{\nu_{21}^{\pm}}{E_2^{\pm}} & -\frac{\nu_{31}^{\pm}}{E_3^{\pm}} \\ -\frac{\nu_{12}^{\pm}}{E_1^{\pm}} & \frac{1}{E_2^{\pm}} & -\frac{\nu_{32}^{\pm}}{E_3^{\pm}} \\ -\frac{\nu_{13}^{\pm}}{E_1^{\pm}} & -\frac{\nu_{23}^{\pm}}{E_2^{\pm}} & \frac{1}{E_3^{\pm}} \end{bmatrix} \quad [17] $$
Limit of the bifurcation domain given by equation [8], are displayed on **Figure 2** for both models. As a remark, the incrementally non-linear model can be seen as incrementally piece-wise linear with an infinity of tensorial zones. Thus, with a numerical effort, an approximation of the bifurcation limit has been displayed. Then, **Figure 3**

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Limit of the bifurcation domain plotted in the 3D stress space for constitutive models of Darve.</figcaption></figure>

presents some instability cones for stress-strain states situated beyond the bifurcation limit according to both models. These stress-strain states have been obtained along the simulation of a drained triaxial test.

Now a physical interpretation is given with the formalism proposed by Nova (Nova, 1994), (Nova, 2004). In fact, what happens when $d^2W$ vanishes along a given loading path is presented. As the positiveness of $d^2W$ depends on the loading rate, stress proportional loading paths are introduced as follows :

$$ \left\{ \begin{array}{l} d\sigma_1 = \text{constant} \quad \text{constant} \in \mathbb{R} \\ d\sigma_1 + R d\sigma_3 = 0 \quad R \in \mathbb{R} - \{0\} \\ d\sigma_2 - R' d\sigma_3 = 0 \quad R' \in \mathbb{R} \end{array} \right. \qquad [18] $$

These paths allow to prospect every loading directions from a given stress-strain state. Thus, we are able to prospect some particular stress paths where the second order work vanishes. The response of such paths is displayed on **Figure 4**. In fact, a previous drained triaxial compression is performed until a stress-strain state situated outside of
<figure><img src="image_6.png" /><figcaption>**Figure 3.** 3D cones of unstable stress direction for a dense sand of Hostun. Figure 3.(a) presents the cones obtained with the octo-linear model. The planes represent the limit between the 8 tensorial zones, the meshes correspond to the analytical solution of Equation [7] and the point clouds the solution obtained with the numerical method. Figure 3.(b) shows results obtained with the non-linear model using the numerical method. $p_0$ is the initial confining pressure, $q = \sigma_1 - \sigma_3$, $\eta = \frac{q}{p}$</figcaption></figure>

the bifurcation domain, then one of the proportional loading (going inside the bifurcation domain) path is followed. To simulate these loading paths, some of the variables have to be controlled. One of the possibilities is as follows :

$$ \begin{bmatrix} d\varepsilon_1 - \frac{d\varepsilon_3}{R} - \frac{R'}{R} d\varepsilon_2 \\ d\sigma_1 + R d\sigma_3 \\ d\sigma_2 - R' d\sigma_3 \end{bmatrix} = \underline{S} \begin{bmatrix} d\sigma_1 \\ \frac{d\varepsilon_3}{R} + \frac{R'}{R} d\varepsilon_2 \\ d\varepsilon_2 \end{bmatrix} \quad [19] $$

For proportional loading paths which show a peak in the $(\varepsilon_1 - \frac{\varepsilon_3}{R} - \frac{R'}{R} \varepsilon_2)$ versus $(\sigma_1)$ plane, following observations can be derived. At this peak, due to the statical conditions imposed by equation [18], we obtain :

$$ \underline{S} \begin{bmatrix} d\sigma_1 \\ \frac{d\varepsilon_3}{R} + \frac{R'}{R} d\varepsilon_2 \\ d\varepsilon_2 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} \quad [20] $$

Equation [20] has no trivial solutions if $det(\underline{S}) = 0$. Then this equation constitutes a generalized failure rule, and the eigenvector corresponding to the vanishing eigenvalue constitutes a generalized flow rule. Furthermore, at this state, $d^2W$ vanishes also. Then it can be concluded that this stress-strain state is situated inside the bifurcation
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Response of proportional stress loading paths obtained with the octo-linear model after a triaxial compression at $p_o = 600$ kPa and $q = 640$ kPa. $R \in [-3.5; -2.5]$ and $R' \in [0.8; 1.2]$. 4 (a), displays the loading paths. On 4 (b), the response path related to the triaxial compression is not displayed. $p = \frac{1}{3} \text{tr}(\underline{\underline{\sigma}})$, $q = \sqrt{\frac{3}{2}} \|\text{dev}(\underline{\underline{\sigma}})\|$, $(\text{dev}(\underline{\underline{\sigma}}) = \underline{\underline{\sigma}} - pI)$.</figcaption></figure>

domain, and that the loading direction belongs to one of the cone presented above (or to the only one unstable loading direction, if the stress-strain state is situated strictly on the bifurcation domain limit).

With such analysis, concept of controllability introduced by Nova (Nova, 1991), (Nova, 1994) is easily understandable. Considering previous loading path which show a peak in $(\varepsilon_1 - \frac{\varepsilon_3}{R} - \frac{R'}{R} \varepsilon_2)$, if $d\sigma_1$ is controlled, this peak can be passed continuously, whereas if $(d\varepsilon_1 - \frac{d\varepsilon_3}{R} - \frac{R'}{R} d\varepsilon_2)$ is controlled, effective collapse of the sample occurs at this peak. More classically, considering a drained triaxial test on dense sample, when axial strain is controlled going through the peak of $q$ is possible, while when axial stress is controlled, effective collapse is noticed at this peak. Furthermore, at such bifurcation point (state for which $d^2W$ vanishes on the considered loading path), Sibille and co-authors (Sibille et al., 2007) have proved with a discrete element model that a “small” mechanical perturbation on the sample induces the effective collapse. Thus, as the second order work vanishes on a given loading path, an unstable state is reached. Effective collapse occurs if particular controlled parameters are used, or if a “small” perturbation is imposed.

To conclude with this section, we present others instability cones on **Figure 5** and **Figure 6** for both models of Darve. In this case, loading paths considered are not any more drained triaxial compressions. In fact, a set of loading paths going towards some given directions located in a stress deviatoric plane are considered. Then, at an arbitrary deviatoric level (situated inside the bifurcation domain), 3D instability cones are plotted. Because, in these conditions, the mean pressure and the second stress invariant are fixed, the influence of the Lode's angle on instability cones can be ob-
served. These figures are 3D, but the point of view is turned along the hydrostatic
axis. That's why the representation seems to be done in the deviatoric plane. Ac-

<figure><img src="image_8.png" /><figcaption>**Figure 5.** 3D cones of unstable stress direction for paths in the deviatoric stress plane with the octo-linear model: $p_o = 600$ kPa, $||dev(\underline{\sigma})|| = 350$ kPa or $q = 428.6$ kPa, $(q = \sqrt{\frac{3}{2}}||dev(\underline{\sigma})||)$. Continuous line represents the bifurcation limit</figcaption></figure>

cording to the Lode's angle, it can be verified that results are 120° periodic. Then, the angle 60° modulo 120° constitutes an axis of symmetry. What is quite interesting on these figures, is to observe that angle 30° modulo 120° seems to show a discontinuity according to instability cones. This angle of 30° corresponds to plane stress conditions. It can be observed for both models that opening of cones are noticeably larger for Lodes angle varying from 30° and 60° than for ones varying from 0° and
<figure><img src="image_9.png" /><figcaption>**Figure 6.** 3D cones of unstable stress direction for paths in the deviatoric stress plane with the non-linear model: $p_o = 600$ kPa, $||dev(\underline{\sigma})|| = 400$ kPa or $q = 489.9$ kPa, ($q = \sqrt{\frac{3}{2}||dev(\underline{\sigma})||}$). Continuous line represents the bifurcation limit</figcaption></figure>

30°. Thus, the three stress invariants plays an important role according to bifurcation, and not only the mean pressure and the deviatoric level.
4. Concluding remarks

In this paper, a bifurcation analysis based on the second order work criterion has been performed to analyse failure in geomaterials. The analysis was restricted to the material point scale. In introduction it has been reminded that this criterion constitutes a lower bound according to the plasticity limit and the localization criterion of Rice. In some particular cases, like the modelling of landslides under very gentle slope, this criterion is necessary in predicting failure. For example, at Petacciato in Italy, a land- slide occurs under a mean slope of 6° without any localisation pattern. The modelling using this criterion has proved to be satisfactory (LESSLOSS-Report, 2004), (Lignon *et al.*, 2008).

In the first part, analytical results valid for any incrementally piece-wise linear
model, have been stated. Existence of instability cones has been proved, and their
equation given in the 3D space of stress principal axes. Then the limit of the bifur-
cation domain has been derived and its equation, which can be solved numerically,
has been established. Finally, it has also been proved that the equation of the bifur-
cation domain limit is not changed when the constitutive tensor *M* is considered or
conversely its inverse *N* in a dual formalism.

In the second part, 3D representation of instability cones as well as of the bifur-
cation domain limit have been given for Darve constitutive models. Then a physical
interpretation of this criterion has been proposed. Noticeably, it has been proved that,
as the second order work vanishes along a given loading path, an extension of the plas-
tic limit criterion and of the flow rule can be defined for proper conjugated variables.
Eventually, the influence of the Lode's angle on the opening of instability cones has
been shown. Thus the third stress invariant plays also an important role when investi-
gating bifurcation in geomaterials.

According to these results, second order work criterion constitutes a good tool to
investigate bifurcation in geomaterials. Noticeably, it opens new perspectives for ex-
perimental testing. All proportional stress loading paths presented in this paper should
be reproducible experimentally. Further more proportional strain loading paths can
also be defined in the same way. Daouadji and co-authors (Darve et al., 2007) have
already investigated some of these paths experimentally. One of the exciting challenge
would be proposing an experimental set-up allowing the description of the bifurcation
domain limit.

5. References

Bigoni D., Hueckel T., “Uniqueness and localization-I. Associative and non-associative elasto-plasticity”, *International Journal of Solids and Structures*, vol. 28(2), p. 197-213, 1991.

Darve F., Chau B., “Constitutive instabilities in incrementally non-linear modelling”, *Constitutive laws for Engineering Materials*, C. S. Desai, p. 301-310, 1987.

Darve F., Flavigny E., Meghachou M., “Yield surfaces and principle of superposition revisited by incrementally non-linear constitutive relations”, *International Journal of Plasticity*, vol. 11(8), p. 927-948, 1995.
Darve F., Laouafa F., “Instabilities in granular materials and application to landslides”, *Mechanics of cohesive-frictional materials*, vol. 5, p. 627-652, 2000.

Darve F., Servant G., Laouafa F., Khoa H.D.V., “Failure in geomaterials : continuous and discrete analyses”, *Computer Methods in Applied Mechanics and Engineering*, vol. 193, p. 3057-3085, 2004a.

Darve F., Sibille L., Daouadji A., Nicot F., “Bifurcations in granular media: macro- and micro-mechanics approaches”, *Comptes Rendus Académie des Sciences Mécanique*, vol. 335, p. 496-515, 2007.

Darve F., Vardoulakis I. (eds), *Degradations and instabilities in geomaterials*, vol. 461 of CISM courses, SPRINGER, 2004b.

Desrues J., Viggiani J., “Strain localization in sand: an overview of the experimental results obtained in Grenoble using stereophotogrammetry”, *International Journal for Numerical and Analytical Methods in Geomechanics*, vol. 28(4), p. 279-321, 2004.

Hill R., “A general theory of uniqueness and stability in elasto-plastic solids”, *Journal of the Mechanics and Physics of Solids*, vol. 6, p. 236-249, 1958.

Hill R., “Eigenmodal deformations in elastic/plastic continua”, *Journal of the Mechanics and Physics of Solids*, vol. 15, p. 371-386, 1967.

Khoa H.D.V., Georgopoulos I.O., Darve F., Laouafa F., “Diffuse failure in geomaterials: Experiments and modelling”, *Computer and Geotechnics*, vol. 33, p. 1-14, 2006.

Lade P.V., “Static Instability and Liquefaction of Loose Fine Sandy Slopes”, *Journal of Geotechnical Engineering*, vol. 118(1), p. 51-71, 1992.

Laouafa F., Darve F., “Modelling of slope failure by a material instability mechanism”, *Computer and Geotechnics*, vol. 29, p. 301-325, 2002.

LESSLOSS-Report, The Petacciato landslide: geological and geotechnical data, Technical report, SGI: Studio Geotechnico Italiano, Campobasso-Italy, 2004.

Lignon S., Laouafa F., Prunier F., Darve F., Khoa H.D.V., “Hydro-mechanical modelling of landslides with a material instability criterion”, *Geotechnique*, 2008. In press.

Mandel J., “Conditions de stabilité et postulat de Drücker”, in J. Kravtchenko, P.M. Sirieys (eds), *Rheology and Soils Mechanics*, Berlin Springer, p. 58-68, 1966.

Mróz Z., “Non-associated flow laws in plasticity”, *Journal de Mécanique*, vol. 2, p. 21-42, 1963.

Mróz Z., “On forms of constitutive laws for elastic-plastic solids”, *Archive of Mechanics Stosovanej*, vol. 18, p. 1-34, 1966.

Nicot F.2008. private communication.

Nova R., “A note on sand liquefaction and soil stability”, 3rd *International Conference on Constitutive Laws for Engineering Materials: Theory and Applications*, Tuscon, USA, p. 153-156, 1991.

Nova R., “Controllability of the incremental response of soil specimens subjected to arbitrary loading programs”, *Journal of the Mechanical behavior of Materials*, vol. 5(2), p. 193-201, 1994.

Nova R., “Controllability of geotechnical testing”, *Failure, Degradation and Instabilities in Geomaterials. Revue française de génie civil*, vol. 8(5-6), p. 613-634, 2004.

Rice J.R., “The localization of plastic deformation”, in W.T. Koiter (ed.), *Theoretical and Applied Mechanics*, p. 207-220, 1976. 14th IUTAM Congress Amsterdam.
Rudnicki J.W., Rice J.R., “Conditions for the localization of deformations in pressure-sensitive dilatant materials.”, *Journal of the Mechanics and Physics of Solids*, vol. 23, p. 371-394, 1975.

Sibille L., Nicot F., Donzé F.V., Darve F., “Material instability in granular assemblies from fundamentally different models”, *International Journal for Numerical and Analytical Methods Geomechanics*, vol. 31(3), p. 457-482, 2007.

Vardoulakis I., Sulem J. (eds), *Bifurcation Analysis in Geomechanics*, Blackie Academic and Professional, 1995.