Bifurcation and generalized mixed loading conditions in geomaterials

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. Copyright © 2010 John Wiley & Sons, Ltd.


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.

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 = 1 − 2 ), whereas the 'volume' is assigned to remain constant ( ε 1 + ε 2 = 0). For this test, the control parameters are: 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: This relation expresses the incremental stress response ( 1 , 2 ) as a function of the incremental strain loading ( ε 1 , ε 2 ); 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: 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: Thus, Equation (1) can be expressed as (see Appendix A) or equivalently 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 = 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 ( C 1 , 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: C 1 = 0 and C 2 = 0. According to Equation (6), a nontrivial incremental response exists if and only if the determinant of the matrix is nil. This condition reads Stress control Strain control 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].
that is, assuming K 22 is not nil: 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 = ( 1 + 2 )/2 continuously decrease and vanish; this is the well-known static liquefaction.

The fundamental role of the determinant of K s
Equation (2) can be inverted as whereas the matrix K is invertible, that is, if det K = 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 s . As geomaterials are known to be nonassociated materials, both matrices K and K s are distinct. In particular, the Bromwich theorem states that the smallest eigenvalue of K s is lower than the real part of the eigenvalues of K (see, for instance, Iordache and Willam [16]). The determinant of K s , 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 s is written as At the deviatoric stress peak, Equation (8) holds and Equation (10) can be rewritten as that is, after rearranging the terms As far as an isochoric triaxial loading path is concerned, the determinant of K s is therefore negative at the loading point where the controllability is lost. In two-dimensional conditions, det K s <0 means that the two eigenvalues have opposite signs. If det K s = 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.

Generalization to proportional strain triaxial paths
The isochoric triaxial loading path is a particular case since we assign the incremental volume, ε 1 + ε 2 , to remain constant. A generalization of this loading path is obtained by prescribing the lateral incremental strain ε 2 to be proportional to the axial incremental strain ε 1 : Proportional strain triaxial paths are controlled by a constant axial strain rate ( ε 1 =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, ε v is plotted versus ε 1 , with ε v = (1−1/ )ε 1 ). The response to this loading path can be investigated by considering the stress term 1 −(1/ ) 2 . For = 1 (isochoric conditions), the deviatoric term q = 1 − 2 is recovered. Figure 3 presents how 1 −(1/ ) 2 evolves in terms of ε 1 , 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 ( ε 1 =constant and ε 1 + ε 2 = 0), the controllability of the test is not lost, even at the ( 1 −(1/ ) 2 ) 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 ( 1 −(1/ ) 2 ) peak:  We assign C 1 as constant and positive, with condition C 2 = 0. According to Equation (3), the response parameters can be chosen as: After rearranging the terms, constitutive relation (3) gives (see Appendix A): If matrix is invertible, Equation (16) is also expressed as: At the 1 −(1/ ) 2 peak (Figure 3), we have C 1 = 0 and C 2 = 0. According to Equation (17), the undefined nontrivial incremental response exists if and only if the determinant of the matrix is nil. This condition reads that is, assuming K 22 is not nil: When condition (19) is met, Equation (16) is no longer valid. A unique incremental response does not exist for the incremental loading ( C 1 , C 2 ). At the ( 1 −(1/ ) 2 ) peak, the controllability of the proportional strain triaxial path is therefore lost when using the control parameters C 1 = 1 −(1/ ) 2 and C 2 = (1/ )ε 1 +ε 2 : this state corresponds to a failure. Interestingly, condition (19) can be explored by noting that the term = K 11 −(1/ )(K 12 + K 21 )+ (1/ 2 )K 22 is a quadratic form with respect to 1/ . Since it can be written as it follows that The quadratic form is therefore associated with the symmetric part K s of matrix K. As a consequence, while K s is definitely positive, the quantity K 11 −(1/ )(K 12 + K 21 )+ (1/ 2 )K 22 is strictly positive for any value of R (the quadratic form is elliptic). When one of the eigenvalues of K s becomes negative, values of exist that make the term K 11 −(1/ )(K 12 + K 21  This is a generalization of the results inferred for the isochoric triaxial path: for the particular value = 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.

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: Matrix K denotes the tangent constitutive operator relating both incremental strain ( ε) and stress ( ) 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 where A and B are two matrices, respectively, of dimension ( p ×n) and ((n − p)×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×2) matrices given by 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: It can be shown (Appendix B) that both control and response parameters can be expressed as where matrix A is decomposed into two block matrices I p, p (identity matrix) and A p,n− p , with the superscripts indicating the dimension of the block matrix: Equations (23a) and (23b) can therefore be rewritten as: Moreover, as (T ε ) −1 = t T , Equation (22) can be equivalently expressed as which, taking Equations (25) and (26) into account, yields: ⎡ where H = T K t T . From Equation (27), it follows that: . . .
and ⎡ 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 R exists, given an incremental loading C.

Controllability analysis
The constitutive matricial equation (36) can be inverted. Indeed, the condition det K>0 implies that det K n− p,n− p = 0. According to Equation (33d), H n− p,n− p = K n− p,n− p , and H n− p,n− p is therefore invertible. Thus, from Equation (35), it follows that ⎡ Combining both Equations (34) and (37) gives which can be expressed as ⎡ For this reason, the loading point considered, at which the controllability of the loading program is lost, is a proper bifurcation point.
Given an integer p, 1 p n, in the following we explore which condition must be fulfilled to ensure the existence of a ( p ×n) matrix A such that det 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: Following Nova's suggestion [12,15], condition (3) can be generalized in an incremental form: As = K ε and C = N R, the second-order work can be regarded as a quadratic form associated with both K s and N s : As a consequence, both N s and K s have the same eigenvalues. However, according to the Bromwich theorem [16,20], stating that the smallest eigenvalue of the symmetric part M s of any square matrix M is lower than any real part of the eigenvalues of M (the inequality is strict when M is nonsymmetric), it follows that the determinant of M s always vanishes before that of M. When the determinant of M is zero, then at least one eigenvalue of M s is strictly negative.
Thus, if the loading program defined by the vector C is not controllable, then det(N) = 0, which implies that at least one eigenvalue of N s (and thus for K s also) is strictly negative.
As a consequence, if there exists one loading program noncontrollable (that is, if an integer p with 1 p n and a ( p ×n) matrix A exist such that det(N) = 0), then K s admits at least one negative eigenvalue.
Conversely, it is assumed that at least one eigenvalue of K s is strictly negative, and at least one eigenvalue is strictly positive. Thus, nonzero vectors x exist so that t xK s x>0, and nonzero vectors y exist so that t yK s y<0. The quadratic form associated with K s , therefore, is neither positive nor negative. The question of whether a matrix A (with p rows and n columns) can be found that ensures det(N) to be nil then arises.
For p = 1, the quantity A 1,n K t A 1,n is a scalar. It follows that det(A 1,n K t A 1,n ) = A 1,n K t A 1,n = A 1,n K s t A 1,n (50) Combining Equations (42), (46) and (50) yields: Since the quadratic form associated with 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 A 1,n is chosen inside the isotropic cone, then A 1,n K s t A 1,n = 0. For this choice of vector t A 1,n , det(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: It is worth noting that when all eigenvalues of K s are strictly positive, all loading programs are controllable. Assuming that prior to any deviatoric loading all eigenvalues of K s are strictly positive, then the existence of a first loading program that is not controllable is observed exactly when det(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(K s ) = 0 corresponds to the boundary of the bifurcation domain.

Illustration of an incrementally nonlinear constitutive relation
3.4.1. The incrementally nonlinear relation. The incremental constitutive equation for rateindependent media relates d and dε by an incrementally nonlinear second-order constitutive relation given by: 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: 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: Finally, using an identification procedure, the following can be obtained: 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].

Instability cones and controllability of mixed loading paths.
An axisymmetric triaxial loading path was simulated after an initial compression at a confining pressure o = 200 kPa. The incrementally nonlinear model was used, calibrated for loose Hostun sand. Then, at a deviatoric stress q = 226 kPa (q = 1 − 3 ), a directional analysis was run: a strain probe ε was imposed in all the directions of the incremental loading space ( ε 1 , ε 2 , ε 3 ) with the same norm (10 −4 ), and the incremental stress response was computed. For each probe direction, the second-order work was assessed: 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 ( 1 = 426 kPa, 2 = 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). For the directions on the cone boundary, the second-order work is nil. These directions therefore belong to the isotropic cone of K s and can be characterized by the ratios [24]: where 2 and 3 satisfy the equation of the isotropic cone f ( 2 , 3 ) = 0. Thus, according to the result inferred in Section 3.3, the vector t A 1,3 = [1, 2 , 3 ] belongs to the isotropic cone of K s . The corresponding loading program defined by the following control parameters, C 2 and C 3 being constant: are no longer controllable at the stress state considered ( 1 = 426 kPa, 2 = 3 = 200 kPa).
To check this, the following proportional strain path is followed after an initial compression at a confining pressure of o = 300 kPa: for the values of 2 and 3 satisfying f ( 2 , 3 ) = 0. Then, to examine the controllability of the loading program defined in Equations (59), the evolution of C 1 = 1 + 2 2 + 3 3 in terms of ε 1 is plotted for certain values of 2 and 3 satisfying f ( 2 , 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 ( 1 = 426 kPa, 2 = 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.

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 i or ε i . A similar notion can be recovered in structural mechanics by introducing additional constraints on a free loaded system: 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 , . . . X n ): for example, a 1 X 1 +···+a n X n = 0. The constrained system can be described by introducing a Lagrangian parameter as follows: More generally, a set of p independent constraints can be imposed: In this case, the constrained system can be described by introducing a set of p Lagrangian parameters i as follows: where A is an (n × p) matrix and = ( 1 , . . ., p ). Equations (65) and (66) can be merged into a single matricial equation, as follows: is a ((n + p)×(n + p)) block matrix. For Equation (67), the existence of another solution that is different from the trivial solution (X = 0 and = 0) requires that the determinant ofK should vanish. Let us assume that K is invertible. Taking advantage of the Schur complement formula (see Appendix D) yields: can be rewritten as: A is an (n × p) matrix. By introducing the ( p × p) matrix H = tÃ KÃ, the vanishing of detK 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.

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 t A 1,n is no longer controllable if the vector t A 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: which, with C 1 = 1 −(1/ ) 2 , C 2 = (1/ )ε 1 +ε 2 , R 1 = ε 1 , and R 2 = 2 , gives: Thus, From Equation (A3): Combining Equations (A4) and (A5) yields: Finally, both Equations (A5) and (A6) can be merged as follows: is invertible, Equation (A7) gives: For = 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 A p, p and A p,n− p (resp. B n− p, p and B n− p,n− p ), with the superscripts indicating the dimension of the block matrix. Thus, adopting an incremental formalism, Equations (22) and (23)  . . .
where both matrices T and T ε below were introduced: Following Nova's suggestion [12,15], condition (3) can be generalized under an incremental form: and then Equation (16) is recovered: