Aspects on the finite-element implementation of the Gurson model including parameter identification

International Journal of Plasticity 15 (1999) 1111±1137

Aspects on the ®nite-element implementation of the Gurson model including parameter identi®cation Rolf Mahnken* Institut fuÈr Baumechanik und Numerische Mechanik, University of Hannover Appelstrasse, 9a, 30167 Hannover, Germany Received in ®nal revised form 10 April 1999

Abstract In this contribution, various aspects on the ®nite-element implementation of the Gurson model are considered. In particular, a linear representation for the plastic potential is used, which shows superior convergence property in the local iteration procedure compared to the original quadratic representation. The formulation of the model is performed in the spatial con®guration based on the multiplicative decomposition of the deformation gradient, and for integration an exponential map scheme is used. A further important aspect is the sensitivity analysis consistent with the underlying integration scheme necessary for minimizing a leastsquares functional for parameter identi®cation by use of a gradient-based optimization algorithm. In a numerical example the local convergence behavior for the two versions of the Gurson model, linear and quadratic are compared. Furthermore material parameters are determined by least-squares minimization based on experimental data obtained for an axisymmetric tensile bar for a ferritic steel. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Multiplicative plasticity; Voids and inclusions; Parameter identi®cation; Optimization; Sensitivity analysis

1. Introduction Metallographic studies (Hancock and Mackenzi, 1976; Thomason, 1990) demonstrate, that ductile damage is basically characterized by three mechanisms of void growth: (i) nucleation of voids due to fracture of particle±matrix interface, failure of * Tel.: +49-511-762-3560; fax: +49-511-762-5496. E-mail address: [email protected] (R. Mahnken). 0749-6419/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S0749-6419(99)00029-7

1112

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

the particle or micro-cracking of the matrix surrounding the inclusion, (ii) growth of voids thus leading to an enlargement of existing holes and (iii) coalescence or microcracks linking of neighboring voids thus leading to vanishing load-carrying capacity of the material, as the void volume fraction approaches unity. A well known approach for modeling the above mentioned mechanisms is formulated by Gurson (1977), where he derived a yield potential for porous plastic materials from simple cell models. Modi®cations of this model have been proposed to improve the predictions at low void volume fractions (Tvergaard, 1981) and to provide a better representation of ®nal void coalescence (Tvergaard and Needleman, 1984). In this respect the model is based on experimental evidence and various applications have been done in the past to model void growth and ductile rupture (see e.g. Tvergaard, 1981; Tvergaard and Needleman, 1984; Bennani et al., 1993; Steglich and Brocks, 1997). The issue of parameter identi®cation and determination of the associated micromechanical material parameters for the Gurson model is still a rather new approach with no generally accepted recommendations. It can be done by combining metallurgical examinations with cell modeling and macroscopic testing results (Steglich and Brocks, 1997). In this work an approach based on macroscopic experimental testing is described. To this end the distance of the simulated data obtained from ®nite element computations to the experimental data is minimized in a least-squares sense by varying (some of) the material parameters. In this way non-uniform distributions of the state variables, such as the void volume fraction, within the specimen are taken into account for the identi®cation process. This approach, including the corresponding sensitivity analysis has been described for problems of plasticity by Mahnken and Stein (1996b) for the geometric linear theory, and in Mahnken and Stein (1997) it has been extended to the geometric nonlinear case based on an algorithmic formulation described by Simo (1992). In the present paper particularly we will resort to the following aspects related to the ®nite-element implementation of the modi®ed Gurson model including parameter identi®cation: 1. Use of a linear formulation for the yield potential which shows better convergence behavior in the local iteration as the quadratic formulation in the original formulation of Gurson. 2. Smooth extension of the cosh(x)-function appearing in the Gurson yield function for larger arguments of x thus avoiding the occurrence of over¯ow during the iterative local iteration. 3. Smoothing of the additional variable f  introduced by Tvergaard and Needleman (1984) thus yielding a smooth derivative for the Jacobian in the local iteration, for the algorithmic tangent modulus in the global equilibrium iteration and the sensitivity operator for the inverse parameter identi®cation algorithm. 4. Application of an implicit algorithm for integration of the constitutive equations as proposed by Steinmann et al. (1993) thus allowing for larger load steps as an explicit algorithm.

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1113

5. Sensitivity analysis of the resulting discretized ®nite element equations consistent with the underlying integration scheme w.r.t. material parameters necessary for using a gradient-based optimization algorithm for parameter identi®cation. An outline of this work is as follows: in Section 2 the constitutive equations of the damage model relative to the spatial con®guration are summarized. In Section 3 an exponential map integrator is used for integration of the resulting rate equations as proposed by Steinmann et al. (1993). The algorithmic tangent operator, needed for the application of a Newton-scheme in the equilibrium iteration, and the sensitivity load term, needed for the application of a gradient-based optimization scheme (e.g. Gauss±Newton method, Quasi-Newton method) for parameter identi®cation, are derived in Section 4. For illustrative purposes two examples are presented in Section 5: ®rstly the convergence behavior of the Newton method is compared for a typical local constitutive problem using the quadratic and the linear formulation of the Gurson plastic potential. Furthermore material parameters are determined by leastsquares minimization based on experimental data for an axisymmetric tensile bar made out of ferritic steel. 1.1. Notations Square brackets [.] are used throughout the paper to denote `function of' in order to distinguish from mathematical groupings with parenthesis (.). 2. Formulation of the model equations 2.1. Notation and kinematics of multiplicative elastoplasticity Let B0  IR3 be the reference con®guration of a continuum body B with smooth boundary @B0 ; I ˆ ‰t0 ; TŠ 2 IR‡ the time interval of interest and K a (vector) space of admissible material parameters. We consider the con®guration ®eld '‰X; t; Š in terms of the three independent variables X  B0 ; t  I and   K. Thus '‰; t; Š : B7!IR3 de®nes a nonlinear deformation map at time t  I for given material parameters   K mapping particles X  B0 to their actual position x  B in the deformed con®guration. The associated nonlinear deformation gradient F ˆ rx '‰X; t; Š with J ˆ detF de®nes a mapping of increments dX  TB0 of a locally de®ned tangent space TB0 associated to the undeformed con®guration to increments dx  TB of a locally de®ned tangent space TB associated to the deformed con®guration. Here and if not stated otherwise also in the subsequent presentation explicit indication of the arguments t and  is omitted. Furthermore we endow the tangent spaces TB0 and TB with co-variant (Riemann) metric tensors G[ ˆ …G] †ÿ1 and g[ ˆ …g] †ÿ1 , respectively. The time and parameter derivatives of the con®guration ®eld are denoted by @t ' and @ ', respectively. For the subsequent representation, we will also introduce

1114

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

kinematic expressions by replacing the time derivative @t …† or the parameter derivative @ …† by its variational counterparts @ …† and @
…†, respectively, where now @ …† and @
…† are shorthand notations for the GaÃteaux derivative at the point x ˆ '…X; t; † in the direction of the virtual variation u and increment for linearization
u, respectively. With the above notation di€erent kinematic variables and associated derivatives can be de®ned (see Mahnken and Stein, 1997). Some examples for velocity gradients, rate of deformation tensors and Lie derivative operators are as follows: d ˆ sym…g[
l†

1: l ˆ rx @t ' 2: l ˆ rx @ '

d ˆ sym…g[
l †

3: l ˆ rx @ '

d ˆ sym…g[
l †

4: l
ˆ rx @
'

d
ˆ sym…g[
l
†

L]t …† ˆ F
@t …Fÿ1
…†]
Fÿt †
Ft

L[t …† ˆ Fÿt
@t …Ft
…†[
F†
Fÿ1 ÿ
L] …† ˆ F
@ Fÿ1
…†]
Fÿt
Ft ÿ
L[ …† ˆ Fÿt
@ Ft
…†[
F
Fÿ1 ÿ
L] …† ˆ F
@ Fÿ1
…†]
Fÿt
Ft ÿ
L[ …† ˆ Fÿt
@ Ft
…†[
F
Fÿ1 ÿ
L]
…† ˆ F
@
Fÿ1
…†]
Fÿt
Ft ÿ
L[
…† ˆ Fÿt
@
Ft
…†[
F Fÿ1

…1†

where …†] and …†[ denote a contra-variant and co-variant tensor object, respectively. The underlying concept of multiplicative elastoplasticity assumes the decomposition F ˆ Fe
Fp

…2†

where Fe and Fp represent the elastic and plastic part of F, respectively, and which implies a plastic intermediate con®guration Bp as macro stress free. In what follows we will also consider the plastic right Cauchy±Green tensor de®ned on the reference con®guration and the elastic left Cauchy±Green tensor de®ned on the actual con®guration, respectively, i.e. Cp ˆ Ftp
G [
Fp ;

be ˆ Fe
G ]
Fte

…3†

where it is assumed, that the tangent space TB 0 on the intermediate con®guration is endowed with a co-variant (Riemann) metric tensor G [ ˆ …G ] †ÿ1 . The above fundamental contra-variant strain measures are related via the pull-back and push-forward operations ÿ1 ÿt Cÿ1 p ˆ F
be
F ;

] ÿ1 t t be ˆ F
Cÿ1 p
F ? Lt be ˆ F
@t Cp
F

…4aÿc†

where L]t is the Lie time-derivative operator de®ned in Eq. (1.1). In what follows we will also exploit the multiplicative split ^ be ˆ J2=3 e be ;

where Je ˆ …detbe †1=2

…5†

such that b^ e and Je represent the isochoric and volumetric part of the elastic deformation, respectively.

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1115

2.2. Elastic constitutive relation The starting point for the derivation of an elastic constitutive relation is the following strain energy function: h i  1 G tr…ln b^ e †2 el ˆ vol ‰Je Š ‡ iso b^ e ; where vol ˆ K…ln Je †2 ; iso ˆ 2 4

…6†

The decoupled form in the above strain energy function is based on the volumetric± isochoric split (5) for the elastic left Cauchy±Green tensor with corresponding bulk modulus K and shear modulus G, respectively. As noted by Steinmann et al. (1993) the use of the spatial, elastic logarithmic Hencky strains lnbe as a strain measure yields some advantage in the ensuing numerical implementation of the model. Next the Kirchho€ stresses are obtained according to  ˆ @be el be with the following result  ˆ vol ‡ iso and vol ˆ m g] ; m ˆ K ln Je iso ˆ G ln b^ e ˆ Gdev ln be ˆ dev :

…7†

A thermodynamic argument for the above result is given by Steinmann et al. (1993) but will not be repeated here for brevity of representation. 2.3. Two representations of the Gurson plastic potential 2.3.1. Quadratic representation A plastic potential for porous plastic materials is derived by Gurson (1977) from simple cell models. Modi®cations of this model have been proposed to improve the predictions at low void volume fractions by Tvergaard (1981) and to provide a better representation of ®nal void coalescence by Tvergaard and Needleman (1984). The original quadratic formulation for the Gurson plastic potential, which also serves as a yield function, often used in the literature is as follows: 1: 2: 3: 4:

2 3 dev ÿ '; where ˆ 2Y2

3q2 m ' ˆ 1 ‡ q3 f 2 ÿ 2q1 f  cosh ;  ˆ 2Y  f if f4f c fˆ fc ‡ k…f ÿ fc † if f > fc Y ˆ Y‰ev Š

…8†

Here m ˆ tr…†=3 is the Cauchy pressure, dev is the deviatoric part of the Cauchy stress tensor , and kk denotes the corresponding inner product norm. The state variable f for the micro-void material is the local micro-void volume fraction

1116

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

fˆ1ÿ

Vm VA

…9†

where VA is the elementary volume of the material, and Vm is the corresponding matrix volume. The parameter f  appearing in (8.3) is a function of the volume void fraction f. This modi®cation was suggested by Tvergaard and Needleman (1984) to provide a better representation of ®nal void coalescence for f > fc . The constants q1 ; q2 and q3 were introduced by Tveergard (1981) to improve the predictions of the original Gurson model at small void fractions. Material hardening is incorporated through Y which is a measure for the average ¯ow strength of the matrix material depending on the strain-like internal variable (the equivalent plastic strain) e . In what follows a combination of a saturation term and a linear hardening term will be used (see e.g. Simo and Hughes, 1998). Y ˆ Y0 ‡ Q‰e Š; where Q‰e Š ˆ q…1 ÿ exp…ÿbe †† ‡ He

…10†

where Y0 ; q; b; H are material parameters. 2.3.2. Linear representation The modi®ed linear formulation used in the sequel of this work is as follows: r q

dev 2

ÿ sign…'† ' YJ; where 1:  ˆ  3 3q2 m  2: ' ˆ 1 ‡ q3 f 2 ÿ 2q1 f  Cosh;  ˆ 2YJ 3: f  ˆ f ‡ …k ÿ 1†…f ÿ fc †T ‰f; fc ;
fc Š 4:

Y ˆ Y‰e Š

…11†

where  ˆ J is the Kirchho€ stress tensor and m ˆ tr…†=3 is the Kirchho€ pressure. The sign-function has been introduced in the above presentation in order to avoid the case ' < 0 for  ˆ 0, which is not admissible in the formulation (8). Then, the equivalence of the two formulations

(8) and (11) can easily be shown for the loading and unloading case 40, since dev is not negative and J and Y are always positive. The threshold function T ‰f; fc ;
fc Š has been introduced into the above formulation (11) in order to obtain a smooth transition of f  ‰ f Š near f ˆ fc and thus to avoid non-existence of the derivative at the critical stage f ˆ fc . A possible choice is the following function, 8 if x4xmin > <0 x ÿ x 1  min T ‰x; xmin ;
xŠ ˆ 1 ‡ sin ÿ if xmin < x < xmin ‡
x > 2
x :2 1 if xmin ‡
x4x …12† and this is illustrated in Fig. 1(a). Alternatively a Hermitian polynomial as introduced by Johansson et al. (1999) can be used. Note, that the above formulation (11.3) reduces to the original formulation (8.3) for
fc ! 0.

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1117

Numerical studies with the cosh function appearing in (8.2) revealed diculties (i.e. over¯ow) for large in-between values of  during the local iteration procedure (described in Section 3.3) although , which is strongly related to the mean-normal stress ratio m =Y, showed relatively small values (usually less than 30) for the ®nal solution points. Therefore, in order to obtain a more stable iteration scheme and still preserving the same results as in the original formulation the cosh-function in (8.2) is replaced by the following Cosh-function: ( cosh x if x4xc 1 …13† Coshx ˆ 2 cosh xc ‡ sinh xc …x ÿ xc † ‡ cosh xc …x ÿ xc † if x > xc 2 From Fig. 1(b) it can be seen that both formulations do coincide for values of jxj4jxc j. This ensures the equivalence of the proposed modi®cation with the original formulation, if the solution for jj of the local iteration is less than jxc j, which in turn implies, that the mean±normal stress ratio m =Y is less than 2xc =3q2 . Thus, e.g. for m =Y < 10, which is a realistic assumption, and q2 ˆ 2 it follows xc < 30. For larger values of jj > jxc j, which may occur during the local iteration procedure the original cosh-function is replaced by a quadratic term arising in the Cosh-function, thus preventing the over¯ow e€ect. Furthermore it is observed, that the Cosh-function shows continuity at the critical point xc up to the second derivative @2x Coshx, thus giving no diculty for the linearization analysis. 2.4. Rate equations In order to describe the change in the material during loading rate equations must be speci®ed. Following Simo and Miehe (1992) an associative ¯ow rule valid for isotropic material behavior is assumed: 1: 2: 3:

1 @ 1 ˆ ln ‡ t1; where ÿ L]t be
bÿ1 e ˆl 2 @ 3 dev

nˆ

dev r 1 ] 3 1 ÿ1 q1 q2 f  @ Cosh t ˆ tr…ÿ Lt bb
be † ˆ l q 2 2 '

…14†

where 1 is a second order unit tensor. The scalar l is the plastic multiplier obtained from the loading and unloading conditions l50; 40; l ˆ 0

…15†

A thermodynamic justi®cation of the above ¯ow rule is also given by Steinmann et al. (1993). Changes in the void volume fraction result from void nucleation as well as from growth of existing voids. The former can be either strain or stress induced, and following

1118

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

Fig. 1. Illustration of (a) Threshold function with smoothing property, (b) Modi®ed Cosh- function (with xc ˆ 5).

Chu and Needleman (1980) only strain induced void nucleation with a normal distribution will be regarded here. The void growth rate is assumed to be proportional to the plastic volume dilatation rate 1:

@t f ˆ @t fnucl ‡ @t fgrowth

2:

@t fnucl ˆ hnucl @t e ; where hnucl ˆ

3:

1 @t fgrowth ˆ …1 ÿ f†tr…ÿ L]t be
bÿ1 e † ˆ …1 ÿ f†t 2

  fN 1 e ÿ " N 2 p exp ÿ … † 2 sN SN 2 …16†

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1119

Furthermore, the e€ective plastic strain rate e is deduced from the equivalence of dissipated power expressions in the matrix and in the material related to the actual con®guration s  ' 2 1 1 ] 1 m t sign…'† ‡ …l † @t e Y…1 ÿ f † ˆ  : …ÿ Lt be
bÿ1 e † ? @ t ev ˆ 3 J 2 1ÿf JY

…17†

2.5. Summary of material parameters We are now in a position to summarize all material constants of the model, characterizing the inelastic behavior, which, apart from the elastic constants K and G, have to be calibrated based on experimental data  ˆ ‰Y0 ; q; b; H; q1 ; q2 ; q3 ; k; fc ; fN ; "N ; sN Št

…18†

In order to be physically meaningful the material parameters are restricted to lower and upper bounds ai ; bi , respectively. These constraints then de®ne the feasible domain K, such that  …19†  2 K  I Rnp ; K :ˆ  : ai 4i 4bi ; i ˆ 1; . . . ; np where np ˆ dim…† ˆ 12 denotes the number of material parameters. 3. Numerical implementation 3.1. Integration scheme In this section, the numerical integration of the evolution equations in the previous section is described over a ®nite time step
t ˆ n‡1t ÿ nt for given initial data n n n ÿ1 q; f; Cp , deformation gradient n‡1 F and Jacobian n‡1 J. Following Steinmann et al. (1993) the starting point for numerical integration of the ¯ow rule (14.1) is its representation in the material setting by use of the relations (4) 1 ÿ1 @ F ÿ @t Cÿ1 p
Cp ˆ lF 2 @

…20†

ÿ1 Then, applying an exponential map integrator to the expression @t Cÿ1 p ˆ ÿ2lF
ÿ1 @=@
F
Cp (Eterovic and Bathe, 1990; Weber and Anand, 1990) and an Euler backward rule to the rate equations (16) and (17) the following incremental objective integration algorithm is obtained

1120

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1: 2: 3:

@n‡1  n‡1 n ÿ1 †
F
Cp @ n‡1 f ˆn f ‡n‡1 hnucl …n‡1 e ÿ ne † ‡ …1 ÿ n‡1f †n‡1 t   2 ' 1 n‡1 m t n‡1 n sign…'† ‡ e ˆ e ‡
l 1 ÿn‡1 f JY 3 n‡1

n‡1 ÿ1 Cÿ1 F exp…ÿ2
l p ˆ

…21†

ÿ
where the relation exp Aÿ1
…†
A ˆ Aÿ1
exp…†
A has been used. From now on, if confusion is out of danger we will neglect the index n ‡ 1 referring to the actual time step. Then, the push-forward relation (4b) applied to (21.1) yields an update for the elastic left Cauchy±Green tensor be ˆ exp…ÿ2
l

@ tr †
b ; @

where

t btr ˆ F
n Cÿ1 p
F

…22†

Due to elastic isotropy be and  commute and due to plastic isotropy  and @=@ commute, and therefore be ;  and btr have identical principal axes. Then the volumetric and the deviatoric part of the spatial logarithmic Hencky strains follow as 1: 2:

ln Je ˆ ln Jtr ÿ t; where Jtr ˆ det…btr †1=2 ln b^ e ˆ ln b^ tr ÿ 2
ln; where b^ tr ˆ Jtrÿ2=3 btr

…23†

Consequently in view of (7) the Kirchho€ stresses are obtained from the relations 1:

 ˆ vol ‡  dev ;

2:

vol ˆ m g] ;

3:

dev ˆ dev;tr ÿ 2G
ln;

where

m ˆ mtr ÿ Kt;

mtr ˆ K ln Jtr

dev;tr ˆ Gdev ln btr

…24†

3.2. Spectral decomposition Upon using a spectral decomposition of the left elastic Cauchy±Green trial tensor btr and recalling, that due to isotropy btr and  commute, we have btr ˆ

3 3 X X 2 …ltr † m ?  ˆ A mA A A Aˆ1

…25a; b†

Aˆ1

2 tr Here …ltr A † and mA ; A ˆ 1; 2; 3 are the eigenvalues and eigenbasis of b , respectively, and A ; A ˆ 1; 2; 3 are the principal values of the Kirchho€ stresses which by use of the vector/matrix notations 2 3 2 3 2 3 2 3 1 1 1 ln ltr 1 1 6 7 6 7 dev 6 6 7 7 "tr :ˆ 4 ln ltr 1 5; I3 :ˆ I3 ÿ 1 1 …26† 2 5; :ˆ 4 2 5; 1 :ˆ 4 1 5; I3 :ˆ 4 3 1 1 3 ln ltr 3

are obtained from the relations

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1:

ˆ vol ‡ dev

2:

vol ˆ m 1;

3:

dev tr

dev ˆ dev;tr ÿ 2G
l; dev;tr ˆ 2GIdev
" ;  ˆ 3

dev



1121

where

m ˆ tr m ÿ Kt;

tr tr m ˆ K1
"

…27†

These sets of equations can be regarded as the counterpart of the relations (24) in the principal directions. Note, that the above structure for the principal Kirchho€ stresses is identical as in the geometric linear theory (see e.g. Simo, 1992), and therefore we can rewrite the representation (27) as 1:

ˆ C"el ;

2:

"el ˆ "tr ÿ
"pl ;

where 1
"pl ˆ
l ‡ t1; C ˆ K1 1 ‡ 2GIdev 3 3

…28†

3.3. Local iteration The four unknowns x ˆ ‰
l; m ; f; e Št appearing in (21) and (27) are obtained from a nonlinear system of equations r ˆ ‰r1 ; r2 ; r3 ; r4 Š ˆ 0, where r

2

dev;tr …Y ‡ Q‰e ŠJ† ‡ p‡ ‰
lŠ 1: r1 ˆ sign…'†

ÿ 2G
l ÿ 3 2: r2 ˆ m ÿ tr m ‡ Kt n 3: r3 ˆ f ÿ f ‡ hnucl …e ÿ ne † ‡ …1 ÿ f†t 0 s 1  ' 2 1  t m @
l A …29† sign…'† ‡ 4: r4 ˆ ne ÿ e ‡ 1ÿf JY 3 Here (29.1) expresses the yield constraint and (29.2) is obtained from the pressure relation (24.2). The conditions (29.3), (29.4) result from the update scheme (21.2), (21.3) for the void volume fraction and the equivalent plastic strain, respectively. In Eq. (29.1) also a penalty term ( 0 if x40 1 2 …30† p‡ ‰xŠ ˆ rp x if x > 0 2 with a large value rp has been introduced, in order to avoid negative values for
l (although this is neither a sucient nor necessary condition). The problem r ˆ 0 can be solved iteratively with a Newton method  ÿ1 x…k‡1† ˆ x…k† ÿ J…k† r…k† ; k ˆ 0; 1; 2; . . . …31† where the Jacobian

1122

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

Jˆ

@r @x

…32†

is required at each iteration step. The expressions for the elements of J can be extracted from the general variation r summarized in Appendix A, and therefore we will not elaborate on more details. 4. Equilibrium problem and associated derivatives 4.1. Weak formulation Denoting ' as the con®guration ®eld at time n‡1 t for given parameters  2 K and using the notation h; i for the L2 dual pairing relative to B0 of functions, vectors or tensor ®elds, the equilibrium problem as the classical weak form of momentum at time n‡1 t with spatial quantities reads Find ' : g‰'Š ˆ h : d i ÿ g ˆ 0 8u for given  2 K

…33†

Here the spatial rate of deformation tensor d induced by the virtual displacement u is de®ned in Eq. (1.3), and g :ˆ hB
ui ‡ hT
ui@ B designates the external part of the weak form for the case of dead loading with dual pairing h; i@ B on the boundary @ B and volume forces B and surface forces T . 4.2. General concept: directional derivative and sensitivity operator Before calculating derivatives of the weak form (33) it is useful to consider the dependencies of some quantities w.r.t. the con®guration ' and the material parameters . For example, from the de®nition of the deformation gradient F ˆ rX ' X;n‡1 t;  we can write F ˆ F‰'‰ŠŠ. Contrary, the plastic part of the deformation gradient at time n t is not dependent on the actual con®guration, and therefore we have n Fp ˆn Fp ‰Š. Both quantities de®ne the left Cauchy±Green trial tensor, and thus btr ˆ btr ‰'‰Š; Š ?  ˆ ‰'‰Š; Š

…34†

where the last relation for the Kirchho€ stresses is due to (24). The above relations motivate the de®nitions of the following operators for any (scalar, vector or tensor-valued) function w‰'‰Š; Š, which is dependent on the material parameters both implicitly via the con®guration ' and explicitly (Mahnken and Stein, 1997; Mahnken, 1999): Firstly, we introduce the standard directional derivative (Gateaux) operator @
w ˆ

dw d

u ˆ w‰'‰Š ‡ "
u; Š "ˆ0 d' d"

…35†

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1123

necessary for linearization of w. Secondly, upon using the notation v ˆ d'=d we de®ne a sensitivity operator 1:

@ w ˆ @' w ‡ @p w;

2:

@' w ˆ

3:

where

 n h io np dw d
v ˆ w '‰Š ‡ "vi ; "ˆ0 d' d"   dw…; † w‰'‰Š; i ‡
i Š ÿ w‰'‰Š; i Š np ˆ lim @p w ˆ
i !0 d
i iˆ1

…36†

which de®nes the total derivative @ w ˆ dw=d. Note, that the term @' w has the same structure as @
w and can thus be obtained from the results for linearization by simply exchanging
u with v . The second term @ w basically excludes the implicit dependence of  via the con®guration at the actual time (or load) step n‡1 t and will subsequently be called partial parameter derivative. Applying the above concept to the Kirchho€ stresses we obtain @
 ˆ L]
 ‡ 2sym…
lt
†; @  ˆ

L] 

‡

2sym…
lt †

‡

@p ;

where

L]
 ˆ 2

where

L] 

@ 1 [ [ : L g ˆ C : d
@g[ 2


@ 1 ˆ 2 [ : L[ g[ ˆ C: d @g 2

…37†

Here the spatial rate of deformation tensors d
and d induced by the linearization displacement
u and con®guration sensitivity v , respectively, are de®ned in Eq. (1). Furthermore, C:ˆ 2@=@g[ is the fourth order algorithmic spatial operator, and as alluded to above, the partial parameter derivative @p  excludes the implicit dependence of  via the con®guration '. Finally, upon using the relations @
d ˆ ÿsym…g[
l
l
† @ d ˆ ÿsym…g[
l
l †

…38†

we obtain the associated derivatives of the weak form as @
g‰@Š ˆ h…C : d
† : d ‡ l
 : d i @ g‰'Š ˆ h…C : d † : d ‡ l  : d ‡ @p  : d i

…39†

The ®rst term is used for computation of an increment
u within a Newton iteration step for iterative solution of the weak form, and the second term is used for computation of the con®guration sensitivity v within an iteration step for minimization of some least-squares functional. The next task consists in determination of the algorithmic spatial tangent operator p C and the partial parameter derivative of Kirchho€ stresses @ , consistent with the integration scheme of Section 3.

1124

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

4.3. Spatial algorithmic tangent operator The spatial algorithmic tangent operator C de®ned in Eq. (37) is obtained from the representation (25b) for the Kirchho€ stress tensor. Taking into account the dependency A ‰"tr ; JŠ the following result is obtained: C

ˆ

3 X 3 X d A Aˆ1 Bˆ1

tr mA mB ‡

d"B

3 3 X dmA X d A m A g[ 2 A ‡ J [ dJ dg Aˆ1 Aˆ1

…40†

An equivalent expression with the dependency A ‰"tr Š, thus yielding the ®rst two terms on the right hand side, is derived by Simo (1992), where also the result for dmA =dg[ is presented. The third term in the above result occurs due to the additional dependency of A on J. The quantities d A =d"tr B and d A =dJ appearing in Eq. (40) are obtained by straightforward di€erentiation of A in Eq. (27). Using the vector notation introduced in Section 3.2 the result is ! d 2G
l dev …2G†2
l d
l d
l

ˆ 2G 1 ÿ

dev tr I3 ‡ dev;tr   ‡ 1 d"tr ÿ 2G d"tr d"tr

d m d"tr d d m d
l 1 ÿ 2G  ˆ dJ dJ dJ ‡1

…41†

Here, the quantities d
l=d"tr ; d m =d"tr ; d
l=dJ; d m =dJ are contained in the following result: we consider the local problem (29) as an implicit function and conclude        tr  tr  dr dr @r @r @r dx dx ˆ ‡ ˆ0 r " ; J; x " ; J ˆ 0 ? ; ; ; d"tr dJ @"tr dJ @x d"tr dJ     …42† dx dx ÿ1 @r @r ; ; ˆ ÿJ ? d"tr dJ @"tr @J   where J is the Jacobian (32). The expressions for @r=@"tr ; @r=@J can be extracted from the general expression for the variation r in the Appendix, and therefore we will not elaborate on more details.

4.4. Partial parameter derivative of Kirchho€ stresses Using the representation (25b) the partial parameter derivative of Kirchho€ stresses appearing in (37)2 consists of two parts (see e.g. Mahnken and Stein, 1997):

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137 3 3 @p  X @p A X @p m A ‡ ˆ mA

A @ @ @ Aˆ1 Aˆ1

1125

…43†

The expressions necessary for determination of @p =@ are summarized in Tables 1 and 2, respectively. In Table 1 those terms are recalled from Mahnken and Stein (1997), which are valid for any model formulated in principal directions. Table 2 contains the expressions resulting from the speci®c model formulation in Section 2. For brevity of presentation quantities relating to the nucleation part of void volume fraction with associated parameters fN ; "N ; sN have been omitted, such that the number of material parameters occurring in the representation (19) is np ˆ 9. From the results of a pre-processing part in Tables 1 and 2, respectively, it can be seen that calculation of @p =@ also involves expressions for the parameter derivative of state variables dn e =d; dn f=d; dn Cÿ1 p =d at the previous time step. Therefore, after having solved the linear problem (39)2 for v , it becomes necessary to determine dn‡1 e =d; dn‡1 f=d; dn‡1 Cÿ1 p =d at the actual time step in a post-processing part, in order to make them available for the next time step. As already mentioned, all expressions of Table 1 are derived and explained in Mahnken and Stein (1997). In what follows, we will brie¯y comment on the expressions of Table 2. Starting point is the representation (28.2) for the principal elastic stretches, which yields the following parameter sensitivity: @"el @"tr @
l @ 1 @t ˆ ‡

‡
l ‡ 1

@ @ @ @ 3 @

…44†

The determination of @
l=@ and @ m =@ is contained in the following result: we consider the local problem (29) as an implicit function and derive   dr @ r @r dx dx @ r ˆ ‡ ˆ0? ˆ ÿJÿ1 r ; x;n e ;n f; "tr ˆ 0 ? d @ @x d d @

…45†

for the total derivative. Here J is the Jacobian (32) and the notation @ r @r @r dn ev @r dn f @r d"tr ‡ n ˆ ‡ n ‡ tr @ @ @ ev d @ f d @" d

…46†

has been used. Then, for the preprocessing step @p A =@ appearing in Eq. (43) is obtained from Eq. (28.1) @p @p "el ˆC @ @

…47†

where @p "el =@ is included in the result (44). More details necessary to obtain @p A =@ are summarized in the preprocessing part of Table 2, where it should be remarked, that ei denotes the ith unit vector.

1126

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

Table 1 Large strain problems formulated in principal directions: Partial parameter sensitivity for the Kirchho€ stresses and the parameter sensitivity of right Cauchy Green plastic strain (material independent part) a. Partial sensitivity for Kirchho€ stresses (pre-processing) Input: nCÿ1 p t @p btr ˆ L] btr ˆ F
@n Cÿ1 p
F 1 1 p tr @p "tr A ˆ ÿ tr
2 mA : @ b ; A ˆ 1; 2; 3 2 l A

for @p A : see material dependant part in Table 2 @p mA ˆ @btr mA : @p btr ; A ˆ 1; 2; 3 (for @btr mA see Miehe 1993) 3

3

Aˆ1

Aˆ1

@p  ˆ  @p A mA ‡  A @p mA b. Parameter sensitivity for the internal strain-like variables (post-processing). Input: n Cÿ1 p ; v @v ÿ1 ? l ˆ @ F
F ? @ J ˆ Jtr…l † @ F ˆ @X tr ] tr t @ b ˆ L  b ‡ 2sym…l
btr †, where L]  btr ˆ F
@n Cÿ1 p
F 1 1 tr @ "tr A ˆ ÿ tr
2 mA : @ b 2 l A

for @ "elA : see material dependent part in Table 2 el2 el @ lel2 A ˆ 2lA @ "A ; A ˆ 1; 2; 3

@ mA ˆ @btr mA : @ btr ; A ˆ 1; 2; 3 (for @btr mA see Miehe, 1993) 3 ÿ
el2 @ be ˆ  @ lel2 A mA ‡ lA @ mA Aˆ1   ÿ1 n ÿ1 ‡ Fÿ1
@ be
Fÿt @ Cÿ1 p ˆ ÿ2sym F
@ F
Cp

For determination of dn‡1 ev =d; dn‡1 f=d; dn‡1 Cÿ1 p =d in Table 2 in the post-processing step ®rstly d
l=d and df=d are calculated from Eq. (45). For determination el of dn‡1 Cÿ1 p =d in Table 1 d" =d has to be computed (see Mahnken and Stein, 1997, Proposition 5.8). The details of the post-processing step are also summarized in Table 2. 5. Numerical examples 5.1. Local iteration problem In this section the convergence behavior of the Newton iteration scheme (31) is compared for the quadratic formulation (8) and the linear formulation (11) of the Gurson plastic potential. To this end a (strain driven) local constitutive problem in the framework of a geometrically linear plane strain theory with sub-incrementation is considered. The following (arbitrary chosen) total strains n‡1 " and initial conditions for n "pl ;n f and n ev are assumed: 2

3 2 3 0:20 0 6 7 6 7 0:21 n‡1 ÿ5 n 7n 607 n " ˆ 0:156 4 ÿ0:10 5; "pl ˆ 4 0 5; f ˆ 2:10 ; ev ˆ 0:1724 ÿ0:26 0

…48†

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1127

Table 2 Large strain multiplicative von Mises elastoplasticity formulated in principal directions: algorithmic constants for the partial parameter sensitivity of the stress variables (pre-processing) and the parameter sensitivity of internal strain-like variables (post-processing)a Input for pre-processing: @n e ; @n f, replace @ "tr by @p "tr from Table 1, part (a) set @ J ˆ 0 Input for post-processing: @n e ; @n f; @ "tr ; @ J from Table 1, part (b) Trial stresses:

  ÿ
dev;tr @ tr ˆ Ktr @tr ˆ 2Gdev @ "tr m " ; @

Normal vector:

  dev;tr 1

dev @  ˆ

dev;tr I3 ÿ   @



Material parameters, Yield function, Threshold function

@ q1 ˆ e5 ; @ q2 ˆ e6 ; @ q3 ˆ e7 ; @ k ˆ e8 ; @ fc ˆ e9 3 2 3 2 0 1 7 6 0 6 qe exp…ÿbe † 7 7 6 6   7 7 6 7 6 0 7 6 6 1 ÿ exp…ÿbe † 7 7 6 7 6 7 6 0 7 6 7 6 e 7 6 7 6 7 6 0 7 6 7 6 @ @ Y ˆ 6 0 T ˆ 7 6  7 7 6 7 6 0 7 6 0 7 6 7 6 7 6 7 6 0 7 6 0 7 6 7 6 7 6 7 6 0 7 6  0 5 4 4 1  f ÿ fc   1 5 cos ÿ ÿ 0
fc 2 2
fc if fc < e < fc ‡
fc else @ T ˆ 0 Further quantities:

3@ q2 m 3q2 m 3q2 m ÿ @ Y ÿ @ J 2 2YJ 2Y J 2YJ2 @ f ˆ @ k… f ÿ fc †T ‰ f; fc ;
fc Š ÿ …k ÿ 1†@ fc T ‰ f; fc ;
fc Š ‡ …k ÿ 1†… f ÿ fc †@ T ‰ f; fc ;
fc Š @  ˆ



‡ 2q3 f  @ f  ÿ 2@ q1 f  Cosh ÿ 2q1 @ f  Cosh ÿ 2q1 f  @ Cosh@  r r 1 3 1 3 1  @ t ˆ ÿt @ ' ‡
l p @ q1 q2 f @ Cosh ‡
l p q1 @ q2 f  @ Cosh 2' 2 ' 2 ' r r 3 1 3 1 ‡
l p q1 q2 @ f  @ Cosh ‡
l p q1 q2 f  @2 Cosh@  2 ' 2 '

@  ' ˆ @  q3 f

2

(continued overpage)

1128

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

Table 2 (continued) Residual r:

@  r1 ˆ v : @ dev;tr ÿ

r r r 2p 2 1 2p 'J@ Y ÿ 'Y@ J p YJ@ ' ÿ 3 32 ' 3

@  r2 ˆ @ m ÿ @ tr m ‡ K@ t @  r3 ˆ @n f ‡ …1 ÿ f†@ t @  r4 ˆ

@n e

! r 1 2 1 m @ t m t m t
l ÿ 2 @ Y ÿ 2 @ J ‡ p @ ' ‡ 1ÿf 32 ' JY JY JY

Solve linear system

@ x ˆ ÿJÿ1 @ÿ  r ? @ t ˆ @ t ‡ @x t@ x Principal elastic stretches

1 @ "el ˆ @ "tr ‡ v @
l ‡
l @ v ‡ 1 @ t 3 If pre-processing then

@p ˆ C@ "el else

n‡1 Store@n‡1 f  e ; @ 

a For brevity of presentation quantities relating to the nucleation part of void volume fraction have been omitted.

The material parameters are chosen as E ˆ 206 000 Mpa,  ˆ 0:3, Y0 ˆ 360:0 Mpa, H ˆ 100 MPa, b ˆ 3:9, q ˆ 436, q1 ˆ 1:5, q2 ˆ 0:5, q3 ˆ 1, i.e. no nucleation part of void volume fraction and no ®nal void coalescence has been considered. The load factor 0.15 appearing in Eq. (48) is divided into 100 sub-increments thus yielding to the results for the von Mises stress and the void volume fraction as shown in Fig. 2. It is noteworthy that identical results for both formulations, i.e. the quadratic formulation (8) and the linear formulation (11) of the Gurson plastic potential are obtained. However, the iteration behavior was di€erent as shown in Table 3: from here the superior convergence behavior for the linear version of the Gurson yield potential can be observed. 5.2. Axisymmetric tensile specimen In this section it is the object to simulate the necking behavior of an axisymmetric tensile bar with geometry as shown in Fig. 3. The material of the specimen is a ferritic steel of the German designation 22 NiMoCr3 7. Five experimental tensile tests were

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1129

Fig. 2. Local iteration problem: von Mises stress and void volume fraction versus load factor. Table 3 Local iteration problem: norm of residual within a Newton iteration for quadratic and linear formulation of the Gurson modela ite 1 2 3 4 5 6 7 8 9 10 11 a

isub=1 quadratic 0.47443 E+05 0.11864 E+05 0.82653 E+05 0.24931 E+04 0.10835 E+03 0.18436 E+00 0.53575 Eÿ06 0.72051 Eÿ11

isub=70 linear 0.55721 E+04 0.27911 Eÿ01 0.78388 Eÿ08 0.55932 Eÿ13

quadratic 0.47750 E+05 0.13672 E+05 0.77185 E+05 0.63410 E+04 0.76053 E+03 0.91718 E+01 0.13842 Eÿ02 0.35815 Eÿ10 0.94389 Eÿ11

isub=100 linear 0.77715 E+04 0.78922 E+01 0.89139 Eÿ03 0.11705 Eÿ08 0.12859 Eÿ12

quadratic 0.48076 E+05 0.30917 E+05 0.17504 E+06 0.11244 E+05 0.79973 E+04 0.30156 E+04 0.40150 E+03 0.94590 E+01 0.56262 Eÿ02 0.20089 Eÿ08 0.67127 Eÿ11

linear 0.95548 E+04 0.87453 E+02 0.36931 E+00 0.27535 Eÿ02 0.23930 Eÿ07 0.24004 Eÿ12

Ite and isub denote the iteration step and the number of sub-increment, respectively.

performed at a temperature of T ˆ 0 C with quasi static displacement controlled loading in axial direction. All specimens showed necking near the middle of the specimen length. Fig. 5 shows experimental results for the load vs. elongation curve of one specimen, which failed in a sudden but completely ductile manner. T is the change of length with respect to a measuring length of L0 =2 ˆ 12:5 mm and in Fig. 6, experimental results for the elongation vs. radius change data of all ®ve (interrupted) tests are shown.

1130

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

Fig. 3. Axisymmetric tensile specimen: Geometry and discretization.

The numerical simulation of the specimen is performed with quadratic nine-node axisymmetric ®nite elements. Due to the symmetry of the structure, only one quarter is discretized as shown in Fig. 3. A small imperfection of
R ˆ 0:00025R0 is implied where R0 ˆ 30 mm is the radius of the specimen, in order to trigger the necking in the middle of the specimen. Integration is performed with the algorithm described in Section 3 in N ˆ 124 load steps. The underlying experimental data for the identi®cation process consist of the necking displacements uj in the middle of the specimen at nd =®ve observation states. Furthermore the total tensile loads Fj at nf ˆ 19 observation states is taken into account. Then the following objective function of least-squares type is examined nf ˆ19 d ˆ5 1 nX 1 X 2 wd …uj ÿ u j † ‡ wf …Fj ÿ F j †2 ! min f…† ˆ  2 jˆ1 2 jˆ1

…49†

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1131

Table 4 Axisymmetric tensile specimen: Starting and obtained values for the material parameters of a ferritic steel

H (N/mm2) q (N/mm2) b fN fc

Starting

Solution

550.0 180.0 20.0 0.004 0.0005

394.3 158.2 22.6 0.00054 0.00040

Fig. 4. Axisymmetric tensile specimen: Value of weighted least-squares functional versus number of iterations.

The values for uj and fj are given in mm and N, respectively, and therefore the weighting factors were chosen as wd ˆ 20 and wf ˆ 1  10ÿ5 . Young's modulus, Poisson's ratio and the initial ¯ow stress were predetermined as E ˆ 196000 N/mm2,  ˆ 0:3 and Y0 ˆ 471:0 N/mm2, respectively. Pre-calculations (not reported here) showed that it is dicult to ®nd a stable and unique set of material parameters characterizing the material. This diculty arises from the incomplete data, and therefore it seems that additional tests are necessary in order to obtain more reliable values for the material parameters, Mahnken and Stein (1997). The following parameters were chosen as constant "N ˆ 0:25; SN ˆ 0:1; f0 ˆ 0:; k ˆ 3:0; q1 ˆ 1:0; q2 ˆ 0:5;
fc ˆ 0:1
10ÿ3 ; q3 ˆ 1. Then ®ve parameters H; q; b; fN ; fc are varied in the optimization process with starting value according to the second

1132

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

Fig. 5. Axisymmetric tensile specimen: Total load versus elongation for experiment and simulation.

Fig. 6. Axisymmetric tensile specimen: Necking displacement versus elongation for experiment and simulation.

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1133

Fig. 7. Axisymmetric tensile specimen: Deformed con®guration at the end of loading.

row in Table 4. For minimization of the least-squares functional (49) an algorithm due to Bertsekas (1982) is used, where the iteration matrix is obtained from the update rule for the BFGS formula. Further details of the algorithm are explained in e.g. Mahnken and Stein (1996a,b). The results of the optimization process for the material parameters are given in the third row of Table 4. These results have been obtained after 34 iteration steps. The corresponding value of the least-squares functional (49) versus the number of iterations is shown in Fig. 4. In Fig. 5 the total load versus the elongation is depicted for both simulation and experiment after the optimization process, thus revealing a very good agreement. In Fig. 7 the deformed con®guration at the end of loading is shown. Finally, in Fig. 8 the equivalent plastic strain and the void volume fraction f at the end of loading are depicted, thus showing the highly nonuniform behavior within the specimen. Furthermore the maximum value of damage is obtained in the center of the specimen thus initiating the cup-cone fracture e€ect, which has also been observed in the experiment.

1134

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

Fig. 8. Axisymmetric tensile specimen: Equivalent plastic strain and void volume fraction at the end of loading thus revealing the highly non-uniform behaviour within the specimen.

6. Summary and conclusions In this contribution various aspects for the ®nite element implementation of the Gurson model were considered. In particular a linear representation for the plastic potential has been used, thus obtaining a better convergence property in the local iteration procedure compared to the original quadratic representation. In order to prevent over¯ow during the local iteration procedure for in-between values, the original cosh-function has been replaced by a modi®ed Cosh-function, which however yields the same solution as the original formulation for moderate values of the mean-normal stress ratio m =Y. Furthermore a threshold function has been introduced which has a smoothing e€ect on the additional variable f introduced by Tvergaard and Needleman (1984). A further important issue was directed to the sensitivity analysis consistent with the underlying exponential map integration scheme. This enables to apply a gradient based optimization scheme for minimizing a least-squares functional as an objective function for parameter identi®cation. The algorithm has been applied in order to obtain material parameters for a ferritic steel based on experimental data for an axisymmetric tensile bar, thus yielding a very good agreement between the simulated data and the experimental data.

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1135

Considering the reliability of the material parameters obtained in this way by simply using macroscopic testing data, it should be kept in mind, that the results may become highly unstable or even non-unique. In Mahnken and Stein (1996a), where this subject is discussed in more detail, it has been pointed out, that possible reasons for the phenomena of instability or non-uniqueness are twofold: ®rstly a de®ciency of the model due to linearities within the model formulation may be present, or secondly a de®ciency of the experiment may occur, in the sense that not all mechanisms intended by the model are properly activated during the testing. For example in the second example of this paper the parameters q1 ; q2 ; q3 were chosen a priori, thus clearly leading to a dependency of some of the other parameters. In this respect we cannot claim uniqueness for the solution of material parameters. Therefore future work should be directed to a more comprehensive study, by combining this purely macroscopic approach with metallurgical examinations in order to take into account more physical insight for some of the parameters and thus leading to more stable solutions.

Acknowledgements The author would like to thank Professor Brocks and Dr.-Ing. Bernauer from the GKSS Forschungszentrum Geesthacht, Germany, for provision of the experimental data of the axisymmetric tensile specimen.

Appendix. Variations of signi®cant quantities In the sequel variations of signi®cant quantities necessary for evaluation of the   Jacobian (32), the quantities @r=@"tr ; @r=@J appearing on the r.h.s. in Eq. (42) necessary for linearization of the weak form and the quantities @ r=@ appearing in (46) and thus necessary for the sensitivity analysis are derived. Variation of Y ˆ Y0 ‡ q…1 ÿ exp…ÿbev †† ‡ Hev Y ˆ Y0 ‡ q…1 ÿ exp…ÿbev †† ‡ qe exp…ÿbe †b ‡ He ‡ …qb exp…ÿbe † ‡ H†e Variation of  ˆ  ˆ

…A1†

3q2 m 2YJ

3q2 m 3q2 m 3q2 m 3q2 m ‡ ÿ Y ÿ J 2 2YJ 2YJ 2Y J 2YJ2

…A2†

1136

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

Variation of f  ˆ f ‡ …k ÿ 1†…f ÿ fc †T ‰ f; fc ;
fc Š f  ˆ f ‡ k…f ÿ fc †T ‰ f; fc ;
fc Š ‡ …k ÿ 1†…f ÿ fc †T ‰ f; fc ;
fc Š ‡ …k ÿ 1†…f ÿ fc †T ‰ f; fc ;
fc Š

…A3†

Variation of ' ˆ 1 ‡ q3 f  2 ÿ 2q1 f  Cosh ' ˆ q3 f 2 ‡ 2q3 f  f  ÿ 2q1 f  Cosh ÿ 2q1 f  Cosh ÿ 2q1 f  @ Cosh r 3 1 p q1 q2 f  @ Cosh Variation of t ˆ
l 2 j'j r 3 1 sign…'† ' t ˆ 
l q  q1 q2 f  @ Cosh ÿ t 2 ' 2'

…A4†

r r 3 1 3 1  q1 q2 f @ Cosh ‡
l q q1 q2 f  @ Cosh ‡
l q 2 ' 2 '

r r 3 1 3 1  ‡
l q  q1 q2 f @ Cosh ‡
l q  q1 q2 f  @2 Cosh 2 ' 2 ' r

2p

Variation of r1 ˆ sign…'† dev;tr ÿ 2G
l ÿ j'jYJ 3   r r q dev;tr :  dev;tr 2 1 2 ' YJ

q  ÿ 2G
l ÿ YJ' ÿ sign…'† r1 ˆ sign…'†

dev;tr 3 2 ' 3

r q 2 sign…'† ' YJ ÿ 3 Variation of r2 ˆ m ÿ tr m ‡ Kt r2 ˆ  m ÿ  tr m ‡ Kt n

…A5†

…A6†

…A7† n

Variation of r3 ˆ f ÿ f ‡ hnucl …e ÿ e † ‡ …1 ÿ f †t …A8† r3 ˆ n f ÿ f ‡ hnucl …e ÿn e † ‡ hnucl …e ÿ n e † ‡ ft ‡ …1 ÿ f †t r   1 2j'j m t sign…'† ‡ Variation of r4 ˆ ne ÿ e ‡
l 1ÿf 3 JY s   2 ' 1 m t  sign…'† ‡ f‡
l r4 ˆ n e ÿ e ‡ JY 3 …1 ÿ f †2 0 1 s r  2' 1 B 2 1 m t m t m t m t C q ' ‡ ‡ ÿ sign…'† ‡
l Y ÿ 2 JA @
l 1ÿf 3 2 ' JY JY JY 2 J Y 3 …A9†

R. Mahnken / International Journal of Plasticity 15 (1999) 1111±1137

1137

References Bennani, B., Picart, P., Oudin, J., 1993. Some basic ®nite element analysis of microvoid nucleation, growth and coalescence. Eng. Comp. 10, 409±421. Bertsekas, D.P., 1982. Projected Newton methods for optimization problems with simple constraints. SIAM J. Con. Opt. 20 (2), 221±246. Chu, C., Needleman, A., 1980. Void nucleation e€ects in biaxially stretched sheets. J. Eng. Mater. Technol. 102, 249. Eterovic, A.L., Bathe, K.J., 1990. A hyperelastic based large strain elasto-plastic constitutive formulation with combined isotropic-kinematic hardening using the logarithmic stress and strain measures. Int. J. Num. Meth. Eng. 30, 1099±1114. Gurson, A.L., 1977. Continuum theory of ductile rupture by void nucleation and growth Ð I. Yield criteria and ¯ow rules for porous ductile media. Engng. Materials Technol. 99, 2±15. Hancock, J.W., Mackenzi, A.C., 1976. On the mechanisms of ductile failure in high-strength steels subjected to multi-axial stress-states. J. Mech. Phy. Solids 24, 147±169. Johansson, M., Mahnken, R., Runesson, K., 1999. Ecient integration for generalized viscoplasticity coupled to damage. Int. J. Num. Meth. Eng. 44, 1727±1747. Mahnken, R. 1999. A comprehensive study of a multiplicative elastoplasticity model coupled to damage including parameter identi®cation. Comp. & Structures, in press. Mahnken, R., Stein, E., 1996a. Parameter identi®cation of viscoplastic models based on analytical derivatives of a least-squares functional and stability investigations. Int. J. Plast. 12 (4), 451±479. Mahnken, R., Stein, E., 1996b. A uni®ed approach for parameter identi®cation of inelastic material models in the frame of the ®nite element method. Comp. Meths. Appl. Mech. Eng. 136, 225±258. Mahnken, R., Stein, E., 1997. Parameter identi®cation for ®nite deformation elasto-plasticity in principal directions. Comp. Meths. Appl. Mech. Eng. 147, 17±39. Miehe, C., 1993. Computation of isotropic tensor functions. Communications in Applied Numerical Methods 9, 889±896. Simo, J.C., 1992. Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the in®nitesimal theory. Comp. Meths. Appl. Mech. Eng. 99, 61±112. Simo, J.C., Hughes, T.J.R., 1998. Computational Inelasticity. Interdisciplinary Applied Mathematics, Vol. 7, Mechanics and Materials. Springer±Verlag. Simo, J.C., Miehe, C., 1992. Associative coupled thermoplasticity at ®nite strains: formulation, numerical analysis and implementation. Comp. Meths. Appl. Mech. Eng. 98, 41±104. Steglich, D., Brocks, W., 1997. Micromechanical modelling of the behaviour of ductile materials including particles. Comp. Mat. Sc. 9, 7±17. Steinmann, P., Miehe, Ch., Stein, E., 1993. Comparison of di€erent ®nite deformation inelastic damage models within multiplicative plasticity for ductile materials. Int. J. Comp. Mech. 13, 458±474. Thomason, P.F., 1990. Ductile Fracture of Metals. Pergamon Press, Oxford. Tvergaard, V., 1981. In¯uence of voids on shear band instabilities under plane strain conditions. Int. J. Fract. 17, 389±407. Tvergaard, V., Needleman, A., 1984. Analysis of the cup-cone fracture in a round tensile bar. Acta Metallurgica 32, 157±169. Weber, G., Anand, L., 1990. Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic±viscoplastic solids. Comp. Meths. Appl. Mech. Eng. 79, 173±202.