- Email: [email protected]
An anisotropic hardening elastoplastic model for clays and sands and its application to FE analysis
Computers and Geotechnics 31 (2004) 37–46 www.elsevier.com/locate/compgeo
An anisotropic hardening elastoplastic model for clays and sands and its application to FE analysis D.A. Sun a
a,*
, H. Matsuoka a, Y.P. Yao b, H. Ishii
c
Department of Civil Engineering, Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya, Japan b Department of Civil Engineering, Beihang University, Beijing, China c Central Japan Railway Company, Nagoya, Japan Received 15 April 2003; received in revised form 18 October 2003; accepted 20 November 2003
Abstract The Sekiguchi–Ohta model is extended to be a unified three-dimensional elastoplastic model for clays, silts and sands by in~ij and a new hardening parameter H . The model can describe the negative and positive troducing a transformed stress tensor r dilatancy of soils with an initially stress-induced anisotropy in three-dimensional (3D) stress. An elastoplastic constitutive tensor is derived for the application of the proposed model to finite element (FE) analyses. A FE analysis example for a test embankment is given and the result demonstrates that the proposed model can predict reasonably the deformation of anisotropically consolidated clay and sand layers under embankment loading. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: Anisotropic; Constitutive model; Clay; Sand
1. Introduction Sekiguchi and Ohta [1] proposed an anisotropic hardening elastoplastic model for clays, which is an extension of the original Cam-clay model [2] to normally K0 -consolidated clays. This model has been widely applied in the finite element computation to predict the stress and deformation of earth structure in engineering practice in Japan [3]. Being the same as the Cam-clay model, the Sekiguchi–Ohta model also uses the Extended Mises criterion as shear yield and failure criteria, which cannot appropriately predict soil behaviors in three-dimensional (3D) stress state. There are some simple and reasonable shear yield criteria or failure criteria such as Matsuoka and NakaiÕs criterion called SMP criterion [4] and LadeÕs criterion [5] for soils in 3D stress. The SMP criterion is one of the best criteria for describing shear yield and failure behavior of soils in 3D
*
Corresponding author. Tel.: +81-52-735-7162; fax: +81-52-7235718. E-mail address: [email protected] (D.A. Sun). 0266-352X/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2003.11.003
stress. In this paper, the Sekiguchi–Ohta model is integrated with the SMP criterion in order to predict reasonably soil behaviors in 3D stress. The SMP criterion is introduced into the Sekiguchi–Ohta model by using a transformed stress method, which has been successfully applied to the integration of the Cam-clay model with the SMP criterion [6]. The transformed stress is deduced from what makes the SMP criterion become a cone with an axis being the space diagonal line in the transformed principal stress space. In this revised Sekiguchi–Ohta model, the SMP criterion is taken as a criterion of shear yield or failure, while the same model parameters are used as in the original Sekiguchi–Ohta model. Another shortcoming of the Sekiguchi–Ohta model is that it cannot predict positive dilatancy occurring in medium-to-dense granular materials such as sand. In this paper, after newly deriving the stress–dilatancy relation of the Sekiguchi–Ohta model, the model integrating the SMP criterion with the Sekiguchi–Ohta model is further extended to be a unified elastoplastic model for clays and sands by introducing a new hardening parameter H, which can describe the negative and positive dilatancy of soils with the initially
38
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
stress-induced anisotropy. Model complexity is controlled through only adding one material parameter to those of the Sekiguchi–Ohta model for practical use. The comparisons with experimental data of anisotropically consolidated-drained triaxial shear tests along different stress paths show that the proposed model gives excellent predictions for normally K0 -consolidated clays and sands. Moreover, in order to apply the presented simple model to practical engineering using FEM, its elastoplastic constitutive tensor is derived in detail. A finite element analysis of a test embankment was performed to demonstrate the usefulness of the model for predicting deformation of anisotropically consolidated clay and sand layers under embankment loading.
2. The Sekiguchi–Ohta model and its stress–dilatancy relation First of all, the Sekiguchi–Ohta model is briefly given for a well understanding of the proposed elastoplastic model. The difference between the Sekiguchi–Ohta model and the original Cam-clay model is only the stress ratio used in the models. In the Sekiguchi–Ohta model, the relative stress ratio g , in which the initial anisotropic stress state is taken into account, is used instead of the stress ratio gð¼ q=pÞ as used in the original Camclay model. The yield function (f ) and the plastic potential function (g) for the Sekiguchi–Ohta model are expressed as follows: kj p g f ¼g¼ ln þ ð1Þ epv ¼ 0; 1 þ e0 p0 M
If gk is defined as Eq. (4), the stress–dilatancy relation in the Sekiguchi–Ohta model can be newly obtained from differentiating Eq. (1) as Eq. (5) gk ¼
3 ðg gij0 Þgij ; 2g ij
gk ¼ M
depv ; depd
Fig. 1. Stress–dilatancy relation in Sekiguchi–OhtaÕs model (K0 ¼ 0:5).
with rij pdij ; p
gij0 ¼
rij0 p0 dij p0
ð3Þ
in which, rij0 is the value of rij at the end of the anisotropic consolidation and dij is KroneckerÕs delta. Whenqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi shearing starts from the isotropic stress g ¼ 3gij gij =2 ¼ q=p because gij0 ¼ 0, i.e., the Sekiguchi–Ohta model becomes the original Cam-clay model. It is worth noticing that the material parameters (k, j, M and m (PoissonÕs ratio)) in the Sekiguchi–Ohta model are the same as those in the Cam-clay model.
ð5Þ
where epd is the plastic deviatoric strain. Eq. (5) is similar in the form to the stress–dilatancy relation in the original Cam-clay model (i.e., g ¼ q=p ¼ M depv =depd ), and can be drawn as shown in Fig. 1 when K0 ¼ 0:5. Since gk will be used in the proposed model, the meaning of gk is explained as follows. Fig. 2 shows the contour lines of gk in the p-plane when K0 ¼ 0:5. From K0 to the point A along the triaxial compression stress path, gk changes from 0:75ð¼ 3ð1 K0 Þ=ð1 þ 2K0 Þ ¼ g0 Þ to M, i.e., gk ¼ g when g > g0 , and from K0 to the origin O along the triaxial compression stress path, gk changes from
where k and j are the slopes of normal consolidation and swelling lines, respectively, M is the slope of the critical state line, e0 is the initial void ratio of soil for p ¼ p0 , epv is the plastic volumetric strain that is used as a hardening parameter, and the relative stress ratio g is defined as [1] rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ðg gij0 Þðgij gij0 Þ g ¼ ð2Þ 2 ij
gij ¼
ð4Þ
Fig. 2. Contour lines of g1 in p-plane (K0 ¼ 0.5).
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
0:75ð¼ g0 Þ to 0, i.e., gk ¼ g when g < g0 . From the origin O to the point B along the triaxial extension stress path, gk changes from 0 to M, i.e., gk ¼ g. In other cases, gk changes from some values in the range of 0:75 0.75 to M when K0 ¼ 0:5. In addition, as is well known, the contour lines of the stress ratio gð¼ q=pÞ are concentric circles with the centers being the origin in the p-plane, and the contour lines of the relative stress ratio g are concentric circles with the centers being the K0 point shown in Fig. 2 in the p-plane. Hence, the stress ratio gk can be explained as an intermediate stress ratio between the relative stress ratio g and the stress ratio g.
3. An anisotropic hardening elastoplastic model for clays and sands As explained before, the Extended Mises criterion with the initial stress-anisotropy (g ¼ constant) is adopted as a shear yield criterion in the Sekiguchi–Ohta model. As is well known, the Extended Mises criterion is not appropriated for describing soil behavior and the SMP criterion is one of the best criteria for describing shear yield and failure behavior of soils in 3D stress. The SMP criterion has been successfully applied to the Cam~ij inclay model by introducing a transformed stress r stead of the conventional stress rij [6]. In the presented study, the model is generalized using the transformed stress based on the SMP criterion. The SMP criterion can be written by [4] I1 I2 ¼ constant; I3
ð6Þ
where I1 , I2 and I3 are the first, second and third stress invariants, respectively. The solid curve in Fig. 3 is the shape of the SMP criterion in the p-plane. The transformed stress is deduced from what makes the SMP curve in the p-plane become a circle with the center being the origin in the transformed p-plane, as shown in Fig. 3. From the geometrical relationship between two points A and A0 shown in Fig. 3, the transformed σ1 (σ1 ) Circle
A
SMP criterion
A'
lθ
θ (θ)
0
σ 2 (σ 2 )
stress can be calculated using the conventional stress as follows: ‘0 ~ij ¼ pdij þ pffiffiffiffiffiffiffiffiffiffi sij ; r skl skl
σ 3 (σ 3)
~-plane (dotted Fig. 3. The SMP criterion in p-plane (solid line) and p line).
ð7Þ
where sij is the deviatoric stress, and ‘0 is written by (Fig. 3) [6] rffiffiffi 2 I1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ‘0 ¼ 2 : ð8Þ 3 3 ðI1 I2 I3 Þ=ðI1 I2 9I3 Þ 1 If we take the place of the stress tensor rij with the ~ij , and introduce a new transformed stress tensor r hardening parameter H instead of the original hardening parameter epv to the Sekiguchi–Ohta model, the following yield function (f ) and plastic potential function (g) can be obtained ! ~p ~g kj f ¼g¼ ln þ H ¼ 0; ð9Þ ~p0 M 1 þ e0 ~ii =3) at the initial shear where p~0 is the value of ~pð¼ r state with the void ratio e0 , and M is explained as the value of q=p at the characteristic state at which the rate of plastic volume change is zero, i.e., depv ¼ 0, for dilative soil such as medium-to-dense sand and at the critical state for contractive soil such as normally consolidated clay in triaxial compression stress [7]. Although the value of M for sand varies with the confining pressure [8], M is usually assumed to be constant for the sake of simplicity and practical use. Moreover, ~ g is defined to be the same in the form as Eq. (2) rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ~g ¼ ð~g ~gij0 Þð~gij ~gij0 Þ ð10Þ 2 ij in which ~ij ~pdij r ~gij ¼ ; ~p
ð11Þ
~gij0 is the value of ~gij at the initial shear state. It can be seen from Eqs. (9) and (10) that the proposed yield surface is anisotropic because the initial stress ratio tensor ~gij0 is fixed. Following a hardening parameter for the isotropic hardening model proposed by authors [9], the increment of the hardening parameter H for the anisotropic hardening is defined as dH ¼
lθ
39
M 4 Mf4 ~g4k p de ; Mf4 M 4 ~g4k v
ð12Þ
where Mf is the value of q=p at failure state (epd ! 1) in triaxial compression stress, and is assumed to be constant for the sake of simplicity and practical use. ~ gk is defined in the same form as Eq. (4). ~gk ¼
3 ð~g ~gij0 Þ~gij : 2~g ij
ð13Þ
40
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
Fig. 4. Proposed stress–dilatancy relation (K0 ¼ 0:5).
Fig. 7. Yield loci in p-plane determined by undrained cyclic loading true triaxial tests on sand (after [11]).
Fig. 5. Proposed yield curve (fine solid line), failure state lines (bold solid line) and characteristic lines (dotted line) in triaxial stress.
By using ~ gk the stress–dilatancy relation in the proposed model is written as follows: ~ gk ¼ M
depv : depd
ð14Þ
Fig. 4 shows the shape of Eq. (14) when K0 ¼ 0:5. The hardening parameter H defined by Eq. (12) is a combination of both plastic volumetric strain epv and plastic shear strain epd , a concept first introduced by
Nova and Wood [10], after taking Eq. (14) into consideration. Since H is a hardening parameter, dH is not less than zero. Taking this into account, the following features of H can be deduced from Eq. (12): (1) When ~ gk ¼ 0 or Mf ¼ M, dH ¼ depv , i.e., the hardening parameter H becomes the plastic volumetric strain epv , which is the same hardening parameter as used in the Sekiguchi– Ohta model and Cam-clay model. (2) When 0 < ~ gk < M, depv > 0; when ~gk ¼ M, depv ¼ 0; and when M < ~ gk < Mf , depv < 0, i.e., the negative and positive dilatancy of soils can be described. In the proposed model, the associated flow rule is ~ij space [6], i.e., adopted in the transformed stress r depij ¼ K
of ; o~ rij
Fig. 6. Yield loci of the proposed model in deviatoric planes.
ð15Þ
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
where K is the proportionality constant, and can be determined from the consistency condition. Fig. 5 shows the proposed yield curve (Eq. (9), fine solid curve), failure state lines (~ gk ¼ Mf , bold solid line) and characteristic lines ~ g ¼ M, dotted line) in triaxial compression and triaxial extension stresses. From Eq. (12), it is known that when Mf ¼ M the failure state line becomes the same as the characteristic state line, and the hardening parameter H becomes the plastic volumetric strain epv , which is widely used as the hardening parameter in the elastoplastic model such as the Cam-clay model and the Sekiguchi–Ohta model for normally consolidated clay. Fig. 6 shows the yield loci of the proposed model in the deviatoric planes. The yield loci in the transformed p-plane (~ p-plane) are concentric circles with the centers being the point K0 as shown in Fig. 6(a), which are the same in the shape as those of the Sekiguchi–Ohta model in the p-plane, while the yield loci become oval-shaped
41
loci in the ordinary p-plane as shown in Fig. 6(b). It is interesting to note that in the p-plane the shape of the yield loci (Fig. 6(b)) is similar to that obtained using a true triaxial test apparatus by Yamada [11], as shown in Fig. 7. The yield points marked with and are determined using the method similar to those by Poorooshasb [12] and Tatsuoka [13]. All specimens were first subjected to a deviatoric stress along y-direction (i.e. h ¼ 0°) from O to A under Constant mean stress p. This stress path corresponds to that in conventional triaxial compression test. After completion of the first cycle along triaxial stress path the specimens were subjected to a loading along the stress path of h ¼ 15° or 30° or 45° to determine the yield stress points marked with from the points where the curvature of shear stress versus shear strain curve is maximum. The method for determining the yield points marked with is similar to the above, but the first experienced deviatoric stress was larger. The different yield loci shown in Fig. 7 are the
5 σ r /σ
a
4 3 r
Comp. R=σ a/σ r=1.5
1 1.5
-10
εa (%)
σ a-σ
2
Proposed S-O model Observed
1
0
-5
0
5
εr (%)
10
p
εv (%)
Ext. 5
Failure
R=1.5
(a) 5 a
4 r
σ r /σ
2 1 2.0
Proposed S-O model Observed
1
-10 ε (%) -5 a
0
σ a-σ
3
Comp. R=2.0
0 5
p
εr (%) 10
Ext.
εv (%)
Failure
5
(b)
R=2.0 5 a
4 r
σ r /σ
2 1 2.5
-10
εa (%)
-5
Proposed S-O model Observed
1 0
σ a-σ
3
5
ε r (%)
0 p
10
Ext.
εv (%)
Failure
5
(c)
Comp. R=2.5
R=2.5
Fig. 8. Comparison between predicted and experimental stress–strain behavior during shearing from anisotropic stresses (experimental data after [14]).
42
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
observed results of sand samples subjected to different values of deviatoric stress in the y-direction. It can be seen from Fig. 6(b) that the yield loci in real stress space are partly non-convex. This is because the convexity of the yield surface is satisfied in the transformed stress space, and the associated flow rule is adopted also in the transformed stress space. That is to say, in the proposed model DruckerÕs postulate is satisfied in the transformed stress space but not in real stress space.
4. Comparisons of model prediction with experimental result The model contains six parameters k, j, M, Mf , m and e0 . In comparison with the Sekiguchi–Ohta model, only one model parameter, Mf , is increased in the proposed model. However, for the soils exhibiting only negative dilatancy, i.e., contractive soils, the model parameters in two models are the same because Mf ¼ M in this case. Fig. 8 shows the comparisons of the predicted results by the Sekiguchi–Ohta model (abbreviated as S–O model in the figure; dotted line) and the proposed model (solid line) with the results of triaxial tests (mark ) on a contractive soil composed of Toyoura sand and Fujinomori clay (experimental data from [14]). Fujinomori clay is a silty clay and has been used as a normally consolidated clay in other soil tests (e.g. [15]). After Toyoura sand was mixed with Fujinomori clay with a ratio of 7:3 by dry weight, the mixture was consolidated up to 49 kPa for triaxial tests. In Fig. 8, ra and rr are the axial and radial stresses, and ea and er are the axial and radial strains in
σ a /σ
triaxial test, respectively. The material parameters used in the model prediction are as follows: k=ð1 þ e0 Þ ¼ 0:021, j=ð1 þ e0 Þ ¼ 0:0034, Mf ¼ M ¼ 1:32 and m ¼ 0:3. The triaxial specimens were first anisotropically consolidated up to p ¼ 980 kPa along the stress paths of ra =rr ¼ 1:5, 2.0 and 2.5, respectively, in triaxial compression and then sheared to triaxial extension up to failure under p ¼ 980 kPa as shown on the right-hand side of Fig. 8. The lefthand side of Fig. 8 shows the predictions and triaxial test results during shearing. It can be seen from Fig. 8 that the proposed model can predict the test results better than the Sekiguchi–Ohta model. Fig. 9 shows the experimental and predicted results of triaxial compression and triaxial extension tests on Toyoura sand from K0 ( ¼ 0.5)-consolidation state (experimental data from [16]). The stress paths tested are from point P at ra =rr ¼ 2 in triaxial compression to triaxial compression or triaxial extension up to failure under p ¼ 196 kPa or rr ¼ 147 kPa, as shown in Fig. 10. The material parameters used in the model prediction are k=ð1 þ e0 Þ ¼ 0:0508, j=ð1 þ e0 Þ ¼ 0:0112, M ¼ 0:95, Mf ¼ 1:66 and m ¼ 0:3. It can be seen from Fig. 9 that the model can properly predict the sand behavior including dilatancy characteristics. Figs. 8 and 9 demonstrate that the proposed model can predict well the stress–strain relation of not only clays but also sands during shearing from the initially anisotropic stress to triaxial compression stress or triaxial extension stress. Therefore, the proposed model is capable of predicting well the stress–strain relations of contractive and dilative soils shearing from an anisotropic stress state to any stress state in 3D stress.
6
σ r /σ
r
6 a
4
2 1 0
Path PA -10
ε r (%)
-5
4
Path PC 5
ε v (%) 4 σ a /σ
ε a (%)
ε a (%)
10 -5
1 2.0
6
σ r /σ
4
-10
ε r (%)
-5 ε v (%)
2 1 0 4
ε r (%)
0 2
ε v (%)
r
Path PB
1
55
6 a
4
Path PD 5
ε a (%)
ε a (%)
10 -5
1 2.0
1 0 ε (%) v 2
ε r (%) 55
Fig. 9. Comparison between predicted and experimental stress–strain behavior of Toyoura sand during shearing from an anisotropic stress (experimental data after [16]).
800
Dijkl ¼ Deijkl Deijmn
43
of of e D =X : o~ rmn orst stkl
ð21Þ
Introducing the elastic constitutive tensor of HookeÕs law for isotropic elasticity into Eqs. (19) and (21), we can obtain X ¼
R= 1
600
A R= 2
failure line
σa (kPa)
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
B 400
k j M3 Mf4 g~4k of of þL 4 2 2 1 þ e0 Mf ðM þ ~gk ÞðM þ ~gk Þ~p orii o~ rjj of of þ 2G ; orij o~ rij
ð22Þ
P
Dijkl ¼ Ldij dkl þ G dik djl þ dil djk 200
D
C
of of dij þ 2G L o~ rmm o~ rij
line failure
!
of of L dkl þ 2G ornn orkl
, X; ð23Þ
0
200
400
600
where L and G are LameÕs constants, which can be written from e ¼ e0 j ln p=p0 in elastic region by
σr (kPa) Fig. 10. Stress paths for Toyoura sand.
5. Elastoplastic constitutive tensor The elastoplastic constitutive tensor for the proposed model to solve elastoplastic boundary value problem by means of the finite element method is derived here. The elastic part of the stress–strain relation can be written in the incremental form as drij ¼ Deijkl deekl ¼ Deijkl ðdekl depkl Þ;
ð16Þ
where Deijkl is the elastic constitutive tensor. ekl and epkl are the total strain tensor and its plastic component, respectively. Because the yield function f is the function of the stress tensor and the hardening parameter, the consistency condition can be written by df ¼
of of dH ¼ 0: drij þ orij oH
ð17Þ
Substituting Eqs. 12,15 and 16 into Eq. (17) gives K¼
of orij
Deijkl dekl X
;
ð18Þ
where X ¼
M 4 Mf4 ~ g4k of of e of dij þ D : 4 4 4 rij orij ijkl o~ rkl Mf M ~ gk o~
ð19Þ
Substituting Eqs. (15) and (18) into Eq. (16), we get a general form of the proposed model drij ¼ Dijkl dekl ;
ð20Þ
where the elastoplastic constitutive tensor is expressed as
G¼
E 3ð1 2mÞð1 þ e0 Þ ¼ p; 2ð 1 þ m Þ 2ð1 þ mÞj
ð24Þ
L¼
E 2 ð1 þ e0 Þ 2 G¼ p G; 3ð1 2mÞ 3 j 3
ð25Þ
where E and m are the elastic modulus and PoissonÕs ratio, respectively.
6. Finite element computation The proposed model is incorporated into the FEM program DACSAR coded by Iizuka and Ohta in Japan [17]. As shown in Figs. 2 and 5, when the stress path is along the K0 -line, the plastic flow at the singular point of non-smooth yield surface is unable to be correctly evaluated due to an unpredictable gradient at the corner. A very simple handling method adopted in the FEM code is that K0 is assumed to be slightly smaller than the stress ratio (r2 =r1 ¼ r3 =r1 ) when the stress state is just on the K0 -line. For the singularity of the Sekiguchi–Ohta model, more precious numerical method has been proposed by Yuttapontada et al. [18] based on KoiterÕs associated flow rule. The governing equations describing the behavior of the soil skeleton coupled with the pore water flow are formulated by Biot [19]. The technique employed in discretizing the pore water flow was originally proposed by Akai and Tamura [20]. As an example of the application of the proposed model to practical engineering, as shown in Fig. 11, a natural ground under embankment loading [17] is analyzed using finite element method. The use of the earth embankment is for carrying a highway. The ground
44
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
36. 00m 10. 26m
15. 48m
10.26m
GL-0.0m GL-1.4m GL-7.8m
AC1 AS
AC2 GL-23.3m Fig. 11. Soil profile of ground and embankment.
X fixed boundary
X fixed boundary
23.3m 10m
Y
X
X, Y fixed boundary
150m Fig. 12. Finite element meshes and boundary conditions.
2
q/p
Mf=1.79 M=1.45
εv (%)
0
5 ε1 (%)
k=ð1 þ e0 Þ j=ð1 þ e0 Þ M Mf m k (m/day)
1
1.30 10 1.30 102 1.40 1.40 0.30 2.0 104
Determination of M and Mf 0 Cc 1+e0
1+e0 As (sand)
0.02 (b)
=0.0217 Measured Predicted
0.01
Cs
9.44 103 4.94 103 1.45 1.79 0.30 1.0 102
10
1
(a)
Table 1 Material parameters used in FE computation Ac2 (clay)
p=392kPa Measured Predicted
1
εv
consists of Ac1 layer (surface clay layer: GL0 )1.4 m), As layer (sand layer: GL-1.4 )7.8 m) and Ac2 layer (uniform marine clay layer, GL-7.8 )23.3m). The ground table is at 0.5 m below the ground surface. The fill material is the dredged sand with a unit weight of 1.73 tf/m3 . Fig. 12 shows the finite element meshes and boundary conditions. Only half section is computed for the symmetry of earth structure. The ground surface and the bottom boundary were treated as drained boundaries, while the two sides of the ground were treated as undrained boundaries. The two sides were restrained in the horizontal direction, while the bottom boundary was fixed in both the vertical and horizontal directions. Table 1 shows the values of the material parameters used in the finite element computation. These values are determined from the laboratory tests on the samples taken from the site of the embankment. In detail, the
=0.0114
100
200
400 p (kPa) Determination of Cc and Cs
1000
Fig. 13. Determination of model parameters for the soil As (experimental data after [21]).
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
0
Predicted Observed
100cm
0 5 10
Fill height (m)
Fig. 15. Surface settlement just after embankment construction.
pw(kPa)
150 100
Depth 9.0m
Predicted Observed
50 0
pw(kPa)
model parameters for the soil Ac2 were determined by using Cc =ð1 þ e0 Þ ¼ 0:3, Cs =ð1 þ e0 Þ ¼ 0:03 and /0 ¼ 34:5°, which are obtained from results of one-dimension consolidation tests and drained triaxial tests on the soil Ac2 [21]. Here, Cc and Cs are compression index and swelling index, respectively; e0 is the initial void ratio and /0 is the effective angle of friction. The model parameters for the soil As are determined from Fig. 13, where the experimental data are from [21] and the solid lines show the model predictions. The value of Mf is determined based on /0 ¼ 43:5° for the soil As that was obtained from the envelope of three MohrÕs stress circles at failure [21]. The values of the coefficient of permeability for the soils As and Ac2 were taken from [17,21], respectively. The surface clay layer (Ac1) is assumed as an elastic material. The observed and computed settlements at the ground surface and depth of 7.8 m (i.e., the top of the soil Ac2) under the center of the embankment are plotted against time during and after embankment construction in Fig. 14. In all the figures, the plot is the observed results and the solid line is the results computed by the finite element method using the proposed
45
200 150 Depth 11.5m 100 50
pw(kPa)
5 10
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1. 6 0
Predicted Observed
At ground surface
100
200
300 Time (day)
0 Predicted Observed
Settlement (m)
0.2 0.4 0.6 0.8 At GL-7.8
(b)
0
200 150 Depth 21.5m 100 50 0
(a)
1
200 150 Depth 15.0m 100 50 0
pw(kPa)
Settlement (m)
Fill height (m)
0
100
200
300 Time (day)
Fig. 14. Loading process and settlements under the center of the embankment.
50
100
150
200
250 300 Time (day)
Fig. 16. Generation and dissipation of excess pore water pressure at different depths of the centerline.
model. Fig. 15 shows the observed and computed ground surface settlements at the time of just completing the embankment construction. The excess pore water pressures induced by embankment loading at different depths of central line are plotted against time during and after embankment construction in Fig. 16. From Figs. 14–16, it can be seen that the settlements and the excess water pressures computed using the proposed model fairly agree with the observed results. However, the computed excess pore water pressure at a depth of 21.5 m is lower than the observed value (Fig. 16). This is due to adopting the treatment of the bottom boundary as a drained boundary in the finite element computation while the true bottom boundary may be of partial drainage.
46
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
7. Conclusions This paper presented an anisotropic hardening model for clays, silts and sands within a simple framework, and the finite element analyses for an embankment on the ground consisting of clays and sands using the proposed model. Several conclusions can be drawn from this study. 1. A unified and simple elastoplastic model for clays and sands was proposed, which is capable of describing deformation and strength characteristics of clays and sands with the initially anisotropic consolidation state in 3D stress. 2. For contractive soils, the model parameters involved in the proposed model are the same in the Cam-clay model or the Sekiguchi–Ohta model. For dilative soil, only one model parameter is added to the model parameters of the Cam-clay model or the Sekiguchi– Ohta model. 3. By using the given elastoplastic constitutive tensor, the proposed model can be easily applied to the finite element computation for earth structure composed of sands and clays. 4. The yield locus for soils subjected to the anisotropic stress should be oval-shaped in the p-plane. Acknowledgements The authors are grateful to Dr. S.H. Liu of Nagoya Institute of Technology and Dr. H.L. Fang of West Japan Engineering Consultants for their help in the revision of the manuscript, and to Prof. A. Iizuka of Kobe University for providing the program DACSAR. References [1] Sekiguchi H, Ohta H. Induced anisotropy and time dependency in clays. In: Murayama S, Schofield AN, editors. Proceedings of the 9th ICSMFE 1977; Specialty Session 9. p. 229–38. [2] Roscoe KH, Schofield AN, Thurairajoh A. Yielding of clay in state wetter than critical. Geotechnique 1963;13(3):211–40. [3] Duncan JM. The role of advanced constitutive relations in practical applications. In: Proceedings of the 13th International Conference on Soil Mechanics and Foundation Engineering, vol. 5; 1994. p. 31–48. [4] Matsuoka H, Nakai T. Stress-deformation and strength characteristics of soil under three different principal stresses. Proc JSCE 1974;232:59–74.
[5] Lade PV, Duncan JM. Elasto-plastic stress–strain theory for cohesionless soil. J Geotech Eng ASCE 1975;101(10):1037–53. [6] Matsuoka H, Yao YP, Sun DA. The Cam-clay models revised by the SMP criterion. Soils Foundations 1999;39(1):81–95. [7] Luong MP. Stress–strain aspects of cohesionless soils under cyclic and transient loading. In: International Symposium on Soils under Cyclic and Transient Loading; 1982. p. 315–24. [8] Lade PV, Ibsaen LB. A study of the phase transformation and the characteristic lines of sand behavior. In: Asaoka et al., editor. International Symposium Deformation and Progressive Failure in Geomechanics; 1997. p. 353–8. [9] Yao YP, Matsuoka H, Sun DA. A unified elastoplastic model for clay and sand with the SMP criterion. In: Vitharana N, Colman R, editors. Proceedings of 8th Australia–New Zealand Conference on Geomechanics, vol. 2; 1999. p. 997–1004. [10] Nova R, Wood DM. A constitutive model for sand in triaxial compression. Int J Numer Analytical Meth Geomech 1979;3(3): 255–78. [11] Yamada K. Deformation characteristics of loose sand under three-dimensional stress. Doctoral thesis at University of Tokyo; 1979 [in Japanese]. [12] Poorooshash HB, Holubec I, Sherbourne AN. Yielding and flow of sand in trixial compression. Can Geotech J 1967;4(4):376–98. [13] Tatsuoka F, Ishihara K. Yielding of sand in triaxial compression. Soils Foundations 1974;14(2):63–76. [14] Ichihara W. Deformation and strength characteristics of intermediate soil under triaxial stress condition. BachelorÕs Thesis at Nagoya Institute of Technology; 1997 [in Japanese]. [15] Nakai T, Matsuoka H, Okuno N, Tsuzuki K. True triaxial tests on normally consolidated clay and analysis of the observed shear behavior using elastoplastic constitutive models. Soils Foundations 1986;26(4):67–78. [16] Funada T, Matsuoka H, Fukumoto S. Deformation characteristics of sand under various stress paths from anisotropic consolidation. In: Proceedings of the 44th Annual Conference of the Japan Society of Civil Engineers, vol. 3; 1989. p. 512–3 [in Japanese]. [17] Iizuka A, Ohta H. A determination procedure of input parameters in elasto-viscoplastic finite element analysis. Soils Foundations 1987;27(3):71–87. [18] Yuttapontada Y, Pipstpongsa T, Iizuka A, Ohta H. FEM for the Sekiguchi–Ohta model consider non-smooth yield surface. In: Proceedings of the 38th Japan National Conference on Geotechnical Engineering, 2003. p. 285-6. [19] Biot MA. General theory of three-dimensional consolidation. J Appl Phys 1941;12:155–64. [20] Akai K, Tamura T. Numerical analysis of multi-dimensional consolidation accompanied with elasto-plastic constitutive equation. Proc JSCE 1978;269:95–104 [in Japanese]. [21] Matsuoka H, Nakai T. Prediction using constitutive equation for soils based on the concept of Spatially Mobilized Plane. Research Report on settlement prediction of ground due to embankment edited by the committee on investigating highway, the Japan Highway Public Corporation 1984. p. 220–42 [in Japanese].
An anisotropic hardening elastoplastic model for clays and sands and its application to FE analysis D.A. Sun a
a,*
, H. Matsuoka a, Y.P. Yao b, H. Ishii
c
Department of Civil Engineering, Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya, Japan b Department of Civil Engineering, Beihang University, Beijing, China c Central Japan Railway Company, Nagoya, Japan Received 15 April 2003; received in revised form 18 October 2003; accepted 20 November 2003
Abstract The Sekiguchi–Ohta model is extended to be a unified three-dimensional elastoplastic model for clays, silts and sands by in~ij and a new hardening parameter H . The model can describe the negative and positive troducing a transformed stress tensor r dilatancy of soils with an initially stress-induced anisotropy in three-dimensional (3D) stress. An elastoplastic constitutive tensor is derived for the application of the proposed model to finite element (FE) analyses. A FE analysis example for a test embankment is given and the result demonstrates that the proposed model can predict reasonably the deformation of anisotropically consolidated clay and sand layers under embankment loading. Ó 2004 Elsevier Ltd. All rights reserved. Keywords: Anisotropic; Constitutive model; Clay; Sand
1. Introduction Sekiguchi and Ohta [1] proposed an anisotropic hardening elastoplastic model for clays, which is an extension of the original Cam-clay model [2] to normally K0 -consolidated clays. This model has been widely applied in the finite element computation to predict the stress and deformation of earth structure in engineering practice in Japan [3]. Being the same as the Cam-clay model, the Sekiguchi–Ohta model also uses the Extended Mises criterion as shear yield and failure criteria, which cannot appropriately predict soil behaviors in three-dimensional (3D) stress state. There are some simple and reasonable shear yield criteria or failure criteria such as Matsuoka and NakaiÕs criterion called SMP criterion [4] and LadeÕs criterion [5] for soils in 3D stress. The SMP criterion is one of the best criteria for describing shear yield and failure behavior of soils in 3D
*
Corresponding author. Tel.: +81-52-735-7162; fax: +81-52-7235718. E-mail address: [email protected] (D.A. Sun). 0266-352X/$ - see front matter Ó 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2003.11.003
stress. In this paper, the Sekiguchi–Ohta model is integrated with the SMP criterion in order to predict reasonably soil behaviors in 3D stress. The SMP criterion is introduced into the Sekiguchi–Ohta model by using a transformed stress method, which has been successfully applied to the integration of the Cam-clay model with the SMP criterion [6]. The transformed stress is deduced from what makes the SMP criterion become a cone with an axis being the space diagonal line in the transformed principal stress space. In this revised Sekiguchi–Ohta model, the SMP criterion is taken as a criterion of shear yield or failure, while the same model parameters are used as in the original Sekiguchi–Ohta model. Another shortcoming of the Sekiguchi–Ohta model is that it cannot predict positive dilatancy occurring in medium-to-dense granular materials such as sand. In this paper, after newly deriving the stress–dilatancy relation of the Sekiguchi–Ohta model, the model integrating the SMP criterion with the Sekiguchi–Ohta model is further extended to be a unified elastoplastic model for clays and sands by introducing a new hardening parameter H, which can describe the negative and positive dilatancy of soils with the initially
38
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
stress-induced anisotropy. Model complexity is controlled through only adding one material parameter to those of the Sekiguchi–Ohta model for practical use. The comparisons with experimental data of anisotropically consolidated-drained triaxial shear tests along different stress paths show that the proposed model gives excellent predictions for normally K0 -consolidated clays and sands. Moreover, in order to apply the presented simple model to practical engineering using FEM, its elastoplastic constitutive tensor is derived in detail. A finite element analysis of a test embankment was performed to demonstrate the usefulness of the model for predicting deformation of anisotropically consolidated clay and sand layers under embankment loading.
2. The Sekiguchi–Ohta model and its stress–dilatancy relation First of all, the Sekiguchi–Ohta model is briefly given for a well understanding of the proposed elastoplastic model. The difference between the Sekiguchi–Ohta model and the original Cam-clay model is only the stress ratio used in the models. In the Sekiguchi–Ohta model, the relative stress ratio g , in which the initial anisotropic stress state is taken into account, is used instead of the stress ratio gð¼ q=pÞ as used in the original Camclay model. The yield function (f ) and the plastic potential function (g) for the Sekiguchi–Ohta model are expressed as follows: kj p g f ¼g¼ ln þ ð1Þ epv ¼ 0; 1 þ e0 p0 M
If gk is defined as Eq. (4), the stress–dilatancy relation in the Sekiguchi–Ohta model can be newly obtained from differentiating Eq. (1) as Eq. (5) gk ¼
3 ðg gij0 Þgij ; 2g ij
gk ¼ M
depv ; depd
Fig. 1. Stress–dilatancy relation in Sekiguchi–OhtaÕs model (K0 ¼ 0:5).
with rij pdij ; p
gij0 ¼
rij0 p0 dij p0
ð3Þ
in which, rij0 is the value of rij at the end of the anisotropic consolidation and dij is KroneckerÕs delta. Whenqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi shearing starts from the isotropic stress g ¼ 3gij gij =2 ¼ q=p because gij0 ¼ 0, i.e., the Sekiguchi–Ohta model becomes the original Cam-clay model. It is worth noticing that the material parameters (k, j, M and m (PoissonÕs ratio)) in the Sekiguchi–Ohta model are the same as those in the Cam-clay model.
ð5Þ
where epd is the plastic deviatoric strain. Eq. (5) is similar in the form to the stress–dilatancy relation in the original Cam-clay model (i.e., g ¼ q=p ¼ M depv =depd ), and can be drawn as shown in Fig. 1 when K0 ¼ 0:5. Since gk will be used in the proposed model, the meaning of gk is explained as follows. Fig. 2 shows the contour lines of gk in the p-plane when K0 ¼ 0:5. From K0 to the point A along the triaxial compression stress path, gk changes from 0:75ð¼ 3ð1 K0 Þ=ð1 þ 2K0 Þ ¼ g0 Þ to M, i.e., gk ¼ g when g > g0 , and from K0 to the origin O along the triaxial compression stress path, gk changes from
where k and j are the slopes of normal consolidation and swelling lines, respectively, M is the slope of the critical state line, e0 is the initial void ratio of soil for p ¼ p0 , epv is the plastic volumetric strain that is used as a hardening parameter, and the relative stress ratio g is defined as [1] rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ðg gij0 Þðgij gij0 Þ g ¼ ð2Þ 2 ij
gij ¼
ð4Þ
Fig. 2. Contour lines of g1 in p-plane (K0 ¼ 0.5).
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
0:75ð¼ g0 Þ to 0, i.e., gk ¼ g when g < g0 . From the origin O to the point B along the triaxial extension stress path, gk changes from 0 to M, i.e., gk ¼ g. In other cases, gk changes from some values in the range of 0:75 0.75 to M when K0 ¼ 0:5. In addition, as is well known, the contour lines of the stress ratio gð¼ q=pÞ are concentric circles with the centers being the origin in the p-plane, and the contour lines of the relative stress ratio g are concentric circles with the centers being the K0 point shown in Fig. 2 in the p-plane. Hence, the stress ratio gk can be explained as an intermediate stress ratio between the relative stress ratio g and the stress ratio g.
3. An anisotropic hardening elastoplastic model for clays and sands As explained before, the Extended Mises criterion with the initial stress-anisotropy (g ¼ constant) is adopted as a shear yield criterion in the Sekiguchi–Ohta model. As is well known, the Extended Mises criterion is not appropriated for describing soil behavior and the SMP criterion is one of the best criteria for describing shear yield and failure behavior of soils in 3D stress. The SMP criterion has been successfully applied to the Cam~ij inclay model by introducing a transformed stress r stead of the conventional stress rij [6]. In the presented study, the model is generalized using the transformed stress based on the SMP criterion. The SMP criterion can be written by [4] I1 I2 ¼ constant; I3
ð6Þ
where I1 , I2 and I3 are the first, second and third stress invariants, respectively. The solid curve in Fig. 3 is the shape of the SMP criterion in the p-plane. The transformed stress is deduced from what makes the SMP curve in the p-plane become a circle with the center being the origin in the transformed p-plane, as shown in Fig. 3. From the geometrical relationship between two points A and A0 shown in Fig. 3, the transformed σ1 (σ1 ) Circle
A
SMP criterion
A'
lθ
θ (θ)
0
σ 2 (σ 2 )
stress can be calculated using the conventional stress as follows: ‘0 ~ij ¼ pdij þ pffiffiffiffiffiffiffiffiffiffi sij ; r skl skl
σ 3 (σ 3)
~-plane (dotted Fig. 3. The SMP criterion in p-plane (solid line) and p line).
ð7Þ
where sij is the deviatoric stress, and ‘0 is written by (Fig. 3) [6] rffiffiffi 2 I1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ‘0 ¼ 2 : ð8Þ 3 3 ðI1 I2 I3 Þ=ðI1 I2 9I3 Þ 1 If we take the place of the stress tensor rij with the ~ij , and introduce a new transformed stress tensor r hardening parameter H instead of the original hardening parameter epv to the Sekiguchi–Ohta model, the following yield function (f ) and plastic potential function (g) can be obtained ! ~p ~g kj f ¼g¼ ln þ H ¼ 0; ð9Þ ~p0 M 1 þ e0 ~ii =3) at the initial shear where p~0 is the value of ~pð¼ r state with the void ratio e0 , and M is explained as the value of q=p at the characteristic state at which the rate of plastic volume change is zero, i.e., depv ¼ 0, for dilative soil such as medium-to-dense sand and at the critical state for contractive soil such as normally consolidated clay in triaxial compression stress [7]. Although the value of M for sand varies with the confining pressure [8], M is usually assumed to be constant for the sake of simplicity and practical use. Moreover, ~ g is defined to be the same in the form as Eq. (2) rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ~g ¼ ð~g ~gij0 Þð~gij ~gij0 Þ ð10Þ 2 ij in which ~ij ~pdij r ~gij ¼ ; ~p
ð11Þ
~gij0 is the value of ~gij at the initial shear state. It can be seen from Eqs. (9) and (10) that the proposed yield surface is anisotropic because the initial stress ratio tensor ~gij0 is fixed. Following a hardening parameter for the isotropic hardening model proposed by authors [9], the increment of the hardening parameter H for the anisotropic hardening is defined as dH ¼
lθ
39
M 4 Mf4 ~g4k p de ; Mf4 M 4 ~g4k v
ð12Þ
where Mf is the value of q=p at failure state (epd ! 1) in triaxial compression stress, and is assumed to be constant for the sake of simplicity and practical use. ~ gk is defined in the same form as Eq. (4). ~gk ¼
3 ð~g ~gij0 Þ~gij : 2~g ij
ð13Þ
40
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
Fig. 4. Proposed stress–dilatancy relation (K0 ¼ 0:5).
Fig. 7. Yield loci in p-plane determined by undrained cyclic loading true triaxial tests on sand (after [11]).
Fig. 5. Proposed yield curve (fine solid line), failure state lines (bold solid line) and characteristic lines (dotted line) in triaxial stress.
By using ~ gk the stress–dilatancy relation in the proposed model is written as follows: ~ gk ¼ M
depv : depd
ð14Þ
Fig. 4 shows the shape of Eq. (14) when K0 ¼ 0:5. The hardening parameter H defined by Eq. (12) is a combination of both plastic volumetric strain epv and plastic shear strain epd , a concept first introduced by
Nova and Wood [10], after taking Eq. (14) into consideration. Since H is a hardening parameter, dH is not less than zero. Taking this into account, the following features of H can be deduced from Eq. (12): (1) When ~ gk ¼ 0 or Mf ¼ M, dH ¼ depv , i.e., the hardening parameter H becomes the plastic volumetric strain epv , which is the same hardening parameter as used in the Sekiguchi– Ohta model and Cam-clay model. (2) When 0 < ~ gk < M, depv > 0; when ~gk ¼ M, depv ¼ 0; and when M < ~ gk < Mf , depv < 0, i.e., the negative and positive dilatancy of soils can be described. In the proposed model, the associated flow rule is ~ij space [6], i.e., adopted in the transformed stress r depij ¼ K
of ; o~ rij
Fig. 6. Yield loci of the proposed model in deviatoric planes.
ð15Þ
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
where K is the proportionality constant, and can be determined from the consistency condition. Fig. 5 shows the proposed yield curve (Eq. (9), fine solid curve), failure state lines (~ gk ¼ Mf , bold solid line) and characteristic lines ~ g ¼ M, dotted line) in triaxial compression and triaxial extension stresses. From Eq. (12), it is known that when Mf ¼ M the failure state line becomes the same as the characteristic state line, and the hardening parameter H becomes the plastic volumetric strain epv , which is widely used as the hardening parameter in the elastoplastic model such as the Cam-clay model and the Sekiguchi–Ohta model for normally consolidated clay. Fig. 6 shows the yield loci of the proposed model in the deviatoric planes. The yield loci in the transformed p-plane (~ p-plane) are concentric circles with the centers being the point K0 as shown in Fig. 6(a), which are the same in the shape as those of the Sekiguchi–Ohta model in the p-plane, while the yield loci become oval-shaped
41
loci in the ordinary p-plane as shown in Fig. 6(b). It is interesting to note that in the p-plane the shape of the yield loci (Fig. 6(b)) is similar to that obtained using a true triaxial test apparatus by Yamada [11], as shown in Fig. 7. The yield points marked with and are determined using the method similar to those by Poorooshasb [12] and Tatsuoka [13]. All specimens were first subjected to a deviatoric stress along y-direction (i.e. h ¼ 0°) from O to A under Constant mean stress p. This stress path corresponds to that in conventional triaxial compression test. After completion of the first cycle along triaxial stress path the specimens were subjected to a loading along the stress path of h ¼ 15° or 30° or 45° to determine the yield stress points marked with from the points where the curvature of shear stress versus shear strain curve is maximum. The method for determining the yield points marked with is similar to the above, but the first experienced deviatoric stress was larger. The different yield loci shown in Fig. 7 are the
5 σ r /σ
a
4 3 r
Comp. R=σ a/σ r=1.5
1 1.5
-10
εa (%)
σ a-σ
2
Proposed S-O model Observed
1
0
-5
0
5
εr (%)
10
p
εv (%)
Ext. 5
Failure
R=1.5
(a) 5 a
4 r
σ r /σ
2 1 2.0
Proposed S-O model Observed
1
-10 ε (%) -5 a
0
σ a-σ
3
Comp. R=2.0
0 5
p
εr (%) 10
Ext.
εv (%)
Failure
5
(b)
R=2.0 5 a
4 r
σ r /σ
2 1 2.5
-10
εa (%)
-5
Proposed S-O model Observed
1 0
σ a-σ
3
5
ε r (%)
0 p
10
Ext.
εv (%)
Failure
5
(c)
Comp. R=2.5
R=2.5
Fig. 8. Comparison between predicted and experimental stress–strain behavior during shearing from anisotropic stresses (experimental data after [14]).
42
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
observed results of sand samples subjected to different values of deviatoric stress in the y-direction. It can be seen from Fig. 6(b) that the yield loci in real stress space are partly non-convex. This is because the convexity of the yield surface is satisfied in the transformed stress space, and the associated flow rule is adopted also in the transformed stress space. That is to say, in the proposed model DruckerÕs postulate is satisfied in the transformed stress space but not in real stress space.
4. Comparisons of model prediction with experimental result The model contains six parameters k, j, M, Mf , m and e0 . In comparison with the Sekiguchi–Ohta model, only one model parameter, Mf , is increased in the proposed model. However, for the soils exhibiting only negative dilatancy, i.e., contractive soils, the model parameters in two models are the same because Mf ¼ M in this case. Fig. 8 shows the comparisons of the predicted results by the Sekiguchi–Ohta model (abbreviated as S–O model in the figure; dotted line) and the proposed model (solid line) with the results of triaxial tests (mark ) on a contractive soil composed of Toyoura sand and Fujinomori clay (experimental data from [14]). Fujinomori clay is a silty clay and has been used as a normally consolidated clay in other soil tests (e.g. [15]). After Toyoura sand was mixed with Fujinomori clay with a ratio of 7:3 by dry weight, the mixture was consolidated up to 49 kPa for triaxial tests. In Fig. 8, ra and rr are the axial and radial stresses, and ea and er are the axial and radial strains in
σ a /σ
triaxial test, respectively. The material parameters used in the model prediction are as follows: k=ð1 þ e0 Þ ¼ 0:021, j=ð1 þ e0 Þ ¼ 0:0034, Mf ¼ M ¼ 1:32 and m ¼ 0:3. The triaxial specimens were first anisotropically consolidated up to p ¼ 980 kPa along the stress paths of ra =rr ¼ 1:5, 2.0 and 2.5, respectively, in triaxial compression and then sheared to triaxial extension up to failure under p ¼ 980 kPa as shown on the right-hand side of Fig. 8. The lefthand side of Fig. 8 shows the predictions and triaxial test results during shearing. It can be seen from Fig. 8 that the proposed model can predict the test results better than the Sekiguchi–Ohta model. Fig. 9 shows the experimental and predicted results of triaxial compression and triaxial extension tests on Toyoura sand from K0 ( ¼ 0.5)-consolidation state (experimental data from [16]). The stress paths tested are from point P at ra =rr ¼ 2 in triaxial compression to triaxial compression or triaxial extension up to failure under p ¼ 196 kPa or rr ¼ 147 kPa, as shown in Fig. 10. The material parameters used in the model prediction are k=ð1 þ e0 Þ ¼ 0:0508, j=ð1 þ e0 Þ ¼ 0:0112, M ¼ 0:95, Mf ¼ 1:66 and m ¼ 0:3. It can be seen from Fig. 9 that the model can properly predict the sand behavior including dilatancy characteristics. Figs. 8 and 9 demonstrate that the proposed model can predict well the stress–strain relation of not only clays but also sands during shearing from the initially anisotropic stress to triaxial compression stress or triaxial extension stress. Therefore, the proposed model is capable of predicting well the stress–strain relations of contractive and dilative soils shearing from an anisotropic stress state to any stress state in 3D stress.
6
σ r /σ
r
6 a
4
2 1 0
Path PA -10
ε r (%)
-5
4
Path PC 5
ε v (%) 4 σ a /σ
ε a (%)
ε a (%)
10 -5
1 2.0
6
σ r /σ
4
-10
ε r (%)
-5 ε v (%)
2 1 0 4
ε r (%)
0 2
ε v (%)
r
Path PB
1
55
6 a
4
Path PD 5
ε a (%)
ε a (%)
10 -5
1 2.0
1 0 ε (%) v 2
ε r (%) 55
Fig. 9. Comparison between predicted and experimental stress–strain behavior of Toyoura sand during shearing from an anisotropic stress (experimental data after [16]).
800
Dijkl ¼ Deijkl Deijmn
43
of of e D =X : o~ rmn orst stkl
ð21Þ
Introducing the elastic constitutive tensor of HookeÕs law for isotropic elasticity into Eqs. (19) and (21), we can obtain X ¼
R= 1
600
A R= 2
failure line
σa (kPa)
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
B 400
k j M3 Mf4 g~4k of of þL 4 2 2 1 þ e0 Mf ðM þ ~gk ÞðM þ ~gk Þ~p orii o~ rjj of of þ 2G ; orij o~ rij
ð22Þ
P
Dijkl ¼ Ldij dkl þ G dik djl þ dil djk 200
D
C
of of dij þ 2G L o~ rmm o~ rij
line failure
!
of of L dkl þ 2G ornn orkl
, X; ð23Þ
0
200
400
600
where L and G are LameÕs constants, which can be written from e ¼ e0 j ln p=p0 in elastic region by
σr (kPa) Fig. 10. Stress paths for Toyoura sand.
5. Elastoplastic constitutive tensor The elastoplastic constitutive tensor for the proposed model to solve elastoplastic boundary value problem by means of the finite element method is derived here. The elastic part of the stress–strain relation can be written in the incremental form as drij ¼ Deijkl deekl ¼ Deijkl ðdekl depkl Þ;
ð16Þ
where Deijkl is the elastic constitutive tensor. ekl and epkl are the total strain tensor and its plastic component, respectively. Because the yield function f is the function of the stress tensor and the hardening parameter, the consistency condition can be written by df ¼
of of dH ¼ 0: drij þ orij oH
ð17Þ
Substituting Eqs. 12,15 and 16 into Eq. (17) gives K¼
of orij
Deijkl dekl X
;
ð18Þ
where X ¼
M 4 Mf4 ~ g4k of of e of dij þ D : 4 4 4 rij orij ijkl o~ rkl Mf M ~ gk o~
ð19Þ
Substituting Eqs. (15) and (18) into Eq. (16), we get a general form of the proposed model drij ¼ Dijkl dekl ;
ð20Þ
where the elastoplastic constitutive tensor is expressed as
G¼
E 3ð1 2mÞð1 þ e0 Þ ¼ p; 2ð 1 þ m Þ 2ð1 þ mÞj
ð24Þ
L¼
E 2 ð1 þ e0 Þ 2 G¼ p G; 3ð1 2mÞ 3 j 3
ð25Þ
where E and m are the elastic modulus and PoissonÕs ratio, respectively.
6. Finite element computation The proposed model is incorporated into the FEM program DACSAR coded by Iizuka and Ohta in Japan [17]. As shown in Figs. 2 and 5, when the stress path is along the K0 -line, the plastic flow at the singular point of non-smooth yield surface is unable to be correctly evaluated due to an unpredictable gradient at the corner. A very simple handling method adopted in the FEM code is that K0 is assumed to be slightly smaller than the stress ratio (r2 =r1 ¼ r3 =r1 ) when the stress state is just on the K0 -line. For the singularity of the Sekiguchi–Ohta model, more precious numerical method has been proposed by Yuttapontada et al. [18] based on KoiterÕs associated flow rule. The governing equations describing the behavior of the soil skeleton coupled with the pore water flow are formulated by Biot [19]. The technique employed in discretizing the pore water flow was originally proposed by Akai and Tamura [20]. As an example of the application of the proposed model to practical engineering, as shown in Fig. 11, a natural ground under embankment loading [17] is analyzed using finite element method. The use of the earth embankment is for carrying a highway. The ground
44
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
36. 00m 10. 26m
15. 48m
10.26m
GL-0.0m GL-1.4m GL-7.8m
AC1 AS
AC2 GL-23.3m Fig. 11. Soil profile of ground and embankment.
X fixed boundary
X fixed boundary
23.3m 10m
Y
X
X, Y fixed boundary
150m Fig. 12. Finite element meshes and boundary conditions.
2
q/p
Mf=1.79 M=1.45
εv (%)
0
5 ε1 (%)
k=ð1 þ e0 Þ j=ð1 þ e0 Þ M Mf m k (m/day)
1
1.30 10 1.30 102 1.40 1.40 0.30 2.0 104
Determination of M and Mf 0 Cc 1+e0
1+e0 As (sand)
0.02 (b)
=0.0217 Measured Predicted
0.01
Cs
9.44 103 4.94 103 1.45 1.79 0.30 1.0 102
10
1
(a)
Table 1 Material parameters used in FE computation Ac2 (clay)
p=392kPa Measured Predicted
1
εv
consists of Ac1 layer (surface clay layer: GL0 )1.4 m), As layer (sand layer: GL-1.4 )7.8 m) and Ac2 layer (uniform marine clay layer, GL-7.8 )23.3m). The ground table is at 0.5 m below the ground surface. The fill material is the dredged sand with a unit weight of 1.73 tf/m3 . Fig. 12 shows the finite element meshes and boundary conditions. Only half section is computed for the symmetry of earth structure. The ground surface and the bottom boundary were treated as drained boundaries, while the two sides of the ground were treated as undrained boundaries. The two sides were restrained in the horizontal direction, while the bottom boundary was fixed in both the vertical and horizontal directions. Table 1 shows the values of the material parameters used in the finite element computation. These values are determined from the laboratory tests on the samples taken from the site of the embankment. In detail, the
=0.0114
100
200
400 p (kPa) Determination of Cc and Cs
1000
Fig. 13. Determination of model parameters for the soil As (experimental data after [21]).
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
0
Predicted Observed
100cm
0 5 10
Fill height (m)
Fig. 15. Surface settlement just after embankment construction.
pw(kPa)
150 100
Depth 9.0m
Predicted Observed
50 0
pw(kPa)
model parameters for the soil Ac2 were determined by using Cc =ð1 þ e0 Þ ¼ 0:3, Cs =ð1 þ e0 Þ ¼ 0:03 and /0 ¼ 34:5°, which are obtained from results of one-dimension consolidation tests and drained triaxial tests on the soil Ac2 [21]. Here, Cc and Cs are compression index and swelling index, respectively; e0 is the initial void ratio and /0 is the effective angle of friction. The model parameters for the soil As are determined from Fig. 13, where the experimental data are from [21] and the solid lines show the model predictions. The value of Mf is determined based on /0 ¼ 43:5° for the soil As that was obtained from the envelope of three MohrÕs stress circles at failure [21]. The values of the coefficient of permeability for the soils As and Ac2 were taken from [17,21], respectively. The surface clay layer (Ac1) is assumed as an elastic material. The observed and computed settlements at the ground surface and depth of 7.8 m (i.e., the top of the soil Ac2) under the center of the embankment are plotted against time during and after embankment construction in Fig. 14. In all the figures, the plot is the observed results and the solid line is the results computed by the finite element method using the proposed
45
200 150 Depth 11.5m 100 50
pw(kPa)
5 10
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1. 6 0
Predicted Observed
At ground surface
100
200
300 Time (day)
0 Predicted Observed
Settlement (m)
0.2 0.4 0.6 0.8 At GL-7.8
(b)
0
200 150 Depth 21.5m 100 50 0
(a)
1
200 150 Depth 15.0m 100 50 0
pw(kPa)
Settlement (m)
Fill height (m)
0
100
200
300 Time (day)
Fig. 14. Loading process and settlements under the center of the embankment.
50
100
150
200
250 300 Time (day)
Fig. 16. Generation and dissipation of excess pore water pressure at different depths of the centerline.
model. Fig. 15 shows the observed and computed ground surface settlements at the time of just completing the embankment construction. The excess pore water pressures induced by embankment loading at different depths of central line are plotted against time during and after embankment construction in Fig. 16. From Figs. 14–16, it can be seen that the settlements and the excess water pressures computed using the proposed model fairly agree with the observed results. However, the computed excess pore water pressure at a depth of 21.5 m is lower than the observed value (Fig. 16). This is due to adopting the treatment of the bottom boundary as a drained boundary in the finite element computation while the true bottom boundary may be of partial drainage.
46
D.A. Sun et al. / Computers and Geotechnics 31 (2004) 37–46
7. Conclusions This paper presented an anisotropic hardening model for clays, silts and sands within a simple framework, and the finite element analyses for an embankment on the ground consisting of clays and sands using the proposed model. Several conclusions can be drawn from this study. 1. A unified and simple elastoplastic model for clays and sands was proposed, which is capable of describing deformation and strength characteristics of clays and sands with the initially anisotropic consolidation state in 3D stress. 2. For contractive soils, the model parameters involved in the proposed model are the same in the Cam-clay model or the Sekiguchi–Ohta model. For dilative soil, only one model parameter is added to the model parameters of the Cam-clay model or the Sekiguchi– Ohta model. 3. By using the given elastoplastic constitutive tensor, the proposed model can be easily applied to the finite element computation for earth structure composed of sands and clays. 4. The yield locus for soils subjected to the anisotropic stress should be oval-shaped in the p-plane. Acknowledgements The authors are grateful to Dr. S.H. Liu of Nagoya Institute of Technology and Dr. H.L. Fang of West Japan Engineering Consultants for their help in the revision of the manuscript, and to Prof. A. Iizuka of Kobe University for providing the program DACSAR. References [1] Sekiguchi H, Ohta H. Induced anisotropy and time dependency in clays. In: Murayama S, Schofield AN, editors. Proceedings of the 9th ICSMFE 1977; Specialty Session 9. p. 229–38. [2] Roscoe KH, Schofield AN, Thurairajoh A. Yielding of clay in state wetter than critical. Geotechnique 1963;13(3):211–40. [3] Duncan JM. The role of advanced constitutive relations in practical applications. In: Proceedings of the 13th International Conference on Soil Mechanics and Foundation Engineering, vol. 5; 1994. p. 31–48. [4] Matsuoka H, Nakai T. Stress-deformation and strength characteristics of soil under three different principal stresses. Proc JSCE 1974;232:59–74.
[5] Lade PV, Duncan JM. Elasto-plastic stress–strain theory for cohesionless soil. J Geotech Eng ASCE 1975;101(10):1037–53. [6] Matsuoka H, Yao YP, Sun DA. The Cam-clay models revised by the SMP criterion. Soils Foundations 1999;39(1):81–95. [7] Luong MP. Stress–strain aspects of cohesionless soils under cyclic and transient loading. In: International Symposium on Soils under Cyclic and Transient Loading; 1982. p. 315–24. [8] Lade PV, Ibsaen LB. A study of the phase transformation and the characteristic lines of sand behavior. In: Asaoka et al., editor. International Symposium Deformation and Progressive Failure in Geomechanics; 1997. p. 353–8. [9] Yao YP, Matsuoka H, Sun DA. A unified elastoplastic model for clay and sand with the SMP criterion. In: Vitharana N, Colman R, editors. Proceedings of 8th Australia–New Zealand Conference on Geomechanics, vol. 2; 1999. p. 997–1004. [10] Nova R, Wood DM. A constitutive model for sand in triaxial compression. Int J Numer Analytical Meth Geomech 1979;3(3): 255–78. [11] Yamada K. Deformation characteristics of loose sand under three-dimensional stress. Doctoral thesis at University of Tokyo; 1979 [in Japanese]. [12] Poorooshash HB, Holubec I, Sherbourne AN. Yielding and flow of sand in trixial compression. Can Geotech J 1967;4(4):376–98. [13] Tatsuoka F, Ishihara K. Yielding of sand in triaxial compression. Soils Foundations 1974;14(2):63–76. [14] Ichihara W. Deformation and strength characteristics of intermediate soil under triaxial stress condition. BachelorÕs Thesis at Nagoya Institute of Technology; 1997 [in Japanese]. [15] Nakai T, Matsuoka H, Okuno N, Tsuzuki K. True triaxial tests on normally consolidated clay and analysis of the observed shear behavior using elastoplastic constitutive models. Soils Foundations 1986;26(4):67–78. [16] Funada T, Matsuoka H, Fukumoto S. Deformation characteristics of sand under various stress paths from anisotropic consolidation. In: Proceedings of the 44th Annual Conference of the Japan Society of Civil Engineers, vol. 3; 1989. p. 512–3 [in Japanese]. [17] Iizuka A, Ohta H. A determination procedure of input parameters in elasto-viscoplastic finite element analysis. Soils Foundations 1987;27(3):71–87. [18] Yuttapontada Y, Pipstpongsa T, Iizuka A, Ohta H. FEM for the Sekiguchi–Ohta model consider non-smooth yield surface. In: Proceedings of the 38th Japan National Conference on Geotechnical Engineering, 2003. p. 285-6. [19] Biot MA. General theory of three-dimensional consolidation. J Appl Phys 1941;12:155–64. [20] Akai K, Tamura T. Numerical analysis of multi-dimensional consolidation accompanied with elasto-plastic constitutive equation. Proc JSCE 1978;269:95–104 [in Japanese]. [21] Matsuoka H, Nakai T. Prediction using constitutive equation for soils based on the concept of Spatially Mobilized Plane. Research Report on settlement prediction of ground due to embankment edited by the committee on investigating highway, the Japan Highway Public Corporation 1984. p. 220–42 [in Japanese].