# Investigating instability in granular materials by means of a microstructural model

Ching S. Chang, Pierre-Yves Hicher, Ali Daouadji

► To cite this version:

Ching S. Chang, Pierre-Yves Hicher, Ali Daouadji. Investigating instability in granular materials by means of a microstructural model. European Journal of Environmental and Civil Engineering, 2009, 13 (2), pp.167-186. 10.1080/19648189.2009.9693098. hal-01006997

HAL Id: hal-01006997

https://hal.science/hal-01006997

Submitted on 8 Apr 2017

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

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

Public Domain
# Investigating instability in granular materials by means of a microstructural model

Ching S. Chang* — Pierre-Yves Hicher** — Ali Daouadji***

* Department of Civil and Environmental Engineering
University of Massachusetts, Amherst, MA 01003, USA

** Institut de Recherche en Génie Civil et Mécanique - UMR CNRS 6183
Ecole Centrale de Nantes, 1, rue de la Noë BP 92101, F-44 321 Nantes Cedex 3, France
pierre-yves.hicher@ec-nantes.fr

*** Laboratoire de Physique et Mécanique des Matériaux – UMR CNRS 7554
Université Paul Verlaine – Metz, Ile du Saulcy, F-57045 Metz Cedex, France

**ABSTRACT.** Under certain loading conditions, instability in loose sand can develop at a shear stress level much lower than the critical state failure line. To analyse these types of problems, we adopt the micromechanical approach developed earlier for the modelling of granular material behaviour. The stress-strain relationship for a granular assembly is determined by integrating the behaviour of the inter-particle contacts in all orientations. The constitutive model is applied to simulate undrained triaxial and constant-q tests on loose sand. Numerical simulations are compared to experimental results in order to evaluate the model's ability to predict the different modes of instability of the granular assembly. With the model we have also analysed instability at inter-particle level to observe how it is connected to instability at the assembly level.

**RÉSUMÉ.** Dans certaines conditions de chargement, des ruptures par instabilité peuvent prendre naissance dans les sables lâches pour des niveaux de contraintes très nettement inférieurs à ceux correspondant à la plasticité parfaite. Pour analyser ce type de problèmes, nous adoptons une approche d'homogénéisation développée précédemment pour la modélisation du comportement des matériaux granulaires. La relation contrainte-déformation macroscopique est obtenue en intégrant le comportement défini au niveau du contact inter-granulaire. Le modèle est utilisé pour simuler des essais non drainés et des essais à déviateur constant sur un sable lâche. Les résultats de la simulation à l'échelle macroscopique sont comparés à des résultats expérimentaux de façon à apprécier la capacité du modèle à prédire la naissance de l'instabilité sur ces chemins particuliers. Une analyse à l'échelle des différents plans de contact inter-granulaire a également été entreprise afin de relier les conditions d'instabilité aux différentes échelles.

**KEYWORDS:** granular material, instability, micromechanics, stress-strain relationship, sand.

**MOTS-CLÉS** : matériau granulaire, instabilité, micromécanique, relation contrainte-déformation, sable.
# 1. Introduction

Instability of granular material is an important topic in geotechnical engineering because it may lead to catastrophic events such as the collapse of earth structures. This phenomenon can be categorised into two major aspects: material instability (also known as intrinsic/constitutive instability) and geometrical instability (see, for example, (Goddard, 2003)). The development of localisation in a specimen of finite size can be studied as a boundary value problem for which initial inhomogeneities, as well as boundary conditions, play an important role. However, it has been shown by Rudnicki and Rice (1975) that the occurrence of localisation can be predicted at a constitutive level. More recently, several authors have demonstrated that other modes of instability can occur within granular materials. Darve and co-authors (2004, 1998, 2002) have studied the mode of instability which they called diffuse failure, using a macroscopic phenomenological approach based on an incrementally non-linear model. Nova (1994) has developed similar concepts and demonstrated that any non-associated elastoplastic model could exhibit domains of instability inside the plastic limit. Experimental evidences support these theoretical approaches. For example, in loose sand under undrained conditions, an unstable condition can be obtained at a low shear stress level and, subsequently, the strength is reduced to almost zero, corresponding to a material state known as static liquefaction.

In this paper, we present a study of material instability, using the micromechanical approach developed by Chang and Hicher (2004). Under this approach, the stress-strain relationship for a granular assembly can be determined by integrating the behaviour of the inter-particle contacts in all orientations, based on a static hypothesis which relates the average stress of the granular assembly to the mean field of particle contact forces. Model simulations are compared to experimental results obtained on loose Hostun sand along two different loading paths: undrained triaxial tests and constant – q tests.

# 2. Experimental evidences of instability in sand

## 2.1. Undrained triaxial tests

Undrained triaxial tests on Hostun sand with various initial relative densities, $D_r$, are presented in Figure 1 (Hicher, 1998). These results show two distinctive trends corresponding to a contractive or a dilative behaviour. For loose sand (small values of $D_r$), a maximum strength is reached in the $q-\varepsilon_1$ plane, and the peak in the stress-strain curve is followed by a rapid decrease of the deviatoric stress down to a minimum strength. This peak corresponds to the development of material instability as will be analysed later on. The minimum strength can be almost zero for a relative density close to zero. This represents the phenomenon called static liquefaction. For medium dense sand (medium values of $D_r$), the tendency of softening still occurs at the beginning of the loading, but it is followed by an increase of the deviatoric stress
up to the ultimate strength, corresponding to the critical state at large deformations. For dense sand (high values of $D_r$), the material is strongly dilatant and no strain softening can be observed. Instead, the deviatoric stress continuously increases up to the ultimate strength achieved at large deformations. These different evolutions can be related to the stress paths followed in a p'-q plane; a continuous decrease of the mean effective stress is observed for strongly contractant materials; a decrease followed by an increase of the mean effective stress up to the critical state is observed for dilatant materials.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Experimental results for undrained triaxial tests on Hostun sand with various densities: a) stress- strain curves, and b) stress paths (Hicher, 1998)</figcaption></figure>

If we examine the position of the peak in the stress plane p'-q, we can see that instability occurs at a stress state below the critical state failure line. Depending on the density of the sand, the position of the peak stress can be more or less distant from the failure line.

## 2.2. Constant-q tests

This type of test consists in shearing the specimen to a prescribed stress ratio along a drained compression triaxial path, and then in decreasing the mean effective stress while keeping constant the deviatoric stress. This stress path can simulate the loading condition of a soil element within a slope when a progressive increase in pore pressure occurs.

Several investigations have demonstrated that instability can occur in loose sand during a constant-q stress path (Sasitharan *et al.*, 1993; Nova and Imposimato, 1997; Gajo *et al.*, 2000; Lade, 2002; Chu and Leong, 2003). Darve *et al.* (2007) have presented similar test results on Hostun sand and observed a sudden collapse of the
sand specimen for stress states located well below the critical state failure line. Typical results of constant-q tests on loose Hostun sand are presented in Figure 2. After an isotropic consolidation stage to a desired initial effective mean pressure p'o, a drained triaxial compression test was applied to the sample up to a prescribed value of the deviatoric stress q. Then, while keeping q constant, a decrease of the mean effective stress p' was applied by increasing the pore water pressure and maintaining constant the total stresses. At a given point of the test, the axial strain rate started to increase very rapidly and the deviatoric stress could no longer be kept constant. The test is no longer controllable (as defined by Nova (1994)) in the sense that the imposed loading program cannot be maintained. As will be shown later, this point corresponds to a loss of stability, since any small change of one control variable, as defined below, will lead to a catastrophic failure.

<figure><img src="image_2.png" /><figcaption>Figure 2. Experimental results of constant-q tests on loose Hostun sand</figcaption></figure>

## 2.3. Analysis of undrained compression and constant-q tests in terms of instability condition

The importance of the stress path can be linked to Hill's sufficient condition of stability (Hill, 1956), which states that a material, progressing from one stress state to another, is stable if the second-order work is strictly positive, i.e.:

$$d^2W = d\sigma_{ij} d\varepsilon_{ij} > 0 \quad [1]$$

Thus, according to Hill's condition, whether a material is stable or not depends not only on the current stress state but also on the direction of the stress increment.

Because $q - \varepsilon_1$ and $\varepsilon_v - \sigma'_3$ are conjugate variables with respect to energy, Equation [1] can be re-arranged to give (Darve et al., 2003):
$$d^2W = dq d\epsilon_1 + d\sigma'_3 d\epsilon_v \quad [2]$$

Thus, as the control variables are dq and dε<sub>v</sub>, the constitutive relation linking the stress increments to the strain increments can be rearrange to give a generalized mixed relation between generalised incremental stresses and generalised incremental strains can be given as:

$$\left\{ \begin{matrix} dq \\ d\epsilon_v \end{matrix} \right\} = \left[ \begin{matrix} E_1 & 2 \frac{E_1}{E_3} v_{13} - 1 \\ 1 - 2v_{31} & \frac{2(1 - v_{33} - 2v_{13}v_{31})}{E_3} \end{matrix} \right] \left\{ \begin{matrix} d\epsilon_1 \\ d\sigma_3 \end{matrix} \right\} = [P] \left\{ \begin{matrix} d\epsilon_1 \\ d\sigma_3 \end{matrix} \right\} \quad [3]$$

where E<sub>i</sub> are the pseudo Young moduli and v<sub>i</sub> are the pseudo Poisson's ratios. As no volumetric variation is allowed during undrained tests, Equation [3] is modified as follows:

$$\left\{ \begin{matrix} dq \\ 0 \end{matrix} \right\} = \left[ \begin{matrix} E_1 & 2 \frac{E_1}{E_3} v_{13} - 1 \\ 1 - 2v_{31} & \frac{2(1 - v_{33} - 2v_{13}v_{31})}{E_3} \end{matrix} \right] \left\{ \begin{matrix} d\epsilon_1 \\ d\sigma_3 \end{matrix} \right\} \quad [4]$$

Equation [4] indicates that instability will occur only if dq = 0 (peak shear stress); the previous equation is therefore:

$$\left\{ \begin{matrix} 0 \\ 0 \end{matrix} \right\} = \left[ \begin{matrix} E_1 & 2 \frac{E_1}{E_3} v_{13} - 1 \\ 1 - 2v_{31} & \frac{2(1 - v_{33} - 2v_{13}v_{31})}{E_3} \end{matrix} \right] \left\{ \begin{matrix} d\epsilon_1 \\ d\sigma_3 \end{matrix} \right\} \quad [5]$$

Thus, instability can take place only if the determinant of the matrix P of the constitutive relation given by Equation [5] becomes equal to zero.

As no variation of the shear stress is imposed during constant-q tests (dq = 0), Equation [3] can be rewritten as:

$$\left\{ \begin{matrix} 0 \\ d\epsilon_v \end{matrix} \right\} = \left[ \begin{matrix} E_1 & 2 \frac{E_1}{E_3} v_{13} - 1 \\ 1 - 2v_{31} & \frac{2(1 - v_{33} - 2v_{13}v_{31})}{E_3} \end{matrix} \right] \left\{ \begin{matrix} d\epsilon_1 \\ d\sigma_3 \end{matrix} \right\} \quad [6]$$

Taking into account Equation [3], instability will occur in this case only if dε<sub>v</sub> = 0, because |dσ'<sub>3</sub>| > 0 is imposed by the loading program. The condition is the same as the one obtained for undrained tests (Equation [4]) and both correspond to:
$$ \det(P) = 1 - 2\nu_{31} + 2(1 - \nu_{33} - \nu_{13})\frac{E_1}{E_3} \quad [7] $$

Equation [7] is verified for non trivial solutions if $det(P) = 0$. As described by Darve et al. (2004), the condition of instability for constant-q tests coincides with the condition for undrained tests, so the stress states in the p'-q plane corresponding to the vanishing of the volumetric strain increment during constant-q tests are the same as the stress states corresponding to the peaks of the stress-strain curves for the undrained tests. Figure 3 presents different results obtained on very loose Hostun sand (Dr = 0). One can see that 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 mobilised friction angle equal to 16°, much lower than the friction angle at critical state equal to 30°.

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Instability line for loose Hostun sand determined from undrained triaxial and constant-q tests</figcaption></figure>

### 3. Stress-strain model based on micromechanical approach

In this Section, the microstructural model developed by Chang and Hicher (2004) is briefly described. In this model, a granular material is viewed as a collection of particles. The deformation of a representative volume of the material is generated by the mobilisation of particle contacts in all orientations. On each contact plane, an auxiliary local coordinate can be established by means of three orthogonal unit vectors $n$, $s$, $t$. The vector $n$ is outward normal to the contact plane. Vectors $s$ and $t$ are on the contact plane.
### 3.1. Inter-particle behaviour

Elastic stiffness: the contact stiffness of a contact plane $k_n^\alpha$ includes normal stiffness, $k_n^\alpha$, and shear stiffness, $k_t^\alpha$. The elastic stiffness tensor is defined by:

$$f_i^\alpha = k_{ij}^{\alpha e} \delta_j^{\alpha e} \qquad [8]$$

which can be related to the contact normal and shear stiffness:

$$k_{ij}^{\alpha e} = k_n^{\alpha} n_i^{\alpha} n_j^{\alpha} + k_t^{\alpha} (s_i^{\alpha} s_j^{\alpha} + t_i^{\alpha} t_j^{\alpha}) \qquad [9]$$

The value of the stiffness for two elastic spheres can be estimated from Hertz-Mindlin's formulation (1953). For sand grains, a revised form was adopted (Chang et al., 1989), given by:

$$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 \qquad [10]$$

where $G_g$ is the elastic modulus for the grains, $f_n$ is the contact force in normal direction, $l$ is the branch length between two particles, $k_{n0}, k_{t0}$ and $n$ are material constants.

*Plastic yield function:* the yield function is assumed to be of the Mohr-Coulomb type, defined in a contact-force space (e.g. $f_n, f_s, f_t$):

$$F(f_i, \kappa) = T - f_n \kappa (\Delta^p) = 0 \qquad [11]$$

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

$$T = \sqrt{f_s^2 + f_t^2} \quad \text{and} \quad \Delta^p = \sqrt{(\delta_s^p)^2 + (\delta_t^p)^2} \qquad [12]$$

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

$$\kappa = \frac{k_{p0} \tan \phi_p \Delta^p}{|f_n| \tan \phi_p + k_{p0} \Delta^p} \qquad [13]$$
*Plastic flow rule:* the plastic sliding often occurs along the tangential direction of the contact plane with an upward or downward movement, thus, shear-induced 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 [14] $$

where the material constant $\phi_0$ can be considered in most cases equal to the inter-particle friction angle $\phi_\mu$. 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 [14]. Thus the flow rule is non-associated.

*Elasto-plastic relationship:* with the elements discussed above, the incremental force-displacement relationship of the inter-particle contact can be obtained. Including both elastic and plastic behaviour, this relationship is given by:

$$ \dot{f}_i^\alpha = k_{ij}^{\alpha p} \dot{\delta}_j^\alpha \quad [15] $$

Detailed expression of the elasto-plastic stiffness tensor can be derived in a straightforward manner from yield function and flow rule, and, therefore, is not given here.

## 3.2. Influence of void ratio on $\phi_p$

The resistance against sliding in a contact plane is dependent on the degree of interlocking by neighbouring particles. The resistance can be related to the packing void ratio $e$ by:

$$ \tan \phi_p = \left( \frac{e_c}{e} \right)^m \tan \phi_\mu \quad [16] $$

where $m$ is a material constant (Biarez and Hicher, 1994) and $e_c$ corresponds to the critical void ratio for a given state of stress. For dense packing, $e_c/e$ is greater than 1 and therefore the apparent inter-particle friction angle $\phi_p$ is greater than the internal friction angle $\phi_\mu$. When the packing structure dilates, the degree of interlocking and the apparent friction angle is reduced, which results in a strain-softening phenomenon. For loose packing, the apparent friction angle $\phi_p$ is smaller than the internal friction angle $\phi_\mu$ and increases during the material contraction.

The critical void ratio $e_c$ is a function of the mean stress applied to the overall assembly and can be written as follows (Schofield, Wroth, 1968):
$$e_c = \Gamma - \lambda \log(p') \quad \text{or} \quad e_c = e_{ref} - \lambda \log\left(\frac{p'}{p_{ref}}\right) \qquad [17]$$

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

### 3.3. Micro-macro relationship

The stress-strain relationship for an assembly can be determined by integrating the behaviour of inter-particle contacts in all orientations. In the integration process, a micro-macro relationship is required. Following the Love-Weber approach, the stress increment can be obtained by the contact forces and branch vectors for all contacts (Christofferson et al., 1981; Rothenburg and Salvadurai, 1981), as follows:

$$\dot{\sigma}_{ij} = \frac{1}{V} \sum_{\alpha=1}^{N} \dot{f}_{j}^{\alpha} l_{i}^{\alpha} \qquad [18]$$

The mean force on the contact plane of each orientation is:

$$\dot{f}_{j}^{\alpha} = \dot{\sigma}_{ij} A_{ik}^{-1} l_{k}^{\alpha} V \qquad [19]$$

where the branch vector $l_k^\alpha$ is defined as the vector joining the centres of two particles, and the fabric tensor is defined as $A_{ik} = \sum_{\alpha=1}^{N} l_i^\alpha l_k^\alpha$.

Using the principle of energy balance, which states that the work done in a representative volume element is equal to the work done on all inter-cluster planes within the element,

$$\sigma_{ij} \dot{u}_{j,i} = \frac{1}{V} \sum_{\alpha=1}^{N} f_{j}^{\alpha} \dot{\delta}_{j}^{\alpha} \qquad [20]$$

From Equations [19] and [20], we obtain the relation between the global strain and the inter-particle displacement:

$$\dot{u}_{j,i} = A_{ik}^{-1} \sum_{\alpha=1}^{N} \dot{\delta}_{j}^{\alpha} l_{k}^{\alpha} \qquad [21]$$

### 3.4. Stress-strain relationship

From Equations [15], [18], and [21], the following relationship between stress increment and strain increment can be obtained:
$$ \dot{u}_{i,j} = C_{ijmp} \dot{\sigma}_{mp} ; \quad \text{where} \quad C_{ijmp} = A_{ik}^{-1} A_{mn}^{-1} V \sum_{\alpha=1}^{N} (k_{jp}^{ep})^{-1} l_k^{\alpha} l_n^{\alpha} \qquad [22] $$

When the contact number $N$ is sufficiently large in an isotropic packing, the summation of the compliance tensor in Equation [22] and the summation of the fabric tensor in Equation [19] can be written in integral form, given by:

$$ C_{ijmp} = A_{ik}^{-1} A_{mn}^{-1} \frac{NV}{2\pi} \int_{0}^{\pi/2} \int_{0}^{2\pi} k_{jp}^{ep} (\gamma, \beta)^{-1} l_k(\gamma, \beta) l_n(\gamma, \beta) \sin \gamma d\gamma d\beta \qquad [23] $$

and:

$$ A_{ik} = \frac{N}{2\pi} \int_{0}^{\pi/2} \int_{0}^{2\pi} l_{i}(\gamma, \beta) l_{k}(\gamma, \beta) \sin \gamma d\gamma d\beta \qquad [24] $$

The integration of Equations [23] and [24] in a spherical coordinate can be carried out numerically by using Gauss integration points over the surface of the sphere.

## 4. Model performances

The model needs a set of input parameters, such as mean particle size, particle stiffness, inter-particle friction, initial porosity, and parameters defining the critical state of the sand. The mean size of the particle for fine Hostun sand is $d = 0.4$ mm. The inter-particle elastic constant $k_{n0}$ is assumed to be equal to 61000 N/mm.

The total number of contacts per unit volume changes during the deformation. Using the experimental data by Oda (1977) for three mixtures of spheres, the total number of contact per unit volume can be approximately related to the void ratio by the following expression:

$$ \frac{N}{V} = \left(\frac{N}{V}\right)_0 \frac{(1+e_0)e_0}{(1+e)e} \qquad [25] $$

where $e_0$ is the initial void ratio of the granular assembly. This equation is used to account for the evolution of the contact number per unit volume. The initial contact number per unit volume can be obtained by matching the predicted and experimentally measured elastic modulus for specimens with different void ratios (Hicher, Chang, 2007).

The value of $k_{t0}/k_{n0}$ is commonly about 0.4, corresponding to a Poisson's ratio for Hostun Sand $\nu = 0.2$ and the exponent $n = 0.5$. From drained triaxial test results, we were able to derive the values of the two parameters corresponding to the position of the critical state in the e-p' plane: $\lambda = 0.2$ and $p_{ref} = 0.01$ MPa for $e_{ref} = e_{max} = 1$. In Equation [16], the value of $m = 1$ was determined from the test
results. The values of $k_{p0}$ are assumed to be same as the elastic stiffness $k_{p0} = k_n$. The set of parameters for fine Hostun sand is presented in Table 1. The model performances will be demonstrated in the following sections by comparing the predicted and measured stress-strain behaviour.

<table><caption>**Table 1.** Model parameters for fine Hostun sand</caption><thead><tr><th>e<sub>ref</sub></th><th>p<sub>ref</sub> (MPa)</th><th>λ</th><th>φ<sub>μ</sub>(°)</th><th>φ<sub>0</sub>(°)</th><th>m</th></tr></thead><tbody><tr><td>1</td><td>0.01</td><td>0.2</td><td>30</td><td>30</td><td>1</td></tr></tbody></table>

## 4.1. Undrained triaxial tests

<figure><img src="image_4.png" /><figcaption>**Figure 4.** Model predictions for undrained triaxial tests on Hostun sand with various relative densities</figcaption></figure>

Figure 4 presents numerical results for undrained triaxial tests on Hostun sand with various initial void ratios, corresponding to relative densities between 0.05 and 1. Both predicted stress-strain curves and stress paths are in agreement with the experimental curves in Figure 1. Results indicate that the model is capable of capturing the general trend observed for contractive and dilative sands. In order to examine the inception of instability, the predicted shear stress and second-order work are plotted against the shear strain in Figure 5 for a test on loose Hostun sand with an initial confining stress $\sigma_{3c} = 300$ kPa. For undrained conditions ($d\epsilon_v = 0$), the second-order work is reduced to $d^2W = dq d\gamma$ where $\gamma$ is the deviatoric strain. As discussed in Section 2.3, since the deviatoric strain increases continuously, $d\gamma$ is always positive and the second-order work can become non-positive if and only if $dq \le 0$ (i.e., decrease in $q$). Figure 5 shows that instability begins at the shear stress peak. The second-order work is positive before the peak
stress. After the peak, the second-order work remains negative and approaches zero at critical state. The numerical simulations can be carried out after the peak because the control variable is the vertical strain $\varepsilon_l$ and not the deviatoric stress $q$.

It should be noted that the undrained condition imposed in the tests is taken into account in the modelling by the condition of no volume change, i.e., isochoric condition. The instability is therefore not triggered by the pore pressure built up, but rather tests on a dry loose sand specimen have also led to the same instability mode (Lanier, Block, 1989).

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Predicted second-order work and stress-strain curves for undrained triaxial test on loose Hostun sand</figcaption></figure>

**4.2. Constant-q tests**

Similarly, we also used the parameters in Table 1 to predict the results of constant-q tests on loose Hostun sand. The predicted and measured results for the confining stress $\sigma_{3c} = 300 \text{ kPa}$ are presented in Figure 6. 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$), according to Equation [1], the second-order work is reduced to $d^2W = dp 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. As for the undrained tests discussed previously, the numerical simulations can be carried out
after the instability condition because the control variable is the mean effective stress $p'$ and not the volume change $\epsilon_v$.

<figure><img src="image_6.png" /><figcaption>**Figure 6.** Comparison of predicted and measured results for constant-q test on Hostun sand</figcaption></figure>

### 4.3. Summary of the condition of instability for loose Hostun sand

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Comparison of predicted and measured instability condition for loose Hostun sand determined from undrained triaxial and constant-q tests</figcaption></figure>
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 presented in Section 2.3 (Figure 7). 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.

### 4.4. Local plane behaviour

In order to gain insight of the material behaviour at the inter-particle level, we selected several contact planes oriented in six different directions indicated by the branch vectors, as shown in Figure 8. The six directions are defined by the angle $\gamma$ between the vector $n$ perpendicular to the contact plane and the vertical axis. A larger value of $\gamma$ indicates a more inclined contact plane.

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Orientations of branch vectors for the inter-particle contact planes selected in the analysis</figcaption></figure>

In a loaded assembly, the forces are different on each contact plane. In order to obtain a more direct comparison between the local and the overall stress-strain behaviour, we retained the following variables, local strain and local stress, instead of inter-particle force and inter-particle displacement. For this purpose, we define a local normal stress $\sigma^\alpha = f_n^\alpha Nl / 3V$ and a local shear stress $\tau^\alpha = f_s^\alpha Nl / 3V$, where $l$ is the branch length and N/V is the total number of contacts per unit volume. The corresponding local normal strain is defined as $\epsilon^\alpha = \delta_n^\alpha / l$ and the local shear strain is defined as $\gamma^\alpha = \delta_r^\alpha / l$.

The local stress paths in the six contact planes are plotted in Figure 9. Since the local stresses are related to the stresses applied to the overall assembly by the static hypothesis, the local stresses on each contact plane can be determined from the applied stresses, even though they are different on each plane. At any given moment, the maximum local shear stress is located in the 45° inclined plane. The less-inclined
planes (oriented more horizontally) carry normal forces larger than the more-inclined planes (oriented more vertically).

<figure><img src="image_9.png" /><figcaption>**Figure 9.** *Comparison of local stress paths for various inter-particle contact planes*</figcaption></figure>

In Figure 9, the circular points on each curve correspond to the moment when the peak stress point occurs in Figure 5. It indicates that all contact planes experience a reduction of the contact shear force after the occurrence of the peak stress. Since the local stress paths for all contact planes display a softening behaviour, it is of interest to know whether all contact planes experience instability at once.

<figure><img src="image_10.png" /><figcaption>**Figure 10.** *Comparison of second-order work in various inter-particle contact planes*</figcaption></figure>

In order to answer this question, we computed the second-order work for each individual plane. Note that in an undrained test, the overall volume change is zero, but the local normal strain does not need to remain constant ($d\epsilon^\alpha \neq 0$). Thus the
second-order work for the $\alpha^{th}$ contact plane, expressed in terms of local stress and local strain, is given by:

$$d^2W^{\alpha} = d\sigma^{\alpha}d\epsilon^{\alpha} + d\tau^{\alpha}d\gamma^{\alpha} \quad [26]$$

The second-order work for the six contact planes is plotted in Figure 10. For the two less-inclined contact planes at 17.9° and 28.3°, the second-order work is always positive, therefore, these two planes remain stable during the entire undrained triaxial test. On the contrary, for the four more-inclined contact planes oriented between 45° and 72°, the second-order work displays negative values but does not become negative simultaneously in each plane. Instability occurs first in the two contact planes with inclination angles 45° and 54.7°. After a few additional load steps, instability occurs subsequently in the contact planes with inclination angles 70.4° and 72°. Thus, even though the local stress paths in all the planes show a softening behaviour, only a part of the contact planes experience instability, but at different times.

<figure><img src="image_11.png" /><figcaption>**Figure 11.** Local stress-strain curves for inter-particle contact planes with inclined angles 17.9° and 28.3°</figcaption></figure>
The assembly instability occurs only when few contact planes (inclined at 45°-54.7°) start to become unstable. Even though the majority of the inter-particle planes are still stable at this time, the unstable condition of the few planes is severe enough for the overall assembly to experience instability. This can be seen from the negative sum of the second-order work for all the planes at this point.

To examine more deeply the behaviour of each individual contact plane, we plotted the local normal stress-strain curves and local shear stress-strain curves for four contact planes. The less-inclined planes (17.9°, 28.3°) are plotted in Figure 11 and the more-inclined planes (54.7°, 70.4°) are plotted in Figure 12. Initially, the specimen is isotropically consolidated to 0.3 MPa and then sheared in undrained condition. In the local normal stress-strain curves in Figures 11 and 12, the part corresponding to the isotropic consolidation is plotted as a dotted line up to 0.3 MPa. From that point on, the curve represents the local response during undrained shear loading. The isotropic consolidation is not displayed in the local shear stress-strain curves. The entire curve represents the local response during the undrained shear loading stage. The four circular points marked on each curve correspond to the four loading stages in Figure 5.

<figure><img src="image_12.png" /><figcaption>**Figure 12.** Local stress-strain curves for inter-particle contact planes with inclined angles 54.7° and 70.4°</figcaption></figure>
It can be observed, in the less-inclined planes, that the local normal strain expands as the normal stress is released. Thus, the second-order work due to the normal stress is positive ($d\sigma^\alpha d\varepsilon^\alpha > 0$). It can also be observed that the local shear stress increases monotonically with the local shear strain until the peak is reached. After the peak, the local shear stress decreases with the local shear strain, producing a positive second-order work ($d\tau^\alpha d\gamma^\alpha > 0$). Thus, the overall second-order work is always positive.

On the plane with an inclination angle of 70.4° (Figure 12), the normal behaviour is similar to the one corresponding to the less-inclined planes, i.e., the local normal strain expands as the normal stress is released ($d\sigma^\alpha d\varepsilon^\alpha > 0$). However, the behaviour is reversed on the plane with an inclination angle of 54.7°, since, when the local normal stress is released, the local normal strain contracts rather than expands, producing large negative second-order work $d\sigma^\alpha d\varepsilon^\alpha < 0$.

The local shear behaviour for the more-inclined planes is also reversed. After the peak stress is reached, the shear stress starts to decrease as the shear strain continues to increase, showing a softening pattern and producing a large negative second-order work $d\tau^\alpha d\gamma^\alpha < 0$. The overall second-order works for the more-inclined planes are negative. The behaviour is dominated by an instable condition.

# 5. Summary and Conclusion

According to Hill's condition, whether a material is stable or not depends not only on the current stress state but also on the direction of the stress increment. Thus, under specific loading conditions, loose sand can remain stable until the critical state failure line is reached. However, under other loading conditions, loose sand can succumb to instability at a shear stress level much lower than the critical state failure line. Therefore, the usual design concept in geotechnical engineering is likely to overestimate the soil strength, which can cause catastrophic collapse of earth slopes.

A micromechanical approach has been adopted for the analysis of this type of instability problem. The model considers the material as an assembly of particles. The stress-strain relationship for the assembly is determined by integrating the behaviour of the inter-particle contacts in all orientations. The inter-particle contact is assumed to have an elasto-plastic behaviour.

The constitutive model has been used to simulate undrained triaxial tests and constant-q tests. The comparison between experimental and predicted results has shown the capability of the model for obtaining the modes of instability at the assembly level. In particular, the position of the instability line, the same for these two types of tests in the p'-q plane, has been predicted with very good accuracy by the model simulations.
Analysis at the inter-particle contact level during an undrained triaxial test has shown that only a few inter-particle planes are unstable at the onset of the assembly instability. Afterwards, strain softening begins with a reduction of the overall shear stress, and progressively, more inter-particle contact planes become unstable.

## 6. References

Biarez J., Hicher, P.Y., *Elementary Mechanics of Soil Behaviour*, Balkema, 1994, p. 208.

Chang C.S., Hicher P.-Y., “An elastic-plastic model for granular materials with microstructural consideration,” *International Journal of Solids and Structures*, vol. 42, 2004, p. 4258-4277.

Chang C.S., Sundaram S.S., Misra A., “Initial moduli of particulate mass with frictional contacts”, *Int. J. for Numerical & Analytical Methods in Geomechanics*, John Wiley & Sons, vol. 13, n° 6, 1989, p. 626-641.

Christofferson J., Mehrabadi M.M., Nemat-Nassar S., “A micromechanical description on granular material behavior”, *ASME Journal of Applied Mechanics*, vol. 48, 1981, p. 339-344.

Chu J., Leong W.K., “Recent progress in experimental studies on instability of granular soil”, *International Workshop on Bifurcations and Instabilities in Geomechanics*, Labuz J.F., Drescher A. (eds.), Swets & Zeitlinger, Lisse, 2003, p. 175-191.

Darve F., Roguiez X., “Homogeneous bifurcation in soils”, *Localization and Bifurcation Theory for Soils and Rocks*, Adachi *et al*. (eds.), Balkema, Rotterdam, 1998, p. 43-50.

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

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

Gajo A., Piffer L., De Polo F., “Analysis of certain factors affecting the unstable behaviour of saturated loose sand”, *Mechanics of Cohesive-Frictional Materials*, vol. 5, 2000, p. 215-237.

Goddard J.D., “Material instability in complex fluids”, *Annual Review of Fluid Mechanics*, 2003, vol. 35, p. 113-133.

Hicher P.Y., “Experimental behavior of granular materials”, *Behavior of Granular Materials*, B. Cambou (ed.), Springer Wien New York, 1998, p. 1-97.

Hicher P-Y., Chang C.S., “An anisotropic non-linear elastic model for particulate materials”, *J. Geotechnical and Geoenvironmental Engrg.*, ASCE, vol. 132, n° 8, 2007.

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

Lade P.V., “Instability, shear banding, and failure in granular materials”, *International Journal of Solids and Structures*, vol. 39, n° 13-14, 2002, p. 3337-3357.

Lanier J., Block J.F., “Essais à volume constant réalisés sur presse tridimensionnelle”, *Greco Geomaterials Report*, 1989, p. 240-243.
Laouafa F., Darve F., “Modelling of slope failure by a material instability mechanism”,
*Computers and Geotechnics*, 2002, vol. 29, n° 4, p. 301-325.

Mindlin R.D., Deresiewicz H., “Elastic spheres in contact under varying oblique forces”,
*ASME Trans. J. Appl. Mech.*, vol. 20, 1953, p. 327-344.

Nova R., “Controllability of the incremental response of soil specimens subjected to arbitrary loading programs”, *J. Mechanical Behaviour of Materials*, 1994, vol. 5, n° 2, p. 193-201.

Nova R., Imposimato S., “Non-uniqueness of the incremental response of soil specimens under true-triaxial stress paths”, *Numerical Models in Geomechanics*, Pietruszczak & Pande (eds.), Balkema, 1997, p. 193-198.

Oda M., “Co-ordination number and its relation to shear strength of granular material”, *Soils and Foundations*, vol. 17, n° 2, 1977, p. 29-42.

Rothenburg L., Selvadurai A.P.S., “Micromechanical definitions of the Cauchy stress tensor for particular media”, *Mechanics of Structured Media*, Selvadurai A.P.S. (eds.), Amsterdam, Elsevier, 1981, 469-486.

Rudnicki J.W., Rice J., “Conditions for the localization of deformation in pressure sensitive dilatant materials”, *International Journal of Solids and Structures*, 1975, vol. 23, p. 371-394.

Sasitharan S., Robertson P.K., Sego D.C., Morgenstern N.R., “Collapse behavior of sand”,
*Canadian Geotechnical Journal*, vol. 30, n° 4, 1993, p. 569-577.

Schofield A.N., Wroth C.P., *Critical State Soil Mechanics*, London, England, McGraw-Hill, 1968.