Formulation of implicit stress integration and consistent tangent modulus for an anisotropic hardening constitutive model

Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

www.elsevier.com/locate/cma

Formulation of implicit stress integration and consistent tangent modulus for an anisotropic hardening constitutive model Seboong Oh a

a,*

, Seung R. Lee

b

Department of Civil Engineering, Yeungnam University, Kyongsan 712-749, Republic of Korea b Department of Civil Engineering, KAIST, Taejon, Republic of Korea Received 15 April 1999; received in revised form 21 March 2001

Abstract An anisotropic hardening constitutive model based on the generalized isotropic hardening (GIH) rule was implemented into the nonlinear ®nite element method. The GIH rule describes discrete formation of homology centers and simultaneous hardening of both inner and outer yield surfaces. Consequently, yielding in the ®eld of reverse loading can be accurately modeled by this rule. In order to preserve the accuracy of nonlinear ®nite element analysis, an implicit stress integration technique employing the generalized trapezoidal rule was applied to the GIH constitutive equation. Furthermore, the consistent tangent modulus was formulated and coded into a nonlinear FEA program, which has the advantage of the quadratic rate of convergence in Newton's iteration. Accuracy and convergence were successfully veri®ed through examples, and thus the anisotropic hardening constitutive model could assess the plastic straining mobilized in overconsolidated state and versatile loading sequences. Ó 2001 Elsevier Science B.V. All rights reserved.

1. Introduction The nonlinear ®nite element solution by Newton type method is widely used to predict responses in geotechnical problems. Such a linearized scheme requires integrating the constitutive equations in order to assure the accuracy and eciency of numerical solutions. The implicit integration method has been regarded as the most appropriate scheme to hold the accuracy of resultant stress, but in order to conserve the quadratic rate of convergence, the tangential stress±strain modulus should be consistent with the linearization [1±6]. Soils undergo plastic strains even under reverse loading conditions, which may be induced by sedimentation history or by loading sequences. Such actual constitutive behavior can be successfully captured by anisotropic hardening models in the context of elasto-plasticity [7±10]. Most anisotropic hardening descriptions are related to the isotropic-kinematic combined hardening rule in which hardening laws have been formulated by adopting an interpolation rule. The authors have also proposed an anisotropic hardening description based on the generalized isotropic hardening (GIH) rule [11]. The anisotropic description of GIH postulates discrete formation of the center of isotropically hardening homology at the stress state of reverse loading events, and simultaneous occurrence of GIH for both inner and outer yield surfaces. Because the GIH model can represent plastic deformations under reverse loading conditions, the plastic behavior of sedimentary soils can be accurately re¯ected in the numerical analysis. *

Corresponding author.

0045-7825/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 1 ) 0 0 2 7 4 - 2

256

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

However, the rigorous anisotropic hardening constitutive models contain multiple hardening variables. Then, robust implementation into a numerical procedure results in a complex mathematical formulation. Therefore, most of the successful anisotropic models used for geomaterials have not yet been facilitated into a numerical analysis such as the ®nite element method. In this study, an implicit stress integration of the GIH model is formulated based on the generalized trapezoidal rule [1]. A tangent modulus is also formulated from the stress integration procedure and coded into a nonlinear analysis program, which is consistent with the linearization of the corresponding nonlinear governing equation [2]. Accuracy and convergence rate are veri®ed by numerical examples. The feasibility of robust implementation for complex anisotropic hardening elasto-plastic constitutive models is presented in this paper.

2. Implicit stress integration by generalized trapezoidal rule The constitutive equation based on the elasto-plasticity is generally of rate form, and its integration is initiated from the following equation: dr ˆ Ce : dee ;

dee ˆ de

dep ;

…2:1†

where Ce is an elastic stress±strain matrix. According to the associated ¯ow rule, the plastic strains are de®ned by dep ˆ d/f;r :

…2:2†

In the above equation d/ is a consistency parameter, and f;r  of =or. Let rn be converged stresses in the previous time step tn . The current stress tensor of kth iteration, rkn‡1 , corresponding to a given strain tensor ekn‡1 ˆ en ‡ Dekn‡1 yields [3] Z tn‡1 rkn‡1 ˆ rn ‡ Ce : dee : …2:3† tn

It is assumed that the elastic modulus Ce is a constant matrix for a given interval, and hence a prescribed trial stress tensor can be de®ned as rtr ˆ rn ‡ Ce : De;

…2:4†

where De  Dekn‡1 . Because the plastic strains involved in elasto-plastic constitutive equations cannot be commonly integrated in a closed form, it is necessary to use a numerical integration scheme. For a closed interval ‰tn ; tn‡1 Š, the integration of (2.2) can be approximated by the generalized trapezoidal rule as follows [1]: Z tn‡1 f;r d/  D/f…1 d†…f;r †n ‡ d…f;r †n‡1 g; …2:5† tn

where d is a trapezoidal parameter, and …f;r †n and …f;r †n‡1 are derivatives of yield functions at each time step. Remark 2.1. If d ˆ 0, the stresses are integrated by an explicit forward Euler method. For the backward Euler method …d ˆ 1†, the closest point projection algorithm is obtained [6]. Oritz and Popov [1] have also reported that stress integrations with d ˆ 1=2 lead to the second-order accuracy. On the other hand, when large strain increments are anticipated, the closest point projection algorithm is optimal. The optimal choice of d depends on the nature of the problem user consideration. Substituting (2.4) and (2.5) into (2.3) results in rkn‡1 ˆ rtr

Ce : …f;r †n‡d ;

…2:6a†

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

…f;r †n‡d ˆ …1

d†…f;r †n ‡ d…f;r †n‡1 :

257

…2:6b†

As the derivative of yield surface in (2.6b) is dependent on the unknown rkn‡1 , the stress integration of (2.6a) should be implicitly performed by an iterative procedure. In order to maintain the quadratic convergence of Newton's method, the tangent modulus must be consistent with the numerical method employed to integrate the constitutive equations. Consistency implies that the stress increment calculated by the tangent modulus operating on the strain increment matches the stress increment calculated by the integration procedure to ®rst-order [3]. The consistent tangent modulus is then derived from the differentiation of stress integration functions (2.6a) and (2.6b) [2] Ckn‡1 ˆ

orkn‡1 : oekn‡1

…2:7†

3. Anisotropic hardening constitutive model based on GIH rule On the basis of the GIH rule, an anisotropic hardening model for reverse loading conditions was proposed in which the center of isotropically hardening homology is selected at the instantaneous stress state of a reverse loading event [11]. The GIH model requires two types of surfaces; e.g., a reference surface and a yield surface, as shown in Fig. 1. Their equations are respectively F ˆ … p

spc †2 ‡

f ˆ …p

a† ‡

2

3 s : s 2 m2

3 …s 2

d 2 pc2 ˆ 0;

b† : …s m2

b†

…3:1a† r2 ˆ 0:

…3:1b†

In the above equations, tensors of stress, strain and center stress are decomposed into volumetric and deviatoric terms as follows: r ˆ p1 ‡ s;

e ˆ 13em 1 ‡ e;

ac ˆ a1 ‡ b:

…3:2†

In the context of soil mechanics, stress terms are positive in compression, and describe e€ective stresses which apply between particle contacts (refer to Refs. [12±14] for tensor notations used in this paper). the variable r is the radius of the yield surface along the p-axis, and pc is the preconsolidation pressure  is a mapping stress onto the reference surface, and which implies the size of reference surface. Further r ac is its center. The yield surface is composed of two half ellipses and constants appearing in (3.1a) and (3.1b) are de®ned as follows, using parameters M; s and c [11]. dˆ1

s;

m ˆ Ms=…1

s† if p P a;

Fig. 1. Anisotropic hardening constitutive model based on GIH rule.

…3:3a†

258

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

d ˆ c ‡ s;

m ˆ Ms=…c ‡ s† if p < a:

…3:3b†

Both the yield surface and the reference surface harden isotropically according to their own hardening rules [11]. The hardening function of the reference surface is described as in the Cam-clay model [14,15], 1‡e pc depm ; k j

dpc ˆ

…3:4†

where e is the void ratio, and k and j are compression and swelling indices, respectively. Assuming that the deviatoric plastic strain is merely dependent on the radius of the yield surface during reverse loading, the hardening function of the yield surface can be de®ned as epd ˆ

Z

depd ˆ

b

a …cr=dpc † ; 2 m h2 1 cr=dpc

…3:5†

where the interval of integration is governed by a reverse loading process. In Eq. (3.5), depd ˆ

p p 2=3kdeij

depkk =3k;

h ˆ

q 2 h2ii =9 ‡ 2…hij hii =3† =3

and h is the center of homology in the stress ®eld normalized by pc . In addition, three parameters a, b and c are needed to describe the anisotropic hardening. Remark 3.1. As shown in Fig. 2(a), the reference surface hardens isotropically during the K0 consolidation process where a deviator stress, q, equals …r1 r3 † under conventional triaxial test conditions. The center of homology h is immediately created at the stress point of abrupt change in the loading direction, such as point A of `swelling' (Fig. 2(b)) and point B0 of `compression after swelling' (Fig. 2(c)). In Fig. 2(c), the stress path directs towards the inside of the previous yield surface, yet a new h is formed within the reference surface. Thereafter yield surfaces harden isotropically respect to h.

Fig. 2. Concept of anisotropic hardening model under triaxial loadings: (a) K0 consolidation, (b) swelling and extension and (c) compression after swelling.

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

259

Remark 3.2. The critical state soil mechanics is the most important background theory of this model, but there are substantial di€erences between the modi®ed Cam-clay model [15] and the GIH model. The modi®ed Cam-clay model has two crucial defects: (1) the undrained shear strength even in the normally consolidated clay is overpredicted and (2) the yielding within the outer yield surface (or reference surface) cannot be modeled which is observed in the overconsolidated or the reversely stressed clays. In the GIH model the shape of reference surface was modi®ed by using the parameter s (see Eqs. (3.3a) and (3.3b)) which controls the shape to predict better the undrained shear strength for the same M. In particular, the anisotropic hardening description can catch the yielding within the reference surface, which is corresponding to the yield surface of modi®ed Cam-clay model. Note that the isotropic hardening Cam-clay model predicts just the elastic response within the reference surface for the behavior of both Figs. 2(b) and (c). Therefore the GIH model has an improved ability to accurately predict the undrained shear strength of normally consolidated clays and to reasonably calculate the plastic strain within the reference surface. Remark 3.3. The nesting yield surface model [7] and the bounding surface model [8] applied a mapping rule (or an interpolation rule) using the mapping stress and its distance from the current stress or any equivalent variable in order to determine the plastic modulus. On the other hand, in the GIH rule, it is not necessary to use the mapping rule for obtaining the plastic modulus, since the size of the yield surface is the only required variable for hardening functions. The discrete formation of h can be regarded as a discontinuous kinematic hardening rule. In virtue of the discontinuous change of h and the GIH concept, the combined hardening phenomena can be simply accounted for without resorting on any complicated kinematic hardening functions.

4. Mathematical formulations of GIH model 4.1. Formulation of stress integrations The elastic bulk and shear moduli, K and G, are linearly dependent on p as in the Cam-clay model according to Kˆ

1‡e p; j

Gˆ

3K…1 2m† ; 2…1 ‡ m†

…4:1†

where m is Poisson's ratio. For simplicity, it is assumed that K and G are constants over ‰tn ; tn‡1 Š such that K ˆ K…tn †. The trial stress in (2.4) is then decomposed as ptr  13tr…rtr † ˆ pn ‡ KDem ;

…4:2a†

str  rtr

…4:2b†

1 tr…rtr †1 3

ˆ sn ‡ 2GDe:

We can now derive integrated stresses from (2.6a) k p  pn‡1 ˆ ptr

KD/…f;p †n‡d ;

…4:3a†

s  skn‡1 ˆ str

2GD/…f;s †n‡d ;

…4:3b†

where the stress derivatives of the yield function are de®ned as …f;p †n‡d ˆ …1

d†…f;p †n ‡ d…f;p †n‡1 ;

…4:4a†

…f;s †n‡d ˆ …1

d†…f;s †n ‡ d…f;s †n‡1 :

…4:4b†

260

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

The integrations of hardening functions are respectively k

pc  …pc †n‡1 ˆ pcn expfKD/…f;p †n‡d g; r 2 p ke k ˆ w…r; pc †; 3

…4:5a† …4:5b†

k b where K ˆ …1 ‡ e†=…k j†, pcn  pc …tn †; ep  …ep †n‡1 ˆ …ep †n ‡ D/…f;s †n‡d and w ˆ a…cr=dpc † =fm2 h2 …1 cr= dpc †g. Therefore (4.3a) and (4.3b) are to be solved implicitly with (4.5a) and (4.5b). Let a residual be de®ned as r ˆ fr1 rT2 r3 r4 r5 gT for a variable vector u ˆ fp sT pc r D/gT such that

r1 ˆ p

ptr ‡ KD/…f;p †n‡d ;

…4:6a†

r2 ˆ s

str ‡ 2GD/…f;s †n‡d ;

…4:6b†

r3 ˆ pc pcn expfhD/…f;p †n‡d g; r 2 p r4 ˆ ke k w…r; pc †; 3 r5 ˆ …p

2

a† ‡

3 …s 2

b† : …s m2

…4:6c† …4:6d† b†

r2 :

…4:6e†

The stress integration is then a problem for solving the following nonlinear system. r…un‡1 † ˆ 0;

…4:7†

where n ‡ 1 is the current step number. Thus a linearization is required and iterations are performed over r0 …uin‡1 †Dui ˆ ri ;

i ui‡1 n‡1 ˆ un‡1

Dui ;

…4:8†

where i is the iteration count under the stress integration. Refer to Appendix A for details of Jacobians r0 in (4.8). The solution procedure is described in Table 1. Remark 4.1. In anisotropic hardening constitutive models, there are usually multiple internal variables. This formulation has two equations of (4.5a) and (4.5b) for the hardening description. In addition, this formulation excludes the loading on reference surfaces since it is similar to a Cam-clay type model. We see such formulations from many sources [4,10,13,14], and hence it is not addressed in this paper. 4.2. Formulation of consistent tangent moduli In order to preserve the quadratic rate of convergence, tangent moduli should be obtained from the stress integration of (2.3). It can be derived from (2.7) [2,13,14], Ckn‡1 ˆ 1

k opn‡1 oskn‡1 ‡ : oekn‡1 oekn‡1

Table 1 Procedure for solving r ˆ 0 1. 2. 3. 4. 5.

Initialize i ˆ 0 Compute ri ˆ r…uin‡1 † If kri k=kr0 k < Tolerance, return; else i Update ui‡1 Dui such that r0 …uin‡1 †Dui ˆ ri n‡1 ˆ un‡1 i i ‡ 1 and go to 2

…4:9†

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

261

Let …†0  o…†kn‡1 =oekn‡1 . The derivatives in (4.9) are then implicitly evaluated from (4.3a) and (4.3b), p0 ˆ …ptr †0 tr 0

s0 ˆ …s †

KD/0 …f;p †n‡d

KD/…f;p †0n‡d ;

2G…f;s †n‡d D/0

0

0

where …ptr † ˆ K1, …str † ˆ 2G…I …f;p †0n‡d ˆ 2d…p0

1 1 3

…4:10b†

1†, and

…f;s †0n‡d ˆ 3d…s0

a0 †;

…4:10a†

2GD/…f;s †0n‡d ;

b0 †=m2 :

…4:11†

The derivatives of center stress tensors in (4.11) are described in Appendix A. Hardening functions (4.5a) and (4.5b) are di€erentiated with respect to ekn‡1 as 0

pc0 ˆ KD/pc …f;p †n‡d ‡ Kpc …f;p †n‡d D/0 ; r o 2 ep n 0 0 : D/ …f;s †n‡d ‡ D/…f;s †n‡d ˆ w;r r0 ‡ w;pc pc0 ; p 3 ke k

…4:12† …4:13†

where w;pc and w;r are also described in Appendix A. Finally, di€erentiating the yield function (3.1b) with respect to ekn‡1 results in f;p …p0

a0 † ‡ f;s : …s0

b0 † ‡ f;r r0 ˆ 0:

…4:14†

Eqs. (4.10a), (4.10b), (4.12), (4.13) and (4.14) include ®ve independent variables of p0 ; s0 ; pc0 ; D/0 and r0 . Thus, we can ®nd p0 and s0 solving these equations, and determine the consistent tangent moduli of (4.9). Details are described in Appendix B. 5. Numerical examples 5.1. Example 1: Conventional triaxial tests for overconsolidated clays In order to perform an accuracy analysis, stress±strain relationships at a point were examined for drained and undrained triaxial compression tests. Undrained triaxial loadings are evaluated by the algorithms presented in Table 1 with volume constraints such that e2 ˆ e3 ˆ e1 =2. For drained traxial tests, a boundary value problem is required to solve the inverse relation en‡1 ˆ e…rn‡1 †. Accordingly, let rn‡1 be applied stresses at a point at tn‡1 and rn‡1 be statically equivalent internal stress. In the drained condition of triaxial loadings the components of rn‡1 are de®ned as r1 ˆ rc ‡ Dr and r2 ˆ r3 ˆ rc where rc is the con®ning cell pressure and Dr (or q) is the principle stress difference. To satisfy the equilibrium equation, rn‡1 ˆ r…en‡1 †;

…5:1†

where en‡1 are the unknowns. A linearization is possible to solve en‡1 by Newton's method such that C kn‡1 Dek ˆ rn‡1

rkn‡1 ;

k k ek‡1 n‡1 ˆ en‡1 ‡ De :

…5:2†

The iteration is performed until a residual k…rn‡1 rkn‡1 †=…rn‡1 rn †k, is less than a tolerance e.g., 1:0e 5. In the tests on Weald clays the cylindrical specimen was loaded to 120 psi by pressure and subsequently unloaded until the con®ning stress of 5 psi. The overconsolidation ratio (OCR) was 24 [16]; the OCR means the ratio of previous maximum pressure to the con®ning stress. In this case, even though the con®ning stress was within the reference surface, i.e., an outer yield surface, large plastic strains actually occurred as the axial stress increased. Table 2 shows the parameter values [11]. The drained stress±strain relationships are shown in Fig. 3, in which axial stress increments are 0.1, 0.5, 1.0 and 2.0 psi. The stresses are very accurately integrated in comparison with the reference case of 0.1 psi, although some errors are shown for q of 8 psi, which approaches peak failure. At this point volumetric strains are also inaccurately calculated. The larger increment induces the smaller strains.

262

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

Table 2 Values of material parameters Parameters

k

j

m

ea

M

sa

ca

ab

b

cc

Values

0.093

0.035

0.25

0.908

0.95

0.45

0.1

0.0004

1.8

1.0

a

For modi®ed Cam-clay, s ˆ 0:5, c ˆ 0:0. Parameters a and b are used for GIH description. c Parameters c affects the yield surface to meet the reference surface before critical state failure if c < 1:0 [11]. b

Fig. 3. Drained triaxial tests of Weald clay: (a) stress±strain relationships and (b) volume change.

Fig. 4. Undrained triaxial tests of Weald clay: (a) stress±strain relationships and (b) e€ective stress paths.

However, other strains do show very accurate volume changes. Further, the measured stress±strain relationships and volume changes are accurately predicted using GIH model, since some plastic strains in overconsolidated states can be appropriately modeled. The undrained stress±strain relationships and e€ective stress paths are shown in Fig. 4. The axial strain increments are 0.1%, 0.5%, 1.0% and 2.0%. Most integrated stresses show very accurate results without visible errors. Fig. 4 shows that only the larger increments, especially at 2%, induces the smaller integrated stresses. From an engineering point of view, it can be concluded that the stress integration algorithm results in an acceptable accuracy.

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

263

For the undrained behavior of Weald clay represented in Fig. 4, the same parameters and initial conditions were used as for the drained tests. Nevertheless, the stress±strain relationships are accurately predicted as shown in Fig. 4(a). In Fig. 4(b), the predicted e€ective stress paths also agree well with the measured behavior. Excessive strains will induce failure at the critical state line (CSL). Consequently, the GIH model consistently predicts both the drained and undrained behaviors of clay. For the drained samples, strains at 4 psi deviatoric stress are summarized in Table 3, showing nearly identical values for all stress increments. Accuracy refers to the rate at which a solution stabilizes as the applied load or the prescribed displacement is divided into subincrements. As stress was divided into subincrements, the solutions stabilized, as shown in Table 3. It can thus be concluded that numerical integrations were accurately performed for all values of d. stress integrations with d ˆ 1=2 led to slightly more accurate solutions than the closest point projection method, although the di€erences were insigni®cant. However, 2 psi incremental stress case with d ˆ 1=2 did not converge at the 8 psi solution, which approaches the failure conditions. This re¯ects the fact that the trapezoidal rule is unconditionally stable only for for d ˆ 1 [1]. The rate of convergence was so rapid that the solutions could converge within 3±4 iterations. If the error is Ek ˆ 10 m in a step, then for the second-order, Ek‡1 ˆ const
10 2m or Ek‡1 ˆ const
Ek2 , such that the number of signi®cant digits is almost doubled in each step. As represented in Fig. 5, the force error norms for 2 psi stress increment show that the convergence rate is asymptotically quadratic. For undrained samples, stresses at 10% axial strains are summarized in Table 4. The solutions also stabilized as the stress was divided into subincrements for any value of d. Stress integration with d ˆ 1=2 led to more accurate solutions than the closest point projection method. Table 3 Strains at q of 4 psi for drained tests Stress increment (psi)

0.1 0.5 1 2 a

d ˆ 0:5

d ˆ 1:0

ea (%)

em (%)

ea (%)

em (%)

1.36 1.37 1.39 1.41a

0.18 0.18 0.19 0.21a

1.36 1.38 1.39 1.42

0.18 0.18 0.19 0.21

Solution diverges at 8 psi.

Fig. 5. Convergence rate for 2 psi stress increment (d ˆ 1:0).

264

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

Table 4 Stresses at axial strain of 10% for undrained tests Stress increment (%)

0.1 0.5 1 2

d ˆ 0:5

d ˆ 1:0

p (psi)

q (psi)

p (psi)

q (psi)

11.92 11.82 11.70 11.49

10.37 10.30 10.20 10.03

11.91 11.77 11.61 11.29

10.38 10.31 10.25 10.14

5.2. Example 2: FE analysis of vertically loaded foundation The implicit stress integration and consistent tangent modulus were coded into a nonlinear FEA program, which solves governing equations by Newton's method. In order to verify the accuracy and eciency of those algorithms, an embedded rigid structure in Weald clay, which was loaded vertically under the plane stain condition was analyzed. The nonlinear FE analysis was performed for a mesh consisting of 4-node quadrilateral elements, as shown in Fig. 6. The convergence criteria for residual forces and energies are de®ned as [13] kRk k=kR0 k 6 eR ;

…5:3a†

jDdk
Rk j=jDd0
R0 j 6 eE ;

…5:3b†

where Rk ˆ …Fext †n‡1 Fint …rkn‡1 †, Fext is external force vector; Fint is internal force vector; and Ddk is search direction. The symbol k k denotes an L2 -vector norm and j j absolute values. Error tolerances were selected as eR ˆ 10 5 and eE ˆ 10 8 . All computations were performed in double precision.

Fig. 6. Finite element mesh of foundation problem.

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

265

Generally, clays have low permeability, and loading situations may rapidly result in insucient drainage; thus an undrained analysis is considered as the more signi®cant situation to design. The undrained problem can be solved by the penalty method, such that a large value of bulk modulus of water penalizes to impose the incompressible constraints. To alleviate mesh locking, we selected a selective reduced integration, or the B (B-bar) method [17]. The initial conditions of soil sedimentation and construction sequences were established as follows [13]: (1) impose the e€ective unit weight of 10.2 kN=m3 , with bulk modulus of water set to zero so that e€ective stresses are generated by gravity forces; (2) preload and remove with a surface traction of 60 kPa, thus forming overconsolidated state of soils simulating sedimentation and erosion of approximately 3m overburden layers; (3) exacavate by removing some clay elements and create a rigid foundation (E ˆ 1012 kPa, unit weight ˆ 25 kN/m3 ); (4) apply unit weight of pore water, 9.8 kN/m3 , setting bulk modulus of water to a large value, 108 ; (5) initialize the displacement vector. Vertical loading was then concentrated on the structure (see Fig. 6). The material parameters of clay deposits are shown in Table 2. Having established the initial conditions, the foundation was loaded by increments of 20, 50, 100, 200 and 500 kN load intensity. The nonlinear relationships of load±displacement are shown in Fig. 7. Errors are insigni®cant, and di€erences are not observable. The values of settlements are shown in Table 5, which shows nearly identical settlements for all load increments. As the applied load is divided into subincrements, the solutions stabilize. Therefore, numerical integrations were performed accurately.

Fig. 7. Load±settlement curves for various loading steps.

Table 5 Settlements and average numbers of iterations Load increment (kN)

20 50 100 200 500

Settlements for the corresponding load (cm)

No. iterations/step

1000 kN

1400 kN

1500 kN

14.749 14.745 14.738 14.725 14.737

36.883 36.879 36.860 36.768 ±

49.762 49.771 49.732 ± 48.419

4.5 5.0 5.0 5.3 6.3

266

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

Furthermore, the rate of convergence is rapid as shown in Table 5. As the rate of convergence is quadratic, the number of signi®cant digits in the error is almost doubled in each step. As shown in Fig. 8, force error norms show that the convergence rate is asymptotically quadratic, since the consistent tangent moduli conserve the eciency of Newton's method. As shown in Fig. 9, the GIH model predicted larger displacements than the Cam-clay model because of the plastic straining considered in the overconslidated states. In particular, the plastic strains could be modeled in the unloading±reloading sequences. Moreover, the reloading induced more plastic strains than the unloading. This can be controlled if the value of parameter c is assumed to be less than 1.0 [11]. The analysis was performed for a very soft clay, and thus resulted in large displacements. Ultimate capacities were estimated as about 1650 and 1700 kN by GIH model and the modi®ed Cam-clay model, respectively. In this problem serviceability was the more important point of consideration. The anisotropic hardening constitutive model assures the reasonable evaluation of serviceability.

Fig. 8. Convergence rate for 500 kN increment (d ˆ 0:5).

Fig. 9. Vertical load±displacement relationships.

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

267

5.3. Example 3: FE analysis of horizontally loaded foundation A rigid structure buried in a soft clay is loaded horizontally under the plane stain condition. Nonlinear analysis was performed for the same material parameters of Table 2 and for the similar mesh shown in Fig. 6 while the symmetric part on the left side was extended from the right side mesh. The initial stresses and construction sequences are similar to the previous example (see Section 5.2) while cyclic loading in the horizontal direction was concentrated at the center of structure. In this example, a drained loading condition was simulated without the process of applying the unit weight of pore water. The deformed mesh around the structure at the load of 250 kN is shown in Fig. 10, where the largest displacement occurred near the structure and the surface was heaving in the side of loading direction. The load±displacement relationships for the horizontal DOF are shown in Fig. 11. As described in Section 5.2, the initial state of soil layer was modeled as the overconsolidated stress state in which the initial stress is within the outer yield surface and hence the isotropic hardening Cam-clay model predicts just an almost linear elastic behavior. However the proposed GIH model can analyze the nonlinear response of the repeatedly loaded ground and also show the hysterisis loop. Each response tends to induce slightly degrading loops.

Fig. 10. Deformed con®guration at the horizontal load of 250 kN.

Fig. 11. Horizontal load±displacement relationship under drained and undrained condition.

268

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

Fig. 12. Horizontal load±displacement with various amplitudes relationship under undrained condition.

In Fig. 11 the drained relationship is more ¯exible than the undrained one where the volume was constrained. Since the drainage condition is dependent on the load amplitude and the loading period for two phase material with soil particle and pore water, it is hard to de®ne the drainage condition prior to the transient analysis or dynamic analysis. Therefore it is found that the mixed formulation is required in order to re¯ect the sti€ness in the soil-structure system, which is dependent on drainage condition, for a rigorous analysis. In Fig. 12, the response for various load amplitude is shown under the undrained condition. The larger load amplitude induced more ¯exible results and load±displacement loops were converged as load cycles continued. These behaviors are similar to those obtained from ®eld load tests and laboratory tests [18]. As a result, the GIH model can simulate the nonlinear behavior under cyclic loading condition for clays. 6. Conclusions In this study, an anisotropic hardening constitutive equation based on the generalized isotropic hardening rule was implicitly integrated by the generalized trapezoidal rule. The consistent tangent moduli were also formulated and coded into a nonlinear FE program. Though triaxial test examples it was veri®ed that the stress integration algorithm is accurate, and that the consistent tangent moduli conserve the quadratic rate of convergence. The accurate and ecient algorithm allow a reasonable analysis of foundation problems using the anisotropic hardening constitutive model; hence, any problems including complex initial conditions and various loading sequences could be analyzed. The GIH model suciently accounted for any plastic straining of clay in reverse loadings via FE solutions. From this study, we conclude that other rigorous anisotropic hardening models might also be implemented in a numerical method, based on an accurate implicit stress integration and ecient consistent tangent moduli. Appendix A. Jacobian matrix for stress integration The derivatives of yield surfaces with respect to stresses and internal variables are derived from (3.1b) as f;p ˆ

f;a ˆ 2…p

f;s ˆ

f;b ˆ 3…s

f;r ˆ

2r:

a†;

…A:1† 2

b†=m ;

…A:2† …A:3†

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

269

The center stress tensor of yield surface ac shows a geometric relationship [11] of ac …r; pc † ˆ ac;pc pc ‡ ac;r r;

…A:4†

where ac;pc ˆ h and ac;r ˆ ac =dpc h=d. T The Jacobians for residual vector r ˆ fr1 rT2 r3 r4 r5 g are derived from di€erentiating (4.6a)±(4.6e) as follows. A.1. Derivatives of r1 The derivatives of r1 are o…f;p †n‡d o…f;p †n‡d ; r1;s ˆ 0; r1;pc ˆ KD/ ; op opc o…f;p †n‡d ; r1;D/ ˆ K…f;p †n‡d ; ˆ KD/ or

r1;p ˆ 1 ‡ KD/ r1;rm

…A:5†

where from (2.6b), o…f;p †n‡d ˆ 2d; op

o…f;p †n‡d ˆ opc

o…f;p †n‡d or

2da;pc ;

2da;r :

…A:6†

A.2. Derivatives of r2 The derivatives of r2 are r2;p ˆ 0; r2;pc

r2;s ˆ I ‡ 2GD/

o…f;s †n‡d ˆ 2GD/ ; opc

o…f;p †n‡d ; os

r2;pc ˆ 2GD/

o…f;s †n‡d ; opc

…A:7†

r2;D/ ˆ 2G…f;s †n‡d ;

where o…f;s †n‡d 3d ˆ 2 I; m os

o…f;s †n‡d ˆ opc

3d b ; m2 ;pc

o…f;s †n‡d ˆ or

3d b : m2 ;r

…A:8†

A.3. Derivatives of r3 The derivatives of r3 are r3;p ˆ r3;r ˆ

o…f;p †n‡d ; op o…f;p †n‡d ; Kpc D/ or Kpc D/

r3;s ˆ 0;

r3;pc ˆ 1

Kpc D/

o…f;p †n‡d ; opc

…A:9†

r3;D/ ˆ Kpc o…f;p †n‡d :

A.4. Derivatives of r4 The derivatives of r4 are r4;p ˆ 0;

r4;s

r 2 okep k ; ˆ 3 os

r4;pc

r 2 okep k ˆ 3 opc

w;pc ;

r4;r

r 2 okep k ˆ 3 or

w;r ;

r4;D/

r 2 okep k ; ˆ 3 oD/ …A:10†

270

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

where epij ˆ …epij †n ‡ D/…of =osij †n‡d okep k ep oep ˆ p : ; os ke k os

o…f;s †n‡d oep ˆ D/ ; os os

…A:11†

okep k ep oep ˆ p : ; opc ke k opc

o…f;s †n‡d oep ˆ D/ ; opc opc

…A:12†

okep k ep oep ˆ p : ; or ke k or

o…f;s †n‡d oep ˆ D/ : or or

…A:13†

From Eq. (3.5), ow a ˆ 2 opc m hpc

b…cr=dpc †

b

…1

…1

b†…cr=dpc †

cr=dpc †

b‡1

b 1

2

b

ow ac b…cr=dpc † ‡ …1 b†…cr=dpc † ˆ 2 : 2 or m hdpc …1 cr=dpc † …A:14†

;

A.5. Derivatives of r5 The derivatives of r5 are r5;p ˆ f;p ;

r5;s ˆ f;s ;

r5;pc ˆ f;a a;pc ‡ f;b : b;pc ;

r5;r ˆ f;a a;r ‡ f;b : b;r ‡ f;r ;

r5;D/ ˆ 0:

…A:15†

Appendix B. Consistent tangent moduli From (4.10a), we obtain a1 p0 ‡ a3 pc0 ‡ a4 r0 ‡ a5 D/0 ˆ K1;

…B:1†

where a0 ˆ 2Kd…D/†;

a1 ˆ 1 ‡ a 0 ;

a3 ˆ

a0 a;pc ;

a4 ˆ

a0 a;r ;

a5 ˆ K…f;p †n‡d :

…B:2†

Eq. (4.10b) is also written as b2 s0 ‡ b3 pc0 ‡ b4 r0 ‡ b5 D/0 ˆ 2G I

1 1 3




1 ;

…B:3†

where b0 ˆ 6Gd…D/†=m2 ;

b2 ˆ 1 ‡ b0 ;

b3 ˆ

b0 b;pc ;

b4 ˆ

b0 b;r ;

b5 ˆ 2G…f;s †n‡d :

…B:4†

From (4.12) c1 p0 ‡ c3 pc0 ‡ c4 r0 ‡ c5 D/0 ˆ 0;

…B:5†

where c0 ˆ Kpc ;

c1 ˆ

2c0 d…D/†;

c3 ˆ 1

c1 a;pc ;

c4 ˆ

c1 a;r ;

c5 ˆ

c0 …f;p †n‡d ;

…B:6†

and from (4.13) d2 : s0 ‡ d3 pc0 ‡ d4 r0 ‡ d5 D/0 ˆ 0;

…B:7†

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

271

where p 6dD/ ; d2 ˆ d0 ep ; m2 kep k r 2 ep : …f;s †n‡d : d5 ˆ 3 kep k

d0 ˆ

d3 ˆ

w;pc

d0 ep : b;pc ;

d4 ˆ

w;r

d0 ep : b;r ; …B:8†

Finally (4.14) becomes e1 p0 ‡ e2 : s0 ‡ e3 pc0 ‡ e4 r0 ˆ 0;

…B:9†

where e1 ˆ f;p ;

e2 ˆ f;s ;

e3 ˆ f;a a;pc ‡ f;b : b;pc ;

e4 ˆ f;a a;r ‡ f;b : b;r ‡ f;r :

…B:10†

Eqs. (B.1), (B.3), (B.5), (B.7) and (B.9) are composed of ®ve independent variables p0 , s0 , pc0 , r0 and D/. We can thus solve a linear system such that 2

a1

0T

6 6 0 b2 I 6 6 6 c1 0T 6 6 0 dT 4 2 e1 eT2

a3

a4

b3

b4

c3

c4

d3

d4

e3

e4

9 8 9 38 K1T > > p0T > a5 > > > > > > > > > > > > 7> > > > > 1 0T > > > > 2G…I 1

1† s b5 7 > > > > 3 = = < < 7 7 0T T ˆ : pc c5 7 0 > > 7> > > > 0T > > > > > > 7 > > r > d5 5> 0T > > > > > > > > > > > > ; : : 0T ; T D/ 0 0

…B:11†

Therefore, the consistent tangent modulus …Ckn‡1 †ijkl in (4.9) can be obtained by solving op=oekl and osij =oekl . The modulus can also be obtained by an explicit form. However, this is a dicult task with rigorous anisotropic constitutive models, and it also requires more expensive computations for inversions than a numerical solution of the linear system.

References [1] M. Oritz, E.P. Popov, Accuracy and stability of integration algorithms for elastoplastic constitutive equations, Int. J. Numer. Meth. Engrg. 21 (1985) 1561±1576. [2] J.C. Simo, R.L. Taylor, Consistent tangent operators for rate-independent elastoplasticity, Comput. Meth. Appl. Mech. Engrg. 48 (1985) 101±118. [3] R.H. Dodds, Numerical techniques for plasticity computations in ®nite element analysis, Comput. Struct. 26 (1987) 767±779. [4] A. Gens, D.M. Potts, Critical state models in computational geomechanics, Engrg. Comput. 5 (1988) 178±197. [5] R.I. Borja, S.R. Lee, R.B. Seed, Numerical simulation of excavation in elasto-plastic soils, Int. J. Numer. Anal. Meth. Geomech. 13 (1989) 231±249. [6] M.A. Cris®eld, Non-linear Finite Element Analysis of Solids and Structures, vol. 1: Essentials, Wiley, New York, 1991. [7] Y.F. Dafalias, Bounding surface plasticity. I: mathematical formulation and hypoplasticity, J. Engrg. Mech. Div., ASCE 112 (1986) 966±987. [8] Z. Mroz, V.A. Noris, O.C. Zienkiwicz, An anisotropic critical state model for soils subject to cyclic loading, Geotechnique 31 (1981) 451±469. [9] S. Somasundaram, C.S. Desai, Modeling and testing for anisotropic behavior of soils, J. Engrg. Mech. Div., ASCE 114 (1988) 1473±1496. [10] S. Oh, An anisotropic hardening constitutive model of clays, Ph.D. Thesis, Department of Civil Engineering, KAIST, Taejon, Korea, 1996 (in Korean). [11] S.R. Lee, S. Oh, An anisotropic hardening constitutive model based on generalized iostropic hardening rule for modelling clay behavior, Int. J. Numer. Anal. Meth. Geomech. 19 (1995) 683±703. [12] S. Oh, S.R. Lee, Nonlinear ®nite element analysis of excavation by anisotropic hardening constitutive elasto-plastic constitutive model, in: Proceedings of the Third Asian-Paci®c Conference on Computational Mechanics, Seoul, 1996, pp. 2023±2028.

272

S. Oh, S.R. Lee / Comput. Methods Appl. Mech. Engrg. 191 (2001) 255±272

[13] R.I. Borja, S.R. Lee, Cam-clay plasticity, part I; Implicit integration of elasto-plastic constitutive relations, Comput. Meth. Appl. Mech. Engrg. 78 (1990) 49±72. [14] S.R. Lee, Non-linear elasto-plastic ®nite element analysis of braced excavations in clays, Ph.D. Thesis, Department of Civil Engineering, Stanford University, Stanford, CA, 1989. [15] K.H. Roscoe, J.B. Burland, On the generalized stress±strain behavior of `wet' clay, in: J. Heyman, Leckie (Eds.), Engineering Plasticity, Cambridge University Press, Cambridge, 1968, pp. 535±609. [16] D.J. Henkel, The e€ects of overconsolidation on the behavior of clays during shear, Geotechnique 6 (1956) 139±1509. [17] T.J.R. Hughes, The Finite Element Method, Prentice-Hall, Englewood Cli€s, NJ, 1980. [18] B.M. Das, Principles of Soil Dynamics, PWS-KENT, 1993.