A Mathematical Approach to Establishing Constitutive Models for Geomaterials

Abstract

The mathematical foundation of the traditional elastoplastic constitutive theory for geomaterials is presented from the mathematical point of view, that is, the expression of stress-strain relationship in principal stress/strain space being transformed to the expression in six-dimensional space. A new framework is then established according to the mathematical theory of vectors and tensors, which is applicable to establishing elastoplastic models both in strain space and in stress space. Traditional constitutive theories can be considered as its special cases. The framework also enables modification of traditional constitutive models.

1. Introduction

The mechanical properties of geomaterials are complex and essential to make a numerical prediction; thus, many researchers have paid attention to constitutive relations of geomaterials. The simplest constitutive model for geomaterials is the elastic model, among which the common nonlinear models are the Cauchy elastic model, the hyperelastic model, and the hypoelastic model. The Cauchy elastic model assumes that the stress (or strain) in the material depends on the current strain (or stress) only, and not on its history. The constitutive equation for the hyperelastic model is established by the strain energy function or complement energy function. The hypoelastic model assumes that the stress state of an elastic material is associated with both the strain state and the stress path. The typical hypoelastic models are the E-μ and E-B models proposed by Duncan et al. [1, 2] and the K-G model [3–5].

According to the experimental results, most deformations of geomaterials are plastic deformations. Therefore, traditional plasticity theory has often been used to establish constitutive models for soil. For example, Drucker et al. [6] described the deformation property of soil by traditional elastoplastic theory and proposed a model with conical yield surface affected by hydrostatic pressure. Roscoe et al. [7] proposed a plastic cap model for normally consolidated clay, which is well known as the Cambridge model. Subsequently, Roscoe and Burland [8] modified the dilatancy equation in the Cambridge model and proposed a modified Cambridge model with elliptical yield surface. Wroth and Bassett [9] and Poorooshasb et al. [10] extended the model to sandy soil, and Yao et al. [11, 12] extended the model to sandy soil and overconsolidated soil. There are plenty of elastoplastic models, such as models with a single yield surface proposed by Desai et al. [13, 14] and Lade et al. [15–18], models with a double yield surface, and three surface models [19]. The concepts of the bounding surface [20–23] and the subloading surface [24, 25], endochronic theory [26], and disturbed states [27] have also been applied to establishing constitutive models for geomaterials.

In this paper, a theoretical framework on establishing constitutive models for geomaterials is proposed, the initial thought of which is provided by the first author in 1988 and in 1990s [28–31], and it has been implemented by some researchers to simulate the behavior of jointed rock masses [32] and soil-structure interface [33].

2. Classical Elastoplastic Theory of Geomaterials

The incremental form of a stress-strain relationship in traditional geomechanics is generally expressed as d𝜎_𝑖𝑗=𝐷_{𝑖𝑗𝑘𝑙}d𝜀_𝑘𝑙. (1) Determining 𝐷_{𝑖𝑗𝑘𝑙} is the major topic for constitutive models of geomaterials. Obviously, 𝐷_{𝑖𝑗𝑘𝑙} can be obtained by fitting experimental data, given that experiments on stress and strain tensors are conducted. However, it is extremely difficult to do so. Therefore, experiments on the stress-strain relationship of geomaterials are usually conducted in principal stress/strain space; that is, only the relationship between the principal stress 𝜎_𝑖(𝑖=1,2,3) and the principal strain 𝜀_𝑖(𝑖=1,2,3), is obtained. To obtain 𝐷_{𝑖𝑗𝑘𝑙}, the constitutive tensor in general coordinate space should be derived from the stress-strain relationship in principal stress/strain space. From a mathematical point of view, these can be treated as the problems of coordinate transformation [34–36].

The relationship between the plastic strain increment and stress increment in principal stress/strain space is defined as d𝜀^𝑝_𝑖3×1=[𝐴]_3×3d𝜎_𝑖3×1,[𝐴]_3×3=𝑎₁₁𝑎₁₂𝑎₁₃𝑎₂₁𝑎₂₂𝑎₂₃𝑎₃₁𝑎₃₂𝑎₃₃, (2a)(2b)where 𝑎_𝑖𝑗 are the functions of total stress 𝜎_𝑖, total strain 𝜀_𝑖(𝑖=1,2,3) or stress path.

When the matrix rank of [𝐴] is 1, or |𝐴|=0, there exists a vector {𝛼₁𝛼₂𝛼₃}^𝑇 and coefficients 𝛽₁,𝛽₂,𝛽₃ to express [𝐴] as [𝐴]=𝛽₁𝛽₂𝛽₃^𝑇𝛼₁𝛼₂𝛼₃, (3) Therefore, substituting (3) into (2a) and (2b) givesd𝜀^𝑝_𝑖=𝛽₁𝛽₂𝛽₃^𝑇𝛼₁𝛼₂𝛼₃d𝜎_𝑖, (4a) that is, d𝜀^𝑝_𝑖=d𝜆𝛽₁𝛽₂𝛽₃^𝑇, (4b)where d𝜆=³_𝑖=1𝛼_𝑖d𝜎_𝑖. (5) According to (4b), d𝜀^𝑝₁∶d𝜀^𝑝₂∶d𝜀^𝑝₃=𝛽₁∶𝛽₂∶𝛽₃. (6)𝛽_𝑖(𝑖=1,2,3) is a function of 𝜎_𝑖 or 𝜀_𝑖. When 𝜷={𝛽₁𝛽₂𝛽₃}^𝑇 is of a field with potential, there is a potential function 𝑄 such that 𝛽_𝑖=𝜕𝑄𝜕𝜎_𝑖. (7) Substituting (7) into (4b), we have d𝜀^𝑝_𝑖=d𝜆𝜕𝑄𝜕𝜎_𝑖. (8) If we assume that d𝜀^𝑝_𝑖 and 𝜎_𝑖 have the same principal directions, the coordinate transformation can be expressed as follows: d𝜀^𝑝_𝑖𝑗=d𝜀^𝑝_𝑖𝜕𝜎_𝑖𝜕𝜎_𝑖𝑗. (9) Substituting (9) into (8) gives d𝜀^𝑝_𝑖𝑗=d𝜆𝜕𝑄𝜕𝜎_𝑖𝑗. (10)

Similarly, from elastic potentials theory, there is a potential function 𝑊 in principal stress space, and 𝜀_𝑖 is defined as 𝜀_𝑖=𝜕𝑊𝜕𝜎_𝑖. (11) If we assume that 𝜎_𝑖 and 𝜀_𝑖 have the same principal direction, the coordinate transformation can be expressed as follows: 𝜀_𝑖𝑗=𝜀_𝑖𝜕𝜎_𝑖𝜕𝜎_𝑖𝑗. (12) Substituting (11) into (12), 𝜀_𝑖𝑗=𝜕𝑊𝜕𝜎_𝑖𝑗. (13)

In conclusion, traditional plastic potential theory corresponds to the case that the matrix [𝐴] in (2a) and (2b) has rank 1, and 𝜷 can be expressed as the gradient vector of a potential function. Based on mathematical principles, a more general potential function-based constitutive framework can be established according to vector field theory and tensor theory as described below.

3. Derivation of Constitutive Framework from Vector Field Theory

Obviously, when the three principal components of the plastic strain increment, d𝜀^𝑝_𝑖(𝑖=1,2,3), are considered to be components of a vector 𝐝𝜺^𝐩, the principal components can be expressed as three linearly independent 3D vectors by a vector fitting method. The gradient vectors of three linearly independent potential functions are selected as the linearly independent vectors.

When 𝐝𝜺^𝐩 is expressed in principal stress space, the coordinate orientations of 𝐝𝜺^𝐩 and 𝜎_𝑖 are the same, and Φ₁,Φ₂,Φ₃ are three linearly independent potential functions in principal stress space, and then the following expression is obtained: d𝜀^𝑝_𝑖=³_𝑘=1𝜆_𝑘𝜕Φ_𝑘𝜕𝜎_𝑖, (14) where 𝜆_𝑘(𝑖=1,2,3) are the coefficients. Suppose that the principal directions of 𝐝𝜺^𝐩 and 𝜎_𝑖 are the same, and we substitute (14) into (9), giving the tensor expression in general coordinate space as d𝜀^𝑝_𝑖𝑗=³_𝑘=1𝜆_𝑘𝜕Φ_𝑘𝜕𝜎_𝑖𝑗. (15)

Similarly, 𝐝𝜺^𝐩 can also be expressed in strain space. Let Ψ₁,Ψ₂,Ψ₃ be three linearly independent potential functions in principal strain space, and the following expression is obtained: d𝜀^𝑝_𝑖𝑗=³_𝑘=1𝜇_𝑘𝜕Ψ_𝑘𝜕𝜀_𝑖𝑗. (16) Define the plastic stress increment as d𝜎^𝑝=𝐷_𝑒d𝜀^𝑝, (17) where [𝐷_𝑒] is the elastic stiffness matrix, and then the expressions in stress space and strain space are d𝜎^𝑝_𝑖𝑗=³_𝑘=1𝛽_𝑘𝜕𝐺_𝑘𝜕𝜎_𝑖𝑗,d𝜎^𝑝_𝑖𝑗=³_𝑘=1𝛼_𝑘𝜕𝐹_𝑘𝜕𝜀_𝑖𝑗. (18)(19)

For the total stress and the total strain, consider the three principal stresses and the three principal strains as vectors in three-dimensional space with the same principal directions, and the following expressions are obtained: 𝜎_𝑖𝑗=³_𝑘=1𝜂_𝑘𝜕𝑊_𝑘𝜕𝜀_𝑖𝑗,𝜀_𝑖𝑗=³_𝑘=1𝜒_𝑘𝜕Ω_𝑘𝜕𝜎_𝑖𝑗, (20) where 𝑊_𝑘,Ω_𝑘(𝑖=1,2,3) are potential functions with linearly independent gradient vectors in strain space and stress space, respectively.

4. Derivation of Constitutive Framework from Tensor Theory

If 𝐴_𝑖𝑗 and 𝐸_𝑖𝑗 are symmetric second-order tensors with the same principal directions, the following equations are obtained according to tensor theory and vector fitting: 𝐴_𝑖𝑗=³_𝑘=1𝜆_𝑘𝜕𝐼_𝐸𝑘𝜕𝐸_𝑖𝑗,𝐸_𝑖𝑗=³_𝑘=1𝛽_𝑘𝜕𝐼_𝐴𝑘𝜕𝐴_𝑖𝑗, (21)(22) where 𝐼_𝐸𝑘(𝑘=1,2,3) are three independent invariants of 𝐸_𝑖𝑗, and 𝐼_𝐴𝑘(𝑘=1,2,3) are three independent invariants of 𝐴_𝑖𝑗. For example, for the stress tensor 𝜎_𝑖𝑗, the three independent invariants can be 𝐼_𝜎1=𝜎_𝑖𝑖, 𝐼_𝜎2=(1/2)𝜎_𝑖𝑗𝜎_𝑗𝑖,  𝐼_𝜎3=(1/3)𝜎_𝑖𝑘𝜎_𝑘𝑛𝜎_𝑛𝑚;  𝐼_𝜎1=𝑝,  𝐼_𝜎2=𝑞,  𝐼_𝜎3=𝜃; or 𝐼_𝜎1=𝜎₁,  𝐼_𝜎2=𝜎₂,  𝐼_𝜎3=𝜎₃, where 𝑝 is the mean stress, 𝑞 is the deviatoric stress, 𝜃 is the Lode’s angle, and 𝜎₁,  𝜎₂, and 𝜎₃ are the three principal stresses.

If 𝐴_𝑖𝑗=d𝜀^𝑝_𝑖𝑗,  𝐸_𝑖𝑗=𝜎_𝑖𝑗,  𝐼_𝜎𝑘(𝑘=1,2,3) are invariants of 𝜎_𝑖𝑗, it can be obtained from (21) that d𝜀^𝑝_𝑖𝑗=³_𝑘=1𝜆_𝑘𝜕𝐼_𝜎𝑘𝜕𝜎_𝑖𝑗, (23) which is equivalent to (15) when Φ_𝑘=𝐼_𝜎𝑘.

If 𝐴_𝑖𝑗=d𝜀^𝑝_𝑖𝑗, 𝐸_𝑖𝑗=𝜀_𝑖𝑗,  𝐼_𝜀𝑘(𝑘=1,2,3) are invariants of 𝜀_𝑖𝑗, then d𝜀^𝑝_𝑖𝑗=³_𝑘=1𝛼_𝑘𝜕𝐼_𝜀𝑘𝜕𝜀_𝑖𝑗. (24) If 𝐴_𝑖𝑗=𝜀_𝑖𝑗,   𝐸_𝑖𝑗=𝜎_𝑖𝑗, then 𝜎_𝑖𝑗=³_𝑘=1𝛽_𝑘𝜕𝐼_𝜀𝑘𝜕𝜀_𝑖𝑗,𝜀_𝑖𝑗=³_𝑘=1𝜆_𝑘𝜕𝐼_𝜎𝑘𝜕𝜎_𝑖𝑗. (25) Equation (25) corresponds to (20), respectively, when 𝑊_𝑘=𝐼_𝜀𝑘,Ω_𝑘=𝐼_𝜎𝑘.

Similarly, if 𝐴_𝑖𝑗=Δ𝜀_𝑖𝑗,  𝐸_𝑖𝑗=Δ𝜎_𝑖𝑗, then according to (21), we have Δ𝜀_𝑖𝑗=³_𝑘=1𝜆_𝑘𝜕𝐼_𝜎𝑘𝜕Δ𝜎_𝑖𝑗. (26)

It should be noted that the derivations from vector field theory and from tensor theory are actually coincident. The potential functions in the derivation from vector field theory are functions of invariants, which degenerate to tensor form after composite derivation.

5. Elastoplastic Matrix in Stress Space

Without considering the effect of the Lode’s angle 𝜃 and rotation of principal stress, the relations for the plastic strain increment and stress increment are expressed asd𝜀^𝑝_𝑣d=𝐴d𝑝+𝐵d𝑞,𝜀^𝑝=𝐶d𝑝+𝐷d𝑞, (27a)(27b)where 𝐴,   𝐵,  𝐶,   𝐷 are parameters.

The following equation is obtained from the potential functions Φ₁=𝑝,   Φ₂=𝑞 in (15) without considering the effect of the Lode’s angle: d𝜀^𝑝_𝑖𝑗=𝜆₁𝜕𝑝𝜕𝜎_𝑖𝑗+𝜆₂𝜕𝑞𝜕𝜎_𝑖𝑗. (28) In matrix form, this is d𝜀^𝑝=𝜆₁𝜕𝑝𝜕𝜎+𝜆₂𝜕𝑞.𝜕𝜎 (29) It is obvious that 𝜆₁=d𝜀^𝑝_𝑣,𝜆₂=d𝜀^𝑝. (30) According to (27a), (27b), (29), and (30), d𝜀^𝑝=𝐴d𝑝+𝐵d𝑞𝜕𝑝+𝜕𝜎𝐶d𝑝+𝐷d𝑞𝜕𝑞.𝜕𝜎 (31) Since d𝑝=𝜕𝑝𝜕𝜎^𝑇{d𝜎},d𝑞=𝜕𝑞𝜕𝜎^𝑇{d𝜎}, (32) it follows thatd𝜀^𝑝=𝐶^𝜎_𝑝},𝐶{d𝜎_𝜎𝑝=𝐴𝜕𝑝𝜕𝜎𝜕𝑝𝜕𝜎^𝑇+𝐵𝜕𝑝𝜕𝜎𝜕𝑞𝜕𝜎^𝑇+𝐶𝜕𝑞𝜕𝜎𝜕𝑝𝜕𝜎^𝑇+𝐷𝜕𝑞𝜕𝜎𝜕𝑞𝜕𝜎^𝑇, (33a)(33b)where [𝐶^𝜎_𝑝] is the plastic compliance matrix.

Since {d𝜀}={d𝜀^𝑒}+{d𝜀^𝑝}, then{d𝜀}=𝐶^𝜎_𝑒𝑝{d𝜎}, (34a) where 𝐶^𝜎_𝑒𝑝=𝐶_𝑒+𝐶^𝜎_𝑝. (34b)[𝐶_𝑒] is the elastic compliance matrix and [𝐶^𝜎_𝑒𝑝] is the elastoplastic compliance matrix.

Therefore, the elastoplastic model in stress space can be established once the parameters 𝐴, 𝐵,  𝐶, and 𝐷 are determined. Note that 𝐴,  𝐵,  𝐶, and 𝐷 are not constants, which evolve with stress and strain, as presented in the following sections.

6. Elastoplastic Matrix in Strain Space

According to (17), the plastic stress increment is defined as {d𝜎^𝑝}=[𝐷_𝑒]{d𝜀^𝑝}, where [𝐷_𝑒] is the elastic matrix and {d𝜀^𝑝} is the plastic strain increment.

𝜀_𝑣,  𝜀, and 𝜓 are three invariants of the strain tensor, where 𝜓 is the strain Lode’s angle. If we ignore the effect of Lode’s angle and take 𝜀_𝑣 and 𝜀 as potential functions, that is, 𝐹₁=𝜀_𝑣,𝐹₂=𝜀 in (19), then d𝜎^𝑝_𝑖𝑗=𝜆₁𝜕𝜀_𝑣𝜕𝜀_𝑖𝑗+𝜆₂𝜕𝜀𝜕𝜀_𝑖𝑗. (35) Written in matrix form, d𝜎^𝑝=𝜆₁𝜕𝜀_𝑣𝜕𝜀+𝜆₂𝜕𝜀.𝜕𝜀 (36) Obviously, 𝜆₁=d𝑝^𝑝,𝜆₂=d𝑞^𝑝. (37) From the definition of (17),d𝑝^𝑝=𝐾_𝑒d𝜀^𝑝_𝑣=𝐾_𝑒d𝜀_𝑣−d𝜀^𝑒_𝑣=𝐾_𝑒d𝜀_𝑣−d𝑝,d𝑞^𝑝=3𝐺_𝑒d𝜀^𝑝=3𝐺_𝑒d𝜀−d𝜀^𝑒=3𝐺_𝑒d𝜀−d𝑞, (38a)(38b)where 𝐾_𝑒,  𝐺_𝑒 are the elastic bulk modulus and shear modulus.

Consideringd𝜀^𝑒_𝑣=1𝐾_𝑒dd𝑝,𝜀^𝑒=13𝐺_𝑒d𝑞, (39a)(39b)substituting (39a) and (39b) into (27a) and (27b) yieldsd𝜀_𝑣=1𝐾_𝑒d+𝐴d𝑝+𝐵d𝑞,1𝜀=𝐶d𝑝+3𝐺_𝑒+𝐷d𝑞. (40a)(40b)It can now be calculated thatd𝑝=𝐴d𝜀_𝑣+𝐵d𝜀,d𝑞=𝐶d𝜀_𝑣+𝐷d𝜀, (41a)(41b)where 1𝐴=|𝐴|1𝐷+3𝐺_𝑒,𝐵𝐵=−|𝐴|,𝐶𝐶=−|𝐴|,1𝐷=|𝐴|1𝐴+𝐾_𝑒,|𝐴|=3𝐷𝐺_𝑒+𝐴𝐾_𝑒+13𝐾_𝑒𝐺_𝑒+(𝐴𝐷−𝐵𝐶). (42) Substituting (41a) and (41b) into (37), (38a), and (38b) gives𝜆₁=d𝑝^𝑝=𝐾_𝑒−𝐴d𝜀_𝑣−𝐵d𝜆𝜀,₂=d𝑞^𝑝=−𝐶d𝜀_𝑣+3𝐺_𝑒−𝐷d𝜀. (43a)(43b)In addition, d𝜀_𝑣=𝜕𝜀_𝑣𝜕𝜀^𝑇{d𝜀},d𝜕𝜀=𝜀𝜕𝜀^𝑇{d𝜀}. (44) Substituting (44) into (43a), (43b), and (36) yieldsd𝜎^𝑝=𝐷^𝜀_𝑝{d𝜀}, (45a) where 𝐷^𝜀_𝑝=𝐾_𝑒−𝐴𝜕𝜀_𝑣𝜕𝜀𝜕𝜀_𝑣𝜕𝜀^𝑇−𝐵𝜕𝜀_𝑣𝜕𝜕𝜀𝜀𝜕𝜀^𝑇−𝐶𝜕𝜀𝜕𝜀𝜕𝜀_𝑣𝜕𝜀^𝑇+3𝐺_𝑒−𝐷𝜕𝜀𝜕𝜕𝜀𝜀𝜕𝜀^𝑇, (45b) that is, 𝐷^𝜀_𝑝=1|𝐴|𝛼𝐴𝜕𝜀_𝑣𝜕𝜀𝜕𝜀_𝑣𝜕𝜀^𝑇+𝐵𝜕𝜀_𝑣𝜕𝜕𝜀𝜀𝜕𝜀^𝑇𝜕+𝐶𝜀𝜕𝜀𝜕𝜀_𝑣𝜕𝜀^𝑇+𝐷𝛼𝜕𝜀𝜕𝜕𝜀𝜀𝜕𝜀^𝑇, (45c)where 𝛼=𝐾_𝑒/3𝐺_𝑒=2(1+𝜇_𝑒)/9(1−2𝜇_𝑒), and 𝜇_𝑒 is the elastic Poisson ratio.

Hence, the total stress increment can be expressed as{d𝜎}=𝐷_𝑒d𝜀^𝑒=𝐷_𝑒{d𝜀}−d𝜀^𝑝=𝐷_𝑒{d𝜀}−d𝜎^𝑝=𝐷_𝑒{d𝜀}−𝐷^𝜀_𝑝}=𝐷{d𝜀_𝜀𝑒𝑝{d𝜀}, (46a) where 𝐷^𝜀_𝑒𝑝=𝐷_𝑒−𝐷^𝜀_𝑝. (46b)The duality of stress and strain is evident in (45a), (45b), (45c), (33a), and (33b).

It should be noted that it is practically impossible to obtain the total strain in soil, and thus the elastoplastic matrix in stress space is more applicable. However, if we could further extend the framework to the space of strain increment, the practicability becomes promising.

7. Relationship with Traditional Elastoplastic Model

7.1. General Form

The elastoplastic compliance matrix of the traditional elastoplastic model is 𝐶^𝜎_𝑒𝑝=𝐶_𝑒+1𝐴_𝐻𝜕𝑔𝜕𝜎𝜕𝑓𝜕𝜎^𝑇, (47) where 𝑓 and 𝑔 are the yield function and plastic potential function, 𝐴_𝐻=−𝜕𝑓/𝜕𝐻{𝜕𝐻/𝜕𝜀^𝑝}^𝑇{𝜕𝑔/𝜕𝜎} is the plastic hardening modulus, and 𝐻 is the hardening parameter.

𝑓 and 𝑔 are usually expressed in terms of the stress invariants, 𝑝,  𝑞,  𝜃. If the effect of the Lode’s angle 𝜃 is not considered, the expression only concerns 𝑝 and 𝑞, that is,𝜕𝑓=𝜕𝜎𝜕𝑓𝜕𝑝𝜕𝑝+𝜕𝜎𝜕𝑓𝜕𝑞𝜕𝑞,𝜕𝜎𝜕𝑔=𝜕𝜎𝜕𝑔𝜕𝑝𝜕𝑝+𝜕𝜎𝜕𝑔𝜕𝑞𝜕𝑞.𝜕𝜎 (48a)(48b)Substituting (48a) and (48b) into (47), we have 𝐶^𝜎_𝑒𝑝=𝐶_𝑒+1𝐴_𝐻×𝜕𝑔𝜕𝑝𝜕𝑓𝜕𝑝𝜕𝑝𝜕𝜎𝜕𝑝𝜕𝜎^𝑇+𝜕𝑔𝜕𝑝𝜕𝑓𝜕𝑞𝜕𝑝𝜕𝜎𝜕𝑞𝜕𝜎^𝑇+𝜕𝑔𝜕𝑞𝜕𝑓𝜕𝑝𝜕𝑞𝜕𝜎𝜕𝑝𝜕𝜎^𝑇+𝜕𝑔𝜕𝑞𝜕𝑓𝜕𝑞𝜕𝑞𝜕𝜎𝜕𝑞𝜕𝜎^𝑇. (49) Comparing this with (33a), (33b), (34a), and (34b), it follows that1𝐴=𝐴_𝐻𝜕𝑔𝜕𝑝𝜕𝑓1𝜕𝑝,𝐵=𝐴_𝐻𝜕𝑔𝜕𝑝𝜕𝑓,1𝜕𝑞𝐶=𝐴_𝐻𝜕𝑔𝜕𝑞𝜕𝑓1𝜕𝑝,𝐷=𝐴_𝐻𝜕𝑔𝜕𝑞𝜕𝑓.𝜕𝑞 (50a) Equation (33a) and (33b) can be seen as a general formula for traditional constitutive models, and (49) is a special form of (33a) and (33b). For associated models, when 𝑓=𝑔, 1𝐴=𝐴_𝐻𝜕𝑓𝜕𝑝²1,𝐵=𝐶=𝐴_𝐻𝜕𝑓𝜕𝑝𝜕𝑓,1𝜕𝑞𝐷=𝐴_𝐻𝜕𝑓𝜕𝑞². (50b)

Most traditional constitutive models are to determine the relationship between 𝑓,  𝑔 and 𝑝, 𝑞, which can be used to calculate the four model parameters 𝐴,  𝐵,  𝐶, and 𝐷 indirectly. It can be seen from (50a) and (50b) that the traditional nonassociated models actually assume that 𝐴𝐷−𝐵𝐶=0. (51a) The associated models still need to satisfy (51a) and also require 𝐵=𝐶. (51b)Rewrite (27a) and (27b) into matrix form d𝜀^𝑝_𝑣d𝜀^𝑝=𝐴𝐵𝐶𝐷d𝑝d𝑞. (52)

Clearly (51a) requires the determinant rank of the coefficient matrix in (52) to be 1 or requires d𝜀^𝑝_𝑣 and d𝜀^𝑝 to be linearly correlated. Equation (51b) additionally requires the coefficient matrix to be symmetric.

The traditional elastoplastic model in stress space can be translated to strain space by the duality of (45a), (45b), (45c), (33a), and (33b). Substituting (50a) and (50b) into (42), (45a), (45b), and (45c) yields the expression in strain space𝐷^𝜀_𝑒𝑝=𝐷_𝑒−1𝛽×𝛼𝜕𝑔𝜕𝑝𝜕𝑓𝜕𝑝𝜕𝜀_𝑣𝜕𝜀𝜕𝜀_𝑣𝜕𝜀^𝑇+𝜕𝑔𝜕𝑝𝜕𝑓𝜕𝑞𝜕𝜀_𝑣𝜕𝜕𝜀𝜀𝜕𝜀^𝑇+𝜕𝑔𝜕𝑞𝜕𝑓𝜕𝑝𝜕𝜀_𝑣𝜕𝜕𝜀𝜀𝜕𝜀^𝑇+1𝛼𝜕𝑔𝜕𝑞𝜕𝑓𝜕𝜕𝑞𝜀𝜕𝜕𝜀𝜀𝜕𝜀^𝑇, (53a) where 𝛽=𝐴_𝐻|𝐴|=13𝐺_𝑒𝜕𝑔𝜕𝑝𝜕𝑓+1𝜕𝑝𝐾_𝑒𝜕𝑔𝜕𝑞𝜕𝑓+𝐴𝜕𝑞3𝐾^𝐻_𝑒𝐺_𝑒. (53b)

Hence, the stiffness matrix of (53a) and (53b) is obtained using 𝑓 and 𝑔 from the traditional constitutive models, and the transformation from stress space to strain space is thus achieved. Obviously, the transformation only changes the mathematical calculation method of the coefficient and elastoplastic matrices and has no influence on a particular model itself or the loading-unloading criterion of the model. Therefore, the transformation is applicable to all traditional elastoplastic models.

7.2. Modified Cambridge Model

In the modified Cambridge model, 𝑞𝑓=𝑔=𝑝+²𝑀²𝑝−𝑝₀𝑒^{(1+𝑒0)/(𝜆−𝜅)𝐻}=0, (54) and so 𝜕𝑓𝜂𝜕𝑝=1−²𝑀²,𝜕𝑓=𝜕𝑞2𝜂𝑀²,𝐴_𝐻=1+𝑒₀𝜂𝜆−𝜅1−²𝑀²𝑝₀𝑒^{(1+𝑒0)/(𝜆−𝜅)𝐻}, (55) where 𝜂=𝑞/𝑝,  𝐻 is the hardening parameter (=𝜀^𝑝_𝑣, for the modified Cambridge model), 𝑀 is the stress ratio at critical state, 𝑝₀ is the initial mean stress, 𝑒₀ is the initial void ratio, 𝜆 is the slope of the normal compression line (NCL), and 𝜅 is the slope of the unloading line. The elastoplastic compliance matrix in stress space is expressed as 𝐶^𝜎_𝑒𝑝=𝐶_𝑒+1𝐴_𝐻×𝜂1−²𝑀²²𝜕𝑝𝜕𝜎𝜕𝑝𝜕𝜎^𝑇+2𝜂𝑀²𝜂1−²𝑀²𝜕𝑝𝜕𝜎𝜕𝑞𝜕𝜎^𝑇+2𝜂𝑀²𝜂1−²𝑀²𝜕𝑞𝜕𝜎𝜕𝑝𝜕𝜎^𝑇+4𝜂²𝑀⁴𝜕𝑞𝜕𝜎𝜕𝑞𝜕𝜎^𝑇. (56) The elastoplastic stiffness matrix in strain space is𝐷^𝜀_𝑒𝑝=𝐷_𝑒−1𝛽×𝐾_𝑒3𝐺_𝑒𝜂1−²𝑀²²𝜕𝜀_𝑣𝜕𝜀𝜕𝜀_𝑣𝜕𝜀^𝑇+2𝜂𝑀²𝜂1−²𝑀²𝜕𝜀_𝑣𝜕𝜕𝜀𝜀𝜕𝜀^𝑇+2𝜂𝑀²𝜂1−²𝑀²𝜕𝜀_𝑣𝜕𝜕𝜀𝜀𝜕𝜀^𝑇+3𝐺_𝑒𝐾_𝑒4𝜂²𝑀⁴𝜕𝜀𝜕𝜕𝜀𝜀𝜕𝜀^𝑇, (57a) where 1𝛽=3𝐺_𝑒𝜂1−²𝑀²²+1𝐾_𝑒4𝜂²𝑀⁴+𝐴_𝐻3𝐺_𝑒𝐾_𝑒. (57b)

A new hardening parameter for the modified Cambridge model was proposed by Yao et al. [12] as𝑀𝐻=d𝐻=^4𝑓−𝜂⁴𝑀⁴−𝜂⁴d𝜀^𝑝_𝑣=1Ωd𝜀^𝑝_𝑣. (58a) in which 𝑀Ω=⁴−𝜂⁴𝑀⁴_𝑓−𝜂⁴, (58b)where 𝑀_𝑓 is the potential failure stress ratio.

Yao improved the modified Cambridge model by replacing 𝐻=𝜀^𝑝_𝑣 with (58a) and (58b), which changes 𝐴_𝐻 in the modified Cambridge model to (1/Ω)𝐴_𝐻. The improved constitutive model is a unified hardening model and is suitable for sandy soil, which actually replaces {d𝜀^𝑝} by Ω{d𝜀^𝑝} with (49), or the following expression with (52): d𝜀^𝑝_𝑣d𝜀^𝑝_𝑣𝑐=d𝜀^𝑝d𝜀^𝑝_𝑐𝑀=Ω=⁴−𝜂⁴𝑀⁴_𝑓−𝜂⁴, (59) where d𝜀^𝑝_𝑣|_𝑐,  d𝜀^𝑝|_𝑐 are the volumetric strain and shear strain, respectively, that are calculated by the modified Cambridge model.

This modification can be further improved. For instance, the volumetric strain and shear strain of triaxial testing are first calculated by the modified Cambridge model, and then the ratio of the volumetric strain and shear strain can be fitted according to the test results, that is d𝜀^𝑝_𝑣d𝜀^𝑝_𝑣𝑐,d=𝜉𝑝,𝑞𝜀^𝑝d𝜀^𝑝_𝑐=𝜁𝑝,𝑞. (60) Therefore, (52) is modified to d𝜀^𝑝_𝑣d𝜀^𝑝=𝜉𝐴𝜉𝐵𝜁𝐶𝜁𝐷d𝑝d𝑞. (61)

𝜉 and 𝜁 can be estimated by polynomial fitting or other fitting methods. Yao’s hardening model is obtained when 𝜉=𝜁=(𝑀⁴−𝜂⁴)/(𝑀⁴_𝑓−𝜂⁴).

The above-mentioned modification method can be extended to other elastoplastic models. The correction coefficients 𝜉 and 𝜁 can be fitted based on the triaxial testing results or other testing results. The obtained matrix can improve the calculation accuracy of existing elastoplastic models or be used to establish new modified models. Note that there is no physical mechanism involved in the framework, and thus the loading-unloading criterion and the evolution of internal state variables of the original model should still be employed in the modified one.

8. Application: A Simple Model

Theoretically, the constitutive model of soil would be established if the parameters 𝐴,  𝐵,  𝐶, and 𝐷 in (27a), (27b), (33a), and (33b) are obtained by experiments such as conventional triaxial test, isotropic compression test, and 𝑝=Const. test.

For example, the equations for the tangent modulus 𝐸_𝑡 and tangential Poisson ratio 𝜇_𝑡 are obtained by fitting the curve from triaxial testing, that is𝐸_𝑡=𝜕𝜎₁−𝜎₃𝜕𝜀₁,𝜇_𝑡=−𝜕𝜀₃𝜕𝜀₁=121−𝜕𝜀_𝑣𝜕𝜀₁. (62a)(62b)

𝐸_𝑡 and 𝜇_𝑡 can be curve fitted by a polynomial or by the formulae in the Duncan-Chang E-μ model [1]. However, two supplementary equations are needed as there are four unknown parameters in (27a) and (27b). Therefore, an isotropic compression test or 𝑝=Const. test should be conducted, or the assumptions 𝐴𝐷−𝐵𝐶=0 and 𝐵=𝐶 are made.

In conventional triaxial test, 𝜎₃=Const,  d𝑝=(1/3)d𝜎₁, d𝑞=d𝜎₁,  d𝜀_𝑣=d𝜀₁+d𝜀₂+d𝜀₃=(1−2𝜇_𝑡)d𝜀₁, and d𝜀=(2/3)(1+𝜇_𝑡)d𝜀₁. According to (62a), (62b), (27a), and (27b) and the assumption in (51a) and (51b), we have𝐾𝐴=^2𝑒𝑝𝜔𝐾,𝐵=𝐶=^𝑒𝑝𝐺_𝑒𝑝𝜔𝐺,𝐷=^2𝑒𝑝𝜔, (63a) where 𝐾_𝑒𝑝=1−2𝜇_𝑡𝐸_𝑡−1−2𝜇_𝑒𝐸_𝑒,𝐺_𝑒𝑝=21+𝜇_𝑡3𝐸_𝑡−21+𝜇_𝑒3𝐸_𝑒,𝜔=𝐺_𝑒𝑝+13𝐾_𝑒𝑝=1𝐸_𝑡−1𝐸_𝑒. (63b)𝐸_𝑒 is the elastic modulus, and the elastic Poisson ratio 𝜇_𝑒 is generally taken as 0.3 for soil. It is obvious that 𝜔>0 is always fulfilled in (63a) and (63b).

By substituting (63a) and (63b) into (33a), (33b), or (35), the elastoplastic compliance matrix in stress space or the stiffness matrix in strain space is obtained. For convenience, we call this model the multiple potential surface model (MPS model).

It should be noted that, contrary to the Duncan-Chang model (DC model), 𝐸_𝑡 and 𝜇_𝑡 in the MPS model are not limited by the generalized Hooke’s law; that is, the MPS model is still available when 𝜇_𝑡>0.5 and the stiffness matrix of the model is not singular. Actually, 𝐸_𝑡 and 𝜇_𝑡 in the new model are not the traditional modulus and Poisson ratio, but just the slope of the curves.

Figure 1 shows the calculation and test results for triaxial testing of Ottawa silica sand conducted by Wu [37]. The unit weight of the sand is 16.8 kN/m³ (=107 pcf). The test is a conventional consolidated-drained triaxial compression test (CD test). The confining pressures were 68.9, 206.7, and 344.5 kPa (= 10, 30, and 50 psi, resp.). During the CD test, confining pressure was firstly applied and then the specimen was consolidated. Deviation stress (𝜎₁−𝜎₃) was applied in the axial direction after consolidation. Variations of deviation stress and volumetric strain versus axial strain can be acquired in the test.

(a)

(b)

Figure 1 

Calculation and test results for triaxial test of Ottawa silica sand.

The calculations were made using the MPS model as well as the Duncan-Chang model, during which 𝐸_𝑡,  𝐸_𝑒, and 𝜇_𝑡 were calculated by the method proposed by Duncan and Chang [1], that is 𝐸_𝑡=𝑅1−^𝑓1−sin𝜙𝜎₁−𝜎₃2𝑐cos𝜙+2𝜎₃sin𝜙²𝐾𝑃_𝑎𝜎₃𝑃_𝑎^𝑛,𝐸_𝑒=𝐾_𝑢𝑟𝑃_𝑎𝜎₃𝑃_𝑎^𝑛, (64)𝜇_𝑡=𝐺−𝐹𝑙𝑔𝜎₃/𝑃_𝑎(1−𝐴^∗)²,𝐴^∗=𝐷𝜎₁−𝜎₃𝐾𝑃_𝑎𝜎₃/𝑃_𝑎^𝑛1−𝑅_𝑓1−sin𝜙𝜎₁−𝜎₃/2𝑐cos𝜙+2𝜎₃,sin𝜙 (65a)(65b)where 𝑐 is the cohesion of the soil, 𝜙 is the friction angle of the soil, 𝑃_𝑎 is the atmospheric pressure, 100 kPa; 𝐾,  𝐾_𝑢𝑟,  𝑛,  𝑅_𝑓,  𝐺, 𝐹, and 𝐷 are parameters.

The parameters in the calculation are taken as 𝑐=0 kPa, 𝜙=38^∘, 𝐾=1116,  𝐾_𝑢𝑟=1500,  𝑛=0.65,  𝑅_𝑓=0.88, 𝐹=0, and 𝐷=0, which are the same for both the DC and MPS models. The value of 𝐺 for the DC model is taken as 0.45 while for the MPS model it is 0.8, which is larger than 0.5. Although the calculation results for deviation stress are identical for the two models, the MPS model can reproduce the dilation of soil. Because of the limitation that 𝜇_𝑡<0.5 in the DC model, the dilatation of soil is not revealed and 𝜀_𝑣<0 is not achieved.

Figure 2 shows the calculation and test results of consolidated-drained triaxial compression test (CD test) of a rock-fill material from Hengshan Dam in China. The unit weight of the material is 20.7 kN/m³. The confining pressures were 300, 500, and 800 kPa, respectively.

(a)

(b)

Figure 2 

Calculation and test results for triaxial test of rockfill material.

The parameters for 𝐸_𝑡 are 𝑐=178 kPa, 𝜙=40.4°, 𝐾=1915, 𝐾_𝑢𝑟=2490,  𝑛=0.18, and 𝑅_𝑓=0.85, which are also the same for both the DC and MPS models. 𝜇_𝑡 in the DC model is still calculated using (65a) and (65b), and 𝐺=0.6,  𝐹=0.37, and 𝐷=0.023 in the calculation, while 𝜇_𝑡 in the MPS model is calculated using the method proposed by Shen and Zhang [38],𝜇_𝑡=12−𝑐_𝑑𝜎₃𝑃_𝑎^𝑛𝑑𝐸_𝑖𝑅_𝑓𝜎₁−𝜎_3𝑓1−𝑅_𝑑𝑅_𝑑×𝑅1−^𝑓𝑆_𝑙1−𝑅_𝑓𝑆_𝑙1−𝑅_𝑑𝑅_𝑑,𝐸_𝑖=𝐾𝑃_𝑎𝜎₃𝑃_𝑎^𝑛, (66a)(66b)in which 𝑆_𝑙=(𝜎₁−𝜎₃)/(𝜎₁−𝜎₃)_𝑓;(𝜎₁−𝜎₃)_𝑓 is the deviation stress at failure, (2𝑐cos𝜙+2𝜎₃sin𝜙)/(1−sin𝜙); 𝑐_𝑑,𝑛_𝑑, and 𝑅_𝑑 are parameters; 𝑐_𝑑=0.000224,  𝑛_𝑑=2.24, and 𝑅_𝑑=0.85 in the calculation.

Again, the calculation results for deviation stress are identical for the two models, while the MPS model reproduces the dilatation of soil giving better results than the DC model. Obviously, more appropriate results of volumetric strain can be acquired using the MPS model if we improve the calculation method of 𝜇_𝑡. However, this is impossible for the DC model due to the limitation that 𝜇_𝑡 cannot exceed 0.5 for the nonlinear elastic model.

9. Conclusions

The main tasks in establishing the constitutive equations for geomaterials are the determination of stress-strain relations in principal stress/strain space and the coordinate transformation of the relationship from principal stress/strain space to general coordinate space. The stress (or strain) and stress increment (or strain increment) in principal stress/strain space is expressed as a vector in a potential field in traditional elastic potential theory and plastic potential theory. However, the vector can be expressed more generally as the gradient vector of linearly independent potential functions. Based on this framework, the traditional models can be easily transformed from stress space to strain space and can be modified in a general way. This framework can also be used to establish new models based on curve fitting. Since it investigates constitutive models from a mathematical point of view independent of the material itself and relevant physical mechanism, the framework can be potentially used in a wider range, not limited to geomaterials. However, the lack of physical insights of materials and constitutive models may also hinder its development, for example, the constitutive model based on the framework may be oversimplified. The loading-unloading criteria as well as the evolutions of internal state variables are not considered in the framework. Therefore, the current framework is mainly useful for modifying the existing models. In the next research, Lode’s angle and noncoaxiality should also be investigated, and more test results are needed to make the verification.

Acknowledgments

The support of the Natural Science Foundation of China (51279085), National Basic Research Program of China (973 Program 2013CB036402), and the State Key Laboratory of Hydroscience and Engineering (sklhse-2012-D-01) is gratefully acknowledged.