# Diffuse failure in geomaterials: Experiments, theory and modelling

A. Daouadji, F. Darve, H. Al Gali, Pierre-Yves Hicher, F. Laouafa, S. Lignon, François Nicot, R. Nova, M. Pinheiro, Florent Prunier, et al.

► To cite this version:

A. Daouadji, F. Darve, H. Al Gali, Pierre-Yves Hicher, F. Laouafa, et al. Diffuse failure in geomaterials: Experiments, theory and modelling. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35 (16), pp.1731-1773. 10.1002/nag.975 . hal-01802801

HAL Id: hal-01802801

https://hal.science/hal-01802801

Submitted on 8 Jan 2020

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
# Diffuse failure in geomaterials: Experiments, theory and modelling

A. Daouadji<fnref n="1" />,*,†, F. Darve<fnref n="2" />, H. Al Gali<fnref n="1" />, P. Y. Hicher<fnref n="3" />, F. Laouafa<fnref n="4" />, S. Lignon<fnref n="5" />,
F. Nicot<fnref n="6" />, R. Nova<fnref n="7" />, M. Pinheiro<fnref n="8" />, F. Prunier<fnref n="2" />, L. Sibille<fnref n="3" /> and R. Wan<fnref n="8" />

## SUMMARY

This paper presents a synthesis of the works performed by various teams from France, Italy and Canada around the question of second-order work criterion. Because of the non-associative character of geomaterials plastic strains, it is now recognized that a whole bifurcation domain exists in the stress space with various possible modes of failure. In a first part these failure modes are observed in lab experimental tests and in discrete element modelling. Then a theoretical study of second-order work allows to establish a link with the kinetic energy, giving a basis to explain the transition from a prefailure (quasi)static regime to a postfailure dynamic regime. Eventually the main features of geomaterials failure are obtained by applying second-order work criterion to five different constitutive rate-independent models—three being phenomenological and two micromechanical. As a whole this paper tries to gather together all the elements for a proper understanding and use of second-order work criterion in geomechanics.

KEY WORDS: failure; diffuse failure; second-order work; bifurcation; nonassociative elastoplasticity; kinetic energy; instability; loss of uniqueness

## 1. INTRODUCTION

The precise definition of failure in materials and its subsequent analysis is an elusive, but yet a very important question in solid mechanics and engineering. Much effort has been expended towards establishing various failure criteria for more than a century now. From a physical point of view, failure appears as a phenomenon associated with the existence of some limit stress states which cannot be passed by the material for any possible monotonous loading path. At such a limit stress state if an additional load is imposed the state of the material changes suddenly with the occurrence of large deformations, cracks, fragmentation,... Indeed it becomes impossible to apply this additional load. This brutal change is called roughly ‘failure’.

<footnote n="1">¹Laboratoire de Physique et Mecanique des Materiaux, Universite Paul Verlaine-Metz, CNRS, Metz, France</footnote>

<footnote n="2">²Laboratoire Sols Solides Structures, UJF-INPG-CNRS, Grenoble, France</footnote>

<footnote n="3">³Institut de Recherche en Genie Civil et Mecanique, Universite de Nantes, UMR CNRS 6183,
Ecole Centrale Nantes, France</footnote>

<footnote n="4">⁴INERIS, Unité Risques Naturels, Ouvrages and Stockages/Direction des Risques du Sol et du Sous-Sol,
Verneuil en Halatte, France</footnote>

<footnote n="5">⁵Groupe de Recherche en Sciences Pour l'Ingénieur, Reims, France</footnote>

<footnote n="6">⁶Cemagref, Geomechanics Group-ETNA, Grenoble, France</footnote>

<footnote n="7">⁷Department of Structural Engineering, Politecnico di Milano, Milan, Italy</footnote>

<footnote n="8">⁸Department of Civil Engineering, University of Calgary, Calgary, Canada</footnote>

<footnote>*Correspondence to: A. Daouadji, Laboratoire de Physique et Mecanique des Materiaux, Universite Paul Verlaine-Metz, CNRS, Metz, France.</footnote>

<footnote>†E-mail: daouadji@univ-metz.fr</footnote>
Based on experimental observations, we can find two broad classes of failure modes that arise due to instabilities that can be either geometric as in column buckling, or material as in constitutive behaviour, respectively. Within the confines of material instability, two major criteria have emerged. The first criterion refers to the vanishing of the determinant of the constitutive tensor which signals homogeneous failure at the plastic limit condition, whereas the second one involves the vanishing of the determinant of the acoustic tensor [1] which coincides with the emergence of plastic strain localization.

For materials with an associative flow rule, as it is generally assumed for metals, the symmetry of the elasto-plastic tensor leads to the compelling fact that the two aforementioned criteria coincide. However, geomaterials are known to be non-associated in character, which leads to the loss of symmetry in the elasto-plastic tensor. Thus, the localization criterion can be met before the plastic limit criterion as has been shown consistently in laboratory experiments, particularly for dense sands and over-consolidated clays [2]. However, undrained triaxial tests on very loose sands, performed in load controlled mode, display a type of failure at peak deviatoric stress that is not described by any of the previously introduced criteria. Because the related failure mode does not display any localization, it has been coined as ‘diffuse failure’ [3] to distinguish it from the localized one. As such, the purpose of this paper is to precisely examine this particular failure mode and the material collapse that ensues from various experimental, analytical, constitutive and numerical frameworks and establish the state-of-the-art on the topic as it is seen by a network of researchers working together since 20 years.

The material instability leads to the loss of uniqueness of the solution of the underlying governing equations, hence to a bifurcation problem. Therefore, it is natural to advocate bifurcation theory as a general framework to analyse all the various kinds of failure. Bifurcation in a mechanical system occurs whenever the state of the system changes suddenly following one of possible multiple branches, which can be either stable or unstable, under continuous variations of state variables. For instance, under properly chosen loading conditions, failure can manifest itself as a sudden transition from a static regime to a dynamic one with an exponential growth of strains. This will be illustrated experimentally as well as numerically in Section 2 with the connection of the phenomenon to diffuse failure. Furthermore, as a property of bifurcated states, the attending failure will depend on perturbations and imperfections in the system as numerically shown in Section 2.2. Because of this dependence on small perturbations, failure can also be viewed as an instability phenomenon in the basic Lyapunov sense [4]. Indeed, for a given bifurcation state on the stress–strain curve, any additional and vanishingly small loading will lead to infinitely large responses. Moreover, the response paths are not unique because they depend on imperfections and not only on state variables characterizing material behaviour.

The failure of materials is certainly complex, but is inarguably related to material instability and a loss of uniqueness problem, although the converse may not be true. The notion of controllability [5] and sustainability of equilibrium states [6] are two concepts that are built around the above-mentioned question to describe failure. Essentially, the former refers to the impossibility to carry out a certain test under a given loading programme according to some control parameters, whereas the latter alludes to a certain stress–strain state for which it is not possible to apply any variations in control parameters without causing collapse. At any rate, these two inspiring concepts are very similar and there is a need to unify them, which is done in Section 3.

It was thus far back mentioned that the peak deviatoric stress in undrained triaxial tests on loose sands could not be described by any ‘classical’ plastic criterion. However, it appears [7] that the nullity of the second-order work criterion (as proposed by authors of [8–13]) is able to describe any kind of quasi-static constitutive instability by taking into account that flutter instabilities are dynamic. For classic elasto-plastic relations, the lowest stress level for which this criterion is fulfilled corresponds to the vanishing of the determinant of the symmetric part of the constitutive matrix related to loading conditions. Thus, this criterion is of particular significance for geomaterials which are non-associated and involve a non-symmetric constitutive matrix. An extensive study of this criterion is proposed in Section 4 of this paper by considering two different approaches. The first approach involves an analytical and numerical analysis of the second-order
work criterion through phenomenological elasto-plastic constitutive relations, whereas the second one considers micro-mechanical models. Both of these approaches allow a proper interpretation of the experimental results presented in Section 2.1 and of the numerical results obtained by direct numerical simulations with a discrete element method given in Section 2.2. The existence of a large bifurcation domain in the stress space and of cones of unstable loading directions for proper control parameters is established from experimental, theoretical and numerical view points. The equation of the bifurcation domain's boundary and of the so-called instability cones is derived in the case of incrementally piecewise constitutive relations. Much emphasis is made on the fact that an essential feature of diffuse failure is the outburst of kinetic energy that ensues with exponentially growing strains and decreasing stresses. The three necessary and sufficient conditions for an effective collapse are thus introduced: (1) the stress state must lie inside the bifurcation domain, (2) the loading directions must be within the instability cone, and (3) the proper control parameters must be set in place. The notions of limit stress states and of flow rules are generalized by considering mixed (i.e. stress-strain) loading conditions. The fact that these mixed loadings are a necessary ingredient of diffuse failure is demonstrated.

## 2. SALIENT ASPECTS OF SOIL BEHAVIOUR AND FAILURE

To base the paper on clear grounds, some failure modes which do not obey neither to a plastic limit condition nor to a localization criterion are first described in this section. In Section 2.1, some typical experimental tests are considered, whereas in Section 2.2 some direct numerical simulations of failures by a discrete element model are presented.

### 2.1. *Experimental evidence of diffuse failure in sand testing*

The experimental programme considered throughout this section consists of two types of tests carried out on the same material: namely, conventional undrained (CU) triaxial compression tests, and constant deviatoric stress drained (CSD) tests. It is recalled that in the CU test series, the sample was isotropically consolidated to the desired confining pressure followed by shearing at constant confining pressure under undrained conditions. On the other hand, in the CSD test series, the specimen follows three stages of loading: (1) an initial isotropic consolidation to the same confining pressure as in the CU test, (2) a drained compression loading at constant confining pressure, and (3) a drained isotropic stress unloading at constant deviatoric stress. In steps (2) and (3), the drainage valve was left open while volume changes, axial displacements and the axial load were continuously measured.

The CU test consists of loading the sample with increasing total stresses along a single monotonic strain- or stress-controlled loading path. In the CSD test, even if the mean effective stress decreases along a path for which the deviatoric stress is kept constant, the stress ratio $\eta = q/p'$ ($q = \sigma_1 - \sigma_3$ and $p' = (\sigma'_1 + 2\sigma'_3)/3$) increases so that also this is actually a loading path. As such, a complex monotonic stress-controlled path is applied to the specimen [14–17].

Generally speaking, as far as granular materials are involved, it is found that loose samples are very prone to diffuse failure (or collapse). Therefore, in the experimental exploration of diffuse failure, it is evident to choose the highest void ratio value as possible for testing. Moreover, a slenderness (height to diameter of specimen) ratio of 1.5 combined with enlarged frictionless top and bottom platens with free ends were used. These are to prevent as much as possible the localization of deformation and maintain a homogeneous deformation.

One of the basic requisites for reproducing a very loose sample in the laboratory is that the material must be uniformly sized. With this property in mind, we chose Hostun sand, a type of sand that is widely tested in France and which contains 99.2% of silica. The main characteristics are summarized in Table I. To perform all series of tests, two Digital Pressure–Volume Controllers (DPVC) were used, one for the cell pressure and the other for the pore water. For the later, a DPVC can impose a volume change and measure the pressure variation or conversely, it can impose a pore pressure and measure the corresponding volume variation. A submersible load cell is used to
<table><caption>Table I. Basic properties of Hostun S28 sand.</caption><thead><tr><th>Type</th><th>Mean size d<sub>50</sub> (mm)</th><th>Cu</th><th>Specific density</th><th>e<sub>max</sub></th><th>e<sub>min</sub></th><th>Ultimate state friction angle (°)</th></tr></thead><tbody><tr><td>Quartz</td><td>0.33</td><td>2</td><td>2.65</td><td>1</td><td>0.656</td><td>32</td></tr></tbody></table>

measure the force and LVDT are used for the axial displacement. Details concerning the sample preparation are given in [17].

### 2.1.1. Conventional Undrained Tests (CU).
Two undrained tests were conducted after an initial isotropic consolidation of $p'_{0}=300$ kPa for a targeted void ratio. The first test is performed in a strain-controlled mode, whereas the second one in a load-controlled mode. This distinction is essential to show the influence of the control parameters, namely the loading method and the isochoric condition, on the macroscopic mechanical behaviour of loose sands.

Figure 1 shows the results of the two undrained tests revealing approximately the same initial stiffness in the $q-\varepsilon_1$ plot. This indicates the reproducibility of both the sample preparation procedure and the test itself. Actually, a more comprehensive study of the repeatability of the tests was undertaken, but is not detailed herein. The two tests give the same peak deviatoric stress at approximately the same axial strain, i.e. around 1%, followed by a decrease in the deviatoric stress. While the strain-controlled test is controllable until very large strains are achieved (more than 15%), in the other case the decrease of the deviatoric stress is sudden in the $q-p'$ plane and accompanied by a jump of axial strain in the $\varepsilon_1-t/t_f$ plane, where the real time $t$ is normalized to the final (end) time $t_f$ of the test. This normalized time is introduced to compare the two tests which have a different duration.

Finally, a plateau is reached in the $q-\varepsilon_1$ stress-strain curve at low stress levels that indicate a partial liquefaction, which can also be seen in the $q-p'$ plot. In this latter plot, the stress ratio reaches a value of 1.32 at the end of the strain-controlled test. We note that this stress ratio corresponds to a mobilized friction angle of 32°, indicating that the material has reached the Mohr-Coulomb failure condition [17].

The load-controlled test differs from the strain-controlled one by the absence of controllability beyond the peak, and as such, neither plateau, nor the plastic limit condition is reached. The stress state corresponding to the peak is in fact a failure state which is well below the Mohr-Coulomb criterion. This distinction is essential for the understanding of the influence of the control parameters, and not the excess pore pressure for triggering collapse. This important fact has been discussed at length in [7, 14, 18–20]. The next section touches further on this point using the CSD test as back drop.

### 2.1.2. Constant deviatoric stress drained tests (CSD).
This type of test consists of loading the specimen to a prescribed deviatoric stress ratio along a drained compression triaxial path, followed by ‘unloading’ the specimen at decreasing mean effective stress under constant deviatoric stress. This later stress path can be regarded as describing the effective stress change in a slope or an embankment undergoing a slow increase in pore water pressure at constant total stress.

In order to check whether instability would occur under drained conditions and how the controlled variables would affect the onset of collapse, two CSD tests, CSD1 and CSD2 (Figure 2), were conducted under the same initial effective pressures of 300 kPa so that a comparison could also be made with the previously described CU tests.

Two methods for achieving the CSD tests after shearing the sample up to a given stress ratio are available. The first method (corresponding to CSD1 test), which is classic, refers to gently increasing the pore pressure while the total stresses $\sigma_1$ and $\sigma_3$ are kept constant. In this work, we propose also a second method (corresponding to CSD2 test) that consists of decreasing total stresses under the constraint that pore pressure remains constant to satisfy drained testing conditions. In fact, the specimen tends to dilate with a concomitant decrease in pore pressure. To keep pore water pressure constant, water has to be injected into the specimen. The volumetric changes that occur during this process were also recorded. The CSD test is thus achieved.
<figure><img src="image_2.png" /><figcaption>Figure 1. Influence of the control parameters on the mechanical behaviour of Hostun loose sand under undrained triaxial compression tests ($p'_0=300$ kPa). No sudden collapse is observed for strain-controlled test whereas it occurs at the peak of the load-controlled test when strain rate increases at collapse with decrease in $q$.</figcaption></figure>

The results of a typical test are depicted in Figures 2 and 3. For the CSD1 test, the specimen was first subjected to a drained triaxial compression test from point 0 to A at a value of 119 kPa, which corresponds to a stress ratio $\eta_A=0.36$ well below the failure line. The volumetric strain increased during the triaxial compression and a tendency to dilation was noted during the subsequent $q$-constant loading. At point A, the deviatoric stress was maintained constant whereas the corresponding axial strain was 0.45% and the volumetric strain was +0.30%. The rate of increase of the pore water pressure of 0.5 kPa per second is due to an injection of water inside the sample at a rate of $0.5 \text{ mm}^3 \text{ per second}$. This increase corresponds to a volumetric strain rate of $-6.55 \times 10^{-6} \text{ s}^{-1}$ (the negative sign corresponds to a dilation of the sample). This stage is labelled as AB in Figure 2.

During this stage, the axial strain increased slightly and then more rapidly up to point B to reach a value of 1.17%. At this point, which corresponds to a stress ratio of $\eta_B=0.73$, the volume variation is imposed to increase slightly and the sample violently collapses. As shown in Figure 2, an exponential increase in the axial strain occurred, corresponding to a jump of the axial strain from 1.17 to 1.88%. The deviatoric stress, axial strain and volumetric strain are given in this figure.

For CSD2 test, a drained triaxial compression test was applied to the sample to a prescribed value of the deviatoric stress of 54 kPa (point A in Figure 2) corresponding to $\eta_A=0.17$. The corresponding axial and volumetric strains are 0.22 and +0.18%, respectively. Then, the radial stress decreased slowly at a constant rate of 2 kPa/min and the volume of the sample changed as the pore water pressure was kept constant.
<figure><img src="image_3.png" /><figcaption>Figure 2. Evolution of the variables *q* and *ε*<sub>v</sub> versus time. The material is stable during AB and the test is controllable. Collapse occurs at point B, the test is no more controllable during BC and the axial strain rate increases.</figcaption></figure>

At point B where η<sub>B</sub>=0.73 and the axial strain is 0.68%, the deviatoric stress cannot be kept constant any longer and consequently the axial strain rate escalated. The test is no longer controllable in the sense that this loading programme (*q* constant and a given volumetric strain) cannot be maintained. As for CSD1, an exponential increase in the axial strain occurred, corresponding to a jump of the axial strain from 0.68 to 1.48%.

Both the methods give the same stress ratio at collapse η<sub>B</sub>=0.73 and the same evolution of the axial strain before collapse. As the material is in equilibrium from A to B, both methods give the same results. However, a peak of the volumetric strain is noticeable in the CSD2 test whereas it is not in the CSD1 test (Figure 3). This is due to the way of using the pore water DPVC. Since the pore water pressure is prescribed constant during the CSD2 test, the volume is free to vary to keep the pore water pressure constant and then the measurements correspond to the volume variation of the specimen. But at point B, the prescribed volume variation does not correspond any more to the volume change of the sample (which is not measured from that point) for CSD1 test.
<figure><img src="image_4.png" /><figcaption>Figure 3. Influence of the control parameters on the mechanical behaviour of Hostun loose sand under q-constant drained tests for two stress ratios $\eta=0.17$ and 0.35. Sudden collapse occurs at the minimum of the volumetric strain.</figcaption></figure>

It will be shown in Section 3 that the point B corresponds to a loss of stability since a small additional increase of one control variable (here the volumetric strain) leads to failure, as it can be observed in Figures 2 and 3, from point B to C. The collapse is then due to the perturbation of the volumetric change for CSD2 while it is due to the imposed volume change rate during CSD1.

The experimental results reported in Figures 2 and 3 show that the curve depicting the evolving volumetric strain passes through a peak. At this peak (point B), the following control parameters can be considered: $q$ (deviatoric stress applied to the specimen) and $\varepsilon_v$ (volumetric strain). The deviatoric stress is maintained constant. At the volumetric strain peak observed experimentally ($d\varepsilon_v=0$), both control parameters are constant. Experimentally, an additional loading was applied by imposing a small value of $d\varepsilon_v$. Under the effect of this small perturbation, the axial strain abruptly increased, leading to the collapse of the sample. A burst of kinetic energy is found. Indeed, the axial strain rate just before collapse is approximately $2 \times 10^{-3} s^{-1}$, whereas it reaches a value of $1.4 s^{-1}$ after the perturbation is applied. The response of the system evolves abruptly from a quasi static regime toward a dynamic one, with no change in the control parameters: this is therefore a proper bifurcation mode.

As mentioned above, the collapse occurs experimentally at a stress ratio $\eta_B \approx 0.73$, which is well below the failure surface (Mohr–Coulomb limit surface) reached at $\eta \approx 1.32$ and very close to the value obtained during CU tests [17]. We should emphasize that these collapses occur without any localization of the deformation. As a consequence, neither Mohr–Coulomb's limit condition nor a localization criterion can provide a proper framework to explain this collapse mechanism.
<figure><img src="image_5.png" /><figcaption>Figure 4. Axisymmetric drained triaxial compression for a confining pressure $\sigma_3 = 100$ kPa.</figcaption></figure>

## 2.2. Numerical evidence of diffuse failure in DEM modelling of granular materials

Direct numerical simulations with the discrete element method (DEM, Cundall and Strack [21]) were carried out with the code SDEC developed by Donzé and Magnier [22]. The three-dimensional (3D) discrete model (Sibille et al., 2007a) is composed of about 10 000 spheres with a continuous size distribution. The interaction between two contacting spheres is modelled by a linear elastic relation (characterized by the stiffness coefficient $k_n$, with $k_n/d_s = 356$ MPa and $d_s$ being the sphere diameter) in the direction normal to the tangent contact plane, and by an elastic perfectly plastic relation (characterized by the stiffness coefficient $k_t$, with $k_t/k_n = 0.42$, and the friction angle $\phi_\mu = 35^\circ$) for the direction included in the tangent contact plane.

Simulations presented in this section were performed with a dense numerical sample (characterized by a void ratio $e=0.618$ and a coordination number $z=4.54$), exhibiting a well-marked dilatant behaviour during an axisymmetric ($\varepsilon_2 = \varepsilon_3$ and $\sigma_2 = \sigma_3$) drained triaxial compression for a confining pressure $\sigma_3 = 100$ kPa (Figure 4).

We analyse the mechanical behaviour of numerical specimens along proportional strain loading paths defined by [23]:

$d\varepsilon_1$ = positive constant

$d\varepsilon_2 = d\varepsilon_3$ (axisymmetric condition) (1)

$d\varepsilon_1 + 2Rd\varepsilon_3 = 0$ ($R$ constant for a given path)

According to the value of $R$, the path is dilatant ($0<R<1$) or contractant ($R>1$); the undrained condition (i.e. isochoric loading path) corresponds to $R=1$. Such loading paths can be seen as a generalization of the undrained tests presented in Section 2.1. Simulations of proportional strain paths, from an isotropic confinement at $\sigma_3 = 100$ kPa, are presented in Figure 5. Since the numerical sample is dense, the isochoric loading path does not lead to a peak of the stress deviator $q = \sigma_1 - \sigma_3$, whereas for more dilatant path (i.e. for lower $R$ values) $q$ passes through a maximum value. A vanishing of stresses (suggesting a liquefaction phenomenon) is even computed for sufficiently low values of $R$ ($R \le 0.825$).

For analyzing the type of test considered in the above, it is important to express the second-order work in terms of the proper energy conjugate variables as suggested by Nova [24], i.e.

$$d_2 W = d\sigma_1 d\varepsilon_1 + 2d\sigma_3 d\varepsilon_3 = d\varepsilon_1(d\sigma_1 - d\sigma_3/R) + d\sigma_3(d\varepsilon_1 + 2Rd\varepsilon_3)/R \quad (2)$$

In Equation (2), we find that the proper response parameters with respect to the control parameters defined in Equation (1) are $\xi_1 = (\sigma_1 - \sigma_3/R)$ and $\xi_2 = \sigma_3/R$. Thus, the deviatoric stress $q$ constitutes a response parameter only when $R=1$. Figure 5(c) shows the evolution of the response parameter $\xi_1$ which passes through a peak for value of $R \le 0.825$. Since the peak at $R=0.825$ corresponds to
<figure><img src="image_6.png" /><figcaption>Figure 5. Simulations of proportional strain loading paths for R = 0.7; 0.8; 0.825; 0.85; 0.875; 0.9; 1.0.</figcaption></figure>

a mobilized friction angle $\phi_m = 15.5^{\circ}$, the numerical specimen fails well below the Mohr–Coulomb failure criterion which refers to a friction angle of $\phi = 24.7^{\circ}$ (which is consistent with the measured value of an assembly of glass ballotini).

Referring far back to the experiments described in Section 2.1, we showed that a diffuse mode of failure well within the Mohr–Coulomb failure criterion was obtained under proper controlled parameters (load controlled). Hence, the next step in this numerical endeavour is to explore possibilities of rapid collapse by changing the load control parameters so that now the proportional strain path is run as follows:

$$ 
\begin{aligned}
d\sigma_1 - d\sigma_3/R &= d\xi_1 = \text{positive constant} \\
d\varepsilon_2 &= d\varepsilon_3 \quad (\text{axisymmetric condition}) \\
d\varepsilon_1 + 2Rd\varepsilon_3 &= d\eta_2 = 0 \quad (R \text{ constant for a given path})
\end{aligned}
\qquad (3) $$

The computational implementation of such mixed loading programme is described in Sibille et al. [25].

Figure 6 shows a comparison of the simulated responses between the two modes of control for R = 0.825. As long as the value of $\xi_1$ is less than the peak value, the response of the numerical sample is very similar irrespective of the control mode. For a control with $\xi_1$, the peak value of $\xi_1$ is slightly exceeded, but finally a decrease is observed (Figure 6(a)), whereas the condition $d\xi_1 > 0$ is still being tried to be fulfilled. When the peak value of $\xi_1$ is approached, the response is characterized by a transition from a quasi-static regime to a dynamical one, as shown by the evolution of the kinetic energy of the numerical sample (i.e. the sum of the kinetic energy of each grain) in Figure 6(b). Owing to the dynamic response, inertial terms are no longer negligible: for instance, when $\sigma_3$ vanishes, $\sigma_1$ is about 23 kPa (Figure 6(c)). Hence, the slight exceeding of the $\xi_1$ peak value could be related to the dynamical response of the sample. Numerical results show that, for a control with the parameters $\xi_1$ and $\eta_2 = \varepsilon_1 + 2R\varepsilon_3$, the peak of $\xi_1$ constitutes a generalized
<figure><img src="image_7.png" /><figcaption>Figure 6. (a), (b) Comparison of simulated responses for a proportion strain loading path controlled by $\varepsilon_1$ or $\xi_1 = \sigma_1 - \sigma_3/R$ and for $R=0.825$; (c) stress evolutions for a control with $\xi_1$.</figcaption></figure>

limit state [23] that cannot be exceeded, or in other terms there is a loss of controllability [16, 24] of the loading programme leading to a sudden failure (as a sudden liquefaction) of the sample.

We now consider two mechanical states at the equilibrium, $ES_a$ and $ES_b$, reached along the proportional strain loading path controlled through parameters defined in Equation (1) for $R=0.825$. These states differ by their relative positions with respect to the peak of $\xi_1$. For instance, $ES_a$ is reached slightly before the peak, whereas $ES_b$ is located slightly after the peak (Figures 5(c) and 6(a)). From $ES_a$ and $ES_b$, the numerical sample is forced to stay at the following equilibrium states:

$$ 
\begin{aligned}
d\xi_1 &= (d\sigma_1 - d\sigma_3)/R = 0 \\
d\varepsilon_2 &= d\varepsilon_3 \quad (\text{axisymmetric condition}) \\
d\eta_2 &= d\varepsilon_1 + 2Rd\varepsilon_3 = 0
\end{aligned}
\qquad (4) $$

Under these conditions, if simulations are run from states $ES_a$ or $ES_b$, the sample does not evolve and stays at the equilibrium states that are imposed. Then, a 'small' perturbation is applied to the sample. This perturbation is realized by imposing, at a given time, a given velocity in a random direction to some 'floating' spheres (simulations are performed without gravity) chosen randomly. The perturbation corresponds to an external supply of kinetic energy of $10^{-4}$ J to the sample. By comparison, the maximum value of kinetic energy computed for a control mode defined by Equation (1) is $1.25 \times 10^{-3}$ J.

Figure 7 presents a comparison of responses of the numerical sample after such a perturbation from equilibrium states $ES_a$ and $ES_b$ controlled by Equation (4). From the equilibrium state $ES_a$, the kinetic energy stays quite small and then vanishes; after some low fluctuations the control parameters $\xi_1$ and $\eta_2$ stay constant and equal to their initial value. Only a slight decrease is observed for the radial stress $\sigma_3$. Thus, the sample reaches a new equilibrium state very close to the initial one. From the equilibrium state $ES_b$ instead, the kinetic energy increases exponentially without vanishing again. On the other hand, we do not succeed in imposing the loading programme defined in Equation (4) and stresses vanish (loss of controllability). Therefore, it seems that, for
<figure><img src="image_8.png" /><figcaption>Figure 7. Simulated responses for a perturbation of the equilibrium states ES<sub>a</sub> and ES<sub>b</sub> controlled by the loading programme defined in Equation (4).</figcaption></figure>

the specific control parameters chosen, the maximum value of ξ₁ constitutes a boundary between
two domains: after the peak sudden failures could easily be triggered by a small perturbation (loss
of sustainability—[6]); whereas, for the same control parameters, failure does not develop from
an equilibrium state localized before the peak.

3. THEORETICAL BACKGROUND

Throughout this section, time and spatial increments of any variable $\psi$ will be distinguished by denoting $d\psi$ as the increment in time of $\psi$ (defined as the product of the particulate derivative $\psi$ by the infinitesimal time increment $dt$) with respect to a given frame, and by denoting $\delta\psi$ as the spatial increment of $\psi$.

3.1. *A general framework: Loss of material stability*

The generally admitted definition of stability comes from Lyapunov [4], which has considered
the influence of small perturbations of the initial positions and celerities of planets on their later
motion around the sun. Applied to continuum mechanics this definition can be written in a proper
mathematical framework under the following form:

A mechanical stress–strain state (for a given material after a given strain history) is called stable if for any positive scalar $\varepsilon$ there is a positive scalar $\eta$ such that for any loading limited in norm by $\eta$, the response will remain limited in norm by $\varepsilon$ ($\eta$ being a function of $\varepsilon$).

This definition implies that all failure states are unstable, since, as recalled in the Introduction,
failure is due to the existence of some limit states that are impossible to be passed without large
strains, cracks, etc. Note that the reciprocal is a false statement: all unstable states do not lead necessarily to failure because of the directional nature of failure as it will be discussed later on.

On the same way failure implies loss of uniqueness (there are an infinity of material responses), the reciprocal being also false (a counterexample is given in the softening regime by the loading-unloading branches).

The classical view to analyse failure is to consider a plastic limit condition, defined as the surface plotted in the six-dimensional (6D) stress space limiting the stress states which can be reached by the material point of a given medium. On such a surface for very small incremental stresses it is possible to obtain large strains.

Thus, for the rate-independent constitutive relation $d\sigma_{ij} = L_{ijkl} d\epsilon_{kl}$ in the 3D space or $d\sigma_{\alpha} = M_{\alpha\beta} d\epsilon_{\beta}$ in the 6D associated space, we have the necessary condition:

$$ \det M = 0 \qquad (5) $$

at least for all incrementally piecewise linear relations.

This equation corresponds to the plastic limit condition and the related relation:

$$ M d\epsilon = 0 \qquad (6) $$

corresponds to the so-called material flow rule (possibly with a vertex for multiple plastic mechanisms or for incrementally non-linear constitutive relations).

However, the experiments showed that some failure modes with some localization patterns (shear banding phenomenon or compaction-dilatation band formation) do not fulfil this criterion.

Taking into account the continuity of the stress vector through the band and the kinematics of the shear band, a proper localization criterion was established [26]:

$$ \det(nLn) = 0 \qquad (7) $$

where $n$ is the unit normal to the band.

However, in Section 2 of this paper, it was shown that some proper failure modes do not fulfil neither the plastic limit condition nor the localization condition, as the initiation of liquefaction of loose undrained sands.

To build another criterion, let us come back to the fact that failure is necessarily an unstable state and let us consider Hill's criterion of instability [8]. Applied to a material point, this criterion states that it is unstable if the deformation can be pursued in a given direction without any input of energy from outside. In that direction the second-order work is taking a negative value. We will see later on that more precisely there is a burst of kinetic energy which is directly linked to these negative values of the second-order work.

For small strains and by neglecting geometrical effects, second-order work writes: $d_2 W = d\sigma : d\epsilon$ in 3D space. It comes in 6D space

$$ d_2 W = d\sigma d\epsilon = {}^t d\sigma d\epsilon = {}^t d\epsilon {}^t M d\epsilon = {}^t d\epsilon M_s d\epsilon \qquad (8) $$

The positivity of second-order work corresponds to the definite positivity of $M$ (or $M_s$).

If we assume that the eigenvalues of $M_s$, assumed to be initially strictly positive, are varying in a continuous manner with the loading parameter, it appears that the positivity of second-order work is equivalent to the positive value of $\det M_s$.

Thus, this new criterion of failure can be written as:

$$ \det M_s \le 0 \qquad (9) $$

Now a clarifying discussion emerges:

For associate materials, the elastoplastic matrix $M$ is symmetric, and conditions (5) and (9) are equivalent.

On the contrary for nonassociate materials (like geomaterials), $M$ is nonsymmetric and according to linear algebra $\det M_s$ will vanish before $\det M$ (according to a loading parameter).
Thus, condition (9) corresponds to a large domain of bifurcations in the stress space whose boundaries are given by the plastic limit condition (det M = 0) and the first vanishing values of second-order work (det M_s = 0).

Moreover, as the second-order work criterion corresponds to negative values of a quadratic form, it has the strange but essential feature to be directional. At a given stress state, it will be positive in some directions, negative in others. We will establish in Section 4.1.3 that in 3D stress space there are some elliptical cones ‘of instability’.

Let us now establish the link between negative values of second-order work and bursts of kinetic energy and for that let us consider an initial configuration C₀ endowed with a volume V₀. After a given loading history, the system evolves to a strained configuration C at time t and occupies a volume V, in equilibrium under a prescribed external loading.

It can be shown [6] that the kinetic energy E_c of the system at the subsequent time t+dt can be expressed following a semi-Lagrangian description as:

$$2E_c(t+dt) = \int_{\Gamma_0} d_2L \, d\Gamma_0 - d_2W \qquad (10)$$

where d₂L = dΠᵢⱼNⱼduᵢ is the incremental work done due to the perturbation of the external control variables on the element boundary, whereas Π denotes the first Piola–Kirchoff stress tensor, N is the current normal to the boundary Γ at the point considered, du is the incremental displacement field, and Γ₀ is the surface bounding the initial volume V₀. Moreover, the second-order work in a semi-Lagrangian description is given by:

$$d_2W = \int_{V_0} d\Pi_{ij} \frac{\partial(du_i)}{\partial X_j} dV_0 \qquad (11)$$

A sudden increase in kinetic energy as signaled by E_c(t+dt)>0 in Equation (10) occurs whenever

$$d_2W < \int_{\Gamma_0} d\Pi_{ij} N_j du_i d\Gamma_0 \qquad (12)$$

In the absence of changes in external boundary loads (for example), the boundary integral in Equation (12) vanishes, and thus:

$$d_2W < 0 \qquad (13)$$

Since the effective collapse of a material specimen implies a burst in kinetic energy, Equations (10), (12) and (13) establish a fundamental connection between failure and the loss of positivity of the second-order work at the material point scale. It follows that exploring failure in homogeneous laboratory tests and interpreting results at the elementary scale can be very insightful. In these lab experimental tests, the loading is applied to the specimen through adequate control parameters while both homogeneous stress and strain fields are assumed. As an example, the subsequent analysis will be restricted to the case of a cubical specimen loaded in axisymmetric conditions, where subscript 1 refers to the vertical direction, and subscript 3 refers to one of the horizontal directions.

Let us first set the incremental external stress applied to each point on the boundary of the specimen as dsᵢ = dΠᵢⱼNⱼ. If both fields ds and du are piecewise on the boundary of the specimen, it follows that:

$$\int_{\Gamma_0} d\Pi_{ij} N_j du_i d\Gamma_0 = A(ds_1 du_1 + 2ds_3 du_3) \qquad (14)$$

where A denotes the cross-sectional area of the cube in the initial configuration. In general, any set of control parameters (C₁ and C₂) is constituted by two independent linear combinations of s₁ and s₃, or u₁ and u₃. Likewise, the related response parameters (R₁ and R₂) are linear combinations of u₁ and u₃ or s₁ and s₃. With the assumption of energy equivalence that:

$$A(ds_1 du_1 + 2ds_3 du_3) = A(dC_1 dR_1 + 2dC_2 dR_2) \qquad (15)$$
<figure><img src="image_9.png" /><figcaption>Figure 8. Illustration of a loss of controllability with a rigid beam.</figcaption></figure>

Equation (10) finally reads in homogeneous conditions as:

$$2E_c(t+dt) = A_1(dC_2dC_2dR_2) - d_2W \quad (16)$$

At this juncture, two lines of reasoning for discussing material stability can be distinguished. The first line of reasoning corresponds to the loss of controllability approach. In the absence of failure, a unique incremental response ($dR_1$ and $dR_2$) exists for a given incremental loading ($dC_1$ and $dC_2$) expressed as:

$$dR_i = M_{ij}dC_j, \quad i, j = 1, 2, 3 \quad (17)$$

with all eigenvalues of rearranged constitutive matrix $M$ being real and positive. When $M$ becomes singular (at least one eigenvalue is nil), the uniqueness of the response is lost. Indeed, the response of the system is no longer governed by Equation (15), but actually follows Equation (14). From a physical view point, the loss of uniqueness corresponds to the fact that the loading path considered is no longer controllable. Interestingly, this approach corresponds to the loss of controllability concept introduced by Nova [24].

A second view point that can be followed in the analysis of material stability consists in investigating under which conditions the system, which is in an equilibrium mechanical state at a given time, can suddenly evolve under the application of a certain infinitesimal perturbation of the control parameters. The existence of both control and response parameters such that the boundary term $A(dC_1dR_1+2dC_2dR_2)$ is zero and a negative second-order work (for example) results into an increase of kinetic energy, as implied by Equation (16). This approach corresponds to the loss of sustainability [27].

Both loss of controllability and loss of sustainability approaches are presented in details in the following.

## 3.2. Loss of controllability approach

Consider the rigid beam of length $l$, shown in Figure 8, that is subject to a distributed load $p$. If $K$ is the value of the torsional stiffness at the hinge, the rotation of the beam $\theta$ is:

$$\theta = \frac{pl^2}{2K} \quad (18)$$

Imagine now that we add a horizontal load $P$, assumed to act in the new geometric configuration. By assuming that the rotation remains small, albeit increasing, the beam rotation is now equal to:

$$\theta = \frac{pl^2}{2(K - Pl)} \quad (19)$$

It is seen that the stiffness of the system is therefore reduced by the effect of the horizontal load. Assume now that we next give a perturbation $dP$ to the horizontal load. The relation between the known applied load perturbation and the consequent rotation is then:

$$dP = \frac{(K - Pl)^2}{2pl^3/2} d\theta \quad (20)$$
As long as $Pl < K$, a unique solution exists to the incremental relation between axial load and rotation, that can be considered as a generalized stress-strain relationship for the discrete system at hand. However, when

$$P = K/l \qquad (21)$$

the term within parentheses in Equation (20), which represents the overall stiffness of the system (elastic plus geometric), vanishes and many values $d\theta$ are possible for zero $dP$, i.e. at constant horizontal load there exist multiple solutions. The load value for which this occurs is known as the latent instability load [28] for the beam and corresponds therefore to the occurrence of a possible bifurcation of the system response. Of course, this solution is approximate, since the hypothesis of small rotations is no longer valid when $P$ approaches $K/l$, but it is anyway a good approximation to the real instability load from an engineering viewpoint.

From the point of view of the experimenter that aims at progressively increasing the loading to investigate the system response, Equation (21) gives the condition at which the controllability of the chosen loading process is lost, because the system is free to choose between infinitely many configurations. We see that the terms loss of controllability of the loading programme, latent instability and point of bifurcation are referred to the same condition, in the light of the example discussed so far, at least.

Consider an element test on a geomaterial specimen in a true triaxial apparatus, where either the principal stresses or principal strains are controlled. If strains are assumed to be small, and the stress-strain state assumed to be uniform within the specimen, the relation between the increments of the external forces acting on the specimen and the corresponding displacements can be expressed directly in terms of stress and strain rates as follows:

$$d\sigma = Dd\epsilon \qquad (22)$$

where a dash indicates effective stresses, $D$ being the stiffness matrix. Stresses are in fact uniquely determined by the equilibrium conditions and strains by compatibility. The controllability of the loading programme, when external forces, i.e. stresses, are controlled, is lost whenever:

$$\det D = 0 \qquad (23)$$

When Equation (23) is fulfilled, there exist infinitely many strain rates at constant stress state and a bifurcation occurs. It is worth noting that such a bifurcation exists only if stresses are controlled. If on the contrary we impose strains,

$$d\epsilon = C d\sigma' \qquad (24)$$

where $C$ is the compliance matrix, when Equation (23) is fulfilled, the determinant of the compliance matrix is still positive and a unique solution exists. In particular, if $d\epsilon$ is chosen in such a way to be proportional to the eigenvector of matrix $D$ (associated with the vanishing eigenvalue),

$$d\sigma' = 0 \qquad (25)$$

and a peak occurs in the stress-strain curve.

Clearly in a test we can control partly strain rates and partly stress rates. Think for instance of a drained constant cell pressure test in a conventional triaxial apparatus, where axial strains are controlled, or of an oedometer test, where the axial load is imposed and the horizontal strain is prevented. We can then write again the constitutive relation given by Equation (22) by partitioning the stress and strain vectors between controlling and response variables as in Equation (15), but by putting at the left-hand side the controlling, i.e. known, variables:

$$\left\{ \begin{matrix} d\sigma_{\alpha} \\ d\epsilon_{\beta} \end{matrix} \right\} = \left[ \begin{matrix} D_{\alpha\alpha} - D_{\alpha\beta}D_{\beta\beta}^{-1}D_{\beta\alpha} & D_{\alpha\beta}D_{\beta\beta}^{-1} \\ -D_{\beta\beta}^{-1}D_{\beta\alpha} & D_{\beta\beta}^{-1} \end{matrix} \right] \left\{ \begin{matrix} d\epsilon_{\alpha} \\ d\sigma_{\beta} \end{matrix} \right\} \quad (\alpha, \beta = 1, 2, 3) \qquad (26)$$
Again, controllability of this loading programme is lost when the determinant of the matrix of Equation (26) vanishes. By exploiting the Schur [29] theorem, it can be proven [5, 24] that this happens when:

$$ \det \mathbf{D}_{\alpha\alpha} = 0 \qquad (27) $$

where $\mathbf{D}_{\alpha\alpha}$ is a principal submatrix of the stiffness matrix $\mathbf{D}$. Whenever a principal minor of $\mathbf{D}$ is nil, therefore, there is a particular loading programme for which the control is lost and a homogeneous bifurcation can occur.

It is important now to investigate for which load levels such a loss of controllability occurs. We note first that if $\mathbf{D}$ is positive definite all its principal minors are positive. Thus, no loss of controllability is possible. Furthermore if $\mathbf{D}$ is symmetric, no minor can be negative if the determinant is positive. For a symmetric matrix, therefore the loss of positive definiteness coincides with the nullity of the determinant of the stiffness matrix, which coincides in turn with the loss of controllability of a fully stress-controlled test. This is the usual condition associated historically with the plastic limit condition. Indeed, for an elastoplastic hardening soil, for instance, Equation (23) is fulfilled when the hardening modulus $H$ becomes zero. If the flow rule is associated, the stiffness matrix is symmetric and therefore no bifurcation can occur for positive values of $H$, i.e. in the hardening regime.

If the flow rule is non-associated however, the stiffness matrix is not symmetric and the loss of positive definiteness:

$$ \det \mathbf{D}_s = 0 \qquad (28) $$

occurs always 'before' the ordinary failure condition, i.e. for a positive hardening modulus [5]. The Ostrowski and Taussky [30] theorem ensures in fact that for a positive semi-definite matrix:

$$ \det \mathbf{D} \geq \det \mathbf{D}_s \qquad (29) $$

Some minors of $\mathbf{D}$ can become zero therefore when the determinant of the stiffness matrix is still positive, i.e. in the hardening regime. For instance, if we decide to follow a load path in which it is controlled the axial strain increment in one principal direction and the principal stress increments in the other two (path a in Figure 9 [16]), loss of controllability occurs after the loss of positive definiteness but before the ordinary failure locus as shown in Figure 9. The constitutive model

<figure><img src="image_10.png" /><figcaption>Figure 9. Loci of loss of controllability under different loading programmes (controllability of path is lost when $C_{11}=0$) (from [16]).</figcaption></figure>
used here for the predictions is the elastoplastic strain-hardening model combining isotropic and kinematic hardening proposed by di Prisco et al. [31].

It is quite interesting to note that, when Equation (28) is met, there is one particular direction in the stress space for which the second-order work is zero. The positive definiteness of the second-order work was proven to be a sufficient condition for stability [8]. Since we have shown that no minor can be nil before the loss of positive definiteness of the stiffness matrix, this result is compatible with the fundamental result obtained by Hill. On the other hand, the analysis so far presented can be pushed a step forward to show that the positive definiteness of the second-order work is also a necessary condition for stability.

In geotechnical tests we often control linear combination of stresses or strains. For instance, in an undrained test on a saturated specimen we impose that the sum of the principal strain increments is zero. On the other hand, controlling the axial load at constant cell pressure implies we control the difference between the axial and radial stress. In general we can therefore define generalized stress and strain variables, e.g.:

$$d\xi = T_{\sigma}d\sigma^{t} \qquad (30)$$

$$d\eta = T_{\epsilon}d\epsilon \qquad (31)$$

that are linked together by a generalized stiffness matrix:

$$\Delta = T_{\sigma}DT_{\epsilon}^{-1} = T_{\sigma}DT_{\sigma}^{t} \qquad (32)$$

The last equality results from the necessary energy equivalence between natural and generalized variables. The work input of the former must be the same of that of the latter in fact.

In principle, we can choose the transformation matrices as we like. In particular we can take:

$$T_{\epsilon} = T_{D}^{t} = T_{\sigma} \qquad (33)$$

where $T_D$ is the orthogonal matrix diagonalizing the symmetric part of the stiffness matrix. If $D_s$ is singular, it is possible to find a particular load programme for which a homogeneous bifurcation can occur. In this case, in fact,

$$\Delta = T_{D}^{t}DT_{D} \qquad (34)$$

and then

$$\Delta = \begin{bmatrix} 0 & e^{T} \\ -e & \Lambda \end{bmatrix} \qquad (35)$$

where the symmetric part of $\Lambda$ is a positive-definite diagonal matrix. Its diagonal terms are the positive eigenvalues of $D_s$, while $e$ is a 2 component vector. If we now take as control variables $\{d\xi_1 d\xi_2^t\}$ the generalized constitutive law can be written as:

$$\begin{Bmatrix} d\xi_1 \\ d\eta_2 \end{Bmatrix} = \begin{bmatrix} e^T \Lambda^{-1} e & e^T \Lambda^{-1} \\ \Lambda^{-1} e & \Lambda^{-1} \end{bmatrix} \begin{Bmatrix} d\eta_1 \\ d\xi_2 \end{Bmatrix} \qquad (36)$$

It is apparent that the determinant of the matrix of Equation (36) is zero, since the first row is a linear combination of the others. Therefore, under constant control variables, infinitely many responses are possible. When the Hill stability condition is violated first, a bifurcation can occur under a particular load path. The positive definiteness of the second-order work is therefore also a necessary condition for uniqueness of the incremental response, in the wider sense discussed above.

It can be shown [32] that some examples of load paths for which the second-order work becomes zero first are special tests in which it is imposed the rate of sample dilation (see Equation (1)).

A special case of this test is the undrained test. Imposimato and Nova [16] showed that the locus for the loss of controllability is a cone in the stress space intermediate between the cone
represented by Equation (28) and that given by Equation (23). In particular, in the axisymmetric triaxial test, it exists a well-known locus of stress deviator peaks [33]. In this case, the loss of controllability cone intersection with the axisymmetric plane is a straight line passing through the origin, which coincides with the so-called Lade's instability line [34]. It is worth mentioning that for this stress level under constant stress state a possible instability in the form of uncontrolled pore pressure generation is possible, as shown by Imposimato and Nova [35].

### 3.3. Loss of sustainability approach

To detect the vanishing of the second-order work, the directional analysis (initially introduced by Gudehus [36] to construct the so-called response-envelopes) is particularly convenient. Starting from a given equilibrium state, incremental stress vector probes $ds$ are imposed along various directions within the Rendulic plane (axisymmetric stress plane) with a fixed norm. Assuming that the equilibrium state considered is not on the plastic limit surface, a unique incremental displacement $du$ exists for any incremental stress loading. This, in fact, corresponds to a quasi-static transition from one equilibrium state to another one. As a consequence, both incremental stress and strain fields are homogeneous within the specimen, so as at any point of the specimen, principal values of both $d\Pi$ and $\nabla u$ verify $d\Pi_i = ds_i$ and $\partial u_i / \partial X_i = du_i / L$, where $L$ is the length of the edge of the specimen. The normalized second-order work can be therefore expressed in terms of the boundary variables as:

$$d_2 W = \frac{ds_1 du_1 + 2ds_3 du_3}{\sqrt{ds_1^2 + 2ds_3^2} \sqrt{du_1^2 + 2du_3^2}} \quad (37)$$

The normalized second-order work is basically a directional quantity that is independent of the norm of the individual vectors. Let us assume that for a given stress probe $ds$ the corresponding second-order work can become negative. Thus, the basic question that arises is to find whether there exists a certain set of control parameters with the related response parameters such that $d_2 W < 0$ and $A(ds_1 du_1 + 2ds_3 du_3) = A(dC_1 dR_1 + 2dC_2 dR_2) = 0$.

Given the directional nature of the second-order work, the incremental stress direction $\alpha_s = \tan^{-1}[ds_1/(\sqrt{2}/ds_3)]$, if it is inside an instability cone, leads to a strictly negative $d_2 W$ for any norm of the stress vector $\{ds_1, \sqrt{2}ds_3\}^T$. Setting $1/R = \tan \alpha_s$, the control parameter $C_1 = s_1 - \sqrt{2}s_3/R$ is introduced. Since:

$$ds_1 du_1 + 2ds_3 du_3 = \left( ds_1 - \frac{\sqrt{2}}{R} ds_3 \right) \delta u_1 + \frac{1}{R} (du_1 + \sqrt{2}R du_3) \sqrt{2} ds_3, \quad (38)$$

it is convenient to define $C_2 = u_1 + \sqrt{2}Ru_3$ as the second control parameter. Once we choose the proper control parameters $C_1$ and $C_2$, we additionally require that they remain constant. Then, $ds_1 du_1 + 2ds_3 du_3 = 0$ and Equation (10) reads:

$$2E_c(t+dt) = -d_2 W \quad (39)$$

As the second-order work associated with this incremental loading programme becomes lower than the boundary integral of Equation (10), then the kinetic energy increases, leading to the collapse of the specimen. As it has been shown in Section 2 in both experimental tests and numerical simulations using a discrete element method, when control parameters are conveniently chosen and maintained constant, an infinitesimal perturbation is sufficient to provoke the collapse of the specimen.

Of course, such a collapse can occur only if at least one incremental stress direction exists leading to the vanishing of the second-order work. The mechanical states satisfying this condition belong to the bifurcation domain. The mechanical states that belong to the limit of the bifurcation domain are such that only a single incremental stress direction exists for which the second-order work vanishes. It is worth noting that the occurrence of failure at any state within the bifurcation domain requires that the control parameters are conveniently chosen and maintained constant at
a bifurcation point, together with the application of a certain infinitesimal perturbation. If these conditions are not fulfilled, no collapse mechanism will be observed.

Finally, a last remark can be added concerning both notions of loss of controllability and sustainability. We assume that the mechanical state considered belongs to the bifurcation domain. The loss of sustainability requires a convenient choice of $C_1$ and $C_2$ under the additional prescription that $dC_1=0$ and $dC_2=0$. Then the application of an infinitesimal perturbation will induce collapse of the specimen. It must be noted that imposing $dC_1=da$ and $dC_2=db$ where $da$ and $db$ are any two infinitesimal values, corresponds to a special class of infinitesimal perturbation of the stationary condition $dC_1=0$ and $dC_2=0$. Following the loss of sustainability approach, the specimen will collapse. Moreover, it can be shown that the generalized tangent constitutive tensor relating both incremental control and response parameters is singular in this case [37]. Thus, in line with the theoretical developments presented in the previous section, this is also a loss of controllability. The main difference between these two notions is that the controllability emphasizes the notion of a loading programme (i.e. a certain path within the mixed strain-stress space), whereas the notion of sustainability is confined to arguments of stability around a mechanical state (i.e. a certain point within the stress space).

# 4. MACRO- AND MICRO-MODELLING OF SOIL FAILURE

## 4.1. Phenomenological approaches

In this section, three constitutive models of increasing complexity are presented to interpret the experimental and theoretical observations of the previous sections. It will be shown that despite the different approaches they all qualitatively predict what is experimentally observed.

### 4.1.1. Non-associated Cam–Clay model.
We herein introduce a basic model for soil behaviour that has the necessary and sufficient ingredients to qualitatively describe the behaviour of loose sands, as observed in triaxial compression tests. Because of the extreme simplicity of the model, possibly the simplest conceivable with the aforementioned characteristics, the conditions for the occurrence of instabilities can be analytically determined. The hierarchy of instability conditions, such as nullity of second-order work, spontaneous generation of pore water pressure or loss of controllability under assigned deformation path can be easily established. This makes the conceptual interpretation of the experimental results shown in Section 2 easier. It should be borne in mind, however, that for a better quantitative agreement with experimental data, reference should be made to more realistic models (e.g. Sinfonietta Classica [38]) even though they in turn do not allow a simple analytical treatment.

It is first assumed that the behaviour of loose sand (or normally consolidated clay) can be described by an elasto-plastic strain hardening model. It is then possible to define a history-dependent domain, within which soil behaviour is considered to be fully reversible. In the so-called triaxial plane $q-p'$ ($q=\sigma_1-\sigma_3$, $p'=(\sigma'_1+2\sigma'_3)/3$), the frontier of this domain, the yield locus, is given by the following equation:

$$f=q+mp'\ln\frac{p'}{p_c}=0 \qquad (40)$$

Equation (40) coincides with the expression of the loading function for the original Cam–Clay model [39], except for the symbol $m$ which is different from the critical state parameter $M$. This is used instead to characterize the plastic potential $g$ as:

$$g=q+Mp'\ln\frac{p'}{p_g}=0 \qquad (41)$$
so that plastic strain rates can be derived as:

$$ \left\{ \begin{matrix} d\epsilon_v^p \\ d\epsilon_d^p \end{matrix} \right\} = \Lambda \left\{ \begin{matrix} \frac{\partial g}{\partial p'} \\ \frac{\partial g}{\partial q} \end{matrix} \right\} \qquad (42) $$

where $\Lambda$, known as plastic multiplier, is a positive quantity when plastic strains occur and is zero upon unloading-reloading ($\epsilon_v = \epsilon_1 + 2\epsilon_3$, $\epsilon_d = 2(\epsilon_1 - \epsilon_3)/3$).

The fact that $m \neq M$ implies that $f \neq g$, i.e. a non-associated flow rule is assumed. Actually, to describe soil behaviour in a realistic way, $m < M$. In Equation (40), $p_c$ is the so-called pre-consolidation pressure, geometrically given by the intersection of the yield locus with the hydrostatic axis, whereas $p_g$ is the corresponding value for the plastic potential. For a given stress state fulfilling Equation (40), these two quantities are connected by the following relationship:

$$ p_g = p' \left( \frac{p'}{p_c} \right)^{-m/M} = p_c e^{((\eta/M) - \eta/m)} \qquad (43) $$

The pre-consolidation pressure depends on the plastic strains experienced by the material. For loose sand, we can assume that only plastic volumetric strain influences the value of $p_c$. Under isotropic loading, a convenient expression for describing the variation of specific volume ($v = 1 + e$) with isotropic pressure is [40],

$$ \ln v = c_0 - \lambda' \ln \frac{p_c}{p_1} \qquad (44) $$

where $c_0$ is a constant. By differentiating Equation (44) and taking account the definitions of specific volume and volumetric strain we have:

$$ \frac{dv}{v} = d\epsilon_v = \lambda' \frac{dp_c}{p_c} \qquad (45) $$

Assuming that upon unloading a similar relation holds true, but with a different compressibility constant $k'$:

$$ d\epsilon_v^e = \kappa' \frac{dp_c}{p_c} \qquad (46) $$

the plastic volumetric strain increment can be obtained by subtracting side by side Equation (46) from Equation (45), so that, eventually, we get:

$$ \frac{\partial p_c}{\partial \epsilon_v^p} = \frac{p_c}{\lambda' - k'} \qquad (47) $$

When plastic strains occur, the loading function $f$ remains constant, so that:

$$ df = \frac{\partial f}{\partial p'} dp' + \frac{\partial f}{\partial q} dq + \frac{\partial f}{\partial p_c} dp_c = 0 \qquad (48) $$

Noting Equations (40), (42), (43) and (47), the plastic multiplier can be calculated as

$$ \Lambda = \frac{(m - \eta)dp' + dq}{H} \qquad (49) $$

where the hardening modulus $H$ is defined as:

$$ H = - \frac{\partial f}{\partial p_c} \frac{\partial p_c}{\partial \epsilon_v^p} \frac{\partial g}{\partial p'} = \frac{mp'}{\lambda' - \kappa'} (M - \eta) \qquad (50) $$
It is apparent that $H$ is positive for $\eta<M$ (hardening), negative for $\eta>M$ (softening) and is zero when $\eta=M$. The constitutive law can then be written as:

$$ \left\{ \begin{matrix} d\epsilon_v \\ d\epsilon_d \end{matrix} \right\} = \left[ \begin{matrix} C_{pp} & C_{pq} \\ C_{qp} & C_{qq} \end{matrix} \right] \left\{ \begin{matrix} dp' \\ dq \end{matrix} \right\} = \left[ \begin{matrix} \frac{k'}{p'} + \left\langle \frac{(\lambda' - k')(m - \eta)}{mp'} \right\rangle & \left\langle \frac{\lambda' - k'}{mp'} \right\rangle \\ \left\langle \frac{(\lambda' - k')(m - \eta)}{mp'(M - \eta)} \right\rangle & \left\langle \frac{\lambda' - k'}{mp'(M - \eta)} \right\rangle \end{matrix} \right] \left\{ \begin{matrix} dp' \\ dq \end{matrix} \right\} \quad (51) $$

The brackets appearing in the matrix of Equation (51) (McCauley brackets) mean that the element of the matrix within parentheses has the value indicated when plastic strains occur ($f=0, df=0$), while it is zero otherwise, i.e. for stress paths causing elastic strains only.

It is apparent from Equation (51) that deviatoric strain rates tend to become infinite when the stress ratio tends to $M$, irrespective of the stress path followed (provided plastic strains take place). The state $\eta=M$ is therefore the so-called limit condition. In fact, for this value of the stress ratio the determinant of the stiffness matrix $D$ vanishes. There exists therefore an eigen solution of the type:

$$ \left\{ \begin{matrix} d\epsilon_v \\ d\epsilon_d \end{matrix} \right\} = \left\{ \begin{matrix} 0 \\ \beta \end{matrix} \right\} \quad (52) $$

where $\beta$ is an arbitrary scalar, that can be unlimitedly large. This means that a load-controlled test cannot be controlled any longer. The ordinary failure condition appears then in this way as a loss of controllability condition for stress controlled tests.

Equation (51) can be easily integrated for assigned loading conditions. It can be easily shown that the essential feature of the behaviour of loose sands can be captured.

In a particular, the strains occurring in an undrained test, starting from an isotropic pressure $p_0$ can be calculated. By imposing the no volume change condition, the effective stress path can be derived from Equation (51) so that the deviatoric strains follow:

$$ \left\{ \begin{matrix} q \\ \epsilon_d \end{matrix} \right\} = \left\{ \begin{matrix} mp' \frac{\lambda'}{\lambda' - \kappa'} \ln \frac{p_0}{p'} \\ \frac{\kappa' \lambda' - \kappa'}{\lambda'} m \ln \frac{M}{M - \eta} \end{matrix} \right\} \quad (53) $$

It is apparent that while the stress ratio monotonically increases with increasing deviatoric strains until the $\eta=M$ condition is asymptotically achieved, the deviator stress reaches a peak and then decreases, as shown in Figure 10. This occurs for:

$$ \eta = \eta_I = \frac{\lambda'}{\lambda' - \kappa'} m \quad (54) $$

Lade [34] introduced the concept of instability line, whose expression is given by Equation (54). Lade showed in fact that, when Equation (54) is fulfilled, an undrained test under load control becomes unstable. This is also clearly predicted by the theory of controllability. If $\eta=\eta_I$, the term of the compliance matrix connecting volumetric strain to isotropic effective stress rates, $C_{pp}$ in Equation (51) is zero. Since the increment of the deviator stress is also zero, it is readily apparent from Equation (51) that in undrained conditions the isotropic effective stress rate can take any (negative) value. The test becomes therefore uncontrollable. In particular, since the total stress is constant, a variation of the isotropic effective stress implies a corresponding (positive) variation of the pore water pressure. This phenomenon, the uncontrolled generation of pore water pressure under constant loading, was experimentally verified [16] and it is clearly associated to an instability condition.

The concept of loss of controllability is more general, however, since the instability in undrained tests is only a special case of possible loss of controllability, as we will see in a clearer way later
<figure><img src="image_11.png" /><figcaption>Figure 10. Predictions of an undrained test: (a) deviatoric stress-strain relationship; (b) stress ratio vs strains $M=1.2$; $m=0.4$; $\lambda'=0.09$; $\kappa'=0.03$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Theoretical predictions in a q constant stress controlled drained test: (a) stress path and (b) volumetric strains.</figcaption></figure>

on. Furthermore, other instability types are possible for lower values of the stress ratio, i.e. ‘before’ Equation (54) is fulfilled.

For the time being, it is interesting to note that in order to observe a peak in the deviator stress in the undrained test it is necessary (but not sufficient) that $m<M$, i.e. the flow rule is non-associated. This is a general result, as shown by Nova [41], and does not depend on the particular model we have chosen here to interpret soil behaviour.

It is furthermore very important to note that the stress state at the instability line is associated to a positive hardening modulus. From Equations (50) and (54) we derive in fact that:

$$H_I = \frac{mp'}{\lambda' - \kappa'}(M - \eta_I) = \frac{mp'}{\lambda' - \kappa'} \left( M - \frac{\lambda'}{\lambda' - \kappa'} m \right) \quad (55)$$

This means that the decrease in the deviator stress after peak is not associated to softening (negative hardening modulus), but it is only a result of the particular load path imposed on the specimen. Indeed, if we change the load programme, e.g. by allowing drainage after the peak has been reached,
the deviator stress increases again and the test is fully controllable up to $\eta=M$, as experimentally shown by di Prisco et al. [42].

Equation (55) is a condition of instability only if drainage is not allowed and load controlled. If we allow for drainage, even by taking constant the axial load and reducing the cell pressure, the test is fully controllable until the usual failure condition $\eta=M$ is achieved because it is fully stress controlled. Strains are in fact given by:

$$ \left\{ \begin{matrix} \varepsilon_v \\ \varepsilon_d \end{matrix} \right\} = \left\{ \begin{matrix} \frac{\lambda' - \kappa'}{m} (\eta - \eta_0) - \lambda' \ln \frac{\eta}{\eta_0} \\ \frac{\lambda' - \kappa'}{M} \left( \frac{M-m}{m} \ln \frac{M-\eta_0}{M-\eta} - \ln \frac{\eta}{\eta_0} \right) \end{matrix} \right\} \quad (56) $$

It is possible to note in Figure 11(b) that while deviatoric strains increase monotonically, volumetric strains are first negative then reach a peak (for $\eta=\eta_I$) (see also Figure 3) and eventually increase up to reaching an horizontal asymptote for $\eta=M$.

A further observation concerns the sign of the second-order work. Since the constitutive law is valid no matter which are the parameters controlled during the test, the value of the second-order work is determined only by the stress state provided no loss of controllability condition is achieved for a given stress path. For instance, at $\eta=\eta_I$, the second-order work is zero for the direction $q$ constant. An instability occurs only if the control is switched, e.g. to undrained conditions, however, while the test is fully controllable if drainage is allowed for, despite the second-order work at that value of the stress ratio is zero in both cases. The nullity of the second-order work is therefore a necessary but not sufficient condition for the occurrence of instability behaviour. This can emerge only if the appropriate variables are controlled.

It is finally possible to determine at what stress level the second-order work may first vanish. As already stated in Section 3, this condition is associated with the nullity of the determinant of the symmetric part of the stiffness matrix, $D_s$. As it can be easily proven (see e.g. [5]) this vanishes at the same time at which the determinant of the symmetric part of the compliance matrix $C_s$ vanishes. From Equation (56) this occurs when:

$$ \frac{1}{p'^2} \left\{ \left[ k' + \frac{(\lambda' - k')(m - \eta)}{m} \right] \left[ \frac{\lambda' - k'}{m(M - \eta)} \right] - \left[ \frac{\lambda' - k'}{2m} \right]^2 \left[ 1 + \frac{m - \eta}{M - \eta} \right] \right\} = 0 \quad (57) $$

The condition described in Equation (57) is fulfilled when:

$$ \eta = \eta_H = M - \frac{\lambda' - \kappa'}{m\kappa'} \left(\frac{M-m}{2}\right)^2 \quad (58) $$

It is obvious that $\eta_H < M$, if the flow rule is non-associated, i.e. the Hill instability condition occurs for positive hardening modulus. It is also easy to show that $\eta_H < \eta_I$ signifies that the so-called instability line is not the locus of the lowest instability condition, a result valid for any elastoplastic model, as shown by Nova [5]. Hill's instability stress ratio can be either lower or greater than $m$, instead, depending on the relative values of $\lambda'$ and $\kappa'$ on the one hand and of $m$ and $M$ on the other.

For $\eta=\eta_H$, the direction along which the second-order work is zero, no matter what are the control variables, is:

$$ dq = \left( \frac{M-m}{2} \right) \left[ 1 - \frac{\lambda' - \kappa'}{\kappa'} \left( \frac{M-m}{2m} \right) \right] dp' = -\alpha dq' \quad (59) $$

Loss of controllability occurs at this stress level whenever we pretend to control the following conjugate generalized stress and strain variables:

$$ \left\{ \begin{matrix} d\eta_2 \\ d\xi_1 \end{matrix} \right\} = \left\{ \begin{matrix} d\varepsilon_v - \alpha d\varepsilon_d \\ \alpha d p' + dq \end{matrix} \right\} = \left\{ \begin{matrix} 0 \\ d\eta_1 \end{matrix} \right\} \quad (60) $$
<figure><img src="image_13.png" /><figcaption>Figure 12. Various instability conditions and region of latent instability.</figcaption></figure>

It is readily apparent that the loading programme defined by Equation (60) is a special case of the R test defined by Darve et al. [23] (see Section 2.2, Equation (4)) with $\alpha=3(R-1)/(2R+1)$.

Figure 12 summarizes the findings of this section. The various instability loci, all straight lines crossing the origin of axes, are plotted together with the failure locus. It is apparent that in the region between $\eta=\eta_H$ and $\eta=M$ it is possible to find a particular load path for which the controllability is lost. The entire region can then be considered as a region of latent instability, which can be put in evidence whenever the appropriate control variables are used.

### 4.1.2. Density, stress and fabric-dependent elasto-plastic model.
We herein turn to an elastoplastic constitutive model that embeds density, stress and fabric level dependencies in order to study diffuse failure in undrained conditions and drained $q$-constant tests under axi-symmetrical conditions as presented earlier in Section 2.

The model is an offshoot of a micro-mechanical analysis of an REV of a granular material undergoing deformation from which a fabric-embedded dilatancy emerges. Apart from fabric dependencies, this dilatancy law is also a function of stress and density levels, pivotal characteristics of granular materials. This so-called stress-density-fabric dilatancy equation enters readily into a classical elasto-plastic model as a flow rule for describing deviatoric material response. Additionally, isotropic (compactive) volumetric response is described by a second flow rule that operates on a cap surface. Hardening rules for both yield surfaces to describe the development of plastic deformations as well as an evolution law to describe the build up of texture (fabric) complete the description of the model. For the sake of simplicity, the model will be presented in the axi-symmetric setting.

We introduce two yield (plastic loading) surfaces $f_c$ and $f_s$ which refer to the isotropic and deviatoric mechanisms, respectively. Thus, in the $p'-q$ plot, expressions of the surfaces are simply:

$$ (\text{cap}) f_c(p', p_c) = p' - p_c \qquad (61) $$

and

$$ (\text{shear}) f_s(p', q, \phi_m) = q - M_{\phi_m} p', \quad M_{\phi_m} = \frac{6 \sin \phi_m}{3 - \sin \phi_m} \qquad (62) $$

where $p_c = p_c(\epsilon_v^p)$ is the pre-consolidation pressure and $\phi_m$ is the mobilized friction angle. It is noted that both $p_c$ and $\phi_m$ are internal variables that control the size of the surfaces during
hardening according to some evolution law that will be described later on. As such, $f_s$ collapses into the limit plastic failure surface $f_f$ at critical state where the friction angle is $\phi_{cv}$.

The plastic flow rule extended to multi-surface plasticity is given as

$$d\epsilon^p = \sum_i \Lambda_i \frac{\partial g_i(\boldsymbol{\sigma}, \zeta)}{\partial \boldsymbol{\sigma}} \qquad (63)$$

where $g$ is a plastic potential function, $\zeta$ is a hardening parameter, $\Lambda$ a positive plastic multiplier, and $i=c, s$ depending on the active plastic mechanism.

The isotropic mechanism referring to plastic compaction is chosen to be simply:

$$g_c(p, \epsilon_v^p) = f_c = p - p_c(\epsilon_v^p) \qquad (64)$$

whereas the deviatoric plastic mechanism is described by the plastic potential $g_s$ as:

$$g_s(p, q, \psi_m) = q - M_{\psi_m} p, \quad M_{\psi_m} = \sin \psi_m \qquad (65)$$

with $\psi_m$ the mobilized dilatancy angle derived from a stress-density-fabric dilatancy equation:

$$\sin \psi_m = \frac{4}{3} \left( \frac{\sin \varphi_m - \sin \varphi_f}{1 - \sin \varphi_m \sin \varphi_f} \right), \quad \sin \varphi_f = \frac{b F_{33} / F_{11} + \gamma^p}{c + \gamma^p} \left( \frac{e}{e_{cr}} \right)^\alpha \sin \varphi_{cv} \qquad (66)$$

where $e$ and $e_{cr}$ are current and critical void ratios, respectively, whereas $b, c$ and $\alpha$ are constants. Fabric information is transmitted through fabric components $F_{11}$ and $F_{33}$ which are projections of the second-order fabric tensor $F_{ij}$ [43, 44] along principal stress directions $\sigma_1$ and $\sigma_3$, respectively. It is evident that both positive and negative rates of dilation can be obtained depending on the relative magnitudes of $\varphi_m$ and $\varphi_f$. In the limit, fabric conditions can be such that a positive rate of dilation is possible even though the current void ratio is higher than critical, i.e. $e > e_{cr}$.

The consolidation (cap) surface $f_c$ evolves according to variations in the consolidation pressure $p_c$ such that:

$$d p_c = \frac{\kappa}{h_c} \left( \frac{\epsilon_{v(c)}^p}{h_c} \right)^{\kappa-1} d\epsilon_{v(c)}^p \qquad (67)$$

where $\kappa$ and $h_c$ are constants and $\epsilon_{v(c)}^p$ is the plastic volumetric strain due to consolidation.

On the other hand, the deviatoric hardening is described by the evolution of $f_s$ such that $\zeta = \varphi_m(\epsilon_d^p)$, i.e.

$$\sin \varphi_m = \frac{\epsilon_d^{p*}}{a + \epsilon_d^{p*}} \left( \frac{e}{e_{cr}} \right)^{-\beta} \sin \varphi_{cv} \qquad (68)$$

where $a$ and $\beta$ are constants and $\epsilon_d^{p*}$ is the nominal plastic deviatoric strain which involves fabric, i.e. $\epsilon_d^{p*} = (2/3 \mathbf{e}^{p*} \cdot \mathbf{e}^{p*})^{1/2}$ with $\mathbf{e}^{p*}$ is the deviatoric part of $\boldsymbol{\varepsilon}^* = \boldsymbol{\varepsilon} \cdot \mathbf{F}^{-1}$ (for $\mathbf{F}$, see Equation (70)). The critical void ratio in turn varies with mean stress $p^*$:

$$e_{cr} = e_{cr0} \exp \left[ - \left( \frac{p^*}{h_{cr}} \right)^{n_{cr}} \right] \qquad (69)$$

where $e_{cr0}$, $n_{cr}$ and $h_{cr}$ are constants, whereas the dependence of fabric is captured by $p^* = \text{trace}(\boldsymbol{\sigma}^*)/3$. As such, the critical state line in the $e$-log $p$ space is non-unique in spite of $\eta$ being constant.

To complete the model, some evolution of fabric during deformation history is prescribed. It is convenient to make fabric change with the level of deviatoric stress level as observed in experiments [45]. Thus, the simplest functional form of this evolution law is given as:

$$\dot{\mathbf{F}} = \chi \dot{\eta}, \quad \eta = s/p \qquad (70)$$

which points to a zero fabric change at critical state ($\eta = \text{constant}$).
<figure><img src="image_14.png" /><figcaption>Figure 13. Numerical results of CSD tests using: (a) direct method through stress-controlled loading programme and (b) forced dilation method through mixed-controlled loading programme. For a fully stress-controlled test (a) the plastic limit condition is reached while for a mixed loading condition (b) it is not possible to pursue the numerical test.</figcaption></figure>

We herein consider a CSD test (see Section 2.1). In line with what was discussed earlier in Section 2.1, we will introduce two ways to carry out numerically such a test following two loading programmes with different control parameters: (a) stress controlled, $d_p'<0$ with $dq=0$ and (b) mixed controlled, $dq=0$ with water injection or volume increase $d\epsilon_v<0$.

The initial mean stress level in all numerical tests is 300 kPa, whereas three distinct levels of deviatoric stresses are investigated, i.e. 50, 100 and 200 kPa. The material studied is a loose sand with an initial void ratio equal to 0.80. Figure 13 shows the numerical simulations of the above tests by integrating the constitutive equations.

The results for the $q$ constant tests under stress-controlled mode are presented in Figure 13(a). In Figure 13(b), we display similar plots for the very same $q$ constant tests simulated under mixed-controlled mode.

The second-order work is computed all along the complete loading paths. Since $dq=0$ and $d_p'<0$ throughout the CSD test, a vanishing second-order work is signalled whenever the volumetric strain reaches a peak ($d\epsilon_v/dp'=0$ since $d_2W=d_p'd\epsilon_v=0$). This condition is indeed verified in the numerical simulations as indicated by small arrows pointing down in Figure 13(a).

Figure 13(a) also shows that the test can proceed past the bifurcation line because the stress-controlled mode allows it. The numerical test eventually stops whenever the effective stress path reaches the plastic limit surface giving way to a localization failure mode. On the other hand, we reveal in Figure 13(b) that when the same test is conducted in mixed loading mode, the solution breaks down as soon as effective stresses reach the bifurcation line previously defined under stress-controlled mode. This is because, in addition to the vanishing of the second-order work, the control variables which are of mixed mode are such that collapse with diffuse failure emerges. It is also interesting to notice that this occurs well below the plastic limit surface. Along with the lines of Section 3.2, the vanishing of the second-order work is a necessary but not sufficient condition for effective diffuse failure occurrence. The proper control variables are therefore also necessary.
<table><caption>Table II. Numerical simulation carried out using the micro-directional model: constitutive parameters and initial conditions.</caption><thead><tr><th>Initial isotropic stress (kPa)</th><th>Initial void ratio</th><th>k<sub>n</sub> (kN m)</th><th>k<sub>t</sub> (kN m)</th><th>φ<sub>μ</sub> (°)</th></tr></thead><tbody><tr><td>125</td><td>0.66</td><td>15 708</td><td>6830</td><td>15</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>Figure 14. Stress-controlled *q* constant test preceded by conventional undrained test up to peak point.</figcaption></figure>

For a further exploration of instability, we performed another series of numerical experiments similar to the CSD tests described earlier, except that the step which corresponds to the drained conventional triaxial loading is now replaced with a CU loading until to the peak deviatoric stress is reached. At this peak, the CSD path is pursued as before. The same material parameters listed in Table II are used. Figure 14 summarizes a few interesting findings. With regards to the CSD phase, we point out that the simulation could be controlled throughout the entire test with stress control variables despite the non-positivity of second-order work from the start of the simulation. This is essentially because det **D** is never zero. However, we see in Figure 14(a) that det (**D**) approaches zero as the effective stress path reaches the plastic limit line. On the other hand, when the forced dilation procedure (with mixed control variables) was attempted to pursue the constant drained phase after the CU phase, the numerical simulations broke down right away. No solution for the *q*-constant path could be obtained regardless of the size of the control parameter increment. These numerical exercises are useful in that they reinforce the main conclusions arrived at in Darve *et al.* [17]. In other words, there certainly exist other failure surfaces within the plastic limit one when examining non-associated materials like geomaterials. In fact, there is a whole domain in axisymmetric stress–strain space where bifurcation, non-uniqueness, instabilities and failure appear and that depends on a series of factors, in special, the stress–strain history, current loading direction and loading mode.

Finally, dependencies of material instability on initial void ratio and fabric are depicted in Figures 15(a), (b), respectively. When looking at the sole effect of initial void ratio, it is found that material instability is encountered earlier for a loose sand than for a dense one as shown in Figure 15(a). Hence, a larger instability domain is expected for a loose sand. Given the importance of fabric, the stability analysis is next repeated by assigning an initial fabric to the specimen and letting it evolve accordingly with deformation history. A fabric tensor **F** was introduced far back in the beginning of this subsection; refer to Equation (70), to describe both inherent and induced anisotropies. Two initial fabric cases are examined for the same initial void ratio: one where a greater number of contact normals is oriented vertically giving rise to a strong fabric in
<figure><img src="image_16.png" /><figcaption>Figure 15. Effect of fabric on the instability response—enlarged bifurcation zone: (a) density level without fabric and (b) consideration of fabric.</figcaption></figure>

the horizontal plane ($\theta=0$) and the other extreme case where the fabric is strong in a vertical plane ($\theta=90^{\circ}$). To reveal the influence of fabric on material instability, Figure 15(b) compares the points of bifurcation at peaks of $\epsilon_v$ along the $q$ constant paths for three different deviatoric stress levels with and without fabric, but at the same initial void ratios of $e_0=0.7$ and 0.8. In the presence of an initial fabric (anisotropy), bifurcation is encountered earlier than when fabric is not considered. As seen in Figure 15(b), the unstable region is now enlarged for the case of the initial void ratio $e_0=0.8$. The anisotropy of the material, as dictated by the grain arrangements, invariably influences the material instability which is herein captured by the model.

### 4.1.3. Incrementally non-linear model of second order.
We next turn to a different class of phenomenological constitutive relations which is the incrementally non-linear ones. A major distinction from the classical elasto-plasticity theory is that this class of constitutive models does not need any assumption to be made on:

1. strain decomposition into an elastic and a plastic part;

2. existence of an elastic limit; and

3. existence of a flow rule.

However, the notion of tensorial zone will be used as was first introduced by Darve and Labanieh [46]. A tensorial zone is a domain of the loading space in which the incremental constitutive relation is linear. From this viewpoint, a classical elasto-plastic model has two tensorial zones: one for loading and the second one for unloading. Beside the introduction of tensorial zones, incrementally non-linear models have been built from the general expression of rate-independent tensorial relations [46, 47] relating stress to strain increments ($d\sigma$ and $d\epsilon$, respectively) as follows:

$$d\epsilon = D_h(u) d\sigma \quad (71)$$

where $D_h$ is a nonlinear operator depending on the previous stress–strain history through state variables and memory parameters $h$ and $u$ is the unit vector in the direction of $d\sigma$: $u = d\sigma / \|d\sigma\|$. An expansion of $D_h(u)$ limited to the second-order terms gives the general expression of incrementally nonlinear constitutive relations of second order, i.e.

$$d\epsilon_{\alpha} = D_{\alpha\beta}^{1} d\sigma_{\beta} + \frac{1}{\|d\sigma\|} D_{\alpha\beta\gamma}^{2} d\sigma_{\beta} d\sigma_{\gamma} \quad (\alpha, \beta, \gamma = 1, \dots, 6) \quad (72)$$
Finally, in order to obtain Darve's constitutive relation, three more hypotheses have to be added:

1. **D** is orthotropic,

2. $d\sigma_{\beta}d\sigma_{\gamma}=0$ for $\beta \neq \gamma$,

3. shear moduli are incrementally linear.

Thus, in principal axes to simplify in this paper, we obtain:

$$ \begin{pmatrix} d\epsilon_1 \\ d\epsilon_2 \\ d\epsilon_3 \end{pmatrix} = \frac{1}{2} [\mathbf{C}^+ + \mathbf{C}^-] \begin{pmatrix} d\sigma_1 \\ d\sigma_2 \\ d\sigma_3 \end{pmatrix} + \frac{1}{2\|d\sigma\|} [\mathbf{C}^+ - \mathbf{C}^-] \begin{pmatrix} (d\sigma_1)^2 \\ (d\sigma_2)^2 \\ (d\sigma_3)^2 \end{pmatrix} \quad (73) $$

where

$$ \mathbf{C}^{\pm} = \begin{bmatrix} \frac{1}{E_1^{\pm}} & -\frac{v_{12}^{\pm}}{E_2^{\pm}} & -\frac{v_{31}^{\pm}}{E_3^{\pm}} \\ -\frac{v_{21}^{\pm}}{E_1^{\pm}} & \frac{1}{E_2^{\pm}} & -\frac{v_{32}^{\pm}}{E_3^{\pm}} \\ -\frac{v_{13}^{\pm}}{E_1^{\pm}} & -\frac{v_{23}^{\pm}}{E_2^{\pm}} & \frac{1}{E_3^{\pm}} \end{bmatrix} \quad (74) $$

Coefficients $E_i^+$ and $v_{ij}^+$ are defined along generalized triaxial loading paths when $d\sigma_i > 0$, whereas $E_i^-$ and $v_{ij}^-$ refer to the same paths when $d\sigma_i < 0$ ($i = 1, 2, 3$). For more details regarding the model, in particular the variation of tangent moduli and Poisson's ratio with stress-strain history, the reader is referred to Darve et al. [47]. In the special one-dimensional case, the relationship in Equation (73) is incrementally piecewise linear (one modulus for loading and another one for unloading). Thus:

$$ d\epsilon = \frac{1}{2} \left( \frac{1}{E^+} + \frac{1}{E^-} \right) d\sigma + \frac{1}{2} \left( \frac{1}{E^+} - \frac{1}{E^-} \right) |d\sigma| \quad (75) $$

By extrapolation, Darve defined an octo-linear relation (eight tensorial zones), that is incrementally piecewise linear. Using previously defined notations, the octo-linear relation is written as follows:

$$ \begin{pmatrix} d\epsilon_1 \\ d\epsilon_2 \\ d\epsilon_3 \end{pmatrix} = \frac{1}{2} [\mathbf{C}^+ + \mathbf{C}^-] \begin{pmatrix} d\sigma_1 \\ d\sigma_2 \\ d\sigma_3 \end{pmatrix} + \frac{1}{2} [\mathbf{C}^+ - \mathbf{C}^-] \begin{pmatrix} |d\sigma_1| \\ |d\sigma_2| \\ |d\sigma_3| \end{pmatrix} \quad (76) $$

which shows explicitly eight tensorial zones for which there are eight linear relations, i.e.

$$ d\epsilon = (C_i)_{i=1,...,8} d\sigma \quad (77) $$

where $(C_i)_{i=1,...,8}$ refer to the matrix **C** with indices (+) are affected to the column ($j$) if $d\sigma_j > 0$, and ($-$) if $d\sigma_j < 0$ ($j = 1, 2, 3$). As for the case of $d\sigma_i = 0$, it can be verified that the relationship established in Equation (73) is continuous [36].

The attractive mathematical simplicity of the octo-linear relation is herein used to explicitly determine Hill's bifurcation domain as well as to derive the analytical equation of the cones containing all directions for which the material is unstable. In the following, all analytical developments are performed using the octo-linear model defined in Equation (76). For the sake of simplicity in the notations, coefficients $E_i^+$, $E_i^-$ and $v_{ij}^+$, $v_{ij}^-$ are not distinguished and are thus written $E_i$ and $v_{ij}$ as analytical developments are assumed to be performed in the proper given tensorial zone. Numerical results are given for both models (i.e. octo-linear and incrementally non-linear) presented above. The analytical developments performed with the octo-linear are general enough to remain true for
every incrementally piecewise linear model. Only the number of tensorial zones has to be adapted (two for classical elasto-plastic relations instead of eight in the present case).

Let us recall **D** as the constitutive matrix that relates dσ to dε, whereas conversely **C** links dε to dσ. If **D** is invertible then **C**=$\mathbf{D}^{-1}$. By splitting matrices **D** and **C** in symmetrical parts **D**<sub>s</sub> and **C**<sub>s</sub> as well as skew symmetrical parts, the second-order work takes the following expression:

$$d_2 W = d\sigma d\varepsilon = d\sigma C_s d\sigma = d\varepsilon D_s d\varepsilon \quad (78)$$

Thus, the sign of the second-order work is clearly linked to the positive definiteness of **C**<sub>s</sub> or **D**<sub>s</sub>. In fact, the positivity of d<sub>2</sub>W and hence of the quadratic forms dσ**C**<sub>s</sub>dσ and dε**C**<sub>s</sub>dε has a geometrical representation in the stress rate and strain rate spaces, respectively. Let us now analyse the condition:

$$d\sigma d\varepsilon = d\sigma C_s d\sigma = 0 \quad (79)$$

Substituting Equation (76) into Equation (79) leads to:

$$\frac{d\sigma_1^2}{E_1} + \frac{d\sigma_2^2}{E_2} + \frac{d\sigma_3^2}{E_3} - \left(\frac{\nu_1^2}{E_1} + \frac{\nu_2^2}{E_2}\right) d\sigma_1 d\sigma_2 - \left(\frac{\nu_3^2}{E_3} + \frac{\nu_2^2}{E_3}\right) d\sigma_2 d\sigma_3 - \left(\frac{\nu_1^3}{E_1} + \frac{\nu_3^3}{E_3}\right) d\sigma_1 d\sigma_3 = 0 \quad (80)$$

which represents an elliptical cone in the stress rate space. Nevertheless, some degenerated form may exist depending on the positivity of the eigenvalues of **C**<sub>s</sub>. Figure 16 presents the four cases that have to be taken into account.

Thus, the directional feature of the second-order work criterion observed by Darve and Laouafa [48] and Laouafa and Darve [49] is exemplified here. It is precisely the condition of det(**C**<sub>s</sub>)≤0 or

<figure><img src="image_17.png" /><figcaption>Figure 16. Geometrical representation of the mathematical solutions of Equation (80). λ<sub>i</sub> are the eigenvalues of **C**<sub>s</sub>.</figcaption></figure>
conversely det($D_s$)≤0 that implies the existence of some loading directions for which bifurcation occurs. In fact, real solutions are the ones (or part of ones) which are fulfilled in the proper tensorial zone. If we call the limit of the bifurcation domain as the set of stress–strain states for which a unique unstable loading direction exists, then this limit can be described by the following condition by assuming that eigenvalues of ($C_s$)$_i$ are evolving continuously with the loading parameter:

$$ \min_{i=1,...,n} (\det(C_s)_i) = 0 \quad \text{with} \quad u_i \subset Z_i \qquad (81) $$

where n is the number of tensorial zones (of the constitutive model), $u_i$ the eigenvector corresponding to the vanishing eigenvalue and $Z_i$ the tensorial zone considered.

To conclude for analytical purposes, it can be proved that [50]:

$$ \det(D_s) = \frac{\det(C_s)}{(\det(C))^2} \qquad (82) $$

which further leads to

$$ \det(D_s) = 0 \Leftrightarrow \det(C_s) = 0 \qquad (83) $$

The consequence of this result is that the bifurcation analysis can be performed without any restriction on the formulation of the constitutive relation. However, this is not the case for the plastic limit state since the condition of det($D$)=0 is not equivalent to det($C$)=0. In the case of an associated material, relation (82) is still valid, but not condition (83), because in this latter case the bifurcation limit coincides with the plasticity limit condition.

The limits of the bifurcation domain given by Equation (83) are illustrated in Figure 17 for both models (octo-linear and non-linear) proposed by Darve. It is noted that the incrementally non-linear model can be seen as incrementally piecewise linear with an infinite number of tensorial zones. Thus, with some numerical developments, an approximation of the bifurcation limit has been calculated [51]. Finally, Figure 18 shows some instability cones as calculated by both models for stress–strain states located beyond the bifurcation limit.

<figure><img src="image_18.png" /><figcaption>Figure 17. Limit of the bifurcation domain plotted in the 3D stress space for constitutive models of Darve.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18. 3D cones of unstable stress directions for a dense sand of Hostun. On the top the cones obtained with the octo-linear model. The planes represent the limit between the eight tensorial zones, the meshes correspond to the analytical solution of Equation (80). On the bottom results obtained with the non-linear model using the numerical method. The point clouds represent the solution plotted on the sphere. $p_o$ is the initial confining pressure, $q = \sigma_1 - \sigma_3$.</figcaption></figure>

## 4.2. Micromechanical approaches

Micro-mechanical approaches offer a paradigm shift in reasoning from classical phenomenological framework. They aim at deriving a constitutive relation between both strain and stress on the macroscopic scale by incorporating information down to micro (a pair of grains) or mesoscopic (cluster of grains) scales. Typically, the mesoscopic scale is characterized by force chains known to play a fundamental role in the stability of granular assemblies [52, 53].

The micromechanical models presented in this section use two well-known schemes to transfer from one scale to another, namely, homogenization and projection (also called localization). Figure 19 shows two equivalent methods that can be used to formulate a micromechanical model. In the micro-directional model ($\mu$-D model) [54] the transfer of scales operates from macroscopic strain tensor to infer particle relative displacements via localization, then on to contact forces through the local contact law, and finally arrive at the macroscopic stress tensor using homogenization. The CH model proposed by Chang and Hicher [55] uses the reverse order starting from the stress tensor to arrive at the strain tensor using the same techniques as in the above.

The two models will be presented in what follows.

### 4.2.1. Micro-directional model ($\mu$-D model)
Throughout this sub-section, time and spatial derivatives of any variable $\psi$ will be distinguished by denoting $d\psi$ the time derivative of $\psi$ (defined as the product of the particulate derivative $\dot{\psi}$ by the infinitesimal time increment $dt$) with respect to a given frame, and by denoting $\delta\psi$ the spatial derivative of $\psi$.

The micro-directional model was initially developed to describe the mechanical behaviour of snow [56–58]. Thus, thanks to the ability of the model to incorporate any local constitutive relation at the microscopic scale, the model was generalized to any type of granular assembly, with a particular emphasis on frictional granular materials [54].

The micro-directional model, which is formulated in a small-strain Eulerian formulation, allows the Cauchy stress tensor $\sigma$ to be related to the small strain tensor $\epsilon$ considering micro-mechanical characteristics. Changes in the fabric of the granular assembly are described by directly modelling the increase or the decrease in the number of contacts along each direction of the physical space.
<figure><img src="image_20.png" /><figcaption>Figure 19. General homogenization/localization scheme relating the incremental stress tensor and the incremental strain tensor.</figcaption></figure>

In this approach, the number of contacts along a given direction is not computed from any fabric tensor, but is related to the normal strain rate along that direction. Following pioneering work based on physical evidence [59], it is observed that the number of contacts increases along contractive directions, whereas it decreases along dilative directions [60]. Therefore, the distribution of contacts is likely to evolve over a general loading path, inducing anisotropy to the fabric. The three basic stages of the homogenization/localization procedure (Figure 19) are very briefly reviewed here (for more details, see [54]):

1. The stress averaging corresponds to the Love-Weber formula [61–64]:

$$ \sigma = \frac{1}{V} \sum_{c=1}^{N_c} \mathbf{f}^c \otimes \ell^c, \qquad (84) $$

where $\ell^c$ is the branch vector joining the centers of particles in contact at $c$, $\mathbf{f}_c$ is the contact force, and the sum is extended to all the $N_c$ contacts occurring in the RVE of volume $V$. The norm of the branch vector $\ell^c$ is assumed to be a constant parameter (equal to the average mean diameter of the grains) whose evolution during loading history is ignored. This ensures that the terms $\mathbf{f}^c$ and $\ell^c$ are uncorrelated. The discrete summation given in Equation (84) can be replaced with a continuous integration over the domain $\Omega$ containing all the contact directions in the physical space, and thus providing the directional character to the model:

$$ \sigma = 2r_g \int_{\Omega} \hat{\mathbf{f}} \otimes \mathbf{n} \omega(\mathbf{n}) d\Omega, \qquad (85) $$

where $\omega(\mathbf{n})$ is the density of contacts along each space direction $\mathbf{n}$, $r_g$ denotes the mean (equivalent) radius of the grains, $\hat{\mathbf{f}}$ is the average of all contact forces $\mathbf{f}^c$ associated with contacts oriented in the direction $\mathbf{n}$, and $d\Omega$ is the elementary solid angle. After differentiation it follows that:

$$ d\sigma = 2r_g \left[ \int_{\Omega} d\hat{\mathbf{f}} \otimes \mathbf{n} \omega d\Omega + \int_{\Omega} \hat{\mathbf{f}} \otimes \mathbf{n} d\omega d\Omega \right]. \qquad (86) $$

2. The kinematical projection rule to express particle displacements (local) in terms of strains (global) is given by:

$$ d\hat{\mathbf{u}}(\mathbf{n}) = 2r_g d\varepsilon \cdot \mathbf{n}, \qquad (87) $$

where $\hat{\mathbf{u}}(\mathbf{n})$ is the directional kinematic variable linked to $\hat{\mathbf{f}}(\mathbf{n})$ along the contact direction $\mathbf{n}$. It is cautioned here that a linear approximation (affine transformation) has been assumed and this may be invalid along dilative directions based on discrete element simulations (see [65]).
3. The local behaviour is described using a standard elastoplastic model relating both the local normal force $f_n^c$ and the local tangential force $f_t^c$ (contained in the tangential plane) to both the local normal relative displacement $u_n^c$ and the local tangential relative displacement $u_t^c$. This constitutive relation can be written in an incremental form as:

$$d f_n^c = k_n d u_n^c \qquad (88)$$

$$d \mathbf{f}_t^c = \vartheta \frac{\mathbf{f}_t^c + k_t \mathbf{d} \mathbf{u}_t^c}{\| \mathbf{f}_t^c + k_t \mathbf{d} \mathbf{u}_t^c \|} - \mathbf{f}_t^c \qquad (89)$$

$$\vartheta = \min[\|\mathbf{f}_t^c + k_t \mathbf{d} \mathbf{u}_t^c\|, \tan \varphi_{\mu}(\mathbf{f}_n^c + k_n \mathbf{d} \mathbf{u}_n^c)]$$

It will be assumed hereafter that Equations (88) and (89) also apply to the directional average variables $\hat{\mathbf{u}}(\mathbf{n})$ and $\hat{\mathbf{f}}(\mathbf{n})$. The contact model based on the Mohr–Coulomb criterion introduces only three parameters: two elastic normal ($k_n$) and tangential ($k_t$) stiffnesses, both constant, and a local friction angle $\varphi_{\mu}$. Other constitutive relations can be envisaged, depending on the material considered. For instance, in the case of snow, a viscoplastic relation was introduced to describe mechanical interactions between ice grains [56, 57].

Finally, it is worth noting that the model can be extended to account for local physical phenomena such as suction due to liquid bridges, cementation due to solid bonds, or time dependence due to viscous behaviour.

Let us now investigate whether the $\mu-D$ model can reveal the occurrence of bifurcation points along any arbitrary loading path. As such, a triaxial loading path is first simulated in drained axisymmetric conditions until various deviatoric stress ratios $\eta = q/p'$ are reached followed by the application of directional stress probes in space according to Gudehus [36]. Incremental stress vectors $d\sigma$ are imposed along all the directions of the incremental stress space and the corresponding incremental strain responses $d\varepsilon$ are computed. Then, for each loading direction, the resulting value of the normalized second-order work $d_2\tilde{W} = d_2W/(\|d\sigma\|\|d\varepsilon\|)$ can be evaluated.

To check whether stress loading directions exist for which the second-order work vanishes, a circular representation is adopted. Basically, for each direction of the loading $\alpha_{\sigma} = \tan^{-1}[d\sigma_1/(\sqrt{2}d\sigma_3)]$, the corresponding points $(\cos\alpha_{\sigma}(d_2W+c), \sin\alpha_{\sigma}(d_2W+c))$ are plotted, where c is a given constant positive parameter. Whenever the second-order work is negative, the corresponding point is located inside the circle of radius c.

The simulations considered herein were run with the set of parameters reported in Table II.

As seen in Figure 20, an instability cone, enclosing all stress directions along which the second-order work is taking negative values, opens up when $\eta$ increases beyond a critical value $\eta_c \approx 2.070$. When $\eta \approx \eta_c$, the cone degenerates more or less into a single direction, which means that the corresponding stress–strain state belongs to the surface of the bifurcation domain. As the stress ratio $\eta$ gradually increases (for $\eta=2.070$ and 2.112), an instability cone clearly develops showing that we are strictly inside the bifurcation domain. In conclusion, at least for the loading path considered, a bifurcation can be detected by using the $\mu-D$ model.

Next, the same issue is explored for an isochoric loading path for which the volume of the sample is preserved while an axial strain rate is imposed. For loose specimens, as shown experimentally in the beginning of this paper, the $q-p'$ curve passes through a peak. Under similar loading conditions, this result is essentially captured by the $\mu-D$ model with the parameters reported in Table II as shown in Figure 21, where this loading path was simulated. At this $q$ peak and isochoric conditions, the second-order work vanishes since $d_2W = dq d\varepsilon_1 = 0$ and $d\varepsilon_1 > 0$ (loading condition).

The corresponding stress–strain state belongs therefore to the bifurcation domain. To know whether this state is located on the surface or strictly within the bifurcation domain, a directional analysis can be carried out. For convenience, incremental strain vectors $d\varepsilon$ are imposed along all the directions of the incremental strain space, and the incremental stress response $d\sigma$ is computed. For each loading direction, the resulting value of the normalized second-order work as previously defined is evaluated and plotted as before.
<figure><img src="image_21.png" /><figcaption>Figure 20. Circular representation of the normalized second-order work using the micro-directional model. The second-order work was computed after an initial axisymmetric drained triaxial loading, for different values of the deviatoric stress ratios. It appears that beyond a critical value of the deviatoric stress ratio, a cone develops gathering directions along which the second-order work is negative or nil.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 21. Simulation of the response of a loose sample along an undrained triaxial loading path in axisymmetric conditions, using the micro-directional model. The *p*–*q* curve passes through a peak.</figcaption></figure>

As seen in Figure 22, an unstable cone that is not reduced to a single direction is recovered. The cone is small, but finite according to the numerical results. This result shows that the deviatoric peak in the *p'*–*q* plane belongs strictly to the bifurcation domain.

The so-called 'Lade's instability line' is therefore also located strictly inside the bifurcation domain.

4.2.2. *CH model*. In this model, we consider a granular material as a collection of particles. The deformation of a representative volume of the material is generated by the mobilization of contact
<figure><img src="image_23.png" /><figcaption>Figure 22. Circular representation of the normalized second-order work using the micro-directional model. The directional analysis was performed exactly at the q peak. A cone gathering directions along which the second-order work is negative or nil is visible even if small.</figcaption></figure>

particles in all orientations. Thus, the stress–strain relationship can be derived as an average of the mobilization behaviour of local contact planes in all orientations. During the integration process, (Figure 19) a relationship is required to link the macro and micro variables. Using the static hypotheses proposed by Liao et al. [66], we obtain the relation between the macro strain and inter-particle displacement:

$$ \dot{u}_{j,i} = A_{ik}^{-1} \sum_{\alpha=1}^{N} \dot{\delta}_{j}^{\alpha} l_{k}^{\alpha} \qquad (90) $$

where $\dot{\delta}_j$ is the relative displacement between two contact particles and the branch vector $l_k$ is the vector joining the centres of two contact particles. It is noted that contact particles include both direct contact and indirect contact of neighbouring particles associated with a Voronoi polyhedron as discussed by Cambou et al. [65]. For convenience, we let $N$ be the total number of contact orientations. The variables $\dot{\delta}_j^\alpha$ and $l_k^\alpha$ are defined, respectively, as the averaged values of $\dot{\delta}_j$ and $l_k$ for all contacts belonging to the $\alpha$th orientation. The fabric tensor in Equation (90) is defined as:

$$ A_{ik} = \sum_{\alpha=1}^{N} l_{i}^{\alpha} l_{k}^{\alpha} \qquad (91) $$

If the magnitudes of the branch vectors $l_i^\alpha$ are different for each orientation $\alpha$, then the fabric tensor $A_{ik}$ is anisotropic. During the deformation of the assembly, the particles are rearranged and the lengths of the branch vectors are changed into a new set of numbers, which is a way to allow the evolution of the strain-induced anisotropy to be considered. One of the major advantages of the model is its capacity to take into account the initial and induced structural anisotropy.

Using both the principle of energy balance and Equation (90), the mean force on the contact plane of each orientation is:

$$ \dot{f}_{j}^{\alpha}=\dot{\sigma}_{i j} A_{i k}^{-1} l_{k}^{\alpha} V \qquad (92) $$

In Equation (92), the stress increment $\dot{\sigma}_{ij}$ can be obtained through the contact forces and branch vectors for contacts in all orientations. The summation of a quantity $F^\alpha$ over all contacts can be
replaced by an integral over all orientations, provided that the quantity can be expressed by a continuous function F(γ, β). In granular packing, this is always the case. Thus,

$$ \sum_{\alpha=1}^{N} F^{\alpha} = \frac{N}{2\pi} \int_{0}^{\pi/2} \int_{0}^{2\pi} E(\gamma, \beta) F(\gamma, \beta) \sin \gamma \, d\gamma \, d\beta \quad (93) $$

$N_c$ is the number of contact orientations. E(γ, β) is the density function of contacts. When the function F(γ, β)=1, the following consistency condition must be satisfied for the density function of contacts:

$$ 1 = \frac{1}{2\pi} \int_{0}^{\pi/2} \int_{0}^{2\pi} E(\gamma, \beta) \sin \gamma \, d\gamma \, d\beta \quad (94) $$

The density function can be obtained from a spherical harmonic expansion of the true distribution of contacts in discrete orientations (for more details concerning the numerical integration see [55]).

As in the (μ–D) model, we consider an elasto-plastic local law relating the local forces $f_i^\alpha$ and the local movements $\delta_i^\alpha$ for a contact plane in the αth orientation. Local forces and displacements can be denoted as follows: $f_i^\alpha = \{f_n^\alpha, f_s^\alpha, f_t^\alpha\}$ and $\delta_i^\alpha = \{\delta_n^\alpha, \delta_s^\alpha, \delta_t^\alpha\}$, where the subscripts n, s, and t represent the components in the three directions of the local coordinate system. The elastic part is derived from Hertz–Mindlin’s formulation [67]. For sand grains, a revised form was adopted [68], giving the following expressions of the normal stiffness $k_n$ and the tangential stiffness $k_t$:

$$ k_n = k_{n0} \left( \frac{f_n}{G_g l^2} \right)^n, \quad k_t = k_{t0} \left( \frac{f_n}{G_g l^2} \right)^n \quad (95) $$

where $G_g$ is the shear modulus for the grains, $k_{n0}$, $k_{t0}$ and n are material constants.

For the plastic part, the yield function is assumed to be of the Mohr–Coulomb type:

$$ F(f_i, \kappa) = T - f_n \kappa(\Delta^p) = 0 \quad (96) $$

where $\kappa(\Delta^p)$ is an isotropic hardening/softening parameter. The shear force T and the rate of plastic sliding $d\Delta^p$ can be defined as:

$$ T = \sqrt{f_s^2 + f_t^2} \quad \text{and} \quad d\Delta^p = \sqrt{(d\delta_s^p)^2 + (d\delta_t^p)^2} \quad (97) $$

Plastic sliding often occurs along the tangential direction of the contact plane with an upward or downward movement, thus shear dilation/contraction takes place. The dilatancy effect can be described by

$$ \frac{d\delta_n^p}{d\Delta^p} = \frac{T}{f_n} - \tan \phi_0 \quad (98) $$

where $\phi_0$ is a material constant which, in most cases, can be considered equal to the internal friction angle $\phi_\mu$.

The hardening function is defined by a hyperbolic curve in $\kappa-\Delta^p$ plane which involves two material constants $\phi_p$ and $k_{p0}$. Thus,

$$ \kappa = \frac{k_{p0} \tan \phi_p \Delta^p}{|f_n| \tan \phi_p + k_{p0} \Delta^p} \quad (99) $$

On the yield surface, under a loading condition, the shear plastic flow is determined by a normality rule applied to the yield function. However, the plastic flow in the direction normal to the contact plane is governed by the stress-dilatancy equation in Equation (98). Thus the flow rule is non-associated.

The peak friction angle, $\phi_p$, on a contact plane is dependent on the degree of interlocking created by neighbouring particles, which can be related to the state of the packing void ratio e by:

$$ \tan \phi_p = \left( \frac{e_c}{e} \right)^m \tan \phi_\mu \quad (100) $$
<table><caption>Table III. Model parameters for Hostun S28 Sand.</caption><thead><tr><th>e<sub>ref</sub></th><th>p<sub>ref</sub> (MPa)</th><th>λ</th><th>φμ (deg.)</th><th>φ<sub>0</sub> (deg.)</th><th>m</th></tr></thead><tbody><tr><td>1</td><td>0.1</td><td>0.2</td><td>32</td><td>32</td><td>1</td></tr></tbody></table>

<figure><img src="image_24.png" /><figcaption>Figure 23. Comparison of predicted (continuous lines) and measured (points) results for undrained triaxial test on Hostun sand (for experimental points, see [17]).</figcaption></figure>

where the internal friction angle $\phi_{\mu}$ and m are material constants and $e_c$ is the void ratio at the critical state, which is in turn a function of the mean stress. The relationship has traditionally been written as:

$$e_c = \Gamma - \lambda \log(p') \quad \text{or} \quad e_c = e_{\text{ref}} - \lambda \log\left(\frac{p'}{p_{\text{ref}}}\right) \qquad (101)$$

where $\Gamma$ and $\lambda$ are two material constants and $p'$ is the mean stress of the packing, and ($e_{\text{ref}}, p_{\text{ref}}$) is a reference point on the critical state line.

In order to illustrate the model's ability to predict the onset of instability in granular materials, we chose to concentrate on two specific loading conditions: undrained triaxial and constant-q tests on loose Hostun sand. The mean size of the particles for fine Hostun sand is $d=0.33$ mm. The inter-particle elastic constant $k_{n0}$ is assumed to be equal to 61 000 N/mm. The value of $k_{t0}/k_{n0}$ is commonly about 0.4, corresponding to Poisson's ratio for Hostun Sand $\nu=0.2$ and the exponent $n=0.5$ (see Equation (95)). The parameters corresponding to the plastic behaviour were determined from drained triaxial tests (Table III).

Figure 23 presents numerical and experimental results for an undrained triaxial test on loose Hostun sand with an initial confining stress $\sigma_{3c}=100$ kPa. Results indicate that the model is capable of capturing the inception of instability, which corresponds, for this type of test, to the q stress peak, as discussed in the previous section. Similar results have been obtained with different initial effective mean stresses.

Similarly, we also used the parameters in Table III to predict the results of constant-q tests on loose Hostun sand. The predicted and measured results for the confining stress $\sigma_{3c}=300$ kPa are presented in Figure 24. The initial part of the $p'-\varepsilon_v$ curve shows that, as the mean stress $p'$ decreases, the volume increases. This trend continues until a certain point where the volume starts to decrease. For constant-q tests ($dq=0$), the second-order work is reduced to $d_2W=d_p d\varepsilon_v$. Since the mean stress is progressively decreased (i.e. $dp<0$), the second-order work becomes negative if and only if $d\varepsilon_v \le 0$ (i.e. the volume contracts). Thus, the onset of instability corresponds to the peak of the $p'-\varepsilon_v$ curve, which is well reproduced by the model simulation.

Instability conditions obtained by numerical simulations of undrained compression and constant-q tests are plotted in the $p'-q$ plane together with the experimental results obtained on loose
<figure><img src="image_25.png" /><figcaption>Figure 24. Comparison of predicted and measured results for constant *q* test on Hostun sand (for experimental points, see [17]).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 25. Comparison of predicted and measured instability condition for loose Hostun sand determined from undrained triaxial and constant-q tests (from [69]).</figcaption></figure>

Hostun sand (Figure 25, see [69]). The condition of instability in the *p'*-q plane is found to be the same for undrained and constant-q triaxial tests, in agreement with the theoretical developments presented above. This condition defines an instability line for a mobilized friction angle equal to 16°, much lower than the friction angle at critical state equal to 30° [70]. One can see that the model is capable of predicting very accurately the condition of instability associated with these two types of tests. The position of the instability line determined by the model simulations is in very good agreement with the position obtained experimentally.

## 5. CONCLUSIONS

Several teams are working in France, Italy and Canada in close connection around the question of the application of second-order work criterion in geomechanics. The first objective of this paper was to propose a synthesis of the main results obtained recently by these teams. All these results converge towards a novel view of failure in geomaterials. Essentially because of the non-associative character of plastic strains in these materials, all conceivable failure states are no more restricted to
a single plastic limit surface. Indeed a whole bifurcation domain in the stress space has emerged, containing various possible failure modes.

First in this paper, the existence of this bifurcation domain was observed from experimental results based on undrained triaxial compressions and on *q*-constant drained paths. Second, direct numerical simulations using a discrete element method have also shown such a failure domain in close qualitative agreement with the lab experiments. The main characteristics of diffuse failure have been discussed: its directional character with respect to loading, the existence of bursts of kinetic energy, the importance of the control variables, and the absence of visible localization pattern.

In order to analyse these failure states strictly inside the plastic limit condition, the proper criterion appears to be the second-order work criterion. This criterion was introduced and analysed from a theoretical point of view and its link with the bursts of kinetic energy, typical of the transition from the quasi-static pre-failure regime to the dynamic post-failure regime, was established. Eventually the equivalence in geomechanics between material instability, loss of controllability and loss of sustainability was clarified. From these theoretical analyses emerge the proof of the existence of a bifurcation domain and the fact that second-order work criterion is the first failure condition to be met along a loading path, if we exclude flutter instabilities. In addition, the second-order work criterion addresses both the plastic limit condition and strain localization, which appear as special cases of the vanishing of second-order work. Thus, the second-order work criterion necessarily plays a central role in geomechanics with respect to safety against failure of geostructures. Finally, the main features of this criterion are discussed. It is directional with respect to loading, while it is also related to mixed modes of loading on homogeneous samples. Thus, the three necessary and sufficient conditions for an effective failure can now be given: (i) the stress state has to belong to the bifurcation domain, (ii) the current loading direction has to be inside the instability cone, and (iii) the loading variables have to be mixed and proper energy conjugate variables. As such, we have here a generalization of the classical conditions of plastic failure where, as it is well known: (i) the stress state is limited to satisfying the plastic limit condition, (ii) the current loading direction has to be inside the plastic loading cone, and (iii) the loading variables have to be controlled solely by stresses.

Moreover, a generalization of the notion of flow rule has also been obtained. The flow rule is classically a homogeneous relation (of degree 1) linking incremental plastic strain components at failure. It corresponds to the eigenvector related to the first vanishing eigenvalue of the tangent elastoplastic matrix. This notion has been generalized into what has been called a 'failure rule', which is a homogeneous relation of the mixed type relating both incremental stresses and strains at failure. The latter relation corresponds to the eigenvector associated to the first vanishing eigenvalue of the tangent generalized constitutive matrix which links mixed loading variables to the mixed conjugate response variables.

In the third (and last) part of this synthesis, these theoretical results are checked by means of various constitutive models. First, three phenomenological elastoplastic relations of increasing complexity are advocated, followed by two micromechanical models. All these constitutive models confirm the observed experimental results as well as the direct discrete numerical ones including theoretical findings, thus constituting altogether a firm basis for a proper use of second-order work criterion in geomechanics.

This synthesis is long enough and any application to geotechnical problems has been excluded from this paper. However, the interested reader will find some illustrations of second-order work applications, for example to landslide modelling, in Darve and Laouafa [48], Laouafa and Darve [49], Lignon et al. [71]. It is noteworthy that Pastor et al. [72] have also observed diffuse failures in their computational modelling.

## REFERENCES

1. Rudnicki JW, Rice JR. Conditions of localization of deformation in pressure-sensitive dilatant materials. *Journal of the Mechanics and Physics of Solids* 1975; **29**:153–172.

2. Vardoulakis I, Sulem J. *Bifurcation Analysis in Geomechanics*. Chapman & Hall Publications: London, 1995.
3. Khoa HDV, Geogopoulos IO, Darve F, Laouafa F. Diffuse failure in geomaterials, experiments and modelling. *Computers and Geotechnics* 2006; **33**:1–14.

4. Lyapunov AM. Problème général de la stabilité des mouvements. *Annales de la Faculté des Sciences de Toulouse* 1907; **9**:203-274.

5. Nova R. Liquefaction, stability, bifurcations of soil via strain-hardening plasticity. In *Numerical Methods for the Localisation and Bifurcation of granular bodies; Proceedings of the International Workshops*, Gdansk, Dembicki E, Gudehus G, Sikora Z (eds). (1989); 117-132.

6. Nicot F, Darve F, Khoa H. Bifurcation, second order-work in granular materials. *International Journal for Numerical and Analytical Methods in Geomechanics* 2007; **31**:1007-1032.

7. Darve F, Chau B. Constitutive instabilities in incrementally non-linear modelling. In *Constitutive Laws for Engineering Materials*, Desai CS (ed.). 1987; 301-310.

8. Hill R. A general theory of uniqueness and stability in elastic-plastic solids. *Journal of Mechanics and Physics of Solids* 1958; **5**:236-249.

9. Mandel J. Conditions de stabilité et postulat de Drucker. In *Rheology and Soil Mechanics*, Kravtchenko J, Sierieys PM (eds). Springer: Berlin, 1966; 58-68.

10. Mroz Z. On forms of constitutive laws for elastic-plastic solids. *Archives of Mechanics* 1966; **18**:1-34.

11. Bazant ZP, Cedolin L. *Stability of Structures: Elastic, Inelastic, Fracture and Damage Theories*. Oxford University Press: New York, 1991.

12. Bigoni D, Hueckel T. Uniqueness localization—I. Associative and nonassociative elastoplasticity. *International Journal of Solids and Structures* 1991; **28**(2):197-213.

13. Petryk H. *Material Instabilities in Elastic and Plasticsolids*. Springer: Berlin, 2000.

14. di Prisco C, Imposimato S. Time dependent mechanical behaviour of loose sands. *Mechanics of Cohesive-Frictional Materials* 1996; **1**:45-73.

15. di Prisco C, Imposimato S. Experimental analysis and theoretical interpretation of triaxial load controlled loose sand specimen collapses. *Mechanics of Cohesive-Frictional Materials* 1997; **2**:93-120.

16. Imposimato S, Nova R. An investigation on the uniqueness of the incremental response of elastoplastic models for virgin sand. *Mechanics of Cohesive-Frictional Materials* 1998; **3**:65-87.

17. Darve F, Sibille L, Daouadji A, Nicot F. Bifurcations in granular media: macro-and micro-mechanics approaches. *Comptes Rendus Mecanique* 2007; **335**:496-515.

18. Sasitharan S, Robertson PK, Sego DC, Morgenstern NR. Collapse behaviour of sand. *Canadian Geotechnical Journal* 1993; **30**:569-577.

19. Skopek P, Morgenstern NR, Robertson PK, Sego DC. Collapse of dry sand. *Canadian Geotechnical Journal* 1994; **31**:1008-1014.

20. Gajo A. The influence of system compliance on collapse of triaxial sand samples. *Canadian Geotechnical Journal* 2004; **41**:257-273.

21. Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. *Geotechnique* 1979; **29**(1):47-65.

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

23. Darve F, Servant G, Laouafa F, Khoa H. Failure in geomaterials: continuous and discrete analyses. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**(27-29):3057-3085.

24. Nova R. Controllability of the incremental response of soil specimens subjected to arbitrary loading programmes. *Journal of the Mechanical Behavior of Materials* 1994; **5**(2):193-201.

25. Sibille L, Nicot F, Donzé F, Darve F. Analysis of failure occurrence from direct simulations. *European Journal of Environmental and Civil Engineering* 2009; **13**:187-201.

26. Rice JR. The localization of plastic deformation. In *Theoretical and Applied Mechanics*, Koiter WT (ed.). North-Holland: Amsterdam, 1976.

27. Nicot F, Darve F. A micro-mechanical investigation of bifurcation in granular materials. *International Journal of Solids and Structures* 2007; **44**:6630-6652.

28. Ziegler H. *Principles of Structural Stability*. Blaisdell Publishing Company: Waltham, 1968.

29. Schur I. Uber Potenzreihen, die im Innern des Einheitskreises beschraenkt sind. *Reine und Angewandte Mathematik* 1917; **147**:205-232.

30. Ostrowski A, Taussky O. On the variation of the determinant of a positive definite matrix. *Nederlandse Akademie van Wetenschappen, Proceedings* 1951; **A54**:333-351.

31. di Prisco C, Nova R, Lanier J. A mixed isotropic-kinematic hardening constitutive law for sand. *Modern Approaches to Plasticity*, Kolymbas D (ed.). Elsevier: Amsterdam, 1993; 83-124.

32. Nova R, Imposimato S. Non-uniqueness of the incremental response of soil specimens under true-triaxial stress paths. In *Proceedings of the VI NUMOG*, Pande, Pietruszczak (eds). Balkema: Montreal, 1997; 193-197.

33. Kramer SL, Seed HB. Initiation of soil liquefaction under static loading conditions. *Journal of Geotechnical Engineering* (ASCE) 1988; **114**(4):412-430.

34. Lade PV. Static instability, liquefaction of loose fine sandy slopes. *Journal of Geotechnical Engineering* (ASCE) 1992; **118**:51-71.

35. Imposimato S, Nova R. Instability of loose sand specimens in undrained tests. In *Proceedings of the 4th International Workshop on Localisation and Bifurcation Theory for Soils and Rocks*, Adachi T, Oka F, Yashima A (eds). Balkema: Rotterdam, 1998; 313-322.
36. Gudehus G. A comparison of some constitutive laws for soils under radially symmetric loading and unloading. In *Third International Conference on Numerical Methods in Geomechanics*, Aachen, Wittke (eds), vol. 4. Balkema Publisher: Rotterdam, 1979; 1309–1324.

37. Nicot F, Sibille L, Darve F. Bifurcation in granular materials: an attempt at a unified framework. *International Journal of Solids and Structures* 2009; **46**:3938–3947.

38. Nova R. Sinfonietta classica: an exercise on classical soil modelling. In *Proceedings, Symposium on Constitutive Equations for Granular Non-cohesive Soils*, Saada CA, Bianchini G (eds). Balkema: Rotterdam, 1988; 501–519.

39. Schofield AN, Wroth CP. *Critical State Soil Mechanics*. McGraw-Hill: New York, 1968.

40. Butterfield R. A natural compression law for soils. *Géotechnique* 1979; **29**:469–480.

41. Nova R. A note on sand liquefaction and soil stability. *Proceedings of the Third International Conference on Constitutive Laws for Engineering Materials: Theory and Applications*, Tucson, AZ, U.S.A., 7–12 January 1991; 153–156.

42. di Prisco C, Matiotti R, Nova R. Theoretical investigation of the undrained stability of shallow submerged slopes. *Géotechnique* 1995; **45**(3):479–496.

43. Wan RG, Guo PJ. Effect of microstructure on undrained behaviour of sands. *Canadian Geotechnical Journal* 2001; **38**:16–28.

44. Wan RG, Guo PJ. Drained cyclic behavior of sand with fabric dependence. *Journal of Engineering Mechanics* 2001; **127**(11):1106–1116.

45. Oda M, Konishi J. Microscopic deformation mechanism of granular material in simple shear. *Soils and Foundations* 1974; **14**:25–38.

46. Darve F, Labanieh S. Incremental constitutive law for sands and clays. Simulations of monotonic and cyclic tests. *International Journal for Numerical and Analytical Methods in Geomechanics* 1982; **6**:243–275.

47. Darve F, Flavigny E, Meghachou M. Yield surfaces and principle of superposition revisited by incrementally nonlinear constitutive relations. *International Journal of Plasticity* 1995; **11**(8):927–948.

48. Darve F, Laouafa F. Instabilities in granular materials and application to landslides. *Mechanics of Cohesive-Frictional Materials* 2005; **8**:627–652.

49. Laouafa F, Darve F. Modelling of slope failure by a material instability mechanism. *Computers and Geotechnics* 2002; **29**(4):301–325.

50. Prunier F, Laouafa F, Lignon S, Darve F. Bifurcation modeling ingecomaterials: from the second order work criterion to spectral analyses. *International Journal for Numerical and Analytical Methods in Geomechanics* 2009; **33**:1169–1202.

51. Prunier F, Nicot F, Darve F, Laouafa F, Lignon S. 3D multi scale bifurcation analysis of granular media. *Journal of Engineering Mechanics* 2009; **135**(6):493–509.

52. Oda M, Konishi J, Nemat-Nasser S. Experimental evaluation of strength of granular materials, effects of particle rolling. *Mechanics of Materials* 1982; **1**:269–283.

53. Radjai F, Roux S, Moreau JJ. Contact forces in a granular packing. *Chaos* 1999; **9**(3):544–550.

54. Nicot F, Darve F. A multiscale approach to granular materials. *Mechanics of Materials* 2005; **37**(9):980–1006.

55. Chang CS, Hicher P-Y. An elasto-plastic model for granular materials with microstructural consideration. *International Journal of Solids and Structures* 2005; **42**(14):4258–4277.

56. Nicot F. Constitutive modelling of a snowcover with a change in scale. *European Journal of Mechanics A/Solids* 2003; **22**(3):325–340.

57. Nicot F. From constitutive modelling of a snowcover to the design of flexible structures. Part I, mechanical modelling. *International Journal of Solids and Structures* 2004; **41**:3317–3337.

58. Nicot F. Constitutive modelling of snow as a cohesive granular material. *Granular Matter* 2004; **6**:47–60.

59. Oda M. The mechanism of fabric changes during compressional deformation of sand. *Soils and Foundations* 1972; **12**:1–23.

60. Calvetti F, Combe G, Lanier J. Experimental micro-mechanical analysis of a 2D granular material, relation between structure evolution and loading path. *Mechanics of Cohesive-Frictional Materials* 1997; **2**:121–163.

61. Love AEH. *A Treatise of Mathematical Theory of Elasticity*. Cambridge University Press: Cambridge, 1927.

62. Weber J. Recherches concernant les contraintes intergranulaires dans les milieux pulvérulents. *Bull. Liaison P. et Ch.* 1966; **20**:1–20.

63. Christoffersen J, Mehrabadi MM, Nemat-Nasser S. A micro-mechanical description of granular material behavior. *Journal of Applied Mechanics* 1981; **48**:339–344.

64. Mehrabadi MM, Oda M, Nemat-Nasser S. On statistical description of stress and fabric in granular materials. *International Journal for Numerical and Analytical Methods in Geomechanics* 1982; **6**:95–108.

65. Cambou B, Chaze M, Dedecker F. Change of scale in granular materials. *European Journal of Mechanics A/Solids* 2000; **19**:999–1014.

66. Liao CL, Chang TP, Young D, Chang CS. Stress–strain relationship for granular materials bases on hypothesis of best fit. *International Journal of Solids and Structures* 1997; **34**(31–32):4087–4100.

67. Mindlin RD. Microstructure in linear elasticity. *Archive for Rational Mechanics and Analysis* 1969; **16**:51–78.

68. Chang CS, Sundaram SS, Misra A. Initial moduli of particulate mass with frictional contacts. *International Journal for Numerical and Analytical Methods in Geomechanics* 1989; **13**(6):626–641.

69. Chang CS, Hicher PY, Daouadji A. Investigating instability in granular materials by means of a micro-structural model. *European Journal of Environmental and Civil Engineering* 2009; **13**(3):167–186.
70. Daouadji A, AlGali H, Darve F, Zeghloul A. Instability in granular materials: an experimental evidence of diffuse mode of failure for loose sands. *Journal of Engineering Mechanics (ASCE)* 2010; **136**(5):10. DOI: 10.1061/_ASCE_EM.1943-7889.101.

71. Lignon S, Laouafa F, Prunier F, Khoa HDV, Darve F. Hydromechanical modelling of landslides with a material instability criterion. *Geotechnique* 2009; **59**(6):513–524.

72. Pastor M, Haddad B, Sorbino G, Cuomo S, Drempetic V. A depth integrated coupled SPH model for flow like landslides and related phenomena. *International Journal for Numerical and Analytical Methods in Geomechanics* 2009; **33**(2):143–172.