- Email: [email protected]
Eulerian constitutive model for finite deformation plasticity with anisotropic hardening
Mechanics of Materials 7 (1989) 279-293 North-Holland
279
EULERIAN CONSTITUTIVE MODEL FOR FINITE DEFORMATION PLASTICITY WITH ANISOTROPIC HARDENING George Z. VOYIADJIS and Peter I. KATTAN Department of Civil Engineering, Louisiana State University, Baton Rouge, 1.,,4 70803, U.S.A. Received 21 July 1988; revised version received 31 August 1988
A constitutive model is presented for finite strain plasticity. The model incorporates both isotropic and kinematic hardening of the Ziegler type. The corotational rate used here is in fine with the theory suggested by Paulun and Pecherski (1985) but not necessarily confined to the yon Mises type yield criterion and the Prager hardening rule. The aspect of integration of the corotational rates is also discussed here. The use of the integration of the material rate of tensors with time as a substitute for the proper integration with time of corotational rates leads to mathematical inconsistencies of the theory of Lie derivatives. The problem of simple shear is investigated and compared with other works.
Introduction The problem of oscillatory stresses produced in simple shear with the use of the Zaremba-Jaumann derivative in conjunction with large deformation kinematic hardening plasticity was pointed out by Nagtegaal and De Jong (1982). In an earlier paper, TruesdeU (1956) noticed the oscillation for simple shear and gave a closed form solution. Later, Dienes (1979) observed the same phenomenon using a hypo-elastic material model. The subject of finite deformation plasticity and in particular the objective of defining an appropriate corotational stress rate and an appropriate constitutive formulation for finite deformation has received wide attention since then. Lee et al. (1983) proposed a modified corotational rate using the spin of the principal direction of a with the largest absolute eigenvalue, a is the deviatoric component of the shift-stress tensor. An alternate approach by Onat (1982, 1984) defines the spatial spin equal to the antisymmetric part of d~'j:ajk multiplied by a constant, d " is the plastic component of the spatial strain rate. The nonoscillatory solution for simple shear is obtained by the proper choice of the constant. Dafalias (1983, 1985) and Loret (1983) obtained similar relations by associating the corotational rate with
the material substructure as defined by Mandel (1973). Mandel (1973) used the triad of director vectors attached to the material substructure and developed the theory of plasticity such that the substructure corotational rate is defined in terms of the spin of the director vectors. He postulated that the constitutive relations require not only the plastic component of the spatial strain rate but also the plastic component of the spatial spin tensor. Dafalias (1985) and Loret (1983) discussed the macroscopic constitutive relations for the plastic spin using the representation theorem for isotropic second-rank antisymmetric tensor valued functions. The importance of the material substructure in defining objective corotational rates is also argued by Pecherski (1985). In inelastic finite deformations of polyerystalline metals, the material moves with respect to the underlying crystal lattice. The lattice itself undergoes elastic deformations and relative rigid-body rotations due to the lattice misorientation (Pecherski, 1985). The work outlined above by Lee et ai. (1983), Onat (1982, 1984), Dafalias (1983, 1985), and Loret (1983), imposes a retardation of the material spin W in order to obtain a non-oscillatory solution for the simple shear problem. The analysis of the solution of the simple shear problem presented
0167-6636/89/$3.50 © 1989, Elsevier Science Publishers B.V. (North-Holland)
280
G.Z. Voyiadis, P.L Kattan / Finite deformation plasticity
by Onat (1982, 1984), Dafalias (1983, 1985), and Loret (1983), results in an unbounded non-oscillatory solution for the shear stress that increases monotonically with increased deformation. Concurrently, the normal stress approaches an asymptotic upper bound. In the case of Lee et al. (1983), both the shear and the normal stresses are unbounded and increase monotonically with increased deformation. We also note that in the case of Dafalias (1983, 1985), Onat (1982, 1984), and Loret (1983), the principal directionsto~ a tend toward the bisector direction of the plane coordinate axes while in the case of Lee et al. (1983), the maximum principal direction of a inclines towards the horizontal axis. The above proposed solutions fail in the proper prediction of the shear stress-shear strain characteristic and the Swift effect in torsion of thin-walled tubes (Paulun and Pecherski, 1985). Other authors have followed different approaches for the proper choice of the corotational objective stress rate. Atluri (1984) based on the idea of a complete hypo-elastic law, modifies the rate of the back stress equation for the case of a rigid-kinematic hardening plastic model. Fressengeas and Molinari (1984), Johnson and Bammann (1984), Simo and Pister (1984), Voyiadjis (1984), Voyiadjis and Kiousis (1987), Moss (1988), and Voyiadjis (1988) have also discussed different aspects of proper choice of the objective stress rate in finite deformation analysis. An ASME publication by William (1984) summarizes the debate on this subject. In this work a constitutive model is presented for finite strain plasticity. The model incorporates both isotropic and kinematic hardening. The modified Prager hardening rule proposed by Ziegler (1959), is used. The corotational rate used here is in line with the theory suggested by Paulun and Pecherski (1985). This theory is not only based on the question of how to avoid the unwanted oscillatory stresses in the simple shear problem, but also refers to the proper constitutive description of anisotropic hardening in the course of finite plastic deformation process. This theory is based on Mandel's (1973, 1981, 1982) theory of material substructure, and the plastic spin expressed by Dafalias (1985) and Loret (1983).
The generalized theory presented here leads to the satisfactory solution of the simple shear traction problem. This is obtained using a numerical solution to solve the governing differential equations. A closed form solution is also obtained for a special case of the differential equations. An alternate numerical integration scheme proposed by Goddard and Miller (1966) is also used here for the solution of the simple shear problem. This scheme is based on integration of the constitutive relations using an appropriate corotational rate. The use of the integration of the material rate of tensors with time as a substitute of the proper integration with time of corotational rates leads to mathematical inconsistencies of the theory of Lie derivatives (Zhong-Heng, 1963, 19841 Yano, 1957). For an isotropic linear hypo-elastic material, for deformations in which the directions of principal stretches remain fixed in the body, the procedure outlined by Ralph and Bathe (Williams. 1984) is equivalent to that of Goddard and Miller (1966). They have presented an equation for the inverse of the Jaumann derivative.
Theoretical formulation
The yield condition used in this work combines both isotropic hardening and kinematic hardening of the type proposed by Ziegler (1959)
(1) A particular form of this yield condition is the generalized von Mises type, with both kinematic and isotropic hardening such that /-
~ ( ~ , - a ~ , ) ( ~ , - ~ , ) - o: + c~ = 0
(2)
In the above expressions rkt is the deviatoric component of the Cauchy stress tensor Okt, a n d akt is the deviatoric component of the shift stress tensor. The constant o corresponds to the initial tensile or compressive yield stress from uniaxial loading, and c is a constant that describes the isotropic component of hardening. In eqns. (1) and (2) i¢, the isotropic hardening parameter, is obtained from the following relation: pp
~¢ = p k t d k t
(3)
G,Z. Voyiadis, P.L Kattan / Finite deformation plasticity
where d~; is the plastic component of the spatial strain rate tensor, and Pkt is a function of the Cauchy stress tensor and the accumulated elastic strain. The superdot in eqn. (3) denotes material time differentiation. Assuming small elastic strains but finite plastic deformations, the spatial rate is decomposed into (Nemat-Nasser, 1979, 1983; Lee, 1981) t
Pl
dkt = dkt + dkt
(4)
where d~t is the elastic component of the spatial strain rate. The above equation will be correct for any amount of elastic strain if the physics of elasto-plasticity is invoked, for example the case of single crystals. A thorough account of this is given by Asaro (1983). An associated flow rule of the type shown below is postulated a;; =
A aokl a/
(5)
in finite plastic deformation analysis. We note in Figure 3a of Paulun and Pecherski (1985) that the shear stress obtained from Paulun and Pecherski (1985), Onat (1982, 1984), Dafalias (1983, 1985) and Loret (1983) for the simple shear problem increases monotonically with increased deformation. On the contrary the shear stress obtained from Lee et al. (1983) approaches an asymptotic upper bound with increased deformation. We also note that the upper bounds for the normal component all of the back stress obtained by Loret (1983), Dafalias (1983, 1985) and Onat (1982, 1984) for the non-oscillatory solution of the problem of simple shear do not seem to be justified physically (Paulun and Pecherski, 1985). The generalization of the modified spin tensor 12 to the three-dimensional case is accomplished through the use of the theory of material substructure as defined by Mandel (1973). Defining 12" such that 12* = w -
where A is a positive scalar function. Following the Prager-Ziegler kinematic hardening rule, the evolution equation for a is given by &kt =/i
(~'k,- ak,)
(6)
where &kt is an objective rate of the back stress that is corotational with the substructure (Paulun and Pecherski, 1985). The objective rate for a is given by (Paulun and Pecherski, 1985). & = a, - 12a + a12
(7)
where 12 is a modified spin tensor. In simple shear 12 causes a retardation to the material spin W which consequently provides a non-oscillatory solution. This is accomplished through the use of an influence function ~o decreasing with the shear strain. The determination of the function ~o is such that the angular velocity ? corresponding to the spin W in simple shear is reduced to the angular velocity 6 of a single material line element. The objective rate of a given by eqn. (7) not only avoids the unwanted oscillatory stress in the simple shear problem, but also provides proper constitutive description of anisotropic hardening
281
(8)
12
we express eqn. (5) as follows:
a = 7, + 12"a - a12"
(9)
where g is the Zaremba-Jaumann derivative given by (10)
7, ffi dt - W a + a W
Using the general form of the yield function (1), the plastic spin may be expressed by the following equation proposed by Dafalias (1985) w " = ~ , ( ~ o - o~) + ~ ( ~ o
- o~)
+ ~ ( ~ 0 ~ - o~o) + ~ , ( ~ o ~ - ~ 0 ~ ) + "% ( o a o 2 - o2ao)
(11)
Using the general case of the yield function (1), we may express 8 f / a o U such that at" _ a/ 0.11 of 0.12 a/ aJ3 aaU aJI aoU + aJ2 a%j + aJ3 aoU
(12)
or
af ; 0%
af af aJ, au + 2 ~ %
af + 3~-~3o..,,.,
(13)
282
G.Z. Voyiadis, P.I. Kattan / Finite deformation plasticity
where J~ = %,, .I2 --- o,j~i, and -/3 = % ~ k O , V Making use of eqns. (13) and (5) we obtain:
The function ~7 is equation (19) for the case of simple shear yields (Paulun and Pecherski, 1985):
( a d " - d"ot ) rl = --~--(1 - ~o)5'
~f oo~)] = A[2~-~f4(ao-o. ) +3-~3(aoo-
or
(14) While retaining the first term only of eqn. (11), it is evident from eqn. (14) that in the absence of J3 terms, the plastic spin may be expressed as follows (Dafalias, 1983, 1985): W " --- f l ( a d " - d " a )
(15)
where 7ll = 2/Z'i ~-~2
(16)
However, eqn. (14) indicates that a particular form of expression (11) such that: W " = ~ , ( a o - oa) + ~13(ao 2 - o2a)
(17)
may be reduced to the form of eqn. (15) provided that the yield function is of such form that constraint (14) is satisfied. Constraint (14) may also be indicative that yield functions containing "/3 terms require that the plastic spin be expressed by eqn. (17) and not by setting ~!3 = 0. That is the 7/3 term needs to be present in the plastic spin expression. In the following work, it is assumed that the yield function does not contain J3 terms but is nevertheless not limited to the von Mises functional form. Equation (15) may be expressed in indicial notation as follows: tP
PP
tt
Wit = fl( a,,,,d,,,t- ak,na.,t)
(18)
where 13 is a scalar function of the isotropic invariants of a and o. The spin fl* is determined from the plastic spin I4,'" and is expressed as follows (Paulun and Pecherski, 1985): It* = n~'
~,2 2 1 + ~,2 ~'
n
y = ~/3-e"
(23)
where e'"
V~"JUklUkl
(24)
We may generalize to three-dimensional analysis the expression for ~ to be (Paulun and Pecherski, 1985)
=} 1 + 3(~") 2
~1 = A~
(26)
such that ~ - - 2~
3(E")2 ~ 2 1 + 3(E") ~
(20)
0 f -0f ] 8°*t 8°*t
(27)
Similarly, using eqn. (5) in expression (20), we eliminate A from eqn. (20) such that: akin
( t r [ ( a k . , d " t - d ~ , ' a . t ) = ] ) '/2
(25)
J
where
(19)
a,7,..,
(22)
where y and 5' are the shear strain and shear strain rate, respectively. The above equation represents the difference between the constant angular speed produced by the material spin W and the angular velocity of the material line element lying initially along the vertical axis. This expression relates to the angular velocity corresponding to the spin 12" and shows to be a function of the accumulated plastic strain "~. Since y may be expressed in terms of the scalar ~", such that
where u,t =
(21)
/ [t
8/ ~Oml
ao.,,
8/ ~Okn ant
ao,.""'
283
G.Z. Voyiadis, P.L Kattan / F i n i t e deformation plasticity
Substituting the expressions for u and 7/ from eqns. (20) and (25) or (28) and (26), respectively, into eqn. (19), the spin D* is determined. Equation (19) may now be expressed as follows: a k l = a k l -- ~Pkmrtmt + 71aknVnl
(29)
We note that in the absence of any plastic deformations, we obtain v ~ Otkl . Olkl
(30)
Equations (15), (18), (20), (21), (22), (25), (27) and (28) have been obtained on the assumption that the yield function (1) has the specific form expressed by eqn. (2). Before deriving the elasto-plastic stiffness matrix, we need to determine the positive parameter /~ associated with the kinematic hardening rule. This parameter is determined assuming that the projection of ,~ on the stress gradient of the yield surface equals to bd":
Of
0om.
bake=&"" Of Of
(31)
Ookt
OOpq OOpq where b is a material parameter. Substituting ,t and d " in eqn. (31) from eqns. (6) and respectively and post-multiplying the resulting pression by Of/OOkl, we obtain the required pression for/2: of of 0a,., 0am.
o/ ( q'pq -- apq )
Ek,m, = X8kA,, + a(Sk..8,. + 8 k : , , )
for (5) exex-
(32)
(35)
and ~ and G are Lam~'s constants. We note that the term - , / ( v a - ap) in eqn. (34) acts a damage effect, during unloading and reloading in the elastic range, due to the accumulated plastic strains. The presence of this damage effect is clearly indicated in the experimental results obtained by Voyiadjis (1984, 1988). We note in this paper the degradation of the elastic stiffness, in uniaxial loading, due to accumulated plastic strains. Eliminating d ' from eqn. (33) through the use of expressions (4) and (5), we obtain:
The parameter A can be calculated from the consistency condition
f(%,, ak,, K)ffi 0
0/
t~ = A b
where the modulus of elasticity Ekt.,. takes the following form:
(37)
The interpretation of eqn. (37) can be made on a consistent mathematical basis with the aid of the theory of Lie derivatives (Zhong-Heng, 1963, 1984; Yano, 1957). Equation (37) is expressed as: +
°akt
~ K= 0
(38)
where the material time derivative of a scalar function is equal to its Lie derivative with respect to the velocity vector v. Making use of equations (3), (6), (4), (5), (32) and (36) in eqn. (38), we obtain an expression for A:
OOvq (39)
In order to obtain the elasto-plastic stiffness tensor, a linear elastic relation is assumed between the corotational stress rate d and the elastic component of the spatial strain rate tensor d ' such that
dkt = ektm.d=.
a/ a%b
~%d
a~/(%/•
Okl = E k l m . d m n -- ~lvkpapl + l]akqPql
(34)
a~
Q = a/ E=b~ a/
(33)
or v
where
a~/)b
a~
al
~au
~au
a/
(¢sh - %,) a%~
(4o)
G.Z. Voyiadis, P.I. Kattan /Finite deformation plasticity
284
Substituting for A, from eqn. (39) into eqn. (36), we obtain the elasto-plastic stiffness tensor which corresponds to the loading function given by expression (1) (provided it does not contain ,/3 terms). The resulting elasto-plastic stiffness tensor is given by Oklmn -~- Eklrn n
1 a/ af Q a~"'-qEpqmnEklijOo,)
(42)
The corotational rates expressed by eqns. (42) and (6) may be integrated numerically using a procedure proposed by Goddard and Miller (1966). Alternatively, by making use of expressions of the type given by relation (7), eqns. (42) and (6) can be expressed as follows: dkl = D k t , . . d m . - Okp~pt q- ~kqOql
7/ I i l
(41)
and ~kl ~ Dktmndmn
~2
:7
Fig. 1. Simple shear problem and material properties of metal. Material properties of aluminum alloy 2024 T4: v = 0.25, E =10,600 ksi, c=138 ksi, b = 73 ksi, o = 56 ksi, k =l and h = At ~ 0.0001.
The deformation gradient F and the velocity gradient L are given by, respectively,
(43)
e--
(44)
~=
1 0 0
kt
O
1 0
0 1
(47)
0
0
(48)
0
0
and akl=~(,Tki--~kl)--akp~plq-
ekqaql
Equation (43) and (44) may now be used in numerical algorithms (finite element methods) to solve boundary value problems.
L =d+ W
The problem of finite simple shear
where d is the rate of deformation (symmetric part of L ) and W is the total spin (antisymmetric part of L ) then we obtain
The problem of simple shear in the x x direction as shown in Fig. 1, is given by the following displacement fidd U! = k t x 2
(45a)
u2 = 0
(45b)
u~ = 0
(45c)
Since,
d=7
v I = kx 2
(46a)
v2 = 0
(46b)
03 = 0
(46c)
0
0
and W=
where k is a constant representing the sheafing strain rate and t is the time. The corresponding velocity field is expressed by
(49)
-1 0
0 0
(51)
Using eqns. (40), (41) and (35) together with eqn.
(50), we obtain dii = kEijl2 - 9--~-k[(,eq--%q)Epq,2E, j k , ( , k , - - a k , )
Q
l
(52)
G.Z. Voyiadis, P.L Kattan / Finite deformation plasticity
together with (54) and (55) and the following expressing for 12
Similarly for & we have 18kGb Q (~',2--f,2)('rkt--fk,)
akt=
(53) ~ = ~ o W = - ~ - [ -01 1 00 1 0 0 0 0
From eqns. (52) and (53), the corotational rates are given explicitly as follows: 36kG 2 = 2Ga- ~('q2 Q
d = 2Gd
and 18kGb
[
TII -- fill
T=
T12 -- f12
0
'/'12 -- f12
1"22 -- f 2 2
0
0
0
1"33 - - f33
36kG2 (~'12 )T+ -~-S ° ~ 0112
a = 18kGb ~ (, ~ 1 : - f 1 2 ) T + - ~ S ~
(55)
a= --ff-(,,2- f, )r
(60)
we obtain the following system of differential equations in 0 and a, respectively
(54)
- al,) r
where
285
(61) (62)
where
]
2012
(56)
S °=
022-011
O22--O11
--2012
0
0
!] (63)
and likewise for S ~.
and Q = 9(b + 2G)('rkt- fkt)(~'kt- ak,) -
Numerical solution of the differential equations of simple shear
(57)
+ 3COeq(,rpq fpq) -
Using the following two equations 0 = # - o~ + ~o
(58)
& = & - f12 + 12a
(59)
The problem of simple shear is first solved numerically by using a Runge-Kutta-Verner fifth and sixth order method for the solution of the
250
to = 0 2
b
~, , , ~ . I
0
I
2
3
4 Y
Fig. 2. Variation of all vs. 7 for different values of to.
0.006895. ~0~
I
I
I
I
5
6
7
8
9
G.Z. Voyiadis, P.I. Kattan /Finite deformation plastwity
286
o~=0.15 600 w=0.20 500
400
~:0.80
co = 0 4 0
300
= 1.00 ZOO
I00 .( I psi - 0 . 0 0 6 8 9 5
o
I
I
I
I
I
2
4
3
M Po)
I
i
I
1
l
5
6
7
8
9
7
Fig. 3. Variation of a12 vs. 7 for different values of ~0.
600
500
400
300
200
tO0
J
O
( Ipui = 0 . 0 0 6 8 9 5 MPo)
I
I
I
I
I
2
3
4
I
I
I
I
I
5 Y
6
7
8
9
Fig. 4. Variation of o12 vs. 7 for different ~. Legend: ~ Dafalias (1985), Loret (1983), Onat (1984); ( ~ Paulun and Pecherski (1985); (~) Dienes (1979)0 Dafalias (1983); and (~) Lee, Mallet and Wertheimer (1983).
G.Z. Voytadis. P.I. Kattan /Finite deformation plasticity
287
300
250
200 0
0
,5o
®
N
I00
50 ( I psi = 0 . 0 0 6 8 9 5 I I
I 2
Z 3
I 4
MPo )
I 5
I 6
I
I
I
7
8
9
7"
Fig. 5. Variation of a]2 vs. y for different ~0. Legend: see Fig. 4.
® 150 140 150
120 I10
I00 9o = ,ac _
80
®
70
60 5o 40 30 20 ~
0
I
2
3
.
.
w
4
5
•
I "Pall
6
Y
Fig. 6. o n versus "r for different expression of ,o. Legend: see Fig. 4.
288
G.Z. Voyiadis. P.I. Kattan / Finite deformation plasticity
governing differential eqns. (61) and (62). The value of ¢ in Figs. 2 and 3 for these equations is assumed to be constant. Several values for ~ are used ranging from 0.05 to 1 in the solution of the problem. In Fig. 2 the normal stress oll is plotted versus the shearing strain Y for different values of ~. It is clear from the figure that for values of greater than 0.2 oscillations in the stress-strain curve occur within this region (up to 1000% strain). We note that as the value of ~ approa¢hes zero, ./., the shearing stress increases mdefimtely as the sheafing strain increases. In Fig. 3, the shear stress ox2 is plotted versus the sheafing strain 7 for different values of ~. This figure also shows that as t~ increases, oscillations occur in the curve as explained earlier. It is interesting to note that for the case of t~ = 1, we obtain the solution for the Jaumann rate which is oscillating as shown in Fig. 3. We note in this •
figure that all values of ~ give similar results up to a strain of about 100%. In Figs. 4 - 6 the sheafing stress or2, the shift component a~2 and the normal stress o H are plotted respectively versus the strain V for different proposed expressions for ~. The numerical solution obtained for o12, o H, a H and a12 is shown as curve (4) in Figs. 7,8,9, and 10 respectively. It is noted that in all the curves in Figs. 4 and 5 oa2 and at2 increase monotonically without bound as 7 increases and no oscillations occur. Nevertheless, the expressions proposed by Dafalias (1985), Loret (1983), and Onat (1984) as well as Paulun and Pecherski (1985) seem to give better correspondence with the results shown in Figs. 7 and 10 for curves 3 and 4. On the other hand, Fig. 6 shows that for the stress o~, the expressions proposed by Dafalias (1983) and Dienes (1979) give results that are closely related to the results
600
500
,=~oo
~oo
®
200
I00 /"
0
( I psi = 0.006895 MPo) I
[
I
I
I
[
I
I
I
i
2
3
4
5
6
7
8
9
Fig. 7. Variation of Ola vs. Y for different solutions. Legend: ( ~ Closed form solution for the special case with w = 0.15; Numerical integration of constitutive equations with (~ - WQ; ( ~ Numerical integration of constitutive equations with Q = R; and ( ~ Numerical integration of differential equation with w = 0.15.
G.Z. Voyiadis, P.L Kattan /Finite deformation plasticity
289
ZSO
® 20C
( I p s i =0.
eD
-I< 15C b I00
®
50
,
0
Fig.
I
2
I
I
I
I
I
3
4
5
6
7
)'
8. Variation o f o i l vs. "f for different solutions. Legend: see
Fig.
J
8
9
7.
® J6oI 150 140 130 120
~
90
~
SO -
7O 6O 5O 40 3O 2O
IO
0
I
2
3
4
)'
5
6
Fig. 9. Variation o f a]; vs. ¥ for different solutions. Legend: see Fig. 7
7
8
9
290
G.Z. Voyiadis, P.L Kattan /Finite deformation plasticn)
25(
®
a~x
~
150
co I 0 0 =to
5O (
10si = 0 . 0 0 6 8 9 5
MPa)
® o
-50
I
2
3 " ~
~
7
8
9
~"
y
FIGURE I0 • VARIATION OF' a IZ vs. 7 FOR DIFFERENT SOLUTIONS Fig. 10. Variation o f a12 vs. 7 for different solutions. Legend: see Fig. 7.
obtained in Fig. 8 by numerically integrating the constitutive equations with the corotational rates given in eqns. (42) and (6). There is a special case of the differential eqns. (61) and (62) where an analytical solution can be obtained. This case is based on the assumption that ~'12= a~2. This assumption is one of two possible ways of achieving 033 = a33 = 0 in eqns. (61) and (62). Next, we investigate the nature of the solution for this case. The other possible alternative of obtaining o 3 = 0 is to assume T33 -- a33. However, this case cannot be solved analytically but can be solved numerically similarly to the solution presented in this section.
Special case of the analytical solution of the differential equations of simple shear For this case the system of eqns. (61) and (62) reduce to an uncoupled system of equations given
by
o 12
=
~
olfo,, /
-1
+
o
0
0"22]
-2
1/ o1: OJlo2z ! (64)
and 0
0t22 ]
0
0
1
at2 ~
- 2
0
0~22/
(65)
We note that in both systems of eqns. (64) and (65), the material parameters c and b do not appear explicitly. They can be incorporated in an appropriate expression for o~. But since ~ is assumed to be constant, this solution does not incorporate the material parameters b and c. F r o m eqns. (64) and (65), we obtain the following two uncoupled second order differential equa-
291
G.Z. Voyiadis, P.L Kattan / Finite deformation plasticity
tions: (66)
O12 q- ~2k2°12 --- 0
Chandra and Mukherjee (1986). This procedure is given by the integral:
a( t ) = Q( t )a( to)Qr( t )
and (67)
fi;12 -I- 602k2¢112 = 0
The solution of equations (66) and (67) is given as follows: + B sin(,~kt)
o,:(t) =A
B
al2(t ) = , 4 c o s ( ~ k t ) + ~
sin((okt)
(68)
+ Q(t)[ft:QX(,)~(~')Q(,)dl"]
Q'r(t) (76)
a( t ) -- O( t ) a ( to)QT( t )
+ Q(t)[ft:Q'r(.r)&(~.)Q(r ) d r ] Q T ( t ) (77)
(69)
where Q is an orthogonal rotation tensor satisfying
where A = o 2(to)
(70)
QQT = QVQ = I
B = kG = ½tok[o22(to) - on(to) ]
(71)
A----- a,2(t0)
(72)
First, the integration in eqns. (76) and (77) is performed using Q = R where R is obtained from the polar decomposition of the deformation gradient F. The expression RTdR transforms the components of d in a local basis that is rotating with respect to the fixed global basis, while R is the measure of this rotation. The integration is then performed with regard to an observer rotating with the basis. Premultiplication by R and postmultiplication by R r of the resulting integration at any time t, gives the Cauchy stress components in the required global basis. Using the polar decomposition theorem and the expression for F in eqn. (47), we obtain:
= ½a~k[a22(to) - an(t0) ]
(73)
and t o is the time parameter indicating the initiation of plastic flow. The corresponding solutions for oil and 0"22 are given by
0"11(t) = A fttok cos(~k~) d~'+ B f t sin(~kr) dT to
to
(74) o22(/) = -0.H(t)
(75)
Similar expressions hold for aH(t ) and a2z(t ). It is noted that this special case of the simple shear problem is similar to that obtained by Paulun and Pecherski (1985).
1 R = Vt4 + k2t 2
(78)
[-2-k/ kt -2 0
0 0
] (79)
0
Alternatively, eqns. (76) and (77) are integrated using Q where Q is obtained from Numerical integration of corotational rates in constitutive equations
The use of integration of material rates for tensors obtained from constitutive equations through the use of corotational rates leads to mathematical inconsistencies of the theory of Lie derivatives (Zhong-Heng, 1963, 1984; Yano, 1957). The integration procedure followed in this work is that proposed by Goddard and Miller (1966) and used later by Reed and Afluri (1985) and by
(~( t) = WQ(t), Q( to) ~ I
(80)
where W is the total spin tensor and is obtained from eqn. (51). Solving the differential equation in expression (80), for Q, we obtain:
[-sin(cos/) ,oos(si/)/i] 0
0
292
G.Z. Voyiadis, P.L Kattan / Finite deformation plasticity
The stresses o12 and O l l are plotted versus the strain in Figs. 7 and 8 respectively. The results are compared for the solution obtained directly by solving the differential equations numerically and the solution obtained by integrating the constitutive equations using an orthogonal rotation tensor. It is clear from Figs. 7 and 8 that using R as an orthogonal rotation tensor gives satisfactory resuits that are comparable with the numerical results obtained from the differential equations using a value for to = 0.15. However, it is noticed that when using Q, obtained from the total spin, the stress-strain curves show oscillations when reaching a strain of about 250%. Thus, it is recommended that R be used in the numerical integration of the constitutive equations. The analytical solution for the shear stress o12 is also shown in Fig. 7 for the special case discussed earher. This solution seems comparable to the solution of the differential equations. The difference in these two solutions is due to the absence of the parameters b and c from the equations used to obtain the analytical solution. In Figs. 9 and 10, the components of the shift in spatial coordinates a H and ~t12 are plotted against the strain. In this case also it is observed that using R gives a more satisfactory solution that is comparable to the numerical solution of the differential equations. Oscillations occur when using Q from the total spin, and the behavior is similar to that of the stresses th~ and o12.
constitutive relations using an appropriate corotational rate. In this case we conclude that the rotation tensor R obtained from the polar decomposition theorem should be used in the integration scheme. For the case when the rotation tensor is obtained from the total spin, it is observed that oscillations occur in the stress-strain curves. A closed form sohttion is also obtained for a special case where T12 = ~t12. This solution reduces to that given by Paulun and Pecherski (1985). It is observed that this solution is independent of the parameters b and c that govern the kinematic and isotropic hardening, respectively. We conclude that the results obtained by using the Runge-Kutta-Verner method are most rehable because they are based on the governing differential equations which reflect the true nature of the finite simple shear problem. Also these results are shown for different proposed expressions for the function to which were suggested by various authors. It should also be noted that all spin tensors are related, and therefore constitutive relations written in terms of one can easily be expressed in terms of others. A comprehensive account of this can be found in Mehrabadi and Nemat-Nasser (1987), and Nemat-Nasser (1983). It is recommended that in finite element implementations, the governing eqns. (43) and (44) should be directly integrated instead of using the corotational rates given by eqns. (42) and (6).
Acknowledgement Conclusion A constitutive model is formulated for finite strain plasticity. The formulation is given in spatial coordinates and incorporates both isotropic and kinematic hardening. As an example, the problem of finite simple shear is solved successfully using two different approaches. The first solution is obtained by directly integrating the governing differential equations numerically using a R u n g e - K u t t a - V e r n e r fifth and sixth order method. In this case we observe that the stress increases monotonically with the increase of the sheafing strain. The second approach is based on the integration of the
The research described in this paper was sponsored by the National Science Foundation under Grant MSM 8800832.
References Asaro, R.J., Advances in Solid Mechanics, Vol. 23, 1983. Atluri, S.N. (1984), On constitutive relations at f'mite strain: Hypo-elasticity and elasto-plasticity with isotropic or kinematic hardening, Comput. Methods App£ Mech. Eng. 43, 137-171. Chandra, A. and S. Mukhftrjee (1986) An analysis of large strain viscoplasticity problems including the effects of induced material anisotropy, J. Appl. Mech. $3, 77-82.
G.Z. Voyiadis, P.I. Kattan / F i n i t e deformation plasticity
Dafalias, Y.F. (1983), Corotational rates for kinematic hardening at large plastic deformations, J. Appl. Mech. 50 561-565. Dafalias, Y.F. (1985), A missing link in the macroscopic constitutive formulation of large plastic deformations,, in: A. Sawczuk and G. Bianchi, eds., Plasticity Toddy, Modellin$ Methods and Applications, proc. International Symposium on Current Trends and Results in Plasticity, CISM, Udine, June 27-30, 1983, Elsevier, London, New York, pp. 135-151. Dafalias, Y.F. (1987), Corotational integrals in constitutive formulations, J. Eng. Mech. ASCE 113 (12), 1967-1973. Dienes, J.K. (1979) On the analysis of rotation an~.steess rate in deforming bodies, Acta Mech. 32 217-232. Fressengeas, C. and A. Molinari (1984), Representations due comportement plastique anisotrope aux grandes deformations, Arch. Mech. 36. Goddasd, J.D. and C. Miller (1966), An inverse of the Janmann derivative and some applications in the rheology of viscoelastic fluids Rheoi. Acta 5, 177-184. Johnson, G.C. and D.J. Bammann (1984), A discussion of stress rates in finite deformation problems, Int. J. Solids Struct. 20, 725-737. Lee, E.H. (1981), Some comments on elastic-plastic analysis, Int. J. Solids Struct. 17, 859-872, 981. Lee, E.H., R.L. Mallet and T.B. Wertheimer (1983), Stress analysis for anisotropic hardening in finite deformation plasticity, J. Appl. Mech. 50, 554-569. Loret, B. (1983), On the effects of plastic rotation in the finite deformation of anisotropic elastoplastic materials, Mech. Mater. 2, 287-304. Mandel, J. (1973), Relations de comportment des milieux elastiques-plastiques et elastiques-viscoplastiques. Notion de Repere Directeur, in: A. Sawczuk, ed., Foundations of Plasticity, Noordhoff, Leiden, pp. 387-400. Mandel, J. (1981), Sur la definition de la vitesse de deformation elastique et sa relation avec la vittesse de constrainte, Int. J. Solids Struct. 17, 873-878. Mandel, J. (1982), Definition d'un repere privilegie pour l'etude des transformations anelastiques du polycristal, J. Mec. Theor. Appl. 1, 1-23. Mehrabadi, M.M. and S. Nemat-Nasser (1987), Some basic kinematical relations for f'mite deformations of continua, Mech. Mater. 6, 127-138. Moss, W.C. (1988), On instabilities in large deformation simple shear loading. Comput. Methods App£ Mech. Eng. 46, 329-338. Nagtegaal, J.C. and J.E. De Jong (1982), Some aspects of non-isotropic work-hardening in finite strain plasticity, in: E.H. Lee and R.L. Mallet, eds., Proc. Workshop on Plasticity of Metals at Finite Strain: Theory, Experiment, and Computation, Stanford University, 1981, Division of Applied Mechanics, Stanford University, Stanford, CA, and De-
293
partment of Mechanical Engineering. Aeronautical Engineering and Mechanics, R.P.I., Troy, NY, pp. 65-102. • Nemat-Nasser, S. (1979), Decomposition of strain measures and their rates in finite deformation elastoplasticity, Int. J. Solids Struct. 15, 155-166. Nemat-Nasser, S. (1983), On finite plastic flow of crystalline solids and geomaterials, J. Appl. Mech. 50, 1114-1126. Onat, E.T., (1982), Representation of inelastic behaviour in the presence of anisotropy and of finite deformations, in: B. Wilshire and R.J. Owen, eds., Recent Advances in Creep and Fracture of Engineering Materials and Structures, Pineridge Press, Swansea, U.K., pp. 231-264. Onat, E.T. (1984), Shear flow of kinematically hardening rigid-plastic materials, in: GJ. Dvorak and R.T. Shield eds., Mechanics of Material Behaviour, Elsevier, Amsterdam, pp. 311-324. Paulun, J.E. and R.B. Pecherski (1985), Study of corotational rates for kinematic hardening in finite deformation plasticity, Arch. Mech. 37. (6) 661-677. Pecherski, R.B. (1985). Discussion of sufficient condition of plastic flow IocaliTation, proc. International Symposium on Current Trends and Results in Plasticity, CISM, Udine, June 27-30, 1983, Eng. Fract. Mech. 21 (the Olszak Memorial Volume), 767-779. Reed, K.W. and S.N. Atluri (1985), Constitutive modeling and computational implementation for finite strain plasticity, Int. J. Plastic 1, 63-87. Simo, J.C. and K.S. Pister, (1984), Remarks on rate constitutive equations for finite deformation problems: Computational implications, Comput. Methods Appl. Eng. 46, 201-215. Voyiadjis, G.Z. (1984), Experimental determination of the material parameters of elasto-plastic, work-hardening metal alloys, Mater. Sci. Eng. J. NG2, 99-107. Voyiadjis, G.Z. and P.D. Kiousis (1987), Stress rate and the Lagrangian formulation of the Finite-strain plasticity for a yon Mises kinematic hardening model, Int. J. Solids Struct. 23, 95-109. Voyiadjis, G.Z. (1988), Degradation of elastic modulus in elasto-plastic coupling with finite strains, Int. J. Plastic. 4, 335-353. William, K.J., ed. (1984), Constitutive Equations: Macro and Computational Aspects, ASME Publication, Book No. G00274, 11, New York. 1984. Yano, K. (1957), The Theory of Lie Derivatives and Its Applications, North-Holland (1959), Amsterdam. Zhong-Heng. G. (1963), Time derivatives of tensor fields in non-linear continuum mechanics, Arch. Mech. Stos. 1 (15) 131-163. Zhong-Hen8, G. and R.H. Dubey (1984), Spins in a deforming continuum, Solids Mech. Arch. 9, 53-61. Ziegler, H. (1959), A modification of Prager's hardening rule, Q. Appl. Math. 17, 55-65.
279
EULERIAN CONSTITUTIVE MODEL FOR FINITE DEFORMATION PLASTICITY WITH ANISOTROPIC HARDENING George Z. VOYIADJIS and Peter I. KATTAN Department of Civil Engineering, Louisiana State University, Baton Rouge, 1.,,4 70803, U.S.A. Received 21 July 1988; revised version received 31 August 1988
A constitutive model is presented for finite strain plasticity. The model incorporates both isotropic and kinematic hardening of the Ziegler type. The corotational rate used here is in fine with the theory suggested by Paulun and Pecherski (1985) but not necessarily confined to the yon Mises type yield criterion and the Prager hardening rule. The aspect of integration of the corotational rates is also discussed here. The use of the integration of the material rate of tensors with time as a substitute for the proper integration with time of corotational rates leads to mathematical inconsistencies of the theory of Lie derivatives. The problem of simple shear is investigated and compared with other works.
Introduction The problem of oscillatory stresses produced in simple shear with the use of the Zaremba-Jaumann derivative in conjunction with large deformation kinematic hardening plasticity was pointed out by Nagtegaal and De Jong (1982). In an earlier paper, TruesdeU (1956) noticed the oscillation for simple shear and gave a closed form solution. Later, Dienes (1979) observed the same phenomenon using a hypo-elastic material model. The subject of finite deformation plasticity and in particular the objective of defining an appropriate corotational stress rate and an appropriate constitutive formulation for finite deformation has received wide attention since then. Lee et al. (1983) proposed a modified corotational rate using the spin of the principal direction of a with the largest absolute eigenvalue, a is the deviatoric component of the shift-stress tensor. An alternate approach by Onat (1982, 1984) defines the spatial spin equal to the antisymmetric part of d~'j:ajk multiplied by a constant, d " is the plastic component of the spatial strain rate. The nonoscillatory solution for simple shear is obtained by the proper choice of the constant. Dafalias (1983, 1985) and Loret (1983) obtained similar relations by associating the corotational rate with
the material substructure as defined by Mandel (1973). Mandel (1973) used the triad of director vectors attached to the material substructure and developed the theory of plasticity such that the substructure corotational rate is defined in terms of the spin of the director vectors. He postulated that the constitutive relations require not only the plastic component of the spatial strain rate but also the plastic component of the spatial spin tensor. Dafalias (1985) and Loret (1983) discussed the macroscopic constitutive relations for the plastic spin using the representation theorem for isotropic second-rank antisymmetric tensor valued functions. The importance of the material substructure in defining objective corotational rates is also argued by Pecherski (1985). In inelastic finite deformations of polyerystalline metals, the material moves with respect to the underlying crystal lattice. The lattice itself undergoes elastic deformations and relative rigid-body rotations due to the lattice misorientation (Pecherski, 1985). The work outlined above by Lee et ai. (1983), Onat (1982, 1984), Dafalias (1983, 1985), and Loret (1983), imposes a retardation of the material spin W in order to obtain a non-oscillatory solution for the simple shear problem. The analysis of the solution of the simple shear problem presented
0167-6636/89/$3.50 © 1989, Elsevier Science Publishers B.V. (North-Holland)
280
G.Z. Voyiadis, P.L Kattan / Finite deformation plasticity
by Onat (1982, 1984), Dafalias (1983, 1985), and Loret (1983), results in an unbounded non-oscillatory solution for the shear stress that increases monotonically with increased deformation. Concurrently, the normal stress approaches an asymptotic upper bound. In the case of Lee et al. (1983), both the shear and the normal stresses are unbounded and increase monotonically with increased deformation. We also note that in the case of Dafalias (1983, 1985), Onat (1982, 1984), and Loret (1983), the principal directionsto~ a tend toward the bisector direction of the plane coordinate axes while in the case of Lee et al. (1983), the maximum principal direction of a inclines towards the horizontal axis. The above proposed solutions fail in the proper prediction of the shear stress-shear strain characteristic and the Swift effect in torsion of thin-walled tubes (Paulun and Pecherski, 1985). Other authors have followed different approaches for the proper choice of the corotational objective stress rate. Atluri (1984) based on the idea of a complete hypo-elastic law, modifies the rate of the back stress equation for the case of a rigid-kinematic hardening plastic model. Fressengeas and Molinari (1984), Johnson and Bammann (1984), Simo and Pister (1984), Voyiadjis (1984), Voyiadjis and Kiousis (1987), Moss (1988), and Voyiadjis (1988) have also discussed different aspects of proper choice of the objective stress rate in finite deformation analysis. An ASME publication by William (1984) summarizes the debate on this subject. In this work a constitutive model is presented for finite strain plasticity. The model incorporates both isotropic and kinematic hardening. The modified Prager hardening rule proposed by Ziegler (1959), is used. The corotational rate used here is in line with the theory suggested by Paulun and Pecherski (1985). This theory is not only based on the question of how to avoid the unwanted oscillatory stresses in the simple shear problem, but also refers to the proper constitutive description of anisotropic hardening in the course of finite plastic deformation process. This theory is based on Mandel's (1973, 1981, 1982) theory of material substructure, and the plastic spin expressed by Dafalias (1985) and Loret (1983).
The generalized theory presented here leads to the satisfactory solution of the simple shear traction problem. This is obtained using a numerical solution to solve the governing differential equations. A closed form solution is also obtained for a special case of the differential equations. An alternate numerical integration scheme proposed by Goddard and Miller (1966) is also used here for the solution of the simple shear problem. This scheme is based on integration of the constitutive relations using an appropriate corotational rate. The use of the integration of the material rate of tensors with time as a substitute of the proper integration with time of corotational rates leads to mathematical inconsistencies of the theory of Lie derivatives (Zhong-Heng, 1963, 19841 Yano, 1957). For an isotropic linear hypo-elastic material, for deformations in which the directions of principal stretches remain fixed in the body, the procedure outlined by Ralph and Bathe (Williams. 1984) is equivalent to that of Goddard and Miller (1966). They have presented an equation for the inverse of the Jaumann derivative.
Theoretical formulation
The yield condition used in this work combines both isotropic hardening and kinematic hardening of the type proposed by Ziegler (1959)
(1) A particular form of this yield condition is the generalized von Mises type, with both kinematic and isotropic hardening such that /-
~ ( ~ , - a ~ , ) ( ~ , - ~ , ) - o: + c~ = 0
(2)
In the above expressions rkt is the deviatoric component of the Cauchy stress tensor Okt, a n d akt is the deviatoric component of the shift stress tensor. The constant o corresponds to the initial tensile or compressive yield stress from uniaxial loading, and c is a constant that describes the isotropic component of hardening. In eqns. (1) and (2) i¢, the isotropic hardening parameter, is obtained from the following relation: pp
~¢ = p k t d k t
(3)
G,Z. Voyiadis, P.L Kattan / Finite deformation plasticity
where d~; is the plastic component of the spatial strain rate tensor, and Pkt is a function of the Cauchy stress tensor and the accumulated elastic strain. The superdot in eqn. (3) denotes material time differentiation. Assuming small elastic strains but finite plastic deformations, the spatial rate is decomposed into (Nemat-Nasser, 1979, 1983; Lee, 1981) t
Pl
dkt = dkt + dkt
(4)
where d~t is the elastic component of the spatial strain rate. The above equation will be correct for any amount of elastic strain if the physics of elasto-plasticity is invoked, for example the case of single crystals. A thorough account of this is given by Asaro (1983). An associated flow rule of the type shown below is postulated a;; =
A aokl a/
(5)
in finite plastic deformation analysis. We note in Figure 3a of Paulun and Pecherski (1985) that the shear stress obtained from Paulun and Pecherski (1985), Onat (1982, 1984), Dafalias (1983, 1985) and Loret (1983) for the simple shear problem increases monotonically with increased deformation. On the contrary the shear stress obtained from Lee et al. (1983) approaches an asymptotic upper bound with increased deformation. We also note that the upper bounds for the normal component all of the back stress obtained by Loret (1983), Dafalias (1983, 1985) and Onat (1982, 1984) for the non-oscillatory solution of the problem of simple shear do not seem to be justified physically (Paulun and Pecherski, 1985). The generalization of the modified spin tensor 12 to the three-dimensional case is accomplished through the use of the theory of material substructure as defined by Mandel (1973). Defining 12" such that 12* = w -
where A is a positive scalar function. Following the Prager-Ziegler kinematic hardening rule, the evolution equation for a is given by &kt =/i
(~'k,- ak,)
(6)
where &kt is an objective rate of the back stress that is corotational with the substructure (Paulun and Pecherski, 1985). The objective rate for a is given by (Paulun and Pecherski, 1985). & = a, - 12a + a12
(7)
where 12 is a modified spin tensor. In simple shear 12 causes a retardation to the material spin W which consequently provides a non-oscillatory solution. This is accomplished through the use of an influence function ~o decreasing with the shear strain. The determination of the function ~o is such that the angular velocity ? corresponding to the spin W in simple shear is reduced to the angular velocity 6 of a single material line element. The objective rate of a given by eqn. (7) not only avoids the unwanted oscillatory stress in the simple shear problem, but also provides proper constitutive description of anisotropic hardening
281
(8)
12
we express eqn. (5) as follows:
a = 7, + 12"a - a12"
(9)
where g is the Zaremba-Jaumann derivative given by (10)
7, ffi dt - W a + a W
Using the general form of the yield function (1), the plastic spin may be expressed by the following equation proposed by Dafalias (1985) w " = ~ , ( ~ o - o~) + ~ ( ~ o
- o~)
+ ~ ( ~ 0 ~ - o~o) + ~ , ( ~ o ~ - ~ 0 ~ ) + "% ( o a o 2 - o2ao)
(11)
Using the general case of the yield function (1), we may express 8 f / a o U such that at" _ a/ 0.11 of 0.12 a/ aJ3 aaU aJI aoU + aJ2 a%j + aJ3 aoU
(12)
or
af ; 0%
af af aJ, au + 2 ~ %
af + 3~-~3o..,,.,
(13)
282
G.Z. Voyiadis, P.I. Kattan / Finite deformation plasticity
where J~ = %,, .I2 --- o,j~i, and -/3 = % ~ k O , V Making use of eqns. (13) and (5) we obtain:
The function ~7 is equation (19) for the case of simple shear yields (Paulun and Pecherski, 1985):
( a d " - d"ot ) rl = --~--(1 - ~o)5'
~f oo~)] = A[2~-~f4(ao-o. ) +3-~3(aoo-
or
(14) While retaining the first term only of eqn. (11), it is evident from eqn. (14) that in the absence of J3 terms, the plastic spin may be expressed as follows (Dafalias, 1983, 1985): W " --- f l ( a d " - d " a )
(15)
where 7ll = 2/Z'i ~-~2
(16)
However, eqn. (14) indicates that a particular form of expression (11) such that: W " = ~ , ( a o - oa) + ~13(ao 2 - o2a)
(17)
may be reduced to the form of eqn. (15) provided that the yield function is of such form that constraint (14) is satisfied. Constraint (14) may also be indicative that yield functions containing "/3 terms require that the plastic spin be expressed by eqn. (17) and not by setting ~!3 = 0. That is the 7/3 term needs to be present in the plastic spin expression. In the following work, it is assumed that the yield function does not contain J3 terms but is nevertheless not limited to the von Mises functional form. Equation (15) may be expressed in indicial notation as follows: tP
PP
tt
Wit = fl( a,,,,d,,,t- ak,na.,t)
(18)
where 13 is a scalar function of the isotropic invariants of a and o. The spin fl* is determined from the plastic spin I4,'" and is expressed as follows (Paulun and Pecherski, 1985): It* = n~'
~,2 2 1 + ~,2 ~'
n
y = ~/3-e"
(23)
where e'"
V~"JUklUkl
(24)
We may generalize to three-dimensional analysis the expression for ~ to be (Paulun and Pecherski, 1985)
=} 1 + 3(~") 2
~1 = A~
(26)
such that ~ - - 2~
3(E")2 ~ 2 1 + 3(E") ~
(20)
0 f -0f ] 8°*t 8°*t
(27)
Similarly, using eqn. (5) in expression (20), we eliminate A from eqn. (20) such that: akin
( t r [ ( a k . , d " t - d ~ , ' a . t ) = ] ) '/2
(25)
J
where
(19)
a,7,..,
(22)
where y and 5' are the shear strain and shear strain rate, respectively. The above equation represents the difference between the constant angular speed produced by the material spin W and the angular velocity of the material line element lying initially along the vertical axis. This expression relates to the angular velocity corresponding to the spin 12" and shows to be a function of the accumulated plastic strain "~. Since y may be expressed in terms of the scalar ~", such that
where u,t =
(21)
/ [t
8/ ~Oml
ao.,,
8/ ~Okn ant
ao,.""'
283
G.Z. Voyiadis, P.L Kattan / F i n i t e deformation plasticity
Substituting the expressions for u and 7/ from eqns. (20) and (25) or (28) and (26), respectively, into eqn. (19), the spin D* is determined. Equation (19) may now be expressed as follows: a k l = a k l -- ~Pkmrtmt + 71aknVnl
(29)
We note that in the absence of any plastic deformations, we obtain v ~ Otkl . Olkl
(30)
Equations (15), (18), (20), (21), (22), (25), (27) and (28) have been obtained on the assumption that the yield function (1) has the specific form expressed by eqn. (2). Before deriving the elasto-plastic stiffness matrix, we need to determine the positive parameter /~ associated with the kinematic hardening rule. This parameter is determined assuming that the projection of ,~ on the stress gradient of the yield surface equals to bd":
Of
0om.
bake=&"" Of Of
(31)
Ookt
OOpq OOpq where b is a material parameter. Substituting ,t and d " in eqn. (31) from eqns. (6) and respectively and post-multiplying the resulting pression by Of/OOkl, we obtain the required pression for/2: of of 0a,., 0am.
o/ ( q'pq -- apq )
Ek,m, = X8kA,, + a(Sk..8,. + 8 k : , , )
for (5) exex-
(32)
(35)
and ~ and G are Lam~'s constants. We note that the term - , / ( v a - ap) in eqn. (34) acts a damage effect, during unloading and reloading in the elastic range, due to the accumulated plastic strains. The presence of this damage effect is clearly indicated in the experimental results obtained by Voyiadjis (1984, 1988). We note in this paper the degradation of the elastic stiffness, in uniaxial loading, due to accumulated plastic strains. Eliminating d ' from eqn. (33) through the use of expressions (4) and (5), we obtain:
The parameter A can be calculated from the consistency condition
f(%,, ak,, K)ffi 0
0/
t~ = A b
where the modulus of elasticity Ekt.,. takes the following form:
(37)
The interpretation of eqn. (37) can be made on a consistent mathematical basis with the aid of the theory of Lie derivatives (Zhong-Heng, 1963, 1984; Yano, 1957). Equation (37) is expressed as: +
°akt
~ K= 0
(38)
where the material time derivative of a scalar function is equal to its Lie derivative with respect to the velocity vector v. Making use of equations (3), (6), (4), (5), (32) and (36) in eqn. (38), we obtain an expression for A:
OOvq (39)
In order to obtain the elasto-plastic stiffness tensor, a linear elastic relation is assumed between the corotational stress rate d and the elastic component of the spatial strain rate tensor d ' such that
dkt = ektm.d=.
a/ a%b
~%d
a~/(%/•
Okl = E k l m . d m n -- ~lvkpapl + l]akqPql
(34)
a~
Q = a/ E=b~ a/
(33)
or v
where
a~/)b
a~
al
~au
~au
a/
(¢sh - %,) a%~
(4o)
G.Z. Voyiadis, P.I. Kattan /Finite deformation plasticity
284
Substituting for A, from eqn. (39) into eqn. (36), we obtain the elasto-plastic stiffness tensor which corresponds to the loading function given by expression (1) (provided it does not contain ,/3 terms). The resulting elasto-plastic stiffness tensor is given by Oklmn -~- Eklrn n
1 a/ af Q a~"'-qEpqmnEklijOo,)
(42)
The corotational rates expressed by eqns. (42) and (6) may be integrated numerically using a procedure proposed by Goddard and Miller (1966). Alternatively, by making use of expressions of the type given by relation (7), eqns. (42) and (6) can be expressed as follows: dkl = D k t , . . d m . - Okp~pt q- ~kqOql
7/ I i l
(41)
and ~kl ~ Dktmndmn
~2
:7
Fig. 1. Simple shear problem and material properties of metal. Material properties of aluminum alloy 2024 T4: v = 0.25, E =10,600 ksi, c=138 ksi, b = 73 ksi, o = 56 ksi, k =l and h = At ~ 0.0001.
The deformation gradient F and the velocity gradient L are given by, respectively,
(43)
e--
(44)
~=
1 0 0
kt
O
1 0
0 1
(47)
0
0
(48)
0
0
and akl=~(,Tki--~kl)--akp~plq-
ekqaql
Equation (43) and (44) may now be used in numerical algorithms (finite element methods) to solve boundary value problems.
L =d+ W
The problem of finite simple shear
where d is the rate of deformation (symmetric part of L ) and W is the total spin (antisymmetric part of L ) then we obtain
The problem of simple shear in the x x direction as shown in Fig. 1, is given by the following displacement fidd U! = k t x 2
(45a)
u2 = 0
(45b)
u~ = 0
(45c)
Since,
d=7
v I = kx 2
(46a)
v2 = 0
(46b)
03 = 0
(46c)
0
0
and W=
where k is a constant representing the sheafing strain rate and t is the time. The corresponding velocity field is expressed by
(49)
-1 0
0 0
(51)
Using eqns. (40), (41) and (35) together with eqn.
(50), we obtain dii = kEijl2 - 9--~-k[(,eq--%q)Epq,2E, j k , ( , k , - - a k , )
Q
l
(52)
G.Z. Voyiadis, P.L Kattan / Finite deformation plasticity
together with (54) and (55) and the following expressing for 12
Similarly for & we have 18kGb Q (~',2--f,2)('rkt--fk,)
akt=
(53) ~ = ~ o W = - ~ - [ -01 1 00 1 0 0 0 0
From eqns. (52) and (53), the corotational rates are given explicitly as follows: 36kG 2 = 2Ga- ~('q2 Q
d = 2Gd
and 18kGb
[
TII -- fill
T=
T12 -- f12
0
'/'12 -- f12
1"22 -- f 2 2
0
0
0
1"33 - - f33
36kG2 (~'12 )T+ -~-S ° ~ 0112
a = 18kGb ~ (, ~ 1 : - f 1 2 ) T + - ~ S ~
(55)
a= --ff-(,,2- f, )r
(60)
we obtain the following system of differential equations in 0 and a, respectively
(54)
- al,) r
where
285
(61) (62)
where
]
2012
(56)
S °=
022-011
O22--O11
--2012
0
0
!] (63)
and likewise for S ~.
and Q = 9(b + 2G)('rkt- fkt)(~'kt- ak,) -
Numerical solution of the differential equations of simple shear
(57)
+ 3COeq(,rpq fpq) -
Using the following two equations 0 = # - o~ + ~o
(58)
& = & - f12 + 12a
(59)
The problem of simple shear is first solved numerically by using a Runge-Kutta-Verner fifth and sixth order method for the solution of the
250
to = 0 2
b
~, , , ~ . I
0
I
2
3
4 Y
Fig. 2. Variation of all vs. 7 for different values of to.
0.006895. ~0~
I
I
I
I
5
6
7
8
9
G.Z. Voyiadis, P.I. Kattan /Finite deformation plastwity
286
o~=0.15 600 w=0.20 500
400
~:0.80
co = 0 4 0
300
= 1.00 ZOO
I00 .( I psi - 0 . 0 0 6 8 9 5
o
I
I
I
I
I
2
4
3
M Po)
I
i
I
1
l
5
6
7
8
9
7
Fig. 3. Variation of a12 vs. 7 for different values of ~0.
600
500
400
300
200
tO0
J
O
( Ipui = 0 . 0 0 6 8 9 5 MPo)
I
I
I
I
I
2
3
4
I
I
I
I
I
5 Y
6
7
8
9
Fig. 4. Variation of o12 vs. 7 for different ~. Legend: ~ Dafalias (1985), Loret (1983), Onat (1984); ( ~ Paulun and Pecherski (1985); (~) Dienes (1979)0 Dafalias (1983); and (~) Lee, Mallet and Wertheimer (1983).
G.Z. Voytadis. P.I. Kattan /Finite deformation plasticity
287
300
250
200 0
0
,5o
®
N
I00
50 ( I psi = 0 . 0 0 6 8 9 5 I I
I 2
Z 3
I 4
MPo )
I 5
I 6
I
I
I
7
8
9
7"
Fig. 5. Variation of a]2 vs. y for different ~0. Legend: see Fig. 4.
® 150 140 150
120 I10
I00 9o = ,ac _
80
®
70
60 5o 40 30 20 ~
0
I
2
3
.
.
w
4
5
•
I "Pall
6
Y
Fig. 6. o n versus "r for different expression of ,o. Legend: see Fig. 4.
288
G.Z. Voyiadis. P.I. Kattan / Finite deformation plasticity
governing differential eqns. (61) and (62). The value of ¢ in Figs. 2 and 3 for these equations is assumed to be constant. Several values for ~ are used ranging from 0.05 to 1 in the solution of the problem. In Fig. 2 the normal stress oll is plotted versus the shearing strain Y for different values of ~. It is clear from the figure that for values of greater than 0.2 oscillations in the stress-strain curve occur within this region (up to 1000% strain). We note that as the value of ~ approa¢hes zero, ./., the shearing stress increases mdefimtely as the sheafing strain increases. In Fig. 3, the shear stress ox2 is plotted versus the sheafing strain 7 for different values of ~. This figure also shows that as t~ increases, oscillations occur in the curve as explained earlier. It is interesting to note that for the case of t~ = 1, we obtain the solution for the Jaumann rate which is oscillating as shown in Fig. 3. We note in this •
figure that all values of ~ give similar results up to a strain of about 100%. In Figs. 4 - 6 the sheafing stress or2, the shift component a~2 and the normal stress o H are plotted respectively versus the strain V for different proposed expressions for ~. The numerical solution obtained for o12, o H, a H and a12 is shown as curve (4) in Figs. 7,8,9, and 10 respectively. It is noted that in all the curves in Figs. 4 and 5 oa2 and at2 increase monotonically without bound as 7 increases and no oscillations occur. Nevertheless, the expressions proposed by Dafalias (1985), Loret (1983), and Onat (1984) as well as Paulun and Pecherski (1985) seem to give better correspondence with the results shown in Figs. 7 and 10 for curves 3 and 4. On the other hand, Fig. 6 shows that for the stress o~, the expressions proposed by Dafalias (1983) and Dienes (1979) give results that are closely related to the results
600
500
,=~oo
~oo
®
200
I00 /"
0
( I psi = 0.006895 MPo) I
[
I
I
I
[
I
I
I
i
2
3
4
5
6
7
8
9
Fig. 7. Variation of Ola vs. Y for different solutions. Legend: ( ~ Closed form solution for the special case with w = 0.15; Numerical integration of constitutive equations with (~ - WQ; ( ~ Numerical integration of constitutive equations with Q = R; and ( ~ Numerical integration of differential equation with w = 0.15.
G.Z. Voyiadis, P.L Kattan /Finite deformation plasticity
289
ZSO
® 20C
( I p s i =0.
eD
-I< 15C b I00
®
50
,
0
Fig.
I
2
I
I
I
I
I
3
4
5
6
7
)'
8. Variation o f o i l vs. "f for different solutions. Legend: see
Fig.
J
8
9
7.
® J6oI 150 140 130 120
~
90
~
SO -
7O 6O 5O 40 3O 2O
IO
0
I
2
3
4
)'
5
6
Fig. 9. Variation o f a]; vs. ¥ for different solutions. Legend: see Fig. 7
7
8
9
290
G.Z. Voyiadis, P.L Kattan /Finite deformation plasticn)
25(
®
a~x
~
150
co I 0 0 =to
5O (
10si = 0 . 0 0 6 8 9 5
MPa)
® o
-50
I
2
3 " ~
~
7
8
9
~"
y
FIGURE I0 • VARIATION OF' a IZ vs. 7 FOR DIFFERENT SOLUTIONS Fig. 10. Variation o f a12 vs. 7 for different solutions. Legend: see Fig. 7.
obtained in Fig. 8 by numerically integrating the constitutive equations with the corotational rates given in eqns. (42) and (6). There is a special case of the differential eqns. (61) and (62) where an analytical solution can be obtained. This case is based on the assumption that ~'12= a~2. This assumption is one of two possible ways of achieving 033 = a33 = 0 in eqns. (61) and (62). Next, we investigate the nature of the solution for this case. The other possible alternative of obtaining o 3 = 0 is to assume T33 -- a33. However, this case cannot be solved analytically but can be solved numerically similarly to the solution presented in this section.
Special case of the analytical solution of the differential equations of simple shear For this case the system of eqns. (61) and (62) reduce to an uncoupled system of equations given
by
o 12
=
~
olfo,, /
-1
+
o
0
0"22]
-2
1/ o1: OJlo2z ! (64)
and 0
0t22 ]
0
0
1
at2 ~
- 2
0
0~22/
(65)
We note that in both systems of eqns. (64) and (65), the material parameters c and b do not appear explicitly. They can be incorporated in an appropriate expression for o~. But since ~ is assumed to be constant, this solution does not incorporate the material parameters b and c. F r o m eqns. (64) and (65), we obtain the following two uncoupled second order differential equa-
291
G.Z. Voyiadis, P.L Kattan / Finite deformation plasticity
tions: (66)
O12 q- ~2k2°12 --- 0
Chandra and Mukherjee (1986). This procedure is given by the integral:
a( t ) = Q( t )a( to)Qr( t )
and (67)
fi;12 -I- 602k2¢112 = 0
The solution of equations (66) and (67) is given as follows: + B sin(,~kt)
o,:(t) =A
B
al2(t ) = , 4 c o s ( ~ k t ) + ~
sin((okt)
(68)
+ Q(t)[ft:QX(,)~(~')Q(,)dl"]
Q'r(t) (76)
a( t ) -- O( t ) a ( to)QT( t )
+ Q(t)[ft:Q'r(.r)&(~.)Q(r ) d r ] Q T ( t ) (77)
(69)
where Q is an orthogonal rotation tensor satisfying
where A = o 2(to)
(70)
QQT = QVQ = I
B = kG = ½tok[o22(to) - on(to) ]
(71)
A----- a,2(t0)
(72)
First, the integration in eqns. (76) and (77) is performed using Q = R where R is obtained from the polar decomposition of the deformation gradient F. The expression RTdR transforms the components of d in a local basis that is rotating with respect to the fixed global basis, while R is the measure of this rotation. The integration is then performed with regard to an observer rotating with the basis. Premultiplication by R and postmultiplication by R r of the resulting integration at any time t, gives the Cauchy stress components in the required global basis. Using the polar decomposition theorem and the expression for F in eqn. (47), we obtain:
= ½a~k[a22(to) - an(t0) ]
(73)
and t o is the time parameter indicating the initiation of plastic flow. The corresponding solutions for oil and 0"22 are given by
0"11(t) = A fttok cos(~k~) d~'+ B f t sin(~kr) dT to
to
(74) o22(/) = -0.H(t)
(75)
Similar expressions hold for aH(t ) and a2z(t ). It is noted that this special case of the simple shear problem is similar to that obtained by Paulun and Pecherski (1985).
1 R = Vt4 + k2t 2
(78)
[-2-k/ kt -2 0
0 0
] (79)
0
Alternatively, eqns. (76) and (77) are integrated using Q where Q is obtained from Numerical integration of corotational rates in constitutive equations
The use of integration of material rates for tensors obtained from constitutive equations through the use of corotational rates leads to mathematical inconsistencies of the theory of Lie derivatives (Zhong-Heng, 1963, 1984; Yano, 1957). The integration procedure followed in this work is that proposed by Goddard and Miller (1966) and used later by Reed and Afluri (1985) and by
(~( t) = WQ(t), Q( to) ~ I
(80)
where W is the total spin tensor and is obtained from eqn. (51). Solving the differential equation in expression (80), for Q, we obtain:
[-sin(cos/) ,oos(si/)/i] 0
0
292
G.Z. Voyiadis, P.L Kattan / Finite deformation plasticity
The stresses o12 and O l l are plotted versus the strain in Figs. 7 and 8 respectively. The results are compared for the solution obtained directly by solving the differential equations numerically and the solution obtained by integrating the constitutive equations using an orthogonal rotation tensor. It is clear from Figs. 7 and 8 that using R as an orthogonal rotation tensor gives satisfactory resuits that are comparable with the numerical results obtained from the differential equations using a value for to = 0.15. However, it is noticed that when using Q, obtained from the total spin, the stress-strain curves show oscillations when reaching a strain of about 250%. Thus, it is recommended that R be used in the numerical integration of the constitutive equations. The analytical solution for the shear stress o12 is also shown in Fig. 7 for the special case discussed earher. This solution seems comparable to the solution of the differential equations. The difference in these two solutions is due to the absence of the parameters b and c from the equations used to obtain the analytical solution. In Figs. 9 and 10, the components of the shift in spatial coordinates a H and ~t12 are plotted against the strain. In this case also it is observed that using R gives a more satisfactory solution that is comparable to the numerical solution of the differential equations. Oscillations occur when using Q from the total spin, and the behavior is similar to that of the stresses th~ and o12.
constitutive relations using an appropriate corotational rate. In this case we conclude that the rotation tensor R obtained from the polar decomposition theorem should be used in the integration scheme. For the case when the rotation tensor is obtained from the total spin, it is observed that oscillations occur in the stress-strain curves. A closed form sohttion is also obtained for a special case where T12 = ~t12. This solution reduces to that given by Paulun and Pecherski (1985). It is observed that this solution is independent of the parameters b and c that govern the kinematic and isotropic hardening, respectively. We conclude that the results obtained by using the Runge-Kutta-Verner method are most rehable because they are based on the governing differential equations which reflect the true nature of the finite simple shear problem. Also these results are shown for different proposed expressions for the function to which were suggested by various authors. It should also be noted that all spin tensors are related, and therefore constitutive relations written in terms of one can easily be expressed in terms of others. A comprehensive account of this can be found in Mehrabadi and Nemat-Nasser (1987), and Nemat-Nasser (1983). It is recommended that in finite element implementations, the governing eqns. (43) and (44) should be directly integrated instead of using the corotational rates given by eqns. (42) and (6).
Acknowledgement Conclusion A constitutive model is formulated for finite strain plasticity. The formulation is given in spatial coordinates and incorporates both isotropic and kinematic hardening. As an example, the problem of finite simple shear is solved successfully using two different approaches. The first solution is obtained by directly integrating the governing differential equations numerically using a R u n g e - K u t t a - V e r n e r fifth and sixth order method. In this case we observe that the stress increases monotonically with the increase of the sheafing strain. The second approach is based on the integration of the
The research described in this paper was sponsored by the National Science Foundation under Grant MSM 8800832.
References Asaro, R.J., Advances in Solid Mechanics, Vol. 23, 1983. Atluri, S.N. (1984), On constitutive relations at f'mite strain: Hypo-elasticity and elasto-plasticity with isotropic or kinematic hardening, Comput. Methods App£ Mech. Eng. 43, 137-171. Chandra, A. and S. Mukhftrjee (1986) An analysis of large strain viscoplasticity problems including the effects of induced material anisotropy, J. Appl. Mech. $3, 77-82.
G.Z. Voyiadis, P.I. Kattan / F i n i t e deformation plasticity
Dafalias, Y.F. (1983), Corotational rates for kinematic hardening at large plastic deformations, J. Appl. Mech. 50 561-565. Dafalias, Y.F. (1985), A missing link in the macroscopic constitutive formulation of large plastic deformations,, in: A. Sawczuk and G. Bianchi, eds., Plasticity Toddy, Modellin$ Methods and Applications, proc. International Symposium on Current Trends and Results in Plasticity, CISM, Udine, June 27-30, 1983, Elsevier, London, New York, pp. 135-151. Dafalias, Y.F. (1987), Corotational integrals in constitutive formulations, J. Eng. Mech. ASCE 113 (12), 1967-1973. Dienes, J.K. (1979) On the analysis of rotation an~.steess rate in deforming bodies, Acta Mech. 32 217-232. Fressengeas, C. and A. Molinari (1984), Representations due comportement plastique anisotrope aux grandes deformations, Arch. Mech. 36. Goddasd, J.D. and C. Miller (1966), An inverse of the Janmann derivative and some applications in the rheology of viscoelastic fluids Rheoi. Acta 5, 177-184. Johnson, G.C. and D.J. Bammann (1984), A discussion of stress rates in finite deformation problems, Int. J. Solids Struct. 20, 725-737. Lee, E.H. (1981), Some comments on elastic-plastic analysis, Int. J. Solids Struct. 17, 859-872, 981. Lee, E.H., R.L. Mallet and T.B. Wertheimer (1983), Stress analysis for anisotropic hardening in finite deformation plasticity, J. Appl. Mech. 50, 554-569. Loret, B. (1983), On the effects of plastic rotation in the finite deformation of anisotropic elastoplastic materials, Mech. Mater. 2, 287-304. Mandel, J. (1973), Relations de comportment des milieux elastiques-plastiques et elastiques-viscoplastiques. Notion de Repere Directeur, in: A. Sawczuk, ed., Foundations of Plasticity, Noordhoff, Leiden, pp. 387-400. Mandel, J. (1981), Sur la definition de la vitesse de deformation elastique et sa relation avec la vittesse de constrainte, Int. J. Solids Struct. 17, 873-878. Mandel, J. (1982), Definition d'un repere privilegie pour l'etude des transformations anelastiques du polycristal, J. Mec. Theor. Appl. 1, 1-23. Mehrabadi, M.M. and S. Nemat-Nasser (1987), Some basic kinematical relations for f'mite deformations of continua, Mech. Mater. 6, 127-138. Moss, W.C. (1988), On instabilities in large deformation simple shear loading. Comput. Methods App£ Mech. Eng. 46, 329-338. Nagtegaal, J.C. and J.E. De Jong (1982), Some aspects of non-isotropic work-hardening in finite strain plasticity, in: E.H. Lee and R.L. Mallet, eds., Proc. Workshop on Plasticity of Metals at Finite Strain: Theory, Experiment, and Computation, Stanford University, 1981, Division of Applied Mechanics, Stanford University, Stanford, CA, and De-
293
partment of Mechanical Engineering. Aeronautical Engineering and Mechanics, R.P.I., Troy, NY, pp. 65-102. • Nemat-Nasser, S. (1979), Decomposition of strain measures and their rates in finite deformation elastoplasticity, Int. J. Solids Struct. 15, 155-166. Nemat-Nasser, S. (1983), On finite plastic flow of crystalline solids and geomaterials, J. Appl. Mech. 50, 1114-1126. Onat, E.T., (1982), Representation of inelastic behaviour in the presence of anisotropy and of finite deformations, in: B. Wilshire and R.J. Owen, eds., Recent Advances in Creep and Fracture of Engineering Materials and Structures, Pineridge Press, Swansea, U.K., pp. 231-264. Onat, E.T. (1984), Shear flow of kinematically hardening rigid-plastic materials, in: GJ. Dvorak and R.T. Shield eds., Mechanics of Material Behaviour, Elsevier, Amsterdam, pp. 311-324. Paulun, J.E. and R.B. Pecherski (1985), Study of corotational rates for kinematic hardening in finite deformation plasticity, Arch. Mech. 37. (6) 661-677. Pecherski, R.B. (1985). Discussion of sufficient condition of plastic flow IocaliTation, proc. International Symposium on Current Trends and Results in Plasticity, CISM, Udine, June 27-30, 1983, Eng. Fract. Mech. 21 (the Olszak Memorial Volume), 767-779. Reed, K.W. and S.N. Atluri (1985), Constitutive modeling and computational implementation for finite strain plasticity, Int. J. Plastic 1, 63-87. Simo, J.C. and K.S. Pister, (1984), Remarks on rate constitutive equations for finite deformation problems: Computational implications, Comput. Methods Appl. Eng. 46, 201-215. Voyiadjis, G.Z. (1984), Experimental determination of the material parameters of elasto-plastic, work-hardening metal alloys, Mater. Sci. Eng. J. NG2, 99-107. Voyiadjis, G.Z. and P.D. Kiousis (1987), Stress rate and the Lagrangian formulation of the Finite-strain plasticity for a yon Mises kinematic hardening model, Int. J. Solids Struct. 23, 95-109. Voyiadjis, G.Z. (1988), Degradation of elastic modulus in elasto-plastic coupling with finite strains, Int. J. Plastic. 4, 335-353. William, K.J., ed. (1984), Constitutive Equations: Macro and Computational Aspects, ASME Publication, Book No. G00274, 11, New York. 1984. Yano, K. (1957), The Theory of Lie Derivatives and Its Applications, North-Holland (1959), Amsterdam. Zhong-Heng. G. (1963), Time derivatives of tensor fields in non-linear continuum mechanics, Arch. Mech. Stos. 1 (15) 131-163. Zhong-Hen8, G. and R.H. Dubey (1984), Spins in a deforming continuum, Solids Mech. Arch. 9, 53-61. Ziegler, H. (1959), A modification of Prager's hardening rule, Q. Appl. Math. 17, 55-65.