Explicit integration scheme for generalized plasticity constitutive models with automatic error control

1222 EXPLICIT INTEGRATION SCHEME FOR GENERALIZED PLASTICITY CONSTITUTIVE MODELS WITH AUTOMATIC ERROR CONTROL MIGUEL M. STICKLE 1 *, PABLO DE LA FUENTE 2 AND CARLOS OTEO 3 1*: Applied Mathematics and Computer Science Department ETSI Caminos, Canales y Puertos Universidad Politécnica de Madrid Avd. Profesor Aranguren s/n, 28040 Madrid, Spain e-mail: miguelstick[email protected]pm.es 2: Continuum Mechanics and Structures Department ETSI Caminos, Canales y Puertos Universidad Politécnica de Madrid Avd. Profesor Aranguren s/n, 28040 Madrid, Spain e-mail: [email protected] 3: Professor on Ground Eng. C / Torpedero Tucumán 26, 28016 Madrid, Spain e-mail: [email protected] Key words: Generalized Plasticity, explicit integration. Abstract. An explicit algorithm for integrating Generalized Plasticity constitutive models is presented. This automatically divides the applied strain increment into subincrements using an estimate of the local error controlling the global integration error in the stress. The algorithm modifies the well known S. W. Sloan substepping scheme to account for Generalized Plasticity constitutive models, in which, unlike Classical Elastoplasticity, the yield surface is not explicitly defined. The integration scheme is described and results are presented for a rigid footing resting on a layer of specific Generalized Plasticity model for sands, in which a hyperelastic formulation is introduced to describe the reversible component of the soil response instead of the hypoelastic approach originally proposed. The explicit algorithm with automatic substepping and error control is shown to be reliable and efficient for these complex constitutive laws. 1 INTRODUCTION Nowadays it is well recognized that the selection of an adequate constitutive model, together with the use of accurate, efficient and robust integration algorithms of the elastoplastic equations, is a key point in finite element analysis of geotechnical problems. As observed by Hughes [1], the integration of the constitutive equations at the local level plays a crucial part in computational plasticity, since it strongly affects the performance of the constitutive equation in actual computations. Implementation of an advanced elastoplastic constitutive model into a finite element XI International Conference on Computational Plasticity. Fundamentals and Applications COMPLAS XI E. Oñate, D.R.J. Owen, D. Peric and B. Suárez (Eds) 1223 Miguel M. Stickle, Pablo De la Fuente, Carlos Oteo. 2 program requires the development of a robust and efficient numerical procedure in order to perform the integration of the constitutive equations along a given loading path. In the context of classical plasticity formulations, in which a yield function is defined in an explicit manner and the enforcement of the consistency condition is a key feature of the integration algorithm, a variety of implicit and explicit integration schemes might be used [2, 3]. On the other hand, based on the assumption that explicit integration schemes for highly non-linear models may potentially lead to inaccuracy and unstable behavior [4], implicit integration algorithms have been considered mostly in the context of non-standard elastoplastic models. The aforementioned assumption cannot be further supported, in the context of classical plasticity formulations, if explicit schemes are endowed with error control techniques [5, 6]. The same situation is observed for Generalized Plasticity based models [7]. The outline of the paper is as follows. We first present the fundamentals of Generalized Plasticity, with particular attention paid to the SandPZ constitutive equations including the modifications in the elastic component introduced by Mira and coworkers in 2009 paper. A novel explicit algorithm for integrating Generalized Plasticity constitutive models is presented in the following section. Finally, Results and conclusions are presented for a rigid footing resting on a sand layer modeled by a SandPZ constitutive relation. 2 GENERALIZED PLASTICITY FRAME WORK. MODIFIED SAND PZ MODEL. The Generalized Plasticity basic idea, introduced by Zienkiewicz and Mroz [8] later extended by Pastor and coworkers [9, 10] and also by Mira and coworkers [4], is that no yield neither plastic potential surface are explicitly defined, but the gradients to the functions themselves. The elastoplastic behavior of the material within the Generalized Plasticity theory is described by the general incremental relationship, ( ) : ep ep ij ijkl kl d Dd d d σε ′′ =⋅ = σDε (1) In which the tangent elastoplastic stiffness four order tensor ep D depends not only on the internal state variables but also on the current effective stress state ′ σ , on the strain-stress history and the direction of the effective stress increments d ′ σ . The dependence of ep D on the direction of d ′ σ is expressed by simply distinguishing between two different loading classes, namely Loading (L) and Unloading (U). Therefore a normalized direction n is defined in the effective stress space for any given ′ σ , determining loading/unloading/neutral loading condition. There are two possibilities for the tangent elastoplastic stiffness tensor in (1) depending on whether loading, ep L D or unloading ep U D is occurring. To guarantee continuity between loading and unloading, ep L D and ep U D are defined as ( ) () [] ( ) () [] 11 11 1 1 ep e L L L ep e U U U H H −− −− = +⋅⊗ = +⋅⊗ D D mn D D mn (2) 1224 Miguel M. Stickle, Pablo De la Fuente, Carlos Oteo. 3 In expression (2), L m and U m are directions of unit norm representing the plastic flow direction in loading (L) and unloading (U) conditions respectively, L H and U H are two scalar functions defined as plastic moduli while e D is the tangent elastic stiffness tensor. By suitable manipulation of (2), ep L D and ep U D can be obtained giving: :: :: :: :: ee ep e L Le LL ee ep e U Ue UU H H ⊗ =− + ⊗ =− + D m nD DD nD m D m nD DD nD m (3) The strain increment d ε can be decomposed into elastic and plastic parts as ep dd d =+ εε ε where ( ) 1 : 1 : for loading 1 : for unloading ee p L L p U U dd dd H dd H − ′ =  ′ = ⋅⊗    ′ = ⋅⊗   εDσ εmnσ εmnσ (4) Therefore, in a Generalized Plasticity approach, the non-linear irreversible behavior of soils can be fully described by simply specifying three directions, , L nm and U m , two scalars, L H and U H and a fourth order tensor e D . Since the hardening moduli L H and U H as well as the plastic flow directions L m and U m are fully determined without reference to any yield surface nor plastic potential, different expressions can be selected for them whether the stress increment implies loading or unloading. Moreover consistency cannot be enforced and the consistency parameter d λ is simply defined as: :: :: e e LU LU d dH λ =+ nD ε nD m (5) Although not explicitly defined, plastic potential and yield surface can be established a posteriori, by integrating LU m and n , respectively. SandPZ model was developed by Pastor and coworkers [9] as a particular type of Generalized Plasticity formulation with the aim of predicting granular soil behavior under both monotonic and cyclic loading. The model assumes an isotropic material response. As a result, the plastic flow direction m , as well as the loading direction n , is expressed in the invariant space defined by ,, pq θ ′ as 1225 Miguel M. Stickle, Pablo De la Fuente, Carlos Oteo. 4 mmm vs ij ij ij pq θ θ σσσ ′ ∂∂∂ =++ ′′′ ∂∂∂ m (6) The value of the coefficients m ,m ,m vs θ are loading class (loading or unloading) dependent. In order to take into account the main features of sand response, i.e. the existence of a critical state condition, dilative response after peak, liquefaction in loose sands, memory of previous stress path, Pastor and coworkers [9] proposed for the plastic modulus L H the following relationship: ( ) 0 L f v s DM H H pH H H H ′ = ⋅⋅ ⋅ + ⋅ (7) Together with () 4 01 0 1 max 1 ; 1 ; exp 1 ; ; 1 1 f f fvs ff g f p Dm s ff H HH MM H dd p M γ α α ηη ββ βξ α α ζη ξ ε ξζ ζα −  = − ⋅ =− = −   +     ′ = ===−⋅     +     ∫∫ In these expressions 001 ,,, H ββγ are constitutive parameters, ξ is the accumulated deviatoric plastic strain and max ζ stands as the maximum value of the mobilized stress function ζ accounting for the soil stress history. In the case of unloading the plastic modulus U H is given by: 0 0 for 1 for 1 u gg Uu uu g Uu u MM HH M HH γ ηη η  =>   =< (8) For where 0 u H is a constitutive parameter and u η , referred as unloading stress ratio, is the stress ratio qp ′ from which unloading takes place. Finally, the PZ model assumes a non-linear elastic response of the soils. As in a large number of constitutive models, the non-linear reversible behavior is described through a hypoelastic approach, in which the tangent bulk modulus K and shear modulus G only depend on the hydrostatic part of the effective stress tensor, according to the following relationships 00 00 , pp KK GG pp ′′ =⋅ =⋅ ′′ (9) Although widely used, one of the major shortcomings of such hypoelastic formulation is that it results in a non-conservative elastic response and energy dissipation over closed stress 1226 Miguel M. Stickle, Pablo De la Fuente, Carlos Oteo. 5 paths [11]. An alternative is to describe the elastic response of soils within a conservative framework adopting the hyperelastic approach based on the existence of an energy potential from which the reversible response can be derived. This naturally leads to a conservative elastic response, guaranteed to obey the First Law of Thermodynamics, and thus avoiding the problems on cycling described above [12, 13] Among the different formulations recently proposed in the geotechnical literature, in this work the hyperelastic approach described by Houslby and coworkers in 2005 (from now on referred as HAR according to authors’ initials) has been adopted to describe the reversible component of the soil response. This Generalized Plasticity PZ model for granular soils has been proposed firstly by Mira and coworkers [4] under its triaxial formulation, extending the model in the present work to deal with a general stress formulation. The stored energy function ϒ of the HAR model in general stress formulation has two different expressions depending on the value assigned to the dimensionless pressure exponent HAR n , which governs the amount of nonlinearity involved in the formulation. For 1 HAR n ≠ , ϒ takes the following form: () ()() ( )( ) 21 0 1 2 HAR HAR nn ea ij HAR HAR HAR HAR pkn kn ευ −− ϒ = ⋅ ⋅− ⋅  ⋅− (10) Where () () () 2 0 2 11 1 11 ee ij ij HAR ee ii jj HAR HAR HAR HAR HAR HAR g ee kn knkn υε ε             =+ ⋅+ + ⋅− ⋅− − , while the asymptotic expression for 1 HAR n = is ( ) () a HAR e ee ii ij ij HAR HAR HAR k g k ee e ij pk e ε ε    ⋅+ ⋅ ⋅ ϒ= ⋅ . , HAR HAR kg are dimensionless constants representing the shear and bulk stiffness factors, respectively, while a p is the atmospheric pressure, adopted as reference stress. The effective stress tensor ′ σ and the tangent elastic tensor e D can be unambiguously determined by taking the first and second order derivatives of (10), obtaining the following expression () 0 2 0 1 12 3 HAR n ij kl e ijkl a HAR HAR HAR HAR ij kl HAR ik jl kl ij a p D p nk k n g pp σσ δδ δδ δδ ′′     =⋅ ⋅ ⋅ + − + −         (11) Where () ( ) 2 0 1 92 HAR HAR mn mn mm nn HAR k n ss pg σσ ⋅− ⋅ ′′ =+ For the present model there are 12 material parameters requiring definition. Generally, all parameters are identified from monotonic and cyclic triaxial tests, though in certain cases some parameters are adopted from previous experiences if full test records are unavailable. 3 EXPLICIT INTEGRATION OF GENERALIZED PLASTICITY MODELS. During a typical step or iteration of an elastoplastic finite element analysis, the forces are applied in increments and the corresponding displacement increments are found from the 1227 Miguel M. Stickle, Pablo De la Fuente, Carlos Oteo. 6 global stiffness equations. Once the nodal displacement increments d u are known, the strain increments at a discrete number of integration points within each element are determined using the strain-displacement relation dd = ε Bu . If the stresses associated with an imposed strain increment cause plastic yielding, it is necessary to solve the system of firs order ordinary differential equations (12)-(13): :: ep ep t ∆  ′==  ∆  ε σDDε ɺɺ (12) p LU d λ =⋅εm ɺ (13) Where :: :: ee LU ep e e LU LU H ⊗ =− + D m nD DD nD m (14) :: :: e e LU LU d dH λ =+ nD ε nD m (15) In these expressions, ′ σ denotes the effective stress tensor, ε the small strain tensor and p ε the plastic strain tensor. The superior dot represents a derivative with respect to time while t ∆ is the time interval over which the external forces have been applied. Since the effective stress and the plastic strains are kwon at the beginning of the time interval, and the known strain rates may be assumed to be constant through the time interval with value t ∆∆ ε , the equations (12) and (13) define an initial value problem. In order to integrate these equations numerically, it is convenient [14] to introduce a pseudotime, T , defined by ( ) 0 T tt t =− ∆ , where 0 t is the time at the start of the load increment, while 0 tt +∆ is the time at the end of the load increment, with 01 T ≤≤ . Since 1 dT dt t =∆ application of the chain rule to ′ ɺ σ and p ɺ ε in (12) and (13) gives :: : :: :: ee LU ep e ee LU e LU LU d dT H λ  ⊗ ′′ = ∆= − ∆=∆ −∆   +  D m nD σD εDεσ Dm nD m (16) p LU d dT λ =∆ ⋅ ε m (17) Where :: :: e e LU LU H λ ∆ ∆= + nD ε nD m (18) Equations (16) and (17) define a classical initial value problem which needs to be 1228 Miguel M. Stickle, Pablo De la Fuente, Carlos Oteo. 7 integrated over the pseudo time interval from 0 T = to 1 T = , where the known values known values in these relations are the imposed strain increments, ∆ ε , together with the effective stresses and plastic strain at the start of the pseudo time increment. The quantities / and LU mn are effective stress functions, while parameter / LU H is a function of both the effective stress and the plastic strain. In order to solve the system of first order ordinary differential equations (16)-(17) Sloan developed a substepping algorithm [5] where the constitutive law is integrated by automatically dividing the strain increment into a number of substeps. An appropriate size for each substep is found through the use of modified Euler or Runge-Kutta-Dormand-Prince formulae, which are specially constructed to provide an estimate of the local error. Later, Sloan and coworkers [6] generalized the 1987 scheme, incorporating new algorithms for handling elastoplastic unloading, computing the yield intersection point, and restoring the stress to the yield surface. Sloan and coworkers schemes were developed exclusively for classical plasticity based models, including classical and generalized critical state models, where non-linear elastic behavior inside the yield surface is exhibited. In all these models the admissible states in the stress space are constrained to lie within the interior or the boundary of the domain explicitly defined by the yield surface. As this is not the case for generalized plasticity based models, Sloan substepping algorithm should be adjusted in order to be able to integrate this kind of models. The proposed integration scheme starts with the known strain increment, ∆ ε , the initial stress 0 σ and initial plastic strain 0 p ε at the start of the increment where 0 T = and 0 tt = . At the end of the integration process the stresses and plastic strains are obtained at the end of the increment where 0 T = and 0 tt = . Consider a pseudo time subincrement in the range 01 n T ≤∆ ≤ and let the subscripts 1 n − and n , denote quantities evaluated at the pseudo times 1 n T − and 1 nn n TT T − = +∆ , respectively. Plastic modulus and plastic flow direction in expressions (16)-(17) are dependent on the direction of the effective stress increments therefore differentiation between the two loading classes should be performed before the proper integration process starts. By means of the strain subincrement nn T∆ =∆ ∆ εε the loading class is firstly established through the following expression ( ) ( ) :: e n nn ′′ ∆ n σDσε (19) If expression (19) is positive an elastoplastic loading process is performed, if negative an elastoplasctic unloading process is implied, while an elastic process comes from a zero value. In the explicit Euler method, the solution for , p ′ σε at the end of the pseudo time step n T ∆ is found from 11 11 nn pp p nn − − ′′ ′ = +∆ = +∆ σσ σ εε ε (20) Where 1229 Miguel M. Stickle, Pablo De la Fuente, Carlos Oteo. 8 ( ) () () 1 11 1 11 1 ,: ,, ep p nn n pp n n n LU n λ −− −− − ′′ ∆= ∆ ′ ∆ =∆ ∆ ⋅ σDσε ε ε σε εmσ (21) A more accurate estimate of the stress and plastic strains at the end of the interval n T ∆ can be found using the modified Euler procedure, which is given by () () 1 12 1 12 1 2 1 2 nn pp p p nn σσ εε − − ′′ ′ ′ = + ∆ +∆ = + ∆ +∆ σσ εε ⌢ ⌢ (22) Where 11 and p ′ ∆∆ σε are obtained from Euler scheme and ( ) () () 2 1 11 1 2 1111 11 ,: ,, ep p p nn n ppp n n n LU n λ −− −− − ′ ′′ ∆ = +∆ +∆ ∆ ′′ ′′ ∆ =∆ +∆ +∆ ∆ ⋅ +∆ σ D σ σε ε ε ε σ σε ε ε m σσ (23) Since the local truncation error [15] in the Euler and modified Euler solutions is ( ) ( ) 23 and OT OT ∆∆ , respectively, the error in n σ and p n ε can be estimated from () () 21 21 1 2 1 2 nn pp pp nn  ′′ ∆ −∆  ′′    −=         ∆ −∆   σσ σσ εε εε ⌢ ⌢ (24) Using any convenient norm, this quantity can be used to compute the relative error measure 21 21 1max , 2 pp np nn R  ′′ ∆ −∆ ∆ −∆  =    εε σσ σε (25) Following 1987 Sloan work, the current strain subincrement is accepted if n R is not greater than some prescribed tolerance, STOL , and rejected otherwise. Regardless of whether the subincrement is accepted or rejected, the next pseudo time step is found from the simple relation 1 nn T qT + ∆ = ⋅∆ (26) where q is chosen so that 1 n R + satisfies the constraint 0.8 , 0.1 1.1 n q STOL R q ≤ ≤≤ (27) 1230 Miguel M. Stickle, Pablo De la Fuente, Carlos Oteo. 9 Two typical controls are finally incorporated. A minimum absolute step size, min T ∆ , and a step size is not allowed to grow immediately after a failed subincrement. 4 RESULTS AND CONCLUSIONS. The behavior of a smooth rigid strip footing resting on an elastoplastic soil mass, governed by the modified SandPZ model presented above, is consider in order to analyze the performance of the proposed integration scheme. Due to the singularity at the edge of the footing and the strong rotation of the principal stresses, this example is a good test for assessing the integration strategy. As loading is prescribed in the form of displacements, an equivalent uniform pressure is found by summing the appropriate nodal reactions. To assess the accuracy of the scheme, an estimate of the stress integration error is found directly from 2 2 ref error ref σ ′′ − =′ σσ σ (28) Where ′ σ are the effective stresses obtained by the proposed integration scheme, ref ′ σ are the reference effective stresses while 2 i is the Euclidean norm. The reference effective stresses are obtained by the explicit Dormand-Prince integration scheme with a stress tolerance of 9 10 STOL − = . Note that the reference stresses provide a very accurate set of stresses for the given mesh and loading sequence and all values are computed at the end of the last load increment. The results for the analyses with 10 load increments of equal size are presented in Table 1. It can be observed how the uniform pressure over the footing after applying 4mm of vertical displacement is similar for all of the specified stress tolerances with values varying by less than 0.6% of the reference pressure. Table 1: Smooth rigid strip footing on SandPZ layer. 10 load steps. Stress Tolerance Equivalent uniform pressure after a vertical displacement of 4mm [N/m 2 ] % of the equivalent pressure obtained under Dormand-Prince integration scheme error σ 2 10 STOL − = 115320 0.14% 3 1.7 10 − ⋅ 4 10 STOL − = 115450 0.02% 4 1.7 10 − ⋅ The error in the computed stresses for the proposed scheme, as defined by equation (28), is less than the integration tolerance for 2 10 STOL − = and within the order of magnitude for 4 10 STOL − = . Therefore the tolerance STOL thus gives a required error control. Figure 1 shows the error spatial distribution induced by the proposed local integration