# A Mathematical Approach to Establishing Constitutive Models for Geomaterials

Authors: Guang-hua Yang; Yu-xin Jie; Guang-xin Li
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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(1)dσij=Dijkldεkl.
Determining Dijkl is the major topic for constitutive models of geomaterials. Obviously, Dijkl 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 σi(i=1,2,3) and the principal strain εi(i=1,2,3), is obtained. To obtain Dijkl, 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(2a){dεip}3×1=[A]3×3{dσi}3×1,(2b)[A]3×3=[a11a12a13a21a22a23a31a32a33],where aij are the functions of total stress σi, total strain εi(i=1,2,3) or stress path.When the matrix rank of[A] is 1, or |A|=0, there exists a vector {α1α2α3}T and coefficients β1,β2,β3 to express [A] as
(3)[A]={β1β2β3}T{α1α2α3},
Therefore, substituting (3) into (2a) and (2b) gives(4a){dεip}={β1β2β3}T{α1α2α3}{dσi},
that is,
(4b){dεip}=dλ{β1β2β3}T,where
(5)dλ=∑i=13αidσi.
According to (4b),
(6)dε1p:dε2p:dε3p=β1:β2:β3.βi(i=1,2,3) is a function of σi or εi. When β={β1β2β3}T is of a field with potential, there is a potential function Q such that
(7)βi=∂Q∂σi.
Substituting (7) into (4b), we have
(8)dεip=dλ∂Q∂σi.
If we assume that dεip and σi have the same principal directions, the coordinate transformation can be expressed as follows:
(9)dεijp=dεip∂σi∂σij.
Substituting (9) into (8) gives
(10)dεijp=dλ∂Q∂σij.Similarly, from elastic potentials theory, there is a potential functionW in principal stress space, and εi is defined as
(11)εi=∂W∂σi.
If we assume that σi and εi have the same principal direction, the coordinate transformation can be expressed as follows:
(12)εij=εi∂σi∂σij.
Substituting (11) into (12),
(13)εij=∂W∂σij.In conclusion, traditional plastic potential theory corresponds to the case that the matrix[A] 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εip(i=1,2,3), are considered to be components of a vector dεp, 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.Whendεp is expressed in principal stress space, the coordinate orientations of dεp and σi are the same, and Φ1,Φ2,Φ3 are three linearly independent potential functions in principal stress space, and then the following expression is obtained:
(14)dεip=∑k=13λk∂Φk∂σi,
where λk(i=1,2,3) are the coefficients. Suppose that the principal directions of dεp and σi are the same, and we substitute (14) into (9), giving the tensor expression in general coordinate space as
(15)dεijp=∑k=13λk∂Φk∂σij.Similarly,dεp can also be expressed in strain space. Let Ψ1,Ψ2,Ψ3 be three linearly independent potential functions in principal strain space, and the following expression is obtained:
(16)dεijp=∑k=13μk∂Ψk∂εij.
Define the plastic stress increment as
(17){dσp}=[De]{dεp},
where [De] is the elastic stiffness matrix, and then the expressions in stress space and strain space are
(18)dσijp=∑k=13βk∂Gk∂σij,(19)dσijp=∑k=13αk∂Fk∂εij.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:(20)σij=∑k=13ηk∂Wk∂εij,εij=∑k=13χk∂Ωk∂σij,
where Wk,Ωk(i=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

IfAij and Eij are symmetric second-order tensors with the same principal directions, the following equations are obtained according to tensor theory and vector fitting:
(21)Aij=∑k=13λ-k∂IEk∂Eij,(22)Eij=∑k=13β-k∂IAk∂Aij,
where IEk(k=1,2,3) are three independent invariants of Eij, and IAk(k=1,2,3) are three independent invariants of Aij. For example, for the stress tensor σij, the three independent invariants can be Iσ1=σii, Iσ2=(1/2)σijσji,  Iσ3=(1/3)σikσknσnm;  Iσ1=p,  Iσ2=q,  Iσ3=θ; or Iσ1=σ1,  Iσ2=σ2,  Iσ3=σ3, where p is the mean stress, q is the deviatoric stress, θ is the Lode’s angle, and σ1,  σ2, and σ3 are the three principal stresses.IfAij=dεijp,  Eij=σij,  Iσk(k=1,2,3) are invariants of σij, it can be obtained from (21) that
(23)dεijp=∑k=13λ-k∂Iσk∂σij,
which is equivalent to (15) when Φk=Iσk.IfAij=dεijp, Eij=εij,  Iεk(k=1,2,3) are invariants of εij, then
(24)dεijp=∑k=13α-k∂Iεk∂εij.
If  Aij=εij,   Eij=σij, then
(25)σij=∑k=13β-k∂Iεk∂εij,εij=∑k=13λ-k∂Iσk∂σij.
Equation (25) corresponds to (20), respectively, when Wk=Iεk,Ωk=Iσk.Similarly, ifAij=Δεij,  Eij=Δσij, then according to (21), we have
(26)Δεij=∑k=13λ-k∂Iσk∂(Δσij).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 as(27a)dεvp=Adp+Bdq,(27b)dε-p=Cdp+Ddq,where A,   B,  C,   D are parameters.The following equation is obtained from the potential functionsΦ1=p,   Φ2=q in (15) without considering the effect of the Lode’s angle:
(28)dεijp=λ1∂p∂σij+λ2∂q∂σij.
In matrix form, this is
(29){dεp}=λ1{∂p∂σ}+λ2{∂q∂σ}.
It is obvious that
(30)λ1=dεvp,λ2=dε-p.
According to (27a), (27b), (29), and (30),
(31){dεp}=(Adp+Bdq){∂p∂σ}+(Cdp+Ddq){∂q∂σ}.
Since
(32)dp={∂p∂σ}T{dσ},dq={∂q∂σ}T{dσ},
it follows that(33a){dεp}=[Cpσ]{dσ},(33b)[Cpσ]=A{∂p∂σ}{∂p∂σ}T+B{∂p∂σ}{∂q∂σ}T+C{∂q∂σ}{∂p∂σ}T+D{∂q∂σ}{∂q∂σ}T,where [Cpσ] is the plastic compliance matrix.Since{dε}={dεe}+{dεp}, then(34a){dε}=[Cepσ]{dσ},
where
(34b)[Cepσ]=[Ce]+[Cpσ].[Ce] is the elastic compliance matrix and [Cepσ] is the elastoplastic compliance matrix.Therefore, the elastoplastic model in stress space can be established once the parametersA, B,  C, and D are determined. Note that A,  B,  C, and D 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σp}=[De]{dεp}, where [De] is the elastic matrix and {dεp} is the plastic strain increment.ε
v,  ε-, 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 εv and ε- as potential functions, that is, F1=εv,F2=ε- in (19), then
(35)dσijp=λ1∂εv∂εij+λ2∂ε-∂εij.
Written in matrix form,
(36){dσp}=λ1{∂εv∂ε}+λ2{∂ε-∂ε}.
Obviously,
(37)λ1=dpp,λ2=dqp.
From the definition of (17),(38a)dpp=Kedεvp=Ke(dεv-dεve)=Kedεv-dp,(38b)dqp=3Gedε-p=3Ge(dε--dε-e)=3Gedε--dq,where Ke,  Ge are the elastic bulk modulus and shear modulus.Considering(39a)dεve=1Kedp,(39b)dε-e=13Gedq,substituting (39a) and (39b) into (27a) and (27b) yields(40a)dεv=(1Ke+A)dp+Bdq,(40b)dε-=Cdp+(13Ge+D)dq.It can now be calculated that(41a)dp=A-dεv+B-dε-,(41b)dq=C-dεv+D-dε-,where
(42)A-=1|A|(D+13Ge),B-=-B|A|,C-=-C|A|,D-=1|A|(A+1Ke),|A|=3DGe+AKe+13KeGe+(AD-BC).
Substituting (41a) and (41b) into (37), (38a), and (38b) gives(43a)λ1=dpp=(Ke-A-)dεv-B-dε-,(43b)λ2=dqp=-C-dεv+(3Ge-D-)dε-.In addition,
(44)dεv={∂εv∂ε}T{dε},dε-={∂ε-∂ε}T{dε}.
Substituting (44) into (43a), (43b), and (36) yields(45a){dσp}=[Dpε]{dε},
where
(45b)[Dpε]=(Ke-A-){∂εv∂ε}{∂εv∂ε}T-B-{∂εv∂ε}{∂ε-∂ε}T-C-{∂ε-∂ε}{∂εv∂ε}T+(3Ge-D-){∂ε-∂ε}{∂ε-∂ε}T,
that is,
(45c)[Dpε]=1|A|{αA{∂εv∂ε}{∂εv∂ε}T+B{∂εv∂ε}{∂ε-∂ε}T=1|A|+C{∂ε-∂ε}{∂εv∂ε}T+Dα{∂ε-∂ε}{∂ε-∂ε}T},where α=Ke/3Ge=2(1+μe)/9(1-2μe), and μe is the elastic Poisson ratio.Hence, the total stress increment can be expressed as(46a){dσ}=[De]{dεe}=[De]({dε}-{dεp})=[De]{dε}-{dσp}=[De]{dε}-[Dpε]{dε}=[Depε]{dε},
where
(46b)[Depε]=[De]-[Dpε].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(47)[Cepσ]=[Ce]+1AH{∂g∂σ}{∂f∂σ}T,
where f and g are the yield function and plastic potential function, AH=-∂f/∂H{∂H/∂εp}T{∂g/∂σ} is the plastic hardening modulus, and H is the hardening parameter.f and g are usually expressed in terms of the stress invariants, p,  q,  θ. If the effect of the Lode’s angle θ is not considered, the expression only concerns p and q, that is,(48a){∂f∂σ}=∂f∂p{∂p∂σ}+∂f∂q{∂q∂σ},(48b){∂g∂σ}=∂g∂p{∂p∂σ}+∂g∂q{∂q∂σ}.Substituting (48a) and (48b) into (47), we have
(49)[Cepσ]=[Ce]+1AH×{∂g∂p∂f∂p{∂p∂σ}{∂p∂σ}T+∂g∂p∂f∂q{∂p∂σ}{∂q∂σ}T×+∂g∂q∂f∂p{∂q∂σ}{∂p∂σ}T+∂g∂q∂f∂q{∂q∂σ}{∂q∂σ}T}.
Comparing this with (33a), (33b), (34a), and (34b), it follows that(50a)A=1AH∂g∂p∂f∂p,B=1AH∂g∂p∂f∂q,C=1AH∂g∂q∂f∂p,D=1AH∂g∂q∂f∂q.
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 f=g,
(50b)A=1AH(∂f∂p)2,B=C=1AH∂f∂p∂f∂q,D=1AH(∂f∂q)2.Most traditional constitutive models are to determine the relationship betweenf,  g and p, q, which can be used to calculate the four model parameters A,  B,  C, and D indirectly. It can be seen from (50a) and (50b) that the traditional nonassociated models actually assume that (51a)AD-BC=0.
The associated models still need to satisfy (51a) and also require
(51b)B=C.Rewrite (27a) and (27b) into matrix form
(52){dεvpdε-p}=[ABCD]{dpdq}.Clearly (51a) requires the determinant rank of the coefficient matrix in (52) to be 1 or requires dεvp and dε-p 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(53a)[Depε]=[De]-1β×{α∂g∂p∂f∂p{∂εv∂ε}{∂εv∂ε}T+∂g∂p∂f∂q{∂εv∂ε}{∂ε-∂ε}T×+∂g∂q∂f∂p{∂εv∂ε}{∂ε-∂ε}T+1α∂g∂q∂f∂q{∂ε-∂ε}{∂ε-∂ε}T},
where
(53b)β=AH|A|=13Ge∂g∂p∂f∂p+1Ke∂g∂q∂f∂q+AH3KeGe.Hence, the stiffness matrix of (53a) and (53b) is obtained using f and g 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,(54)f=g=p+q2M2p-p0e(1+e0)/(λ-κ)H=0,
and so
(55)∂f∂p=1-η2M2,∂f∂q=2ηM2,AH=1+e0λ-κ(1-η2M2)p0e(1+e0)/(λ-κ)H,
where η=q/p,  H is the hardening parameter (=εvp, for the modified Cambridge model), M is the stress ratio at critical state, p0 is the initial mean stress, e0 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
(56)[Cepσ]=[Ce]+1AH×{(1-η2M2)2{∂p∂σ}{∂p∂σ}T×+2ηM2(1-η2M2){∂p∂σ}{∂q∂σ}T×+2ηM2(1-η2M2){∂q∂σ}{∂p∂σ}T×+4η2M4{∂q∂σ}{∂q∂σ}T{(1-η2M2)2{∂p∂σ}{∂p∂σ}T}.
The elastoplastic stiffness matrix in strain space is(57a)[Depε]=[De]-1β×{Ke3Ge(1-η2M2)2{∂εv∂ε}{∂εv∂ε}T×+2ηM2(1-η2M2){∂εv∂ε}{∂ε-∂ε}T×+2ηM2(1-η2M2){∂εv∂ε}{∂ε-∂ε}T×+3GeKe4η2M4{∂ε-∂ε}{∂ε-∂ε}TKe3Ge(1-η2M2)2},
where
(57b)β=13Ge(1-η2M2)2+1Ke4η2M4+AH3GeKe.A new hardening parameter for the modified Cambridge model was proposed by Yao et al. [12] as(58a)H=∫dH=∫Mf4-η4M4-η4dεvp=∫1Ωdεvp.
in which
(58b)Ω=M4-η4Mf4-η4,where Mf is the potential failure stress ratio.Yao improved the modified Cambridge model by replacingH=εvp with (58a) and (58b), which changes AH in the modified Cambridge model to (1/Ω)AH. The improved constitutive model is a unified hardening model and is suitable for sandy soil, which actually replaces {dεp} by Ω{dεp} with (49), or the following expression with (52):
(59)dεvpdεvp|c=dε-pdε-p|c=Ω=M4-η4Mf4-η4,
where dεvp|c,  dε-p|c 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(60)dεvpdεvp|c=ξ(p,q),dε-pdε-p|c=ζ(p,q).
Therefore, (52) is modified to
(61){dεvpdε-p}=[ξAξBζCζD]{dpdq}.ξ and ζ can be estimated by polynomial fitting or other fitting methods. Yao’s hardening model is obtained when ξ=ζ=(M4-η4)/(Mf4-η4).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.


## 7.1. General Form

The elastoplastic compliance matrix of the traditional elastoplastic model is(47)[Cepσ]=[Ce]+1AH{∂g∂σ}{∂f∂σ}T,
where f and g are the yield function and plastic potential function, AH=-∂f/∂H{∂H/∂εp}T{∂g/∂σ} is the plastic hardening modulus, and H is the hardening parameter.f and g are usually expressed in terms of the stress invariants, p,  q,  θ. If the effect of the Lode’s angle θ is not considered, the expression only concerns p and q, that is,(48a){∂f∂σ}=∂f∂p{∂p∂σ}+∂f∂q{∂q∂σ},(48b){∂g∂σ}=∂g∂p{∂p∂σ}+∂g∂q{∂q∂σ}.Substituting (48a) and (48b) into (47), we have
(49)[Cepσ]=[Ce]+1AH×{∂g∂p∂f∂p{∂p∂σ}{∂p∂σ}T+∂g∂p∂f∂q{∂p∂σ}{∂q∂σ}T×+∂g∂q∂f∂p{∂q∂σ}{∂p∂σ}T+∂g∂q∂f∂q{∂q∂σ}{∂q∂σ}T}.
Comparing this with (33a), (33b), (34a), and (34b), it follows that(50a)A=1AH∂g∂p∂f∂p,B=1AH∂g∂p∂f∂q,C=1AH∂g∂q∂f∂p,D=1AH∂g∂q∂f∂q.
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 f=g,
(50b)A=1AH(∂f∂p)2,B=C=1AH∂f∂p∂f∂q,D=1AH(∂f∂q)2.Most traditional constitutive models are to determine the relationship betweenf,  g and p, q, which can be used to calculate the four model parameters A,  B,  C, and D indirectly. It can be seen from (50a) and (50b) that the traditional nonassociated models actually assume that (51a)AD-BC=0.
The associated models still need to satisfy (51a) and also require
(51b)B=C.Rewrite (27a) and (27b) into matrix form
(52){dεvpdε-p}=[ABCD]{dpdq}.Clearly (51a) requires the determinant rank of the coefficient matrix in (52) to be 1 or requires dεvp and dε-p 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(53a)[Depε]=[De]-1β×{α∂g∂p∂f∂p{∂εv∂ε}{∂εv∂ε}T+∂g∂p∂f∂q{∂εv∂ε}{∂ε-∂ε}T×+∂g∂q∂f∂p{∂εv∂ε}{∂ε-∂ε}T+1α∂g∂q∂f∂q{∂ε-∂ε}{∂ε-∂ε}T},
where
(53b)β=AH|A|=13Ge∂g∂p∂f∂p+1Ke∂g∂q∂f∂q+AH3KeGe.Hence, the stiffness matrix of (53a) and (53b) is obtained using f and g 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,(54)f=g=p+q2M2p-p0e(1+e0)/(λ-κ)H=0,
and so
(55)∂f∂p=1-η2M2,∂f∂q=2ηM2,AH=1+e0λ-κ(1-η2M2)p0e(1+e0)/(λ-κ)H,
where η=q/p,  H is the hardening parameter (=εvp, for the modified Cambridge model), M is the stress ratio at critical state, p0 is the initial mean stress, e0 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
(56)[Cepσ]=[Ce]+1AH×{(1-η2M2)2{∂p∂σ}{∂p∂σ}T×+2ηM2(1-η2M2){∂p∂σ}{∂q∂σ}T×+2ηM2(1-η2M2){∂q∂σ}{∂p∂σ}T×+4η2M4{∂q∂σ}{∂q∂σ}T{(1-η2M2)2{∂p∂σ}{∂p∂σ}T}.
The elastoplastic stiffness matrix in strain space is(57a)[Depε]=[De]-1β×{Ke3Ge(1-η2M2)2{∂εv∂ε}{∂εv∂ε}T×+2ηM2(1-η2M2){∂εv∂ε}{∂ε-∂ε}T×+2ηM2(1-η2M2){∂εv∂ε}{∂ε-∂ε}T×+3GeKe4η2M4{∂ε-∂ε}{∂ε-∂ε}TKe3Ge(1-η2M2)2},
where
(57b)β=13Ge(1-η2M2)2+1Ke4η2M4+AH3GeKe.A new hardening parameter for the modified Cambridge model was proposed by Yao et al. [12] as(58a)H=∫dH=∫Mf4-η4M4-η4dεvp=∫1Ωdεvp.
in which
(58b)Ω=M4-η4Mf4-η4,where Mf is the potential failure stress ratio.Yao improved the modified Cambridge model by replacingH=εvp with (58a) and (58b), which changes AH in the modified Cambridge model to (1/Ω)AH. The improved constitutive model is a unified hardening model and is suitable for sandy soil, which actually replaces {dεp} by Ω{dεp} with (49), or the following expression with (52):
(59)dεvpdεvp|c=dε-pdε-p|c=Ω=M4-η4Mf4-η4,
where dεvp|c,  dε-p|c 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(60)dεvpdεvp|c=ξ(p,q),dε-pdε-p|c=ζ(p,q).
Therefore, (52) is modified to
(61){dεvpdε-p}=[ξAξBζCζD]{dpdq}.ξ and ζ can be estimated by polynomial fitting or other fitting methods. Yao’s hardening model is obtained when ξ=ζ=(M4-η4)/(Mf4-η4).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 parametersA,  B,  C, and D in (27a), (27b), (33a), and (33b) are obtained by experiments such as conventional triaxial test, isotropic compression test, and p=Const. test.For example, the equations for the tangent modulusEt and tangential Poisson ratio μt are obtained by fitting the curve from triaxial testing, that is(62a)Et=∂(σ1-σ3)∂ε1,(62b)μt=-∂ε3∂ε1=12(1-∂εv∂ε1).E
t and μt 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 p=Const. test should be conducted, or the assumptions AD-BC=0 and B=C are made.In conventional triaxial test,σ3=Const,  dp=(1/3)dσ1, dq=dσ1,  dεv=dε1+dε2+dε3=(1-2μt)dε1, and dε-=(2/3)(1+μt)dε1. According to (62a), (62b), (27a), and (27b) and the assumption in (51a) and (51b), we have(63a)A=Kep2ω,B=C=KepGepω,D=Gep2ω,
where
(63b)Kep=1-2μtEt-1-2μeEe,Gep=2(1+μt)3Et-2(1+μe)3Ee,ω=Gep+13Kep=1Et-1Ee.Ee is the elastic modulus, and the elastic Poisson ratio μe 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),Et and μt in the MPS model are not limited by the generalized Hooke’s law; that is, the MPS model is still available when μt>0.5 and the stiffness matrix of the model is not singular. Actually, Et and μt in the new model are not the traditional modulus and Poisson ratio, but just the slope of the curves.Figure1 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/m3 (=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 (σ1-σ3) was applied in the axial direction after consolidation. Variations of deviation stress and volumetric strain versus axial strain can be acquired in the test.Calculation and test results for triaxial test of Ottawa silica sand.
(a)
(b)The calculations were made using the MPS model as well as the Duncan-Chang model, during whichEt,  Ee, and μt were calculated by the method proposed by Duncan and Chang [1], that is
(64)Et=[1-Rf(1-sinϕ)(σ1-σ3)2ccosϕ+2σ3sinϕ]2KPa(σ3Pa)n,Ee=KurPa(σ3Pa)n,(65a)μt=G-Flg(σ3/Pa)(1-A*)2,(65b)A*=D(σ1-σ3)KPa(σ3/Pa)n[1-Rf(1-sinϕ)(σ1-σ3)/(2ccosϕ+2σ3sinϕ)],where c is the cohesion of the soil, ϕ is the friction angle of the soil, Pa is the atmospheric pressure, 100 kPa; K,  Kur,  n,  Rf,  G, F, and D are parameters.The parameters in the calculation are taken asc=0 kPa, ϕ=38°, K=1116,  Kur=1500,  n=0.65,  Rf=0.88, F=0, and D=0, which are the same for both the DC and MPS models. The value of G 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 μt<0.5 in the DC model, the dilatation of soil is not revealed and εv<0 is not achieved.Figure2 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/m3. The confining pressures were 300, 500, and 800 kPa, respectively.Calculation and test results for triaxial test of rockfill material.
(a)
(b)The parameters forEt are c=178 kPa, ϕ=40.4°, K=1915, Kur=2490,  n=0.18, and Rf=0.85, which are also the same for both the DC and MPS models. μt in the DC model is still calculated using (65a) and (65b), and G=0.6,  F=0.37, and D=0.023 in the calculation, while μt in the MPS model is calculated using the method proposed by Shen and Zhang [38],(66a)μt=12-cd(σ3Pa)ndEiRf(σ1-σ3)f1-RdRd×(1-RfSl1-RfSl1-RdRd),(66b)Ei=KPa(σ3Pa)n,in which Sl=(σ1-σ3)/(σ1-σ3)f;(σ1-σ3)f is the deviation stress at failure, (2ccosϕ+2σ3sinϕ)/(1-sinϕ); cd,nd, and Rd are parameters; cd=0.000224,  nd=2.24, and Rd=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μt. However, this is impossible for the DC model due to the limitation that μt 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.



---

*Source: 739068-2013-06-22.xml*