INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) NUMERICAL ASPECTS OF A NON-PROPORTIONAL CYCLIC PLASTICITY MODEL UNDER PLANE STRESS CONDITIONS STEFAN HARTMANN1;∗ , MARC KAMLAH 2 AND ANDREAS KOCH 3 2 1 Institut f ur Mechanik; Universitat Gh Kassel; Monchebergstrae 7; 34109 Kassel; Germany Forschungszentrum Karlsruhe; Institut fur Materialforschung II; Postfach 3640; 76021 Karlsruhe; Germany 3 Institut f ur Mechanik (Bauwesen); Universitat Stuttgart; Pfaenwaldring 7; 70550 Stuttgart; Germany This paper is dedicated to Professor Peter Haupt on the occasion of his 60th birthday ABSTRACT The objective of this paper is the investigation of the inuence of a material model representing nonproportional loading conditions with respect to cyclic plasticity phenomena on several plane boundary-value problems and the development of the corresponding stress algorithm. This material model, developed by Haupt and Kamlah,1 contains linear isotropic elastic behaviour, a von Mises yield function, an associated ow rule, non-linear kinematic hardening of Armstrong and Frederick type, a modied arc-length representation considering cyclic plasticity phenomena and the inclusion of non-proportional hardening eects. These rate-independent constitutive equations are based on the assumption of small strains. The boundary-value problem will be solved by the nite element method including investigations of a semi-analytical computation of the consistent tangent operator. Concluding examples will show non-proportional hardening eects as well as inhomogenization phenomena stated by Luhrs and Haupt 2 for specimen under uniaxial cyclic loading. ? 1998 John Wiley & Sons, Ltd. KEY WORDS: nite elements; cyclic plasticity; non-radiality; non-linear kinematic hardening; stress algorithm INTRODUCTION The detailed representation of material phenomena in metal plasticity involves rather complicated constitutive models. The increase in the number of evolution equations leads to the diculty of developing simple and cheap stress algorithms in the context of nite element analysis. This is reinforced by the application of cyclic plasticity. A widely used rate-independent plasticity model for small strains deals with a linear isotropic elasticity relation, a von Mises ow function, an associated ow rule and special evolution equations for the hardening phenomena. The consideration of cyclic processes leads to the incorporation of more complicated kinematic hardening equations. We choose a kinematic hardening rule proposed by Armstrong and Frederick3 which has been enhanced by a transformed arc-length by Haupt et al.4 This model has the advantage that only one further scalar evolution equation has to be introduced to represent most uniaxial cyclic plasticity phenomena such as cyclic hardening and softening as well as, for example, the independence of the mean strains and cyclic mean stress relaxation. The basic idea consists of ∗ Correspondence to: Stefan Hartmann, Institut f ur Mechanik, Universitat Gh Kassel, Monchebergstrae 7, 34109 Kassel, Germany. E-mail: hart@ifm.maschinenbau.uni-kassel.de CCC 0029–5981/98/081477–22$17.50 ? 1998 John Wiley & Sons, Ltd. Received 4 September 1997 Revised 6 February 1998 1478 S. HARTMANN, M. KAMLAH AND A. KOCH the introduction of a generalized arc-length. The implementation of this model in a nite element code was performed by Luhrs.5 Several materials exhibit a signicantly enlarged stress response to non-proportional biaxial loading histories compared to a uniaxial loading history of the same equivalent strain amplitude (e.g. References 6 –8). In the past, investigations concerning the experiments and the constitutive modelling of biaxial loading cases were enforced. A very simple approach describing non-proportional loading conditions is due to Haupt and Kamlah,1 where two further evolution equations are introduced: the rst one governs a measure of non-proportionality based on the angle between the current backstress and their corresponding rate introduced by Benallal et al.9 The other internal variable represents the evolution of the yield stress, i.e. isotropic hardening. This model has the advantage of employing no additional memory surfaces and thus the evolution of the internal variables is continuous. In this paper we present the implementation of this material model into a nite element code based on the assumption of plane stress conditions. The applied implicit stress algorithm, which is an extension of the algorithms proposed in Hartmann et al.10 and Luhrs and Haupt,2 can be reduced to the solution of three non-linear equations. These equations are solved by a Newtonlike method. The corresponding consistent tangent operator, which is needed to achieve quadratic convergence of the global equilibrium iteration, might be computed either fully numerically or semi-analytically. To show the inuence of a material model including non-proportional hardening eects, several examples are investigated. The rst example deals with inhomogenization eects in uniaxially loaded specimens under cyclic loading conditions (cf. Luhrs and Haupt) and the other describes the benchmark example of the plane strip with a circular hole. CONSTITUTIVE MODELLING Usually, constitutive equations are motivated by considering uniaxial stress–strain curves of tensile or torsional tests. The three-dimensional generalization of these material models is then carried out by tensors and tensor invariants retaining the structure of the uniaxial equations. Since the tests of Lamba and Sidebottom6 it is known that a dierent hardening behaviour occurs during non-proportional biaxial or multiaxial loading conditions. Most of these experiments are carried out by tension=torsion loadings of thin-walled cylindrical specimens. Several works dealt with the modelling of these hardening eects (see, for example, References 7, 8, 11–15 and the literature cited there). A review has been given by Ohno.16 Most of these models are based on von Mises plasticity or Perzyna-Chaboche-type viscoplasticity models, which are extended by further evolution equations. The geometrically linear model developed by Haupt and Kamlah1 is built up as follows: we assume an additive decomposition of the linearized strain tensor E = 12 (grad u+grad T u), u = u(x; t) representing the displacement eld, into elastic strains E e and plastic strains E p : E = Ee + Ep (1) The elastic strains are dened by the linear, isotropic elasticity relation T = K(tr E e )I + 2GEDe = CE e : T denotes the stress tensor and C = KI ⊗ I + 2G I − 13 I ⊗ I ? 1998 John Wiley & Sons, Ltd. (2) (3) Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) NUMERICAL ASPECTS OF A CYCLIC PLASTICITY MODEL 1479 the fourth-order elasticity tensor. K = E=(3(1 − 2)) denes the bulk modulus and G = E=(2(1 + )) the shear modulus. E and are the Young’s modulus and the Poisson ratio. tr A = A · I denes the trace operator. The point between two second-order tensors symbolizes the inner product, A · B = aij bij . The superscript D describes the deviator operator, AD = A − 13 (tr A)I. I stands for the identity tensor of second order and I for the fourth-order counterpart. With the introduction of the von Mises yield function f(T; X; k) := 12 (T − X)D · (T − X)D − 13 k 2 ; (4) yield and loading conditions are usually dened: elasticity f¡0 or f = 0 and ḟĖp = 0 ¡0 elastoplasticity f = 0 and ḟĖp = 0 ¿0 (5) ḟĖ denes the time derivative of the yield function f at constant internal variables. Since p all evolution equations depend on the plastic strain rate Ėp it is equivalent to demand Ėp = 0. X describes the backstress tensor which corresponds to the kinematic hardening behaviour and k represents the isotropic hardening parameter. The plastic strains evolve in the case of plastic loading according to the normality rule Ėp = (T − X)D = N k(T − X)D k (6) with the normal tensor N= (T − X)D k(T − X)D k (7) p and the norm kAk = tr (AT A). The plastic strains exhibit plastic incompressibility, tr Ėp = 0, caused by their deviatoric evolution. The superscript T denotes the transposition. is the plastic multiplier, which can be calculated from the consistency condition (ḟ = 0). The kinematic hardening rule, attributed to Reference 3, Ẋ = cĖp − bżX (8) is modied by Haupt et al.4 by an arc-length z dened by ż = with the accumulated plastic strain r ṡ = ṡ 1 + ap 2 Ėp · Ėp = 3 (9) r 2 3 (10) c; b and a are material parameters and p is a further internal variable. The introduction of the modied arc-length z makes it possible to describe essential uniaxial properties under cyclic loading conditions, namely cyclic hardening and softening eects, independence of the mean strain and cyclic mean stress relaxation. This is guaranteed by a further evolution equation for the parameter p, ṗ = ? 1998 John Wiley & Sons, Ltd. ṡ (kXk − p) sp (11) Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1480 S. HARTMANN, M. KAMLAH AND A. KOCH The quantity p can be interpreted as a measure of the plastic strain amplitude.4; 17 sp denes a material parameter. The extension of this model describing non-proportional hardening behaviour is proposed in Haupt and Kamlah.1 To this end, two further constitutive equations are introduced: the rst one, governs the measure of the non-proportionality with the material parameters s and ˆ : ˆ ˙ = ṡ (| sin |ˆ − ) = ṡ ((1 − cos2 )=2 − ) s s (12) Thus, is calculated from the history of the non-proportionality angle introduced by Benallal et al.:9 cos = X·Ẋ kXkkẊk (13) The current values of measures the average non-proportionality of the recent process history. In the case of initial loading (X = 0) is assumed to be zero. The quantity inuences in the proposed model the isotropic hardening behaviour 1 k̇ = ṡ k0 − k ; k(0) = k0 (14) 1 + p in which k describes the second internal variable. k0 denes the initial uniaxial yield stress. Note, that no isotropic hardening takes place under uniaxial loading conditions. The whole model contains the material parameters E and for the elasticity relation, c, b, a, sp and k0 for the proportional part of the elastoplastic material model and s , ˆ , and for non-proportional loading conditions. These equations could be integrated by evaluating the consistency condition ḟ = 0 to calculate the plastic multiplier occurring in (6) or (10). In the context of the nite element method, using a return mapping scheme, another procedure is favourable. WEAK FORMULATION Implicit non-linear nite element formulations resulting from a weak formulation, e.g. the principle of virtual work or other displacement driven formulations like the strain projection method (see, for example, Hughes18 ) or the enhanced assumed strain method (see Simo and Rifai19 ), are solved mostly in two nested procedures. The outer algorithm iteratively solves the nodal displacements under the agency of external loads. Usually, a Newton–Raphson method is applied. The inner procedure computes the stresses and the concerning internal variables by the so-called stress algorithms. To demonstrate the aforementioned procedure, we choose for simplicity the principle of virtual displacements. The external loads will be applied in a nite number of load steps, denoted by an index n characterizing the time t (n) . The dierence of the virtual work of the external loads and the stresses at time t (n) will be assumed to vanish (u(n) ; u) = 0 with (n) (u ; u) := ? 1998 John Wiley & Sons, Ltd. Z V E · T (n) Z dV − A u · t (15) (n) Z dA − V u · k(n) dV (16) Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) NUMERICAL ASPECTS OF A CYCLIC PLASTICITY MODEL 1481 for a material body which is subjected to the equilibrium condition. T(n) describes the stress state at time t (n) , E = 12 (grad u+grad T u) are the corresponding virtual strains, which are calculated from the symmetric part of the gradient of the virtual displacements u (u = 0 on Au where the displacements u(n) = u(x; t (n) ) are known, u(n) = u on Au ), t(n) represents the surface tractions on the material surface A (A = A ∪ Au , but A ∩ Au = ∅) and k(n) denes the specic volume force with the mass density . The external loads are applied in a nite number of steps, for two reasons: rstly, the Newton– Raphson method for the equilibrium iteration is only locally convergent and one has to be close by the solution. The second motive arises when considering the time integration of the evolution equations, since the accuracy of the time integration procedures depends on the step-size t = t (n+1) − t (n) . To this end, we look for the equilibrium state at time t (n+1) (u(n+1) ; u) = 0 (17) on the basis of the solution of equation (15). The application of the Newton–Raphson method and the nite element specic discretization lead to a linear system of equations in each iteration step denoted by the index i: K iT U = F(n+1) − R i (18) R K iT = V BT CiL B dV stands for the tangential stiness matrix, U = Ui+1 − Ui represents the increR R ment of the nodal displacements, F(n+1) = A NT t (n+1) dA + V NT k (n+1) dV denes the external R forces and R i = V BT Ti dV the residual force vector. These matrices are assembled from the element matrices usually written with a superscript e. Here, we omit this index for the sake of brevity. CiL represents the so-called consistent tangent operator resulting from a correct derivation of the derivative CiL = dTi dEi (19) (see References 20–22 and 10) and Ti are the stresses computed by the stress algorithm. N represents the matrix of the shape functions and B stands for the strain–displacement matrix, Ei = B Ui . The dimension of the aforementioned matrices depends on the element formulation. In the three-dimensional case the stresses and strains are stored in six-dimensional vectors, Ti ∈ R6 and Ei ∈ R6 , and the consistent tangent operator in a quadratic matrix, CiL ∈ R6×6 . In the case of plane stress conditions the matrices have the dimensions Ti ∈ R3 , Ei ∈ R3 and CiL ∈ R3×3 . STRESS ALGORITHM UNDER PLANE STRESS CONDITIONS To investigate the inuence of the proposed model on several boundary conditions, we choose plane stress conditions. This approach states that the out-of-plane normal and shear stresses in that direction are zero (33 = 31 = 23 = 0). This assumption denes constraints for the resulting boundary-value problem, which are incorporated in advance in the representation of the tensors in the constitutive equations. Both Simo and Taylor 23 and Simo and Govindjee 24 investigate their resulting return mapping algorithm used in the context of von Mises plasticity with kinematic hardening of a special type, i.e.: Ẋ = c(s)Ėp . The special structure of this hardening rule is exploited in the resulting stress algorithm based on a spectral decomposition representation. However, the ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1482 S. HARTMANN, M. KAMLAH AND A. KOCH application of this idea does not lead to an eective algorithm in the context of the constitutive framework of this paper. In a rst step we reformulate the constitutive equations (1) – (14), thereby obtaining the following matrix representation under plane stress conditions: T = C E e = C{E − E p } Ė p = P 1 N ( Ẋ = cN − r ṗ = (20) (21) r 2 b X 3 1 + ap ) 21 (kXkPSC − p) 3 sp (22) (23) r 2 1 ˆ ((1 − cos2 )=2 − ) 3 s r 1 2 k̇ = k0 − k 3 1 + p ˙ = f(S; X; k) = 12 {S − X}T P 2 {S − X} − 13 k 2 = 12 nT P 2 n − 13 k 2 cos = with 11 T = 22 ; 12 XT P 2 Ẋ kXkPSC kẊkPSC X11 11 e11 p11 X = X22 ; E = 22 ; E e = e22 ; E p = p22 X12 12 e12 p12 211 − 22 s11 1 n −11 + 222 ; n = S − X; N = q S = P 3 T = s22 = 3 T n P 2n s12 312 1 0 2G 1 0 C= 1− 1− 0 0 2 2 −1 0 1 2 1 0 1 −1 ; 1 2 0; 2 0 1 P = P = P1 = 2 3 3 0 0 3 2 0 0 2 ? 1998 John Wiley & Sons, Ltd. (24) (25) (26) (27) (28) (29) (30) (31) Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1483 NUMERICAL ASPECTS OF A CYCLIC PLASTICITY MODEL Single-underlined characters denote column vectors and those with double-underlining matrices. Scalar values are not underlined. In the same manner we distinguish parantheses: ( ) are used for scalars, { } for column vectors and [ ] for matrices. in the usual manner, ij = 2ij The shear strains ij , i 6= j, are related to the tensor components √ for i 6= j (see equation (28)3;4;5 ). The norm of a tensor A, kAk = A · A, has to be adjusted to the vector representation, occuring under plane stress conditions (PSC), q (32) kAk 7→ kAkPSC = AT P 2 A; in which the inner product A · A changes to AT P 2 A. The constraint 33 = 0 leads to e33 = −=(1 − )(e11 + e22 ) and the evolution equations for the plastic strains as well as the backstresses yield p33 = −(p11 + p22 ) and X33 = −(X11 + X22 ), respectively. Therefore, the strains 33 = e33 + p33 are known if the plastic strains p11 ; p22 and the stresses 11 ; 22 have been calculated: 33 = (−(11 + 22 ) − (1 − 2)(p11 + p22 ))(1 − )−1 . According to the well-known elastic predictor–plastic corrector scheme the plastic strains are held constant and the deviatoric trial state, using elasticity relation (20), T S = P 3 T T = P 3 C{Ei − E(n) p } (33) is computed. Inserting this stress state in the yield function (26), f(T S; X(n) ; k (n) ), yields an algorithmic counterpart of the loading criterion (5). If f(T S; X(n) ; k (n) )60 holds, the stress state lies in the elastic region and the new state variables for given strains Ei are known: Ti = T T; Eip = E(n) p ; Xi = X(n) ; pi = p(n) ; i = (n) ; k i = k (n) (34) If the trial stress state corresponds to f(T S; X(n) ; k (n) )¿0, the evolution equations of the internal variables are integrated with a Backward–Euler step under the condition f(Si ; Xi ; k i ) = 0 (this is equivalent to the method applied to dierential algebraic equations, in which f = 0 denes the algebraic constraint): Si = P 3 Ti = P 3 C{Ei − Eip } i i Eip = E(n) p + P 1N ( Xi = X(n) + i cN i − r i p =p (n) + r i r (35) 2 b Xi 3 1 + api ) 21 (kXi kPSC − pi ) 3 sp (n) = + i f(Si ; Xi ; k i ) = 12 {Si − Xi }T P 2 {Si − Xi } − 13 k i T 2 = 12 ni P 2 ni − 13 k i = 0 ? 1998 John Wiley & Sons, Ltd. (37) (38) 2 1 ˆ ((1 − cos2 i )=2 − i ) 3 s r 1 2 i (n) i i k =k + k0 − k 3 1 + pi i i (36) (39) (40) 2 (41) Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1484 S. HARTMANN, M. KAMLAH AND A. KOCH Here, we introduce the algorithmic counterpart of the plastic multiplier i := i t and the dierence tensor ni := Si − Xi . The cosine of the angle i is dened by the time derivative of the backstress tensor Ẋ, which is replaced by the increment Ẋ = {Xi − X(n) }=t, T i cos = Xi P 2 {Xi − X(n) } kXi kPSC kXi − X(n) kPSC (42) Equations (35)–(41) represent a system of thirteen non-linear equations with thirteen unknowns (Si , Eip Xi , pi , i , k i , i ). To transfer this set of equations to a more eective algorithm several conversions can be carried out, which lead to a reduction to only three non-linear equations. For the non-linear kinematic hardening model this conversion is known10; 25; 26 and leads to one non-linear equation to calculate i . For the extension of the model to cyclic processes, i.e. the introduction of pi in the stress algorithm, Luhrs and Haupt2 yielded a further equation, dominated by pi . The enhancement of the model taking into account non-proportionality exhibits a further equation, calculating the yield stress k i . The process of reducing the set of equations will be shown next. In a rst step the plastic strains (36) are inserted in the elasticity relation (35), Si = T S − i Ĉ N i (43) Ĉ = P 3 C P 1 (44) Xi = i {X(n) + i cN i } (45) with Then we recast equation (37) to obtain with i = (i ; pi ) = 1+ q 1 2 i 3 b=(1 + api ) Equations (43) and (45) are inserted in the dierence tensor Si − Xi r 2 i i kN Si − Xi = T S − i X(n) − i [Ĉ + i c1]N i = 3 where we have exploited the yield condition (41) with denitions (29)2; 3 r q 2 i i i i i i i i i k S − X = n = N niT P 2 ni = N kS − X kPSC = N 3 (46) (47) (48) Now, we obtain from equation (47) the normal tensor Ni = i −1 T { S − i X(n) } with r i ? 1998 John Wiley & Sons, Ltd. i = Ĉ + i i c+ 2 i k 3 (49) ! 1 (50) Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1485 NUMERICAL ASPECTS OF A CYCLIC PLASTICITY MODEL which is introduced in relation (45), Xi = i {X(n) + i c i −1 T { S − i X(n) }} (51) If we now calculate the norm of the tensor (1=i )Xi − X(n) = i cN i , keeping in mind that kN i kPSC = 1, the rst scalar equation arises: 1 i X − X(n) = i c (52) i PSC with 1 i X − X(n) = i c i i −1 T { S − i X(n) } (53) or g (i ; pi ; k i ) = 0 (54) i −1 T (55) with g (i ; pi ; k i ) := i c − ki c { S − i X(n) }kPSC The aforementioned second and third equations result from rearranging of equations (38) and (40) gp (i ; pi ; k i ) = 0 i i (56) i gk ( ; p ; k ) = 0 with i i i i gp ( ; p ; k ) := p − p(n) + q (57) i 2 i 3 ( =sp )kX kPSC 1+ q (58) 2 i 3 =sp q k (n) + 23 i k0 q gk (i ; pi ; k i ) := k i − 1 + 23 i =(1 + pi i ) (59) where the unknown i is calculated from equation (39), q ˆ (n) + 23 i =s (1 − cos2 i )=2 q i = 1 + 23 i =s (60) and is introduced in equation (57). Here, we use cos i from equation (42) and kXi kPSC from equation (51). Remark. The inversion of the matrix i ? 1998 John Wiley & Sons, Ltd. i = is simple to perform since 0 11 12 0 12 22 0 0 33 i has the structure Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1486 S. HARTMANN, M. KAMLAH AND A. KOCH with the inverse i −1 = 11 1 22 − 2 12 − 22 0 12 − 0 11 0 0 12 2 11 22 − 12 : 33 Remark. The three non-linear equations (54), (56) and (57) with the unknowns i , pi and k i have to be solved. We choose a Newton-like single-step algorithm (see Schwarz27 ). For a given non-linear system of equations, fj (x1 ; : : : ; xn ) = 0; i = 1; : : : ; n, this scheme is given by xj(l+1) = xj(l) − (l+1) (l) fj (x1(l+1) ; : : : ; xj−1 ; xj ; : : : ; xn(l) ) (l+1) (l+1) (l) @ ; : : : ; xj−1 ; xj ; : : : ; xn(l) )|xj = x(l) @xj fj (x1 j ; j = 1; : : : ; n l symbolizes the iteration index. This algorithm has the advantage of needing only the derivatives @g =@i , @gp =@pi and @gk =@k i and the disadvantage of only linear convergence. Remark. The algorithm applied to the three-dimensional case is quite similar. Merely equation (43) and therefore equation (47) is much more simple, since the elasticity matrix Ĉ reduces to the scalar 2G, i.e. the inversion of the matrix i in equation (49) is trivial. Once i , pi and k i have been calculated, the normal tensor N i is determined by equation (49). Then the stresses, using equations (35) and (43), Si = P−1 {T S − i Ĉ N i } Ti = P−1 3 3 i i = C{Ei − E(n) p } − C P1 N (61) and the internal variables Xi from equation (45) as well as i from equation (60) have to be evaluated. Table I summarizes the stress algorithm. CONSISTENT TANGENT OPERATOR The Newton–Raphson method carries out the equilibrium iteration and requires a consistently derived functional matrix K iT , which is assembled by using the so-called consistent tangent operator (19). This means that the equations determining the stresses, i.e. equation (61), have to be dierentiated with respect to the current strains Ei : dTi d i i {C{Ei − E(n) p } − C P 1N } i = dE dEi " " T i ## di i i dN − = C I − P1 N dEi dEi CiL = (62) (63) In this expression we have to calculate the two derivatives {di =dEi } ∈ R3 and [dN i =dEi ] ∈ R3×3 . {di =dEi } is derived from the implicit function theorem applied to the set of equations (54), (56) ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1487 NUMERICAL ASPECTS OF A CYCLIC PLASTICITY MODEL Table I. Stress algorithm (1) Calculate elastic predictor T S = P 3 C{Ei − E(n) p } (2) Proof yield function 2 IF f(T S; X(n) ; k (n) ) = 12 {T S − X(n) }T P 2 {T S − X(n) } − 13 k (n) ≤ 0 THEN Ti = P−1 T S, Xi = X(n) , pi = p(n) , i = (n) , k i = k (n) 3 (3) ELSE, calculate i , pi and k i from q −1 i i i T i (n) 2 i { S − X } =0 g ( ; p ; k ) = c − Ĉ + c + 3 k 1 PSC p2 i (n) i i i i i gp (i ; pi ; k i ) = pi − i i i i gk ( ; p ; k ) = k − + p =sp kX kPSC 1+ k (n) + 1+ 3 p 2 p2 p2 3 3 =0 i =sp i k0 =0 i 3 1+pi i with i (equation (46)), i (equation (60)), cos i (equation (42)) and Xi from equation (51) (4) Renew N i (equation (49)) and calculate Ti = P−1 Si = P−1 {T S − i Ĉ N i } 3 3 Xi = i {X(n) + i cN i } i = (n) + p2 3 ˆ i =s (1 − cos2 i )=2 1+ p2 3 i =s and (57), g (i (Ei ); pi (Ei ); k i (Ei ); Ei ) = 0 gp (i (Ei ); pi (Ei ); k i (Ei ); Ei ) = 0 (64) gk (i (Ei ); pi (Ei ); k i (Ei ); Ei ) = 0 since the values i , pi and k i depend implicitly on the current strains Ei . Therefore, the chain rule can be applied ( j = ; p; k): dgj dEi T = @gj @i ? 1998 John Wiley & Sons, Ltd. di dEi T + @gj @pi dpi dEi T + @gj @k i dk i dEi T + @gj @Ei T = 0T (65) Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1488 S. HARTMANN, M. KAMLAH AND A. KOCH This is equivalent to the linear system of equations with several right-hand sides T i T @g @i @gp @i @gk @i @g @pi @gp @pi @gk @pi @g @k i @gp @k i @gk @k i @g d @Ei dEi T T @gp i dp =− i @Ei ; dE i T @g T dk (66) k dEi @Ei as proposed by Tamme22 or Hartmann et al.,10 from which the unknown vectors {di =dEi }, {dpi =dEi } and {dk i =dEi } can be calculated. There are several possible solution strategies: (1) Analytical derivation of the coecient matrix as well as the right-hand side vector and the analytical calculation of the unknowns (This will be the most time consuming option for the analyst) (2) Analytical derivation of the coecient matrix as well as the right-hand side vector and the numerical calculation of the unknowns by solving the set of linear equations (This semi-analytical method has been applied in our analyses) (3) Numerical derivation of the coecient matrix as well as the right-hand side vector performed by numerical dierentiation and the numerical calculation of the unknowns by solving the set of linear equations (this purely numerical procedure turned out to consume approximately 50 per cent more computation time than the 2nd alternative during the subsequently investigated uniaxial specimens). The further primary unknown matrix [dN i =dEi ] can be derived by a straightforward dierentiation of equation (49). The dierentiation with respect to the current strains Ei is given by q i T i T i 2 i d dN 3b i −1 (n) I N i d Ĉ + i c 1 − i P = C − X − 3 1 + api dEi dEi dEi ( i i 2 i T r i T )# 2 dk dp i + N abc + (67) 1 + api 3 dEi dEi with di dEi r = 2 2 i b 3 1 + api ai 1 + api dpi dEi − di dEi (68) The calculation of matrix (67) is derived in Appendix I. Remark. The calculation of the coecients @gj =@, j = ; p; k, = i ; pi ; k i is straightforward. There may occur in radial processes singular expressions, since 1 − cos2 i occurs in the denominator. This problem is avoided by numerical perturbation. NUMERICAL EXAMPLES The aforementioned applied nite element formulation based on the principle of virtual displacements and the derived stress algorithm has been implemented in the nite element program FEAP.28 ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) NUMERICAL ASPECTS OF A CYCLIC PLASTICITY MODEL 1489 Figure 1. Tension-compression test specimen (geometric data), boundary conditions of a quarter (coordinates point I (0·158,0·247), point II (5·842,10·747)), discretization with 360 four-noded isoparametric elements (414 nodes, 764 unknown DOF) The calculations are performed on an IBM RISC 6000 (380) computer with a UNIX-type operating system. We investigate the inuence of non-proportionality occuring in uniaxial loaded specimens. To this end, we choose two examples: a tension-compression test specimen tapered in the middle region and a plane membrane with a circular hole. Tension-compression test specimen In the rst example, we simulate a plane tension-compression test with a specimen tapered to the middle (see Figure 1(a)). Experimental investigations of material behaviour by such specimens relies often on the assumption of a homogeneous middle region. However, Luhrs5 as well as Luhrs and Haupt2 showed that for an axial loaded axisymmetric specimen the assumed homogeneity may get lost during cyclic processes. This was shown experimentally and by numerical simulations using the material model (1)–(11) enhanced by viscous modelling with a Perzyna–Chaboche-type viscoplasticity model (see Chaboche29 ). The inhomogenization of the thinner region is caused by the transition from the thicker to the thinner part of the sample and evolves during cyclic plastic processes. The change of strain distribution depends in the experiments essentially on the investigated material and in the simulations on the constitutive model. In Figure 1(c), the ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1490 S. HARTMANN, M. KAMLAH AND A. KOCH Table II. Material parameters E MPa 21 0000 k0 c b a sp s ˆ = 0·3 MPa 180 MPa 50 000 = 830 MPa−1 0·005 = 0·01 MPa−1 0·01 = 50 = 0·01 = 1 Figure 2. Axial-strain distribution after the 1st, 7th, 13th, and 20th cycle without non-proportional isotropic hardening ( = 0) discretization with four-noded isoparametric elements is shown and in Figure 1(b) the boundary conditions concerned are depicted. The rst calculation was performed with the material parameters shown in Table II, where we omit the non-proportional inuences ( = 0, see equations (12) – (14)). The material parameters are taken from Reference 1 and refer to stainless steel XCrNi18.9. The parameters a, sp , and s are slightly modied to reinforce some cyclic phenomena occurring in the subsequent calculations. The displacement driven load (maximal displacement amplitude d = 0·3 mm) was applied over 20 cycles with 200 load steps per cycle. Figure 2 shows that no essential inhomogeneity occurs in the slender part of the sample and after only a few cycles the strain distribution is stationary. This is conrmed by Figure 3 in which the axial-strains versus the axial stresses of points I and II ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) NUMERICAL ASPECTS OF A CYCLIC PLASTICITY MODEL 1491 Figure 3. Axial-strain versus axial-stresses at point I and point II (‘proportional’ model ( = 0)) (their location is shown in Figure 1) are depicted. The magnitude of the axial-strains and axialstresses at points I and II are nearly identical and the saturation level is reached very quickly. This result is observed in the plane strain case too (private communications with G. Luhrs), i.e. that the inhomogenization eects shown in References 5 and 2 are initiated from the state variables characterizing the circumferential direction. However, if we take the non-proportional hardening model into account (, and s dened in Table II) the simulation answers are quite dierent. Now, the axial strains are inhomogeneously distributed in the slender part of the sample and it seems that no saturation of the strain distribution takes place until the 20th cycle (see Figure 4). Furthermore, it seems that the inhomogenization in the slender region continues (see Figure 5), which can be observed at the hysteresis loops of point I. This process is not monotonic, i.e. the mean strains increase and decrease. A measure of non-proportionality is characterized by the isotropic hardening function (25). In Figure 6 the values of the ‘yield stresses’ k are depicted for the 2nd, 3rd, 4th and 20th cycle, thereby supporting the statement of inhomogenization. The resulting deformations of the specimens shown in Figure 7 are qualitatively dierent, but quantitatively of similar order. Strip with a circular hole The next example treats the well-known benchmark example of the plane membrane with a circular hole. The discretization is shown in Figure 8. The load was applied horizontally in 20 cycles with a maximal amplitude of d = 0·3 mm with 150 load increments per cycle. Of course, the sample strain distribution is in general inhomogeneous. In Table II the corresponding material parameters are listed. If we choose in a rst calculation = 0, i.e. no non-proportional hardening occurs, the equivalent strains of cycle 1, 7, 13 and 20 are depicted in Figure 9. These strains are greater than those of the non-proportional model shown in Figure 10, as is to be expected. Although the non-proportional ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1492 S. HARTMANN, M. KAMLAH AND A. KOCH Figure 4. Axial-strain distribution after the 1st, 7th, 13th, and 20th cycle (‘non-proportional’ model) Figure 5. Axial-strain versus axial-stresses at point I and point II (‘non-proportional’ model) ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) NUMERICAL ASPECTS OF A CYCLIC PLASTICITY MODEL 1493 Figure 6. Isotropic hardening distribution in the 2nd, 3rd, 4th, and 20th cycle eects appear to be nearly saturated, i.e. the evolution of additional changes in the strain distribution seems to be nished, the isotropic hardening distribution shows additional small changes between the 7th and the 20th cycle, i.e. the whole process apparently continues (Figure 11). If we look at the displacements of the boundary of the hole in radial direction (pressure mode in the 20th cycle), we see that there are quite dierent deformations (see Figure 12). The horizontal displacements (’ = 0) dier about approximately 0·01 mm. The vertical displacements are nearly zero. Obviously, the non-proportional hardening model shows a stier response. CONCLUSIONS In this paper we apply the return mapping algorithm to a non-proportional plasticity model proposed by Haupt and Kamlah.1 The return mapping algorithm leads to the solution of three non-linear equations. Under plane stress conditions it can be shown that for a tension-compression specimen tapered towards the middle, the inhomogenization process dealt with by Luhrs and Haupt2 can only be simulated by consideration of evolution equations representing non-proportional hardening eects. However, in other inhomogeneous, deformed structures, taking the non-proportional hardening model into account, only very small dierences are apparent, compared to the propor? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1494 S. HARTMANN, M. KAMLAH AND A. KOCH Figure 7. Deformed meshes after the 20th cycle of the proportional and the non-proportional model (magnied 50 times) Figure 8. Mesh of the plane strip with a circular hole. Boundary conditions of a quarter, discretization with 200 four-noded isoparametric elements (231 nodes, 429 unknown DOF) ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) NUMERICAL ASPECTS OF A CYCLIC PLASTICITY MODEL 1495 Figure 9. Equivalent strains after the 1st, 7th, 15th and 20th cycle with the proportional model ( = 0) Figure 10. Equivalent strains after the 1st, 7th, 15th and 20th cycle with the non-proportional model tional model. Consequently, this study contributes to a more general experience: the response of an initially homogeneous structure is highly sensitive to details of the constitutive model, while the behaviour of an initially inhomogeneous structure is dominated by some basic features of the constitutive model (like the existence of a yield stress). APPENDIX I To form the derivative [dN i =dEi ] we apply the Gateaux-derivative to the normal tensor (49) i dN DEi N i (Ei )[H] = H dEi ( i T ) d i −1 i i −1 (n) T i (n) = DEi (E )[H] { S − X } + (69) P 3C − X H dEi ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1496 S. HARTMANN, M. KAMLAH AND A. KOCH Figure 11. Isotropic hardening distribution of a part of the strip (50 × 50 mm) in the 2nd, 3rd, 7th, and 20th cycle Figure 12. Radial displacements at the boundary of the hole after the 20th cycle of the proportional and the non-proportional model Here, we have dierentiated the trial stresses T S from equation (33) DEi T S(Ei )[H] = P 3 C H (70) and look for the Gateaux-derivatives of i and the inverse of i dened by equations (46) and (50), respectively. q The dierential of i = (i ; pi ) = (1 + 23 bi =(1 + api ))−1 is achieved by applying the chain rule: with r r 2 2 @(i ; pi ) 2 bi 2 i @(i ; pi ) i ab =− ; = (71) @i 3 1 + api @pi 3 1 + api ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) NUMERICAL ASPECTS OF A CYCLIC PLASTICITY MODEL we get ( T T T ) @(i ; pi ) dpi di @(i ; pi ) di H= + H DEi ((E ); p(E ))[H] = @i @pi dEi dEi dEi (r i i )T 2 2 i b ai dp d = H − i i i 3 1 + ap 1 + ap dE dEi i 1497 i The dierentiation of the inverse of DEi i −1 i (72) , (Ei )[H] = − i −1 DEi i (Ei )[H] i −1 leads to the calculation of ( ( i T r i T )) T T 2 dk d di di i i i +I +c + H DEi (E )[H] = Ĉ 3 dEi dEi dEi dEi (73) (74) If we insert this equation in the expression, DE i i −1 (Ei )[H] {T S − i X(n) } = − i −1 DEi i (Ei )[H] N i where we have used equation (49), after a few calculations we eventually arrive at "" ! # r T 2 bi i di i −1 i i −1 i T i (n) i (E )[H] { S − X } = − Ĉ + c 1 − I N DEi 3 1 + api dEi ( r 2 i T r i T )# i i 2 2 dk dp H + N i abc + i 3 1 + api 3 dEi dE (75) ACKNOWLEDGEMENTS We would like to express our deep gratitude to Professor Dr. Peter Haupt for his generous support and for teaching us his approach to the fascinating eld of mechanics. REFERENCES 1. P. Haupt and M. Kamlah, ‘Representation of cyclic hardening and softening properties using continuous variables’, Int. J. Plasticity, 11, 267–291 (1995). 2. G. Luhrs and P. Haupt, ‘On the evolution of inhomogeneous deformations in test specimens under cyclic loading conditions’, in: D.R.J. Owen, E. Onate and E. Hinton, eds., ‘Computational plasticity. Fundamentals and applications’, COMPLAS V Int. Conf. on Computational Plasticity, Barcelona, 1997, pp. 945–950. 3. P. J. Armstrong and C. O. Frederick, ‘A mathematical representation of the multiaxial Bauschinger eect’, General Electricity Generating Board, Report No.: RD=B=N731, Berceley Nuclear Laboratories 1966. 4. P. Haupt, M. Kamlah and C. Tsakmakis, ‘Continuous representation of hardening properties in cyclic plasticity’, Int. J. Plasticity, 8, 803 – 817 (1992). 5. G. Luhrs, Randwertaufgaben der Viskoplastizitat. Modellierung, Simulation und Vergleich mit experimentellen Daten aus zyklischen Prozessen und Umformvorgangen, Doctoral Thesis (Report No. 1=1997) at the Institute of Mechanics, University of Kassel, Germany, 1997. 6. H. S. Lamba and O. M. Sidebottom, ‘Cyclic plasticity for nonproportional paths: part 1 — cyclic hardening, erasure of memory, and subsequent strain hardening experiments’, J. Engng. Mater. Technol., 100, 96 –103 (1978). ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998) 1498 S. HARTMANN, M. KAMLAH AND A. KOCH 7. E. Tanaka, S. Murakami and M. Ooka, ‘Eects of plastic strain amplitudes on non-proportional cyclic plasticity’, Acta Mech., 57, 167–182 (1985). 8. A. Benallal and D. Marquis, ‘Constitutive equations for nonproportional cyclic elasto-viscoplasticity’, J. Engng. Mater. Technol., 109, 326 –336 (1987). 9. A. Benallal, P. LeGallo and D. Marquis, ‘The role of mean strain on the stress response in nonproportional cyclic plasticity’, in: A.S. Khan and M. Tokuda eds., Advances in Plasticity 1989, Proc. Plasticity ’89, Oxford, 1989, pp. 203 –206. 10. S. Hartmann, G. Luhrs and P. Haupt, ‘An ecient stress algorithm with applications in viscoplasticity and plasticity’, Int. J. Numer. Methods Engng., 40, 991–1013 (1997). 11. D. L. McDowell, ‘A Two surface model for transient nonproportional cyclic plasticity, Part 1 developement of appropriate equations’, J. Appl. Mech., 52, 298 –308 (1985). 12. A. Benallal and D. Marquis, ‘Eects of non-proportional loadings in cyclic elasto-viscoplasticity: experimental, theoretical and numerical aspects’, Engng. Comput., 5, 241–247 (1988). 13. A. Benallal, P. LeGallo and D. Marquis, ‘An experimental investigation of cyclic hardening of 316 stainless steel and of 2024 aluminium alloy under multiaxial loadings’, Nucl. Engng. Des., 114, 345 –353 (1989). 14. J. C. Moosbrugger and D. L. McDowell, ‘On a class of kinematic hardening rules for nonproportional cyclic plasticity’, J. Engng. Mater. Technol., 111, 87–98 (1989). 15. O. T. Bruhns, T. Lehmann and A. Pape, ‘On the description of transient cyclic hardening behaviour of mild steel CK 15, Int. J. Plasticity, 8, 331–359 (1992). 16. N. Ohno, ‘Recent topics in constitutive modelling of cyclic plasticity and viscoplasticity’, Appl. Mech. Rev., 43 (11), 283 –295 (1990). 17. M. Kamlah, Zur Modellierung des Verfestigungsverhaltens von Materialien mit statischer Hysterese im Rahmen der phanomenologischen Thermomechanik, Doctoral Thesis (Report No. 3=1994) at the Institute of Mechanics of the University of Kassel, Germany, 1994. 18. T. J. R. Hughes, The Finite Element Method, Prentice-Hall, Englewood Clis, NJ, 1987. 19. J. C. Simo and M. S. Rifai, ‘A class of mixed assumed strain methods and the method of incompatible modes’, Int. J. Numer. Methods Engng., 29, 1595 –1638 (1990). 20. J. C. Simo and R. L. Taylor, ‘Consistent tangent operators for rate-independent elastoplasticity’, Comput. Methods Appl. Mech. Engng., 48, 101–118 (1985). 21. H. J. Braudel, M. Abouaf and J. L. Chenot, ‘An implicit and incremental formulation for the solution of elastoplastic problems by the nite element method’, Comput. Struct., 22, 801– 814 (1986). 22. A. Tamme, ‘A unied approach to modelling and numerical solution of coupled eld problems in nonlinear solid mechanics’, in: D. Besdo and E. Stein, eds., Finite Inelastic Deformations — Theory and Applications, IUTAM Symposium Hannover=Germany 1991, Springer, Berlin, 1992, pp. 93 –106. 23. J. C. Simo and R. L. Taylor, A return mapping algorithm for plane stress elastoplasticity’, Int. J. Numer. Methods Eng., 22, 649 – 670 (1986). 24. J. C. Simo and S. Govindjee, ‘Exact closed-form solution of the return mapping algorithm in plane stress elastoviscoplasticity’, Engng. Comput., 5, 254 –258 (1988). 25. S. Hartmann and P. Haupt, ‘Stress computation and consistent tangent operator using non-linear kinematic hardening models’, Int. J. Numer. Methods Engng., 36, 3801–3814 (1993). 26. F. Auricchio and R. L. Taylor, ‘Two material models for cyclic plasticity: non-linear kinematic hardening and generalized plasticity’, Int. J. Plasticity, 11, 65 –98 (1995). 27. H. R. Schwarz, Numerische Mathematik, B.G. Teubner Verlag Stuttgart, Germany, 1986. 28. O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, Vols 1 and 2, 4th edn., McGraw-Hill, London, 1989 and 1991. 29. J. L. Chaboche, ‘Time independent constitutive theories for cyclic plasticity’, Int. J. Plasticity, 2, 149 –188 (1986). ? 1998 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 42, 1477–1498 (1998)

1/--страниц