- Email: [email protected]
A hierarchical single yield surface for frictional materials
Computers and Geotechnics 36 (2009) 960–967
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
A hierarchical single yield surface for frictional materials G. Mortara * Department of Mechanics and Materials, University of Reggio Calabria, Via Graziella, Feo di Vito, 89122 Reggio Calabria, Italy
a r t i c l e
i n f o
Article history: Received 23 January 2009 Received in revised form 2 March 2009 Accepted 4 March 2009 Available online 8 April 2009 Keywords: Constitutive modelling Three-dimensional stress Yield surface Shear strength Frictional materials
a b s t r a c t This paper extends the formulation of a recent yield and failure criterion [Mortara G. A new yield and failure criterion for geomaterials. Géotechnique 2008;58(2):125–32] by modifying its shape in the deviatoric representation. Such a modification improves the performance of the surface introducing a constitutive parameter associated to the variation of the Lode angle. The criterion is validated through the comparison with experimental data obtained from multiaxial tests on soils, rocks and concrete. A special form of the criterion will be also derived for numerical applications in order to remove singularities of the previous criterion. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction Elastoplastic models of geomaterials need the definition of yield and potential surfaces under three-dimensional stress conditions. The Mohr–Coulomb failure condition is often used in geomechanics even though it has corners that complicate the implementation in numerical codes. For this reason, over the years, smooth curves close to the Mohr–Coulomb criterion have been formulated [14,10,22]. Recently Mortara [18] formulated a new yield and failure criterion that generalises the criteria by Matsuoka and Nakai [14] and Lade and Duncan [10] and can include as well the criteria by Lade [8], Kim and Lade [6] and Houlsby [4]. The criterion can take into account tensile strength and different shapes in the meridian sections ranging from linear and curvilinear cones to closed domains. All surfaces have the characteristics to change their shape in the deviatoric plane for changing values of the friction angle /0 . In particular the continuous transition from circular to triangular shapes is observed for increasing values of / 0 as for the well known criteria by Matsuoka and Nakai [14] and Lade and Duncan [10]. The update version shown in this paper adds the possibility to modify the deviatoric shape of the criterion through the introduction of a correcting function dependent on Lode’s angle. Some changes are made in the expression of the previous criterion in order to avoid some drawbacks. Later, some comparisons will be presented between the criterion and experimental data based on triaxial and biaxial test on soils, rocks and concrete.
* Tel.: +39 0965 875 271; fax: +39 0965 875 201. E-mail address: [email protected] 0266-352X/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2009.03.007
2. General form of the criterion The criterion assumes the following general form:
f ¼ In1 Im m¼ 2 nI 3 ¼ 0
3n 2
ð1Þ
where n is the parameter determining the continuous transition of the surface shape in the deviatoric plane. The Matsuoka and Nakai [14] failure surface is obtained for n = 1 while for n = 3 criterion by Lade and Duncan [10] is obtained. I1 , I2 and I3 are the first, the second and the third stress invariants of the translated principal i ¼ ri þ t being t the absolute value of the isotropic tensile stresses r strength:
1 þ r 2 þ r 3 I1 ¼ r
1r 2r 3r 2 þ r 3 þ r 1 I2 ¼ r
1r 2r 3 I3 ¼ r
ð2Þ
Finally, n is the hierarchical function ruling both the size and shape of the criterion. Its general form is given by:
n ¼ nh þ an1 n2 n3
a ¼ n0 nh
ð3Þ
where n0 is the basic size and shape function dependent on the triaxial compression friction angle /00 and the parameter n; constant nh is the value that In1 Im 2 =I 3 takes in isotropic conditions; n1 and n2 are two correcting functions generating different shapes in the deviatoric planes; n3 is a correcting function introduced to provide the possibility of new shapes in deviatoric sections. The definition of the above mentioned functions is listed below.
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
n0 ¼
ð3 sin /00 Þn ð3 þ sin /00 Þm
ð4Þ
ð1 þ sin /00 Þð1 sin /00 Þ2m
nh ¼ 33m c1 I0 n1 ¼ I1 n2 ¼ 1 sgnðc2 Þ
ð5Þ ð6Þ I1 I0
!c 2 ð7Þ
3 n3 ¼ 1 þ e 1 cos h 2
ð8Þ
When c2 = 0 the geometrical meaning of constant I0 is the following: the curved surface obtained for c1 > 0 and the straight surface obtained for c1 = 0 intersect for I1 ¼ I0 . When c2 > 0, instead, constant I0 represents the value of I1 > 0 for which the surface crosses the hydrostatic axis. The Lode’s angle h in Eq. (8) is such that:
cos 3h ¼
pffiffiffiffiffiffi 27 J 3 2 J3
ð9Þ
where the deviatoric invariants J and J3 are defined as follows:
J¼
pffiffiffiffiffiffiffiffi 1 J2 J 2 ¼ I2 I21 3
1 2 3 J 3 ¼ I3 I2 I1 þ I 3 27 1
ð10Þ
Definition (9) provides h = 0° for triaxial compression and h = 60° for triaxial extension. In its general form the definition of the criterion needs seven parameters, namely /00 (or a), n, t, I0, c1, c2 and e. However, the basic form of the criterion, which is represented by a cone pinned at the origin of the principal stress space, requires only the two parameters /00 and n. The others can be introduced hierarchically depending on the experimental data that have to be represented. Introduction of the curvature in the meridian sections is obtained through parameters c1 and I0. Closed surfaces are instead obtained by choosing a non zero value for parameter c2. The geometrical meaning of the constants appearing in the yield function expression is reported in Mortara [18]. In order to show the different surfaces that can be obtained for different parameters, in Fig. 1 the three-dimensional representation of two surfaces is reported where the effect of taking c1 > 0 and c2 = 0 (Fig. 1a), and c2 > 0 (Fig. 1b) is clearly shown.
m 1 2 1 1 3 f ¼ In1 J 2 þ I21 n pffiffiffiffiffiffi cos 3hJ 3 I1 J 2 þ I1 ¼ 0 3 3 27 27
961
ð11Þ
Eq. (11) permits to observe that a third degree equation is obtained both for the Matsuoka and Nakai [14] and the Lade and Duncan [10] criteria. A third degree equation provides one or three real roots and for n = 1 or n = 3 Eq. (11) results in one negative and two positive values of J. In particular, for the Matsuoka and Nakai criterion as for the Lade and Duncan criterion it is possible to find stress states for which f < 0 even for points outside the current elastic domain [18]. This problem, particularly important for numerical applications, is common to different yield criteria [3]. For any non integer value pffiffiffi n, instead, Eq. (11) is not defined for J > J u ¼ I1 = 3, namely outside a circular cone in the pffiffiffi general stress space pffiffiffi having origin for I1 ¼ 0 and radius J u ¼ I1 = 3. Condition J I1 = 3 prevents the possibility to find more than a positive solution to f = 0 for stress states outside the current elastic domain. However taking non integer values for n does not avoid the possibility to find stress states such that J > Ju. In order to avoid discontinuities for the yield surface some comments should be made on the singularities of Eq. (1). First of all, if 1 > 0 and r 2 ¼ r 3 ¼ 0 it results f = 0 we apply Eq. (1) for a > 0, r independently of the value of / 0 given that this stress state corresponds to I1 > 0, I2 ¼ I3 ¼ 0. Nevertheless, this stress state should correspond to / 0 = 90°. If we limit ourselves to study the behaviour of frictional materials for /00 /00L being /00L a value close to 90° (e.g. /00L ¼ 89 ) we can remove pitfalls of Eq. (1), or equivalently of Eq. (11), for all values of n. For this purpose, Eq. (11) is rewritten in the form:
f ¼ f1 f2 nf3 þ I1N þ I2N ¼ 0
ð12Þ
where
n f1 ¼ hI1 i þ 1 sgnhI1 i 1 þ sgnhI1 i m 1 f2 ¼ hJ 2 þ I21 i þ 1 sgnhI2 i 3 f3 ¼
J3 J 3 þ 1 sgnJ
3. Special form of the criterion
I1N ¼
I1 þ jI1 j 2
Criterion (1) can be also expressed in terms of invariants I1 , J and h as follows
I2N ¼
I2 þ jI2 j 2
1 1 J 3 hI1 iJ 2 þ hI1 i3 3 27
Fig. 1. Effect of parameter c1 (a) and parameter c2 (b) in the principal stress space.
ð13Þ
962
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
The Macaulay brackets hi applied to any scalar x gives hxi ¼ ðx þ jxÞ=2. Variable J* in Eq. (13) is given by
J ¼ J L hJ L Ji
ð14Þ
Variable JL is the limit value for invariant J related to the chosen value of /00L that clearly depends on the Lode’s angle. For the sake of simplicity its value is found analysing the Lode dependence of the Mohr–Coulomb criterion for /00 ¼ /00L :
b 2 sin /00L 3 sin /00L J L ¼ pffiffiffi b¼ 0 hI 1 i 3 þ sin /00L 3b cos h þ ð2 bÞ sin h 3 sin /0L
/0h ¼ arcsin ð15Þ
Constant b is the aspect ratio, namely the ratio between the value of invariant J for h = 60° and h = 0° [2]. Variables I1N and I2N , instead, represent the value of f for negative values of I1 and I2 , respectively. For such values, in fact, I1N and I2N are equal to jI1 j and jI2 j, respectively and they are the only non null terms in Eq. (12) while for positive values of I1 and I2 they do not influence the value of f. The new form is completed by new expressions of functions n1 and n2. As can be noted analysing their expression, function (6) is discontinuous when I1 ¼ 0 like function (7) when I1 ¼ 0 and c2 = 0. For this reason the expressions of n1 and n2 are modified as follows:
n1 ¼
I0 I1L þ hI1 I1L i
!c 1 n2 ¼ 1 sgnðc2 Þ
I1L þ hI1 I1L i I0
!c 2
ð16Þ where I1L is an arbitrary small value of invariant I1 . In the updated form (12) no stress states can be found giving f < 0 outside the current elastic domain and the function is defined for all values of I1 , J, h and n. In order to understand the importance of form (12), Table 1 reports some example of stress states for which Eq. (11) fails contrarily to Eq. (12). 4. New shapes in the deviatoric sections: effect of function n3 Function n3 provides the possibility to modify the surface shape in the deviatoric sections. The simple idea to improve the capability of the new criterion in the deviatoric plane is to generate a surface with shape varying with the Lode’s angle h. For this purpose the correcting function n3, which is a trigonometric function of period 2p/3, is used. The role of parameters n and e is similar, given that their effect is to enlarge the deviatoric section with continuity from values of h ranging from 0° (triaxial compression) to 60° (triaxial extension). However, as shown in Fig. 2, they provide different shapes in the deviatoric section. Fig. 2a reports five functions drawn for /00 ¼ 30 , t = 100, e = 0 and n = 2, 4, 6, 8, 10 while Fig. 2b reports five functions drawn for /00 ¼ 30 , t = 100, n = 2 and e = 0, 0.2, 0.4, 0.6, 0.8. As it can be observed the surfaces for increasing value of n tend to become circular while those obtained for increasing value of e tend to assume a sharper profile in the deviatoric plane as the Lode’s angle tends to approach triaxial extension conditions. Fig. 3 shows the different shapes produced
Table 1 Values of Eqs. (11) and (12) for a > 0, t = c1 = c2 = e = 0 and for some values of n, r1, r2 and r3. n
r1
r2
r3
Eq. (11)
Eq. (12)
60 or >3 decimal <3 P3 Odd integer All but odd integer
0 <0 >0 >0 >0 >0
0
0
Undefined Undefined 0 Undefined <0 Undefined
0 I1 >0 >0 I2 I2
r1
r1
0 0 r1/2 r1/2
0 0 r1/2 r1/2
by parameters n and e in the biaxial plane. Parameters are the same as those of Fig. 2 and the figures represent the sections of the 3D surfaces with the plane r3 = 0. It can be still observed the different effects that parameters n and e produce on the biaxial plane especially as the Lode’s angle approaches 60°. Variation of parameters n and e influences the value of the friction angle /0h , which is related to the current state of principal stress through the relation:
r1 r3 r1 þ r3
ð17Þ
Fig. 4 shows the variation of /0h with the Lode’s angle h of seven curves: six are obtained for /00 ¼ 30 and different values of n and e while curve C0 refers to the Mohr–Coulomb criterion that, as well known, assumes the constant value /0h ¼ /00 . The maximum difference /0h /00 among the curves related to Eq. (1) and the Mohr–Coulomb curve is not obtained for the same value of the Lode’s angle: for e = 0, and irrespectively on the value /00 , such difference is attained for h<30° if n<3 and for h > 30° if n > 3. For the Lade and Duncan criterion, instead, the maximum difference /0h /00 is always obtained for h = 30°. The effect of parameter e on the shape in the deviatoric plane can be discussed by analysing curves C1, C3 and C4. Curve C1 is related to the Matsuoka–Nakai criterion while curve C3 is related to n = 1 and e = 0.5. The difference in terms of friction angle of these curves is evident and their maximum difference is attained for h = 60°. Curve C4 is related to e = 0 and n = 4.36 that is the value of n that assures the same value of shear strength for triaxial extension conditions of curve C3. It is evident how, although this two curves attain the same value of /0h for h = 60°, they produce different shapes in the deviatoric plane. Similar considerations follow for curves C2, C5 and C6. 5. Derivation of parameters and comparisons with experimental data In this section the proposed criterion is compared with experimental data related to triaxial and biaxial tests on soils, rocks and concrete. The values of parameters used to match the experimental curves are reported in Table 2. Concerning the derivation of parameters, when the experimental data show that a closed surface is necessary for their representation, it is necessary to specify a value for parameter c2. To derive parameters of the yield surface from triaxial compression or extension tests one can observe that Eq. (1) can be rewritten in the form
fg ¼ F
Rcp1 1 Rcp2
a¼0
ð18Þ
where
27 F ¼ 2 cos 3h 27
1 3
m g2 27
g3 13 g2 þ 1
nh
g¼
q p0 þ t
Rp ¼
p0 þ t p0 þ t
ð19Þ
For the stress state where the deviatoric stress is maximum (q = qg) condition ofg/og = 0 must be verified and this allows us to derive parameter c1
c1 ¼
Rcp2 gg @F c2 F @g 1 Rcp2
ð20Þ
Once the coordinates of the point in the q–p0 plane are known, it is possible to solve Eq. (18) that provides a value of a for a chosen value of c2. Fig. 5 shows the effect of increasing values of parameter c2 on the shape of the curves in the q–p0 plane for n = 1, t=p00 ¼ 0:1, qg =p00 ¼ 0:6 and p0g =p0 ¼ 0:7 (gg = 0.75). The inner curve is related to the value c2 = 2.649 that is the minimum value of c2 assuring
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
963
Fig. 2. Different shapes in the deviatoric sections: (a) increasing values of n for e = 0; and (b) increasing value of e for n = 2.
Fig. 3. Shapes of the criterion in the biaxial plane: (a) effect of parameter n; and (b) effect of parameter e.
c1 > 0. The set of shapes achievable in meridian planes by varying the values of parameters c1 and c2 makes the model similar to existing yield surfaces and plastic potentials in the literature [20,19,11,12,1,3]. An advantage of the proposed model resides in the opportunity to avoid the use of tension cut-off that instead is necessary in models like for example the modified Cam Clay [21]. On the other hand a possible drawback of the model is that, for a given set of parameters, the value of the mean stress p0 ¼ pg for which the deviatoric stress is maximum is not constant but increases with h. Such change increases for increasing values of gg. The procedure summarised by Eqs. (18) and (20) is applied to derive the parameters related to triaxial compression tests on carbonate sand reported by [5]. In Fig. 6 two sets of data are interpreted, the first one related to uncemented carbonate sand and the second related to carbonate sand with C = 20% of cementation. It is observed that the model curves provide a good interpretation of the experimental data. The results shown in Figs. 5 and 6 highlight the flexibility of the proposed model in the meridian representation. The derivation of parameter e is simple once the parameters related to the meridian representation of the surface are known. Con-
stant e can be easily found from Eq. (11) by providing the value of Je namely the value that invariant J assumes for h = 60°.
0 1
m 2 n 1 2 1 B I1 J e þ 3 I1 C e¼ n A1 @ an1 n2 p2ffiffiffiffi J3e 13 I1 J2e þ 271 I31 h 27
ð21Þ
In Fig. 7 the comparison among experimental true triaxial test data on dense Santa Monica Beach sand [9] and the criteria by Lade [8], Matsuoka and Nakai [14] and Eq. (1) with e > 0 is shown. The comparison is done both in terms of octahedral shear stresses (Fig. 7a) and Mohr–Coulomb friction angles (Fig. 7b). As it can be observed, Lade’s model overestimates the octaedral stresses, and then the friction angles, in the midrange of h values. Lade [9] imputed such overestimation to shear banding occurring in the sand specimens. Matsuoka and Nakai model [14] matches quite well the experimental data for h < 30° but underestimates the remaining stress points. The value of e for the new model curve is obtained from (21) for Je = 118 kPa. As one can observe the scatter between new model and experimental points is limited for all values of the Lode’s angle. Fig. 8 shows the experimental data on dolomite performed by Mogi [15] are compared with the new criterion and the models
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
f
964
Fig. 4. Influence of parameters n and e on the variation of the friction angle with the Lode’s angle for /00 ¼ 30 .
Table 2 Values of the parameters for the new criterion reported in Figs. 6–10. Figures
Experimental data
n
t (kPa)
a
c1
c2
I0 (kPa)
e
6 6 7 8 9a 9b 10a 10b
Islam et al. [5] (C = 0%) Islam et al. [5] (C = 20%) Lade [9] Mogi [15] Mogi [16] Mogi [16] Lee et al. [13] Lee et al. [13]
1.00 1.00 0.90 2.90 4.10 3.90 2.70 2.70
0:00 p0e 0:05 p0e 0.00 7150 9000 6700 3000 4000
10.00 13.98 8.25 9.0 108 1.25 105 5.70 104 1.10 105 1.40 105
0.795 0.563 0.144 1.800 2.000 2.000 1.150 1.150
0.100 0.500 0.000 0.000 0.000 0.000 0.000 0.000
1:00 p0e 1:00 p0e 100 100 100 100 100 100
0.000 0.000 0.665 0.800 0.900 0.900 0.450 0.450
Fig. 5. Effect of parameter c2 on the shape of the criterion in the q–p0 plane.
by Kim and Lade [6] and Mortara [18]. The new model, thanks to the introduction of parameter e, matches very well the data in the Rendulic plane (Fig. 8a) while the other two models are less
Fig. 6. Derivation of parameters from experimental data in the q–p0 plane.
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
965
Fig. 7. Comparison among true triaxial test data on dense Santa Monica Beach sand [9] and three failure criteria: (a) deviatoric plane; and (b) friction angle versus Lode’s angle.
Fig. 8. Comparison among three models and experimental data on dolomite performed by Mogi [15]: (a) Rendulic plane; and (b) biaxial plane.
satisfactory in extension conditions. The positive effect of function n3 is evident also in Fig. 8b that reports the experimental data and model curves in biaxial stress conditions. The curve related to the previous model [18] was obtained for n = 3.9 and improved the simulation by Kim and Lade (n = 3) but underestimating the uniaxial compressive strength. The new model curve, obtained for n = 2.9 and e = 0.8, improves the agreement with the experimental data for all values of the ratio r1/r2 providing at the same time a better estimation of the uniaxial compressive strength. Fig. 9 reports the comparison between experimental data by Mogi [16,17] and model. In particular, Fig. 9a shows the comparison between biaxial tests on Solenhofen limestone for three constant values of r3, while in Fig. 9b the comparison deals with tests on Mizuno trachyte for four constant values of r3. In both
cases the comparison is fairly good except for the model curve for r3 = 100 MPa in Fig. 9b that underestimates the values of experimental data. Finally, Fig. 10, deals with the comparison of the model with two sets of biaxial data on concrete performed for different biaxial stress ratios r1/r2 reported by Lee et al. [13]. The tests deal with two different values of the uniaxial compressive strength namely rc = 30.3 MPa and rc = 39.0 MPa reported in Fig. 10a and b, respectively. The experimental data show trends already observed in previous works like the one by Kupfer et al. [7]. These trends highlight that for stress ratio conditions ranging from uniaxial compression and uniaxial extension (0° 6 h 6 60°) the experimental envelopes show concavity. However, the comparison between data and model is satisfactory for both sets of data. It should be noted that only
966
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
Fig. 9. Comparison between new model and experimental data in the biaxial plane: (a) Solenhofen limestone [16]; and (b) Mizuno trachyte [17].
Fig. 10. Comparison between new model and experimental data on concrete by Lee et al. [13].
parameters t and a have been changed in the model curves reported in Fig. 10. 6. Conclusions This paper has dealt with the formulation of a new isotropic yield and failure criterion for frictional materials in terms of the three stress invariants. Starting from the recent model by Mortara [18], which formulation generalises well known criteria, a more general formulation has been derived by adding a function that increases the performance of the criterion in the deviatoric represen-
tation. The introduction of the new function that modifies the deviatoric section of the criterion has revealed an increasing capability to model the experimental data in both the deviatoric and the biaxial plane. The paper has reported the procedure to derive the constitutive parameters both in the deviatoric and meridian planes. The criterion has been validated using experimental data related to biaxial and triaxial tests on soils, rocks and concrete, showing an overall good attitude in modelling the behaviour of frictional materials. A special form of the criterion has been derived in order to remove all discontinuities observed in the previous formulation [18] mak-
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
ing the new form suitable for the implementation in numerical codes. References [1] Aubertin M, Li L. A porosity-dependent inelastic criterion for engineering materials. Int J Plasticity 2004;20:2179–208. [2] Bardet JP. Lode dependences for isotropic pressure-sensitive materials. ASME, J Appl Mech 1990;57:498–506. [3] Bier W, Hartmann S. A finite strain constitutive model for metal powder compaction using a unique and convex single surface yield function. Eur J Mech A/Solids 2006;25:1009–30. [4] Houlsby GT. A general failure criterion for frictional and cohesive materials. Soils Found 1986;26(2):97–101. [5] Islam MK, Carter JP, Airey DW. Comparison of the yield locus and stressdilatancy function of some critical state constitutive models with experimental data for carbonate sand. J Inst Civil Eng, India 2004;84:267–74. [6] Kim MK, Lade PV. Modeling rock strength in three dimensions. Int J Rock Mech Min Sci Geomech Abstr 1984;21(1):21–33. [7] Kupfer H, Hilsdorf HK, Rusch H. Behaviour of concrete under biaxial stress. ACI J 1969;66(8):656–66. [8] Lade PV. Elastoplastic stress–strain theory for cohesionless soil with curved yield surfaces. Int J Solids Struct 1977;13:1019–35. [9] Lade PV. Instability, shear banding and failure in granular materials. Int J Solids Struct 2002;39:3337–57. [10] Lade PV, Duncan M. Elastoplastic stress–strain theory for cohesionless soil. ASCE, J GeoTech Eng Div 1975;101(GT10):1037–53.
967
[11] Lagioia R, Nova R. An experimental and theoretical study of a behaviour of calcarenite in triaxial compression. Geotechnique 1995;40(4):633–48. [12] Lagioia R, Puzrin AM, Potts DM. A new versatile expression for yield and plastic potential surfaces. Comput GeoTech 1996;19(3):171–91. [13] Lee SK, Song YC, Han SH. Biaxial behavior of plain concrete of nuclear containment building. Nucl Eng Des 2004;227:143–53. [14] Matsuoka H, Nakai T. Stress-deformation and strength characteristics of soil under three different principal stresses. Proc Japan Soc Civil Eng 1974;232:59–70. [15] Mogi K. Effect of the intermediate principal stress on rock failure. J Geophys Res 1967;72(20):5117–31. [16] Mogi K. Effect of the triaxial stress system on the failure of dolomite and limestone. Tectonophysics 1971;11:111–27. [17] Mogi K. Fracture and flow of rocks under high triaxial compression. J Geophys Res 1971;76(5):1255–69. [18] Mortara G. A new yield and failure criterion for geomaterials. Géotechnique 2008;58(2):125–32. [19] Nova R. Sinfonietta classica: an example on classical soil modelling. In: Saada, Bianchini, editors. Constitutive equations for granular non–cohesive soils. Cleveland, Rotterdam: Balkema; 1988. [20] Nova R, Wood DM. An experimental program to define the yield function for sand. Soils Found 1978;18(4):77–86. [21] Roscoe KH, Burland JB. On the generalized stress–strain behaviour of ‘wet’ clay. In: Heyman J, Leckie FA, editors. Engineering plasticity. Cambridge University Press; 1968. p. 535–609. [22] Zienkiewicz OC, Pande GN. Some useful forms of isotropic yield surfaces for soil and rock mechanics. In: Gudehus G, editor. Numerical methods in geomechanics, Vol. 5. J. Wiley & Sons; 1977. p. 179–90 [chapter 5].
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
A hierarchical single yield surface for frictional materials G. Mortara * Department of Mechanics and Materials, University of Reggio Calabria, Via Graziella, Feo di Vito, 89122 Reggio Calabria, Italy
a r t i c l e
i n f o
Article history: Received 23 January 2009 Received in revised form 2 March 2009 Accepted 4 March 2009 Available online 8 April 2009 Keywords: Constitutive modelling Three-dimensional stress Yield surface Shear strength Frictional materials
a b s t r a c t This paper extends the formulation of a recent yield and failure criterion [Mortara G. A new yield and failure criterion for geomaterials. Géotechnique 2008;58(2):125–32] by modifying its shape in the deviatoric representation. Such a modification improves the performance of the surface introducing a constitutive parameter associated to the variation of the Lode angle. The criterion is validated through the comparison with experimental data obtained from multiaxial tests on soils, rocks and concrete. A special form of the criterion will be also derived for numerical applications in order to remove singularities of the previous criterion. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction Elastoplastic models of geomaterials need the definition of yield and potential surfaces under three-dimensional stress conditions. The Mohr–Coulomb failure condition is often used in geomechanics even though it has corners that complicate the implementation in numerical codes. For this reason, over the years, smooth curves close to the Mohr–Coulomb criterion have been formulated [14,10,22]. Recently Mortara [18] formulated a new yield and failure criterion that generalises the criteria by Matsuoka and Nakai [14] and Lade and Duncan [10] and can include as well the criteria by Lade [8], Kim and Lade [6] and Houlsby [4]. The criterion can take into account tensile strength and different shapes in the meridian sections ranging from linear and curvilinear cones to closed domains. All surfaces have the characteristics to change their shape in the deviatoric plane for changing values of the friction angle /0 . In particular the continuous transition from circular to triangular shapes is observed for increasing values of / 0 as for the well known criteria by Matsuoka and Nakai [14] and Lade and Duncan [10]. The update version shown in this paper adds the possibility to modify the deviatoric shape of the criterion through the introduction of a correcting function dependent on Lode’s angle. Some changes are made in the expression of the previous criterion in order to avoid some drawbacks. Later, some comparisons will be presented between the criterion and experimental data based on triaxial and biaxial test on soils, rocks and concrete.
* Tel.: +39 0965 875 271; fax: +39 0965 875 201. E-mail address: [email protected] 0266-352X/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2009.03.007
2. General form of the criterion The criterion assumes the following general form:
f ¼ In1 Im m¼ 2 nI 3 ¼ 0
3n 2
ð1Þ
where n is the parameter determining the continuous transition of the surface shape in the deviatoric plane. The Matsuoka and Nakai [14] failure surface is obtained for n = 1 while for n = 3 criterion by Lade and Duncan [10] is obtained. I1 , I2 and I3 are the first, the second and the third stress invariants of the translated principal i ¼ ri þ t being t the absolute value of the isotropic tensile stresses r strength:
1 þ r 2 þ r 3 I1 ¼ r
1r 2r 3r 2 þ r 3 þ r 1 I2 ¼ r
1r 2r 3 I3 ¼ r
ð2Þ
Finally, n is the hierarchical function ruling both the size and shape of the criterion. Its general form is given by:
n ¼ nh þ an1 n2 n3
a ¼ n0 nh
ð3Þ
where n0 is the basic size and shape function dependent on the triaxial compression friction angle /00 and the parameter n; constant nh is the value that In1 Im 2 =I 3 takes in isotropic conditions; n1 and n2 are two correcting functions generating different shapes in the deviatoric planes; n3 is a correcting function introduced to provide the possibility of new shapes in deviatoric sections. The definition of the above mentioned functions is listed below.
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
n0 ¼
ð3 sin /00 Þn ð3 þ sin /00 Þm
ð4Þ
ð1 þ sin /00 Þð1 sin /00 Þ2m
nh ¼ 33m c1 I0 n1 ¼ I1 n2 ¼ 1 sgnðc2 Þ
ð5Þ ð6Þ I1 I0
!c 2 ð7Þ
3 n3 ¼ 1 þ e 1 cos h 2
ð8Þ
When c2 = 0 the geometrical meaning of constant I0 is the following: the curved surface obtained for c1 > 0 and the straight surface obtained for c1 = 0 intersect for I1 ¼ I0 . When c2 > 0, instead, constant I0 represents the value of I1 > 0 for which the surface crosses the hydrostatic axis. The Lode’s angle h in Eq. (8) is such that:
cos 3h ¼
pffiffiffiffiffiffi 27 J 3 2 J3
ð9Þ
where the deviatoric invariants J and J3 are defined as follows:
J¼
pffiffiffiffiffiffiffiffi 1 J2 J 2 ¼ I2 I21 3
1 2 3 J 3 ¼ I3 I2 I1 þ I 3 27 1
ð10Þ
Definition (9) provides h = 0° for triaxial compression and h = 60° for triaxial extension. In its general form the definition of the criterion needs seven parameters, namely /00 (or a), n, t, I0, c1, c2 and e. However, the basic form of the criterion, which is represented by a cone pinned at the origin of the principal stress space, requires only the two parameters /00 and n. The others can be introduced hierarchically depending on the experimental data that have to be represented. Introduction of the curvature in the meridian sections is obtained through parameters c1 and I0. Closed surfaces are instead obtained by choosing a non zero value for parameter c2. The geometrical meaning of the constants appearing in the yield function expression is reported in Mortara [18]. In order to show the different surfaces that can be obtained for different parameters, in Fig. 1 the three-dimensional representation of two surfaces is reported where the effect of taking c1 > 0 and c2 = 0 (Fig. 1a), and c2 > 0 (Fig. 1b) is clearly shown.
m 1 2 1 1 3 f ¼ In1 J 2 þ I21 n pffiffiffiffiffiffi cos 3hJ 3 I1 J 2 þ I1 ¼ 0 3 3 27 27
961
ð11Þ
Eq. (11) permits to observe that a third degree equation is obtained both for the Matsuoka and Nakai [14] and the Lade and Duncan [10] criteria. A third degree equation provides one or three real roots and for n = 1 or n = 3 Eq. (11) results in one negative and two positive values of J. In particular, for the Matsuoka and Nakai criterion as for the Lade and Duncan criterion it is possible to find stress states for which f < 0 even for points outside the current elastic domain [18]. This problem, particularly important for numerical applications, is common to different yield criteria [3]. For any non integer value pffiffiffi n, instead, Eq. (11) is not defined for J > J u ¼ I1 = 3, namely outside a circular cone in the pffiffiffi general stress space pffiffiffi having origin for I1 ¼ 0 and radius J u ¼ I1 = 3. Condition J I1 = 3 prevents the possibility to find more than a positive solution to f = 0 for stress states outside the current elastic domain. However taking non integer values for n does not avoid the possibility to find stress states such that J > Ju. In order to avoid discontinuities for the yield surface some comments should be made on the singularities of Eq. (1). First of all, if 1 > 0 and r 2 ¼ r 3 ¼ 0 it results f = 0 we apply Eq. (1) for a > 0, r independently of the value of / 0 given that this stress state corresponds to I1 > 0, I2 ¼ I3 ¼ 0. Nevertheless, this stress state should correspond to / 0 = 90°. If we limit ourselves to study the behaviour of frictional materials for /00 /00L being /00L a value close to 90° (e.g. /00L ¼ 89 ) we can remove pitfalls of Eq. (1), or equivalently of Eq. (11), for all values of n. For this purpose, Eq. (11) is rewritten in the form:
f ¼ f1 f2 nf3 þ I1N þ I2N ¼ 0
ð12Þ
where
n f1 ¼ hI1 i þ 1 sgnhI1 i 1 þ sgnhI1 i m 1 f2 ¼ hJ 2 þ I21 i þ 1 sgnhI2 i 3 f3 ¼
J3 J 3 þ 1 sgnJ
3. Special form of the criterion
I1N ¼
I1 þ jI1 j 2
Criterion (1) can be also expressed in terms of invariants I1 , J and h as follows
I2N ¼
I2 þ jI2 j 2
1 1 J 3 hI1 iJ 2 þ hI1 i3 3 27
Fig. 1. Effect of parameter c1 (a) and parameter c2 (b) in the principal stress space.
ð13Þ
962
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
The Macaulay brackets hi applied to any scalar x gives hxi ¼ ðx þ jxÞ=2. Variable J* in Eq. (13) is given by
J ¼ J L hJ L Ji
ð14Þ
Variable JL is the limit value for invariant J related to the chosen value of /00L that clearly depends on the Lode’s angle. For the sake of simplicity its value is found analysing the Lode dependence of the Mohr–Coulomb criterion for /00 ¼ /00L :
b 2 sin /00L 3 sin /00L J L ¼ pffiffiffi b¼ 0 hI 1 i 3 þ sin /00L 3b cos h þ ð2 bÞ sin h 3 sin /0L
/0h ¼ arcsin ð15Þ
Constant b is the aspect ratio, namely the ratio between the value of invariant J for h = 60° and h = 0° [2]. Variables I1N and I2N , instead, represent the value of f for negative values of I1 and I2 , respectively. For such values, in fact, I1N and I2N are equal to jI1 j and jI2 j, respectively and they are the only non null terms in Eq. (12) while for positive values of I1 and I2 they do not influence the value of f. The new form is completed by new expressions of functions n1 and n2. As can be noted analysing their expression, function (6) is discontinuous when I1 ¼ 0 like function (7) when I1 ¼ 0 and c2 = 0. For this reason the expressions of n1 and n2 are modified as follows:
n1 ¼
I0 I1L þ hI1 I1L i
!c 1 n2 ¼ 1 sgnðc2 Þ
I1L þ hI1 I1L i I0
!c 2
ð16Þ where I1L is an arbitrary small value of invariant I1 . In the updated form (12) no stress states can be found giving f < 0 outside the current elastic domain and the function is defined for all values of I1 , J, h and n. In order to understand the importance of form (12), Table 1 reports some example of stress states for which Eq. (11) fails contrarily to Eq. (12). 4. New shapes in the deviatoric sections: effect of function n3 Function n3 provides the possibility to modify the surface shape in the deviatoric sections. The simple idea to improve the capability of the new criterion in the deviatoric plane is to generate a surface with shape varying with the Lode’s angle h. For this purpose the correcting function n3, which is a trigonometric function of period 2p/3, is used. The role of parameters n and e is similar, given that their effect is to enlarge the deviatoric section with continuity from values of h ranging from 0° (triaxial compression) to 60° (triaxial extension). However, as shown in Fig. 2, they provide different shapes in the deviatoric section. Fig. 2a reports five functions drawn for /00 ¼ 30 , t = 100, e = 0 and n = 2, 4, 6, 8, 10 while Fig. 2b reports five functions drawn for /00 ¼ 30 , t = 100, n = 2 and e = 0, 0.2, 0.4, 0.6, 0.8. As it can be observed the surfaces for increasing value of n tend to become circular while those obtained for increasing value of e tend to assume a sharper profile in the deviatoric plane as the Lode’s angle tends to approach triaxial extension conditions. Fig. 3 shows the different shapes produced
Table 1 Values of Eqs. (11) and (12) for a > 0, t = c1 = c2 = e = 0 and for some values of n, r1, r2 and r3. n
r1
r2
r3
Eq. (11)
Eq. (12)
60 or >3 decimal <3 P3 Odd integer All but odd integer
0 <0 >0 >0 >0 >0
0
0
Undefined Undefined 0 Undefined <0 Undefined
0 I1 >0 >0 I2 I2
r1
r1
0 0 r1/2 r1/2
0 0 r1/2 r1/2
by parameters n and e in the biaxial plane. Parameters are the same as those of Fig. 2 and the figures represent the sections of the 3D surfaces with the plane r3 = 0. It can be still observed the different effects that parameters n and e produce on the biaxial plane especially as the Lode’s angle approaches 60°. Variation of parameters n and e influences the value of the friction angle /0h , which is related to the current state of principal stress through the relation:
r1 r3 r1 þ r3
ð17Þ
Fig. 4 shows the variation of /0h with the Lode’s angle h of seven curves: six are obtained for /00 ¼ 30 and different values of n and e while curve C0 refers to the Mohr–Coulomb criterion that, as well known, assumes the constant value /0h ¼ /00 . The maximum difference /0h /00 among the curves related to Eq. (1) and the Mohr–Coulomb curve is not obtained for the same value of the Lode’s angle: for e = 0, and irrespectively on the value /00 , such difference is attained for h<30° if n<3 and for h > 30° if n > 3. For the Lade and Duncan criterion, instead, the maximum difference /0h /00 is always obtained for h = 30°. The effect of parameter e on the shape in the deviatoric plane can be discussed by analysing curves C1, C3 and C4. Curve C1 is related to the Matsuoka–Nakai criterion while curve C3 is related to n = 1 and e = 0.5. The difference in terms of friction angle of these curves is evident and their maximum difference is attained for h = 60°. Curve C4 is related to e = 0 and n = 4.36 that is the value of n that assures the same value of shear strength for triaxial extension conditions of curve C3. It is evident how, although this two curves attain the same value of /0h for h = 60°, they produce different shapes in the deviatoric plane. Similar considerations follow for curves C2, C5 and C6. 5. Derivation of parameters and comparisons with experimental data In this section the proposed criterion is compared with experimental data related to triaxial and biaxial tests on soils, rocks and concrete. The values of parameters used to match the experimental curves are reported in Table 2. Concerning the derivation of parameters, when the experimental data show that a closed surface is necessary for their representation, it is necessary to specify a value for parameter c2. To derive parameters of the yield surface from triaxial compression or extension tests one can observe that Eq. (1) can be rewritten in the form
fg ¼ F
Rcp1 1 Rcp2
a¼0
ð18Þ
where
27 F ¼ 2 cos 3h 27
1 3
m g2 27
g3 13 g2 þ 1
nh
g¼
q p0 þ t
Rp ¼
p0 þ t p0 þ t
ð19Þ
For the stress state where the deviatoric stress is maximum (q = qg) condition ofg/og = 0 must be verified and this allows us to derive parameter c1
c1 ¼
Rcp2 gg @F c2 F @g 1 Rcp2
ð20Þ
Once the coordinates of the point in the q–p0 plane are known, it is possible to solve Eq. (18) that provides a value of a for a chosen value of c2. Fig. 5 shows the effect of increasing values of parameter c2 on the shape of the curves in the q–p0 plane for n = 1, t=p00 ¼ 0:1, qg =p00 ¼ 0:6 and p0g =p0 ¼ 0:7 (gg = 0.75). The inner curve is related to the value c2 = 2.649 that is the minimum value of c2 assuring
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
963
Fig. 2. Different shapes in the deviatoric sections: (a) increasing values of n for e = 0; and (b) increasing value of e for n = 2.
Fig. 3. Shapes of the criterion in the biaxial plane: (a) effect of parameter n; and (b) effect of parameter e.
c1 > 0. The set of shapes achievable in meridian planes by varying the values of parameters c1 and c2 makes the model similar to existing yield surfaces and plastic potentials in the literature [20,19,11,12,1,3]. An advantage of the proposed model resides in the opportunity to avoid the use of tension cut-off that instead is necessary in models like for example the modified Cam Clay [21]. On the other hand a possible drawback of the model is that, for a given set of parameters, the value of the mean stress p0 ¼ pg for which the deviatoric stress is maximum is not constant but increases with h. Such change increases for increasing values of gg. The procedure summarised by Eqs. (18) and (20) is applied to derive the parameters related to triaxial compression tests on carbonate sand reported by [5]. In Fig. 6 two sets of data are interpreted, the first one related to uncemented carbonate sand and the second related to carbonate sand with C = 20% of cementation. It is observed that the model curves provide a good interpretation of the experimental data. The results shown in Figs. 5 and 6 highlight the flexibility of the proposed model in the meridian representation. The derivation of parameter e is simple once the parameters related to the meridian representation of the surface are known. Con-
stant e can be easily found from Eq. (11) by providing the value of Je namely the value that invariant J assumes for h = 60°.
0 1
m 2 n 1 2 1 B I1 J e þ 3 I1 C e¼ n A1 @ an1 n2 p2ffiffiffiffi J3e 13 I1 J2e þ 271 I31 h 27
ð21Þ
In Fig. 7 the comparison among experimental true triaxial test data on dense Santa Monica Beach sand [9] and the criteria by Lade [8], Matsuoka and Nakai [14] and Eq. (1) with e > 0 is shown. The comparison is done both in terms of octahedral shear stresses (Fig. 7a) and Mohr–Coulomb friction angles (Fig. 7b). As it can be observed, Lade’s model overestimates the octaedral stresses, and then the friction angles, in the midrange of h values. Lade [9] imputed such overestimation to shear banding occurring in the sand specimens. Matsuoka and Nakai model [14] matches quite well the experimental data for h < 30° but underestimates the remaining stress points. The value of e for the new model curve is obtained from (21) for Je = 118 kPa. As one can observe the scatter between new model and experimental points is limited for all values of the Lode’s angle. Fig. 8 shows the experimental data on dolomite performed by Mogi [15] are compared with the new criterion and the models
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
f
964
Fig. 4. Influence of parameters n and e on the variation of the friction angle with the Lode’s angle for /00 ¼ 30 .
Table 2 Values of the parameters for the new criterion reported in Figs. 6–10. Figures
Experimental data
n
t (kPa)
a
c1
c2
I0 (kPa)
e
6 6 7 8 9a 9b 10a 10b
Islam et al. [5] (C = 0%) Islam et al. [5] (C = 20%) Lade [9] Mogi [15] Mogi [16] Mogi [16] Lee et al. [13] Lee et al. [13]
1.00 1.00 0.90 2.90 4.10 3.90 2.70 2.70
0:00 p0e 0:05 p0e 0.00 7150 9000 6700 3000 4000
10.00 13.98 8.25 9.0 108 1.25 105 5.70 104 1.10 105 1.40 105
0.795 0.563 0.144 1.800 2.000 2.000 1.150 1.150
0.100 0.500 0.000 0.000 0.000 0.000 0.000 0.000
1:00 p0e 1:00 p0e 100 100 100 100 100 100
0.000 0.000 0.665 0.800 0.900 0.900 0.450 0.450
Fig. 5. Effect of parameter c2 on the shape of the criterion in the q–p0 plane.
by Kim and Lade [6] and Mortara [18]. The new model, thanks to the introduction of parameter e, matches very well the data in the Rendulic plane (Fig. 8a) while the other two models are less
Fig. 6. Derivation of parameters from experimental data in the q–p0 plane.
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
965
Fig. 7. Comparison among true triaxial test data on dense Santa Monica Beach sand [9] and three failure criteria: (a) deviatoric plane; and (b) friction angle versus Lode’s angle.
Fig. 8. Comparison among three models and experimental data on dolomite performed by Mogi [15]: (a) Rendulic plane; and (b) biaxial plane.
satisfactory in extension conditions. The positive effect of function n3 is evident also in Fig. 8b that reports the experimental data and model curves in biaxial stress conditions. The curve related to the previous model [18] was obtained for n = 3.9 and improved the simulation by Kim and Lade (n = 3) but underestimating the uniaxial compressive strength. The new model curve, obtained for n = 2.9 and e = 0.8, improves the agreement with the experimental data for all values of the ratio r1/r2 providing at the same time a better estimation of the uniaxial compressive strength. Fig. 9 reports the comparison between experimental data by Mogi [16,17] and model. In particular, Fig. 9a shows the comparison between biaxial tests on Solenhofen limestone for three constant values of r3, while in Fig. 9b the comparison deals with tests on Mizuno trachyte for four constant values of r3. In both
cases the comparison is fairly good except for the model curve for r3 = 100 MPa in Fig. 9b that underestimates the values of experimental data. Finally, Fig. 10, deals with the comparison of the model with two sets of biaxial data on concrete performed for different biaxial stress ratios r1/r2 reported by Lee et al. [13]. The tests deal with two different values of the uniaxial compressive strength namely rc = 30.3 MPa and rc = 39.0 MPa reported in Fig. 10a and b, respectively. The experimental data show trends already observed in previous works like the one by Kupfer et al. [7]. These trends highlight that for stress ratio conditions ranging from uniaxial compression and uniaxial extension (0° 6 h 6 60°) the experimental envelopes show concavity. However, the comparison between data and model is satisfactory for both sets of data. It should be noted that only
966
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
Fig. 9. Comparison between new model and experimental data in the biaxial plane: (a) Solenhofen limestone [16]; and (b) Mizuno trachyte [17].
Fig. 10. Comparison between new model and experimental data on concrete by Lee et al. [13].
parameters t and a have been changed in the model curves reported in Fig. 10. 6. Conclusions This paper has dealt with the formulation of a new isotropic yield and failure criterion for frictional materials in terms of the three stress invariants. Starting from the recent model by Mortara [18], which formulation generalises well known criteria, a more general formulation has been derived by adding a function that increases the performance of the criterion in the deviatoric represen-
tation. The introduction of the new function that modifies the deviatoric section of the criterion has revealed an increasing capability to model the experimental data in both the deviatoric and the biaxial plane. The paper has reported the procedure to derive the constitutive parameters both in the deviatoric and meridian planes. The criterion has been validated using experimental data related to biaxial and triaxial tests on soils, rocks and concrete, showing an overall good attitude in modelling the behaviour of frictional materials. A special form of the criterion has been derived in order to remove all discontinuities observed in the previous formulation [18] mak-
G. Mortara / Computers and Geotechnics 36 (2009) 960–967
ing the new form suitable for the implementation in numerical codes. References [1] Aubertin M, Li L. A porosity-dependent inelastic criterion for engineering materials. Int J Plasticity 2004;20:2179–208. [2] Bardet JP. Lode dependences for isotropic pressure-sensitive materials. ASME, J Appl Mech 1990;57:498–506. [3] Bier W, Hartmann S. A finite strain constitutive model for metal powder compaction using a unique and convex single surface yield function. Eur J Mech A/Solids 2006;25:1009–30. [4] Houlsby GT. A general failure criterion for frictional and cohesive materials. Soils Found 1986;26(2):97–101. [5] Islam MK, Carter JP, Airey DW. Comparison of the yield locus and stressdilatancy function of some critical state constitutive models with experimental data for carbonate sand. J Inst Civil Eng, India 2004;84:267–74. [6] Kim MK, Lade PV. Modeling rock strength in three dimensions. Int J Rock Mech Min Sci Geomech Abstr 1984;21(1):21–33. [7] Kupfer H, Hilsdorf HK, Rusch H. Behaviour of concrete under biaxial stress. ACI J 1969;66(8):656–66. [8] Lade PV. Elastoplastic stress–strain theory for cohesionless soil with curved yield surfaces. Int J Solids Struct 1977;13:1019–35. [9] Lade PV. Instability, shear banding and failure in granular materials. Int J Solids Struct 2002;39:3337–57. [10] Lade PV, Duncan M. Elastoplastic stress–strain theory for cohesionless soil. ASCE, J GeoTech Eng Div 1975;101(GT10):1037–53.
967
[11] Lagioia R, Nova R. An experimental and theoretical study of a behaviour of calcarenite in triaxial compression. Geotechnique 1995;40(4):633–48. [12] Lagioia R, Puzrin AM, Potts DM. A new versatile expression for yield and plastic potential surfaces. Comput GeoTech 1996;19(3):171–91. [13] Lee SK, Song YC, Han SH. Biaxial behavior of plain concrete of nuclear containment building. Nucl Eng Des 2004;227:143–53. [14] Matsuoka H, Nakai T. Stress-deformation and strength characteristics of soil under three different principal stresses. Proc Japan Soc Civil Eng 1974;232:59–70. [15] Mogi K. Effect of the intermediate principal stress on rock failure. J Geophys Res 1967;72(20):5117–31. [16] Mogi K. Effect of the triaxial stress system on the failure of dolomite and limestone. Tectonophysics 1971;11:111–27. [17] Mogi K. Fracture and flow of rocks under high triaxial compression. J Geophys Res 1971;76(5):1255–69. [18] Mortara G. A new yield and failure criterion for geomaterials. Géotechnique 2008;58(2):125–32. [19] Nova R. Sinfonietta classica: an example on classical soil modelling. In: Saada, Bianchini, editors. Constitutive equations for granular non–cohesive soils. Cleveland, Rotterdam: Balkema; 1988. [20] Nova R, Wood DM. An experimental program to define the yield function for sand. Soils Found 1978;18(4):77–86. [21] Roscoe KH, Burland JB. On the generalized stress–strain behaviour of ‘wet’ clay. In: Heyman J, Leckie FA, editors. Engineering plasticity. Cambridge University Press; 1968. p. 535–609. [22] Zienkiewicz OC, Pande GN. Some useful forms of isotropic yield surfaces for soil and rock mechanics. In: Gudehus G, editor. Numerical methods in geomechanics, Vol. 5. J. Wiley & Sons; 1977. p. 179–90 [chapter 5].