3D bifurcation analysis in geomaterials

ABSTRACT In this paper, a study of failure in geomaterials is proposed using the second order work criterion and a phenomenological approach. In a first part, an analytical investigation of this criterion is proposed. General 3D equation of instability cones as well as of the 3D bifurcation domain limit are given for every incrementally piece-wise linear constitutive model. In the second part, these instability cones and bifurcation domain are displayed for constitutive models of Darve. Then a physical interpretation of these results is proposed. Noticeably, it is proved that, when the second order work vanishes along a given loading path, an extension a generalized flow rule can be defined for proper conjugated variables. Finally we conclude that this approach constitutes a good starting point when investigating bifurcation in geomaterials, and opens new horizons in experimental testing.


3D bifurcation analysis in geomaterials
farid.laouafa@ineris.frABSTRACT.In this paper, a study of failure in geomaterials is proposed through the analysis of the second order work criterion.The analysis is restricted to the material point scale.Then a phenomenological approach based on "macroscopic" constitutive models is adopted.In a first part, an analytical investigation of this criterion is proposed.General 3D equation of instability cones as well as of the 3D bifurcation domain limit are given for every incrementally piece-wise linear constitutive model.In the second part, these instability cones and bifurcation domain are displayed for constitutive models of Darve.Then a physical interpretation of these results is proposed.Noticeably, it is proved that, when the second order work vanishes along a given loading path, an extension of the notions of failure limit condition as well as of flow rule can be defined for proper conjugated variables.It is also shown, that the Lode's angle (and not only the mean pressure associated with a deviatoric level) is also important according to bifurcation.Finally we conclude that this approach constitutes a good starting point when investigating bifurcation in geomaterials, and opens new horizons in experimental testing.
Methods which allow to analyse such kind of failures, are generally based on loss of uniqueness of particular constitutive equation solutions.In that sense, criteria issue from these analyses are particular type of bifurcation criteria by opposition to stability criteria.Motions equations are not investigated, and as a consequence, stability of an equilibrium state (or even of a trajectory of a solid) is not studied.Nevertheless, loss of uniqueness of such equation solutions characterize the occurrence of a discontinuity in the response mode by keeping continuous loading parameters.This is the basic feature of a bifurcation.
Limit of plasticity and localization criterion of Rice (Rice, 1976), (Rudnicki et al., 1975), are the most common bifurcation criteria used by engineers.Nevertheless, as already written above, strain localization can appear before the plasticity limit condition, but in addition, diffuse failure modes which are neither described by the plasticity condition, nor by the Rice criterion can occur.The Hill's sufficient condition of stability seems describing accurately these different failure modes.Furthermore, Bigoni and Hueckel (Bigoni et al., 1991) have proved that this condition is reached before Rice condition and limit plasticity one.
This condition, stated by Hill (Hill, 1958), reads the following expression: with V the volume of the body at stage t, s i j the transpose of the first Piola-Kirchoff tensor, and d ∂u j ∂x i the deformation velocity.In the following of this paper, the analysis is performed only at the material point scale.Moreover, small strain assumption is done and the geometrical effects are neglected.Thus the following condition is derived: This is the so called second order work criterion.σ i j and ε i j denotes respectively, the components of the Cauchy stress tensor, and of the small strain tensor.It is worth noting that those quantities are linked by their constitutive relation.In the following of the paper, bifurcation in geomaterials is investigated through this criterion.
In a first part, some analytical developments and results about the second order work are proposed for the class of incrementally piece-wise linear constitutive models.
In fact, this class of models covers the majority of those proposed in the literature for rate independent materials.Secondly, a physical interpretation of these results is given, and numerical illustrations are presented with the constitutive models of Darve (Darve et al., 1995).

Some analytical results
In this section, an analytical analysis of relation [2] is proposed.To do so, following notations are used.M is the constitutive matrix which links dσ to dε dσ = Mdε .
N is the one which links dε to dσ.If M is invertible, N = M −1 .Finally, we define du, dv the vectors of mixed linear combinations of dσ and dε and S the constitutive matrix (built with the help of M, orN) which links du and dv.This formulation allows the description of any loading path on an homogeneous sample.For example, if we want to describe an undrained triaxial compression test 1 driven with the axial force, the following problem can be built: with dq = dσ 1 − dσ 3 , and dε v = (dε 1 + 2dε 3 ).In that way, du = t dq, dε v , and dv = t [dε 1 , dσ 3 ] are conjugated variables according the second order work.By splitting matrices M and N in a symmetrical part M s and N s and a skew symmetrical part, the second order work takes the following expression : Thus, the sign of the second order work is strongly linked to the positiveness of det N s or det M s .Positiveness of the quadratic form t dσN s dσ gives a geometrical representation of the positiveness of d 2 W in the stress rate space, while t dεM s dε in the strain rate space.Let us now analyse the condition : For any incrementally piece-wise linear constitutive relation, N can take the following form in a given tensorial zone 2 and in principal axes: 1. assuming the fact that the pore water and each grain of the skeleton are incompressible, which is the case for the considered pressures 2. A tensorial zone is a domain of the loading space in which the incremental constitutive relation is linear.A classical elasto-plastic model has two tensorial zones: one for plastic loading and the second one for elastic unloading.
Hence, equation [5] reads : This equation, is the one of an elliptical cone in the stress rate space.Nevertheless, some degenerated form can exist depending on the positiveness of the N s eigenvalues.
If we assume that all eigenvalues are positive at the virgin state, and are evolving continuously with the loading parameter, four cases have to be taken into account.
3) λ 1 < 0, λ 2 > 0 : det N s < 0. Equation [7] admits an elliptical cone as solution.Second order work vanishes for loading directions upon the cone, is already negative for loading directions included inside the cone and positive elsewhere.4) λ 1 < 0, λ 2 = 0, λ 3 > 0 : det N s = 0. Equation [7] admits two secant planes as solution.Second order work vanishes on these planes, is positive in two of the subdivisions of the stress rate space, and negative in the two others.However, this solution has never been reached with constitutive models used by co-authors.
The Figure 1 illustrates the previous discussion.Thus, the directional feature of this criterion observed by Darve and Laouafa (Darve et al., 2000), (Laouafa et al., 2002) is proved here.This means that, as det N s ≤ 0 or conversely det M s ≤ 0, it exists some loading directions for which bifurcation occurs.As a remark, this directional feature is stronger for non elastic behaviour.In fact, real solutions are the ones (or part of ones) of conditions (2) to (4) which are fulfilled in the proper tensorial zone.If we call the limit of the bifurcation domain, the set of stress-strain states for which a unique unstable loading direction exists, this limit can be describe by the following condition : with n, the number of tensorial zones (of the constitutive model), u i the eigenvector corresponding to the vanishing eigenvalue, and Z i the tensorial zone considered.
To conclude with this section, we will prove that : as equation [4] impose it.In fact, we have stated that : which is equivalent to : Nicot (Nicot, 2008) proposed the following proof of relation [11].Considering the left hand side of equation [11], we have: and by considering the right hand side of equation [11], it follows: thus, as det A = det t A for any matrix A, it follows : This result means that bifurcation analysis can be performed without any restriction on the use of the constitutive relation, whereas it is not the case for the plasticity limit condition.det M = 0 is not equivalent to det N = 0.In the case of associated material, relation [10] is still valid, but condition [9] not any more, because in this case the bifurcation limit coincides with the plasticity limit condition.

Physical interpretation
In this section, numerical results, as well as a physical interpretation of analytical developments presented in the previous part are given.
Numerical results, are displayed with the constitutive models of Darve (Darve et al., 1995).Without going into details of these models, we just recall that they are not based on the classical concepts of the elasto-plasticity.Decomposition of the strain in an elastic and plastic part is not assumed, and no plastic potentials are defined.According to the first model, the non linear relation which links strain rate to stress rate is directly described by an incrementally non-linear relation.The second model is a simplification of the first one, and becomes incrementally piece-wise linear with eight tensorial zones.In principal axes, these models are written as follows : with Limit of the bifurcation domain given by equation [8], are displayed on Figure 2 for both models.As a remark, the incrementally non-linear model can be seen as incrementally piece-wise linear with an infinity of tensorial zones.Thus, with a numerical effort, an approximation of the bifurcation limit has been displayed.Then, Figure 3 presents some instability cones for stress-strain states situated beyond the bifurcation limit according to both models.These stress-strain states have been obtained along the simulation of a drained triaxial test.Now a physical interpretation is given with the formalism proposed by Nova (Nova, 1994), (Nova, 2004).In fact, what happens when d 2 W vanishes along a given loading path is presented.As the positiveness of d 2 W depends on the loading rate, stress proportional loading paths are introduced as follows : These paths allow to prospect every loading directions from a given stress-strain state.Thus, we are able to prospect some particular stress paths where the second order work vanishes.The response of such paths is displayed on Figure 4.In fact, a previous drained triaxial compression is performed until a stress-strain state situated outside of shows results obtained with the non-linear model using the numerical method.p o is the initial confining pressure, q = σ 1 − σ 3 , η = q p the bifurcation domain, then one of the proportional loading (going inside the bifurcation domain) path is followed.To simulate these loading paths, some of the variables have to be controlled.One of the possibilities is as follows : For proportional loading paths which show a peak in the ) plane, following observations can be derived.At this peak, due to the statical conditions imposed by equation [18], we obtain : Equation [20] has no trivial solutions if det S = 0. Then this equation constitutes a generalized failure rule, and the eigenvector corresponding to the vanishing eigenvalue constitutes a generalized flow rule.Furthermore, at this state, d 2 W vanishes also.
Then it can be concluded that this stress-strain state is situated inside the bifurcation  the response path related to the triaxial compression is not displayed.p = 1 3 tr σ , domain, and that the loading direction belongs to one of the cone presented above (or to the only one unstable loading direction, if the stress-strain state is situated strictly on the bifurcation domain limit).With such analysis, concept of controllability introduced by Nova (Nova, 1991), (Nova, 1994) is easily understandable.Considering previous loading path which show a peak in ε 1 − ε 3 R − R R ε 2 , if dσ 1 is controlled, this peak can be passed continuously, whereas if dε 1 − dε 3 R − R R dε 2 is controlled, effective collapse of the sample occurs at this peak.More classically, considering a drained triaxial test on dense sample, when axial strain is controlled going through the peak of q is possible, while when axial stress is controlled, effective collapse is noticed at this peak.Furthermore, at such bifurcation point (state for which d 2 W vanishes on the considered loading path), Sibille and co-authors (Sibille et al., 2007) have proved with a discrete element model that a "small" mechanical perturbation on the sample induces the effective collapse.Thus, as the second order work vanishes on a given loading path, an unstable state is reached.Effective collapse occurs if particular controlled parameters are used, or if a "small" perturbation is imposed.
To conclude with this section, we present others instability cones on Figure 5 and Figure 6 for both models of Darve.In this case, loading paths considered are not any more drained triaxial compressions.In fact, a set of loading paths going towards some given directions located in a stress deviatoric plane are considered.Then, at an arbitrary deviatoric level (situated inside the bifurcation domain), 3D instability cones are plotted.Because, in these conditions, the mean pressure and the second stress invariant are fixed, the influence of the Lode's angle on instability cones can be ob-served.These figures are 3D, but the point of view is turned along the hydrostatic axis.That's why the representation seems to be done in the deviatoric plane.Ac- 489.9 kPa, q = 3 2 dev σ .Continuous line represents the bifurcation limit 30 • .Thus, the three stress invariants plays an important role according to bifurcation, and not only the mean pressure and the deviatoric level.

Concluding remarks
In this paper, a bifurcation analysis based on the second order work criterion has been performed to analyse failure in geomaterials.The analysis was restricted to the material point scale.In introduction it has been reminded that this criterion constitutes a lower bound according to the plasticity limit and the localization criterion of Rice.In some particular cases, like the modelling of landslides under very gentle slope, this criterion is necessary in predicting failure.For example, at Petacciato in Italy, a landslide occurs under a mean slope of 6 • without any localisation pattern.The modelling using this criterion has proved to be satisfactory (LESSLOSS-Report, 2004), (Lignon et al., 2008).
In the first part, analytical results valid for any incrementally piece-wise linear model, have been stated.Existence of instability cones has been proved, and their equation given in the 3D space of stress principal axes.Then the limit of the bifurcation domain has been derived and its equation, which can be solved numerically, has been established.Finally, it has also been proved that the equation of the bifurcation domain limit is not changed when the constitutive tensor M is considered or conversely its inverse N in a dual formalism.
In the second part, 3D representation of instability cones as well as of the bifurcation domain limit have been given for Darve constitutive models.Then a physical interpretation of this criterion has been proposed.Noticeably, it has been proved that, as the second order work vanishes along a given loading path, an extension of the plastic limit criterion and of the flow rule can be defined for proper conjugated variables.Eventually, the influence of the Lode's angle on the opening of instability cones has been shown.Thus the third stress invariant plays also an important role when investigating bifurcation in geomaterials.
According to these results, second order work criterion constitutes a good tool to investigate bifurcation in geomaterials.Noticeably, it opens new perspectives for experimental testing.All proportional stress loading paths presented in this paper should be reproducible experimentally.Further more proportional strain loading paths can also be defined in the same way.Daouadji and co-authors (Darve et al., 2007) have already investigated some of these paths experimentally.One of the exciting challenge would be proposing an experimental set-up allowing the description of the bifurcation domain limit.

Figure 2 .
Figure 2. Limit of the bifurcation domain plotted in the 3D stress space for constitutive models of Darve.

Figure 3 .
Figure 3. 3D cones of unstable stress direction for a dense sand of Hostun.Figure 3.(a) presents the cones obtained with the octo-linear model.The planes represent the limit between the 8 tensorial zones, the meshes correspond to the analytical solution of Equation [7] and the point clouds the solution obtained with the numerical method.Figure 3.(b) shows results obtained with the non-linear model using the numerical method.p o is the initial confining pressure, q = σ 1 − σ 3 , η = q