Abaqus
Abaqus
Abaqus
Colby C. Swan 4120 Seamans Center for Engineering Arts & Sciences Department of Civil and Environmental Engineering Center for Computer-Aided Design The University of Iowa Iowa City, Iowa 52242, USA E-mail: colby-swan@uiowa.edu Phone : 1(319) 335-5831 Young-Kyo Seo, Post-Doctoral Research Associate Department of Civil and Environmental Engineering Computational Solid Mechanics Laboratory The University of Iowa Iowa City, Iowa 52242, USA E-mail: yseo@diego.ecn.uiowa.edu Phone : 1(319) 335-6455
Abstract
One of the primary strengths of elasto-plastic cap models is their ability to capture the gross inelastic coupling between deviatoric and volumetric behaviors of many porous media. While numerical integration algorithms for these models have been presented in the literature, performing implicit analysis of earthen systems using cap models remains a challenging endeavor. One of the difculties associated with most isotropic cap models is that the three independent surfaces comprising the yield surface intersect non-smoothly. It is shown here that the elastoplastic tangent operators at the corner points on such yield surfaces are singular, giving rise to potential numerical difculties. To address this issue, a novel, three-surface elasto-plastic cap model in which the three surfaces intersect smoothly is introduced and developed here. The rate form of the model constitutive equations are rst presented, followed by an unconditionally stable integration algorithm with expressions for consistent tangent operators. Sample computations demonstrating the very good performance of the model in slope stability and bearing capacity problems are also presented. Key words: cap models; elasto-plasticity; soil models; computational plasticity;
Table of Contents
Page
1. 2. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . ASSESSMENT OF A NON-SMOOTH CAP MODELS . . . . . . . . . .
. . . 1 . . . 4
2.1 2.2
3.
Basic Forms and Rate Equations . . . . . . . . . . . . . . . . . . 4 The Problem: Singular Tangent Operators at Corner Points . . . . . . . 7
. . . 9
Stress Updates and Active Yield Surface Determination . . . . . . . Case 1 Integration Algorithm . . . . . . . . . . . . . . . . . . Case 2 Integration Algorithm . . . . . . . . . . . . . . . . . . Case 3 Integration Algorithm . . . . . . . . . . . . . . . . . . Consistent Tangent Operators . . . . . . . . . . . . . . . . . .
. .
EXAMPLE COMPUTATIONS . . . . . . . . . . . . . . . . . . . .
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . .
ii
Page
List of Figures
1. 2. 3. 4. 5. 6. 7. 8. Non-smooth, three-surface, two-invariant cap model yield surface with two corner points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Smooth, three-surface two-invariant yield function for cap model . . . . . . .
10 17 17 28 29 31 32
tr 1
0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
tr 1
Single element hydrostatic compression test . . . . . . . . . . . . . . . . Four-element limit state analysis computation . . . . . . . . . . . . . . .
Mesh and results of rigid foundation bearing capacity computations on loose and dense sandy soils . . . . . . . . . . . . . . . . . . . . . . . . . Mesh and results of slope stability computations for loose and dense sandy soils . .
iii
Page
1. 2. 3. 4. 5. 6. 7. 8. 9. 10.
Box 2: Algorithm for determination of yield surface activity . . . . . . . . . Box 3: Iterative closest point return map to the Drucker-Prager surface . . . . . Box 4: Fully implicit return map for CASE 2
. . . . . . . . . . . . . . .
and
Box 6: Return map algorithm for CASE 3 . . . . . . . . . . . . . . . . . Table 1: Material parameters used in hydrostatic compression test . . . . . . . Table 2: Material parameters used in 4-element limit analysis test . . . . . . . Table 3: Material parameters used in bearing capacity computations . . . . . . Table 4: Material parameters used in slope stability computations . . . . . . .
iv
13 16 19 22 23 24 28 29 30 32
. . . . . . . . . . . .
1.
stress-strain models for soils that provide both adequate physical representation of observed mechanical behaviors, and also sound numerical performance in implicit computational analyses. While physical realism in soil models is clearly necessary, it alone does not permit them to be successfully employed in limit state computational analysis of earthen structures. For successful usage in implicit analysis of geomechanical structures, soil models must have numerical implementations that are also be complete, stable, and continuously differentiable as will be discussed below. This will allow for efcient simulation of the inelastic load-deformation response of earthen structures with relatively rapid convergence, and thus enhanced computational efciency. Elastoplastic soil models are a class of material models that can provide some degree of both physical realism and potentially sound numerical performance. One of the benets of such models is that their underlying theory has become increasingly well established, having beneted from preceding developments in metal plasticity theory. One of the pioneering extensions of metal plasticity theory to soil plasticity was performed by Drucker and Prager [1,2] when they extended the von Mises yield criterion to account for connement strengthening of granular media. The result was their new yield criterion, which is now commonly called the Drucker-Prager criterion. In 1957, Drucker et al. [3] further proposed that the volumetric plasticity behavior of soils might be successfully modeled with a strain a hardening compression cap surface that closes off the open end of the Drucker-Prager failure envelope. While preceding soil models and yield criteria such as Mohr-Coulomb had captured frictional connement strengthening behavior, this strain hardening model was one of the rst that attempted to couple the deviatoric and volumetric deformation behaviors of granular media. Since then, a variety of alternative strain hardening plasticity models have been developed [4-7], but in this work, attention will be conned to Drucker-Prager models with hardening compression cap surfaces. In 1971, DiMaggio and Sandler [8] proposed a specic elasto-plastic cap model (Figure 1) and numerical implementation numerical implementation. Since the model was apparently intended to be used in explicit nite difference and/or nite element computations, the implementation assumed extremely small strain increments, and the treatment did not deal with material tangent
stiffness operators. While cap models such as this were developed initially for sands, they have also been successfully employed with other materials such as clays, and concrete [9-11]. These cap models are based on classical isotropic elasto-plasticity theory, and couple the Drucker-Prager failure envelope and with a hardening compression cap surface. Since they are generally used with associated ow rules, they can permit dilatancy (that is, incremental plastic volume increase under shear loading) when loading occurs on the Drucker-Prager failure envelope and the tension cap surface, and plastic compaction (incremental plastic volume decrease) when loading occurs on the compression cap surface. In the family of cap models advanced by DiMaggio and Sandler, the hardening cap is an elliptical surface with a constant ratio of major to minor radius, and it intersects
the failure envelope in a non-smooth fashion. The compression cap surface moves along the
axis with incremental changes in plastic volumetric strain. Under dilatation, the compression cap moves inward while contracting, and under plastic compaction, it moves outward while expanding.
()
Figure 1: Non-smooth, three-surface, two-invariant cap model yield surface with two corner regions.
The intersection point, or the so-called corner point, at which the Drucker-Prager envelope and the moving compression cap intersect has long been recognized as point that can give rise to both numerical and constitutive instabilities [12]. One of the difculties stems from the fact that the original hardening law [8] for such cap models is posed in a form that gives rise to one-to-one
||s||
Tension Cutoff
I1
in a form
similar to that expressed below in Equation (8). However, if such a one-to-one hardening law is actually used, then unexpected and perhaps undesired softening behavior may occur when the stress point is at the compression corner point [13,14]. Motivated by a desire to control and eliminate unwanted softening behavior at the compression corner point, Sandler and Rubin [13] proposed a modied discontinuous hardening law which prevented softening material responses. Resende and Martin [14] subsequently showed, however, that a discontinuous hardening law can lead to an incomplete model in the sense that there is a region in strain rate space at the compressive singular corner point which cannot be covered by any possible stress rate. The basic reason for this lack completeness was traced to the discontinuity of the hardening law at the compression corner. Using an elegant integration theory for non-smooth multi-surface plasticity [15], based on closest point projection return mapping algorithms [16] Simo et al. [17] then developed a robust numerical implementation of the cap model in terms of algorithms for determining yield activity based on the Kuhn-Tucker conditions. As part of their new implementation, Simo et al. [17] proposed a modied continuous hardening law which precluded softening and avoided the completeness issue raised by Resende and Martin [14]. They formulated consistent algorithmic loading/unloading conditions and closest point projection return mapping algorithms for all possible modes of yield surface activity. In a subsequent work [18], Hofstetter, Simo and Taylor (1993) modied the cap models hardening law such that it was associated, thus resulting in symmetrical consistent tangent operators. Nevertheless, even with the vital improvements advanced in [15-18], problems remain with the non-smooth, three-surface cap model. As will be shown below, one of the primary remaining difculties is that the material tangent operators in both the compression corner and tension corner regions are singular in that they do not provide any bulk/volumetric material stiffness. This is completely unrealistic, and can lead to difculties in structural analysis of soil systems. While there are a number of ad-hoc ways to deal with this singularity of tangent operators at the corner points (such as the introduction of articial bulk stiffness in the corner regions, or resorting to usage of visco-plasticity as proposed in [15] for general non-smoothess), these approaches are less than satisfactory. The intent of this manuscript is, therefore, to introduce and fully develop a novel, smooth cap model wherein the three surfaces (a Drucker-Prager failure envelope; a hardening circular
compression cap; and a xed circular tension cap) intersect in a smooth fashion such that there are no corner regions. This smoothness therefore precludes many of the historical difculties associated with the corner regions including: material softening response, lack of completeness, and singularity of material tangent operators. The model features essentially the same physical degree of realism as the preceding cap models, but its numerical performance characteristics are signicantly enhanced. In Section 2 of this manuscript, a non-smooth cap model is briey introduced, and the difculty with singular tangent operators at the corner regions is highlighted. Section 3 then presents a novel smooth cap model along with a detailed integration algorithm and expressions for consistent tangent operators. The integration algorithm presented is based on a Backward Euler integration of the rate constitutive equations, which give rise to an elasticpredictor, plastic-corrector (return map) stress update algorithm. Differentiation of the incremental stress update algorithms provides expressions for the so-called consistent tangent operators which facilitate good convergence characteristics in implicit structural analysis of soil structures. In Section 4, the excellent performance of the model is rst demonstrated on a number of simple material test type computations, and then on full scale slope stability and bearing capacity type problems.
2.
2.1
For simplicity and clarity, the elastoplastic constitutive equations of a non-smooth three surface cap model are considered here in a small deformation framework. These equations can be straightforwardly extended to a large-deformation framework as necessary. The small strain tensor admits the usual additive elastic, plastic decomposition as follows: = +
where
and
are the total, elastic, and plastic strain vectors, respectively. The elastic response
of the material is assumed to be characterized by a constant isotropic tensor C = such that the incremental stress response of the material is given by: = C:
(1)
1 + 2 I dev
(2)
In stress space, the elastic domain is bounded by three distinct yield surfaces which are functions
s = Idev : ). The three surfaces comprising the yield surface intersect in a non-smooth manner as shown in Figure 1.
= 1 2 3 are: (3)
The hardening law for this model derives from the fact that the volumetric crush curve (plastic
volumetric strain
plastic volumetric strain in the soil (or porous medium) as measured from a virgin completely unloaded state;
2
represents the maximum possible plastic volumetric strain for the medium, with
1
8 @
(absolute) of
at which
8 94
the reference state being the materials virgin unloaded state; and
where ( )
axis;
ref 1
5 64
2 3
the left (
versus
1)
and cohesion . The aspect ratio of the elliptical compression . The cap is permitted to translate along the
1
, the constants
and
yield surfaces
= 0 and
and
' )!
' (
0,
0,
0,
0, and
'
'
&
)=
( )
$ %
! "
( 1) =
and
are: exp
1
3(
)=
2(
)=
1(
)=
( 1)
0 ) 0
0. The
(i.e.
(4) (5)
(6) (7)
dimensionless constant providing the ratio of major to minor radii of the elliptical compression cap
( ) ( )
The ow rule for this non-smooth model is associated, and since multiple surfaces are potentially active at any given instant, it takes Koiters generalized form:
proportional to the instantaneous magnitudes of the plastic deformation processes at a xed material point with respect to each of the three yield functions. Loading and unloading criteria are specied by the Karesh-Kuhn-Tucker conditions as
'
0;
0;
In accordance with the KKT conditions, this elasto-plastic constitutive model poses six different different possibilities: 0. The stress point lies inside of yield surface and 0, = 1 2 3. In this case the material response is incrementally elastic; 1. Loading is occuring on surface 1, such that 1 = 0 and 1 0; 2. Loading is occuring on surface 2, such that 2 = 0 and 2 0; 3. Loading is occuring on surface 3, such that 3 = 0 and 3 0; 4. Surfaces 1 and 2 are simultaneously active 1 = 2 = 0 and 1 0, 2 0; 5. Surfaces 1 and 3 are simultaneously active 1 = 3 = 0 and 1 0, 3 0; In essence, this model is comprised of ve elasto-plastic sub-cases as listed above. Those that can create some difculty are the corner cases 4 and 5 as briey detailed below.
=0
and
(10)
(11)
(12)
(13)
2 3
where
=1
then
( )=
modulus
( ) for
as follows: ] (9) ( )
To illustrate the difculty associated with the corner cases of non-smooth cap models, the return mapping algorithms for the corner case 4 is briey formulated below. Upon differentiating the case 4 stress update algorithm, it is shown that the associated consistent tangent operator obtained is singular in that it provides no volumetric stiffness for the material. This can lead to singular (rank decient) global tangent stiffness matrix operators, and severe computational difculties when attempting to solve global force balance nite element equations. While band-aid type remedies to this problem are available, the difculty can be avoided altogether by resorting to a cap model in which the yield surfaces intersect in a smooth fashion. Such a model is introduced in the following section. If a Backward Euler integration algorithm is applied to the elastoplastic rate constitutive equations laid out above, then a stress update algorithm with an elastic predictor and a plastic corrector (return mapping) is the result. The elastic predictor assumes incrementally elastic behavior, and leads to a trial stress state as follows:
If the elastic predictor lies outside of the elastic domain, then one or more of the intersecting yield surfaces will be active. If it is assumed that the elastic predictor lies in the attractor region for the corner point at which the compression cap and failure envelope intersect (i.e. the compression corner region of Figure 1) then the return map (plastic corrector) brings the stress state back to the indicated corner point. Specically,
+1 +
where (
+1
+1 )
+1
= Idev : (
the deviatoric plastic strain increment. The total plastic strain increment associated with the return
+1
+1
+1
= 1 +1
+ 2 +1
map is
+1 1
tr
+1
tr
tr
+ C:
+1
(14) (15)
+1
(16)
+1 )
is
(17)
The magnitudes of the plastic deformation proceesses 1 +1 and 2 +1 are computed by enforcing ( 1)
+1
= 0 and ( 2 )
+1
= 0, which provides
Hence the return map to the corner point is very simple and involves no iterations. However, when one differentiates the stress increment produced by this integration algorithm with respect to the driving strain increment, a singular consistent tangent operator is obtained. Differentiating Eq. (16) [and the subsidiary terms in Eqs. (17)-(19)] with respect to
the following consistent tangent operator when surfaces 1 and 2 are simultaneously active:
tr +1 +1
+1
where n
+1
is the deviatoric component of the unit normal to the yield surface at the converged
stress point. Clearly, this tangent material stiffness operator has no bulk stiffness, and thus it is singular. A similarly singular tangent operator with no bulk stiffness occurs at the corner point between the tension cutoff surface and the failure envelope. It is the authors experience that such singular material tangent operators can give rise to rank decient global tangent stiffness matrices when implicit nonlinear nite element calculations are performed with non-smooth cap models. There are a number of band-aid type solutions that can be employed to deal with this problem. Two in particular are: i introduction of an articial (and small) bulk stiffness to the material tangent dened in Equation (20); and ii introduction of visco-plasticity to the model, which allows the tangent operators in corner regions to remain non-singular as long as the inviscid limit is not approached. Neither of these band-aid type approaches deal with the root cause of the problem, however, which is the existence of corner points on the yield surface. Furthermore, the introduction of the articial bulk stiffness in the corner regions can lead to a tangent operator that is inconsistent with the stress update algorithm, and thus result in slow convergence behavior in solving nonlinear global force balance equations. The introduction of viscoplastic behaviors to avoid singularity of tangent operators in
&
=2
[Idev
+1 ]
+1
2 ( 1 +1 + 2 +1 )
+1
2 +1 =
+1
tr ( 1)
1 +1 =
tr ( 1)
+1 +1 +1
(18) (19)
+1
provides
(20)
the corner regions makes the structural analysis problem somewhat more involved, since for each structural load increment applied, a number of time steps must be permitted to compute the timedependent deformation response of the structure. Furthermore, the rate dependence introduced with viscoplasticity can obscure more physically based rate effects such as those associated with pore pressure diffusion phenomena. A more fundamental approach to this problem of singular tangent operators in corner regions is therefore proposed below, with the introduction of a smooth cap model having no corner regions.
3.
3.1
To deal with the difculties associated with a non-smooth yield surface, a yield surface (Figure 2) made up of three smoothly intersecting yield functions is introduced here. The form of the
are respectively called the Drucker-Prager envelope function, the compression cap function, and the tension cap function:
1(
where
=s
and
concert with the non-smooth cap model and are dened here as
1
( 1) =
)=
( )
)2
( )
( 1) =
and
exp(
1)
( )
3(
)=
( 1)
2(
)=
)=
( 1)
0 ) 0
)(
and
which
10
||||
()
IC1()
I1
IT1
envelope, and a moving delimiting point between the Drucker-Prager envelope and the compression cap.
primary differences between the yield functions for the smooth cap model and those for preceding generations of non-smooth cap models are: a. The failure envelope function features a saturation of strength gains associated with increased conning stresses (or 1 ). This feature is very important when the model is to be used in limit state analysis of geo-structures, as in [20,21]. b. The compression cap surface is no longer elliptical, but rather strictly circular in the - 1 space representation. As in the non-smooth model, the cap surface can translate left (outward) on the 1 axis with plastic compression of the medium, and right (inward) on the 1 axis with plastic dilatation. d. The radius ( ) of the compression cap is such that the it maintains smooth tangency with the exponential failure envelope. Hence the radius ( ) is the minimal distance from the current centerpoint ( 0) of the compression cap, to the failure envelope surface 2 = 0. Due to the nonlinearity of the failure envelope surface, a closed form expression for ( ) is not available, although it can be determined algorithmically as described below in Box 1. e. The tension surface is no longer a planar cutoff, but rather circular cap. The center of the tension cap resides at 1 = 0, and the radius of the surface is a constant which is determined (Box 1) based on , , and .
delimiting points [
and
( )] as well as
( ) denote, respectively, a xed delimiting point between the tension cap and the Drucker-Prager
denotes the xed radius of the tension cap. Specic algorithms for determination of these are presented in the following subsection. The
and
and
axis at
= .
11
Since the Mohr-Coulomb yield criterion has a long history of usage in classical soil mechanics, geotechnical engineers often think of soil shear strength characteristics in terms of the Mohr0 1
1(
Coulomb cohesion
with tension and compression caps, and with saturating frictional effects is employed since: 1. It captures the saturation of soil strength with increasing effective conning stresses; and 2. It is a completely smooth model having no corner regions. The Drucker-Prager model does not suffer from the non-smooth corner regions that generally afict MohrCoulomb type soil models [22]. So that the results of the present work can be placed in context and even compared with results from classical geotechnical analysis methods wherein the Mohr-Coulomb failure criterion is routinely considered, the failure envelope in Equations (21) and (24) is re-written by taking its Taylor series
1(
)=
which is simply a variation of the linear Drucker-Prager yield criterion. A correspondence can
and
example, translations from Mohr-Coulomb parameters to linear Drucker-Prager parameters have been provided by Chen and Saleeb (1982) as:
1 2
These equations can be inverted to provide a translation from linear Drucker-Prager envelope parameters to Mohr-Coulomb which will be used in the following section.
8 0
tan
4 1 + tan2 3 2
1 2
4 1 + 3 tan2
4 3 1 + 3 tan2
=
1 2
2 tan
1 2
where
)=
1+
1
24
&
(i.e.
1
expansion about
= 0:
1
1)
1)
(28)
(29)
(30)
12
While this generalized form of the ow rule would appear to imply that more than one yield function can be active at a given instant, this will not occur with this model due to the smoothness of the yield surface of the cap model. The Karesh-Kuhn-Tucker loading/unloading conditions and the plastic consistency conditions for the smooth model are expressed in a manner identical to that for the non-smooth cap model in Equations (12) and (13). In rate form, the non-associated hardening
( )=
kinematic hardening law is employed with this model, the rate form of which is
3.2
With the non-smooth cap model, the minor radius of the elliptical compression cap surface can
the cap surface smoothly intersects with the exponential Drucker-Prager envelope. For a linear Drucker-Prager failure envelope, expressions are available for the cap surface radius as a function of
circular compression cap is determined such that when the cap is centered at (
Idev
as ( ) =
= min 0
. The parameter ( ) still denes the intersection of the cap with hydrostatic stress
where
and
are material parameters which retain the same meaning as they have in the noncan take is limited to 0, such that
5 64 2 5 64 @
where
( ) is the tangent hardening modulus for the cap parameter which is still given by [ ( )] ( ) (33)
( ) ( )
(31)
(32)
(34)
13
[Seo, 1998]. However, for the nonlinear exponential failure envelope described in Equation (24)
are required, an algorithm for determination of both are provided below in Box 1. The essential
End while
2
( ) 2
[2
serves as the delimiting point between the tension cap and the Drucker-Prager failure envelope.
This point is easily found by applying the algorithm of Box 1 with the value
= 0.
3.3
The basic problem of integrating the elastoplastic constitutive equations at a xed material point can be stated as follows. On the time interval of interest [0
], it is assumed that at
For the tension cap surface of xed radius , it is necessary to nd the abscissa value
8
!
8
( )=
+(
( )= (
)+(
)2
1 2
+1
=2
) + 2(
+1
=(
= 2(
While(abs( )
)]
=2
( 1 ) + 2(
Initialize: k = 0;
( )=
) .
1 2
= ; ) (B1.1)
= 0, then
)=
( 1) + (
1
= 0. Hence the
)2 is the
(B1.5) (B1.6)
which
algorithm and expressions for consistent tangent operators that follow, values for
( ), and
14
+1 ],
is accomplished here using well-established operator-splitting, elastic-predictor, plastic-corrector methods. Briey, the updated stress can be written as follows:
+1
+C:
+1
+C:
+1
C:
+1
+1
incremental material response which is fully elastic, leading to a the so-called incrementally elastic trial stress predictor computed as
= =
where
+1
+1
Taking the trace and deviatoric parts of the elastic trial stress leads to:
( 1 )tr +1 , ( 2 )tr +1 , ( 3 )tr +1 are then computed based on the trial stress and hardening parameters. If
( 1 )tr +1
TOL1 or ( 2 )tr +1
TOL2 or ( 3 )tr +1
outside of the elastic domain, and a plastic correction must be performed. The manner in which the plastic correction is computed depends upon which one of the yield constraints is in fact active. For the general case of non-smooth plasticity models, determination of the so-called active set of yield contraints can be a challenging task as demonstrated by Simo et al [55]. Even for this smooth plasticity model, great caution is needed in determining which, if any, of the three yield constraints will be active at any given instant. The algorithm used for
+1
and
+1
tr
remain unchanged
+1
= s + 2 Idev
+1
( 1)
+1
= ( 1) + 9
+1
(37 ) (37 )
+1 1
+2
tr
+ C:
+1
+1
+1
+1
+1
+1
variables: that is
] the total and plastic strain tensors are known as are the stress and hardening
+1
over
(35)
(36 ) (36 )
15
determining yield function activity with the proposed smooth cap model is shown in Box 2. The underlying concept for the algorithm of Box 2 is that since the Drucker-Prager surface envelopes both the tension and compression cap surfaces, when it is violated , then the second and third yield
surfaces could in fact be active. Ultimately, however, only one will indeed be active. Due to the
smoothness of the yield surface and the nonlinearity of the hardening law for
a priori criterion for determining which of the three yield constraints will be active based on the elastic trial stress state. Since a determination of active yield constraints is not possible from the trial stress, a determination of activity must be determined from updated stress states. Consequently, the closest point projection return map of the elastic trial stress to the Drucker-Prager surface is performed and an
update of consistent with the return map is performed (leading to a value denoted by
case1 +1 ).
returned stress point lies in the updated domain of the Drucker-Prager surface, then that constraint will indeed be active. However, if the return point should lie in the updated domain of either the compression cap or the tension cap, they would represent the active constraint. These ideas
then the return point lies in the domain of the tension cap surface, and so the tension cap surface will be active which means that the return map must be re-computed using the Case 3 integration
the compression cap is active and the return map must be re-computed using the Case 2 integration algorithm. It should be further noted, corresponding to Figure 4, that it is possible for the elastic trial
+ 2 3 1
it must be determined whether Case 2 or Case 3 is active. A criterion based on the location of
algorithm.
+1
with respect to
and
stress state
0 and
tr
0. In this case,
case1 +1 )
case1 +1 ),
there is no simple
If the
then
16
case1 +1
goto Case 2.
Elseif ( 1 )case1 +1
goto Case 2.
Elseif (
tr 3
TOL3 and
tr 1
goto Case 3. Else Elastic predictor lies in elastic domain. Endif Endif Box 2. Determination of yield surface activity.
3.4
surface
If (
tr 2
TOL2 and
tr 1
tr
If ( 1 )case1 +1
case1 +1 )
If (
tr 1
TOL1 )
+1 ))
based on
tr 1,
tr 2,
tr 3
tr
+1
and
tr
+1
+1
tr
tr
+1 .
+1 .
is returned to the
17
|||| (,I ) 1 (,I ) 1
tr n+1
||||
n+1
case1
n+1
()
3a
I1C
I1
I1T
I1
T 1
3b
(,I ) 1 (,I ) 1
tr n+1
case1
||||
n+1
()
case1
I1C
I1
I1T
3c
n+1
Figure 3: 3a) Trial stress point that violates 1 0 and thus violates 2 0 and 3 0 also; 3b) Case 1 return point to the Drucker-Prager surface but lying in the tension cap domain; and 3c) Case 1 return point to the Drucker-Prager surface lying in the compression cap domain.
|||| (,I1)
tr
||||
n+1
(,I1) () n
I1C
tr
n+1
4a
I1
I1T
I1
IT
1
4b
||||
(,I1) () n
I1C
tr
n+1 I1T
I1
4c
Figure 4: 4a) Trial stress point that violates 2 0 and lies in the domain of the compression cap surface; 4b) Trial stress point that violates 3 0 and lies in the domain of the tension cap; and 4c) Trial stress point that violates both 2 0 and 3 0 and still lies in the elastic domain.
18
+1
1.
Using a Backward Euler integration rule which gives rise to plastic-correction of the elastic-
+1
+1
+1
1 +1 n
The deviatoric and trace portions of the updated stress can be written as
+1
+1
) 1 +1
(42)
To complete the stress update, the plasticity consistency parameter 1 +1 is computed such that the stresses lie on the Drucker-Prager surface
1 +1 is provided in Box 3. Once 1 +1 is computed and the stresses have been updated, then the
+1
+1 )
+1 )
3 1 +1 (
+1 )
+1
(43)
+1 .
( 1)
+1
= ( 1)
+1
+9
( 1)
1 +1
(41 )
+1
=s
2 1 +1 n
C:
+1
+3
+1
+1
tr
tr
2 1 +1 n
1 +1
where n =
+1 +1
+1 +1
is the normalized deviatoric component of the normal vector to the with respect to
+1
& $
= 1 +1
= 1 +1 n
+1
(38)
(39)
(40)
(41 )
19
+1
+1
+1
3.5
When the algorithm of Box 2 indicates that the compression cap surface is active (Case 2), then the elastic stress predictor must be returned to the compression cap surface, which will generally translate and grow/shrink during the return map process. (Note: If the portion of the compression
( ).
Backward Euler integration of the Case 2 associated ow rule gives a plastic strain increment
+1
+1
+1
& $
= 2 +1
= 2 +1 2
as follows
and
( 1)
+1
+1
+1
((
+1 ) +1 )
update (
+1 ) +1
such that (
+1 ) +1
+1 )
9 ( 1 +1 )
+1
((
+1 ) +1 )
(2 +
)( 1 +1 )
+1
( 1 +1 )
= ( 1 +1 )
( 1) ( 1)
( 1) =
(2 +
+1 1
9 1 9 1 +1
While (abs( ( 1 ) )
( 1) =
+1
((
+1 )
Initialize: k = 0; ( 1 +1 ) = 0; (
+1 )
=(
+1 )
+1
+1
+1
(B3.1)
(44)
20
+1
Decomposing the updated stress into its deviatoric and bulk components provides:
+1 +1
+1
+1
+1
Taking the difference between Eqs. (46a) and (47b) yields the following useful results: 1 + 2(2 + )
+1
The objective of the plastic correction stress update algorithm for Case 2 is thus to solve for
the value of the generalized plastic consistency parameter 2 +1 , as well as Equations (46b), (47a) and the plastic consistency condition, namely,
Hence, the Case 2 return map involves solving a system of three highly nonlinear equations [(46b),
+1
and
system of equations can be reliably obtained, a robust and fully implicit algorithm is developed below. Since a sequential linearization algorithm is used to solve the system, the total derivative of with respect to
+1
+1 2 +1
2 4(2 + ) +1 1 + 2(2 + ) 2 +1
2 +1
2 2 +1
+1 2 +1
1 + 2(2 + ) 2 +1
2 +1
1 +1
+1 )
=0
+1 .
+1
1 + 2(2 +
) 2 +1
+1
+1
+1
+1
(48)
(49)
(50)
(51)
+1
+2
+1
2 +1
(47 )
+1
3 (
2 +1 ) +1
( 1)
= ( 1 )tr +1 + 9
2 +1
=s
4 2 +1
+1
+3
+1
tr
2 +1 1
4 2 +1 s
(45)
(46 ) (46 )
21
+1
and
where:
and
Utilizing Cramers Rule, the desired derivatives can be straightforwardly obtained as follows: ( ( (
+1
11
22
12
21 )
With all the preceding derivatives in hand, the fully implicit closest point projection algorithm for determination of 2 +1 is shown in Box 4. Once a trial value of 2 +1 is obtained (in the algorithm of Box 4), then
is nontrivial, even with a xed value of 2 +1 . Dening a residual vector based on Eqs. (46b) and (47a) as follows,
a Newtons method based on iterative linearization of the residual function is presented in Box 5
for updating
and .
r(
)=
( 1)
+1
( 1 ) +1 9 2 +1 + 3 ( + 1) 2 +1
(55)
updated. Due to the nonlinearity of the update Equations (46b) and (47a), the update of
and
must be
1
and
&
11
22 11
12
21 ) 21 )
(54 )
1 2 +1
22
12 )
(54 )
F=
&
9 3
+1
2 1
+1
1 + 3
(49 )
A=
2 1
2 +1
2 +1
21
22
2 +1
11
12
2 +1
+1
(52)
(53 )
22
update
compute ( 2 ) =
+1 +1
by Eq. (49)
+1
End while
given 2 +1 , ( 1 ) return
+1 ,
+1 ,
iteratively update ( 1 )
+1 +1
and
+1
+1
by Eq. (48)
+1
2 2
( 2)
+1 +1
= ( 2)
compute
While (
TOL2 ) then
+1
+1 +1
compute
2(
+1
( 1)
+1
tr 2
tr
(see Box 5)
( 1)
+1
= ( 1)
+1
+1
+1
+1
+1 )
=0
+1
(B4.1)
23
=0
let
= ( 1)
+1
and
+1
While r(x ) x =x
+1
TOL) then
update r(x = +1
+1
) by Eq. (55)
End while
Remark 1: For each iteration cycle of the algorithm in Box 4, once 2 is updated, exact updates of both 1 and are coupled and nonlinear, requiring a sub-iteration loop.
3.6
When the algorithm of Box 2 indicates that the tension cap surface is active (Case 3), or when the return point for Case 2 lies in the domain of the tension cap, then the stress must be returned to the tension cap surface. Since the tension cap has a xed radius, the Case 3 return mapping algorithm is considerably simpler than that for the compression cap (Case 2). Integration of the associated Case 3 ow rule gives the following plastic strain increment,
+1
+1
+1
+2
3 +1
+1
+3
tr
3 +1 1
4 3 +1
+1
+1
+1
+1
& $
= 3 +1
= 3 +1 2
r x
and
(B5.1)
and
for CASE 2.
(56)
(57 ) (57 )
24
Decomposition of the corrected stresses and back stresses into deviatoric and bulk components gives
+1
( 1 )tr +1 = 1 + 18 3 +1 The remaining objective of the return map algorithm is to enforce the plastic consistency condition 1 + 2(2 + ) 3 +1
2 2
The return map algorithm for Case 3 is necessarily iterative, and the Newtons Method algorithm of Box 3 is very efcient and robust. Once the Case 3 return map is successfully completed, the
=0
update
+1 +1
compute ( 3 ) =
+1
End while
+1
+1
update ( 1 )
+1 +1
by Eq. (59)
+1
3 3
( 3)
+1 +1
= ( 3)
+1
+1
+1 ,
compute
2 4(2 + ) 1+2(2 + ) 3
9 ( 1 9
)2 3
While (( 3 )
+1
TOL3 ) then
compute ( 3 )
+1
3(
+1
= ( 3 )tr +1 =
( 1)
( 1)
+1
= ( 1)
+1
+1
+1 +1 )
+1
+1
( 3)
1 +1 )
=0
+1
(59)
(B6.1) (B6.2)
+1
( 1)
18
= ( 1 )tr +1 + 9
3 +1
= ( 1 )tr +1
3 +1
1 + 2(2 +
) 3 +1
+1
(58 )
+1
+1
(58 )
25
In modern computational plasticity, it is now recognized [19] that in order to achieve the asymptotically quadratic rate of force-balance convergence that is theoretically possible with global Newton-Raphson force balance iterations, material tangent operators that are consistent with the implemented (discrete) form of the constitutive models must be utilized. This leads to so-called consistent tangent operators which generally differ quite signicantly from the continuum elastoplastic tangent moduli which can be derived from the rate form of the constitutive equations and the plastic consistency condition. Since the derivation of expressions for consistent tangent operators is conceptually straightforward [19], albeit algebraically complex, expressions for the consistent tangent operators for the three subcases of the smooth cap model are presented in the following sub-sections.
3.7.1
3.7.2
where =1 9
1 + 3 2 +1
3 2 +1
2 1
1 + 3
2 1
2 +1
(62 )
+1
2 1
27
( 2 +1 )2
= Celastic
(1
+1 )
1) +
(1
1)
+1 )
+1
+1
[I
+1 1 +1
2 1
2 1
+1
(2 +
)+
9 ( )2 1 9 1 +1
+1
=C
+1 ]
elastic
+1
+1
2 n
2 n
3 1 9 1 +1
3 1 9 1 +1
(60)
9 Idev
(61)
(62 )
26 (62 ) 2 +1 +
1
Inspection of the above consistent tangent operator proves it to be non-symmetrical arising from the terms with coecients
and
6.
cap parameter
symmetric consistent tangent only for associated ow rules and associated hardening laws. [Simo and Hughes, 1998] Hofstetter et al [18] developed an associated hardening law for the cap parameter and found that this did indeed lead to a symmetrical expression for the Case 2 consistent tangent
response characteristics of the model. Specically, in the implementation presented in [18], the Drucker-Prager envelope and the tension cap yield functions were not expressed as functions of , and consequently, loading on these yield surfaces, which results in plastic dilatance, does not lead to the usual retraction of the cap. For this reason, the associated hardening law proposed in [18] has not been adopted here. The consistent tangent operator expression above can be symmetrized and utilized. While the symmetrized consistent tangent operator is not precisely consistent with the Case 2 integration algorithm, it still gives much better performance in nonlinear nite computations than the classical elastic-plastic continuum tangent operator. It is proposed here that the consistent tangent operator be symmetrized by introducing the following modied constants:
8 2 2 +1 1 + 2(2 + ) 2 +1
2 2 +1
&
(62 ) (62 )
16
1 + 4(2 + ) 2 +1 (1 + 2(2 + ) 2 +1 )2
2 2 +1
=9
2 3 +1 +
2 4 +1
2 1
2 2 +1
(62 )
(62 )
12 1 + 2(2 + ) 2 +1
2 2 +1
12 1 + 2(2 + ) 2 +1
2 1
1 2 +1
(62 )
2 4 +1
(62 )
2 1
27
2 +1
10
+1
3.7.3
4.
EXAMPLE COMPUTATIONS
4.1
This simple test is designed to show the pressure versus volumetric strain behavior of the smooth cap model. A single nite element in Figure 5a is subject to a stress-controlled hydrostatic loading test, and the strain response of the element is computed and displayed in Figure 5b. To demonstrate the rate of convergence achieved with the consistent tangent operator expressions provided above, the convergence behavior for a single load-step during this test is shown in Figure 5, and shows an asymptotically quadratic rate of convergence. The cap model material parameters used in this test are listed in Table 1.
2 4 +1 (2 + ) 1+2(2 + ) 3 +1
9 ( 1 9
)2 3 +1
+1
3 +1
+1
4 1+2(2 + ) 3 +1
3 1 9
4 1+2(2 + ) 3 +1
3 1 9
+1
= Celastic
8 2 3 +1 9 2 Idev + 1 + 2(2 + ) 3 +1 1 9
3 +1 1 3 +1
+1 )
1)
+1 )
9 Idev
3 +1
= Celastic
(1 5
( 6
1) + (1 7
&
1 2
2 4 +1
where
2 +1
(64)
(65)
(66)
=9 7
2 3 +1
5 6
(63 )
= = 5 6
12 1 + 2(2 +
10 ) 2 +1
(63 )
28
10 8 Pressure (MPa) 6 4 2 1m 0 P z y x fixed P free to move x residual norm free to move x,y 1m 1m 0 0.005 0.01 volumetric strain 0.015
104 106 6
3 2 iteration number
Figure 5: 5a) Single element mesh for hydrostatic loading test; 5b) pressure versus volumetric strain response of the material model; and 5c) nite element force balance convergence characteristics of the model during a representative load step, for both the consistent tangent operator and the continuum tangent operator.
Material Parameter
Value 170.0 MPa 210.0 MPa -1.0 KPa 3.86 KPa 0.21 0.01
8
1.2E-6 Pa
29
While the preceding model demonstrated the performance of the model under purely hydrostatic loading, this simple four-element test is designed to show the models performance under combined deviatoric and compressive loading. In this test computation, the loading shown on the four-element mesh (Figure 6a) is increased until the limit state of the model is found using methods outlined in [20]. The computed load-displacement response of the model is shown in Figure 6b, and the material model parameters utilized are listed below in Table 2. In this test, which involved plastic loading at the crown of the cap at the limit state, very good convergence behavior was achieved.
Load
0 1
Displacement (m *104)
1m
1m
Figure 6: 6a) Simple plane-strain four-element mesh with applied loading and restraints; and 6b) the computed load-deformation response provided by the smooth cap model.
Parameter
Parameter
4 2
3.2E-7 Pa
30
In this relatively simple test problem, the bearing capacity of a shallow strip footing foundation
of width
= 2 24
solutions are compared with analytical results. The bearing capacity computations were performed using symmetry conditions, and the mesh utilized 600 bilinear continuum elements (Figure 7a).
The problem was solved twice, once for a loose sandy soil (
10
remaining material model parameters used for the loose and dense soils are identical as shown in Table 3. Before the rigid foundation loading was applied to the soil using displacement control, the soil deposits were rst allowed to consolidate under their own self-weight. Once the soils were properly stressed in this manner, the loose soil was normally consolidated due to the chosen value
of
0, 0.
while the dense soil would be considered over-consolidated, also due to the chosen value The computed foundation load versus displacement responses for these two soil conditions
of
Terzaghis theory with the assumption of a general shear failure. For the case of the loose sandy
0.1
0.1
Table 3. Sand-like material parameters used in Bearing Capacity Tests. Corresponding to = 0 15,
= 18 89 .
0 102
which would be computed using Terzaghis model with the assumption of a local shear
Value (Dense Sand) 2000 kg/m3 40.0 MPa 100.0 MPa -10.0 MPa 3.86 KPa 0.15 3.2E-8 Pa
1
are shown in Figure 7b. For the dense sandy soil, a bearing capacity of 0 268
was computed,
resting on a sandy soil is computed using the proposed cap model, and the
31
20m
2.0
2.5
3.0
30m
Load (MN)
7b) Computed load-displacement responses for loose sand (dashed) and compacted sand (solid).
Figure 7: 7a) Mesh of bilinear continuum plane-strain nite elements used in strip footing bearing capacity computations; and 7b) the computed load-displacement response for a rigid footing on a loose soil (dashed curve), and on a compacted soil (solid curve).
4.4 Slope Stability Analysis Computations
In these computations, an earthen slope model comprised of a sandy soil is analyzed for stability using the methods proposed in [20], which involve increasing the gravitational loading on the slope model until a failure mechanism develops, and the slope model can take no further loading. The mesh used to model the slope is shown in Figure 8 and contains 1130 bilinear continuum planestrain nite elements. The nite slope shown has a height of 30m and a repose angle of 29.98 . For a loose sandy soil and a dense sandy soil whose parameters are shown in Table 4, the computed stability factors of safety for this model were 0.95 and 0.98, respectively.
32
Figure 8: Mesh of 1130 bilinear continuum plane-strain nite elements used in slope-stability analysis.
Parameter
Value (Dense Sand) 1800 kg/m3 1.2 GPa 3.0 GPa -10.0MPa 30.0 Pa 0.2 0.07
8
0.07
5.
SUMMARY AND CLOSURE A smooth, three-surface cap model has been presented here along with sample computations
which demonstrate its very good performance characteristics. The model retains many of the positive physical attributes of preceding non-smooth cap models, but avoids many of the numerical difculties associated with corner points in those models. The numerical integration algorithms presented here deal with determination of active yield surfaces, iterative return mapping algorithms, and consistent tangent operator expressions. The rate of convergence achieved in the sample implicit nite element computations presented was typically asymptotically quadratic.
3.2E-8 Pa
4 2
33
Elastoplastic cap models of the type presented here are quite useful in modeling ductile soil behaviors. To capture softening behaviors in porous media which are also of considerable interest, cap models such as this one can be straightforwardly coupled with continuum damage mechanics models.
6.
ACKNOWLEDGEMENTS This research was funded by a University of Iowa Old Gold Fellowship, and by a research
grant from the Whitaker Foundation 96-0636. This support is gratefully acknowledged.
7.
BIBLIOGRAPHY
[1] Drucker, D. C. and Prager, W. Soil Mechanics and plastic analysis or limit design Q. Appl. Math., 10(2), 1952, pp. 157-165 [2] Drucker, D. C. Limit analysis of two- and three-dimensional soil mechanics problems J. Mech. Phys. Solids, 1, 1953, pp. 217-226 [3] Drucker, D. C., Gibson, R. E. and Henkel, D. J. Soil Mechanics and Work-Hardening Theories of Plasticity, Transactions, ASCE, Vol. 122, 1957, pp. 1692-1653 [4] Roscoe, K. H., Schoeld, A. N. and Worth, C. P. On the yielding of soils Geotechnique, 8(1), 1958, pp. 22-53 [5] Roscoe, K. H., Schoeld, A. N. and Thurairajah, A. Yielding of clays in state wetter than critical Geotechnique, 13(3), 1963, pp. 211-240 [6] Schoeld, A. N. and Wroth, C. P. Critical State Soil Mechanics McGraw-Hill, New York NY, 1968, pp. 310 [7] Burland, J. B. The yielding and dilation of clay Correspondence Geotechnique, 15(2), 1965, pp 211-214 [8] DiMaggio, F. L., and Sandler, I. S. Material models for granular soils, J. of Engng Mech., ASCE, Vol. 97. No. EM3. June 1971, pp. 935-950 [9] Chen, W. F. Plasticity in Reinforced Concrete McGraw-Hill, New York, 1982. [10] Chen, W.F. and Mizuno, E. Nonlinear analysis in soil mechanics: theory and implementation, Elsevier, 1990. [11] Desai, C. S. and Siriwardane Constitutive laws for engineering materials Prentice-Hall, 1984. [12] Bathe, K.J., Snyder, M.D., Cimento, A.P. and Rolph, W.D. On some current procedures and difculties in nite element analysis of elasto-plastic response, Computers & Structures, 12, 1980, pp. 607-624. [13] Sandler, I. S. and Rubin, D. An Algorithm and a modular subroutine for the cap model, Intl. J. Numer. Analy. Meth. Geomech., 3, 1979 pp. 173-186 [14] Resende, L. and Martin, J. B. Formulation of Drucker-Prager cap model J. Eng. Mech., 111 No. 7., 1985 pp. 855-881 [15] Simo, J.C., Kennedy, J.G. and Govindjee, S. Non-Smooth multisurface viscoplasticity: Loading/unloading conditions and numerical algorithms International Journal for Numerical Methods in Engineering, 26 2161-2185 (1988). [16] Ortiz, M. and Simo, J.C. An analysis of a new class of integration algorithms for elastoplastic constitutive relations, International Journal for Numerical Methods in Engineering 23, 1986, 353-366. [17] Simo, J.C., J-W. Ju, K.S. Pister and R.L. Taylor, Assessment of cap model: consistent return algorithms and rate-dependent extension, Journal of Engineering Mechanics, 114, 1988, 191-218. [18] G. Hofstetter, J.C. Simo and R.L. Taylor, A modied cap model: closest point solution algorithms, Computers & Structures 48-2 203-214 (1993).
34
[19] J. C. Simo and R. L. Taylor, Consistent tangent operators for rate-independent elastoplasticity, Computer Methods in Applied Mechanics and Engineering, 48 1985, pp. 101-118 [20] C.C. Swan and Y.-K. Seo, Limit state analysis of earthen slopes using dual continuum/FEM approaches, International Journal for Numerical and Analytical Methods in Geomechanics, 23 1359-1371 (1999). [21] Y.-K. Seo and C.C. Swan, Stability analysis of embankments on saturated soil deposits using elasto-plastic porous medium models, Journal of Getoechnical and Geoenvironmental Engineering 127-5 436-445 (2001). [22] O.C. Zienkiewicz, C. Humpheson and R.W. Lewis, Associated and non-associated viscoplasticity and plasticity in soil mechanics, Geotechnique 25 671-689 (1975).