# Experimental and numerical investigation of diffuse instability in granular materials using a microstructural model under various loading paths

Ali Daouadji, Pierre-Yves Hicher, Mohamad Jrad, Beena Sukumaran, Salim Belouettar

► To cite this version:

Ali Daouadji, Pierre-Yves Hicher, Mohamad Jrad, Beena Sukumaran, Salim Belouettar. Experimental and numerical investigation of diffuse instability in granular materials using a microstructural model under various loading paths. Geotechnique, 2013, 63 (5), pp.368-381. 10.1680/geot.10.P.121. hal-01501440

HAL Id: hal-01501440

https://hal.univ-lorraine.fr/hal-01501440

Submitted on 18 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
# Experimental and numerical investigation of diffuse instability in granular materials using a microstructural model under various loading paths

A. DAOUADJ, P.-Y. HICHER, M. JRAD, B. SUKUMARAN and S. BELOUETTAR

Depending on stress state and loading conditions, instability in loose sand can develop at a shear stress level much lower than the Mohr-Coulomb failure criterion. This type of failure mode, very different from strain-softening in dense sand, can lead to catastrophic collapse of earth structures. An appropriate constitutive model is therefore needed that can capture the correct failure modes and accurately predict the stress-strain response for both dense and loose sands. In this paper, a micromechanical approach developed earlier for the modelling of granular material behaviour is adopted. The constitutive model is utilised to simulate undrained triaxial, constant-q and proportional strain tests on loose sand, and numerical results are compared with experimental results. The model is capable of accurately reproducing the macroscopic experimental tests, and is also able to provide microstructural information. This enables the observation of instability at the particle scale, and provides a better understanding of the linkage between particle-scale instability and instability at the assembly level.

KEYWORDS: failure; liquefaction; numerical modelling

## INTRODUCTION

Instability of granular materials can be categorised into two main aspects: material instability (also known as intrinsic/ constitutive instability) and structural instability (e.g. Goddard, 2003). Localised failure 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 & Rice (1975) that the occurrence of localised failures can be predicted at a constitutive level. More recently, several authors have demonstrated that other modes of instability, such as diffuse failures, can occur within granular materials. Diffuse instability occurs when no strain localisation is noticed in the sample. Although strain localisation occurs even in loose samples (Mokni & Desrues, 1999), it appears at an axial strain of about 5% during a biaxial test, whereas the diffuse mode of instability occurs at about 1% axial strain during a biaxial or triaxial test under undrained conditions. Darve & Roguiez (1998), Laouafa & Darve (2002), Darve et al. (2004, 2007), Daouadji & Hicher (2010), Daouadji et al. (2011) and Nicot et al. (2011) have studied this mode of instability using a macroscopic phenomenological approach based on incrementally linear and non-linear models.

Nova (1994) has developed similar concepts, and has demonstrated that any non-associated elasto-plastic model could exhibit domains of instability inside the plastic limit. In this paper, a study of material instability, using a micro-mechanical approach developed by Chang & Hicher (2004) is presented. Under this framework, the stress-strain relationship for a granular assembly is determined by integrating

the behaviour of the interparticle 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.

Experimental evidence supports 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*. The same can be said for constant-q and proportional strain tests. The former stress path can simulate the loading condition of a soil element on sloping ground when a progressive increase in pore pressure occurs, and the latter allows the influence of drainage conditions on the mechanical behaviour of saturated sands to be studied.

Numerical results obtained using the micromechanical-based model are compared with experimental results obtained on loose Hostun sand along three different loading paths: undrained triaxial tests, constant-q and proportional strain tests. It will be shown that the model is capable of accurately reproducing the observed mechanical behaviour. The analysis of the behaviour of loose sand using the proper control parameter is adequate to study the possible instabilities that may occur during loading. In addition, the local response is also studied along various orientations of branch vectors using the microstructural model.

## STRESS-STRAIN MODEL BASED ON A MICROMECHANICAL APPROACH

The microstructural CH model developed by Chang & Hicher (2004) is briefly described. In this model, a granular material is viewed as a collection of particles. The deformation of a representative elemental volume of the material is obtained 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* and *t*. The vector *n* is outward normal to the contact plane, and vectors *s* and *t* are on the contact plane.
Interparticle behaviour

Elastic stiffness. The contact stiffness of a contact plane α includes normal stiffness, k_n^α, and shear stiffness, k_t^α. The elastic stiffness tensor is defined by

$$ \mathbf{f}_i^\alpha = k_{ij}^{\alpha e} \delta_j^{\alpha e} \qquad (1) $$

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

$$ 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 (2) $$

The value of the stiffness for two elastic spheres can be estimated from the Hertz–Mindlin formulation (Mindlin & Deresiewicz, 1953). For sand grains, a revised form was adopted (Chang et al., 1989), which is given by

$$ k_n = k_{n0} \left( \frac{f_n}{G_g l^2} \right)^n \qquad (3) $$

$$ k_t = k_{t0} \left( \frac{f_n}{G_g l^2} \right)^n \qquad (4) $$

where G_g is the elastic modulus for the grains, f_n is the contact force in the normal direction, l is the branch length between two particles, and 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) as

$$ F(f_r, \kappa) = f_r - f_n \kappa (\delta_r^p) = 0 \qquad (5) $$

where κ(δ_r^p) is a hardening/softening parameter. The shear force f_r and the rate of plastic sliding δ_r^p are defined as

$$ f_r = \sqrt{f_s^2 + f_t^2} \qquad (5a) $$

and

$$ \delta_r^p = \sqrt{(\delta_s^p)^2 + (\delta_t^p)^2} \qquad (5b) $$

The hardening function is defined by a hyperbolic curve in the κ–δ_r^p plane, where φ_p is the peak friction angle, and k_{p0} is the plastic modulus.

$$ \kappa = \frac{k_{p0} \tan \phi_p \delta_r^p}{|f_n| \tan \phi_p + k_{p0} \delta_r^p} \qquad (6) $$

Plastic flow rule. Plastic sliding occurs in the contact plane, accompanied by an upward or downward movement that provokes a shear-induced dilation/contraction of the whole specimen. The dilatancy effect can be described by

$$ \frac{d\delta_n^p}{d\delta_r^p} = \frac{f_r}{f_n} - \tan \phi_0 \qquad (7) $$

where the material constant φ_0 is equal to the interparticle friction angle φ_μ, for most cases. 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 (equation (7)). Thus the flow rule is non-associated.

Elasto-plastic relationship. With the elements discussed above, the incremental force–displacement relationship of the interparticle contact can be obtained. Including both elastic and plastic behaviour, this relationship is given by

$$ \dot{f}_i^\alpha = k_{ij}^{\alpha p} \delta_j^\alpha \qquad (8) $$

A detailed expression for the elasto-plastic stiffness tensor can be derived in a straightforward manner from the yield function and flow rule, and therefore is not provided here. More details can be obtained from earlier publications (Chang & Hicher, 2004; Hicher & Chang, 2007).

Influence of void ratio on φ_p

The resistance against sliding in a contact plane is dependent on the degree of interlocking with 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 \qquad (9) $$

where m is a material constant (Biarez & 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 interparticle friction angle φ_p is greater than the internal friction angle, φ_μ. When the packing structure dilates, the degree of interlocking and the apparent friction angle are reduced, which results in a strain-softening phenomenon. For loose packing, the apparent friction angle φ_p is smaller than the internal friction angle φ_μ, and increases during 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 (Schofield & Wroth, 1968)

$$ e_c = \Gamma - \lambda \log p' \qquad (10a) $$

or

$$ e_c = e_{\text{ref}} - \lambda \log \left( \frac{p'}{p_{\text{ref}}} \right) \qquad (10b) $$

where Γ and λ 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.

Micro–macro relationship

The stress–strain relationship for an assembly can be determined by integrating the behaviour of interparticle 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 from the contact forces and branch vectors for all contacts (Christofferson et al., 1981; Rothenburg & Selvadurai, 1981), as

$$ \dot{\sigma}_{ij} = \frac{1}{V} \sum_{\alpha=1}^{N} f_j^\alpha I_i^\alpha \qquad (11) $$

The mean force on the contact plane for each orientation is

$$ \dot{f}_j^\alpha = \dot{\sigma}_{ij} \mathbf{A}_{ik}^{-1} I_k^\alpha V \qquad (12) $$

where the branch vector I_k^α is defined as the vector joining the centres of two particles, and the fabric tensor is defined as
$$ \mathbf{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 intercluster 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 (13) $$

From equations (12) and (13), the relation between the global strain and interparticle displacement is obtained

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

### Stress-strain relationship

From equations (8), (11) and (14), the relationship between stress increment and strain increment can be obtained as

$$ \dot{u}_{i,j} = \mathbf{C}_{ijmp} \dot{\sigma}_{mp} \qquad (15) $$

where

$$ \mathbf{C}_{ijmp} = \mathbf{A}_{ik}^{-1} \mathbf{A}_{mn}^{-1} V \sum_{\alpha=1}^{N} (k_{jp}^{\alpha p})^{-1} l_{k}^{\alpha} l_{n}^{\alpha} $$

When the contact number N is sufficiently large in an isotropic packing, the summation of the compliance tensor in equation (15) and the summation of the fabric tensor in equation (12) can be written in integral form as

$$ \mathbf{C}_{ijmp} = \mathbf{A}_{ik}^{-1} \mathbf{A}_{mn}^{-1} \frac{NV}{2\pi} \int_{0}^{\pi/2} \times \int_{0}^{2\pi} k_{jp}^{\alpha p} (\gamma, \beta)^{-1} l_{k}(\gamma, \beta) l_{n}(\gamma, \beta) \sin\gamma \, d\gamma \, d\beta \qquad (16) $$

and

$$ \mathbf{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 (17) $$

The integration of equations (16) and (17) in a spherical coordinate system can be carried out numerically by using Gauss integration points over the surface of the sphere.

## EXPERIMENTAL EVIDENCE OF INSTABILITY IN SAND

### Undrained triaxial tests

Undrained triaxial compression tests on Hostun sand with various initial relative densities, $D_r$, are presented in Fig. 1 (Hicher, 1998). These results show two distinctive trends corresponding to a contractive or dilative behaviour. For loose sand (small values of $D_r$), a maximum peak strength is reached in the $q-\epsilon_1$ plane, which is then 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 demonstrated later. 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 to softening still occurs after the peak, but it is followed by an increase of the deviatoric stress up to the ultimate strength, corresponding to the state at large deformations.

<figure><img src="image_1.png" /><figcaption>Fig. 1. Experimental results for undrained triaxial tests on Hostun sand with various densities during displacement-controlled tests: (a) stress-strain curves; (b) stress paths (redrawn after Hicher, 1998)</figcaption></figure>

(Carrera et al., 2011). For dense sand (high values of $D_r$) the material is strongly dilative, 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 observed in the $p'-q$ plane: a continuous decrease of the mean effective stress is observed for strongly contractive materials, and a decrease followed by an increase of the mean effective stress up to the critical state is observed for dilatant materials.

The position of the peak in the stress plane $p'-q$ 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 distant from the failure line for loose sand and closer to the failure line for denser sand.

### Constant-q tests

The constant-q test consists of shearing the specimen to a prescribed stress ratio along a drained compression triaxial path, and then decreasing the mean effective stress while keeping the deviatoric stress constant. This stress path can simulate the loading condition of a soil element on sloping ground 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 (for a brief review see Darve et al., 2007). The authors have presented experimental results on Hostun sand, and have 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 sand (CSD) are presented in Fig. 2. After an isotropic consolidation stage to a desired initial effective mean pressure $p'_0$, a drained triaxial compression test is performed on the sample at two prescribed values of deviatoric stress q (CSD1 and CSD2). Then, while keeping q constant, the mean effective stress $p'$ is decreased by increasing the pore water pressure while maintaining the vertical and radial total stresses applied at the boundaries of the sample constant. This pore water increase is due to an imposed volume change on the sample, so the test is controlled by the deviatoric stress q and the volumetric strain $\varepsilon_v$. At a given point during this test stage, the axial strain rate starts to increase very rapidly, illustrated by a jump in axial strain as seen in Fig. 3, at which point the deviatoric stress can no longer be kept constant. The test, corresponding to a given constant value of q, is no longer controllable (as defined by Nova, 1994) in the sense that the imposed loading programme cannot be maintained. It will be shown in the next subsection that this point corresponds to a loss of stability, since any small change of one control variable will lead to a catastrophic failure. At collapse, the axial strain is 1.27% for CSD2, whereas this value is only 0.78% for CSD1.

### Proportional strain tests

Proportional strain path tests are performed to study the influence of drainage conditions on the mechanical behaviour of saturated loose sands (Chang et al., 2009). These strain paths reproduce the in situ conditions of a soil element when subjected to an excessive, partial or expansive drainage (Sivathayalan & Logeswaran, 2007) condition.

The loading programme for such tests is

$$ \left\{ \begin{array}{l} d\varepsilon_1 > 0 \\ d\varepsilon_2 = d\varepsilon_3 \\ \alpha = \frac{d\varepsilon_v}{d\varepsilon_1} \end{array} \right. \qquad (18) $$

where axial compression and volume outflow are positive. A positive value of $\alpha$ corresponds to a contractive test; a negative value corresponds to an expansive or dilative test, and a nil value is the conventional strain-controlled undrained test. From the definition of $\alpha$, it is clear that $\alpha = 1$ corresponds to a one-dimensional test.

<figure><img src="image_2.png" /><figcaption>Fig. 2. Experimental loading paths imposed during constant-q tests on loose Hostun sand: $p'_0 = 300$ kPa, $q_1 = 55$ kPa (CSD1), $q_2 = 120$ kPa (CSD2)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 3. Experimental results of constant-q tests on loose Hostun S28 sand: $p'_0 = 300$ kPa, $q_1 = 55$ kPa (CSD1), $q_2 = 120$ kPa (CSD2). Onset of collapse highlighted by the black dot</figcaption></figure>

Four values of $\alpha$ were chosen: 0.5, 0.26, 0 and -0.47. The first two values induce contractive behaviour, the third value induces an isochoric condition, and the fourth imposes dilative behaviour. These four applied conditions are shown in Fig. 4.

The stress-strain curves of loose Hostun sand subjected to these four proportional strain tests preceded by isotropic consolidation to an effective mean stress of 100 kPa are shown in Fig. 5. The undrained test ($\alpha = 0$) presents a peak in the $q-\varepsilon_1$ and $q-p'$ planes, which occurs at an axial strain of about 0.4%. This peak is followed by a total loss of effective stress and a collapse of the sample. Such a collapse was also observed for $\alpha = -0.47$, as a total loss of effective stress was noticed. However, the values of the corresponding deviatoric stress and axial strain are lower than observed during the undrained tests. Indeed, in the case of such an expansive or dilative test, the axial strain corresponding to the peak is about 0.1%.

A peak is also noticeable for $\alpha = 0.26$, but no collapse was observed, as the material recovers strength after a decrease of the deviatoric stress. It is noteworthy that the curves corresponding to $\alpha = 0.5$ do not present any peak. For this test, one should expect no instability during this
<figure><img src="image_4.png" /><figcaption>Fig. 4. Four proportional strain paths applied to loose Hostun sand: two tests are contractive (α = 0.26 and 0.5), one is isochoric (α = 0) and one is expansive (α = -0.47)</figcaption></figure>

loading path for Hostun loose sand. However, the analysis in terms of stability will give us more insight into the mechanical behaviour of loose Hostun sand during proportional strain tests.

## ANALYSIS OF EXPERIMENTAL RESULTS IN TERMS OF STABILITY CONDITION

### Analysis of undrained compression and constant-q tests

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

$$ 
\begin{aligned}
d^2 W &= d\sigma'_{ij} d\varepsilon_{ij} \\
&= dq d\varepsilon_d + dp' d\varepsilon_v \\
&> 0
\end{aligned} 
\quad (19) $$

where under axisymmetric conditions, and for saturated sands, the deviatoric stress, $q = \sigma_1 - \sigma_3$; the effective mean stress, $p' = \frac{1}{3}(\sigma'_1 + 2\sigma'_3)$; and the deviatoric strain, $\varepsilon_d = \frac{2}{3}(\varepsilon_1 - \varepsilon_3)$.

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 (19) can be rearranged to give (Darve et al., 2004)

$$ d^2 W = dq d\varepsilon_1 + d\sigma'_3 d\varepsilon_v \quad (20) $$

Thus, as the control variables during the load-controlled tests are $dq$ and $d\varepsilon_v$ (as mentioned earlier, a volume of water needs to be injected inside the sample in order to reach the targeted pore pressure value), the constitutive relation linking the stress increments to the strain increments can be rearranged to give a generalised mixed relation between generalised incremental stresses and generalised incremental strains as

$$ 
\begin{aligned}
\left\{ \begin{matrix} dq \\ d\varepsilon_v \end{matrix} \right\} &= [\mathbf{P}] \left\{ \begin{matrix} d\varepsilon_1 \\ d\sigma'_3 \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\varepsilon_1 \\ d\sigma'_3 \end{matrix} \right\}
\end{aligned} 
\quad (21) $$

where $E_i$ are pseudo Young's moduli and $v_i$ are pseudo Poisson's ratios. As no volumetric strain is allowed during undrained tests, equation (21) indicates that instability will occur only if $dq = 0$, which corresponds to the peak shear stress (Fig. 1). Instability can only occur if the determinant of the matrix $[\mathbf{P}]$ of the relation given by equation (21) reduces to zero.

$$ \det(\mathbf{P}) = 1 - 2v_{31} + 2(1 - v_{33} - v_{13})\frac{E_1}{E_3} = 0 \quad (22) $$

On the other hand, as no variation of the shear stress is imposed during constant-q tests ($dq = 0$), instability occurs during this test only if $d\varepsilon_v = 0$, as $d\varepsilon_1 > 0$ and $d\sigma'_3 < 0$ are imposed by the loading programme (Figs 3 and 6). This condition coincides with that obtained for undrained tests

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

<figure><img src="image_6.png" /><figcaption>Fig. 5. Experimental results on Hostun sand: deviatoric stress plotted against (a) axial strain and (b) effective mean stress for four proportional strain tests</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Experimental results of constant-q tests on loose Hostun S28 sand: $p'_{0} = 300$ kPa, $q_{1} = 55$ kPa (CSD1), $q_{2} = 120$ kPa (CSD2). Vanishing of the second-order work corresponds to the black dots</figcaption></figure>

(equation (22)). In other words, 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 that corresponding to the deviatoric stress peak of the undrained

tests. This condition corresponds to a mobilised friction angle equal to 16°, which is much lower than the friction angle at critical state, equal to 31° (Daouadji *et al.*, 2010). It should be emphasised that during the displacement-
controlled tests under undrained conditions, the tests are kinematically controlled so that large deformations can be induced. Therefore mixed control variables have to be used to trigger diffuse instabilities in granular materials with an effective collapse of the sample and no localisation pattern. The experimental second-order work is computed using equations (19) and (20) and is plotted against axial strain for CSD1 and CSD2. It is clearly shown that the loss of controllability coincides with the vanishing of the second-order work, which corresponds to a diffuse failure.

### Analysis of proportional strain tests

Darve et al. (2004) have used an equivalent definition of the loading programme when numerically checking the occurrence of instability conditions using an incrementally non-linear constitutive relation. It is expressed as

$$ \left\{ \begin{array}{l} d\varepsilon_1 > 0 \\ d\varepsilon_2 = d\varepsilon_3 \\ d\varepsilon_1 + 2Rd\varepsilon_3 = 0 \end{array} \right. \qquad (23) $$

where $\alpha$ and $R$ are related by $\alpha = 1 - (1/R)$.

As in the previous subsection, equation (19) can be rewritten using the conjugate variables with respect to energy as

$$ d^2 W = \left( d\sigma'_{1} - \frac{d\sigma'_{3}}{R} \right) d\varepsilon_{1} + \left( d\varepsilon_{1} + 2Rd\varepsilon_{3} \right) \frac{d\sigma'_{3}}{R} \qquad (24) $$

As, during the proportional strain tests, $d\varepsilon_1 + 2Rd\varepsilon_3$ is nil, the second-order work is therefore given by

$$ d^2 W = \left( d\sigma'_{1} - \frac{d\sigma'_{3}}{R} \right) d\varepsilon_{1} \qquad (25) $$

which does not correspond to $dq \cdot d\varepsilon_1$ when $R \neq 1$.

As, during the proportional strain test, the axial strain is increased (equations (18) and (23)), equation (25) vanishes only if $d\sigma'_{1} - (d\sigma'_{3}/R)$ is nil, that is, when $\sigma'_{1} - (\sigma'_{3}/R)$ reaches a maximum as $d\varepsilon_1 + 2Rd\varepsilon_3 = 0$ during the test.

The second-order work is computed using the experimental data, and is plotted in Fig. 7 against the axial strain. In this figure the deviatoric stress $q$, the control parameter $d\sigma'_{1} - (d\sigma'_{3}/R)$ and the experimental second-order work $d^2 W$ are plotted against axial strain for $R = 0.68$, 1, 1.36 and 2 using equation (19) and equation (25), and both give exactly the same result.

A peak of the deviatoric stress exists for the first three values of $R$, but not for the last one, and the common interpretation is that there is no instability when $R = 2$. However, using equation (19) and equation (25), $d^2 W$ vanishes during the four tests, indicating that instability occurs during the expansive or compressive tests up to $R = 2$. It is found that the vanishing point occurs before the deviatoric stress peak for $R = 0.68$, at the peak for $R = 1$, and after the peak $R = 1.36$ and 2. Interestingly, the experimental second-order work always vanishes at

<figure><img src="image_8.png" /><figcaption>Fig. 7. Experimental results during four proportional triaxial tests on Hostun sand for (a) $R = 0.68$, (b) $R = 1$, (c) $R = 1.36$ and (d) $R = 2$, for which second-order work vanishes respectively after peak, at peak and before peak of deviatoric stress. It always vanishes at peak of control parameter even if no peak is noticeable for $R = 2$ in $q-p'$ plane</figcaption></figure>
the peak of the control parameter even for R = 2 (Fig. 7). Therefore this quantity is a good indicator of instabilities in sands, and the classical analysis based on the presence of a peak in the deviatoric stress must be re-examined in terms of the presence of a peak of the control parameter.

## MODEL PERFORMANCES

### Numerical results of undrained triaxial tests

Figure 8 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 Fig. 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 Fig. 9 for a test on loose Hostun sand with an initial confining stress $p'_0 = 300$ kPa. For undrained conditions ($d\epsilon_v = 0$), the second-order work is reduced to $d^2W = dq \cdot d\epsilon_d$, where $d\epsilon_d$ is the deviatoric strain. Since the deviatoric strain increases continuously, $d\epsilon_d$ is always positive, and the second-order work can become non-positive if and only if $dq \le 0$ (i.e. a decrease in $q$). Fig. 9 shows that instability begins at the peak shear stress. 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 $\epsilon_1$ and not the deviatoric stress $q$.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Model predictions for undrained triaxial tests on Hostun sand with various relative densities</figcaption></figure>

### Numerical results of constant-q tests

The parameters in Table 1 were used to predict the results of constant-q tests on loose Hostun sand. These parameters were determined by Chang and Hicher for earlier studies on Hostun sand (Chang & Hicher, 2004; Hicher & Chang, 2007). The predicted and measured results for the confining stress $p'_0 = 300$ kPa are presented in Fig. 10. The initial part of the $p'-\epsilon_v$ curve shows that, as the mean stress $p'$ decreases, the volume strain increases. This trend continues until a certain point, when the volumetric strain starts to decrease. For constant-q tests ($dq = 0$), according to equation (19), the second-order work is reduced to $d^2W = dp' \cdot d\epsilon_v$. Since the mean stress progressively decreases (i.e. $dp' < 0$), the second-order work becomes negative or nil if and only if $d\epsilon_v \le 0$ (i.e. if the volume contracts). Thus the onset of instability corresponds to the peak of the $p'-\epsilon_v$ curve, which is well reproduced by the model simulation. Using an approach similar to that used for the undrained tests discussed in the previous subsection, 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$.

### Numerical results of proportional strain tests

Using the same sets of model parameters given in Table 1, the model was used to predict the four proportional strain tests described earlier. For that purpose, equation (23) is used and the predicted results obtained for the confining stress $p'_0 = 100$ kPa are presented in Fig. 11. It is found that the model is capable of reproducing the trends of the mechanical behaviour of sands even under such a loading programme. Indeed, the peaks of the deviatoric stresses are determined accurately for $R = 0.68$, 1, 1.36, and no peak is reached for $R = 2$. Moreover, the axial strains corresponding to these maximum values are well simulated.

Using the same approach as described earlier, the deviatoric stress, the control parameter and the second-order work are plotted against the axial strain for the four values of R (Fig. 12). With the numerical results, it is also found that the axial strain corresponding to the vanishing of $d^2W$ is located before the peak for $R = 0.68$, at the peak for $R = 1$, and after the peak $R = 1.36$ and 2. Interestingly, the numerical second-order work always vanishes at the peak of the control parameter, even for $R = 2$.

### Summary of condition of instability for loose Hostun sand

Instability conditions obtained by numerical simulations of undrained compression and constant-q tests are plotted in the $p'-q$ plane (Fig. 13) together with the experimental results presented earlier (Figs 1 and 2). One can see that the model is capable of accurately predicting the condition of instability associated with these two types of test. The position of the instability line determined by the model simulations is in very good agreement with those obtained experimentally. Thus both loading programmes can be seen as a validation for the study of diffuse instabilities using different loading paths, as they need to fulfil the same condition (equation (22)) for its onset.

It should be emphasised that the $q-p'$ plane is relevant to observing the occurrence of instability in undrained triaxial compression tests, which correspond to a value of $R = 1$. However, for expansive or compressive loadings, the control parameter is $d\sigma'_1 - (d\sigma'_3/R)$, so the analysis has to be performed in the $[\sigma'_1 - (\sigma'_3/R)] - \epsilon_1$ plane.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Predicted second-order work and stress–strain curves for undrained triaxial test on loose Hostun sand ($p'_0 = 300$ kPa)</figcaption></figure>

<table><caption>Table 1. Model parameters for fine Hostun sand (CH model)</caption><thead><tr><th>e<sub>ref</sub></th><th>p<sub>ref</sub>: MPa</th><th>λ</th><th>φ<sub>μ</sub>: degrees</th><th>φ<sub>0</sub>: degrees</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>

LOCAL PLANE BEHAVIOUR

In order to gain insight into the material behaviour at the interparticle level, several contact planes oriented in six different directions indicated by the branch vectors have been selected, as shown in Fig. 14. The six directions are defined by the angle $\theta$ between the vector $n$ perpendicular to the contact plane and the vertical axis. A larger value of $\theta$ indicates a more inclined contact plane.

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

The local stress paths in the six contact planes are plotted in Fig. 15. 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 along the 45° inclined plane. The less-inclined

planes (oriented more horizontally) carry normal forces larger than the more-inclined planes (oriented more vertically).

In Fig. 15, the points indicated on each curve correspond to the stress state when the peak stress shown in Fig. 9 occurs during the undrained test. This 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.

In order to answer this question, the second-order work for each individual plane is computed. 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^a \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^2 W^a = d\sigma^a d\epsilon^a + d\tau^a d\gamma^a \quad (26)$$

The second-order work for the six contact planes is plotted in Fig. 16. 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. By contrast, for the four more-inclined contact planes oriented between 45° and 72°, the second-order work displays negative values, but this does not happen simultaneously along all planes. Instability occurs first in the two contact planes with inclination angles of 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 few of the contact planes experience instability, and at different load steps.
<figure><img src="image_11.png" /><figcaption>Fig. 10. Comparison of predicted and measured results for constant-q test on Hostun sand ($p'_{0} = 300$ kPa)</figcaption></figure>

The assembly instability occurs only when a few contact planes (inclined at 45°–54.7°) start to become unstable. Even though the majority of the interparticle planes are still stable at this time, the unstable condition along these few planes is severe enough for the overall assembly to experience instability. This can be seen from the negative cumulative sum of the second-order work for all the planes at this point.

To examine the behaviour of each individual contact plane more closely, the local normal stress–strain curves and local shear stress–strain curves for four contact planes are plotted. The less-inclined planes (17.9°, 28.3°) and the more-inclined planes (54.7°, 70.4°) are plotted in Fig. 17. Initially, the specimen is isotropically consolidated to 300 kPa and then sheared in undrained condition. In the local normal stress–strain curves in Fig. 10, the isotropic consolidation phase is plotted as a dotted line up to 300 kPa. 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 shown in Fig. 9.

It can be observed, in the less-inclined planes, that the

<figure><img src="image_12.png" /><figcaption>Fig. 11. Numerical results on Hostun sand: deviatoric stress plotted against (a) axial strain and (b) effective mean stress for four proportional strain tests</figcaption></figure>

local normal strain expands as the normal stress is released. Thus the second-order work due to the normal stress is positive ($d\sigma^{a}d\varepsilon^{a} > 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^{a}d\gamma^{a} > 0$). Thus the overall second-order work is always positive.

On the plane with an inclination angle of 70.4° (Fig. 17), the normal behaviour is similar to that corresponding to the less-inclined planes: that is, the local normal strain expands as the normal stress is released ($d\sigma^{a}d\varepsilon^{a} > 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^{a}d\varepsilon^{a} < 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^{a}d\gamma^{a} < 0$). The overall second-order work for the more-inclined planes is negative. The behaviour is therefore dominated by an unstable condition.
<figure><img src="image_13.png" /><figcaption>Fig. 12. Numerical results during four proportional triaxial tests on Hostun sand for (a) R = 0.68, (b) R = 1, (c) R = 1.36 and (d) R = 2 for which second-order work vanishes respectively after peak, at peak and before peak of deviatoric stress, but always at peak of control parameter</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Comparison of predicted and measured onset of instability for loose Hostun sand determined from undrained triaxial and constant-q tests for various initial effective mean stresses</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Orientations of branch vectors for interparticle contact planes selected in analysis</figcaption></figure>

## SUMMARY AND CONCLUSION

According to Hill's condition, material stability is dependent 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 below the critical state failure line. Therefore the usual design concept used in geotechnical engineering is likely to overestimate the soil strength for such loading conditions, which can lead to catastrophic collapse of earth slopes.

Experimental evidence from undrained triaxial tests and
<figure><img src="image_16.png" /><figcaption>Fig. 15. Comparison of local stress paths for various interparticle contact planes</figcaption></figure>

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

<figure><img src="image_18.png" /><figcaption>Fig. 16. Comparison of second-order work in various interparticle contact planes</figcaption></figure>

constant-q tests has been presented to demonstrate diffuse failures in loose sands. Moreover, proportional strain tests have also been carried out to highlight the influence of the loading paths under compressive, isochoric and expansive conditions. For such test conditions, it is found that the second-order work vanishes after the peak, at the peak, or even before the peak. However, the second-order work vanishes exactly at the peak of the relevant control parameter. Moreover, it is clearly shown that, even if no deviatoric peak is noticeable during a compressive test, the second-order work may vanish, and will correspond to the peak of the control parameter.

A micromechanical approach has also been adopted for 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 interparticle contacts in all orientations. The interparticle contact is assumed to be elasto-plastic.

This constitutive model has been used to simulate undrained triaxial tests at various relative densities, and constant-q and strain-proportional tests on loose Hostun sand. The comparison between experimental and predicted results has shown that the model is capable of capturing diffuse modes of instability at the assembly level for sands. Specifically, the position of the instability line has also been predicted with very good accuracy by the model simulations, and has a much shallower slope than the critical state line.

The behaviour at the particulate level is also examined using the microstructural model along various orientations of branch vectors. Analysis at the interparticle contact level during an undrained triaxial test has shown that only a few interparticle planes are unstable at the onset of instability. Afterwards, strain-softening begins with a reduction in the overall shear stress, and progressively more interparticle contact planes become unstable. This demonstrates that not all the particles in the assembly are unstable at the onset of instability in the assemblage of particles.

### ACKNOWLEDGEMENTS

This work has been partly done during a Mobility In of Dr Daouadji (FNR/11/AM2c/26) founded by Fonds National de la Recherche Luxembourg (FNR), which is gratefully acknowledged. This material is based upon work supported by the National Science Foundation under Grant CMMI-0959415, which is gratefully acknowledged by Beena Sukumaran.

### NOTATION

<table><tr><td><b>A</b><sub>ik</sub>, <b>A</b><sub>mn</sub></td><td>fabric tensor</td></tr><tr><td><b>C</b><sub>ijmp</sub></td><td>compliance tensor</td></tr><tr><td><i>D</i><sub>t</sub></td><td>relative density</td></tr><tr><td><i>E</i><sub>i</sub></td><td>pseudo Young's moduli</td></tr><tr><td><i>e</i></td><td>packing void ratio</td></tr><tr><td><i>e</i><sub>c</sub></td><td>critical void ratio</td></tr><tr><td><i>e</i><sub>ref</sub></td><td>reference void ratio</td></tr><tr><td><b><i>f</i></b><sup><i>a</i></sup><sub><i>i</i></sub></td><td>elastic stiffness tensor</td></tr><tr><td><b><i>f</i></b><sup><i>a</i></sup><sub><i>j</i></sub></td><td>mean force on contact plane</td></tr><tr><td><i>f</i><sub>n</sub></td><td>contact force in normal direction</td></tr><tr><td><i>f</i><sub>t</sub></td><td>shear force</td></tr><tr><td><i>f</i><sub>s</sub></td><td>contact force in direction <i>s</i> in sliding plane</td></tr><tr><td><i>f</i><sub>t</sub></td><td>contact force in direction <i>t</i> in sliding plane</td></tr><tr><td><i>G</i><sub>g</sub></td><td>elastic modulus for grains</td></tr><tr><td><i>k</i><sup><i>ae</i></sup><sub>ij</sub>, <i>k</i><sup><i>ab</i></sup><sub>jp</sub></td><td>component of the elastic local tensor in the plane of orientation <i>a</i></td></tr><tr><td><i>k</i><sub>n</sub>, <i>k</i><sub>n0</sub></td><td>material constants</td></tr><tr><td><i>k</i><sup><i>a</i></sup><sub>n</sub></td><td>normal stiffness</td></tr><tr><td><i>k</i><sub>p0</sub></td><td>plastic modulus</td></tr><tr><td><i>k</i><sub>t</sub>, <i>k</i><sub>t0</sub></td><td>material constants</td></tr><tr><td><i>k</i><sup><i>a</i></sup><sub>t</sub></td><td>shear stiffness</td></tr><tr><td><i>l</i>, <i>l</i><sub>i</sub>, <i>l</i><sub>k</sub>, <i>l</i><sub>n</sub></td><td>branch length between two particles</td></tr><tr><td><i>I</i><sup><i>a</i></sup><sub>k</sub></td><td>branch vector joining centres of two particles</td></tr><tr><td><i>m</i></td><td>material constant</td></tr><tr><td><i>N</i></td><td>contact number</td></tr><tr><td><b><i>n</i></b>, <i>s</i>, <i>t</i></td><td>orthogonal unit vectors</td></tr><tr><td><i>p</i></td><td>mean stress</td></tr><tr><td><i>p</i>'</td><td>effective mean stress</td></tr><tr><td><i>p</i>'<sub>0</sub></td><td>initial effective mean pressure; initial confining stress</td></tr><tr><td><i>p</i><sub>ref</sub></td><td>reference mean stress</td></tr><tr><td><i>q</i></td><td>deviatoric stress</td></tr><tr><td><i>q</i><sub>1</sub>, <i>q</i><sub>2</sub></td><td>particular values of <i>q</i></td></tr><tr><td><i>R</i></td><td>constant value during a proportional strain test</td></tr></table>
<figure><img src="image_19.png" /><figcaption>Fig. 17. Local shear stress against local normal stress curves for interparticle contact planes with inclined angles (a) 7.9°; (b) 28.3°; (c) 54.7°; (d) 70.4°</figcaption></figure>
$\dot{u}_{j,i}$ strain increment

$V$ volume

$d^2W$ second-order work

$\alpha$ contact plane in the orientation $\alpha$

$\beta$ angle in local coordinate at the inter-particle contact

$\Gamma$ material constant

$\gamma$ angle in local coordinate at the inter-particle contact

$\gamma^a$ local shear strain

$\delta^a_j$ component of the elastic displacement

$\delta^p_r$ rate of plastic sliding

$\delta^p_s$ component of the tangential plastic displacement in direction $s$

$\delta^p_t$ component of the tangential plastic displacement in direction $t$

$\delta^a_n$ normal plastic displacement

$\delta^a_r$ rate of plastic sliding on $\alpha$ plane

$\varepsilon_d$ deviatoric strain

$\varepsilon_v$ volumetric strain

$\varepsilon^a$ local normal strain

$\varepsilon_1$ axial strain; vertical strain

$\varepsilon_2, \varepsilon_3$ radial (lateral) strains

$\theta$ angle between the vector $n$ perpendicular to the contact plane and the vertical axis

$\kappa(\delta^p_r)$ hardening/softening parameter

$\lambda$ material constant

$v_i$ pseudo Poisson's ratios

$\dot{\sigma}_{ij}$ stress increment

$\sigma^a$ local normal stress

$\sigma_1, \sigma_3$ axial and radial total stress respectively

$\sigma'_1, \sigma'_3$ axial and radial effective stress respectively

$\tau^a$ local shear stress

$\phi_p$ peak friction angle

$\phi_\mu$ interparticle friction angle

$\phi_0$ material constant

## REFERENCES

Biarez, J. & Hicher, P. Y. (1994). *Elementary mechanics of soil behaviour*, p. 208. Rotterdam, the Netherlands: Balkema.

Carrera, A., Coop, M. & Lancellotta, R. (2011). Influence of grading on the mechanical behaviour of Stava tailings. *Géotechnique* **61**, No. 11, 935–946, http://dx.doi.org/10.1680/geot.9. P.009.

Chang, C. S. & Hicher, P.-Y. (2004). An elastic-plastic model for granular materials with microstructural consideration. *Int. J. Solids Struct.* **42**, No. 14, 4258–4277.

Chang, C. S., Sundaram, S. S. & Misra, A. (1989). Initial moduli of particulate mass with frictional contacts. *Int. J. Numer. Anal. Methods Geomech.* **13**, No. 6, 626–641.

Chang, C. S., Hicher, P.-Y. & Daouadji, A. (2009). Investigating instability in granular materials by means of a micro-structural model. *Eur. J. Environ. Civ. Engng* **13**, No. 3, 167–186.

Christofferson, J., Mehrabadi, M. M. & Nemat-Nassar, S. (1981). A micromechanical description on granular material behavior. *ASME J. Appl. Mech.* **48**, No. 2, 339–344.

Daouadji, A. & Hicher, P. Y. (2010). An enhanced constitutive model for crushable granular materials. *Int. J. Numer. Anal. Methods Geomech.* **34**, No. 6, 555–580.

Daouadji, A., Al Gali, H., Darve, F. & Zeghloul, A. (2010). Instability in granular materials: an experimental evidence of diffuse mode of failure for loose sands. *J. Engng Mech. ASCE* **136**, No. 5, 575–588.

Daouadji, A., Darve, F., Al Gali, H., Hicher, P. Y., Laouafa, F., Lignon, S., Nicot, F., Nova, R., Pinheiro, M., Prunier, F., Sibille, L. & Wan, R. (2011). Diffuse failure in geomaterials: experiments, theory and modelling. *Int. J. Numer. Anal. Methods Geomech.* **35**, No. 16, 1731–1773.

Darve, F. & Roguiez, X. (1998). Homogeneous bifurcation in soils. *Proc. 4th Int. Workshop on Localization and Bifurcation Theory for Soils and Rocks, Gifu*, pp. 43–50. Rotterdam, the Netherlands: Balkema.

Darve, F., Servant, G., Laouafa, F. & Khoa, H. D. V. (2004). Failure in geomaterials: continuous and discrete analyses. *Comput. Methods Appl. Mech. Engng* **193**, No. 27–29, 3057–3085.

Darve, F., Sibille, L., Daouadji, A. & Nicot, F. (2007). Bifurcations in granular media: macro-and micro-mechanics approaches. *C. R. Mécanique* **335**, No. 9, 496–515.

Goddard, J. D. (2003). Material instability in complex fluids. *Ann. Rev. Fluid Mech.* **35**, 113–133.

Hicher, P. Y. (1998). Experimental behaviour of granular materials. In *Behaviour of granular materials* (ed. B. Cambou), pp. 1–97. Wien, Austria: Springer.

Hicher, P.-Y. & Chang, C. S. (2007). Anisotropic nonlinear elastic model for particulate materials. *J. Geotech. and Geoenviron. Engng ASCE* **132**, No. 8, 1052–1061.

Hill, R. (1958). A general theory of uniqueness and stability in elasto-plastic solids. *J. Mech. Phys. Solids* **6**, 236–249.

Laouafa, F. & Darve, F. (2002). Modelling of slope failure by a material instability mechanism. *Comput. Geotech.* **29**, No. 4, 301–325.

Mindlin, R. D. & Deresiewicz, H. (1953). Elastic spheres in contact under varying oblique forces. *ASME Trans. J. Appl. Mech.* **20**, 327–344.

Mokni, M. & Desrues, J. (1999). Strain localization measurements in undrained plane-strain biaxial tests on Hostun RF sand. *Mech. Cohes.-Frict. Mater.* **4**, No. 4, 419–441.

Nicot, F., Daouadji, A., Laouafa, F. & Darve, F. (2011). Second-order work and kinetic energy in granular media. *Granular Matter* **13**, No. 1, 19–28.

Nova, R. (1994). Controllability of the incremental response of soil specimens subjected to arbitrary loading programs. *J. Mech. Behav. Mater.* **5**, No. 2, 193–201.

Rudnicki, J. W. & Rice, J. (1975). Conditions for the localization of deformation in pressure sensitive dilatant materials. *J. Mech. Phys. Solids* **23**, No. 6, 371–394.

Rothenburg, L. & Selvadurai, A. P. S. (1981). Micromechanical definitions of the Cauchy stress tensor for particular media. In *Mechanics of structured media* (ed. A. P. S. Selvadurai), pp. 469–486. Amsterdam, the Netherlands: Elsevier.

Sivathayalan, S. & Logeswaran, P. (2007). Behaviour of sands under generalized drainage boundary conditions. *Can. Geotech. J.* **44**, No. 2, 138–150.

Schofield, A. N. & Wroth, C. P. (1968). *Critical state soil mechanics*. London, UK: McGraw-Hill.