# Bifurcation and generalized mixed loading conditions in geomaterials

François Nicot, Noël Challamel, Jean Lerbet, Florent Prunier, Félix Darve

► To cite this version:

François Nicot, Noël Challamel, Jean Lerbet, Florent Prunier, Félix Darve. Bifurcation and generalized mixed loading conditions in geomaterials. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35 (13), pp.1409-1431. 10.1002/nag.959 . hal-00767565

HAL Id: hal-00767565

https://hal.science/hal-00767565

Submitted on 25 Jun 2018

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Bifurcation and generalized mixed loading conditions
in geomaterials

François Nicot¹,*,†, Noel Challamel², Jean Lerbet³, Florent Prunier⁴ and Félix Darve⁵

¹Cemagref, Unité de Recherche Erosion Torrentielle Neige et Avalanches, Grenoble, France

²INSA, Laboratoire de Génie Mécanique et Génie Civil, Rennes, France

³Université d'Evry Val d'Essonne, UFR Sciences et Technologie, Evry, France

⁴INSA, Lyon, France

⁵Laboratoire Sols Solides Structures, UJF-INPG-CNRS, Grenoble, France

The concept of bifurcation has been widely discussed, based, for example, on homogeneous tests in soil mechanics. Essentially two-dimensional (axisymmetric) conditions have been considered. This paper aims at deriving a bifurcation criterion in geomechanics, independently of the problem's dimension. This criterion is related to the vanishing of the determinant of the symmetric part of the material's (or the discrete system's) constitutive (or stiffness) matrix. The derivation is essentially based on the notion of loading parameters (controlling the loading). Basically, a bifurcation occurs when the existence of a unique solution to the quasistatic problem, involving a given set of loading parameters, is lost. Interestingly, the criterion is shown to be independent of the choice of loading parameters. As a possible extension, the context of structure mechanics is considered, and the close connection between both soil and structure analysis is discussed.

KEY WORDS: bifurcation; constitutive uniqueness; controllability; sustainability; second-order work; loading parameters

# 1. INTRODUCTION

During the last few decades, soil mechanics problems have benefitted from the extensive development of complex incremental constitutive laws, which have been incorporated in most finite element-based codes. These constitutive laws have been calibrated with elementary soil testing, based on the usual control parameters. These control parameters are mainly lateral pressure and axial displacement on the cylinder and are finally reduced to a few basic parameters. From theoretical and experimental points of view, it seems that the specific mode of testing the geomaterial has not been exhaustively explored, compared with how thoroughly efficient constitutive laws have been developed. However, the relevant choice of the control parameters, displacement-based or force-based parameters, is probably a determinant factor of the limit behavior of the soil core. The aim of the present paper is mainly to show that the choice of the control parameters in geotechnical testing is crucial to characterizing a material's strength.

<footnote>*Correspondence to: François Nicot, Cemagref, Unité de Recherche Erosion Torrentielle Neige et Avalanches, Grenoble, France.</footnote>

<footnote>†E-mail: francois.nicot@cemagref.fr</footnote>
Incremental inelastic problems may lose the uniqueness of the evolution problem for large
strains. This loss of uniqueness is generally coupled with a number of specific theoretical and
numerical obstacles. The bifurcation problem can generally be associated with the loss of stability
of the fundamental homogeneous path, eventually leading to the failure state.

In classical views, failure corresponds to a limit state. At this state, it has been shown that certain
loading conditions can no longer be applied. In fact, if an infinitesimal increase in the loading
parameters is prescribed, then the nature of the response abruptly changes from a quasistatic
mode to a dynamical mode. Since there is a discontinuous change in the response of the system
under continuously evolving loading parameters, failure states are also related to a bifurcation
problem. This point has been clearly evidenced in both experimental and numerical investigations.
Indeed, by considering a discrete element method, it was shown that a very small disturbance
applied to a granular assembly in a failure state is sufficient to provoke the sudden collapse of the
specimen [1, 2]. Failure states are also described as unsustainable [3–5].

At a failure state, different response paths can be followed thereafter, depending on the loading
conditions and various imperfections. For a convenient choice of loading parameters and under
the effect of disturbances, the response of the system becomes dynamical. The dynamical regime
is detected by an exponential increase in the system's strain rate, together with the occurrence
of sudden outbursts in kinetic energy. At this point, the static problem no longer has a unique
solution from a mathematical point of view. Hence, failure also implies losing the uniqueness of the
quasistatic problem (see, for instance, Vardoulakis et al. [6, 7], Petrik [8], Bigoni and Hueckel [9],
Bigoni [10], Chambon and Caillerie [11]). However, the reverse is not true. In geomechanics, loss
of uniqueness is also described as a loss of controllability of a certain loading path [12], pointing
out the basic role played by the loading (or control) parameters adopted [13, 14]. This approach
was broadly discussed in two-dimensional (or axisymmetric) contexts. This is, for example, the
relevant context to investigate homogeneous laboratory tests, such as the undrained triaxial tests.

In this paper, the notion of loss of uniqueness (or loss of controllability) is revisited, by gener-
alizing to any dimension of the loading space. For this purpose, the constitutive relation under a
general incremental form is investigated with any *n* degree of freedom (*n*≤6), and the existence of
a unique incremental response, given an incremental loading (defined with specific control param-
eters), is queried. This generalization, based on splitting the loading parameters between kinematic
parameters (involving strain variables) and static parameters (involving stress variables), highlights
the basic role played by the sign of the eigenvalues of the symmetric part of the tangent stiffness
matrix. Moreover, by computing the determinant of the so-called generalized constitutive matrix,
relating the loading variables to the response variables, a relation was inferred that is submitted
to the possible field application of structure mechanics. This last aspect is discussed at the end of
the paper.

Throughout the paper, all stress states considered are strictly within the plastic limit surface, so
that the determinant of the constitutive matrix is strictly positive and has never vanished throughout
the loading history. In addition, only rate-independent materials are considered.

Two-order tensors are represented by bold capital letters (A), whereas vectors are denoted by
X. The summation convention on repeated indices will be employed. For any (one- or two-order)
tensor A, A denotes the transpose tensor.

2. A REVIEW OF THE CONTROLLABILITY APPROACH

2.1. Introduction

In the past decade, Nova has introduced the notion of laboratory test controllability, querying
whether or not a given loading program, at any loading state, could be pursued [12]. The notion
of controllability is intimately related to that of loading (or control) parameters, associated with
the manner in which the loading is applied to the boundary of the specimen considered [13, 14].
To exemplify, the case of the isochoric biaxial test can be commented. In this test, a deviatoric
stress rate is imposed (δq = const, with q = σ₁ − σ₂), whereas the ‘volume’ is assigned to remain
constant ($\delta\epsilon_1 + \delta\epsilon_2 = 0$). For this test, the control parameters are:

$$C_1 = \sigma_1 - \sigma_2 \qquad (1a)$$

$$C_2 = \epsilon_1 + \epsilon_2 \qquad (1b)$$

One control parameter is expressed with stress components, the second one with strain components. This is therefore a mixed control program.

To pursue the notion of controllability further, it is useful to introduce the constitutive relation, in the rate form:

$$\begin{bmatrix} \delta\sigma_1 \\ \delta\sigma_2 \end{bmatrix} = \mathbf{K} \begin{bmatrix} \delta\epsilon_1 \\ \delta\epsilon_2 \end{bmatrix} \qquad (2)$$

This relation expresses the incremental stress response ($\delta\sigma_1, \delta\sigma_2$) as a function of the incremental strain loading ($\delta\epsilon_1, \delta\epsilon_2$); $\mathbf{K}$ is the tangent stiffness matrix. This relation can be rearranged so as to express two proper response parameters ($R_1, R_2$) in terms of the control parameters ($C_1, C_2$). Following Nova's approach, both control and response parameters are related by the following condition, stating that both control and response parameters must be conjugated with respect to energy:

$$C_i R_i = \sigma_i \epsilon_i \qquad (3)$$

where the terms $R_i$ are linear combinations of either stress or strain components.

In the case of the isochoric triaxial loading path, from Equation (1), Equation (3) yields:

$$R_1 = \epsilon_1 \qquad (4a)$$

$$R_2 = \sigma_2 \qquad (4b)$$

Thus, Equation (1) can be expressed as (see Appendix A)

$$\begin{bmatrix} \delta R_1 \\ \delta R_2 \end{bmatrix} = \frac{1}{K_{11} + K_{22} - K_{12} - K_{21}} \begin{bmatrix} 1 & K_{22} - K_{12} \\ K_{21} - K_{22} & K_{11}K_{22} - K_{12}K_{21} \end{bmatrix} \begin{bmatrix} \delta C_1 \\ \delta C_2 \end{bmatrix} \qquad (5)$$

or equivalently

$$\begin{bmatrix} \delta C_1 \\ \delta C_2 \end{bmatrix} = \frac{1}{K_{22}} \begin{bmatrix} K_{11}K_{22} - K_{12}K_{21} & K_{12} - K_{22} \\ K_{22} - K_{21} & 1 \end{bmatrix} \begin{bmatrix} \delta R_1 \\ \delta R_2 \end{bmatrix} \qquad (6)$$

This new constitutive relation (5) between both control and response parameters is defined if and only if the quantity $K_{11} + K_{22} - K_{12} - K_{21}$ does not vanish. While $K_{11} + K_{22} - K_{12} - K_{21} \neq 0$, the test can be pursued. The loading program is applicable, and the test is reputed to be controllable.

On the contrary, if $K_{11} + K_{22} - K_{12} - K_{21} = 0$, Equation (5) is no longer valid. A unique incremental response does not exist for the incremental loading ($\delta C_1, \delta C_2$). The controllability of the test is lost.

For loose specimens, it is shown that the deviatoric stress reaches a peak (Figure 1). At this peak, both incremental deviatoric stress and volume are nil: $\delta C_1 = 0$ and $\delta C_2 = 0$. According to Equation (6), a nontrivial incremental response exists if and only if the determinant of the matrix

$$\mathbf{N} = \frac{1}{K_{22}} \begin{bmatrix} K_{11}K_{22} - K_{12}K_{21} & K_{12} - K_{22} \\ K_{22} - K_{21} & 1 \end{bmatrix}$$

is nil. This condition reads

$$\frac{K_{11} + K_{22} - K_{12} - K_{21}}{K_{22}} = 0 \qquad (7)$$
<figure><img src="image_2.png" /><figcaption>Figure 1. Typical undrained biaxial behavior of a loose sand. At the q-peak, following the axial control parameters, the specimen can liquefy (strain control) or can give rise to a sudden failure (stress control). From Darve et al. [15].</figcaption></figure>

that is, assuming $K_{22}$ is not nil:

$$K_{11} + K_{22} - K_{12} - K_{21} = 0 \qquad (8)$$

Both conditions

$$N \begin{bmatrix} \delta R_1 \\ \delta R_2 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$$

and det $K=0$ imply that nontrivial incremental responses exist but are undefined. As a consequence, at the deviatoric stress peak, the controllability of the test is lost. This state corresponds to a failure for this choice of control parameters. It is worth noting that the choice of the control parameters is crucial, with respect to the controllability of the test. Indeed, if the isochoric test is controlled by the axial strain (a constant axial strain rate is imposed), the deviatoric stress also reaches a peak (for loose specimens), but the test can be pursued beyond this peak. Both deviatoric stress and mean effective pressure $p'=(\sigma_1+\sigma_2)/2$ continuously decrease and vanish; this is the well-known static liquefaction.

## 2.2. The fundamental role of the determinant of $K_s$

Equation (2) can be inverted as

$$\begin{bmatrix} \delta\epsilon_1 \\ \delta\epsilon_2 \end{bmatrix} = K^{-1} \begin{bmatrix} \delta\sigma_1 \\ \delta\sigma_2 \end{bmatrix} \qquad (9)$$

whereas the matrix $K$ is invertible, that is, if det $K \neq 0$. The condition det $K=0$ precisely corresponds to the plastic limit condition. When this condition is not met, a unique incremental strain response exists for any incremental stress loading. In geomaterials, the plastic limit condition is usually given by the Mohr–Coulomb criterion. Experimental results run in the laboratory show that the Mohr–Coulomb criterion is not reached during the isochoric triaxial test. For loose specimens, the stress state at the deviatoric stress peak is far from the Mohr–Coulomb straight line.

The eigenvalues of $K$ are strictly positive at the end of the initial isotropic compression, and it can be assumed that they evolve continuously over any subsequent loading. Moreover, the vanishing of the determinant of $K$ occurs along the Mohr–Coulomb line. At the deviatoric stress
peak, which is far from the Mohr-Coulomb straight line, the determinant of **K** is therefore not nil. It is strictly positive. Thus, a criterion to detect potential loss of controllability should not be based on the matrix **K**. In fact, it is relevant to consider the symmetric part of **K**, denoted as **K**<sub>s</sub>. As geomaterials are known to be nonassociated materials, both matrices **K** and **K**<sub>s</sub> are distinct. In particular, the Bromwich theorem states that the smallest eigenvalue of **K**<sub>s</sub> is lower than the real part of the eigenvalues of **K** (see, for instance, Iordache and Willam [16]). The determinant of **K**<sub>s</sub>, which is strictly positive after a given isotropic compression, may therefore vanish before the determinant of **K** along a given loading path.

The determinant of **K**<sub>s</sub> is written as

$$ \det \mathbf{K}_s = K_{11} K_{22} - \left( \frac{K_{12} + K_{21}}{2} \right)^2 \qquad (10) $$

At the deviatoric stress peak, Equation (8) holds and Equation (10) can be rewritten as

$$ \det \mathbf{K}_s = K_{11} K_{22} - \left( \frac{K_{11} + K_{22}}{2} \right)^2 \qquad (11) $$

that is, after rearranging the terms

$$ \det \mathbf{K}_s = - \left( \frac{K_{11} - K_{22}}{2} \right)^2 \qquad (12) $$

As far as an isochoric triaxial loading path is concerned, the determinant of **K**<sub>s</sub> is therefore negative at the loading point where the controllability is lost. In two-dimensional conditions, det **K**<sub>s</sub> < 0 means that the two eigenvalues have opposite signs. If det **K**<sub>s</sub> = 0, one of the two eigenvalues is nil. In the following section, we attempt to generalize this result by considering more general proportional strain loading paths.

## 2.3. Generalization to proportional strain triaxial paths

The isochoric triaxial loading path is a particular case since we assign the incremental volume, δε<sub>1</sub> + δε<sub>2</sub>, to remain constant. A generalization of this loading path is obtained by prescribing the lateral incremental strain δε<sub>2</sub> to be proportional to the axial incremental strain δε<sub>1</sub>:

$$ \delta\varepsilon_1 + \lambda\delta\varepsilon_2 = 0 \qquad (13) $$

Proportional strain triaxial paths are controlled by a constant axial strain rate (δε<sub>1</sub> = constant), together with condition (13). Different proportional strain paths can be considered, for different values of λ coefficient, as shown in Figure 2 (in this figure, ε<sub>v</sub> is plotted versus ε<sub>1</sub>, with ε<sub>v</sub> = (1 - 1/λ)ε<sub>1</sub>).

The response to this loading path can be investigated by considering the stress term σ<sub>1</sub> - (1/λ)σ<sub>2</sub>. For λ = 1 (isochoric conditions), the deviatoric term q = σ<sub>1</sub> - σ<sub>2</sub> is recovered. Figure 3 presents how σ<sub>1</sub> - (1/λ)σ<sub>2</sub> evolves in terms of ε<sub>1</sub>, by considering a piecewise incrementally linear constitutive relation [15]. Interestingly, it can be observed that the curve passes through a peak for sufficiently small values of λ. As a consequence, a similar reasoning can be applied, as previously done for isochoric triaxial loading paths.

If the test is fully strain-controlled (δε<sub>1</sub> = constant and δε<sub>1</sub> + λδε<sub>2</sub> = 0), the controllability of the test is not lost, even at the (σ<sub>1</sub> - (1/λ)σ<sub>2</sub>) peak. As under isochoric conditions, stresses continuously decrease until they vanish (liquefaction). Thus, the following mixed loading control is proposed to query the controllability of the test at the (σ<sub>1</sub> - (1/λ)σ<sub>2</sub>) peak:

$$ C_1 = \sigma_1 - \frac{1}{\lambda}\sigma_2 \qquad (14a) $$

$$ C_2 = \frac{1}{\lambda}\varepsilon_1 + \varepsilon_2 \qquad (14b) $$
<figure><img src="image_3.png" /><figcaption>Figure 2. Proportional strain triaxial loading paths. From Darve et al. [15].</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Proportional strain triaxial loading paths. From Darve et al. [15].</figcaption></figure>

We assign $C_1$ as constant and positive, with condition $\delta C_2=0$.

According to Equation (3), the response parameters can be chosen as:

$$R_1 = \varepsilon_1 \qquad (15a)$$

$$R_2 = \sigma_2 \qquad (15b)$$

After rearranging the terms, constitutive relation (3) gives (see Appendix A):

$$\begin{bmatrix} \delta R_1 \\ \delta R_2 \end{bmatrix} = \frac{1}{K_{11} + \frac{1}{\lambda^2} K_{22} - \frac{1}{\lambda} (K_{12} + K_{21})} \begin{bmatrix} 1 & \frac{1}{\lambda} (K_{22} - \lambda K_{12}) \\ \frac{1}{\lambda} (\lambda K_{21} - K_{22}) & K_{11} K_{22} - K_{12} K_{21} \end{bmatrix} \begin{bmatrix} \delta C_1 \\ \delta C_2 \end{bmatrix} \qquad (16)$$
If matrix

$$ \begin{bmatrix} 1 & \frac{1}{\lambda}(K_{22} - \lambda K_{12}) \\ \frac{1}{\lambda}(\lambda K_{21} - K_{22}) & K_{11}K_{22} - K_{12}K_{21} \end{bmatrix} $$

is invertible, Equation (16) is also expressed as:

$$ \begin{bmatrix} \delta C_1 \\ \delta C_2 \end{bmatrix} = \frac{1}{K_{22}} \begin{bmatrix} K_{11}K_{22} - K_{12}K_{21} & \frac{1}{\lambda}(\lambda K_{12} - K_{22}) \\ \frac{1}{\lambda}(K_{22} - \lambda K_{21}) & 1 \end{bmatrix} \begin{bmatrix} \delta R_1 \\ \delta R_2 \end{bmatrix} \quad (17) $$

At the $\sigma_1 - (1/\lambda)\sigma_2$ peak (Figure 3), we have $\delta C_1 = 0$ and $\delta C_2 = 0$. According to Equation (17), the undefined nontrivial incremental response exists if and only if the determinant of the matrix

$$ N = \frac{1}{K_{22}} \begin{bmatrix} K_{11}K_{22} - K_{12}K_{21} & \frac{1}{\lambda}(\lambda K_{12} - K_{22}) \\ \frac{1}{\lambda}(K_{22} - \lambda K_{21}) & 1 \end{bmatrix} $$

is nil. This condition reads

$$ \frac{K_{11} - \frac{1}{\lambda}(K_{12} + K_{21}) + \frac{1}{\lambda^2}K_{22}}{K_{22}} = 0 \quad (18) $$

that is, assuming $K_{22}$ is not nil:

$$ K_{11} - \frac{1}{\lambda}(K_{12} + K_{21}) + \frac{1}{\lambda^2}K_{22} = 0 \quad (19) $$

When condition (19) is met, Equation (16) is no longer valid. A unique incremental response does not exist for the incremental loading ($\delta C_1, \delta C_2$). At the ($\sigma_1 - (1/\lambda)\sigma_2$) peak, the controllability of the proportional strain triaxial path is therefore lost when using the control parameters $C_1 = \sigma_1 - (1/\lambda)\sigma_2$ and $C_2 = (1/\lambda)\epsilon_1 + \epsilon_2$: this state corresponds to a failure.

Interestingly, condition (19) can be explored by noting that the term $\chi = K_{11} - (1/\lambda)(K_{12} + K_{21}) + (1/\lambda^2)K_{22}$ is a quadratic form with respect to $1/\lambda$. Since it can be written as

$$ \chi = \begin{bmatrix} 1 & -\frac{1}{\lambda} \end{bmatrix} \begin{bmatrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{bmatrix} \begin{bmatrix} 1 \\ -\frac{1}{\lambda} \end{bmatrix} \quad (20) $$

it follows that

$$ \chi = \begin{bmatrix} 1 & -\frac{1}{\lambda} \end{bmatrix} \mathbf{K} \begin{bmatrix} 1 \\ -\frac{1}{\lambda} \end{bmatrix} = \begin{bmatrix} 1 & -\frac{1}{\lambda} \end{bmatrix} \mathbf{K}_s \begin{bmatrix} 1 \\ -\frac{1}{\lambda} \end{bmatrix} \quad (21) $$

The quadratic form is therefore associated with the symmetric part $\mathbf{K}_s$ of matrix $\mathbf{K}$.

As a consequence, while $\mathbf{K}_s$ is definitely positive, the quantity $K_{11} - (1/\lambda)(K_{12} + K_{21}) + (1/\lambda^2)K_{22}$ is strictly positive for any value of $R$ (the quadratic form is elliptic). When one of the eigenvalues of $\mathbf{K}_s$ becomes negative, values of $\lambda$ exist that make the term $K_{11} - (1/\lambda)(K_{12} + K_{21}) + (1/\lambda^2)K_{22}$ vanish. For such values, vector $[1 \ -1/\lambda]$ belongs to the isotropic cone of the quadratic form.‡

<footnote>‡The isotropic cone of a quadratic form contains all the vectors vanishing the quadratic form.</footnote>
This is a generalization of the results inferred for the isochoric triaxial path: for the particular value $\lambda=1$, condition (8) is recovered.

In two-dimensional conditions, given a loading point, after a given loading history, the existence of one negative eigenvalue (the second one being nil or positive) therefore stands as a necessary and sufficient condition so that a certain loading mode (defined with proper control parameters) exists that cannot be applied. This loading mode is not controllable at that loading point.

This criterion, involving the eigenvalues of $K_s$, is intrinsic. However, this result was derived in a particular two-dimensional context. The purpose of the next section is to generalize this result for any dimension of the loading space.

### 3. AN ATTEMPT AT GENERALIZATION

#### 3.1. Generalized constitutive relation

In geomechanics, for rate-independent materials, stress rates are generally related to strain rates through a constitutive relation that can be written in an incremental form as:

$$ \begin{bmatrix} \delta\sigma_1 \\ \vdots \\ \delta\sigma_n \end{bmatrix} = \mathbf{K} \begin{bmatrix} \delta\epsilon_1 \\ \vdots \\ \delta\epsilon_n \end{bmatrix} \qquad (22) $$

Matrix $K$ denotes the tangent constitutive operator relating both incremental strain ($\delta\epsilon$) and stress ($\delta\sigma$) vectors. $n$ is the dimension of the problem. Equation (22) is related to the material point scale (in the sense of continuum mechanics), and $n$ ranges from 2 to 6. However, for the sake of generality, this section will assume no restriction on $n$. Thereafter, it will be assumed that the mechanical state considered is strictly within the plastic limit surface, in the hardening regime, so that $\det K > 0$.

Given a loading point, after a given loading history, we query whether a certain loading mode (defined with proper control parameters) exists that cannot be applied. This loading mode will be reputed as not controllable at that loading point.

For this purpose, it is useful to generalize the notion of control and response variables, used in the previous section. The loading can be described by a mixed control vector $C$, whose components are composed of linear combinations of both stress and strain components.

Thus, vector $C$ can be expressed as

$$ C_i = A_{ij}\sigma_j \quad \text{for } i=1, \dots, p \quad \text{(static control parameters)} \qquad (23a) $$

$$ C_i = B_{ij}\epsilon_j \quad \text{for } i=p+1, \dots, n \quad \text{(kinematic control parameters)} \qquad (23b) $$

where $A$ and $B$ are two matrices, respectively, of dimension $(p \times n)$ and $((n-p) \times n)$. The $n$ components $C_i$ are referred to as the control parameters.

For two-dimensional proportional strain triaxial loading discussed in the previous section, $A$ and $B$ are both $(1 \times 2)$ matrices given by $A=[1 -1/\lambda]$ and $B=[1 \ \lambda]$.

Many examples of mixed loading involving both static and kinematic variables exist in the field of geomechanics: in an oedometric test, a lateral rigid confinement is prescribed, whereas an axial stress is imposed; in a triaxial test, the lateral confining pressure is assigned, whereas an axial strain rate is imposed. In large-scale geologic topics, such mixed loadings are also relevant: see, for example, Leroy and Molinari [17].

Without altering the generality of the problem, the response vector $R$ can be considered to be composed of stress or strain components:

$$ R_i = \epsilon_i \quad \text{for } i=1, \dots, p \qquad (24a) $$

$$ R_i = \sigma_i \quad \text{for } i=p+1, \dots, n \qquad (24b) $$
It can be shown (Appendix B) that both control and response parameters can be expressed as

$$ \begin{bmatrix} \delta C_1 \\ \vdots \\ \delta C_p \\ \delta R_{p+1} \\ \vdots \\ \delta R_n \end{bmatrix} = \mathbf{T}^{\sigma} \begin{bmatrix} \delta \sigma_1 \\ \vdots \\ \delta \sigma_p \\ \delta \sigma_{p+1} \\ \vdots \\ \delta \sigma_n \end{bmatrix} \qquad (25) $$

$$ \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_p \\ \delta C_{p+1} \\ \vdots \\ \delta C_n \end{bmatrix} = \mathbf{T}^{\varepsilon} \begin{bmatrix} \delta \varepsilon_1 \\ \vdots \\ \delta \varepsilon_p \\ \delta \varepsilon_{p+1} \\ \vdots \\ \delta \varepsilon_n \end{bmatrix} \qquad (26) $$

with

$$ \mathbf{T}^{\sigma} = \begin{bmatrix} \mathbf{I}^{p, p} & \mathbf{A}^{p, n-p} \\ \mathbf{0}^{n-p, p} & \mathbf{I}^{n-p, n-p} \end{bmatrix} \qquad (27) $$

$$ \mathbf{T}^{\varepsilon} = ({}^{t}\mathbf{T}^{\sigma})^{-1} = \begin{bmatrix} \mathbf{I}^{p, p} & \mathbf{0}^{p, n-p} \\ -{}^{t}\mathbf{A}^{p, n-p} & \mathbf{I}^{n-p, n-p} \end{bmatrix} \qquad (28) $$

where matrix A is decomposed into two block matrices $\mathbf{I}^{p, p}$ (identity matrix) and $\mathbf{A}^{p, n-p}$, with the superscripts indicating the dimension of the block matrix: $\mathbf{A}=[\mathbf{I}^{p, p}, \mathbf{A}^{p, n-p}]$. It can be noted that $\det(\mathbf{T}^{\sigma})=\det(\mathbf{I}^{p, p})\det(\mathbf{I}^{n-p, n-p})=1$.

Equations (23a) and (23b) can therefore be rewritten as:

$$ C_i = \sigma_i + \mathbf{A}_{i,1} \sigma_{p+1} + \cdots + \mathbf{A}_{i,n-p} \sigma_n \quad \text{for } i=1, \dots, p \qquad (29a) $$

$$ C_i = -\mathbf{A}_{1,i-p}\varepsilon_1 - \cdots - \mathbf{A}_{p,i-p}\varepsilon_p + \varepsilon_i \quad \text{for } i=p+1, \dots, n \qquad (29b) $$

Moreover, as $(\mathbf{T}^{\varepsilon})^{-1}={}^{t}\mathbf{T}^{\sigma}$, Equation (22) can be equivalently expressed as

$$ \mathbf{T}^{\sigma} \begin{bmatrix} \delta \sigma_1 \\ \vdots \\ \delta \sigma_n \end{bmatrix} = \mathbf{T}^{\sigma} \mathbf{K} {}^{t}\mathbf{T}^{\sigma} \mathbf{T}^{\varepsilon} \begin{bmatrix} \delta \varepsilon_1 \\ \vdots \\ \delta \varepsilon_n \end{bmatrix} \qquad (30) $$
which, taking Equations (25) and (26) into account, yields:

$$ 
\begin{bmatrix}
\delta C_1 \\
\vdots \\
\delta C_p \\
\delta R_{p+1} \\
\vdots \\
\delta R_n
\end{bmatrix}
= \mathbf{H}
\begin{bmatrix}
\delta R_1 \\
\vdots \\
\delta R_p \\
\delta C_{p+1} \\
\vdots \\
\delta C_n
\end{bmatrix}
\qquad (31) $$

where $\mathbf{H} = \mathbf{T}^{\sigma} \mathbf{K}^{t} \mathbf{T}^{\sigma}$.

From Equation (27), it follows that:

$$ \mathbf{H} = \mathbf{T}^{\sigma} \mathbf{K}^{t} \mathbf{T}^{\sigma} = \begin{bmatrix} \mathbf{I}^{p, p} & \mathbf{A}^{p, n-p} \\ \mathbf{0}^{n-p, p} & \mathbf{I}^{n-p, n-p} \end{bmatrix} \mathbf{K} \begin{bmatrix} \mathbf{I}^{p, p} & \mathbf{0}^{n-p, p} \\ \mathbf{t}^{A} & \mathbf{I}^{n-p, n-p} \end{bmatrix} \qquad (32) $$

Expressing matrix $\mathbf{K}$ as a block matrix

$$ \begin{bmatrix} \mathbf{K}^{p, p} & \mathbf{K}^{n-p, p} \\ \mathbf{K}^{p, n-p} & \mathbf{K}^{n-p, n-p} \end{bmatrix} $$

and developing the matrix products in Equation (32) finally yield

$$ \mathbf{H} = \begin{bmatrix} H^{p, p} & H^{n-p, p} \\ H^{p, n-p} & H^{n-p, n-p} \end{bmatrix} $$

with

$$ \mathbf{H}^{p, p} = \mathbf{K}^{p, p} + \mathbf{A}^{p, n-p} \mathbf{K}^{n-p, p} + \mathbf{K}^{p, n-p} \mathbf{t} \mathbf{A}^{p, n-p} + \mathbf{A}^{p, n-p} \mathbf{K}^{n-p, n-p} \mathbf{t} \mathbf{A}^{p, n-p} \quad (33a) $$

$$ \mathbf{H}^{p, n-p} = \mathbf{K}^{p, n-p} + \mathbf{A}^{p, n-p} \mathbf{K}^{n-p, n-p} \qquad (33b) $$

$$ \mathbf{H}^{n-p, p} = \mathbf{K}^{n-p, p} + \mathbf{K}^{n-p, n-p} \mathbf{t} \mathbf{A}^{p, n-p} \qquad (33c) $$

$$ \mathbf{H}^{n-p, n-p} = \mathbf{K}^{n-p, n-p} \qquad (33d) $$

As a consequence, it can be deduced that

$$ \begin{bmatrix} \delta C_1 \\ \vdots \\ \delta C_p \end{bmatrix} = \mathbf{H}^{p, p} \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_p \end{bmatrix} + \mathbf{H}^{p, n-p} \begin{bmatrix} \delta C_{p+1} \\ \vdots \\ \delta C_n \end{bmatrix} \qquad (34) $$

and

$$ \begin{bmatrix} \delta R_{p+1} \\ \vdots \\ \delta R_n \end{bmatrix} = \mathbf{H}^{n-p, p} \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_p \end{bmatrix} + \mathbf{H}^{n-p, n-p} \begin{bmatrix} \delta C_{p+1} \\ \vdots \\ \delta C_n \end{bmatrix} \qquad (35) $$

When matrix $\mathbf{H}^{p, p}$ is invertible, Equations (34) and (35) can be merged into the following matricial relation:

$$ \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_n \end{bmatrix} = \begin{bmatrix} (\mathbf{H}^{p, p})^{-1} & -(\mathbf{H}^{p, p})^{-1} \mathbf{H}^{p, n-p} \\ \mathbf{H}^{n-p, p} (\mathbf{H}^{p, p})^{-1} & \mathbf{H}^{n-p, n-p} - \mathbf{H}^{n-p, p} (\mathbf{H}^{p, p})^{-1} \mathbf{H}^{p, n-p} \end{bmatrix} \begin{bmatrix} \delta C_1 \\ \vdots \\ \delta C_n \end{bmatrix} \qquad (36) $$
where tensor

$$ M = \begin{bmatrix} (\mathbf{H}^{p,p})^{-1} & -(\mathbf{H}^{p,p})^{-1}\mathbf{H}^{p,n-p} \\ \mathbf{H}^{n-p,p}(\mathbf{H}^{p,p})^{-1} & \mathbf{H}^{n-p,n-p} - (\mathbf{H}^{n-p,p}(\mathbf{H}^{p,p})^{-1})\mathbf{H}^{p,n-p} \end{bmatrix} $$

relating both control and response variables, is referred to as the generalized constitutive tensor. Equation (36) is a generalization of Equation (5) or (16) to any loading space dimension $n$. For $n=2$, Equation (16) can be recovered from Equation (36) (see Appendix C).

It is worth noting that no restriction was assigned to infer Equation (36); the validity of this result therefore holds for any rate-independent incrementally piece-wise linear constitutive relation considered and for any dimension of the loading space.

The next section explores the conditions in which the existence of a unique incremental response $\delta \underline{R}$ exists, given an incremental loading $\delta \underline{C}$.

## 3.2. Controllability analysis

The constitutive matricial equation (36) can be inverted. Indeed, the condition $\det \underline{K}>0$ implies that $\det \underline{K}^{n-p,n-p} \neq 0$. According to Equation (33d), $\underline{H}^{n-p,n-p} = \underline{K}^{n-p,n-p}$, and $\underline{H}^{n-p,n-p}$ is therefore invertible. Thus, from Equation (35), it follows that

$$ \begin{bmatrix} \delta C_{p+1} \\ \vdots \\ \delta C_n \end{bmatrix} = -(\underline{H}^{n-p,n-p})^{-1} \underline{H}^{n-p,p} \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_p \end{bmatrix} + (\underline{H}^{n-p,n-p})^{-1} \begin{bmatrix} \delta R_{p+1} \\ \vdots \\ \delta R_n \end{bmatrix} \quad (37) $$

Combining both Equations (34) and (37) gives

$$ \begin{bmatrix} \delta C_1 \\ \vdots \\ \delta C_p \end{bmatrix} = \underline{H}^{p,p} \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_p \end{bmatrix} - \underline{H}^{p,n-p} (\underline{H}^{n-p,n-p})^{-1} \underline{H}^{n-p,p} \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_p \end{bmatrix} + \underline{H}^{p,n-p} (\underline{H}^{n-p,n-p})^{-1} \begin{bmatrix} \delta R_{p+1} \\ \vdots \\ \delta R_n \end{bmatrix} \quad (38) $$

which can be expressed as

$$ \begin{bmatrix} \delta C_1 \\ \vdots \\ \delta C_p \end{bmatrix} = (\underline{H}^{p,p} - \underline{H}^{p,n-p} (\underline{H}^{n-p,n-p})^{-1} \underline{H}^{n-p,p}) \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_p \end{bmatrix} + \underline{H}^{p,n-p} (\underline{H}^{n-p,n-p})^{-1} \begin{bmatrix} \delta R_{p+1} \\ \vdots \\ \delta R_n \end{bmatrix} \quad (39) $$

Finally, by merging Equations (37) and (39), it follows that

$$ \begin{bmatrix} \delta C_1 \\ \vdots \\ \delta C_n \end{bmatrix} = \begin{bmatrix} \underline{H}^{p,p} - \underline{H}^{p,n-p} (\underline{H}^{n-p,n-p})^{-1} \underline{H}^{n-p,p} & \underline{H}^{p,n-p} (\underline{H}^{n-p,n-p})^{-1} \\ -(\underline{H}^{n-p,n-p})^{-1} \underline{H}^{n-p,p} & (\underline{H}^{n-p,n-p})^{-1} \end{bmatrix} \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_n \end{bmatrix} \quad (40) $$
The controllability of the loading program defined with control parameters $C_i$, $i=1, \dots, n$, requires that given $\delta C_i$, $\forall i=1, \dots, n$, a unique solution to Equation (40) exists. This requirement is fulfilled when the determinant of the matrix

$$ \mathbf{N} = \begin{bmatrix} \mathbf{H}^{p,p} - \mathbf{H}^{p,n-p} (\mathbf{H}^{n-p,n-p})^{-1} \mathbf{H}^{n-p,p} & \mathbf{H}^{p,n-p} (\mathbf{H}^{n-p,n-p})^{-1} \\ -(\mathbf{H}^{n-p,n-p})^{-1} \mathbf{H}^{n-p,p} & (\mathbf{H}^{n-p,n-p})^{-1} \end{bmatrix} $$

is not nil. Taking advantage of the Schur complement formula recalled in Appendix D gives

$$ \det \mathbf{N} = \det(\mathbf{H}^{n-p,n-p})^{-1} \det(\mathbf{H}^{p,p} - \mathbf{H}^{p,n-p} (\mathbf{H}^{n-p,n-p})^{-1} \mathbf{H}^{n-p,p} + \dots + \mathbf{H}^{p,n-p} (\mathbf{H}^{n-p,n-p})^{-1} \mathbf{H}^{n-p,n-p} (\mathbf{H}^{n-p,n-p})^{-1} \mathbf{H}^{n-p,p}) \quad (41) $$

that is

$$ \det \mathbf{N} = \frac{\det(\mathbf{H}^{p,p})}{\det(\mathbf{K}^{n-p,n-p})} \quad (42) $$

Thus, the condition $\det \mathbf{N} \neq 0$ finally reads:

$$ \det(\mathbf{H}^{p,p}) \neq 0 \quad (43) $$

When $\det(\mathbf{H}^{p,p})=0$, an infinity of vectors

$$ \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_n \end{bmatrix} $$

satisfy Equation (36) or Equation (40) for a given vector

$$ \begin{bmatrix} \delta C_1 \\ \vdots \\ \delta C_n \end{bmatrix} $$

For this reason, the loading point considered, at which the controllability of the loading program is lost, is a proper bifurcation point.

### 3.3. Bifurcation limit equation

Since $\det(\mathbf{H}^{p,p})$ plays a basic role in the controllability of a given loading program, it is very useful to analyze conditions in which the determinant of $\mathbf{H}^{p,p}$ vanishes.

Equation (33a) is rearranged as follows:

$$ \mathbf{H}^{p,p} = \mathbf{I}^{p,p} \mathbf{K}^{p,p} \mathbf{I}^{p,p} + \mathbf{A}^{p,n-p} \mathbf{K}^{n-p,p} \mathbf{I}^{p,p} + \mathbf{I}^{p,p} \mathbf{K}^{p,n-p} \mathbf{A}^{p,n-p} + \mathbf{A}^{p,n-p} \mathbf{K}^{n-p,n-p} \mathbf{A}^{p,n-p} \quad (44) $$

Equation (44) can therefore be expressed as a product of block matrices, as follows:

$$ \mathbf{H}^{p,p} = [\mathbf{I}^{p,p} \quad \mathbf{A}^{p,n-p}] \begin{bmatrix} \mathbf{K}^{p,p} & \mathbf{K}^{n-p,p} \\ \mathbf{K}^{n-p,p} & \mathbf{K}^{n-p,n-p} \end{bmatrix} \begin{bmatrix} \mathbf{I}^{p,p} \\ \mathbf{A}^{p,n-p} \end{bmatrix} \quad (45) $$

Recalling that $\mathbf{A} = [\mathbf{I}^{p,p} \quad \mathbf{A}^{p,n-p}]$ and

$$ \mathbf{K} = \begin{bmatrix} \mathbf{K}^{p,p} & \mathbf{K}^{n-p,p} \\ \mathbf{K}^{p,n-p} & \mathbf{K}^{n-p,n-p} \end{bmatrix} $$

finally gives:

$$ \mathbf{H}^{p,p} = \mathbf{A} \mathbf{K}^{t} \mathbf{A} \quad (46) $$
Given an integer $p$, $1 \le p \le n$, in the following we explore which condition must be fulfilled to ensure the existence of a ($p \times n$) matrix $\mathbf{A}$ such that $\det \mathbf{H}^{p,p} = 0$.

For this purpose, it is useful to introduce the second-order work $W_2$, defined as the scalar product between both incremental strain and stress [18]. Ignoring geometrical effects [19], $W_2$ is expressed as:

$$W_2 = \delta\sigma_i \delta\epsilon_i \qquad (47)$$

Following Nova's suggestion [12, 15], condition (3) can be generalized in an incremental form:

$$\delta C_i \delta R_i = \delta \sigma_i \delta \epsilon_i \qquad (48)$$

As $\delta\underline{\sigma} = \mathbf{K}\delta\underline{\epsilon}$ and $\delta\underline{C} = \mathbf{N}\delta\underline{R}$, the second-order work can be regarded as a quadratic form associated with both $\mathbf{K}_s$ and $\mathbf{N}_s$:

$$W_2 = {}^t\delta\underline{\epsilon}\mathbf{K}_s\delta\underline{\epsilon} = {}^t\delta\underline{R}\mathbf{N}_s\delta\underline{R} \qquad (49)$$

As a consequence, both $\mathbf{N}_s$ and $\mathbf{K}_s$ have the same eigenvalues.

However, according to the Bromwich theorem [16, 20], stating that the smallest eigenvalue of the symmetric part $\mathbf{M}_s$ of any square matrix $\mathbf{M}$ is lower than any real part of the eigenvalues of $\mathbf{M}$ (the inequality is strict when $\mathbf{M}$ is nonsymmetric), it follows that the determinant of $\mathbf{M}_s$ always vanishes before that of $\mathbf{M}$. When the determinant of $\mathbf{M}$ is zero, then at least one eigenvalue of $\mathbf{M}_s$ is strictly negative.

Thus, if the loading program defined by the vector $\underline{C}$ is not controllable, then $\det(\underline{N})=0$, which implies that at least one eigenvalue of $\mathbf{N}_s$ (and thus for $\mathbf{K}_s$ also) is strictly negative.

As a consequence, if there exists one loading program noncontrollable (that is, if an integer $p$ with $1 \le p \le n$ and a ($p \times n$) matrix $\mathbf{A}$ exist such that $\det(\underline{N})=0$), then $\mathbf{K}_s$ admits at least one negative eigenvalue.

Conversely, it is assumed that at least one eigenvalue of $\mathbf{K}_s$ is strictly negative, and at least one eigenvalue is strictly positive. Thus, nonzero vectors $\underline{x}$ exist so that ${}^t\delta\underline{x}\mathbf{K}_s\delta\underline{x} > 0$, and nonzero vectors $\underline{y}$ exist so that ${}^t\delta\underline{y}\mathbf{K}_s\delta\underline{y} < 0$. The quadratic form associated with $\mathbf{K}_s$, therefore, is neither positive nor negative. The question of whether a matrix $\mathbf{A}$ (with $p$ rows and $n$ columns) can be found that ensures $\det(\underline{N})$ to be nil then arises.

For $p=1$, the quantity $\mathbf{A}^{1,n}\mathbf{K}^t\mathbf{A}^{1,n}$ is a scalar. It follows that

$$\det(\mathbf{A}^{1,n}\mathbf{K}^t\mathbf{A}^{1,n}) = \mathbf{A}^{1,n}\mathbf{K}^t\mathbf{A}^{1,n} = \mathbf{A}^{1,n}\mathbf{K}_s^t\mathbf{A}^{1,n} \qquad (50)$$

Combining Equations (42), (46) and (50) yields:

$$\det(\underline{N}) = \frac{\mathbf{A}^{1,n}\mathbf{K}_s^t\mathbf{A}^{1,n}}{\det\mathbf{K}^{n-1,n-1}} \qquad (51)$$

Since the quadratic form associated with $\mathbf{K}_s$ is neither positive nor negative, the isotropic cone (containing all the vectors vanishing the quadratic form) is not reduced to the nil vector. If the vector ${}^t\mathbf{A}^{1,n}$ is chosen inside the isotropic cone, then ${}^t\mathbf{K}_s^t\mathbf{A}^{1,n} = 0$. For this choice of vector ${}^t\mathbf{A}^{1,n}$, $\det(\underline{N})=0$, and the corresponding loading program is no longer controllable. It can be noted that this loading program is defined by the following set of control parameters:

$$C_1 = \sigma_1 + \mathbf{A}_{1,2}\sigma_2 + \cdots + \mathbf{A}_{1,n}\sigma_n \quad (\mathbf{A}_{1,1} = 1) \qquad (52a)$$

$$C_i = -\mathbf{A}_{1,i}\varepsilon_1 + \varepsilon_i \quad \text{for } i=2, \dots, n \qquad (52b)$$

It is worth noting that when all eigenvalues of $\mathbf{K}_s$ are strictly positive, all loading programs are controllable. Assuming that prior to any deviatoric loading all eigenvalues of $\mathbf{K}_s$ are strictly positive, then the existence of a first loading program that is not controllable is observed exactly when $\det(\mathbf{K}_s)=0$.

In this section, it was specified that the loading point at which the controllability of a loading program is lost is a proper bifurcation point. The set of mechanical states for which a bifurcation can occur is denoted by the bifurcation domain. Thus, the bifurcation limit equation $\det(\mathbf{K}_s)=0$ corresponds to the boundary of the bifurcation domain.
### 3.4. Illustration of an incrementally nonlinear constitutive relation

#### 3.4.1. The incrementally nonlinear relation.
The incremental constitutive equation for rate-independent media relates $d\underline{\sigma}$ and $d\underline{\epsilon}$ by an incrementally nonlinear second-order constitutive relation given by:

$$d\epsilon_{ij} = M_{ijkl}^{1} d\sigma_{kl} + \frac{1}{\|d\underline{\sigma}\|} M_{ijklmn}^{2} d\sigma_{kl} d\sigma_{mn} \quad (53)$$

Assuming in addition that relation (53) is orthotropic, that the shear moduli are incrementally linear, and that there are no ‘crossed’ terms in relation (53), in principal we obtain incremental stress and strain axes:

$$\begin{bmatrix} d\epsilon_1 \\ d\epsilon_2 \\ d\epsilon_3 \end{bmatrix} = A_h \begin{bmatrix} d\sigma_1 \\ d\sigma_2 \\ d\sigma_3 \end{bmatrix} + \frac{1}{\|d\underline{\sigma}\|} B_h \begin{bmatrix} (d\sigma_1)^2 \\ (d\sigma_2)^2 \\ (d\sigma_3)^2 \end{bmatrix} \quad (54)$$

Let us now introduce ‘generalized’ triaxial paths (two constant lateral stresses in fixed principal axes), ‘generalized’ Young moduli, and ‘generalized’ Poisson ratios along these paths: $E_i = (\partial\sigma_i / \partial\epsilon_i)_{\sigma_j, \sigma_h}$ ($\sigma_j, \sigma_k$ constant lateral stresses) and $V_i^j = -(\partial\epsilon_j / \partial\epsilon_i)_{\sigma_j, \sigma_k}$. By distinguishing triaxial compressions ($d\sigma_i > 0$, index ‘+’) from triaxial extensions ($d\sigma_i < 0$, index ‘-’) in axial direction $i$, one can introduce both matrices $H^+$ and $H^-$ defined by:

$$H^+ = \begin{bmatrix} \frac{1}{E_1^+} & -\frac{V_2^{1+}}{E_2^+} & -\frac{V_3^{1+}}{E_3^+} \\ -\frac{V_1^{2+}}{E_1^+} & \frac{1}{E_2^+} & -\frac{V_3^{2+}}{E_3^+} \\ -\frac{V_1^{3+}}{E_1^+} & -\frac{V_2^{3+}}{E_2^+} & \frac{1}{E_3^+} \end{bmatrix} \quad \text{and} \quad H^- = \begin{bmatrix} \frac{1}{E_1^-} & -\frac{V_2^{1-}}{E_2^-} & -\frac{V_3^{1-}}{E_3^-} \\ -\frac{V_1^{2-}}{E_1^-} & \frac{1}{E_2^-} & -\frac{V_3^{2-}}{E_3^-} \\ -\frac{V_1^{3-}}{E_1^-} & -\frac{V_2^{3-}}{E_2^-} & \frac{1}{E_3^-} \end{bmatrix} \quad (55)$$

Finally, using an identification procedure, the following can be obtained:

$$\bar{A} = \frac{1}{2}(\bar{H}^+ + \bar{H}^-) \quad \text{and} \quad \bar{B} = \frac{1}{2}(\bar{H}^+ - \bar{H}^-) \quad (56)$$

Relations (54), (55), and (56) define the incrementally nonlinear model considered throughout this section. Greater detail can be found in Darve and Labanieh [21], Darve [22], and Darve et al. [23].

#### 3.4.2. Instability cones and controllability of mixed loading paths.
An axisymmetric triaxial loading path was simulated after an initial compression at a confining pressure $\sigma_o = 200$ kPa. The incrementally nonlinear model was used, calibrated for loose Hostun sand. Then, at a deviatoric stress $q = 226$ kPa ($q = \sigma_1 - \sigma_3$), a directional analysis was run: a strain probe $\delta\underline{\epsilon}$ was imposed in all the directions of the incremental loading space ($\delta\epsilon_1, \delta\epsilon_2, \delta\epsilon_3$) with the same norm $(10^{-4})$, and the incremental stress response $\delta\underline{\sigma}$ was computed. For each probe direction, the second-order work was assessed:

$$W_2 = \delta\sigma_1\delta\epsilon_1 + \delta\sigma_2\delta\epsilon_2 + \delta\sigma_3\delta\epsilon_3 \quad (57)$$

It is shown that there exists a set of directions along which the second-order work is negative. These directions are contained in a cone, the contour of which is given in Figure 4. Thus, the stress state considered ($\sigma_1 = 426$ kPa, $\sigma_2 = \sigma_3 = 200$ kPa) belongs to the bifurcation domain. In this state, $\det(K_s) < 0$ (one eigenvalue is negative, the two other eigenvalues are positive), where $K_s$ is the symmetric part of the tangent stiffness operator associated with the constitutive relation (54).
<figure><img src="image_5.png" /><figcaption>Figure 4. Intersection of the instability cone with the unit sphere in the incremental strain space.</figcaption></figure>

For the directions on the cone boundary, the second-order work is nil. These directions therefore belong to the isotropic cone of Ks and can be characterized by the ratios [24]:

$$ \lambda_2 = \frac{\delta\varepsilon_2}{\delta\varepsilon_1} \quad \text{and} \quad \lambda_3 = \frac{\delta\varepsilon_3}{\delta\varepsilon_1} \qquad (58) $$

where $\lambda_2$ and $\lambda_3$ satisfy the equation of the isotropic cone $f(\lambda_2, \lambda_3)=0$.

Thus, according to the result inferred in Section 3.3, the vector $^tA^{1,3}=[1, \lambda_2, \lambda_3]$ belongs to the isotropic cone of Ks. The corresponding loading program defined by the following control parameters, $C_2$ and $C_3$ being constant:

$$ C_1 = \sigma_1 + \lambda_2\sigma_2 + \lambda_3\sigma_3 \qquad (59a) $$

$$ C_2 = -\lambda_2\varepsilon_1 + \varepsilon_2 \qquad (59b) $$

$$ C_3 = -\lambda_3\varepsilon_1 + \varepsilon_3 \qquad (59c) $$

are no longer controllable at the stress state considered ($\sigma_1=426$ kPa, $\sigma_2=\sigma_3=200$ kPa).

To check this, the following proportional strain path is followed after an initial compression at a confining pressure of $\sigma_o=300$ kPa:

$$ \delta\varepsilon_1 = \text{const} \qquad (60a) $$

$$ -\lambda_2\delta\varepsilon_1 + \delta\varepsilon_2 = 0 \qquad (60b) $$

$$ -\lambda_3\delta\varepsilon_1 + \delta\varepsilon_3 = 0 \qquad (60c) $$

for the values of $\lambda_2$ and $\lambda_3$ satisfying $f(\lambda_2, \lambda_3)=0$.

Then, to examine the controllability of the loading program defined in Equations (59), the evolution of $C_1=\sigma_1+\lambda_2\sigma_2+\lambda_3\sigma_3$ in terms of $\varepsilon_1$ is plotted for certain values of $\lambda_2$ and $\lambda_3$ satisfying $f(\lambda_2, \lambda_3)=0$. As seen in Figure 5, all the curves pass through a peak. In addition, Figure 6 shows that the peaks are reached more or less at the same stress state, very close to that considered for the directional analysis ($\sigma_1=426$ kPa, $\sigma_2=\sigma_3=200$ kPa).

According to the results obtained in the previous sections, the existence of a unique incremental response to the loading program defined by Equations (59) is lost at the $C_1$ peak.

## 4. APPROACH BY THE LAGRANGE MULTIPLIER

In Section 3, the notion of generalized control and response parameters was introduced for homogeneous problems. Control parameters are built as linear combinations of stress or strain components. The loading is exerted through external constraints defined from linear combinations of $\sigma_i$ or $\varepsilon_i$. A similar notion can be recovered in structural mechanics by introducing additional constraints on a free loaded system:

$$ KX=0 \qquad (61) $$
<figure><img src="image_6.png" /><figcaption>Figure 5. Existence of $C_1$ peaks for different values of $\lambda_2$ and $\lambda_3$ such that vectors $^tA^{1,3}=[1, \lambda_2, \lambda_3]$ belong to the isotropic cone of $K_s$.</figcaption></figure>

where **K** is the linearized generalized stiffness matrix, and **X** is the perturbation around the equilibrium solution of the structural problem.

A tangent constraint can be defined as a linear equation between kinematic variables ($X_1, \dots, X_n$): for example, $a_1X_1+\dots+a_nX_n=0$. The constrained system can be described by introducing a Lagrangian parameter $\lambda$ as follows:

$$ K\underline{X} + \lambda\underline{A} = \underline{0} \qquad (62) $$

$$ ^tA\underline{X} = \underline{0} \qquad (63) $$

More generally, a set of $p$ independent constraints can be imposed:

$$ A_{1i}X_1 + \dots + A_{ni}X_n = \underline{0}, \quad i=1, \dots, p \qquad (64) $$

In this case, the constrained system can be described by introducing a set of $p$ Lagrangian parameters $\lambda_i$ as follows:

$$ K\underline{X} + A\underline{\lambda} = \underline{0} \qquad (65) $$

$$ ^tA\underline{X} = \underline{0} \qquad (66) $$

where **A** is an $(n \times p)$ matrix and $\underline{\lambda}=(\lambda_1, \dots, \lambda_p)$.

Equations (65) and (66) can be merged into a single matricial equation, as follows:

$$ \begin{bmatrix} \mathbf{K} & \mathbf{A} \\ ^t\mathbf{A} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \underline{\mathbf{X}} \\ \underline{\lambda} \end{bmatrix} = \underline{0} \qquad (67) $$

where

$$ \tilde{\mathbf{K}} = \begin{bmatrix} \mathbf{K} & \mathbf{A} \\ ^t\mathbf{A} & \mathbf{0} \end{bmatrix} $$

is a $((n+p) \times (n+p))$ block matrix.
<figure><img src="image_7.png" /><figcaption>Figure 6. $\sigma_1$, $\sigma_2$, and $\sigma_3$ versus $C_1=\sigma_1+\lambda_2\sigma_2+\lambda_3\sigma_3$. At the maximum value of $C_1$, $\sigma_2$ and $\sigma_3$ are close to 200 kPa (between 180 and 220 kPa) and $\sigma_1$ is close to 426 kPa (between 420 and 440 kPa).</figcaption></figure>

For Equation (67), the existence of another solution that is different from the trivial solution ($\underline{X}=\underline{0}$ and $\underline{\lambda}=\underline{0}$) requires that the determinant of $\tilde{K}$ should vanish. Let us assume that $K$ is invertible. Taking advantage of the Schur complement formula (see Appendix D) yields:

$$ \det \tilde{K} = \det K \det(\underline{0} - \underline{1}A(K)^{-1}A) \quad (68) $$

Noting that $\underline{1}A(K)^{-1}A = \underline{1}A(\underline{1}K)^{-1}\underline{1}K(K)^{-1}A = \underline{1}\tilde{A}\underline{1}K\tilde{A} = \underline{1}(\tilde{A}K\tilde{A})^{-1}A$, Equation (68) can be rewritten as:

$$ \frac{\det \tilde{K}}{\det K} = (-1)^p \det(\tilde{A}K\tilde{A}) \quad (69) $$

$\tilde{A}$ is an $(n \times p)$ matrix. By introducing the $(p \times p)$ matrix $H=\tilde{A}K\tilde{A}$, the vanishing of $\det \tilde{K}$ is related to the vanishing of $\det H$.

Thus, we formally arrive exactly at the same problem as that considered in Section 3.2. Although the physical meaning of the matrix is different (strictly speaking, $K$ is not a constitutive matrix and depends on the equilibrium configuration and forces required for this configuration to exist), the same conclusion thereby holds true: the vanishing of $\det H$ implies that the symmetric part $K_s$ admits two eigenvalues of opposite signs.
In conclusion, as far as structure mechanics is concerned, det($K_s$)≤0 is a necessary condition for the existence of loadings (parameterized by A) associated with several equilibrium configurations (bifurcation condition).

As in soil mechanics, where the existence of a regular matricial relation between both control and response parameters was investigated, the vanishing of det $K_s$ plays a fundamental role in structural mechanics to detect a bifurcation point. As such, the notion of control parameters (characterizing the loading) is replaced with the notion of constraints (applied to the kinematic variables of the problem). This question was recently considered in a series of papers [25, 26]. The use of the second-order work criterion to detect structure collapse will be investigated in a forthcoming paper.

## 5. CONCLUDING REMARKS

This paper has investigated the question of bifurcation in soil mechanics, with a possible extension in structural mechanics. The singularity of the symmetric part of the tangent stiffness matrix was shown to play a fundamental role. In soil mechanics, assuming a material point problem, the tangent stiffness matrix corresponds to the constitutive relation, whereas in structural mechanics the stiffness matrix encompasses both the constitutive behavior of the materials and the balance equations.

In soil mechanics, the loading is usually applied by means of specific control variables, and the response is characterized by response parameters. In this paper, we have investigated the conditions in which, at a given mechanical state, the existence and uniqueness of the incremental response of the system are lost given any incremental loading: for a problem of dimension $n$, the loading is composed of $p$ static variables (linear combinations of the stress components) and of $n-p$ kinematic variables (linear combinations of the strain components). It was shown that the determinant of the symmetric part of the tangent stiffness matrix must vanish: $K_s$ admits two eigenvalues of opposite signs. In addition, for the case $p=1$ and when $K_s$ admits one negative eigenvalue and one positive eigenvalue, the loading program associated with the vector $^tA^{1,n}$ is no longer controllable if the vector $^tA^{1,n}$ is chosen inside the isotropic cone. This result was checked in three-dimensional conditions by considering an incrementally nonlinear constitutive relation.

Interestingly, this approach was extended to the case of structural mechanics by introducing the notion of kinematic constraints. The existence of kinematic constraints prescribed to a structure can be regarded in structural mechanics as a notion similar to that of the control (or loading) parameter. Generalizing the results inferred in two-dimensional conditions [25], an $n$-dimension problem was considered, with the adjunction of $p$ kinematic constraints. The uniqueness of the static problem's solution was shown to be lost when the determinant of the so-called stiffness matrix $K_s$ has vanished.

# APPENDIX A

Equation (2) can be transformed to make both control and response variables appear:

$$ \begin{bmatrix} 1 & -\frac{1}{\lambda} \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \delta\sigma_1 \\ \delta\sigma_2 \end{bmatrix} = \begin{bmatrix} 1 & -\frac{1}{\lambda} \\ 0 & \lambda \end{bmatrix} \begin{bmatrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{bmatrix} \begin{bmatrix} 1 & 0 \\ \frac{1}{\lambda} & 1 \end{bmatrix}^{-1} \begin{bmatrix} 1 & 0 \\ \frac{1}{\lambda} & 1 \end{bmatrix} \begin{bmatrix} \delta\epsilon_1 \\ \delta\epsilon_2 \end{bmatrix} \quad (A1) $$

which, with $C_1 = \sigma_1 - (1/\lambda)\sigma_2$, $C_2 = (1/\lambda)\epsilon_1 + \epsilon_2$, $R_1 = \epsilon_1$, and $R_2 = \sigma_2$, gives:

$$ \begin{bmatrix} \delta C_1 \\ \delta R_2 \end{bmatrix} = \begin{bmatrix} K_{11} - \frac{1}{\lambda}(K_{12} + K_{21}) + \frac{1}{\lambda^2}K_{22} & K_{12} - \frac{1}{\lambda}K_{22} \\ K_{21} - \frac{1}{\lambda}K_{22} & K_{22} \end{bmatrix} \begin{bmatrix} \delta R_1 \\ \delta C_2 \end{bmatrix} \quad (A2) $$
Thus,

$$ \delta C_1 = \left( K_{11} - \frac{1}{\lambda}(K_{12} + K_{21}) + \frac{1}{\lambda^2}K_{22} \right) \delta R_1 + \left( K_{12} - \frac{1}{\lambda}K_{22} \right) \delta C_2 \quad (A3) $$

$$ \delta R_2 = \left( K_{21} - \frac{1}{\lambda} K_{22} \right) \delta R_1 + K_{22} \delta C_2 \quad (A4) $$

From Equation (A3):

$$ \delta R_1 = \frac{1}{K_{11} - \frac{1}{\lambda}(K_{12} + K_{21}) + \frac{1}{\lambda^2}K_{22}} \left( \delta C_1 - \left( K_{12} - \frac{1}{\lambda}K_{22} \right) \delta C_2 \right) \quad (A5) $$

Combining Equations (A4) and (A5) yields:

$$ \delta R_2 = \frac{1}{K_{11} - \frac{1}{\lambda}(K_{12} + K_{21}) + \frac{1}{\lambda^2}K_{22}} \left( \left( K_{21} - \frac{1}{\lambda}K_{22} \right) \delta C_1 + (K_{11}K_{22} - K_{12}K_{21}) \delta C_2 \right) \quad (A6) $$

Finally, both Equations (A5) and (A6) can be merged as follows:

$$ \begin{bmatrix} \delta R_1 \\ \delta R_2 \end{bmatrix} = \frac{1}{K_{11} - \frac{1}{\lambda}(K_{12} + K_{21}) + \frac{1}{\lambda^2}K_{22}} \begin{bmatrix} 1 & \frac{1}{\lambda}K_{22} - K_{12} \\ K_{21} - \frac{1}{\lambda}K_{22} & K_{11}K_{22} - K_{12}K_{21} \end{bmatrix} \begin{bmatrix} \delta C_1 \\ \delta C_2 \end{bmatrix} \quad (A7) $$

If matrix

$$ \begin{bmatrix} 1 & \frac{1}{\lambda}K_{22} - K_{12} \\ K_{21} - \frac{1}{\lambda}K_{22} & K_{11}K_{22} - K_{12}K_{21} \end{bmatrix} $$

is invertible, Equation (A7) gives:

$$ \begin{bmatrix} \delta C_1 \\ \delta C_2 \end{bmatrix} = \frac{1}{K_{22}} \begin{bmatrix} K_{11}K_{22} - K_{12}K_{21} & \frac{1}{\lambda}(\lambda K_{12} - K_{22}) \\ \frac{1}{\lambda}(K_{22} - \lambda K_{21}) & 1 \end{bmatrix} \begin{bmatrix} \delta R_1 \\ \delta R_2 \end{bmatrix} \quad (A8) $$

For $\lambda=1$, Equations (5) and (6) are readily obtained from Equations (A7) and (A8).

## APPENDIX B

Starting from Equations (22) and (23), matrix **A** (resp. **B**) can be decomposed into two block matrices $\mathbf{A}^{p,p}$ and $\mathbf{A}^{p,n-p}$ (resp. $\mathbf{B}^{n-p,p}$ and $\mathbf{B}^{n-p,n-p}$), with the superscripts indicating the dimension of the block matrix. Thus, adopting an incremental formalism, Equations (22) and (23) can be rewritten as:

$$ \begin{bmatrix} \delta C_1 \\ \vdots \\ \delta C_p \\ \delta R_{p+1} \\ \vdots \\ \delta R_n \end{bmatrix} = \mathbf{T}^{\sigma} \begin{bmatrix} \delta \sigma_1 \\ \vdots \\ \delta \sigma_p \\ \delta \sigma_{p+1} \\ \vdots \\ \delta \sigma_n \end{bmatrix} \quad (B1) $$
$$ \begin{bmatrix} \delta R_1 \\ \vdots \\ \delta R_p \\ \delta C_{p+1} \\ \vdots \\ \delta C_n \end{bmatrix} = \mathbf{T}^{\mathbf{\epsilon}} \begin{bmatrix} \delta \epsilon_1 \\ \vdots \\ \delta \epsilon_p \\ \delta \epsilon_{p+1} \\ \vdots \\ \delta \epsilon_n \end{bmatrix} \qquad (B2) $$

where both matrices $\mathbf{T}^{\sigma}$ and $\mathbf{T}^{\epsilon}$ below were introduced:

$$ \mathbf{T}^{\sigma} = \begin{bmatrix} \mathbf{A}^{p,p} & \mathbf{A}^{p,n-p} \\ \mathbf{0}^{n-p,p} & \mathbf{I}^{n-p,n-p} \end{bmatrix} \qquad (B3) $$

$$ \mathbf{T}^{\epsilon} = \begin{bmatrix} \mathbf{I}^{p,p} & \mathbf{0}^{p,n-p} \\ \mathbf{B}^{n-p,p} & \mathbf{B}^{n-p,n-p} \end{bmatrix} \qquad (B4) $$

Following Nova's suggestion [12, 15], condition (3) can be generalized under an incremental form:

$$ \delta C_i \delta R_i = \delta \sigma_i \delta \epsilon_i \qquad (B5) $$

By combining Equations (B1) and (B2), it follows that

$$ \delta C_i \delta R_i = [\epsilon_1 \ \dots \ \delta \epsilon_p \ \delta \epsilon_{p+1} \ \dots \ \delta \epsilon_n]^{\mathrm{t}\mathbf{T}^{\epsilon}} \mathbf{T}^{\sigma} \begin{bmatrix} \delta \sigma_1 \\ \vdots \\ \delta \sigma_p \\ \delta \sigma_{p+1} \\ \vdots \\ \delta \sigma_n \end{bmatrix} \qquad (B6) $$

and according to Equation (B5),

$$ \delta C_i \delta R_i = [\delta \epsilon_1 \ \dots \ \delta \epsilon_p \ \delta \epsilon_{p+1} \ \dots \ \delta \epsilon_n] \begin{bmatrix} \delta \sigma_1 \\ \vdots \\ \delta \sigma_p \\ \delta \sigma_{p+1} \\ \vdots \\ \delta \sigma_n \end{bmatrix} \qquad (B7) $$

for any $\delta\sigma_i$ and $\delta\epsilon_i$, for $i=1, \dots, n$.
Equations (B6) and (B7) readily imply that

$$ \mathbf{T}^{\sigma \mathrm{t} \mathbf{T}^{\epsilon}} = \mathbf{I} \quad \text{or equivalently} \quad (\mathbf{T}^{\epsilon})^{-1} = \mathrm{t}^{\mathrm{T}^{\sigma}} \qquad (B8) $$
Moreover, from Equations (B3) and (B4), it can be written that

$$
\begin{aligned}
\mathbf{T}^{\sigma} \mathbf{T}^{\epsilon} &= \begin{bmatrix} \mathbf{A}^{p,p} & \mathbf{A}^{p,n-p} \\ \mathbf{0}^{n-p,p} & \mathbf{I}^{n-p,n-p} \end{bmatrix} \begin{bmatrix} \mathbf{I}^{p,p} & \mathbf{t}\mathbf{B}^{n-p,p} \\ \mathbf{0}^{p,n-p} & \mathbf{t}\mathbf{B}^{n-p,n-p} \end{bmatrix} \\
&= \begin{bmatrix} \mathbf{A}^{p,p} & \mathbf{A}^{p,p} \mathbf{t}\mathbf{B}^{n-p,p} + \mathbf{A}^{p,n-p} \mathbf{t}\mathbf{B}^{n-p,n-p} \\ \mathbf{0}^{n-p,p} & \mathbf{t}\mathbf{B}^{n-p,n-p} \end{bmatrix}
\end{aligned}
$$

Taking advantage of Equation (B8), $T^{\sigma}T^{\epsilon} = I$, and it can finally be derived that

$$ \begin{bmatrix} \mathbf{A}^{p,p} & \mathbf{A}^{p,p} \mathbf{t} \mathbf{B}^{n-p,p} + \mathbf{A}^{p,n-p} \mathbf{t} \mathbf{B}^{n-p,n-p} \\ \mathbf{0}^{n-p,p} & \mathbf{t} \mathbf{B}^{n-p,n-p} \end{bmatrix} = \begin{bmatrix} \mathbf{I}^{p,p} & \mathbf{0}^{p,n-p} \\ \mathbf{0}^{n-p,p} & \mathbf{I}^{n-p,n-p} \end{bmatrix} \quad (B9) $$

that is:

$$ \mathbf{A}^{p,p} = \mathbf{I}^{p,p} \quad (B10) $$

$$ \mathbf{B}^{n-p,n-p} = \mathbf{I}^{n-p,n-p} \quad (B11) $$

$$ \mathbf{A}^{p,n-p} = -\mathbf{t}\mathbf{B}^{n-p,p} \quad (B12) $$

Finally, both matrices $T^{\sigma}$ and $T^{\epsilon}$ are expressed as:

$$ T^{\sigma} = \begin{bmatrix} \mathbf{I}^{p,p} & \mathbf{A}^{p,n-p} \\ \mathbf{0}^{n-p,p} & \mathbf{I}^{n-p,n-p} \end{bmatrix} \quad (B13) $$

$$ T^{\epsilon} = (\mathbf{t}T^{\sigma})^{-1} = \begin{bmatrix} \mathbf{I}^{p,p} & \mathbf{0}^{p,n-p} \\ -\mathbf{t}\mathbf{A}^{p,n-p} & \mathbf{I}^{n-p,n-p} \end{bmatrix} \quad (B14) $$

## APPENDIX C

For $n=2$ and $p=1$, $\mathbf{A}^{p,n-p}$ is a scalar. Let $\mathbf{A}^{p,n-p} = -(1/\lambda)$. Equation (36) reads

$$ \begin{bmatrix} \delta R_1 \\ \delta R_2 \end{bmatrix} = \begin{bmatrix} (\mathbf{H}^{p,p})^{-1} & -(\mathbf{H}^{p,p})^{-1}\mathbf{H}^{p,n-p} \\ \mathbf{H}^{n-p,p}(\mathbf{H}^{p,p})^{-1} & \mathbf{H}^{n-p,n-p} - \mathbf{H}^{n-p,p}(\mathbf{H}^{p,p})^{-1}\mathbf{H}^{p,n-p} \end{bmatrix} \begin{bmatrix} \delta C_1 \\ \delta C_2 \end{bmatrix} \quad (C1) $$

with

$$ \mathbf{H}^{p,p} = K_{11} - \frac{1}{\lambda} K_{21} - \frac{1}{\lambda} K_{12} + K_{22} \frac{1}{\lambda^2} \quad (C2a) $$

$$ \mathbf{H}^{p,n-p} = K_{12} - \frac{1}{\lambda} K_{22} \quad (C2b) $$

$$ \mathbf{H}^{n-p,p} = K_{21} - \frac{1}{\lambda} K_{22} \quad (C2c) $$

$$ \mathbf{H}^{n-p,n-p} = K_{22} \quad (C2d) $$

Thus,

$$ \begin{aligned} \mathbf{H}^{n-p,n-p} - \mathbf{H}^{n-p,p} (\mathbf{H}^{p,p})^{-1} \mathbf{H}^{p,n-p} &= (\mathbf{H}^{p,p})^{-1} (\mathbf{H}^{n-p,n-p} \mathbf{H}^{p,p} - \mathbf{H}^{n-p,p} \mathbf{H}^{p,n-p}) \\ &= (\mathbf{H}^{p,p})^{-1} \det(\mathbf{H}) \quad (C3) \end{aligned} $$
As $\mathbf{H}=\mathbf{T}^{\sigma}\mathbf{K}^{\prime}\mathbf{T}^{\sigma}$, with $\det(\mathbf{T}^{\sigma})=1$, it follows that $\det\mathbf{H}=\det\mathbf{K}=K_{11}K_{22}-K_{12}K_{21}$.
Thus,

$$ \begin{bmatrix} \delta R_1 \\ \delta R_2 \end{bmatrix} = \begin{bmatrix} (\mathbf{H}^{p,p})^{-1} & -(\mathbf{H}^{p,p})^{-1}\mathbf{H}^{p,n-p} \\ (\mathbf{H}^{p,p})^{-1}\mathbf{H}^{n-p,p} & (\mathbf{H}^{p,p})^{-1}\det\mathbf{H} \end{bmatrix} \begin{bmatrix} \delta C_1 \\ \delta C_2 \end{bmatrix} \quad (C4) $$

and then Equation (16) is recovered:

$$ \begin{bmatrix} \delta R_1 \\ \delta R_2 \end{bmatrix} = \frac{1}{K_{11} - \frac{1}{\lambda}(K_{12} + K_{21}) + \frac{1}{\lambda^2}K_{22}} \begin{bmatrix} 1 & \frac{1}{\lambda}(K_{22} - \lambda K_{12}) \\ \frac{1}{\lambda}(\lambda K_{21} - K_{22}) & K_{11}K_{22} - K_{12}K_{21} \end{bmatrix} \begin{bmatrix} \delta C_1 \\ \delta C_2 \end{bmatrix} \quad (C5) $$

# APPENDIX D

Schur complement formula: For any block matrix

$$ \mathbf{N} = \begin{bmatrix} \mathbf{A} & \mathbf{C} \\ \mathbf{B} & \mathbf{D} \end{bmatrix} $$

if **A**, **B**, **C**, and **D** are, respectively, $(p \times p)$, $(p \times q)$, $(q \times p)$, and $(q \times q)$ matrices, and **A** is invertible, then $\det \mathbf{N} = \det \mathbf{D} \det(\mathbf{A} - \mathbf{B}\mathbf{D}^{-1}\mathbf{C})$.

# ACKNOWLEDGEMENTS

The French Federative Research Structure VOR (Natural Hazards and Structure Vulnerability) is gratefully acknowledged by the authors. The European project LESSLOSS, and the French national (ANR) projects SIGMA, STABROCK, and Snow-White also supported this research.

# REFERENCES

1. Sibille L, Nicot F, Donze F, Darve F. Material instability in granular assemblies from fundamentally different models. *International Journal for Numerical and Analytical Methods in Geomechanics* 2007; **31**:457–481.

2. Sibille L, Donzé F, Nicot F, Chareyre B, Darve F. Bifurcation detection and catastrophic failure. *Acta Geotecnica* 2007; **3**(1):14–24.

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

4. Nicot F, Darve F, Khoa HDV. Bifurcation and second-order work in geomaterials. *International Journal for Numerical and Analytical Methods in Geomechanics* 2007; **31**:1007–1032.

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

6. Vardoulakis I, Goldscheider M, Gudehus G. Formation of shear bands in sand bodies as a bifurcation problem. *International Journal for Numerical and Analytical Methods in Geomechanics* 1978; **2**(2):99–128.

7. Vardoulakis I, Sulem J. Bifurcation Analysis in Geomechanics. Chapman and Hall: London, 1995.

8. Petryk H. Theory of bifurcation and instability in time-independent plasticity. In *Bifurcation and Stability of Dissipative Systems*, Nguyen QS (ed.). CISM Courses and Lecturers, vol. 327. Springer: Berlin, 1993; 95–152.

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

10. Bigoni D. Bifurcation and instability of non-associative elastoplastic solids. In *Material Instabilities in Elastic and Plastic Solids*, Petryk H (ed.). CISM Courses and Lecturers, vol. 414. Springer: Berlin, 2000; 1–52.

11. Chambon R, Caillerie D. Existence and uniqueness theorems for boundary value problems involving incrementally nonlinear models. *International Journal of Solids and Structures* 1999; **36**:5089–5099.

12. Nova R. Controllability of the incremental response of soil specimens subjected to arbitrary loading programs. *Journal of Mechanical behavior of Materials* 1994; **5**(2):193–201.

13. Runesson K, Mroz Z. A note on nonassociated plastic flow rules. *International Journal of Plasticity* 1989; **5**(6):639–658.

14. Klisinski M, Mroz Z, Runesson K. Structure of constitutive equations in plasticity for different choices of state and control variables. *International Journal of Plasticity* 1992; **8**(3):221–243.
15. Darve F, Servant G, Laouafa F, Khoa HDV. Failure in geomaterials, continuous and discrete analyses. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**:3057–3085.

16. Iordache MM, Willam KJ. Localized failure analysis in elastoplastic Cosserat continua. *Computer Methods in Applied Mechanics and Engineering* 1998; **151**(3-4):559–586.

17. Leroy Y, Molinari A. Stability of steady states in shear zones. *Journal of the Mechanics and Physics of Solids* 1992; **40**:181–212.

18. Hill R. A general theory of uniqueness and stability in elastic–plastic solids. *Journal of the Mechanics and Physics of Solids* 1958; **6**:236–249.

19. Petryk H. Theory of material instability in incrementally nonlinear plasticity. In *CISM Material Instabilities in Elastic and Plastic Solids*, Petrik H (ed.), vol. 414. Springer: Berlin, 1999; 261–331.

20. Ishaq M. Sur les spectres des matrices. *Séminaire Dubreil, Algèbre et Théorie des Nombres* 1955; **9**:1–14.

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

22. Darve F. The expression of rheological laws in incremental form and the main classes of constitutive equations. In *Geomaterials Constitutive Equations and Modelling*, Darve F (ed.). Taylor & Francis: London, 1990; 123–148.

23. Darve F, Flavigny E, Meghachou M. Constitutive modeling and instabilities of soil behaviour. *Computers and Geotechnics* 1995; **17**:203–224.

24. Prunier F, Nicot F, Darve F, Laouafa F, Lignon S. 3D multiscale bifurcation analysis of granular media. *Journal of Engineering Mechanics (ASCE)* 2009; **135**(6):1–17.

25. Challamel N, Nicot F, Lerbet J, Darve F. On the stability of non-conservative elastic systems under mixed perturbations. *European Journal of Environmental and Civil Engineering* 2008; **13**(3):347–367.

26. Challamel N, Nicot F, Lerbet J, Darve F. A theoretical link between second-order work and stability under kinematics constraints. *Engineering Structures* 2010; in press.