Constitutive Modeling of Soils and Rocks

Preface to the English Edition

The French version of this book appeared in 2002 as part of the Material Mechanics and Engineering series. The objective of this book was to create as complete as possible a corpus of knowledge and methods in this field.

In designing this book on the mechanical behavior of soils and rocks, we gathered together a number of internationally known specialists, who each brought a significant contribution to the knowledge of the experimental behavior of these materials, as well as their constitutive modeling. Our goal was to cover as far as possible the theories at the basis of the different approaches of modeling, and also to address the most recent advances in the field.

In translating this book into English, we hope to make available to a wider scientific and engineering public the approaches and school of thought which have dominated the field of geomaterial mechanics in France over the past few decades. We have put together present-day knowledge of mechanical behavior and their theoretical bases in order to construct an original, analytical framework which, we hope, will give readers a useful guide for their own research. Most of the chapters have been updated in order to include the most recent findings on the respective topics.

Finally, we wish to dedicate this book to the memory of Professor Jean Biarez, who not only played a ground-breaking role in the history of soil mechanics in France, but remains a source of inspiration to many of us today.

Pierre-Yves Hicher

Jian-Fu Shao

Preface to the French Edition

Soils and rocks possess a number of similar characteristics: both are highly heterogenous materials formed by natural grains. This alone gives them certain rheological features which distinguish them from other solid materials, such as a strongly non-linear character, a behavior which depends on the mean stress and shearing which induces volume variations, often dilatancy, which leads to unassociated plastic strains.

Soils and rocks can be studied at different scales. At the scale of one or several grains (from μm to cm), we can examine the discrete phenomena which govern the interactions between grains. They can be described using micro-mechanical models or analyzed in order to better understand the material behavior at a larger scale, typically the size of the material specimen: this approach corresponds to passing from a discontinuous to an equivalent continuous medium. Even though the size of the latter can vary, it has to be sufficiently large (typically from 1 cm to 1 dm) compared to the size of the material discontinuities in order to be representative of the equivalent continuous medium, whose behavior can be modeled by using certain concepts of continuous medium mechanics which ignore the notion of scaling in its basic equations.

However, some phenomena, such as the development of defects or cracks within the material specimen, are located at an intermediary scale, called the meso scale. It is thus necessary, in a constitutive model for continuous medium, to use scaling techniques in order to take into account these intermediary scales. This approach, still recent but potentially strong, can also be adapted to change the scale from the material specimen to the in situ soil or rock masses in geotechnical work modeling.

The constitutive models developed to describe the mechanical behaviors at the macroscopic scale can be roughly classified into two categories: those adapted to the behavior of ductile materials and those adapted to the behavior of fragile materials. The first category corresponds mainly to sandy or clayey soils, but also to soft rocks subjected to high confining stresses. The second category corresponds mainly to hard rocks, but also to certain soft rocks and highly overconsolidated clays subjected to small confining stresses. In ductile materials, the non-linear behavior is essentially due to irreversible grain displacements, which leads to a more or less significant hardening and to a pore volume change which induces volume changes at the scale of the specimen. In fragile materials, the non-linear behavior is due to the development of cracks, whose size may vary and whose direction depends on the principal stress directions.

In order to model ductile behaviors, plasticity (elastoplasticity or viscoplasticity) has shown to be an operational framework and the large majority of the constitutive models for soils and certain soft rocks belong to this category. However, for noncohesive soils in particular, the difficulty of characterizing an elastic domain, determining the plastic mechanisms (potential and yield surface) experimentally, has led to the development of specific constitutive models, whose structure can be defined as incrementally non-linear.

In order to model fragile behaviors, the damage mechanics framework has been used to propose constitutive models adapted to describing irreversible phenomena linked to the deterioration of certain physical properties. In particular, they can take into account a large amount of rock properties: irreversible strains, dilatancy, induced anisotropy, hysteresis loop during loading-unloading due to opening and closing of mesocracks and frictional mechanisms along closed mesocracks.

In intermediary materials, the non-linear behavior can be due to microstructural changes, associating damage and hardening phenomena. Models coupling plasticity and damage have been developed to take into account this type of behavior.

After a general presentation of the constitutive models and their internal structures, each chapter will give a brief description of the different approaches mentioned above by focusing on a given class of materials. The first three chapters are devoted to the elastoplasticity theory applied to soils and soft sedimentary rocks. An alternative approach is then presented by means of the so-called incrementally non-linear models. The time-effect in clayey soils is analyzed in the framework of viscoplasticity. The behavior of hard rocks is then studied in Chapters 8 and 9, through the use of the damage theory at different scales. The modeling of the poromechanical behavior is also introduced in order to take into account the hydromechanical coupling in saturated porous rocks.

As the validity of any given model lies in its capacity to reproduce the observed material characteristics, the authors have placed the experimental data, obtained mainly from laboratory testing on intact soil and rock samples, under special consideration. The final chapter is devoted to parameter identification procedures. This is an important topic when dealing with natural materials because, each site being different from another, accurate parameter identification is essential for the quality of geotechnical work calculations, which is the final goal of this modeling approach.

Pierre-Yves Hicher

Jian-Fu Shao

Chapter 1

The Main Classes of Constitutive Relations ¹

1.1. Introduction

The study of the mechanical behavior of solid materials and its description by constitutive relations was for many years developed within the framework of isotropic linear elasticity characterized by Hooke’s law, plasticity characterized by the Von Mises, Tresca and Mohr-Coulomb criteria, and viscosity characterized in the linear case by Newton’s law. However, since the end of the 1960s, the development of more powerful numerical methods such as the finite element method and the use of high-performance computers has revived the study of material behavior, as it became possible to take into account a more realistic viscoelastoplastic modeling, albeit at the expense of much more complex formalisms.

Inside the three sets of equations defining a continuous medium mechanics problem, i.e. general equations (conservation equations), constitutive laws and boundary conditions, constitutive laws correspond to the more difficult part, particularly since the general framework in which the constitutive equations are inscribed remains often numerically imprecise. It is the comprehension of the absence of physical laws in this domain which gradually changed the designation of constitutive laws to constitutive models. The latter corresponds better to the objective of giving a mathematical form to the mechanical properties of materials, whose complexity has been demonstrated by the diversity of the experimental results.

During the last 30 years, a large variety of constitutive models have been developed and many workshops organized all over the world have shown that it is important for developers as well as users of models to be able to obtain guiding ideas and a general framework of analysis. The objective of this chapter is to try to formulate both of these.

This general framework will be more useable if it can be unified, and we intend to show that it can be applied to elastoplasticity as well to viscoplasticity or damage theory. We thus invite the reader to a wide presentation of constitutive relations for solid materials.

Two preliminary comments need to be made. First, we should explain why the chapter covers rheology in an incremental form. Two main reasons have made such an incremental presentation indispensable. The first is physical and is linked to the fact that, as soon as some plastic irreversibility is mobilized within the material, the global constitutive functional, which relates the stress state σ(t) at a given time t to the strain state ε(t) history up to this time, is in principle very difficult to formulate explicitly as this functional is singular at all stress-strain states (or more precisely non-differentiable, as will be shown). An incremental formulation enables us to avoid this fundamental difficulty. The second reason is numerical and stems from the fact that material behavior, and usually also the modeling of engineering works, exhibits many non-linearity sources which imply that the associated boundary value problem must be solved by successive steps linked to increments of loading at the boundary. Therefore, such finite element codes need to express the constitutive relations incrementally.

Our second comment concerns the use of incremental stress and strain rather than the stress and strain rates. Here also, it is the physical nature of the phenomena which determines our choice: in elastoplasticity, and more generally for all non-viscous behaviors, physical time does not play any role and, as a consequence, the derivatives with the physical time have no real meaning. Therefore, the incremental form appears to be intrinsically significant and can in fact be attached straightforwardly to the rate: the incremental strain is the product of the strain rate with the time increment, while the incremental stress is the product of Jaumann’s derivative of the stress tensor with the time increment. It is, however, incorrect to speak of stress and strain increments, since the incremental strain (for example) corresponds to a small strain variation only in the case of a sufficiently small strain.

This chapter begins with a traditional presentation of the rheological functional. We will show the limits of the functional expression and overcome this limitation by establishing the incremental rheological formalism. First, we will cover the case of non-viscous materials. The notion of tensorial zones will allow us to present the different classes of non-viscous models. Then, we will come back to the general case by considering models which take into account any kind of irreversibility.

1.2. The rheological functional

The basic concepts of continuous medium mechanics are taken for granted. The tangent linear transformation, characterized by the matrix of the gradient of the material particle positions, is assumed to describe correctly the material geometric deformation, even if some theories, called second gradient theories, consider that this first order approximation by the tangent linear transformation from the positions at a given time to the actual positions is not sufficient, and subsequently introduce second order terms [MUH 91]. We also assume that the constitutive law of a material element does not depend on the neighboring elements (some theories called non-local theories consider that the behavior of a basic material particle depends on a finite deformation field around that particle [PIJ 87]). These two hypotheses define a specific class of materials called simple media [TRU 74] for which we will develop a theoretical analysis.

The starting point of rheology is thus based upon a principle of determinism, which can be expressed as follows: if a given loading path is applied to a material sample, the material response is determined and unique, i.e., the principle of determinism applies only in conditions where there is uniqueness of the rheological response. Passing through a bifurcation point gives several possible responses. The choice of one of these responses is guided by existing imperfections which are not taken into account in the description of the material mechanical state or in the mode of loading application (control in force or in displacement, for example).

The first expression of the principle of determinism is obtained by stating that stress state σ(t) at a given time t is a functional of the history of the tangent linear transformation up to this time t. This implies that it is necessary to know the entire loading path in order to deduce the associated response path.

From a mathematical point of view, this is stated by the existence of a rheological functional F:

(1.1)

where E(t) is the strain part of the tangent linear transformation E at time t, also called the deformation gradient. Deformation gradient E is the Jacobian matrix of position f(X,t) of material point X at time t. The existence of such a functional, and not a function, is related to an essential physical characteristic: for irreversible behaviors, knowing strain ε(t) at time t does not enable us to determine the stress, and vice versa. For example, we can think of viscous or plastic materials where a given level of stress can be related to an infinite number of different strain states.

Since this chapter is devoted to the study of fundamental properties of constitutive relations, the general properties of rheological functional F need to be examined:

– isotropy of F: due to the principle of isotropy of space, F has to be an isotropic function of ε (if σ, ε and the internal variables are subjected to an equal rotation, F remains identical);

– non-linearity of F: the hypothesis of linearity for F is expressed by:

In such a case, the material response to a sum of histories will simply be equal to the sum of the responses to each history. This constitutes Boltzmann’s principle and is the basis of linear viscoelasticity theory, but it is not at all valid in the general case, as in elastoplasticity, for example, where, when we double the strain for example, the stress is obviously not doubled, due to the non-linear behavior.

In the general case, F must be studied in the framework of non-linear functionals:

  F is furthermore non-differentiable as soon as there is some plastic irreversibility. Owen and Williams (1969) showed in fact that the assumptions of non-viscosity and differentiability of stress functional F imply that there is no internal dissipation.

In other words, a non-viscous material whose constitutive functional is differentiable is necessarily elastic. Basically, this is due to the fact that, in plasticity, the tangent loading modulus is not equal to the tangent unloading modulus. Therefore, if we want to describe the behavior of anelastic materials by using a stress-strain relationship, this relation must be formulated using a non-linear and non-differentiable functional:

– degeneration of F: the only case of degeneration of functional F into a function corresponds to elasticity (possibly non-linear and anisotropic), where there is a one-to-one mapping between stress and strain.

Finally, if we want to describe irreversible behavior, we have to consider a nonlinear, non-differentiable functional, which, mathematically, is very difficult to use. We therefore need to study constitutive relations using an incremental formulation rather than a global one.

1.3. Incremental formulation of constitutive relations

We shall now introduce an incremental formulation using a second statement from the principle of determinism. The second principle of determinism, which can be called in the small to be distinguished from the first principle in the large, is obtained by stating that a small load applied during a time increment dt induces a small uniquely determined response.

As stated previously, this principle applies only if the uniqueness of the incremental constitutive relation is maintained. For bifurcation cases, the principle is no longer valid and the choice of the bifurcated branch will depend on the boundary conditions and material imperfections. In addition, the principle assumes implicitly that the loading rate is kept constant during the time increment (even if it can vary from one increment to another), which excludes dynamic loads due to shocks.

We denote by dε = D dt the incremental strain tensor of order two equal to the product of the second-order strain rate tensor, D (symmetric part of transformation rate L.: L = E-1 where E is here the deformation gradient) and time increment dt, and by dσ = σˆdt the incremental stress tensor, equal to the product of an objective time derivative of Cauchy stress tensor σ and dt.

(1.2)

What are the properties of this tensorial function Fh?

The first comment concerns the fact that Fh depends on the previous stress-strain history. This history is generally characterized by some scalar and tensorial variables denoted by h which will appear as parameters in the previous relation. These parameters describe, as far as possible, the actual deformed state of the solid. According to various constitutive theories, they are sometimes called memory variables, hardening parameters, internal variables, etc.

Secondly, Fh must satisfy the objectivity principle, which means that Fh must be independent of any observer movement relative to the solid. Thus, Fh is an isotropic function of all its arguments: dε, dσ and also the state tensorial variables, which characterize its presently deformed state. However, if the material is anisotropic insofar as its mechanical properties are concerned, then Fh is an anisotropic function of dε and dσ.

Finally, Fh is essentially a non-linear function as long as there is some plastic irreversibility. If Fh is linear, we can write:

which is the general form of viscoelastic laws where M is the fourth-order elastic tensor and C the second-order creep rate tensor of the material.

This property of non-linearity for Fh is directly linked to the non-differentiability of rheological functional F, the property of differentiability of F being equivalent to the linearity of Fh.

In conclusion, relation (1.2) corresponds to the general incremental form of the constitutive relations. We will now distinguish between viscous and non-viscous materials in order to represent this incremental form more precisely.

1.4. Rate-independent materials

For non-viscous materials, the loading rate (characterized by time gradation on a loading path) has no influence on material constitutive behavior: a given loading path, followed at any given rate, gives the same response path. In other words, the behavior class considered is rate-independent. This restriction of the constitutive law implies that constitutive function Fh, which relates dε and dσ, is independent of time increment dt, during which the incremental loading is applied. Therefore, Fh is independent of dt and we can write:

(1.3)

or

(1.4)

The possibility of inversing G or H is linked to the uniqueness of the constitutive relations. This question will not be studied here; for more details see [DAR 94, DAR 95a].

From a mathematical point of view this independence of non-viscous behaviors on loading rates implies the following identity:

thus

(1.5)

which states that if the stress rate is multiplied by any positive scalar, the strain rate response is also multiplied by the same scalar.

This is the first property of G: G is a homogenous function of degree 1 in dσ with respect to the positive values of the multiplying parameter. This homogenity property must not be confused with that of positively homogenous functions, which is given by;

In addition to this property of homogenity of degree 1, functions G and H, as we have seen in general for function Fh, are non-linear and anisotropic in dσ (or in dε).

1.4.1. Non-linearity of G and H

If, in relation (1.3), dε is the response to an incremental loading dσ, the response to an incremental loading – dσ, following dσ, is not equal to – dε, because plastic irreversibility or damage takes place in the material. Therefore, G and H are necessarily non-linear functions of dσ and dε respectively, which implies that the principle of incremental superposition cannot be rigorously verified, except within the elastic domain, or more generally within a domain of incremental linearity of the constitutive model. Calculus shows, however, that the principle of incremental superposition can be roughly verified along step-wise paths, approaching a given loading path [DAR 95b].

1.4.2. Anisotropy of G and H

Following the same reasoning as for function Fh, we can deduce that G and H are anisotropic functions of dσ and dε respectively.

This anisotropy is directly linked to the geometrical meso-structure of the material, which is gradually modified by the strain (particularly irreversible) history. We have seen that this history can be characterized by scalar and tensorial state parameters.

In simple cases, this anisotropy is directly imposed by the choice of these state parameters. If we consider, for example, only scalar memory parameters (such as void ratio), defined independently of any frame, based on the objectivity principle functions G and H will be isotropic functions, which is not supported by experiments.

If we add one single tensor variable (such as the stress tensor) to these scalar memory parameters, G and H are orthotropic functions of dσ and dε respectively, the orthotropy axes being identified with the principal stress or strain axes. In this case, it means that G and H are invariant by symmetry with respect to any plane containing two principal stress or strain directions.

In the more general case of state variables with at least two second-order non-commutating tensorial variables, anisotropy is not defined. Orthotropy thus becomes a constitutive assumption, which must be considered as an approximation of the real behavior of the material for classes of loading in which stress and strain principal axes rotate.

1.4.3. Homogenity of degree 1 of G and H

Having described the three main properties of G, we will now focus on the first (homogenity of degree 1) to see the mathematical consequences of such a property. Let us for this purpose recall Euler’s identity for homogenous regular functions of degree 1, by writing it for a function of two variables:

(1.6)

where partial derivatives δf/δx and δf/δy are homogenous functions of degree 0.

In formulating constitutive relations, it is often more convenient to replace stress tensor σ and strain tensor ε, which are second order, by two vectors of IR⁶ defined in a six-dimensional related space. In this space, vectorial function G is written:

with summation on index β.

The partial derivatives of a homogenous function of degree 1 are homogenous functions of degree 0. Therefore, functions ∂Gα/∂(dσβ) depend only on the direction of dσ, characterized by the unit vector:

with:

Finally, we obtain:

(1.7)

or:

(1.8)

Equations (1.7) and (1.8) are the general expressions for all rate-independent constitutive relations. Constitutive tensors M and N also depend on state variables and memory parameters, which characterize the loading history. These two matrices are the gradient matrices of non-linear functions G and H, respectively. In that sense, they can be considered as tangent constitutive tensors and are therefore uniquely defined. However, it is possible to construct from them an infinite number of secant constitutive tensors by adding to the M or N lines the components of any unit vector perpendicular to dσ or dε, respectively.

Relations (1.7) and (1.8) will now allow us to propose a classification of all the existing rate-independent constitutive relations with respect to their intrinsic structure.

1.5. Notion of tensorial zones

First of all, we need to define the notion of a tensorial zone [DAR 82]. We will call a tensorial zone any domain in the incremental loading space on which the restriction of G or H is a linear function. In other words, the relationship between dε and dσ in a given tensorial zone is incrementally linear. If we denote the tensorial zone being considered as Z, the following definition implies:

In zone Z, the constitutive relation is characterized by a unique tensor Mz. If u belongs to Z, any vector collinear to u also belongs to Z for all real positive values. Therefore, a zone is defined by a set of half-infinite straight lines, whose apex is the same and is at the origin of the incremental loading space. Tensorial zones thus comprise adjacent hypercones, whose common apex is this origin. What does the constitutive relation become on the common boundary of two (or several) adjacent tensorial zones? If Mz1 and Mz2 are constitutive tensors attached, respectively, to tensorial zones Z1 and Z2, they must obviously satisfy the condition of continuity of the response to any loading direction u which belongs to the boundary between Z1 and Z2:

(1.9)

Relation (1.9) can be called a continuity condition for zone change. This condition prohibits, in particular, an arbitrary choice of the constitutive tensors in two adjacent tensorial zones.

Furthermore, we will see that conventional elastoplastic relations satisfy this condition by means of the consistency condition. This is also the case for damage models when they are built in a rigorous manner. On the other hand, hypoelastic models do not necessarily fulfil this condition, which has to be verified a posteriori. It has been proven that this is not the case for some of these models [GUD 79].

The response-envelopes, as proposed by Gudehus [GUD 79], constitute geometrical diagrams which completely characterize a constitutive relation at a given stress strain state after a given strain history. At this state, all the incremental loadings, having the same norm but oriented in all directions, are considered and all the incremental responses are plotted. The extremities of the response vectors form a hypersurface which is called the response-envelope. Figure 1.1 gives an example of an elastoplastic model in axisymmetric condition: the continuity of the response at the boundary of two tensorial zones appears well fulfilled. Figure 1.2 gives an example of a model with discontinuities, whereas Figure 1.3 corresponds to a continuous non-elastoplastic model.

In fact, the number of tensorial zones characterizes how a given model describes the irreversibility due to plasticity or damage, and the directional change of behavior, i.e. how constitutive tensor M (or N) evolves with the direction of loading u (or v). More precisely, the number of tensorial zones of a given constitutive model is an intrinsic criterion, which fully represents the model structure. Therefore, we have chosen this criterion to classify, in the next section, the different rate-independent constitutive models.

1.6. The main classes of rate-independent constitutive relations

1.6.1. Constitutive relations with one tensorial zone

The first class of relations that we are going to look at is related to the simplest assumption that there is only one tensorial zone. Therefore:

Therefore:

or:

The behavior is therefore entirely reversible (except possibly in the case of the existence of memory parameters h, but this corresponds rather to an artefact in hypoelasticity). As there is a unique linear relationship between dε and dσ (incremental linearity), we have here in this first class all the elastic laws, isotropic or anisotropic, linear or non-linear (in this last case, M and N depend on the actual state of stress).

The best way to reproduce an elastic behavior (without any internal dissipation) in a rigorous manner is to introduce an elastic potential V defined by:

As V is an exact total differential, we obtain the following expression:

Therefore:

using Schwarz’s identity. As a consequence, matrices M and N are symmetric and tensors C and D, defined by

have major symmetries

In the general case of non-linear elasticity, the existence of a potential also implies conditions of integrability [LOR 85], which have to be satisfied by the components of M and N. All these laws are called hyperelastic, while, in the absence of a potential, they are called hypoelastic. The hypoelastic models generate energy dissipation, and should thus not be used in practice, the behavior represented by these models being poorly identified.

1. 6. 2. Constitutive relations with two tensorial zones

In the presence of two tensorial zones, we can call one the loading zone and the other the unloading zone. We thus define two different behaviors (two different constitutive tensors), one representing the loading condition, and the other the unloading condition. Each matrix is attached to a different tensorial zone, these two tensorial zones being separated by a hyperplane in dσor dε space. A loading-unloading criterion, a linear and homogenous inequation in dσ or dε, allows us to discriminate between the two behaviors. The hyperplane equation corresponds, by construction, to the zero value of the loading-unloading criterion. The continuity condition at the crossing of the hyperplane gives a link to the two constitutive tensors and the hyperplane equation:

(1.10)

Figure 1.1. Response envelopes [GUD 79] in axisymmetric conditions for an elastoplastic material with two tensorial zones (characterized by fs = ±1) for three different stress levels. The continuity of the response envelopes is verified

Numerous constitutive models follow these general rules and are therefore based on the definition of two tensorial zones. Their formalisms are basically similar, even if, sometimes, the detailed equations do not clearly show their fundamentally bilinear structure. These models are divided into three different families: elastoplastic models with one plastic potential, hypoelastic models with a unique loading-unloading criterion and damage laws. We will now examine them successively.

1.6.2.1. Elastoplastic models with one plastic potential

The first assumption concerns the additive decomposition of the incremental strain into an elastic part (reversible) and a plastic part (irreversible):

(1.11)

The plastic deformations exist only beyond a given limit surface, the elastic limit, which depends on the loading history and evolves due to the hardening created by plastic strains, as has been shown experimentally. Its equation is given by:

(1.12)

The loading condition is obtained by writing that the incremental stress is directed outwards from the elastic limit. The unloading condition is obtained if the incremental stress is directed inwards. It follows that:

(1.13)

The equation of the hyperplane, the border between the two zones in the dσ space, is thus given by:

When the elastic limit is reached, the direction of the incremental plastic strain is given by the flow rule which is often specified in terms of a plastic potential g(σ) as:

(1.14)

where dλ is an arbitrary scalar, whose value is determined by the consistency rule which mathematically expresses that,