International Journal of Plasticity 25 (2009) 1833–1878
Contents lists available at ScienceDirect
International Journal of Plasticity journal homepage: www.elsevier.com/locate/ijplas
A large deformation theory for rate-dependent elastic–plastic materials with combined isotropic and kinematic hardening David L. Henann, Lallit Anand * Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
a r t i c l e
i n f o
Article history: Received 24 August 2008 Received in final revised form 25 November 2008 Available online 7 December 2008
Keywords: Large deformation viscoplasticity Isotropic hardening Kinematic hardening Time-integration procedure Finite elements
a b s t r a c t We have developed a large deformation viscoplasticity theory with combined isotropic and kinematic hardening based on the dual decompositions F ¼ Fe Fp [Kröner, E., 1960. Allgemeine kontinuumstheorie der versetzungen und eigenspannungen. Archive for Rational Mechanics and Analysis 4, 273–334] and Fp ¼ Fpen Fpdis [Lion, A., 2000. Constitutive modelling in finite thermoviscoplasticity: a physical approach based on nonlinear rheological models. International Journal of Plasticity 16, 469–494]. The elastic distortion Fe contributes to a standard elastic free-energy wðeÞ , while Fpen , the energetic part of Fp , contributes to a defect energy wðpÞ – these two additive contributions to the total free energy in turn lead to the standard Cauchy stress and a back-stress. Since Fe ¼ FFp1 and Fpen ¼ Fp Fp1 dis , the evolution of the Cauchy stress and the backstress in a deformation-driven problem is governed by evolution equations for Fp and Fpdis – the two flow rules of the theory. We have also developed a simple, stable, semi-implicit timeintegration procedure for the constitutive theory for implementation in displacement-based finite element programs. The procedure that we develop is ‘‘simple” in the sense that it only involves the solution of one non-linear equation, rather than a system of non-linear equations. We show that our time-integration procedure is stable for relatively large time steps, is first-order accurate, and is objective. Ó 2008 Elsevier Ltd. All rights reserved.
1. Introduction In classical small deformation theories of metal plasticity, when attempting to model the Bauschinger phenomenon (Bauschinger, 1886) and other phenomena associated with cyclic loading, kine* Corresponding author. Tel.: +1 617 253 1635. E-mail address:
[email protected] (L. Anand). 0749-6419/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijplas.2008.11.008
1834
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
matic strain hardening is often invoked by introducing a symmetric and deviatoric stress-like tensorial internal variable Tback called the back-stress, which acts to oppose the Cauchy stress T in the formulation of an appropriate yield criterion and a corresponding flow rule for the material. The simplest model for the evolution of such a back-stress is Prager’s linear kinematic-hardening rule (Prager, 1949):
T_ back ¼ BE_ p
so that Tback ¼ BEp ;
ð1:1Þ p
where B > 0 is a back-stress modulus and E is the plastic stain tensor in the standard decomposition E ¼ Ee þ Ep of the infinitesimal strain tensor E into elastic and plastic parts. A better description of kinematic hardening is given by an evolution equation for the back-stress Tback proposed initially by Armstrong and Frederick (1966):1
T_ back ¼ BE_ p cTback e_ p ;
def with e_ p ¼ jE_ p j;
ð1:2Þ
where the term ðcTback e_ p Þ with c > 0 represents a dynamic recovery term, which is collinear with the back-stress Tback and is proportional to the equivalent plastic strain rate e_ p . For monotonic uniaxial loading, the evolution of back-stress, instead of being linear as prescribed by (1.1), is then exponential with a saturation value. The Armstrong–Frederick model for kinematic hardening is widely used as a simple model to describe cyclic loading phenomena in metals. For a discussion of the Armstrong–Frederick model and its various extensions to model cyclic plasticity in the small deformation regime see Lemaitre and Chaboche (1990) and the recent review by Chaboche (2008). While there have been numerous attempts to extend the Armstrong–Frederick rule to large deformations, a suitable extension is still not widely agreed upon. Dettmer and Reese (2004) have recently reviewed and numerically analyzed several previous models for large deformation kinematic hardening of the Armstrong–Frederick type. One class of theories (cf., Dettmer and Reese, 2004, for a list of references to the literature) follow Prager and Armstrong and Frederick, and introduce a stress-like internal variable to model kinematic hardening; these theories require suitable frame-indifferent evolution equations to generalize an evolution equation resembling (1.2) to large deformations. An alternative to such theories are theories based on a proposal by Lion (2000), who in order to account for energy storage mechanisms associated with plastic flow, introduced a kinematic constitutive assumption that the plastic distortion Fp in the standard Kröner (1960) elastic–plastic decomposition F ¼ Fe Fp of the total deformation gradient, may be further multiplicatively decomposed into energetic and dissipative parts Fpen and Fpdis , respectively, as Fp ¼ Fpen Fpdis (cf., Lion’s Fig. 5).2 The elastic distortion Fe contributes to a standard elastic free-energy wðeÞ , while Fpen contributes a defect energy wðpÞ – these two additive contributions to the total free energy in turn lead to the standard Cauchy stress and a back-stress. Since Fe ¼ FFp1 and Fpen ¼ Fp Fp1 dis , the evolution of the Cauchy stress and the back-stress in a deformation-driven problem is governed by evolution equations for Fp and Fpdis – the two flow rules of the theory. An important characteristic of such a theory is that since Fp is invariant under a change in frame, then so also are Fpen and Fpdis , and the resulting evolution equations for the back-stress are automatically frame-indifferent, and do not require special considerations such as those required when generalizing equations of the type (1.2) to proper frame-indifferent counterparts. Examples of kinematic-hardening constitutive theories based on Lion-type Fp ¼ Fpen Fpdis decomposition may be found in Dettmer and Reese (2004), Menzel et al. (2005), Håkansson et al. (2005), Wallin and Ristinmaa (2005), and Vladimirov et al. (2008). While Dettmer and Reese (2004) in Section 4.3 of their paper, and more recently Vladimirov et al. (2008) in Section 2.2 of their paper have presented a generalization of the Armstrong–Frederick type of kinematic-hardening rate-independent plasticity theory using the Lion decomposition Fp ¼ Fpen Fpdis , we find their presentation rather brief, and (in our view) not sufficiently expository. The primary focus of the papers by Dettmer and Reese and Vladimirov et al., was the development of suitable time-integration algorithms, and they have developed several fully-implicit time-integration procedures for
1
Also recently published as Frederick and Armstrong (2007). Lion uses the notation FM F, Fe Fe , Fin Fp , Fine Fpen , Finp Fpdis . Based on a statistical-volume-averaging argument for the response of a polycrystalline aggregate, Clayton and McDowell (2003, Eq. (13)) have also introduced a Lion-type decomposition ~p ; they call F ~i F ~i Fp a ‘‘meso-incompatibility tensor.” Fp ¼ F en 2
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1835
their model; one such integration procedure involves the solution of a system of 13 highly non-linear equations – cf., Eq. (53) of Vladimirov et al. (2008). Motivated by the work of these authors, it is the purpose of this paper to: (1) Develop a thermodynamically consistent rate-dependent plasticity theory with combined isotropic and kinematic hardening of the Armstrong–Frederick type, based on the dual decompositions F ¼ Fe Fp and Fp ¼ Fpen Fpdis . The limit m ! 0, where m is a material rate-sensitivity parameter in our theory, corresponds to a rate-independent model. Importantly, we develop the theory based on the principle of virtual power and carefully lay down all assumptions and specializations that we have adopted so that it may in the future be possible to generalize the theory with other, more flexible isotropic and kinematic-hardening models. (2) Develop a simple, stable, semi-implicit time-integration procedure for our rate-dependent constitutive theory for implementation in displacement-based finite element programs. The procedure that we develop is ‘‘simple” in the sense that it involves the solution of only one stiff nonlinear equation,3 rather than a system of non-linear equations. We show that our time-integration procedure is stable for relatively large time steps, is first-order accurate, and is objective.
2. Kinematics We consider a homogeneous body B identified with the region of space it occupies in a fixed reference configuration, and denote by X an arbitrary material point of B. A motion of B is then a smooth oneto-one mapping x ¼ vðX; tÞ with deformation gradient, velocity and velocity gradient given by4
F ¼ rv;
_ v ¼ v;
_ 1 : L ¼ grad v ¼ FF
ð2:1Þ
To model the inelastic response of materials, we assume that the deformation gradient F may be multiplicatively decomposed as (Kröner, 1960)
F ¼ Fe Fp :
ð2:2Þ
As is standard, we assume that
J ¼ det F > 0; and consistent with this we assume that def
J e ¼ det Fe > 0; e
def
J p ¼ det Fp > 0;
ð2:3Þ
p
so that F and F are invertible. Here, suppressing the argument t: Fp ðXÞ represents the local inelastic distortion of the material at X due to a ‘‘plastic mechanism.” This local deformation carries the material into – and ultimately ‘‘pins” the material to – a coherent structure that resides in the structural space5 at X (as represented by the range of Fp ðXÞ); Fe ðXÞ represents the subsequent stretching and rotation of this coherent structure, and thereby represents the corresponding ‘‘elastic distortion.”
3 With reference to the power-law model (9.35) in Section 9.3, the stiffness of the equations depends on the strain-ratesensitivity parameter m, and the stiffness increases to infinity as m tends to zero, the rate-independent limit. For small values of m special care is required to develop stable constitutive time-integration procedures; cf., e.g., Lush et al. (1989). 4 Notation: We use standard notation of modern continuum mechanics. Specifically: r and Div denote the gradient and divergence with respect to the material point X in the reference configuration; grad and div denote these operators with respect to the point x ¼ vðX; tÞ in the deformed body; a superposed dot denotes the material time-derivative. Throughout, we write Fe1 ¼ ðFe Þ1 , Fp> ¼ ðFp Þ> , etc. We write trA, symA, skwA, A0 , and sym0 A, respectively, for the trace, symmetric, skew, deviatoric, and symmetric–deviatoric parts of a tensor A. Also, the inner product of tensors A and B is denoted by A : B, and the magnitude of A pffiffiffiffiffiffiffiffiffiffiffi by jAj ¼ A : A. 5 Also sometimes referred to as the ‘‘intermediate” or ‘‘relaxed” local space at X.
1836
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
We refer to Fp and Fe as the inelastic and elastic distortions. By (2.1)3 and (2.2),
L ¼ Le þ Fe Lp Fe1 ;
ð2:4Þ
Le ¼ F_ e Fe1 ;
ð2:5Þ
with
Lp ¼ F_ p Fp1 :
As is standard, we define the total, elastic, and plastic stretching and spin tensors through
9 D ¼ symL; W ¼ skwL; > = e e e e D ¼ symL ; W ¼ skwL ; > ; Dp ¼ symLp ; Wp ¼ skwLp ;
ð2:6Þ
so that L ¼ D þ W, Le ¼ De þ We , and Lp ¼ Dp þ Wp . The right and left polar decompositions of F are given by
F ¼ RU ¼ VR;
ð2:7Þ
where R is a rotation (proper-orthogonal tensor), while U and V are symmetric, positive-definite tensors with
U¼
pffiffiffiffiffiffiffiffi F> F;
V¼
pffiffiffiffiffiffiffiffi FF> :
ð2:8Þ
Also, the right and left Cauchy–Green tensors are given by
C ¼ U2 ¼ F> F;
B ¼ V2 ¼ FF> :
ð2:9Þ e
p
Similarly, the right and left polar decompositions of F and F are given by
Fe ¼ Re Ue ¼ Ve Re ; e
Fp ¼ Rp Up ¼ Vp Rp ;
p
e
e
ð2:10Þ p
p
where R and R are rotations, while U , V , U , V are symmetric, positive-definite tensors with
Ue ¼
pffiffiffiffiffiffiffiffiffiffiffiffi Fe> Fe ;
Ve ¼
pffiffiffiffiffiffiffiffiffiffiffiffi Fe Fe> ;
Up ¼
pffiffiffiffiffiffiffiffiffiffiffiffi Fp> Fp ;
Vp ¼
pffiffiffiffiffiffiffiffiffiffiffiffi Fp Fp> :
ð2:11Þ
Also, the right and left elastic Cauchy–Green tensors are given by
Ce ¼ Ue2 ¼ Fe> Fe ;
Be ¼ Ve2 ¼ Fe Fe> ;
ð2:12Þ
and the right and left plastic Cauchy–Green tensors are given by
Cp ¼ Up2 ¼ Fp> Fp ;
Bp ¼ Vp2 ¼ Fp Fp> :
ð2:13Þ
2.1. Incompressible, irrotational plastic flow We make two basic kinematical assumptions concerning plastic flow: (i) First, we make the standard assumption that plastic flow is incompressible, so that
J p ¼ det Fp ¼ 1 and tr Lp ¼ 0:
ð2:14Þ
Hence, using (2.2) and (2.14)1 ,
J e ¼ J:
ð2:15Þ
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1837
(ii) Second, from the outset we constrain the theory by limiting our discussion to circumstances under which the material may be idealized as isotropic. For isotropic elastic–viscoplastic theories utilizing the Kröner decomposition it is widely assumed that the plastic flow is irrotational in the sense that6
Wp ¼ 0:
ð2:16Þ p
p
Then, trivially, L D and
F_ p ¼ Dp Fp ;
with tr Dp ¼ 0:
ð2:17Þ
Thus, using (2.1), (2.4), (2.5), and (2.17), we may write (2.4) for future use as
grad v ¼ Le þ Fe Dp Fe1 ;
with trDp ¼ 0:
ð2:18Þ
2.2. Decomposition of Fp into energetic and dissipative parts Next following Lion (2000), in order to account for energy storage mechanisms associated with plastic flow, we introduce an additional kinematic constitutive assumption that Fp may be multiplicatively decomposed into an energetic Fpen and a dissipative part Fpdis as follows:
Fp ¼ Fpen Fpdis ;
det Fpen ¼ det Fpdis ¼ 1:
ð2:19Þ
Fpdis ðXÞ
Fpen p
We call the range of the local substructural space. Thus maps material elements from the substructural space to the structural space (which is the range of F ðXÞ). Lion (2000) calls the substructural space as the ‘‘intermediate configuration of kinematic hardening.” A ‘‘physical interpretation” of such an ‘‘intermediate configuration” is not completely clear7; here we simply treat it as a mathematical construct resulting from Lion’s presumed decomposition (2.19). Then, from (2.5)
Lp ¼ Lpen þ Fpen Lpdis Fp1 en
ð2:20Þ
with p1 Lpen ¼ F_ pen Fen
with trLpen ¼ 0;
p1 Lpdis ¼ F_ pdis Fdis
with trLpdis ¼ 0;
ð2:21Þ
and corresponding stretching and spin tensors through
6 A note on role of the plastic spin in isotropic solids: There are numerous publications discussing the role of plastic spin in theories based on the Kröner decomposition (cf., e.g., Dafalias, 1998). Based on the classical principle of frame-indifference (cf., Section 3), Green and Naghdi (1971) (cf., also Casey and Naghdi, 1980), introduced the notion of a change in frame of the intermediate structural space. This notion leads to transformation laws for Fe and Fp of the form
p F
¼ QFp ;
e F
¼ Fe Q > ;
in which Q ðX; tÞ is an arbitrary time-dependent rotation of the structural space. Unfortunately, Green and Naghdi viewed structural frame-indifference as a general principle; that is, a principle that stands at a level equivalent to that of conventional frameindifference. This view has been refuted by many workers (cf., e.g., Dashner, 1986, and the references therein). While we agree with the view that structural frame-indifference is not a general principle, this hypothesis may be used to represent an important facet of the behavior of a large class of polycrystalline materials based on the Kröner decomposition. Indeed, for polycrystalline materials without texture the structural space is associated with a collection of randomly oriented lattices, and hence the evolution of dislocations through that space should be independent of the frame with respect to which this evolution is measured. Using the notion of a change in frame of the structural space, Gurtin and Anand (2005, p. 1714) have shown that within a framework for isotropic plasticity based on the Kröner decomposition (which we follow here), one may assume, without loss in generality, that the plastic spin vanishes Wp ¼ 0. Here, we adopt the assumption Wp ¼ 0 from the outset; we do so based solely on pragmatic grounds: when discussing finite deformations for isotropic materials a theory without plastic spin is far simpler than one with plastic spin. 7 However, see Clayton and McDowell (2003) for a ‘‘statistical-volume-averaging” argument for the response of a polycrystalline aggregate.
1838
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
)
Dpen ¼ symLpen ;
Wpen ¼ skwLpen ;
Dpdis ¼ symLpdis ;
Wpdis ¼ skwLpdis ;
ð2:22Þ
so that Lpen ¼ Dpen þ Wpen and Lpdis ¼ Dpdis þ Wpdis . Paralleling the assumption of plastic irrotationality for isotropic materials, Wp ¼ 0, we assume that
Wpdis ¼ 0: Then, trivially
ð2:23Þ
Lpdis
¼
Dpdis ,
F_ pdis ¼ Dpdis Fpdis
and
with trDpdis ¼ 0:
ð2:24Þ
Since
; Wp ¼ Wpen þ skw Fpen Dpdis Fp1 en
ð2:25Þ
in order to enforce the constitutive constraint Wp ¼ 0, we require that both
Wpen ¼ 0;
ð2:26Þ
¼ 0: skw Fpen Dpdis Fp1 en
ð2:27Þ
and
Also, (2.20) reduces to
; Dp ¼ Dpen þ sym Fpen Dpdis Fp1 en
with trDpen ¼ trDpdis ¼ 0:
ð2:28Þ
The right and left polar decompositions of Fpen and Fpdis are given by
Fpen ¼ Rpen Upen ¼ Vpen Rpen ; where
Rpen
Upen
and
Rpdis
Fpdis ¼ Rpdis Updis ¼ Vpdis Rpdis ;
are rotations, while
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi p ¼ Fp> en Fen ;
Vpen
Upen ,
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ Fpen Fp> en ;
Vpen ,
Updis
Updis ,
Vpdis
ð2:29Þ are symmetric, positive-definite tensors with
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p ¼ Fp> dis Fdis ;
Vpdis ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Fpdis Fp> dis :
ð2:30Þ
Also, the corresponding right and left Cauchy–Green tensors are given by p> p Cpen ¼ Up2 en ¼ Fen Fen ;
Cpdis
¼
Up2 dis
¼
p Fp> dis Fdis ;
p p> Bpen ¼ Vp2 en ¼ Fen Fen ; p p> Bpdis ¼ Vp2 dis ¼ Fdis Fdis :
ð2:31Þ ð2:32Þ
3. Frame-indifference Changes in frame (observer) are smooth time-dependent rigid transformations of the Euclidean space through which the body moves. We require that the theory be invariant under such transformations, and hence under transformations of the form
v ðX; tÞ ¼ Q ðtÞðvðX; tÞ oÞ þ yðtÞ;
ð3:1Þ
with Q ðtÞ a rotation (proper-orthogonal tensor), yðtÞ a point at each t, and o a fixed origin. Then, under a change in frame, the deformation gradient transforms according to
F ¼ QF;
ð3:2Þ
and hence
C ¼ C that is; C is invariant;
ð3:3Þ
also F_ ¼ Q F_ þ Q_ F, and by (2.1)3 ,
L ¼ QLQ > þ Q_ Q > :
ð3:4Þ
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1839
Thus,
W ¼ QWQ > þ Q_ Q > :
D ¼ QDQ > ; e p
ð3:5Þ
e p
Moreover, ðF F Þ ¼ Q ðF F Þ, and therefore, since observers view only the deformed configuration,
Fe ¼ QFe ;
Fp ¼ Fp ;
and thus Fp is invariant:
ð3:6Þ
Hence, by (2.5),
Le ¼ QLe Q > þ Q_ Q > ;
ð3:7Þ
Lp ; Dp ; and Wp are invariant:
ð3:8Þ
and
Further, by (2.10),
Fe ¼ Re Ue ! QFe ¼ QRe Ue ; Fe ¼ Ve Re ! QFe ¼ QVe Q > QRe ; and we may conclude from the uniqueness of the polar decomposition that
Re ¼ QRe ;
Ve ¼ QVe Q > ; e
Ue ¼ Ue :
ð3:9Þ
e
Hence, from (2.12), B and C transform as
Be ¼ QBe Q >
and Ce ¼ Ce ;
ð3:10Þ
p
further, since F is invariant,
Bp and Cp are also invariant:
ð3:11Þ
p
Finally, since F is invariant under a change in frame,
Fpen and Fpdis are invariant; as are all of the following quantities based on Fpen and Fpdis :
Upen ;
Cpen ;
Vpen ;
Bpen ;
Rpen ;
Updis ;
Cpdis ;
Vpdis ;
Bpdis ;
Rpdis :
4. Development of the theory based on the principle of virtual power Following Gurtin and Anand (2005), the theory presented here is based on the belief that the power expended by each independent ‘‘rate-like” kinematical descriptor be expressible in terms of an associated force system consistent with its own balance. However, the basic ‘‘rate-like” descriptors, namely, v, Le , Dp , Dpen , and Dpdis are not independent, since by (2.18) and (2.28) they are constrained by
grad v ¼ Le þ Fe Dp Fe1 ;
trDp ¼ 0;
ð4:1Þ
and
; Dp ¼ Dpen þ sym Fpen Dpdis Fp1 en
trDpen ¼ trDpdis ¼ 0;
ð4:2Þ
and it is not apparent what forms the associated force balances should take. It is in such situations that the strength of the principle of virtual power becomes apparent, since the principle of virtual power automatically determines the underlying force balances.
1840
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
We denote by Pt an arbitrary part (subregion) of the deformed body with n the outward unit normal on the boundary @Pt of Pt . The power expended on Pt by material or bodies exterior to Pt results from a macroscopic surface traction tðnÞ, measured per unit area in the deformed body, and a macroscopic body force b, measured per unit volume in the deformed body, each of whose working accompanies the macroscopic motion of the body; the body force b presumed to account for inertia; that is, granted the underlying frame is inertial,
_ b ¼ b0 qv;
ð4:3Þ
with b0 the noninertial body force and q the mass density in the deformed body. We therefore write the external power as
Wext ðPt Þ ¼
Z
tðnÞ vda þ
@Pt
Z
b vdv:
ð4:4Þ
Pt
We assume that power is expended internally by an elastic stress Se power-conjugate to Le , a plastic stress Tp power-conjugate to Dp , an energetic stress Sen conjugate to Dpen , and a dissipative stress Tdis p1 conjugate to symðFpen Dpdis Fen Þ. Since Dp , Dpen , and symðFpen Dpdis Fp1 en Þ are symmetric and deviatoric, we assume that Tp , Sen , and Tdis are symmetric and deviatoric. We write the internal power as
Wint ðPt Þ ¼
Z
Se : Le þ J 1 Tp : Dp þ Sen : Dpen þ Tdis : symðFpen Dpdis Fp1 dv; en Þ
Pt
¼
Z
p> p Se : Le þ J 1 Tp : Dp þ Sen : Dpen þ sym0 ðFp> dv: en Tdis Fen Þ : Ddis
Pt
The term J arises because Tp : Dp þ Sen : Dpen þ Tdis : sym Fpen Dpdis Fp1 is measured per unit volume in en the structural space, but the integration is carried out within the deformed body. Defining a new stress measure 1
def
p> Sdis ¼ sym0 ðFp> en Tdis Fen Þ;
ð4:5Þ
the internal power may more compactly be written as
Wint ðPt Þ ¼
Z
Pt
Se : Le þ J 1 Tp : Dp þ Sen : Dpen þ Sdis : Dpdis dv:
ð4:6Þ
4.1. Principle of virtual power Assume that, at some arbitrarily chosen but fixed time, the fields v, Fe , Fpen (and hence F, Fp , Fpdis ) are known, and consider the fields v, Le , Dp , Dpen , and Dpdis as virtual velocities to be specified independently ~e , ~, L in a manner consistent with the constraints (4.1) and (4.2). That is, denoting the virtual fields by v ~ pen , and D ~ p to differentiate them from fields associated with the actual evolution of the body, we ~ p, D D dis require that
~ p Fe1 ; ~ ¼ L~e þ Fe D grad v and
~ p ¼ 0; trD
~p ¼ D ~ p p1 ; ~ p þ sym Fp D D en en dis Fen
ð4:7Þ
~ p ¼ trD ~ p ¼ 0: trD en dis
ð4:8Þ
More specifically, we define a generalized virtual velocity to be a list
~ p; D ~p ; D ~ p Þ; ~; L~e ; D V ¼ ðv en dis consistent with (4.7) and (4.8). Then, writing
Wext ðPt ; VÞ ¼ Wint ðPt ; VÞ ¼
Z Z
~ da þ tðnÞ v @Pt
Pt
Z
~dv; bv Pt
~e þ J 1 Tp : D ~ p þ Sen : D ~ p þ Sdis : D ~p Se : L dv; en dis
) ð4:9Þ
1841
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
respectively, for the external and internal expenditures of virtual power, the principle of virtual power is the requirement that the external and internal powers be balanced. That is, given any part Pt ,
Wext ðPt ; VÞ ¼ Wint ðPt ; VÞ for all generalized virtual velocities V:
ð4:10Þ
4.1.1. Consequences of frame-indifference We assume that the internal power Wint ðPt ; VÞ is invariant under a change in frame and that the virtual fields transform in a manner identical to their nonvirtual counterparts. Then given a change in frame, invariance of the internal power requires that
Wint ðPt ; V Þ ¼ Wint ðPt ; VÞ;
ð4:11Þ
Wint ðPt ; V Þ
Pt
where and represent the region and the internal power in the new frame. In the new frame Pt transforms rigidly to a region Pt , the stresses Se , Tp , Sen , and Sdis transform to Se , Tp , Sen , and ~e transforms to Sdis , while L
~ e Q > þ Q_ Q > ; L~e ¼ Q L ~ p, D ~ pen , and D ~ p are invariant. Thus, under a change in frame Wint ðPt ; VÞ transforms to and D dis
Wint ðPt ; V Þ ¼ Wint ðPt ; V Þ Z ~ e Q > þ Q_ Q > þ J 1 Tp : D ~ p þ S : D ~ p þ S : D ~p dv ¼ Se : Q L en en dis dis P Z t ~ p þ S : D ~ p þ S : D ~p ¼ dv: ðQ > Se Q Þ : L~e þ Se : ðQ_ Q > Þ þ J 1 Tp : D en en dis dis Pt
Then (4.11) implies that
Z
~ p þ S : D ~ p þ S : D ~p ðQ > Se Q Þ : L~e þ Se : ðQ_ Q > Þ þ J 1 Tp : D dv en en dis dis Pt Z ~ p þ Sen : D ~ p þ Sdis : D ~p dv; ¼ Se : L~e þ J 1 Tp : D en dis
ð4:12Þ
Pt
or equivalently, since the part Pt is arbitrary,
~ p þ S : D ~ p þ S : D ~p ðQ > Se Q Þ : L~e þ Se : ðQ_ Q > Þ þ J 1 Tp : D en en dis dis ~ p þ Sen : D ~ p þ Sdis : D ~p : ¼ Se : L~e þ J 1 Tp : D en dis
ð4:13Þ
Also, since the change in frame is arbitrary, if we choose it such that Q is an arbitrary time-independent rotation, and that Q_ ¼ 0, we find from (4.13) that
h i ~ p þ S Sdis : D ~ p ¼ 0: ~ p þ S Sen : D ðQ > Se Q Þ Se : L~e þ J 1 ðTp Tp Þ : D en en dis dis
~e , D ~ p, D ~ pen , and D ~ p , we find that the stress Se transforms according to Since this must hold for all L dis
Se ¼ QSe Q > ;
ð4:14Þ
p
while T , Sen , and Sdis are invariant
Tp ¼ Tp ;
Sen ¼ Sen ;
Sdis ¼ Sdis :
ð4:15Þ
Next, if we assume that Q ¼ 1 at the time in question, and that Q_ is an arbitrary skew tensor, we find from (4.13), using (4.14), (4.15) that
Se : Q_ ¼ 0; or that the tensor Se is symmetric,
Se ¼ Se> :
ð4:16Þ e
Thus, the elastic stress S is frame-indifferent and symmetric.
1842
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
4.1.2. Macroscopic force balance In applying the virtual balance (4.10) we are at liberty to choose any V consistent with the con~ is arbitrary, and straints (4.7) and (4.8). Consider a generalized virtual velocity with V for which v ~ pen ¼ D ~ p 0, so that ~p ¼ D D dis
~: L~e ¼ grad v
ð4:17Þ
For this choice of V, (4.10) yields
Z
~da þ tðnÞ v
@Pt
Z
~dv ¼ bv
Pt
Z
~dv: Se : grad v
ð4:18Þ
Pt
Further, using divergence theorem, (4.18) yields
Z
~da þ ðtðnÞ Se nÞ v
@Pt
Z
~dv ¼ 0: ðdivSe þ bÞ v
ð4:19Þ
Pt
~, standard variational arguments yield the traction condition Since (4.19) must hold for all Pt and all v
tðnÞ ¼ Se n;
ð4:20Þ
and the local force balance
divSe þ b ¼ 0:
ð4:21Þ e
This traction condition and force balance and the symmetry and frame-indifference of S are classical conditions satisfied by the Cauchy stress T, an observation that allows us to write def
T ¼ Se
ð4:22Þ
and to view
T ¼ T>
ð4:23Þ
as the standard macroscopic Cauchy stress and (4.21) as the local macroscopic force balance. Since we are working in an inertial frame, so that (4.3) is satisfied, (4.21) reduces to the local balance law for linear momentum,
_ divT þ b0 ¼ qv;
ð4:24Þ
with b0 the noninertial body force. 4.1.3. Microscopic force balances To discuss the microscopic counterparts of the macroscopic force balance, first consider a generalized virtual velocity with
~ p 0; ~ 0 and D v dis
ð4:25Þ
~ p arbitrarily, and let also, choose the virtual field D
~ p Fe1 L~e ¼ Fe D
~p ¼ D ~ p; and D en
ð4:26Þ
consistent with (4.7) and (4.8). Thus
~ p: T : L~e ¼ ðFe> TFe> Þ : D
ð4:27Þ
Next, define an elastic Mandel stress by def
Me ¼ JFe> TFe> ;
ð4:28Þ
which in general is not symmetric. Then, on account of our choice (4.25)1 the external power vanishes, so that, by (4.10), the internal power must also vanish, and satisfy
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
Z
Wint ðPt ; VÞ ¼
1843
~ p dv ¼ 0: J 1 ðTp þ Sen Me Þ : D
Pt
~ p , a standard argument Since this must be satisfied for all Pt and all symmetric and deviatoric tensors D yields the first microforce balance
ð4:29Þ Next consider a generalized virtual velocity with
~ p 0; ~ 0 and D v
ð4:30Þ
~ p arbitrarily, and let also, choose the virtual field D dis
~ p ¼ sym Fp D ~ p p1 D en en dis Fen
ð4:31Þ
consistent with (4.8). Thus
~ p ¼ ðFp> Sen Fp> Þ : D ~p : Sen : D en en en dis
ð4:32Þ
Next, define a new stress measure by def
p> Mpen ¼ Fp> en Sen Fen ;
ð4:33Þ Mpen
which in general is also not symmetric. For lack of better terminology, we call a plastic Mandel stress. Then, on account of our choice (4.30)1 the external power vanishes, so that, by (4.10), the internal power must also vanish, and satisfy
Z
Wint ðPt ; VÞ ¼
Pt
~ p dv ¼ 0: Sdis Mpen : D dis
~ p , we obtain the secSince this must be satisfied for all Pt and all symmetric and deviatoric tensors D dis ond microforce balance
ð4:34Þ
5. Free-energy imbalance We consider a purely mechanical theory based on a second law requiring that the temporal increase in free energy of any part Pt be less than or equal to the power expended on Pt . The second law therefore takes the form of a dissipation inequality
Z
_ wJ 1 dv 6 Wext ðPt Þ ¼ Wint ðPt Þ;
ð5:1Þ
Pt
with w the free energy, measured per unit volume in the structural space. Since J 1 dv with J ¼ det F represents the volume measure in the reference configuration, and since Pt deforms with the body,
Z Pt
Z _ wJ 1 dv ¼
_ 1 dv: wJ
Pt
Thus, since Pt is arbitrary, we may use (4.6), (4.22), and (4.23) to localize (5.1); the result is the local free-energy imbalance
1844
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
w_ JT : De Tp : Dp Sen : Dpen Sdis : Dpdis 6 0:
ð5:2Þ
We introduce a new stress measure def
Te ¼ JFe1 TFe> ;
ð5:3Þ e
which is symmetric, since T is symmetric; T represents a second Piola stress with respect to the intermediate structural space. Note that by using the definition (5.3), the Mandel stress defined in (4.28) is related to the stress measure Te by
Me ¼ Ce Te :
ð5:4Þ e
Next, differentiating (2.12)1 results in the following expression for the rate of change of C ,
C_ e ¼ Fe> F_ e þ F_ e> Fe : Hence, since Te is symmetric,
Te : C_ e ¼ 2Te : ðFe> F_ e Þ ¼ 2ðFe Te Fe> Þ : ðF_ e Fe1 Þ ¼ 2ðFe Te Fe> Þ : De ; or, using (5.3) we obtain
JT : De ¼
1 e _e T :C : 2
ð5:5Þ
Use of (5.5) in (5.2) allows us to express the free-energy imbalance in an alternate convenient form as
1 w_ Te : C_ e Tp : Dp Sen : Dpen Sdis : Dpdis 6 0: 2
ð5:6Þ
Finally, we note that w is invariant under a change in frame since it is a scalar field, and on account of the transformation rules discussed in Section 3, the transformation rules (4.14) and (4.15), and the definitions (4.28), (5.3), (4.33), and (5.4), the fields
Ce ;
Cpen ;
Bpen ;
Dp ;
Dpen ;
Dpdis ;
Te ;
Tp ;
Sen ;
Sdis ;
Me ;
and Mpen
ð5:7Þ
are also invariant. 6. Constitutive theory To account for the classical notion of isotropic strain hardening we introduce a positive-valued scalar internal state-variable S, which has dimensions of stress. We refer to S as the isotropic deformation resistance. Since S is a scalar field, it is invariant under a change in frame. Next, guided by the free-energy imbalance (5.6), and by experience with previous constitutive theories, we assume the following special set of constitutive equations:
9 ðeÞ ðCe Þ þ w ðpÞ ðBp Þ; > w¼w en > > > > e ðCe Þ; > Te ¼ T > > > p = Sen ¼ Sen ðB Þ; en
p ðDp ; SÞ; Tp ¼ T dis ðDp ; dp ; SÞ; Sdis ¼ S dis
p S_ ¼ hðd ; SÞ;
ð6:1Þ
> > > > > > > > > ;
where p def
d ¼ jDp j
ð6:2Þ p
is the scalar flow rate corresponding to D . Note that since w is the free energy per unit volume of the ðpÞ on Bp , because Ce and Bp are squared stretch ðeÞ is chosen to depend on Ce and w structural space, w en en like measures associated with the intermediate structural space. Also note that we have introduced a
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1845
p
possible dependence of Sdis on d because for some materials we anticipate that Dpdis (and thereby Sdis ) p may be constrained by the magnitude of d . ðeÞ represents the standard energy associated with intermolecular interactions, and the The energy w ðpÞ may be ðpÞ is a ‘‘defect-energy” associated with plastic deformation. For metallic materials w energy w attributed to an energy stored in the complex local microstructural defect state, which typically includes dislocations, sub-grain boundaries, and local incompatibilities at second-phase particles and ðpÞ leads to the developgrain boundaries. At the macroscopic continuum level, the ‘‘defect-energy” w ment and evolution of the energetic stress Sen , and this allows one to phenomenologically account for strain-hardening phenomena commonly called kinematic hardening. Isotropic hardening is modeled by the evolution (6.1)6 of the isotropic deformation resistance S. Finally, note that on account of the transformation rules listed in the paragraph containing (5.7) and since S is also invariant, the constitutive equations (6.1) are frame-indifferent. 6.1. Thermodynamic restrictions From (6.1)1
ðeÞ ðCe Þ ðpÞ ðBp Þ @w @w en : C_ e þ : B_ pen ; w_ ¼ e @Bpen @C and, using8
B_ pen ¼ Dpen Bpen þ Bpen Dpen ; and the symmetry of
Bpen
ð6:3Þ
p and @ w=@B en ,
ðpÞ Bp ðpÞ ðBp Þ p p @w @w en en : B_ pen ¼ : Den Ben þ Bpen Dpen p p @Ben @Ben ðpÞ p @ w ðBen Þ p B : Dpen : ¼2 en @Bpen 0 Hence, satisfaction of the free-energy imbalance (5.6) requires that the constitutive equations (6.1) satisfy
ðpÞ p ðeÞ ðCe Þ 1 e e @w en ðBp Þ 2 @ w ðBen Þ Bp T ðC Þ : C_ e þ S : Dpen e p en en 2 @Ben @C 0 dis ðDp ; dp ; SÞ : Dp P 0; p ðDp ; SÞ : Dp þ S þT dis
dis
ð6:4Þ
and hold for all arguments in the domains of the constitutive functions, and in all motions of the body. Thus, sufficient conditions that the constitutive equations satisfy the free-energy imbalance are that9 (i) the free-energy determine the stresses Te and Sen via the stress relations
ðeÞ e e ðCe Þ ¼ 2 @ w ðC Þ ; T @Ce ðpÞ p en ðBp Þ ¼ 2 @ w ðBen Þ Bp : S p en en @Ben 0
ð6:5Þ ð6:6Þ
(ii) the plastic distortion-rates Dp and Dpdis satisfy the dissipation inequality
Tp : Dp þ Sdis : Dpdis P 0:
ð6:7Þ
We assume henceforth that (6.5) and (6.6) hold in all motions of the body. We assume further that the material is strictly dissipative in the sense that Recall that Wpen ¼ 0, cf., (2.26). We content ourselves with constitutive equations that are only sufficient, but generally not necessary for compatibility with thermodynamics. 8 9
1846
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
Tp : Dp > 0 whenever Dp –0;
ð6:8Þ
Sdis : Dpdis > 0 whenever Dpdis –0:
ð6:9Þ
and
Some remarks: It is possible to enhance the theory by introducing another positive-valued dimen ðpÞ ðjÞ, sionless strain-like scalar internal variable j which contributes an additional defect energy w iso p which evolves as j evolves according to an evolution equation of the type j_ ¼ hðd ; jÞ. This will lead def ðpÞ to a scalar internal stress KðjÞ ¼ @ w iso ðjÞ=@ j, and the dissipation inequality (6.7) will be modified as
Tp : Dp þ Sdis : Dpdis K j_ P 0:
ð6:10Þ
ðpÞ ð w iso
Although a term such as jÞ might be needed for a more complete characterization of the ‘‘stored energy of cold work,” and a more precise determination of dissipation and thereby heat generation (cf., e.g., Rosakis et al., 2000; Ristinmaa et al., 2007), here, since our interest is in an isothermal theory, we refrain from introducing such additional complexity. However, we do wish to emphasize that our internal variable S, which represents an isotropic deformation resistance to plastic flow, is in no sense ðpÞ ðjÞ – it is introduced on a purely pragmatic pheintroduced as being derived from a defect energy w iso nomenological basis. Note that Ristinmaa et al. (2007) consider only ‘‘isotropic” hardening theories, while the theory under consideration here is for combined isotropic and kinematic hardening, and ðpÞ ðBp Þ already accounts for some aspects of the ‘‘stored energy of cold work.” the term w en Vladimirov et al. (2008) have considered a combined isotropic and kinematic-hardening theory with an additional defect energy wiso ðjÞ, and write their corresponding stress as @wiso ðjÞ=@ j ¼ ðRðjÞÞ as the strain-hardening contribution to the isotropic ðRðjÞÞ KðjÞ. They associate the stress pffiffiffiffiffiffiffiffi p deformation resistance, and take j_ ð 2=3Þd , so that j represents an accumulated plastic strain. It is our opinion that for materials which show strong isotropic strain hardening, the theory of Reese and co-workers will grossly underestimate the amount of dissipation. This is the major reason why we do not follow the approach of Reese and co-workers, but directly introduce our stress-like internal variable S to represent isotropic strain hardening. 6.2. Isotropy The following definitions help to make precise our notion of an isotropic material (cf., Anand and Gurtin, 2003): þ
(i) Orth = the group of all rotations (the proper-orthogonal group); (ii) the symmetry group GR , is the group of all rotations of the reference configuration that leaves the response of the material unaltered; (iii) the symmetry group GS at each time t, is the group of all rotations of the structural space that leaves the response of the material unaltered; (iv) the symmetry group GSS at each time t, is the group of all rotations of the substructural space that leaves the response of the material unaltered. Note that there is a slight difficulty in attaching a ‘‘physical meaning” to a rotation in GSS . Since the substructural space itself is a mathematical construct defined as the range of the linear map Fpdis , a rotation Q 2 GSS should be considered as a rotation in a ‘‘thought experiment” involving the substructural space; such a rotation may never be ‘‘physically attainable.” We now discuss the manner in which the basic fields transform under such transformations, granted the physically natural requirement of invariance of the internal power or equivalently, the requirement that
Te : C_ e ;
Tp : Dp ;
Sen : Dpen ;
and Sdis : Dpdis be invariant:
6.2.1. Isotropy of the reference space Let Q be a time-independent rotation of the reference configuration. Then F ! FQ , and hence
ð6:11Þ
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
Fp ! Fp Q ;
1847
and Fe and Fpen are invariant; and hence Ce and Cpen are invariant;
Fpdis ! Fpdis Q
ð6:12Þ so that, by (2.5), (2.12), and (2.21),
C_ e ; Dp ; Dpen ; and Dpdis are invariant: We may therefore use (6.11) to conclude that
Te ; Tp ; Sen ; and Sdis are invariant:
ð6:13Þ
Thus, with reference to the special constitutive functions (6.1), we note the these functions are not affected by rotations Q 2 GR . 6.2.2. Isotropy of the structural space Next, let Q , a time-independent rotation of the intermediate structural space, be a symmetry transformation. Then F is unaltered by such a rotation, and hence
Fe ! Fe Q ;
Fp ! Q > Fp ;
Fpen ! Q > Fpen ; and Fpdis is invariant;
ð6:14Þ
and also
Ce ! Q > Ce Q ;
C_ e ! Q > C_ e Q ;
Dpen ! Q > Dpen Q ;
Dp ! Q > Dp Q ;
Bpen ! Q > Bpen Q ;
and Dpdis is invariant:
ð6:15Þ
Then (6.15) and (6.11) yield the transformation laws
Te ! Q > Te Q ;
Tp ! Q > Tp Q ;
Sen ! Q > Sen Q ;
and Sdis is invariant:
ð6:16Þ
Thus, with reference to the constitutive equations (6.1) we conclude that
ðeÞ ðQ > Ce Q Þ; ðeÞ ðCe Þ ¼ w w p ðpÞ ðpÞ ðQ > Bp Q Þ; ðB Þ ¼ w w en en > e e e ðQ > Ce Q Þ; Q T ðC ÞQ ¼ T en ðQ > Bp Q Þ; en ðBp ÞQ ¼ S Q >S en
en
p ðDp ; SÞQ ¼ T p ðQ > Dp Q ; SÞ; Q >T
9 > > > > > > = > > > > > > ;
ð6:17Þ
must hold for all rotations Q in the symmetry group GS at each time t, while the constitutive functions dis ðDp ; dp ; SÞ, and hðdp ; SÞ are invariant to such rotations. S dis We refer to the material as one which is structurally isotropic if þ
GS ¼ Orth ;
ð6:18Þ
so that the response of the material is invariant under arbitrary rotations of the intermediate structural space at each time t. Henceforth, we assume that (6.18) holds. In this case, the response functions en , and T e , S p must each be isotropic. ðeÞ , w ðpÞ , T w 6.2.3. Isotropy of the substructural space Finally, let Q , a time-independent rotation of the substructural space, be a symmetry transformation. Then F, Fe , and Fp are unaltered by such a rotation, and hence
Fpen ! Fpen Q ;
Fpdis ! Q > Fpdis ;
ð6:19Þ
and also
Dpdis ! Q > Dpdis Q ;
and Dpen is invariant:
ð6:20Þ
Then (6.20) and (6.11) yield the transformation laws
Sdis ! Q > Sdis Q ;
and Sen is invariant:
ð6:21Þ
1848
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
Thus, with reference to the constitutive equations (6.1) we conclude that p
p
dis ðQ > Dp Q ; d ; SÞ dis ðDp ; d ; SÞQ ¼ S Q >S dis dis
ð6:22Þ
must hold for all rotations Q in the symmetry group GSS at each time t, while the other constitutive functions are invariant to such rotations. We refer to the material as one which is substructurally isotropic if þ
GSS ¼ Orth ;
ð6:23Þ
so that the response of the material is also invariant under arbitrary rotations of the substructural dis is an isotrospace at each time. Henceforth, we assume that (6.23) also holds at each time, so that S pic function. Thus, by assumption, all the tensorial response functions in (6.1) are isotropic. We confine our attention to materials which may be adequately defined by such isotropic functions, and refer to such materials as isotropic. 6.3. Consequences of isotropy of the free energy ðeÞ ðCe Þ is an isotropic function of Ce , it has the representation Since w
~ ðeÞ ðICe Þ; ðeÞ ðCe Þ ¼ w w
ð6:24Þ
where
ICe ¼ I1 ðCe Þ; I2 ðCe Þ; I3 ðCe Þ is the list of principal invariants of Ce . Thus, from (6.5),
~ ðeÞ e ðCe Þ ¼ 2 @ w ðICe Þ ; Te ¼ T @Ce
ð6:25Þ
e ðCe Þ is an isotropic function of Ce . Then since the Mandel stress is defined by (cf., (5.4)) and T
Me ¼ Ce Te ; e ðCe Þ is isotropic, we find that Te and Ce commute, and T
Ce Te ¼ Te Ce ;
ð6:26Þ e
and hence that the elastic Mandel stress M is symmetric. Further the defect free energy has a representation
~ ðpÞ ðI p Þ where I p ¼ I1 ðBp Þ; I2 ðBp Þ; I3 ðBp Þ ; ðpÞ ¼ w w Ben Ben en en en
ð6:27Þ
and this yields that
~ ðpÞ ðI p Þ @w Ben @Bpen
Bpen
ð6:28Þ
is a symmetric tensor, and (6.6) reduces to
en ðBp Þ ¼ 2 Sen ¼ S en
~ ðpÞ ðI p Þ @w Ben @Bpen
! Bpen
; 0
a symmetric and deviatoric tensor. Also the plastic Mandel stress defined in (4.33) is
Mpen
¼
2Fp> en
~ ðpÞ ðI p Þ @w Ben @Bpen
! Bpen
Fp> en ; 0
ð6:29Þ
1849
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
" ¼ 2 Fp> en " ¼ 2 Fp> en
~ ðpÞ ðI p Þ @w Ben @Bpen ~ ðpÞ ðI p Þ @w Ben @Bpen
!
#
Bpen Fp> en !
# Fpen
; 0
;
ð6:30Þ
0
~ ðpÞ ðI p Þ=@Bp is symmetric, this stress is symmetric and deviatoric. and since @ w Ben en 7. Flow rules Henceforth, in accord with common terminology we call the symmetric and deviatoric tensor Sen a back-stress, and we denote an effective deviatoric Mandel stress by
Meeff
0
def
¼ Me0 Sen :
ð7:1Þ
Then, upon using the constitutive relation (6.1)4 and the first microforce balance (4.29), together with symmetry of the Mandel stress discussed above, a central result of the theory is the first flow rule
Meeff
0
p ðDp ; SÞ: ¼T
ð7:2Þ
Next, on account of the symmetry and deviatoric nature of the plastic stress Mpen (cf., (6.30)), the second microforce balance (4.34) and the constitutive equation (6.4)5 yield the second flow rule of the theory: p
dis ðDp ; d ; SÞ: Mpen ¼ S dis
ð7:3Þ
We now make two major assumptions concerning the plastic flow for isotropic materials: (1) Codirectionality hypotheses: Recall p
d ¼ jDp j;
p
and let ddis ¼ jDpdis j;
also let def
Np ¼
Dp p d
def
and Npdis ¼
Dpdis p ddis
ð7:4Þ
denote the directions of plastic flow whenever Dp –0 and Dpdis –0. Then the dissipation inequality (6.7) may be written as
dis ðdp ; Np ; dp ; SÞ : Np dp P 0: p ðdp ; Np ; SÞ : Np dp þ S T dis dis dis dis
ð7:5Þ
Guided by (7.5), we assume henceforth that the dissipative flow stress Tp is parallel to and points in the same direction as Np so that
p ðdp ; Np ; SÞ ¼ Yðdp ; Np ; SÞNp ; T
ð7:6Þ
where p p ðdp ; Np ; SÞ : Np Yðd ; Np ; SÞ ¼ T
ð7:7Þ
represents a scalar flow strength of the material. Similarly we assume that dissipative flow stress Sdis is parallel to and points in the same direction as Npdis so that
dis ðdp ; Np ; dp ; SÞ ¼ Y dis ðdp ; Np ; dp ; SÞNp ; S dis dis dis dis dis
ð7:8Þ
where p p dis ðdp ; Np ; dp ; SÞ : Np Y dis ðddis ; Npdis ; d ; SÞ ¼ S dis dis dis
ð7:9Þ
represents another scalar flow strength of the material. We refer to the assumptions (7.6) and (7.8) as the codirectionality hypotheses.10 10
These assumptions corresponds to the classical notion of maximal dissipation in Mises-type theories of metal plasticity.
1850
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
(2) Strong isotropy hypotheses: p We also assume that the scalar flow strength Yðd ; Np ; SÞ is independent of the flow direction Np , so that p
Yðd ; SÞ:
ð7:10Þ
Similarly, we also assume that the scalar flow strength flow direction Npdis , so that p
p p Y dis ðddis ; Npdis ; d ; SÞ
is independent of the
p
Y dis ðddis ; d ; SÞ:
ð7:11Þ
We refer to the assumptions (7.10) and (7.11) as the strong isotropy hypotheses. Thus, using (7.10) and (7.6), the flow rule (7.2) reduces to, p
ðMeeff Þ0 ¼ Yðd ; SÞNp ;
ð7:12Þ
which immediately gives
Np ¼
ðMeeff Þ0 ; jðMeeff Þ0 j
ð7:13Þ
and p
jðMeeff Þ0 j ¼ Yðd ; SÞ:
ð7:14Þ
jðMeeff Þ0 j
p
When and S are known, (7.14) serves as an implicit equation for the scalar flow rate d . Next, using (7.11) and (7.8), the flow rule (7.3) reduces to, p
p
Mpen ¼ Y dis ðddis ; d ; SÞNpdis ;
ð7:15Þ
which gives
Npdis ¼
Mpen ; jMpen j
ð7:16Þ
and p
p
jMpen j ¼ Y dis ðddis ; d ; SÞ:
ð7:17Þ
jðMeeff Þ0 j
p
When and S are known, (7.17) serves as an implicit equation for the scalar flow rate ddis . Finally, using (7.6), (7.10), and (7.14), as well as (7.8), (7.11), and (7.17) the dissipation inequality (6.7) reduces to p
p
jðMeeff Þ0 jd þ jMpen jddis P 0:
ð7:18Þ
8. Summary of the constitutive theory In this section we summarize our constitutive theory. (1) Free energy.
~ ðeÞ ðICe Þ þ w ~ ðpÞ ðI p Þ; w¼w Ben
ð8:1Þ e
where ICe and IBpen are the lists of the principal invariants of C and (2) Cauchy stress.
~ ðeÞ ðICe Þ @w def T ¼ J 1 Fe Te Fe> where Te ¼ 2 : @Ce
Bpen ,
respectively.
ð8:2Þ
(3) Mandel stresses. Back-stress. The elastic Mandel is given by
Me ¼ Ce Te ;
ð8:3Þ
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1851
and is symmetric. The back-stress is given by
~ ðpÞ ðI p Þ @w Ben
Sen ¼ 2
@Bpen
! Bpen
ð8:4Þ
; 0
and is symmetric deviatoric. The stress difference
Meeff ¼ Me Sen
ð8:5Þ
is called an effective elastic Mandel stress. The plastic Mandel stress is given by
"
! ~ ðpÞ ðI p Þ @w Ben
Mpen ¼ 2 Fp> en
@Bpen
# Fpen
ð8:6Þ
; 0
and is symmetric and deviatoric. (4) Flow rules. The evolution equation for Fp is
F_ p ¼ Dp Fp ;
ð8:7Þ
p
with D given by p
Dp ¼ d Np ;
Np ¼
ðMeeff Þ0 ; jðMeeff Þ0 j
ð8:8Þ p
where the scalar flow rate d is obtained by solving the scalar strength relation p
jðMeeff Þ0 j ¼ Yðd ; SÞ; for given for Fpdis is
ðMeeff Þ0
ð8:9Þ p
and S, where Yðd ; SÞ is a rate-dependent flow strength. The evolution equation
F_ pdis ¼ Dpdis Fpdis ; with
Dpdis
ð8:10Þ
given by
p
Dpdis ¼ ddis Npdis ;
Npdis ¼
Mpen ; jMpen j
ð8:11Þ
p
where the scalar flow rate ddis is obtained by solving the scalar strength relation p
p
jMpen j ¼ Y dis ðddis ; d ; SÞ; Mpen ,
ð8:12Þ
p
for given d , and S, where (5) Evolution equation for S.
p p Y dis ðddis ; d ; SÞ
is another rate-dependent flow strength.
p S_ ¼ hðd ; SÞ:
ð8:13Þ
The evolution equations for Fp , Fpdis , and S need to be accompanied by initial conditions. Typical initial conditions presume that the body is initially (at time t ¼ 0, say) in a virgin state in the sense that
FðX; 0Þ ¼ Fp ðX; 0Þ ¼ Fpdis ðX; 0Þ ¼ 1; e p
p
so that by F ¼ F F and F ¼
Fpen Fpdis
SðX; 0Þ ¼ S0 ð¼ constantÞ; e
we also have F ðX; 0Þ ¼ 1 and
Fpen ðX; 0Þ
ð8:14Þ ¼ 1.
9. Specialization of the constitutive equations In this section, based on experience with existing recent theories of viscoplasticity with isotropic and kinematic hardening, we specialize our constitutive theory by imposing additional constitutive assumptions.
1852
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
9.1. Free-energy wðeÞ The spectral representation of Ce is
Ce ¼
3 X
xei rei rei ; with xei ¼ ke2 i ;
ð9:1Þ
i¼1
where ðre1 ; re2 ; re3 Þ are the orthonormal eigenvectors of Ce and Ue , and ðke1 ; ke2 ; ke3 Þ are the positive eigenvalues of Ue . Instead of using the invariants ICe , the free-energy wðeÞ for isotropic materials may be alternatively expressed in terms of the principal stretches as
ðeÞ ðke ; ke ; ke Þ: wðeÞ ¼ w 1 2 3
ð9:2Þ
Then, by the chain-rule and (8.2)2 , the stress Te is given by
Te ¼ 2
3 3 ðeÞ ðke ; ke ; ke Þ ðeÞ ðke ; ke ; ke Þ @ke X ðeÞ ðke ; ke ; ke Þ @ xi X @w @w 1 @w 1 2 3 1 2 3 i 1 2 3 ¼2 : e e e ¼ e @ki ki @kei @C @Ce @C i¼1 i¼1
ð9:3Þ
Assume that the squared principal stretches xei are distinct, so that the xei and the principal directions rei may be considered as functions of Ce ; then
@ xei ¼ rei rei ; @Ce
ð9:4Þ
and, granted this, (9.4) and (9.3) imply that
Te ¼
3 ðeÞ ðke ; ke ; ke Þ X 1 @w 1 2 3 rei rei : e k @kei i i¼1
ð9:5Þ
Further, from (8.2)1 , 1 e e e>
T¼J F T F
1
e
e e
e
¼J R U T U R
e>
1
¼J R
e
3 X i¼1
kei
! ðeÞ ðke ; ke ; ke Þ @w e e 1 2 3 ri ri Re> : @kei
ð9:6Þ
Next, since Me ¼ Ce Te (cf., (8.3)), use of (9.1) and (9.5) gives the Mandel stress as
Me ¼
3 X
kei
i¼1
ðeÞ ðke ; ke ; ke Þ @w 1 2 3 rei rei : @kei
ð9:7Þ
Let def
Ee ¼ ln Ue ¼
3 X
Eei rei rei ;
ð9:8Þ
i¼1
denote the logarithmic elastic strain with principal values def
Eei ¼ ln kei ;
ð9:9Þ
and consider an elastic free-energy function of the form
ðeÞ ðEe ; Ee ; Ee Þ; ðeÞ ðke ; ke ; ke Þ ¼ w w 1 2 3 1 2 3
ð9:10Þ
so that, using (9.7),
Me ¼
3 ðeÞ ðEe ; Ee ; Ee Þ X @w 1 2 3 rei rei : @Eei i¼1
ð9:11Þ
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1853
We consider the following simple generalization of the classical strain-energy function of infinitesimal isotropic elasticity which uses a logarithmic measure of finite strain11
ðeÞ ðEe Þ ¼ GjEe j2 þ 1 KðtrEe Þ2 ; w 0 2
ð9:12Þ
G > 0 and K > 0;
ð9:13Þ
where
are the shear modulus and bulk modulus, respectively. Then, (9.11) gives
Me ¼ 2GEe0 þ KðtrEe Þ1
ð9:14Þ
and on account of (9.6), (9.7), and (9.14),
T ¼ J 1 Re Me Re> :
ð9:15Þ
9.2. Free-energy wðpÞ The spectral representation of Bpen is
Bpen ¼
3 X
bi l i l i ;
with bi ¼ a2i ;
ð9:16Þ
i¼1
where ðl1 ; l2 ; l3 Þ are the orthonormal eigenvectors of Bpen and Vpen , and ða1 ; a2 ; a3 Þ are the positive eigenvalues Vpen . Instead of using the invariants IBpen , the free-energy wðpÞ may alternatively be expressed as
ðpÞ ða1 ; a2 ; a3 Þ: wðpÞ ¼ w
ð9:17Þ
Then, by the chain-rule 3 3 ðpÞ ða1 ; a2 ; a3 Þ X ðpÞ ða1 ; a2 ; a3 Þ @ai ðpÞ ða1 ; a2 ; a3 Þ @bi @w @w 1X 1 @w ¼ : p p ¼ @a 2 a @ai @B @Bpen @Ben i i en i¼1 i¼1
ð9:18Þ
Assume that bi are distinct, so that the bi and the principal directions li may be considered as functions of Bpen . Then,
@bi ¼ li li ; @Bpen
ð9:19Þ
and, granted this, (9.18) implies that 3 ðpÞ ða1 ; a2 ; a3 Þ 1 X ðpÞ ða1 ; a2 ; a3 Þ @w 1 @w ¼ li li : p 2 a @ai @Ben i i¼1
ð9:20Þ
Also, use of (9.16) and (9.20) in (8.4) gives the deviatoric back-stress as
Sen ¼
3 X i¼1
ðpÞ ða1 ; a2 ; a3 Þ @w ai li li @ai
! :
ð9:21Þ
0
Let def
Epen ¼ ln Vpen ¼
3 X
ln ai li li :
ð9:22Þ
i¼1
11 Vladimirov et al. (2008) employ a Neo-Hookean free-energy function modified for elastic compressibility (cf., their Eq. (28)). Here, as an alternate, we use the classical strain-energy function for infinitesimal elasticity, but employ the large deformation Hencky-logarithmic elastic strain Ee . Anand (1979, 1986) has shown that the simple free-energy function (9.12) is useful for moderately large elastic stretches. Additionally, use of the logarithmic strain helps in developing implicit time-integration procedures based on the exponential map (cf., Appendix A).
1854
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
denote a logarithmic energetic strain with principal values ðln ai Þ. Note that since trVpen ¼ a1 a2 a3 ¼ 1,
trðEpen Þ ¼ ln a1 þ ln a2 þ ln a3 ¼ lnða1 a2 a3 Þ ¼ 0; and hence the strain tensor
Epen
ð9:23Þ
is traceless, Next, consider the following simple defect energy:
h i ðpÞ ðEp Þ ¼ BjEp j2 ¼ B ðln a1 Þ2 þ ðln a2 Þ2 þ ðln a3 Þ2 ; w en en
ð9:24Þ
with B > 0, which, since Epen is deviatoric, is analogous to (9.12). Then using (9.21), (9.22), and (9.23)
Sen ¼ 2BEpen ;
ð9:25Þ
we call the positive-valued constitutive parameter B the back-stress modulus. Then the plastic Mandel stress (4.33) becomes
p> p> p p p1 Mpen ¼ Fp> Rpen ; en Sen Fen ¼ R en V en ð2BEen ÞV en or p Mpen ¼ Rp> en Sen R en :
ð9:26Þ
9.3. Strength functions First consider the strength relation p
jðMeeff Þ0 j ¼ Yðd ; SÞ;
ð9:27Þ
appearing in (8.9). Define an equivalent shear stress by
1
s def ¼ pffiffiffi jðMeeff Þ0 j;
ð9:28Þ
2
and an equivalent shear strain rate by
pffiffiffi
pffiffiffi
p mp def ¼ 2d ¼ 2jDp j;
ð9:29Þ
respectively. Using the definitions (9.28) and (9.29), we rewrite the strength relation (9.27) as12
s ¼ Yðmp ; SÞ;
ð9:30Þ
and assume that the shear flow strength Yðmp ; SÞ has the simple form
Yðmp ; SÞ ¼ gðmp ÞS;
ð9:31Þ
where S now represents an isotropic flow resistance in shear, and
gðmp Þ is a strictly increasing function of
gð0Þ ¼ 0;
mp :
ð9:32Þ
p
We refer to the dimensionless function gðm Þ as the rate-sensitivity function. Next, by (9.32) the function gðmp Þ is invertible, and the inverse function f ¼ g 1 is strictly increasing and hence strictly positive for all nonzero arguments; further f ð0Þ ¼ 0. Hence the relation (9.30) may be inverted to give an expression
mp ¼ g 1
s s f S S
ð9:33Þ
for the equivalent plastic shear strain rate. An example of a commonly used rate-sensitivity function is the power-law function
gðmp Þ ¼
12
p m
m m0
ð9:34Þ
;
We absorb inconsequential factors of
pffiffiffi 2 in writing (9.30).
1855
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
where m > 0, a constant, is a rate-sensitivity parameter and m0 > 0, also a constant, is a reference flowrate. The power-law function allows one to characterize nearly rate-independent behavior, for which m is very small. Further, granted the power-law function (9.34), the expression (9.33) has the specific form
mp ¼ m0
m1 s : S
ð9:35Þ
Next consider the strength relation p
p
jMpen j ¼ Y dis ðddis ; d ; SÞ;
ð9:36Þ
appearing in (8.12). Define another equivalent shear stress by
1
sen def ¼ pffiffiffi jMpen j;
ð9:37Þ
2
and another equivalent shear strain rate by
pffiffiffi
pffiffiffi
p mpdis def ¼ 2ddis ¼ 2jDpdis j;
ð9:38Þ
respectively. Using the definitions (9.37) and (9.38), we rewrite the strength relation (9.36) as
sen ¼ Y dis ðmpdis ; mp ; SÞ: A specialization of Y dis ðm ening is
Y dis ðmpdis ; mp ; SÞ ¼
ð9:39Þ
p dis ;
p
m ; SÞ which leads to Armstrong and Frederick (1966) type kinematic hard
B
cm
13
p
mpdis ;
ð9:40Þ
where B is the back-stress modulus (cf., (9.25)), and c P 0 is a dimensionless constant. Using (9.40), (9.39) becomes
sen ¼
B
mpdis ;
cmp
ð9:41Þ
so that the term ðB=cmp Þ represents a ‘‘pseudo”-viscosity (Dettmer and Reese, 2004). The relation (9.41) may be inverted to give
mpdis ¼
sen
B
cmp ;
ð9:42Þ
which shows that mpdis ¼ 0 when either c ¼ 0 or when mp ¼ 0. Thus note that the constitutive assumption (9.42) constrains the theory such that Dpdis ¼ 0 whenever Dp ¼ 0, and this in turn will lead to no change in the back-stress Sen when Dp ¼ 0, that is when there is no macroscopic plastic flow. 9.4. Evolution equation for S The evolution equation for the isotropic deformation resistance is specialized in a rate-independent form as
S_ ¼ hðSÞmp
ð9:43Þ
with hðSÞ a hardening function, which we take to be given by (Brown et al., 1989):
( hðSÞ ¼
13
a h0 1 SS
for S0 6 S 6 S ;
0
for S P S ;
Again absorbing inconsequential factors of
pffiffiffi 2 in writing (9.39).
ð9:44Þ
1856
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
where S , a, and h0 are constant moduli with S > S0 , a P 1, and h0 > 0. The hardening function (9.44) is strictly decreasing for S0 6 S 6 S and vanishes for S P S . 9.5. Dissipation inequality Finally, using (9.28), (9.29), (9.37), (9.38), and (9.42), the dissipation inequality (7.18) reduces to
c s þ ðsen Þ2 mp P 0:
ð9:45Þ
B
9.6. Summary of the specialized constitutive model In this section, we summarize the specialized form of our theory, which should be useful in applications. The theory relates the following basic fields: x ¼ vðX; tÞ, F ¼ rv; J ¼ det F > 0, F ¼ Fe Fp , Fe ; J e ¼ det Fe ¼ J > 0, Fp ; J p ¼ det Fp ¼ 1, Fe ¼ Re Ue , Dp ¼ symðF_ p Fp1 Þ, Fp ¼ Fpen Fpdis , Fpen ; det Fpen ¼ 1, Fpdis ; det Fpdis ¼ 1, Fpen ¼ Vpen Rpen , Þ, Dpdis ¼ symðF_ pdis Fp1 dis T ¼ T> , w, S > 0;
motion; deformation gradient; elastic–plastic decomposition of F; elastic distortion; inelastic distortion; polar decomposition of Fe ; plastic stretching corresponding to Fp ; elastic–plastic decomposition of Fp ; energetic part of Fp ; dissipative part of Fp ; polar decomposition of Fpen ; plastic stretching corresponding to Fpdis ; Cauchy stress; free-energy density per unit volume of intermediate structural space; scalar internal variable.
(1) Free energy. We consider a separable free energy
w ¼ wðeÞ þ wðpÞ : With
Ue ¼
3 X
kei rei rei ;
ð9:46Þ
ð9:47Þ
i¼1
denoting the spectral representation of Ue , and with
Ee ¼
3 X
ln kei rei rei ;
ð9:48Þ
i¼1
denoting an elastic logarithmic strain measure, we adopt the following special form for the freeenergy wðeÞ :
1 wðeÞ ¼ GjEe0 j2 þ KðtrEe Þ2 ; 2
ð9:49Þ
where
G > 0;
and K > 0;
ð9:50Þ
are the shear modulus and bulk modulus, respectively. Further, with
Vpen ¼
3 X i¼1
ai li li ;
ð9:51Þ
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1857
denoting the spectral representation of Vpen , and with
Epen ¼
3 X
ln ai li li ;
trEpen ¼ 0;
ð9:52Þ
i¼1
we adopt a free-energy wðpÞ of the form
wðpÞ ¼ BjEpen j2 ;
ð9:53Þ
where the positive-valued parameter
B P 0;
ð9:54Þ
is a back-stress modulus. (2) Cauchy stress. Mandel stresses. Back-stress. Effective stress. Corresponding to the special free-energy functions considered above, the Cauchy stress is given by
T ¼ J 1 Re Me Re> ;
ð9:55Þ
where
Me ¼ 2GEe0 þ KðtrEe Þ1;
ð9:56Þ
is an elastic Mandel stress. The symmetric and deviatoric back-stress is defined by
Sen ¼ 2BEpen ;
ð9:57Þ p
and the driving stress for D is the effective stress given by
ðMeeff Þ0 ¼ Me0 Sen :
ð9:58Þ
The corresponding equivalent shear stress is given by
1
s ¼ pffiffiffi jðMeeff Þ0 j: 2
ð9:59Þ
The driving stress for Dpdis is the stress measure p Mpen ¼ Rp> en Sen R en :
ð9:60Þ
The corresponding equivalent shear stress is given by
1
sen ¼ pffiffiffi jMpen j: 2
ð9:61Þ
(3) Flow rules. The evolution equation for Fp is
9 F_ p ¼ Dp Fp ; Fp ðX; 0Þ ¼ 1; > > > e = ðM Þ Dp ¼ mp 2effs 0 ; > > > ; mp ¼ m0 sS 1=m ;
ð9:62Þ
where m0 is a reference plastic strain rate with units of 1=time, and m is a strain-rate-sensitivity parameter. The evolution equation for Fpdis is
9 F_ pdis ¼ Dpdis Fpdis ; Fpdis ðX; 0Þ ¼ 1; > > = p M Dpdis ¼ mpdis 2senen ; > > ; mpdis ¼ c sBen mp ;
where c P 0 is a dimensionless constant.
ð9:63Þ
1858
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
(4) Evolution equation for the isotropic-hardening variable S. The internal variable S is taken to obey the evolution equation
S_ ¼ hðSÞmp ; SðX; 0Þ ¼ S0 ; ( a for S0 6 S 6 S ; h0 1 SS hðSÞ ¼ 0 for S P S ;
9 > > = > > ;
ð9:64Þ
where S , a, and h0 are constant moduli with S > S0 , a P 1, and h0 > 0.
10. Some numerical solutions In Appendix A, we develop a semi-implicit time-integration procedure, and obtain an (approximate) algorithmic tangent for the special constitutive theory summarized in the previous section. The time-integration procedure and the associated algorithmic tangent have been implemented in the commercial finite element program Abaqus/Standard (2008) by writing a user material subroutine (UMAT). In this section, we investigate some salient properties of the model and the accompanying time-integration procedure. 10.1. Simple tension/compression First, in order to demonstrate that the constitutive theory captures the various isotropic and kinematic-hardening effects, we carried out a numerical simulation of a symmetric strain-cycle in simple tension and compression. The material parameters used in our numerical study are:
G ¼ 80 GPa; m0 ¼ 0:001 s1 ; S0 ¼ 100 MPa; B ¼ 40 GPa;
K ¼ 175 GPa; m ¼ 0:02; h0 ¼ 1250 MPa; a ¼ 1 S ¼ 250 MPa; c ¼ 400;
ð10:1Þ
the value of m ¼ 0:02 is intended to represent a nearly rate-independent material. The axial strain-cycle was imposed between true strain limits of ¼ 0:02 at an absolute axial true strain rate of j_ j ¼ 0:01 s1 ; a fixed time step Dt ¼ 0:01 s was used in the numerical simulations. The calculations were performed for: (i) no hardening, h0 ¼ 0 and B ¼ 0; (ii) isotropic hardening only, h0 –0 and B ¼ 0; (iii) kinematic hardening only, h0 ¼ 0 and B–0; and finally (iv) combined kinematic and isotropic hardening, h0 –0 and B–0. Plots of the axial stress r versus the axial strain are shown in Fig. 1(a). With no hardening, we obtain an elastic-perfectly-plastic response, as expected. When isotropic hardening only is allowed, we observe an increase in the flow stress, and the stress level at which plastic flow recommences on strain reversal is equal in magnitude to the stress level from which the reversal in strain was initiated. In the case of kinematic hardening only, we again observe an increase in the flow stress, but this time the magnitude of the stress level at which plastic flow recommences on strain reversal is substantially smaller than the stress level from which the reversal in strain direction was initiated; a clear manifestation of the Bauschinger effect. When isotropic hardening is combined with kinematic hardening, we see the continual evolution in the flow stress associated with combined isotropic and kinematic hardening during forward straining, and a clear Bauschinger effect due to the kinematic hardening upon strain reversal. 10.2. Simple shear Next, we concentrate on numerical solutions related to homogeneous simple shear. With respect to a rectangular Cartesian coordinate system with origin o and orthonormal base vectors fei ji ¼ 1; 3g, a simple shearing motion is described by
x ¼ X þ ðmtÞX 2 e1 ;
ð10:2Þ
with m a shear strain rate, and C ¼ mt is the amount of shear. Throughout this section we use the material parameters summarized in (10.1).
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
a
1859
600
Axial stress, σ (MPa)
400 200 0
No hardening Isotropic only Kinematic only Combined
−200 −400 −600 −0.02
b
−0.01
0
Axial strain, ε
0.01
0.02
300
Stress, T
12
(MPa)
200 100 No hardening Isotropic only Kinematic only Combined
0 −100 −200 −300 0
0.01
0.02
0.03
Amount of shear, Γ
0.04
0.05
Fig. 1. (a) Comparison of axial stress r versus axial strain for various types hardening in a symmetric strain-cycle simulation in simple tension and compression. (b) Comparison of shear stress T 12 versus shear strain C for various types of hardening in reversed simple shear.
The specific goals of our numerical solutions related to simple shear are as follows: To further demonstrate that the constitutive theory captures the major features of isotropic and kinematic hardening in the case of simple reversed shear. To investigate the convergence and stability properties of the time-integration scheme in simple reversed shear. To demonstrate that the constitutive equations do not suffer from oscillations in the stress response at large shear strains. To demonstrate that the constitutive equations and the time-integration procedure are indeed objective.
1860
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
First, we wish to further demonstrate that the constitutive theory captures the various isotropic and kinematic-hardening effects. To this end, we consider simple shear, as described above, for a time of 5 s at a shear strain rate of m ¼ 0:01 s1 for a total shear strain C ¼ mt ¼ 0:05, and then at a shear strain rate of m ¼ 0:01 s1 for another 5 s to complete one reversal of strain, while using a fixed time step Dt ¼ 0:01 s. As for simple tension/compression, the calculations were performed for: (i) no hardening, h0 ¼ 0 and B ¼ 0; (ii) isotropic hardening only, h0 –0 and B ¼ 0; (iii) kinematic hardening only, h0 ¼ 0 and B–0; and finally (iv) combined kinematic and isotropic hardening, h0 –0 and B–0, and the results are shown in Fig. 1(b). Again, we observe the important qualitative features of each case in simple reversed shear, i.e. elastic-perfectly-plastic response in the case of no hardening, a continual increase in the flow stress in cases involving isotropic hardening, and a clear Bauschinger effect in cases involving kinematic hardening. Next, we investigated the stability and convergence properties of our time-integration procedure. Since our integration procedure is semi-implicit rather than fully-implicit, one might be concerned that there may be a restrictive time step constraint for stability. To ascertain any time step restrictions, we performed a convergence test. Reversed simple shear to a maximum shear strain of C ¼ 0:05 at a shear strain rate of jmj ¼ 0:01 s1 was considered while increasingly coarsening the time step. The specific time steps Dt of
1 103 ;
1 102 ;
2:5 102 ;
1 101 ;
2:5 101 ;
and 1 s
were considered; these time steps translate to shear strain increments of
1 105 ;
1 104 ;
2:5 104 ;
1 103 ;
2:5 103 ;
and 1 102 ;
respectively. Fig. 2 shows the stress–strain curves corresponding to these cases. For increasing time steps, no evidence of instability is observed, but of course the accuracy of the solution degrades, especially in the regions of elastic–plastic transitions. The order of convergence is demonstrated in Fig. 3; since no exact solution is available, the solution obtained using a time step of Dt ¼ 1 103 s was treated as the ‘‘exact” solution and used to calculate the pointwise error for the cases corresponding to the larger time steps. The error E was then obtained by taking the L-infinity norm of the pointwise error. As seen in Fig. 3, the time-integration procedure is first-order accurate. Thus, while no stability problems are encountered, the accuracy of the integration procedures deteriorates, especially when large time increments are taken during periods of rapidly changing elastic–plastic transitions.
300
(MPa)
200
Stress, T
12
100 Δ t = 1×10−3 s −2
s
−1
s
0
Δ t = 1×10
−100
Δ t = 1×10 Δt=1s
−200 −300
0
0.01
0.02
0.03
Amount of shear, Γ
0.04
0.05
Fig. 2. Comparison of solutions for reversed simple shear for various time steps.
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1861
8.5
8
1
log(E)
1 7.5
7
6.5
6 −2.5
−2
−1.5
−1
log(Δ t)
−0.5
0
0.5
Fig. 3. Convergence diagram for the time-integration procedures demonstrating first-order accuracy.
We wish next to demonstrate that the constitutive equations do not produce shear-stress oscillations at large strains – a feature which has plagued numerous previous large deformation constitutive theories for kinematic hardening. To this end, we performed a monotonic simple shear simulation to a very large final shear strain of C ¼ 10, at a shear strain rate of m ¼ 0:01 s1 using a fixed time step Dt ¼ 0:1 s. The shear T 12 and normal components, T 11 and T 22 , of the Cauchy stress T are plotted in Fig. 4(a). Clearly there are no oscillations in the T 12 versus C response; the shear stress rises and saturates (because of the chosen value of the material parameters) at a shear strain of C 1:5. We have also carried out a similar calculation with all other material parameters the same as before, but by setting the dynamic recovery parameter c ¼ 0; this corresponds to the case of ‘‘linear” kinematic hardening. The variation of the shear and normal components of the Cauchy stress T for this case are plotted in Fig. 4(b). Here, since there is no dynamic recovery for the kinematic hardening, the stress levels for this physically unrealistic case get quite large, but the T 12 versus C curve still does not show any oscillations; the peak in this curve is due to the use of a logarithmic energetic strain lnðVpen Þ and not due to any oscillatory behavior in the evolution of the back-stress. Finally, in order to demonstrate the frame-indifference of the constitutive equations and numerical algorithm, we have carried out a numerical calculation for simple shear with superposed rigid rotation, which is described by Weber et al. (1990):
x ¼ Q ðtÞ½X þ ðmtÞX 2 e1 ;
ð10:3Þ
Q ðtÞ ¼ ðe1 e1 þ e2 e2 Þ cosðxtÞ þ ðe1 e2 e2 e1 Þ sinðxtÞ þ e3 e3 ;
ð10:4Þ
with
a rotation about the e3 -axis. In performing our numerical calculation, we used the material parameters listed in (10.1), and the calculation was performed with x ¼ 0:1p radians per second and m ¼ 0:01 per second, for t 2 ½0; 20 , so that h xt 2 ½0; 2p and C mt 2 ½0; 0:2 . The calculation was carried by using a fixed time step of Dt ¼ 0:1 s, which corresponds to a shear strain increment of DC ¼ 0:001 and a rotation increment of Dh ¼ 3:6 . Snapshots of the initial and deformed geometry at a few representative stages of the simple shear plus superposed rotation are shown in Fig. 5. Let the solution for the stress corresponding to the motion (10.3) be denoted by T ðtÞ. The result of the shear component of the ‘‘unrotated” Cauchy stress ðQ ðtÞ> T ðtÞQ ðtÞÞ versus the amount of shear C, obtained from the numerical calculation is shown in Fig. 6. Also shown in this figure is T 12 versus C for the motion (10.3) with Q ðtÞ ¼ 1. The two shear stress versus shear strain curves are indistinguishable from each other; the
1862
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
a
400 350
Stress (MPa)
300 250 200
T12
150
T
T
11 22
100 50 0 −50 0
b
2
4
6
8
10
4
6
8
10
Amount of shear, Γ
5
1.5
x 10
T
12
1
T11 T
Stress (MPa)
22
0.5 0 −0.5 −1 −1.5 0
2
Amount of shear, Γ
Fig. 4. Simple shear to large shear strain (a) with dynamic recovery c–0; and (b) without dynamic recovery, c = 0. Note that there are no shear-stress oscillations in either case.
excellent agreement between the two calculations verifies the objectivity of the constitutive equations and time-integration procedure. 10.3. Cyclic loading of a curved bar Next, in order to exercise the algorithmic tangent in a case of inhomogeneous deformation, we consider the cyclic deformation of a curved bar of circular cross-section. A finite element mesh for the curved bar is shown in Fig. 7(a). The bar, which has a diameter of 50 mm, curves through a radius of 100 mm for 90 , and then extends straight for another 100 mm. The mesh consists of 2184 Abaqus-C3D8R elements. The nodes on the face denoted as ‘‘fixed”-face in Fig. 7(a) have prescribed null displacements, u1 ¼ u2 ¼ u3 ¼ 0, while the boundary conditions on the face denoted as the ‘‘moved”-face were specified as
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1863
Fig. 5. Simple shear with superposed rigid rotation. The deformed geometry at each stage is shown with a solid line, while the dashed line represents the initial geometry.
u1 ¼ u3 ¼ 0;
and u2 ðtÞ ¼ a sinðxtÞ;
with a ¼ 20 mm and x ¼ 0:2p s1 :
ð10:5Þ
The calculation was run to a final time of t f ¼ 30 s, so that the u2 -displacement on the ‘‘moved”-face goes through three complete cycles. The material properties used in the calculation were the same as those in the previous calculations (10.1). The deformed mesh at the point of maximum deflection is also shown in Fig. 7(a). In this (and subsequent) multi-element calculations, a variable time-stepping integration procedure was employed. As in the paper by Lush et al. (1989), the calculation was performed with a slightly modified version of the Abaqus/Standard static analysis procedure. The modification was introduced to enhance the automatic time-stepping procedure in Abaqus to control the accuracy of the constitutive time-integration. This was done by using as a control measure the maximum equivalent plastic shear strain increment
Dcpmax ¼ ðDtÞðmpnþ1 Þmax occurring at any integration point in the model during an increment. The value of Dcpmax was kept close to a specified nominal value Dcps . From Fig. 2, we note that the accuracy of the solution begins to
1864
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
350
Shear and rotation Shear only
Stress, T12 (MPa)
300 250 200 150 100 50 0
0
0.05
0.1
Amount of shear, Γ
0.15
0.2
Fig. 6. Comparison of results for T 12 versus C curves from shear with superposed rotation against shear with no superposed rotation.
degrade for shear strain increments of around 1 103 ; accordingly, we chose Dcps to be half this value, 0:5 103 . The automatic time-stepping algorithm operated to keep the ratio
R
Dcpmax Dcps
ð10:6Þ
close to 1.0 by adjusting the size of the time increments.14 As the curved bar is cyclically displaced in the ðe2 ; e3 Þ-plane at the ‘‘moved-face”, the section of the bar near the ‘‘fixed”-face is subjected to reversed-torsion, while the part of the bar towards the ‘‘moved”-face is subjected to reversed-bending, resulting in states of shear, tension, and compression in various parts of the body. With increasing displacement magnitude at the ‘‘moved”-face, inhomogeneous plastic deformation initiates and spreads in various disparate regions of the curved bar. The integrated reaction force in the 2-direction on the ‘‘moved”-face is plotted against the prescribed u2 -displacement in Fig. 7(b). This resulting overall cyclic load–displacement curve is the result of the complex evolution of the internal variables controlling the combined isotropic and kinematic hardening at the various sections of the bar which are undergoing plastic deformation. The overall load versus displacement response is smooth, and one can clearly see the effects of both kinematic hardening in the prominent Bauschinger effect and isotropic hardening in the continual increase in the magnitude of the load with cyclic deformation. This example shows that our semi-implicit integration procedure, when coupled with our heuristic automatic time-stepping algorithm and our approximate algorithmic tangent, is effective in
14 After an equilibrium solution for a time increment Dtn ¼ t nþ1 tn was found, the value of the ratio R was checked to determine whether this solution would be accepted or not. If R was greater than 1.25, then the solution was rejected, and a new time increment was used that was smaller by the factor ð0:85=RÞ. If R 6 1:25, then the solution was accepted, and the value of R was used to determine the first trial size of the next time increment. The following algorithm was used:
if if if
0:8 < R 6 1:25 then Dt nþ1 ¼ Dtn =R; 0:5 < R 6 0:8 then Dt nþ1 ¼ 1:25Dt n ; R 6 0:5 then Dt nþ1 ¼ 1:50Dtn :
ð10:7Þ
Note that the measure Dcpmax was allowed to exceed the user specified value Dcps by up to 25%. This was done to avoid having to recalculate increments that came out just slightly above the specified nominal value, but were otherwise essentially acceptable.
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1865
a
b
150 100
Load (kN)
50 0 −50 −100 −150
−20
−10
0
10
20
Displacement, u , (mm) 2
Fig. 7. (a) Undeformed and deformed meshes for the cyclic loading of a curved bar. (b) Resulting tip load versus tip displacement.
obtaining an accurate, stable, and efficient solution in an inhomogeneous large deformation implicit finite element calculation. 10.4. Bending and spring-back of an aluminum sheet As our final numerical example, we consider a simple plane-strain sheet forming operation in which an aluminum sheet is plastically bent about a mandrel of fixed radius and then unloaded. We use the recent experimental data of Cao et al. (2008) for AA6111-T4 sheets to estimate the material parameters in our model. The fit of the model to the experimental data is shown in Fig. 8, and the estimated material parameters are listed in (10.8):
G ¼ 25:5 GPa;
K ¼ 68 GPa;
m0 ¼ 0:001 s1 ;
m ¼ 0:02;
S0 ¼ 67:55 MPa; h0 ¼ 1200 MPa; a ¼ 1; S ¼ 156:7 MPa; B ¼ 1:09 GPa;
c ¼ 61:08:
ð10:8Þ
1866
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
400 Experiment Model
True axial stress (MPa)
300 200 100 0 −100 −200 −300 −400
0
0.02
0.04
0.06
0.08
0.1
0.12
True axial strain Fig. 8. Fit of the model to experimental data for AA6111-T4 sheet material. The data is from Cao et al. (2008); note that the kink in the reversed loading data for the smallest strain amplitude is actually reported by these authors, and is probably an artifact of the experimental difficulties associated with performing reversed loading experiments on sheet materials.
The initial configuration for our plane-strain sheet forming simulation is shown in Fig. 9(a). We consider a 1 mm thick sheet with a length of 75 mm, which is to be to bent to a radius of 50 mm between a pair of matched-rigid dies; accordingly, the radius of the top and bottom dies is prescribed to be 50.5 mm and 49.5 mm, respectively. Due to the symmetry of the problem, we consider only half of the geometry with suitable boundary conditions at the symmetry-plane. The sheet is modeled using a mesh consisting of 1880 Abaqus-CPE4H plane-strain elements with seven elements through the sheet thickness, and the dies are modeled using rigid surfaces. Contact between the sheet and the dies was modeled as frictionless. In the first step of the sheet forming simulation, the top die is moved downward to bend the sheet completely around the bottom die. Fig. 9(b) shows the geometry of the initial and fully bent geometries of the sheet. In the second step, the top die is moved back upwards, and the sheet is allowed to spring-back. The unloaded sprung-back geometry is also shown in Fig. 9(b). It is widely believed that plasticity models which do not account for kinematic hardening tend to incorrectly predict the amount of spring-back in simulations of sheet-metal forming operations (cf., e.g., Zhao and Lee, 2001). Thus, our large deformation theory which not only includes both kinematic and isotropic hardening but also accounts for large rotations, should be useful in improved modeling of spring-back phenomena in sheet forming operations.
11. Concluding remarks In this paper, we have formulated a large deformation constitutive theory for combined isotropic and kinematic hardening based on the dual decomposition F ¼ Fe Fp and Fp ¼ Fpen Fpdis . We now show that a similar kinematic-hardening theory may also be constructed without using the decomposition Fp ¼ Fpen Fpdis . Recall the relation (6.3) for the evolution of Bpen ,
B_ pen ¼ Dpen Bpen þ Bpen Dpen ; and the kinematical equation (2.28)
ð11:1Þ
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1867
a
b
5
x2 (mm)
0 −5 Undeformed Fully bent Sprung−back
−10 −15 −20
0
10
20
x (mm)
30
40
Fig. 9. (a) Initial configuration for plane-strain matched-die sheet bending simulations. (b) Undeformed and fully bent sheet geometries, together with the sprung-back geometry.
: Dp ¼ Dpen þ sym Fpen Dpdis Fp1 en
ð11:2Þ
The relation (11.2) gives
p1 Dpen ¼ Dp sym Fpen Dpdis Fen ; p mdis p1 ¼ Dp ðusing ð9:63ÞÞ; sym Fpen Mpen Fen en 2s p mdis p p1 ðusing ð9:60ÞÞ; sym Fpen Rp> ¼ Dp en Sen R en Fen 2sen p mdis p1 ¼ Dp ; sym Vpen Sen Ven 2sen p mdis ðusing ð9:57Þ and Epen ð1=2Þ ln Bpen Þ; sym Vpen ðB ln Bpen ÞVp1 ¼ Dp en 2sen
ð11:3Þ
1868
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
p Bmdis ðln Bpen Þ; en 2s 1 p ¼ Dp cm ðln Bpen Þ ðusing ð9:42ÞÞ: 2 ¼ Dp
ð11:4Þ
Using (11.4) in (11.1) gives
1 p 1 B_ pen ¼ Dp cm ðln Bpen Þ Bpen þ Bpen Dp cmp ðln Bpen Þ ; 2 2 or
B_ pen ¼ Dp Bpen þ Bpen Dp cmp Bpen ðln Bpen Þ:
ð11:5Þ
Thus, if we set
A Bpen ; we obtain the evolution equation
A_ ¼ Dp A þ ADp cmp Aðln AÞ:
ð11:6Þ p
¼ Fpen Fpdis wðpÞ ðBpen Þ,
Therefore, instead of using the decomposition F and developing a theory which gives rise to one may develop an alternate theory based a back-stress Sen because of a defect free-energy ðpÞ ðAÞ, and an evolution equation on an internal variable A, with a corresponding defect free-energy w for A obeying (11.6); indeed, this is the track taken recently by Anand et al. (2008). However, we do note that the development of an implicit time-integration procedure for a constitutive theory based on A is not as straightforward as the one developed here for a theory utilizing the Fp ¼ Fpen Fpdis decomposition. Finally, as is abundantly clear from the extensive literature on hardening models for metal plasticity (cf., e.g., Chaboche, 2008), the simple theory with combined isotropic and kinematic hardening developed in this paper is only foundational in nature, and there are numerous specialized enhancements/modifications to the theory that need to be incorporated in order to match actual experimental data for different metals; we leave such work to the future. Acknowledgements This work was supported by the National Science Foundation under Grant Nos. CMS-0555614 and CMMI-062524. Appendix A In this Appendix we develop a semi-implicit time-integration procedure, and obtain an (approximate) algorithmic tangent for the theory formulated in this paper. The time-integration procedure and the associated tangent are intended for use in the context of ‘‘implicit” finite element procedures. We have implemented the procedure described below in the commercial finite element program Abaqus/ Standard (2008) by writing a user material subroutine (UMAT), and it is using this finite element program that we have carried out the calculations whose results are reported in the main body of the paper.
A.1. Time-integration procedure The constitutive time-integration problem is that given
fTn ; Fpn ; Sn ; ðSen Þn ; ðFpdis Þn g; as well as Fn and Fnþ1 ; at time t n , we need to calculate
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1869
fTnþ1 ; Fpnþ1 ; Snþ1 ; ðSen Þnþ1 ; ðFpdis Þnþ1 g at time tnþ1 ¼ tn þ Dt. The evolution equation for F_ p ¼ Dp Fp is integrated by means of an exponential map (Weber and Anand, 1990)
Fpnþ1 ¼ exp DtDpnþ1 Fpn ; the inverse of
Fp1 nþ1
¼
Fpnþ1
Fp1 n
^ p ðððMe Þ Þ ; Snþ1 Þ; with Dpnþ1 ¼ D eff nþ1 0 nþ1
ð12:1Þ
is then
exp DtDpnþ1 :
e
p1
Next, using F ¼ FF given by
ð12:2Þ
and (12.1)1 , the elastic deformation gradient at the the end of the step is
Fenþ1 ¼ Fetrial exp DtDpnþ1 ;
ð12:3Þ
where def
Fetrial ¼ Fnþ1 Fp1 n
ðtrial values are evaluated with plastic flow frozenÞ
ð12:4Þ
e
is a trial value of F at the end of the step. The elastic right Cauchy–Green tensor and its trial value at the end of the step are e Cenþ1 ¼ Fe> nþ1 Fnþ1 ;
e Cetrial ¼ Fe> trial Ftrial :
ð12:5Þ
Thus, using (12.3),
Cenþ1 ¼ exp DtDpnþ1 Cetrial exp DtDpnþ1 :
ð12:6Þ
To proceed further we make our first approximation: (A1) To first order in Dt, we approximate
: exp DtDpnþ1 ¼ 1 DtDpnþ1 :
ð12:7Þ
Hence, (12.6) becomes
: Cenþ1 ¼ 1 DtDpnþ1 Cetrial 1 DtDpnþ1 ; : ¼Cetrial DtDpnþ1 Cetrial DtCetrial Dpnþ1 ; : e1 p ¼Cetrial 1 DtCtrial Dnþ1 Cetrial DtDpnþ1 :
ð12:8Þ
e e p eG 1 Next, consider the term Ce1 trial Dnþ1 Ctrial ; with Etrial ¼ 2 ðCtrial 1Þ denoting a Green strain corresponde e eG ing to Ctrial , we have Ctrial ¼ 1 þ 2Etrial , and hence
: e1 p p eG Ctrial Dnþ1 Cetrial ¼ 1 2EeG trial Dnþ1 1 þ 2Etrial ; : eG p eG p eG ¼ Dpnþ1 þ 2Dpnþ1 EeG trial 2Etrial Dnþ1 4Etrial Dnþ1 Etrial :
ð12:9Þ
We now make our second approximation: (A2) We assume that jEeG trial j 1, so that (12.9) may be approximated as
: e1 p Ctrial Dnþ1 Cetrial ¼ Dpnþ1 :
ð12:10Þ
Using (12.10), we find that (12.8) reduces to
: Cenþ1 ¼ Cetrial 1 2DtDpnþ1 :
ð12:11Þ
Further, using the symmetry of Cenþ1 we note that the symmetric tensors Cetrial and 1 2DtDpnþ1 commute, and hence share the same principal directions. Thus taking the logarithm of (12.11) we obtain
1870
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1 : Eenþ1 ¼Eetrial þ ln 1 2DtDpnþ1 ; 2
ð12:12Þ
where
Eenþ1 ¼
1 ln Cenþ1 2
and Eetrial ¼
1 ln Cetrial : 2
ð12:13Þ
Next, we make our third approximation: (A3) Recalling the series representation of the logarithm of a tensor, we assume that
: ln 1 2DtDpnþ1 ¼ 2DtDpnþ1 for small DtDpnþ1 . Thus, under the assumptions (A1)–(A3) of small Dt and small jEeG trial j, (12.6) yields the important relation
: Eenþ1 ¼Eetrial DtDpnþ1 :
ð12:14Þ
Next, let
1 def C ¼ 2G I 1 1 þ K1 1 3
ð12:15Þ
denote the elasticity tensor, and let
Menþ1 ¼ C½Eenþ1 and Metrial ¼ C½Eetrial ;
ð12:16Þ
respectively, denote the elastic Mandel stress at the end of the time step, as well as its trial value. Then, operating on (12.14) by C gives
Menþ1 ¼ Metrial DtC½Dpnþ1 :
ð12:17Þ
Subtracting ðSen Þnþ1 from the left-hand side and
ðSen Þtrial ðSen Þn ;
15
ðtrial values are evaluated with plastic flow frozenÞ
from the right-hand side of (12.17), and writing def
ðMeeff Þtrial ¼ Metrial ðSen Þtrial ;
ð12:18Þ
ðMeeff Þnþ1 ¼ ðMeeff Þtrial DtC½Dpnþ1 :
ð12:19Þ
gives
Next, from (9.62) we have
1 Dpnþ1 ¼ pffiffiffi mpnþ1 Npnþ1 ; 2
Npnþ1 ¼
ððMeeff Þnþ1 Þ0 pffiffiffi ; nþ1 2s
1
snþ1 ¼ pffiffiffi jððMeeff Þnþ1 Þ0 j; 2
ð12:20Þ
where
pffiffiffi
mpnþ1 ¼ 2jDpnþ1 j:
ð12:21Þ
Using (12.15) and (12.20) in (12.19) we obtain
ðMeeff Þnþ1 ¼ ðMeeff Þtrial Since
Npnþ1
pffiffiffi 2GðDt mpnþ1 ÞNpnþ1 :
ð12:22Þ
is deviatoric, the deviatoric and spherical parts of (12.22) give
15 Subtracting ðSen Þn rather than ðSen Þnþ1 from the right-hand side of (12.17) to arrive at (12.19), is a critical step in our timeintegration procedure. It allows us to decouple the update of Fp and Fpdis into to two separate, yet otherwise implicit updates.
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
ðMeeff Þnþ1
0
¼ ððMeeff Þtrial Þ0
1871
pffiffiffi 2GðDt mpnþ1 ÞNpnþ1 ;
trðMeeff Þnþ1 ¼ trðMeeff Þtrial :
ð12:23Þ
Using (12.20)2 , (12.23)1 may be arranged as
pffiffiffi nþ1 þ GðDtmpnþ1 Þ Npnþ1 ¼ ððMeeff Þtrial Þ0 : 2 s
ð12:24Þ
Next, defining def
Nptrial ¼
ððMeeff Þtrial Þ0 pffiffiffi ; trial 2s
1
strial def ¼ pffiffiffi jððMeeff Þtrial Þ0 j;
ð12:25Þ
2
(12.24) may be written as
snþ1 þ GðDtmpnþ1 Þ Npnþ1 ¼ strial Nptrial ;
ð12:26Þ
which immediately gives
Npnþ1 ¼ Nptrial ; snþ1 þ GðDtmpnþ1 Þ ¼ strial :
ð12:27Þ Npnþ1
Thus the direction of plastic flow at the end of the step is determined by the trial direction of plastic flow Nptrial . Next, the implicit form of the flow strength relation (9.30) is
snþ1 ¼ Y mpnþ1 ; Snþ1 :
ð12:28Þ
We integrate the evolution equation for S in a semi-implicit fashion as
Snþ1 ¼ Sn þ hðSn ÞðDtmpnþ1 Þ:
ð12:29Þ
Since the hardening function h typically does not change rapidly, we do not expect the semi-implicit integration of the evolution equation for S to have a significant effect on the stability of our algorithm. Using (12.29) in (12.28), we have
snþ1 ¼ Y mpnþ1 ; Sn þ hðSn ÞðDtmpnþ1 Þ :
ð12:30Þ
Finally, using (12.30) in (12.27)2 gives the following implicit equation for
trial GðDt mpnþ1 Þ Y mpnþ1 ; Sn þ hðSn ÞðDtmpnþ1 Þ ¼ 0: Gðmpnþ1 Þ ¼ s
mpnþ1 : ð12:31Þ
Once a solution for pnþ1 elastic Mandel stress Menþ1
is obtained by solving (12.31), Snþ1 is easily evaluated by using (12.29). The is obtained from (12.17) using
1 Dpnþ1 ¼ pffiffiffi mpnþ1 Nptrial ; 2
ð12:32Þ
m
as
Menþ1 ¼ Metrial
pffiffiffi 2GðDtmpnþ1 ÞNptrial :
ð12:33Þ
p
The update for F is obtained from (12.1)1 using the result (12.32):
1 Fpnþ1 ¼ exp pffiffiffi ðDt mpnþ1 ÞNptrial Fpn : 2
ð12:34Þ
Next, calculating
Fenþ1 ¼ Fnþ1 Fp1 nþ1 ;
ð12:35Þ Fenþ1
Renþ1 .
and performing a polar decomposition of gives and (9.55) gives the update for the Cauchy stress as > e e e Tnþ1 ¼ J 1 nþ1 ðRnþ1 ÞðMnþ1 ÞðR nþ1 Þ
Finally, using
with J nþ1 ¼ det Fnþ1 :
Renþ1 ,
and the relations (12.33)
ð12:36Þ
1872
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
We turn next to updating Sen and Fpdis . The evolution equation for F_ pdis ¼ Dpdis Fpdis is also integrated by using an exponential map and (9.63) as
ðFpdis Þnþ1 ¼ exp DtðDpdis Þnþ1 ðFpdis Þn
with ðDpdis Þnþ1 ¼
c mp ðMp Þ ; 2B nþ1 en nþ1
ð12:37Þ
the inverse of ðFpdis Þnþ1 is then
p 1 p ðFpdis Þ1 nþ1 ¼ ðFdis Þn exp DtðDdis Þnþ1 :
Hence, using
Fpen
ðFpen Þnþ1
¼
¼
Fp Fp1 dis ,
ðFpen Þtrial
ð12:38Þ p
the energetic part of F at the end of the step is given by
exp DtðDpdis Þnþ1 ;
ð12:39Þ
where def
ðFpen Þtrial ¼ Fpnþ1 ðFpdis Þ1 n is
a trial value of Fpen at The tensors ðFpen Þnþ1
ð12:40Þ
the end of the step. and ðFpen Þtrial admit the polar decompositions
ðFpen Þnþ1 ¼ ðRpen Þnþ1 ðUpen Þnþ1 ;
ðFpen Þtrial ¼ ðRpen Þtrial ðUpen Þtrial :
ð12:41Þ
Using (12.41) in (12.39) and rearranging, we obtain
ðRpen Þnþ1 ðUpen Þnþ1 exp DtðDpdis Þnþ1 ¼ ðRpen Þtrial ðUpen Þtrial :
ð12:42Þ
Next, from (9.57) and (9.60) p p> p p p Mpen ¼ Rp> en Sen R en ¼ R en ð2B ln V en ÞRen ¼ 2B ln Uen :
Thus the principal directions of those ðUpen Þnþ1 . Hence
ðMpen Þnþ1 ,
ð12:43Þ
and hence from (12.37)2 those of
ðDpdis Þnþ1 ,
ðUpen Þnþ1 exp DtðDpdis Þnþ1 is symmetric:
are the same as
ð12:44Þ
Then, because of the uniqueness of the polar decomposition theorem
ðRpen Þnþ1 ¼ ðRpen Þtrial ; ðUpen Þnþ1 exp DtðDpdis Þnþ1 ¼ ðUpen Þtrial :
ð12:45Þ ð12:46Þ
Eq. (12.46) implies that ðUpen Þnþ1 and ðUpen Þtrial have the same principal directions. Thus taking the logarithm on both sides and rearranging we have
lnðUpen Þnþ1 ¼ lnðUpen Þtrial DtðDpdis Þnþ1 :
ð12:47Þ
Multiplying (12.47) through by 2B and using (12.43) gives
ðMpen Þnþ1 ¼ ðMpen Þtrial 2BDtðDpdis Þnþ1 ;
def
where ðMpen Þtrial ¼ 2B lnðUpen Þtrial :
Finally, using (9.60), (12.37)2 , (12.45), and pre-multiplying (12.48) by ðRpen Þ> trial gives
ðSen Þnþ1 ¼ ðSen Þtrial cðDtmpnþ1 ÞðSen Þnþ1 ;
def
ðRpen Þtrial
ð12:48Þ
and post-multiplying by
where ðSen Þtrial ¼ 2B lnðVpen Þtrial :
ð12:49Þ
Thus, the back-stress Sen is updated as
! 1 ðSen Þtrial : 1 þ cðDtmpnþ1 Þ
ð12:50Þ
ðMpen Þnþ1 ¼ ðRpen Þ>trial ðSen Þnþ1 ðRpen Þtrial ;
ð12:51Þ
ðSen Þnþ1 ¼ Correspondingly,
use of which in (12.37) provides the update ðFpdis Þnþ1 for Fpdis .
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1873
Remark. Due to the use of the exponential map in integrating the evolution equations for Fp and Fpdis , the constraint of plastic incompressibility is exactly maintained by our time-integration procedure. This is easily verified by recalling the identity detðexp AÞ ¼ expðtrAÞ and recognizing the deviatoric nature of Dp and Dpdis in (12.1) and (12.37), respectively. A.1.1. Summary of the time-integration procedure Given: fTn ; Fpn ; Sn ; ðSen Þn ; ðFpdis Þn g, as well as Fn and Fnþ1 , at time tn . Calculate: fTnþ1 ; Fpnþ1 ; Snþ1 ; ðSen Þnþ1 ; ðFpdis Þnþ1 g at time t nþ1 ¼ t n þ Dt. Step 1. Calculate the trial elastic deformation gradient
Fetrial ¼ Fnþ1 Fp1 n :
ð12:52Þ
Step 2. Perform the polar decomposition
Fetrial ¼ Retrial Uetrial :
ð12:53Þ
Step 3. Calculate the trial elastic strain
Eetrial ¼ ln Uetrial :
ð12:54Þ
Step 4. Calculate the trial Mandel stress and the trial effective Mandel stress
Metrial ¼ C½Eetrial ; ðMeeff Þtrial ¼ Metrial ðSen Þn :
ð12:55Þ ð12:56Þ
Step 5. Calculate the trial mean normal pressure, the deviatoric part of the trial effective Mandel stress, the trial equivalent shear stress, and the trial direction of plastic flow
1 trial ¼ trðMeeff Þtrial ; p 3 trial 1; ððMeeff Þtrial Þ0 ¼ ðMeeff Þtrial þ p rffiffiffi 1 strial ¼ ððMeeff Þtrial Þ0 ; 2 e ððM Þ Þ ffiffiffi trial 0 : Nptrial ¼ peff trial 2s
ð12:57Þ ð12:58Þ ð12:59Þ ð12:60Þ
mpnþ1 by solving
Step 6. Calculate
trial GðDt mpnþ1 Þ Y Gðmpnþ1 Þ ¼ s
mpnþ1 ; Sn þ hðSn ÞðDtmpnþ1 Þ ¼ 0:
ð12:61Þ
p
Step 7. Update D
Dpnþ1 ¼
pffiffiffiffiffiffiffiffiffiffi p ð1=2Þmnþ1 Nptrial :
ð12:62Þ
p
Step 8. Update F
Fpnþ1 ¼ exp DtDpnþ1 Fpn :
ð12:63Þ
Step 9. Update the Mandel stress Me
Menþ1 ¼ Metrial
pffiffiffi 2GðDtmpnþ1 ÞNptrial :
ð12:64Þ
e
Step 10. Update F
Fenþ1 ¼ Fnþ1 Fp1 nþ1 :
ð12:65Þ
Step 11. Perform the polar decomposition
Fenþ1 ¼ Renþ1 Uenþ1 :
ð12:66Þ
1874
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
Step 12. Update the Cauchy stress
J nþ1 ¼ det ðFnþ1 Þ;
ð12:67Þ
> e e e Tnþ1 ¼ J 1 nþ1 ðR nþ1 ÞðMnþ1 ÞðR nþ1 Þ :
ð12:68Þ
Step 13. Update S
Snþ1 ¼ Sn þ hðSn ÞðDt mpnþ1 Þ:
ð12:69Þ
Step 14. Calculate ðFpen Þtrial
ðFpen Þtrial ¼ Fpnþ1 ðFpdis Þ1 n :
ð12:70Þ
Step 15. Perform the polar decomposition
ðFpen Þtrial ¼ ðVpen Þtrial ðRpen Þtrial :
ð12:71Þ
Step 16. Calculate the trial back-stress
ðSen Þtrial ¼ 2B lnðVpen Þtrial :
ð12:72Þ
Step 17. Update the back-stress Sen
ðSen Þnþ1 ¼
! 1 ðSen Þtrial : 1 þ cðDtmpnþ1 Þ
ð12:73Þ
Step 18. Update the plastic Mandel stress Mpen
ðMpen Þnþ1 ¼ ðRpen Þ>trial ðSen Þnþ1 ðRpen Þtrial : Step 19. Update
ðDpdis Þnþ1 ¼
ð12:74Þ
Dpdis
c mp ðMp Þ : 2B nþ1 en nþ1
ð12:75Þ
Step 20. Update Fpdis
ðFpdis Þnþ1 ¼ exp DtðDpdis Þnþ1 ðFpdis Þn :
ð12:76Þ
A.2. Jacobian matrix In typical ‘‘implicit” finite element procedures utilizing a Newton-type iterative solution method, one needs to compute an algorithmically consistent tangent, often called the Jacobian matrix. We obtain an estimate for our Jacobian matrix below. From the outset we note that Jacobian matrices are used only in the search for the global finite element solution, and while an approximate Jacobian might affect the rate of convergence of the global iteration scheme, it will not impair the accuracy of our constitutive time-integration algorithm. Consider the Cauchy stress at the end of the increment:
Tnþ1 ¼ ðdet Fenþ1 Þ1 Renþ1 Menþ1 ðRenþ1 Þ> :
ð12:77Þ
Then with D denoting a variation, (12.77) gives
DTnþ1 ðDRenþ1 ÞðRenþ1 Þ> Tnþ1 þ Tnþ1 ðDRenþ1 ÞðRenþ1 Þ> þ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} incremental elastic rotation
¼
incremental elastic rotation
ðdet Fenþ1 Þ1 Renþ1 DMenþ1 ðRenþ1 Þ> :
e1 tr DFenþ1 Fnþ1 |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
Tnþ1
incremental elastic vol: change
ð12:78Þ
of the incremental elastic rotaWe assume that the variations and tr DFenþ1 Fe1 nþ1 tion and volume change are reasonably well-estimated by a commercial finite element program, such
ðDRenþ1 ÞðRenþ1 Þ>
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1875
as Abaqus/Standard (2008), in which we have implemented our constitutive model by writing a user material subroutine (UMAT).16. Thus, here we concentrate on evaluating the variation DMenþ1 in (12.78), which is computed from
DMenþ1 ¼ C DEetrial ;
ð12:79Þ
where the fourth-order tensor
@Menþ1 @Eetrial
def
C¼
ð12:80Þ
is the important constitutive contribution to the global Jacobian matrix. Recall that the time-integration procedure gives
Menþ1 ¼ Metrial
pffiffiffi 2GðDtmpnþ1 ÞNptrial :
ð12:81Þ
Using (12.27)2 , (12.81) may be written as
Menþ1 ¼ Metrial þ
pffiffiffi nþ1 s trial ÞNptrial : 2ðs
ð12:82Þ
By the product rule
C¼
pffiffiffi nþ1 @ s trial @Nptrial @Metrial pffiffiffi @s p þ 2 ð s s Þ þ N 2 : nþ1 trial trial @Eetrial @Eetrial @Eetrial @Eetrial
ð12:83Þ
Since
Metrial ¼ C Eetrial ;
ð12:84Þ
we have
@Metrial ¼ C: @Eetrial
ð12:85Þ
Next, since
Nptrial ¼
ððMeeff Þtrial Þ0 pffiffiffi ; trial 2s
ð12:86Þ
we have
@Nptrial ¼ @Eetrial
rffiffiffi trial 1 1 @ððMeeff Þtrial Þ0 ððMeeff Þtrial Þ0 @ s e : e 2 trial 2 s strial @Etrial @Etrial
ð12:87Þ
Now,
ðMeeff Þtrial ¼ Metrial ðSen Þn ¼ C Eetrial ðSen Þn 1 ¼ 2G I 1 1 þ K1 1 Eetrial ðSen Þn : 3
ð12:88Þ
Hence,
e 1 ððMeeff Þtrial Þ0 ¼ 2G I 1 1 Etrial ðSen Þn ; 3
16
Abaqus uses a hypoelastic constitutive equation for the stress of the type
e p _ _ T ¼ TW TþTWe ¼ TWTþTW |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} ¼ C½DD : for Wp ¼0
ð12:89Þ
1876
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
and therefore
@ððMeeff Þtrial Þ0 1 ¼ 2G I 1 1 : e 3 @Etrial
ð12:90Þ
Further,
strial
rffiffiffi 1 ððMeeff Þtrial Þ0 : ¼ 2
ð12:91Þ
Therefore, by the chain-rule
rffiffiffi > 1 @ððMeeff Þtrial Þ0 @jððMeeff Þtrial Þ0 j e e 2 @E @ððMeff Þtrial Þ0 rffiffiffi trial > 1 1 ððMeeff Þtrial Þ0 2G I 1 1 ; ¼ 2 3 jððMeeff Þtrial Þ0 j
pffiffiffi 1 ððMeeff Þtrial Þ0 ¼ 2G I 1 1 3 jððMeeff Þtrial Þ0 j pffiffiffi ððMe Þ Þ eff trial 0 ¼ 2G jððMeeff Þtrial Þ0 j pffiffiffi ¼ 2GNptrial :
trial @s ¼ @Eetrial
ð12:92Þ ð12:93Þ ð12:94Þ ð12:95Þ ð12:96Þ
Thus, using (12.90) and (12.96) in (12.87) we obtain
@Nptrial 2G ¼ pffiffiffi @Eetrial trial 2s
1 I 1 1 Nptrial Nptrial ; 3
ð12:97Þ
and substituting (12.85) and (12.97) in (12.83), we have
C ¼Cþ
1 I 1 1 Nptrial Nptrial 3 pffiffiffi nþ1 pffiffiffi p @s 2 2GNtrial : þ Nptrial e @Etrial
snþ1 1 2G strial
ð12:98Þ
trial , we have, using (12.96) Noting that Eetrial enters the equations through s
nþ1 @ s nþ1 @ s trial @ s nþ1 pffiffiffi p @s ¼ ¼ 2GNtrial : e e trial @Etrial @ strial @Etrial @ s
ð12:99Þ
Finally, substituting (12.99) into (12.98), using (12.15), and rearranging, we have
~ 2G snþ1 @ snþ1 Np Np C¼C trial strial @ strial trial
ð12:100Þ
snþ1 ~ I 1 1 1 þ K1 1 with G ~ def ~ ¼ 2G C G: ¼ 3 strial
ð12:101Þ
where
Thus, it remains to determine the derivative
nþ1 @s trial @s nþ1 ; (12.27)2 , and Snþ1 , (12.69), as in (12.100). First, we rewrite the updates for s
; Snþ1 snþ1 ; ¼ Sn þ DthðSn Þf Snþ1
snþ1 ¼ strial GDtf Snþ1
snþ1
ð12:102Þ
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
1877
respectively, where f is the flow function of (9.33). Using the definitions def
X¼
snþ1 Snþ1
;
def
trial g; and Y ¼ f s
ð12:103Þ
the system of equations (12.102) may be written as
X ¼ GðX; YÞ;
ð12:104Þ
8 9 > = < strial GDtf sSnþ1 > nþ1 : G¼ > ; : Sn þ DthðSn Þf sSnþ1 > nþ1
ð12:105Þ
where
Differentiating (12.104) with respect to Y (at the solution point) we obtain
@X @GðX; YÞ @GðX; YÞ @X ¼ þ ; @Y @Y @X @Y
ð12:106Þ
from which
@X @GðX; YÞ ¼ A ; @Y @Y
ð12:107Þ
where17
A ¼
1 @GðX; YÞ I : @X
ð12:108Þ
Straightforward calculations show
2
GDt@@fs
@GðX; YÞ 6 tnþ1 ¼4 @f @X DthðSn Þ@ s
tnþ1
@f GDt @S
3 tnþ1
@f DthðSn Þ@S
7 5;
and
1 @GðX; YÞ ¼ : @Y 0
ð12:109Þ
tnþ1
Thus, we obtain that
nþ1 @s ¼ A11 ; trial @s
ð12:110Þ
substitution of which in (12.100) completes our estimate for the Jacobian matrix C. References Abaqus/Standard, 2008. SIMULIA, Providence, RI. Anand, L., 1979. On H. Hencky’s approximate strain-energy function for moderate deformations. ASME Journal of Applied Mechanics 46, 78–82. Anand, L., 1986. Moderate deformations in extension-torsion of incompressible isotropic elastic materials. Journal of the Mechanics and Physics of Solids 34, 293–304. Anand, L., Gurtin, M.E., 2003. A theory of amorphous solids undergoing large deformations, with applications to polymeric glasses. International Journal of Solids and Structures 40, 1465–1487. Anand, L., Ames, N.M., Srivastava, V., Chester, S., 2008. A thermo-mechanically-coupled theory for large deformations of amorphous polymers. Part 1: Formulation. International Journal of Plasticity, doi:10.1016/j.ijplas.2008.11.004. Armstrong, P.J., Frederick, C.O., 1966. A mathematical representation of the multiaxial Bauschinger effect. Report RD/B/N731, CEGB, Central Electricity Generating Board, Berkeley, UK. Bauschinger, J., 1886. Über die Veränderung der Position der Elastizitätsgrenze des Eisens und Stahls durch Strecken und Quetschen und durch Erwärmen und Abkühlen und durch oftmals wiederholte Beanspruchungen. Mitteilung aus dem Mechanisch-technischen Laboratorium der Königlichen polytechnischen Hochschule in München 13, 1–115. Brown, S.B., Kim, K.H., Anand, L., 1989. An internal variable constitutive model for hot working of metals. International Journal of Plasticity 5, 95–130.
17
Here I is the second-order identity matrix.
1878
D.L. Henann, L. Anand / International Journal of Plasticity 25 (2009) 1833–1878
Cao, J., Lee, W., Cheng, H.S., Seniw, M., Wang, H.-P., Chung, K., 2008. Experimental and numerical investigation of combined isotropic-kinematic hardening behavior of sheet metals. International Journal of Plasticity, 10.1016/j.ijplas.2008.04.007. Casey, J., Naghdi, P.M., 1980. A remark on the use of the decomposition F ¼ Fe Fp in plasticity. Journal of Applied Mechanics 47, 672–675. Chaboche, J.L., 2008. A review of some plasticity and viscoplasticity constitutive theories. International Journal of Plasticity 24, 1642–1693. Clayton, J., McDowell, D., 2003. A multiscale multiplicative decomposition for elastoplasticity of polycrystals. International Journal of Plasticity 19, 1401–1444. Dafalias, Y., 1998. Plastic spin: necessity or redundancy? International Journal of Plasticity 14, 909–931. Dashner, P.A., 1986. Invariance considerations in large strain elasto-plasticity. Journal of Applied Mechanics 53, 55–60. Dettmer, W., Reese, S., 2004. On the theoretical and numerical modelling of Armstrong–Frederick kinematic hardening in the finite strain regime. Computer Methods in Applied Mechanics and Engineering 193, 87–116. Frederick, C.O., Armstrong, P.J., 2007. A mathematical representation of the multiaxial Bauschinger effect. Materials at High Temperatures 24, 1–26. Green, A.E., Naghdi, P.M., 1971. Some remarks on elastic–plastic deformation at finite strain. International Journal of Engineering Science 9, 1219–1229. Gurtin, M.E., Anand, L., 2005. The decomposition F ¼ Fe Fp , material symmetry, and plastic irrotationality for solids that are isotropic-viscoplastic or amorphous. International Journal of Plasticity 21, 1686–1719. Håkansson, P., Wallin, M., Ristinmaa, M., 2005. Comparison of isotropic hardening and kinematic hardening in thermoplasticity. International Journal of Plasticity 21, 1435–1460. Kröner, E., 1960. Allgemeine kontinuumstheorie der versetzungen und eigenspannungen. Archive for Rational Mechanics and Analysis 4, 273–334. Lemaitre, J., Chaboche, J.L., 1990. Mechanics of Solid Materials. Cambridge University Press, Cambridge. Lion, A., 2000. Constitutive modelling in finite thermoviscoplasticity: a physical approach based on nonlinear rheological models. International Journal of Plasticity 16, 469–494. Lush, A.M., Weber, G., Anand, L., 1989. An implicit time-integration procedure for a set of internal variable constitutive equations for isotropic elasto-viscoplasticity. International Journal of Plasticity 5, 521–549. Menzel, A., Ekh, M., Runesson, K., Steinmann, P., 2005. A framework for multiplicative elastoplasticity with kinematic hardening coupled to anisotropic damage. International Journal of Plasticity 21, 397–434. Prager, W., 1949. Recent developments in the mathematical theory of plasticity. Journal of Applied Physics 20, 235–241. Ristinmaa, M., Wallin, M., Ottosen, N.S., 2007. Thermodynamic format and heat generation of isotropic hardening plasticity. Acta Mechanica 194, 103–121. Rosakis, P., Rosakis, A.J., Ravichandran, G., Hodowany, J., 2000. A thermodynamic internal variable model for the partition of plastic work into heat and stored energy in metals. Journal of the Mechanics and Physics of Solids 48, 581–607. Vladimirov, I.N., Pietryga, M.P., Reese, S., 2008. On the modelling of non-linear kinematic hardening at finite strains with application to springback – comparison of time integration algorithms. International Journal for Numerical Methods in Engineering 75, 1–28. Wallin, M., Ristinmaa, M., 2005. Deformation gradient based kinematic hardening model. International Journal of Plasticity 21, 2025–2050. Weber, G., Anand, L., 1990. Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic–viscoplastic solids. Computer Methods in Applied Mechanics and Engineering 79, 173–202. Weber, G.G., Lush, A.M., Zavaliangos, A., Anand, L., 1990. An objective time-integration procedure for isotropic rate-independent and rate-dependent elastic–plastic constitutive equations. International Journal of Plasticity 6, 701–744. Zhao, K.M., Lee, J.K., 2001. Material properties of aluminum alloy for accurate draw-bend simulation. Journal of Engineering Materials and Technology 123, 287–292.
