close

Вход

Забыли?

вход по аккаунту

?

441

код для вставкиСкачать
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)
Документ
Категория
Без категории
Просмотров
4
Размер файла
289 Кб
Теги
441
1/--страниц
Пожаловаться на содержимое документа