- Email: [email protected]
Semi-implicit finite strain constitutive integration of porous plasticity models
Finite Elements in Analysis and Design 104 (2015) 41–55
Contents lists available at ScienceDirect
Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel
Semi-implicit finite strain constitutive integration of porous plasticity models P. Areias a,d,n,1, T. Rabczuk b, J. César de Sá c a
Department of Physics, University of Évora, Colégio Luís António Verney, Rua Romão Ramalho, 59, 7002-554 Évora, Portugal Institute of Structural Mechanics, Bauhaus-University Weimar, Marienstraße 15, 99423 Weimar, Germany c Mechanical Engineering Department, Faculty of Engineering, University of Porto, Rua Dr. Roberto Frias, s/n, 4200-465 Porto, Portugal d ICIST, Instituto Superior Técnico, Lisboa, Portugal b
art ic l e i nf o
a b s t r a c t
Article history: Received 11 December 2014 Received in revised form 16 May 2015 Accepted 19 May 2015 Available online 12 June 2015
Two porous plasticity models, Rousselier and Gurson–Tvergaard–Needleman (GTN), are integrated with a new semi-implicit integration algorithm for finite strain plasticity. It consists of using relative Green– Lagrange during the iteration process and incremental frame updating corresponding to a polar decomposition. Lowdin's method of orthogonalization is adopted to ensure incremental frameinvariance. In addition, a smooth replacement of the complementarity condition is used. Since porous models are known to be difficult to integrate due to the combined effect of void fraction growth, stress and effective plastic strain evolution, we perform a complete assessment of our semi-implicit algorithm. Semi-implicit algorithms take advantage of different evolution rates to enhance the robustness in difficult to converge problems. A detailed description of the constitutive algorithm is performed, with the key components comprehensively exposed. In addition to the fully detailed constitutive algorithms, we use mixed finite strain elements based on Arnold's MINI formulation. This formulation passes the inf– sup test and allows a direct application with porous models. Isoerror maps for two common initial stress states are shown. In addition, we extensively test the two models with established benchmarks. Specifically, the cylindrical tension test as well as the butterfly shear specimen are adopted for validation. A 3D tension test is used to investigate mesh dependence and the effect of a length scale. Results show remarkable robustness. & 2015 Elsevier B.V. All rights reserved.
Keywords: Constitutive integration Porous plasticity Finite strains Semi-implicit Löwdin's method
1. Introduction For common metals at room temperature, void nucleation and particle debonding, void growth and coalescence are the mechanisms for crack growth from the process region [10]. Detailed multiple-scale analysis of these mechanisms require robust polycrystalline plasticity and representation of dislocation dynamics and grain boundary barriers (Roters et al. [30] present a comprehensive review). Therefore, the so-called phenomenological models (e.g. [18]) are considered a cost-effective tool for predicting damage and fracture in metals when overall quantities, like maximum load and dissipated energy, are of interest. The now well-established GTN model is based on the original yield surface by Gurson [16] and was subsequently modified by Tvergaard and Needleman [34,35]. The void growth part of the GTN model is based on the micromechanics of the ductile process. Phenomenological models formulated based on thermodynamical n
Corresponding author. 1 Researcher ID: A-8849-2013 http://www.researcherid.com/rid/A-8849-2013 Web-of-Science search: areias, p*
http://dx.doi.org/10.1016/j.finel.2015.05.005 0168-874X/& 2015 Elsevier B.V. All rights reserved.
principles have also been proposed such as the Lemaitre model [18] and its subsequent modifications. Among this class of constitutive models, the Rousselier model [31,4] is interesting from a computational perspective as it is simpler and requires the specification of fewer parameters than the GTN model, while possessing the ability to closely match the predictions of the latter by determination of its parameters. It adheres to the idea that the occurrence of localization, crack initiation and propagation emerges as a direct consequence of strain softening due to void growth. Hence, unlike micro-mechanics models of porous metal plasticity, the Rousselier model does not directly model the mechanism of coalescence. However, two properties are worth mentioning: 1. The Rousselier model predicts void fraction growth in pure shear. 2. The Rousselier yield surface has an isochoric vertex causing a non-zero plastic volumetric term component. As a more complete alternative to the Rousselier model, it was recently (see [22]) created an extended version of the GTN model, which includes the shear void fraction as an independent history variable. We here focus on variants of the GTN and Rousselier
42
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
models to assess the isoerror maps and robustness with a new semi-implicit integration algorithm. To compare the GTN and the Rousselier models, we use the following procedure 1. The GTN model uses the same coalescence approach as the classical Tvergaard–Needleman technique (cf. [35]). 2. Both void fraction and damage variables are integrated in closed-form using the Lambert W and the ERF functions. Concerning the semi-implicit integration method with porous plasticity models, some considerations are worth mentioning: 1. Due to their complementarity form (see, e.g. [20]), finite strain plasticity problems often exhibit convergence difficulties for large strain values and severe sensitivity to step size [3]. We have two measures to attenuate these difficulties: the use of a smoothed complementarity condition and the removal of iterative rotation matrices from the constitutive laws in finite strains. 2. It is computationally convenient that both hyperelastic and finite-strain elasto-plastic laws are implemented in a unique algorithm and applicable to any discretization scheme. 3. Classical assumed-strain elements typically do not directly provide the deformation gradient (e.g. [9]), necessary for many constitutive formulations. An estimated deformation gradient can be calculated from the polar decomposition if an approximate rotation matrix is available. Since a constitutive frame is adopted, the rotation matrix is obtained from this frame in two distinct configurations. Compatible frames use Löwdin's method [19]. 4. Since the seminal contributions of Weber and Anand [36] and Simo [32], Kirchhoff stress tensors (i.e. τ ¼ J σ ) are frequently employed in the yield functions. This is a computational convenience, as commonly adopted elasto-plastic and hyperelastic models are often quasi-incompressible. This is not the case of porous plasticity (cf. [4]). We therefore avoid that approach. 5. Semi-implicit formulations, where certain quantities are fixed in the flow vector (but not the flow vector itself), cf. [23,8] can lead to substantial savings in constitutive integration. We further extend the semi-implicit algorithm presented in [3] to achieve very large time steps. Consistent linearization of integrated form of objective rates is intricate and computationally expensive, making it a possible candidate for the explicit integration part of the semi-implicit scheme. However, rigid body motions are exactly accounted for.
t
t ¼σ v u¼u
u
on Γ ta
2.1. Integration algorithm: equilibrium and relative Green–Lagrange strain Cauchy equations of equilibrium for a rotated reference configuration are obtained from the corresponding spatial equilibrium (derivations for the latter are shown in Ogden [27]). Using standard notation (cf. [27]) we write the spatial equilibrium equations as ð1Þ
with the Cauchy tensor components σij (i; j ¼ 1; 2; 3). In (1) i is the direction index and j is the facet index. The components of the body force vector are identified as bi. In (1), coordinates xaj are the spatial, or deformed, coordinates of a given point under consideration. In addition, the following natural and essential boundary conditions
Γta is the ð2Þ
on Γ a
u
ð3Þ
where t is the known stress vector on Γ v is the outer normal and u is the known displacement field on Γua. It is assumed that (1) and (2)– (3) are satisfied for a time parameter t a A ½0; T with T being the total time of observation and for a point with position xa A Ωa belonging to the deformed position domain at the time of analysis. Equilibrium configuration corresponds to the domain Ωa and is identified by the subscript a. In tensor notation, Eq. (1) can be presented as t a,
∇ σT þ b ¼ 0
ð4Þ
with ∇ ¼ ∂=∂xa being the spatial gradient operator. After multiplication _ integration in the deformed configuration Ωa by the velocity field u, and application of integration by parts component-wise, we obtain the _ int is the internal and W _ ext is the external following power form (W power): Z Z Z σ : L dΩa ¼ b u_ dΩa þ t u_ dΓ a ð5Þ Ωa Ωa Γ ta |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl {zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl } _ _ W int W ext _ where L is the velocity gradient: L ¼ ∂u=∂x a . L can be decomposed into a symmetric part (D, called strain rate) and a skew-symmetric part (W , called vorticity), L ¼ D þW . We also introduce the Jacobian J ¼ det F and F is the deformation gradient. The use of equivalence _ int and external W _ ext power for a continuum has between internal W been used to obtain objective rates of the Cauchy stress σ . In the corotational case, we can rotate the Kirchhoff stress and the strain rate to obtain Z _ int ¼ RT0a τR0a : RT0a DR0a dΩ0 ð6Þ W Ω0
From each element nodal arrangement, we use three directions, where e 1a , e 2a and e 3a to define R0a as R0a ¼ ½e 1a ; e 2a ; e 3a
ð7Þ
where each column is a unit vector and e ia e ja ¼ δij . Section 2.3 describes the calculation of R0a . We then use the following convention for writing tensors: the bold symbol indicates tensor components. Therefore, the frame where the tensor components are defined is identified as a superscript. For example, the secondorder tensor S ab with the equilibrium configuration Ωa and reference configuration Ωb has the following components in frame c:S cab . If a global frame 0 is adopted, component transformation between frame 0 and frame c follows: S cab ¼ Rc0 S 0ab RTc0
2. Constitutive integration in finite strains
∂σ ij þ bi ¼ 0 ∂xaj
hold on each part of the boundary Γ a ¼ Γ a [ Γ a where natural boundary and Γua is the essential boundary:
ð8Þ
Note that this approach, although related, does not correspond to the rotated Algorithm of Areias and Belytschko (cf. [2]). Generalizing (8), change of basis2 from d to c results in the following transformation: S cab ¼ Rcd S dab RTcd
ð9Þ
In addition, change in reference configuration reads S 0ac
¼ Rbc S 0ab RTbc
ð10Þ
Using (10) and transformation (9) for coinciding reference configuration and frame, S cac ¼ S bab
2
This corresponds to a change of frame, denoted by a superscript.
ð11Þ
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
From the relation between the time derivative right Cauchy– Green tensor C and the strain rate D, we obtain C_ ¼ 2F T DF
ð12Þ
Simo [33, p. 281, Eq. 8.1.13] derived the following one-step scheme for integrating the strain rate in frame 0, producing the following “strain” from integration of D: h i e0ab;β ffi 12 F Tbβ F Taβ F aβ I F bβ ð13Þ where β is an arbitrary configuration. Using the previous power conjugacy (6), the internal power can be written in the initial configuration Ω0 as Z _ int ¼ S 0a0 : D0a0 dΩ0 W Ω0
where τ and Using a hypoelastic relation, stress updating is performed as follows: S 0a0
¼ RT0a 0a R0a
D0a0
0 ¼ RT0a da R0a .
S 0a0 ¼ S 0b0 þ ΔS a ðRβ0 e0ab;β RTβ0 Þ
S bab ¼ S bbb þ ΔS a ðRβ0 e0ab;β RTβ0 Þ
where Sba0 ¼ Voigt½S ba0 . We note that, omitting indices a, b, and 0, V S ðFÞ can be written as 2
F 211 6 2 6 F 12 6 6 2 6F V S ðFÞ ¼ 6 13 6 F 11 F 12 6 6 4 F 11 F 13 F 12 F 13
F 221
F 231
2F 21 F 11
2F 31 F 11
2F 31 F 21
F 222 F 223
F 232 F 233
2F 22 F 12
2F 32 F 12
2F 32 F 22
2F 23 F 13
2F 33 F 13
2F 33 F 23
F 21 F 22 F 21 F 23
F 31 F 32 F 31 F 33
F 21 F 12 þF 11 F 22 F 21 F 13 þF 11 F 23
F 31 F 12 þ F 11 F 32 F 31 F 13 þ F 11 F 33
F 22 F 23
F 32 F 33
F 22 F 13 þF 12 F 23
F 32 F 13 þ F 12 F 33
S bab ¼ S bbb þ ΔS a ðebab Þ
ð14Þ
with Simo integration scheme reducing to h i ebab ffi 12 Rb0 F Tab F ab I RTb0
ð15Þ
Ωb in (14) as
1 b 1 1 S ¼ S bbb þ ΔS a ðebab Þ J b0 ab J b0 J b0 |fflffl{zfflffl}S ⋆b |fflffl{zfflffl}σ b
ð16Þ
b
Local constitutive quantities, as well as quadratic strain incre Local frames a and b, and therefore the rotation tensor are explicitly integrated.
Given F bab , ebab (both in frame b), R0b and R0a (both in frame 0) and Tb and Ta Recovered from storage F bb0 , σ 0b , ebb0 (σ b is stored in frame 0 for purpose of representation) Represent Cauchy stress in frame b σ b ¼ V RT σ 0 b
Relevant Jacobian determinants
ð17Þ
then we obtain the Cauchy stress in configuration as 1 a 1 J S ¼ S b ¼ b0 S ⋆b J a0 aa J a0 ab J a0 ab
Ωa and frame a ð18Þ
If ΔS a ðebab Þ is homogeneous of degree one, then it follows that b ~b S ⋆b ab ¼ σ b þ ΔS a e ab
0b
b
J b0 ¼ det F bb0 J ab ¼ det F bab Rab ¼ RT0a R0b
Relative rotation and total deformation gradient update
1 b e J b0 ab
ð19Þ
For hyperelastic models, we require the total strain, which is the strain from configuration Ωa with respect to configuration Ω0. Using frame b, we have the following update formula for the Green–Lagrange strain: b b eba0 ¼ ebb0 þ F bT b0 eab F b0
ð20Þ S ⋆b ab
We can observe that is work-conjugate to eba0 in configuration Ωb. Global Cauchy stress can be obtained as σ 0b ¼ R0b σ bb RT0b . In frame a, eaa0 reads eaa0 ¼ RTba eba0 Rba
Thermal effect in the relative strain eb○ ¼ eb þ αðT a T b ÞI ab ab b○ eba0 ¼ ebb0 þ V E F bT b0 eab
Total strain update
b○ Corrected relative strain for UL e~ ab ¼ J1 eb○ b0 ab b ~ b○ and Qab along with (UL) Determine S⋆b ab ¼ σ b þ ΔS a e ab
S a ∂ΔS a 1 ∂Q ab ab , ∂T a , J ~ b○ and ∂Q sensitivity quantities Cab ¼ J1 ∂Δ~ b○ ∂T a b0 ∂e ab b0 ∂e ab (TL) Determine Sba0 eba0 and Qab, along with sensitivity ∂Sb
quantities Ca0 ¼ ∂ea0 b ,
∂Sba0 ∂Q ab ∂T a , ∂eb a0
and
∂Q ab ∂T a
b b 1 S⋆b ab ¼ J b0 V S F b0 Sa0 (TL) Determine relative sensitivities C ¼ 1 V F b C V ðF bT Þ a0 E ab b0 b0 J b0 S b ∂S⋆b ∂Sa0 b 1 ab ∂T a ¼ J b0 V S F b0 ∂T a ∂Q ab ab ¼ ∂Q V E F bT b0 ∂eb ∂eb a0
(TL)Determine relative stresses
where b e~ ab ¼
S
F ba0 ¼ F bab F bb0
1 b S J b0 ab
σ aa ¼ S ⋆a aa ¼
7 7 7 7 7 7 F 31 F 22 þ F 21 F 32 7 7 7 F 31 F 23 þ F 21 F 33 5 F 32 F 23 þ F 22 F 33
and, for the strain quantities, V E ðFÞ ¼ V TS ðFÞ. The full algorithm is given in Algorithm 1. In summary:
Defining S ⋆b ab ¼
3
Algorithm 1. General relative Lagrangian formulation (Voigt notation adopted).
In the case β ¼ b, the result is
ab
b For symmetric tensors, such as S ⋆b ab and eab , it is preferable to use the Voigt notation. Upright bold symbols denote Voigt form of symmetric tensors: 1 b b ⋆b S⋆b V S F b0 Sa0 ð22Þ ab ¼ Voigt½S ab ¼ J b0
ments are implicitly integrated.
where ΔS a is a function of the strain e0ab;β . This term ΔS a is called the constitutive part of the stress. However, since S 0a0 ¼ S bab and S 0b0 ¼ S bbb , we have
We now insert the Jacobian for configuration
43
ð21Þ
a0
ab
α
∂Sab ∂Sab ∂T a ∂T a þ C ab I eaa0 ¼ V E ðRab Þeba0 F aa0 ¼ Rab F ba0 RTab
Update temperature sensitivity Determine strain in frame a
Determine deformation gradient in frame a Determine Cauchy stress in frame 0 σ 0 ¼ 1 V S ðR0a ÞS⋆b ab a J ab
Store F aa0 , σ 0a , eaa0 Return Qab, S⋆b ab , C ab ,
∂S⋆b ab ∂Q ab ∂T a , ∂eb ab
and
∂Q ab ∂T a
44
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
2.2. Local constitutive system
We determine the relevant constitutive quantities by solving a general nonlinear constitutive system. Omitting the frame superscript, it reduces to the root finding for the following nonlinear system: ! zfflfflfflfflffl}|fflfflfflfflffl{Governing φ Sab ; Q ab ; eKab ; eab ; Ta ; χa ¼0 |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}Unknown |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl}Known |{z}Internalvar: ð23Þ where χ a is the internal variable vector. For the present application, it is defined as ( )
χa ¼
εp fa
ð24Þ
where εp is the effective plastic strain and fa is the void fraction variable. Newton iteration for (23) provides, in general, the solution for the constitutive unknowns Sab , Qab, and χ a . For the present application, system (23) has Q ab ¼ 0, T a ¼ 0, eKab ¼ ∅. A smoothed version of (23) is introduced by using the Chen– Mangasarian replacement functions, which depend on a parameter Error such that (cf. [3]) lim φError ¼ φ
Error-0
ð25Þ
Newton iteration on (23), using φError as a replacement provides the following scheme: ( ) h ∂φ i ΔSab ∂φError Error ð26Þ ¼ φError ∂Sab ∂χ a |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}J Δχ a After achieving a solution using the Newton scheme (26), we calculate the sensitivities from the following equation: ( ) dSab ∂φError deab ð27Þ ¼ J 1 dχ a ∂eab Specifically (with further details being provided in [3]), the present elasto-plastic system is therefore given as 8 9 1 ΔS ab nΔγ > > < eab Clinear =
ð28Þ φError ðS ab ; eab ; χ a Þ ¼ μ⋆ Δγ μ⋆ Δγ þ ϕ Error > > : ; Δχ a Δγωðχ a Þ where Δγ is the plastic multiplier increment, μ⋆ is a dimensional parameter described in [3] and the smooth ramp function of Chen and Mangasarian [14] is used for the third equation, which depends on Error. Function ϕ in (28) is the yield function. In addition, the flow vector n is also required for the flow rule, and it is defined as n¼
∂ϕ ∂Sab
ð29Þ
Fig. 1 shows the effect of Error in the satisfaction of the complementarity condition (see also [3]). In (28), function ωðχ a Þ is the internal variable evolution function, such that _ ðSab ; χ a Þ χ_ a ¼ γω
ð30Þ
The motivation for a unique treatment for all yield functions is that any particular yield function is inserted by means of ϕ, n and derivatives of n. This is accomplished with Mathematica and the AceGen add-on [29,17]. Functions hi and hiError are respectively given by hi ¼
þ j j 2
ð31Þ
Fig. 1. Replacement of μ⋆ Δγ μ⋆ Δγ þ ϕ by μ⋆ Δγ 〈μ⋆ Δγ þ ϕ〉Error as a function of a Error parameter (μ⋆ ¼ 1 is depicted).
hiError ¼ þ
Error lnð2Þ ln 1 þ exp lnð2Þ Error
ð32Þ
where hiError is an approximation of hi. 2.3. Löwdin frame With the objective of uniquely defining the element frame compatible with the concept of polar decomposition, use is made of the symmetric orthogonalization formula by Per-Olov Löwdin (cf. [19]). With that purpose, we resort to the singular value decomposition (SVD). Given a real 3 3 matrix J, there exists a decomposition with the form J ¼ U ΣV T
ð33Þ
where U is a 3 3 unit matrix, U U ¼ I, Σ is a 3 3 diagonal matrix and V is a 3 3 unit matrix with V T V ¼ I. Columns of U are normalized eigenvectors of JJ T and columns of V are normalized eigenvectors of J T J. The diagonal matrix Σ contains the singular values of J in its principal diagonal. Singular values are nonnegative reals. Columns of U are the left singular vectors and columns of V are the right singular vectors. The non-zero singular values of J are the square roots of the non-zero eigenvalues of J T J. In the present case columns of J form a basis. Löwdin's formula provides the following orthogonal frame: T
R ¼ UV T
ð34Þ
It was proved by Carlson and Keller [12] that R obtained by this formula (34) is the proper orthogonal matrix whose columns have the minimum Euclidean distance to the corresponding columns of J. A more concise formula for R is R ¼ JV Σ
1
VT
ð35Þ
and, since Σ is a diagonal matrix, we have the following index notation: 1 Rij ¼ pffiffiffiffiJ ik V kl V jl
λl
ð36Þ
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
45
where λl, l ¼ 1; …; n, are the eigenvalues of J T J and columns of V are the corresponding eigenvectors. In our case, and for a tetrahedron, we have the following specialization: 1
R0a ¼ J a V a Σa V Ta
ð37Þ
where J a is the element Jacobian for configuration Ωa, V a are the eigenvectors of J Ta J a . For the configuration Ωb, we have R0b ¼ J b V b σ b 1 V Tb
ð38Þ
The following notation is used for the Jacobian matrices: Ja ¼
∂xa ∂ξ
ð39Þ
Jb ¼
∂xb ∂ξ
ð40Þ
where ξ is the vector of parent-domain coordinates. The motivation for the use of R0a and R0b to obtain a rotation can be summarized as follows: In problems where local frames are required (anisotropic laws, fracture, shells), rotation can be obtained by multiplication Rab ¼ RT0a R0b . For uniformity, we use the Löwdin local frames in all problems. 2.4. The extended GTN model and its integration The version of the Gurson–Tvergaard–Needleman model introduced by Reis et al. [28] is rewritten here in non-dimensional form: " !# 3ST I6 Tdev S q2 ST I3 2 ϕGTN y; S; f ⋆ ¼ 1 þ q3 f ⋆ 2q1 f ⋆ cosh 2y 2y2 ð41Þ The effective void fraction f ⋆ is a function of the void fraction f . We use the nucleation part of the void fraction following the work of Chu and Needleman [15] using the normal distribution: " # f 1 εp εN 2 f_ N ¼ pNffiffiffiffiffiffi exp ε_ p ð42Þ 2 SN SN 2π where fN, SN and εN are constitutive properties for the nucleation (N) part of the void fraction law. The classical void growth term of Gurson (cf. [16]) is adopted: f G ¼ 1 expð εv Þ
ð43Þ
where εv is the trace of deformation. In addition, the shear part of Xue and Nahshon (cf. [37,24,25]) are combined to obtain the following term: 0 1 2 3q6 6 1=3 B det ½S dev C 1=3 2 fS ¼ ð44Þ @1 54 3 Af εpa 2 π ST I6 Tdev S Including the initial void fraction, f0, then the total f follows the following closed-form laws from the above classical relations: f εN ε εp pffiffiffi ERF N pffiffiffi f ¼ f 0 þ N ERF ð45Þ þ 1 expð εv Þ 2 SN 2 SN 2 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} f 0 1 2 3q6 6 1=3 B det ½S dev C 1=3 2 þ @1 54 3 Af εpa 2 π ST I6 Tdev S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}f S
ð46Þ
The additional property q6 controls the void growth in shear. A depiction of (41) with is shown in Fig. 2. In contrast with the Rousselier model the GTN model is symmetric with respect to the
Fig. 2. Extended GTN yield function (y ¼ q1 ¼ q2 ¼ q3 ¼ 1) with p ¼ ST I3 =3 and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 T 2S I6 Tdev S.
q¼
trace of S. The effect of increasing f can be observed as a uniform reduction in radius of the yield surface. The error function ERFðÞ results from closed-form integration of the Gaussian curve and is available in several open-source libraries and is defined as Z 2 2 ERFðÞ ¼ pffiffiffiffi e t dt ð47Þ
π
0
Two distinct types of effective plastic strain are adopted,3 and εp. The latter being given as
ε_ p ¼
γ_ ST n ð1 f Þy ab
εv
ð48Þ
Integration of constitutive quantities follow the unconditionally stable backward-Euler scheme:
εpa ¼ εpb þ
Δγ
ST n ð1 f Þy ab
ð49Þ
In addition, we have the trace of the flow vector, which is used to obtain εv: ε_ v ¼ γ_ nT I3 ð50Þ It is straightforward to obtain the critical ff, corresponding to collapse, which is given by: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q1 þ q21 q3 ff ¼ ð51Þ q3 The effective void fraction f ⋆ is calculated according to the additional property fc, using the Tvergaard–Needleman modification (cf. [35]): 8 f; f of c > > > > < f f f c f f c ; f c rf rf m ð52Þ f⋆ ¼ fcþ f m f c > > > > :ff; f 4f m The effective void fraction (52) is used in the yield function, and is not updated as is the case of f. 3 No power-consistent closed-form solution exists for the effective plastic strain in the classical Gurson model.
46
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
Fig. 3. Rousselier yield function (y ¼ σ 1 ¼ D ¼ 1) with p ¼ ST I3 =3 and q ¼
2.5. The Rousselier model and its integration
Table 1 Constitutive properties for the two porous models (Rousselier and extended GTN).
The Rousselier [31] model for porous metals is characterized by the following yield function in non-dimensional form: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " # T 3 T σ1 f ⋆ ST I3 2S I6 Tdev S þ D exp ϕR y; S; f ⋆ ¼ 1 ð53Þ y 3 1 f ⋆ σ 1 1f ⋆ y in which f ⋆ is the effective void fraction in configuration Ωa and y is the yield stress. In addition to the hardening characteristics reflected in y, the only additional properties are σ1, D and the initial porosity f0. The property σ1 can be estimated from the original yield stress, y0 ¼ yð0Þ and the ultimate tensile strength, ymax ¼ maxεp y as σ estimate ¼ 13 y0 þ ymax ð54Þ 1 Note that ymax must be provided by the actual hardening curve in the porous case. The effective void fraction f ⋆ depends on the corresponding constitutive void fraction f. A depiction of (53) is shown in Fig. 3. Note that the Rousselier model is asymmetric with respect to the trace of S. The effect of increasing f can be observed. Void fraction follows from volume change and is calculated exactly as in the previous section. Fracture value ff coincides with the zero stress condition fulfilling the condition ϕ ¼ 0: ff ¼
y Dσ 1
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 T 2S I6 Tdev S.
ð55Þ
Property
Rousselier
Extended GTN
Elastic properties Hardening law Initial void fraction Void fraction properties
E, ν y yðεp Þ f0 σ1 D
Nucleation
fN SN εN q6 fc fm
E, ν y yðεp Þ f0 q1 q2 q3 fN SN εN q6 fc fm
Shear Coalescence
where Wð) is the Lambert W function. The effective plastic strain rate is given by this power equivalence (see also [4]):
γ_ ε_ p ¼ STab n
ð58Þ
y
and
ε_ v ¼
γ_
ð1 f Þy
STab I3
nT I3
ð59Þ
Integration of constitutive quantities follow the unconditionally stable backward-Euler scheme:
Δγ
It can be observed that D indirectly controls the void fraction rate. We can force the growth near the critical value of void fraction to model coalescence as in the GTN model, but with a slight difference (see [35]): 8 f; f of c > > > > < f f f c f f c ; f c r f rf m ð56Þ f⋆ ¼ fcþ f m f c > > > > :ff; f 4f m
Time stepping is adapted so that Δεp ¼ εpn þ 1 εpn is kept below 5%. This measure is necessary for accuracy reasons in problems with localization.
with fm, fc and ff being additional properties. An mentioned, initial void fraction f0 must also be provided. The equivalent plastic strain increment is computed based on the power equivalence relation: 1 y ε_ p ð57Þ STab ε_ p ¼ ð1 f Þ y 3σ 1 W Df ⋆ exp 3 3σ 1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}y
For Error ¼ 1 10 3 we test the two porous plasticity models, along with the classical von-Mises yield criterion. Properties are shown in Table 1 for the two presented models. Von-Mises is a particularization. Isoerror maps are shown in Fig. 4 for the three yield functions. We use the properties shown in Table 2 for the runs of the isoerror maps. Plane strain condition is assumed (4 direct stress components
εpa ¼ εpb þ
y
STab n
ð60Þ
2.6. Summary of properties for porosity models and isoerror maps
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
47
Fig. 4. Isoerror maps for von-Mises, Rousselier and GTN yield functions (a single step is used).
and 1 shear stress component). For each strain increment a single step is used. The following two conclusions are drawn:
yield function is not smooth for all stress states. It is also acknowledgedly difficult to integrate (cf. [4]).
Very small integration errors are observed for the von-Mises
yield function. This is well known since in a simpler case an exact solution exists. The Rousselier yield function exhibits a larger error than the GTN yield function. This is now relatively well known since this
2.7. Fitting of Rousselier properties to the GTN model To estimate the properties for the Rousselier model from known properties of the GTN model, we can obtain D, σ1 from q1, q2 and q3
48
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
Table 2 Isoerror maps: relevant properties for two prototype materials (consistent units). Property
Rousselier
Property
Extended GTN
E ν y f0
72 400 0.33 352
E ν y f0
72 400 0.33 352
2
2 10 500 106
σ1 D
2 10 1 0.042 0.1 0.2 1 0.06 0.10
fN SN εN q6 fc fm
q1 q2
2 10 1.5 1.0
q3 fN SN εN q6 fc fm
2.25 0.042 0.1 0.2 1 0.06 0.10
[5], see also Bathe [7] and Cao [11], is based on a two-field formulation where:
Pressure is linearly interpolated using the corner nodes. An internal shape function, named a bubble, enriches the velocity or displacement fields.
2
Cauchy stress is calculated from the constitutive stress S⋆b ab as
σ aa ¼
1 ⋆b S J ab ab
ð66Þ
, Cauchy pressure is obtained as T S⋆b I3 ab pa ¼ 3J ab
ð67Þ
We can therefore write S⋆b ab in Voigt form as a sum of deviatoric and pressure terms:
Initial stress, Case I 352 σxx σyy 352 Initial stress, Case II σxx 352 σyy 0
⋆b S⋆b ab ¼ J ab pa I3 þ Tdev Sab
Table 3 Properties for 2024-T351 aluminum (consistent units are adopted). Rousselier
Extended GTN
E ν y
72 400 0.33 352 þ 668ε0:63302 p
E ν y
72 400 0.33 352 þ 668ε0:63302 p
f0 σ1 D fN SN εN q6 fc fm
0.0 174.51 1.07291 0.042 0.1 0.2 1 0.06 0.99
f0 q1 q2 q3 fN SN εN q6 fc fm
0.0 1.5 1.0 2.25 0.042 0.1 0.2 1 0.06 0.99
where Tdev is the 22 13 3 6 1 2 6 3 3 6 6 1 13 6 Tdev ¼ 6 3 60 0 6 6 0 40 0
0
following sparse matrix: 3 13 0 0 0 7 13 0 0 0 7 7 7 2 0 0 07 3 7 0 1 0 07 7 7 0 0 1 05 0 0 0 1
ð68Þ
ð69Þ
In terms of power balance, we use the following relation, where ⋆b ~ It corresponds to a S~ ab depends on the independent pressure p. classical two-field variational principle: 0 T 1 Z Z S⋆b I3 C ab ⋆b T b B _ ext ð70Þ ¼W S~ ab e_ ab dΩb þ @J ab p~ þ Ap~_ dΩb 3 Ωb Ωb |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}W_ int where the relative Jacobian Jab is used to ensure correct volume calculation. We note that the product J ab p~ cannot be used as an unknown field. In (70), we have the following quantities:
at a given void-fraction level fc. Using least-squares, we solve the following least-squares minimization problem Cf. Table 3: 2 2 i 1h ð61Þ qGTN qRousselier þ pGTN pRousselier min σ 1 ;D 2
⋆b ~ 3 þTdev S⋆b S~ ab ¼ J ab pI ab
ð71Þ
Discretization follows the standard MINI formulation: uðξÞ ¼
5 X
N K ξ uK
ð72Þ
K¼1
where
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qGTN ¼ y 1 þf c q3 2q1
ð62Þ
qRousselier ¼ ð1 f c Þ y σ 1 Df c 2
2y 1 þ f c q3 pGTN ¼ arccosh q2 2f c q1 pRousselier ¼ 3ð1 f c Þσ 1 log
ð63Þ
! ð64Þ
y
Df c σ 1
ð65Þ
Independent pressure p~ is interpolated using the corner nodes: ~ ξÞ ¼ pð
4 X
NK ξ p~ K
ð73Þ
K¼1
N1 ðξÞ ¼ 1 ξ1 ξ2 ξ3
ð74aÞ
N 2 ðξ Þ ¼ ξ 2
ð74bÞ
N 3 ðξ Þ ¼ ξ 3
ð74cÞ
N 4 ðξ Þ ¼ ξ 1
ð74dÞ
3. Element technology for quasi-incompressible problems 3.1. Tetrahedron Quasi-incompressible constitutive laws severely reduce the performance of displacement-based elements. Usually, mixed elements are a solution for avoiding locking in quasiincompressible problems. The low-order MINI element by Arnold
The bubble function, for tetrahedra, is given by N5 ðξÞ ¼ ξ1 ξ2 ξ3 1 ξ1 ξ2 ξ3
ð5Þ
For the calculation of the stiffness matrix, the variation of (70) is required. Not all quantities are determined by hand-derivation, and we use Mathematica [29] with the AceGen (cf. [17]) add-on to
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
49
Fig. 5. Axisymmetric representation of tension specimen and relevant contour plots for the GTN and Rousselier models.
calculate some derivatives. Using (70), we obtain Z ⋆b T ⋆b T b _ int ¼ dS~ ab e_ ab þ S~ ab debab dΩb dW Ωb
T 3 dS⋆b I3 7 ab 6 þ 4dJ ab p~ þ J ab dp~ þ 5p~_ dΩb 3 Ωb Z
⋆b
~ 3 ~ 3 J ab dpI dS~ ab ¼ Tdev Cab debab dJ ab pI ð75Þ
2
Specifically, the terms and dJ ab are determined by AceGen. In terms of rotation tensor, we use the following process to obtain e 1 , e 2 and e 3 .
ð76Þ
3.2. Dimensional reduction
ð77Þ
Dimensional reduction for plane-strain requires careful consideration of out-of-plane stress, which is in general non-null. We therefore explicitly enforce the plane-strain condition by setting 13 and 23 components of both stress and strain as zero, but
where the following notation was used: b dS⋆b ab ¼ C ab deab
ð78Þ
debab
50
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
Fig. 6. Void nucleation and growth in the GTN and Rousselier models (10 000 triangles).
h ⋆b i including the out-of-plane stress S~ ab
33
as active. Therefore, for
plane strain, we only omit out-of-plane shear components of stress and strain. In the axisymmetric case, we use reduced integration, which is sufficient to remove locking in the nearincompressible case. 4. Numerical tests All examples are run in the SIMPLAS framework, cf. [1] with the present formulations. MINI elements were implemented to comply with the Löwdin frames and the relative Green–Lagrange strains. Internal degrees-of-freedom of mixed elements are condensed out. 4.1. Assessment of mesh size sensitivity using classical tension specimens An axisymmetric idealization of a simple tension test with smooth geometry (cf. [6] for a complete testing program) is made, with the relevant dimensions shown in Fig. 5. The corresponding properties are shown in Table 1. Also shown in Fig. 5 are the contour plots for the effective plastic strain and void fraction in both models. Rousselier yield function produces earlier failure due to a steeper void fraction growth in the early nucleation stages. We can observe that no hourglass or mesh instabilities are visible and contour plots are relatively smooth. We compare the void fraction evolution for the two models, separating the nucleation and void growth terms, in Fig. 6. It can
Fig. 7. Aluminum specimen, GTN model, mesh size sensitivity: effective plastic strain, void fraction and reaction–displacement results. Step size Δv ¼ 0:09 mm.
be observed that the growth term in Rousselier model4 is much quicker to develop and tends to localize strain in a sharper form (see also Fig. 5). Evolution of void fraction and effective plastic strain for the GTN yield function is shown in Fig. 7. Reaction–displacement results are also shown in the same figure. We obtain two conclusions:
Effective plastic strain and void fraction evolution are close (however not perfectly close in the case of the effective plastic strain) to the ones reported by Malcher et al. [21]. Reaction– displacement results are very similar to the reported experimental curve (cf. [21])
4 Note that properties have been fitted so that intersections in p- and q-axes for both models coincide at f ¼ f c .
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
51
Fig. 9. Butterfly specimen: relevant dimensions, coarse mesh (4483 nodes and 16 473 MINI tetrahedra) and fine mesh (6351 nodes and 24 132 MINI tetrahedra). One symmetry plane (zz) is taken into account.
Fig. 8. Aluminum specimen, Rousselier model, step size sensitivity: effective plastic strain, void fraction and reaction–displacement results. Mesh with 10 000 triangles and 5151 nodes.
Some mesh sensitivity is observed, however it can be considered mild and acceptable. Concerning the Rousselier model, evolution of void fraction and effective plastic strain is shown in Fig. 8, with mild mesh sensitivity. In terms of reaction–displacement results, results are very robust, with near-absence of mesh sensitivity. 4.2. Butterfly test This test, introduced by Bai and Wierzbicki [6] is required to assess the performance in shear of both yield functions. A description of this test is shown in Fig. 9 with the corresponding geometrical data. Both [21,22] perform extensive testing with this test to calibrate the constitutive properties. In this work, and for this test, Aluminum 2024-T351 is used and both mesh size and time step size sensitivities are inspected.
For the coarse mesh (due to size constraints), we show a sequence of relevant contour plots in Fig. 10. Since the MINI mixed element technology is employed, severe deformations are possible without volumetric locking. In addition, we note that a shear band is formed at the specimen weak region. Effect of step size is shown in Fig. 11 for the GTN model and in Fig. 12 for the Rousselier model. Only mild dependence on step size is observed. A complete depiction of the dependence on mesh size is shown in Fig. 13 where again only slight dependence is observed. We note that, if stronger localization is required, then the mesh size has a serious effect in the results and regularization will be required. Mesh size effect on the effective plastic strain and void fraction evolution at the monitored node are shown in Fig. 14. In this case some effect is noted, which is not sufficient to detract from the overall very robust results. 4.3. Mesh dependence and regularization: 3D tension test To inspect the mesh dependence with full localization (until collapse), we use the geometry from Norris tension test [26] (see also [32]). Fig. 15 shows the tension test specimen with effective
52
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
Fig. 10. Butterfly specimen (coarse mesh): sequence of contour plots for the GTN model (effective plastic strain and void fraction).
Fig. 11. GTN model: comparison with experimental results and numerical results by Malcher et al. [22]. Effect of step size.
plastic strain and void fraction contour plots. Two meshes are used: a coarse one with 3175 nodes and a finer one with 5409 nodes (the latter is used in the contour plots). A fixed displacement step of Δw ¼ 0:02 is adopted for this problem. We use the constitutive properties described in Table 4 which produce full localization. Regularization is performed by introducing a lengthscale lc ratio in Eq. (42) with εp -lc εp , and Eq. (43) with εv -lc εv
Fig. 12. Rousselier model: comparison with experimental results. Effect of step size.
for the GTN model and similarly for the Rousselier model. This is sufficient to produce nearly mesh-insensitive results. The lengthscale ratio is defined as follows: lc ¼
min le lf
ð79Þ
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
53
Fig. 13. Effect of mesh-size (von-Mises, GTN and Rousselier).
Fig. 15. 3D tension test: relevant geometrical data and contour plot for the GTN model.
Table 4 Properties for the 3D tension test (consistent units are adopted). Rousselier
Fig. 14. (a) Effective plastic strain and (b) void fraction evolution, comparison with Malcher et al. [21,22] and the model of Bai and Wierzbicki (cf. [6]).
where le is the edge length and lf is the fixed length size, a constitutive property (given in Table 4). Relatively smooth effective plastic strain and void fraction contour plots are obtained. In terms of displacement/reaction results for the GTN problem, we test the unregularized (lc ¼ 1) and regularized cases with the two meshes (see Fig. 16 for the GTN case). The length scale is sufficient to obtain insensitive results. For the Rousselier model, see Fig. 16, the same conclusion is obtained.
Extended GTN
lf E ν y
0.35 72 400 0.33 352 þ 668ε0:63302 p
lf E ν y
0.35 72 400 0.33 352 þ 668ε0:63302 p
f0 σ1 D fN SN εN q6 fc fm
0.002 20.4 9.53 0.08 0.01 0.05 0.5 0.06 0.2
f0 q1 q2 q3 fN SN εN q6 fc fm
0.002 1.5 1.0 2.25 0.08 0.01 0.05 0.5 0.06 0.2
5. Conclusions We presented a general semi-implicit constitutive algorithm and specialized it for the integration of porous plasticity yield
54
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
shown in the 3D tension test, where both models were tested with two meshes. In conclusion, for the numerical tests, we assessed the mesh-size independence and the choice of yield function in the 2D tension test, butterfly test and 3D tension test. Very robust results are obtained and the algorithm allows severe strains with accompanying mesh distortions in shear without failure of convergence. In particular, for relatively low values of void fraction, results are in general insensitive to mesh and step sizes without regularization. For fully localized problems, a simple length scale is sufficient to remove mesh dependence near the collapse.
References
Fig. 16. Tension test: comparison of two meshes in non-regularized and regularized cases. (a) GTN model (both non-regularized and regularized cases). (b) Rousselier model (regularized cases).
functions. New contributions and distinctive features of the proposed algorithm are:
The combination of relative Green–Lagrange strains with the
exact corotational method. Quadratic strains are used during iteration whereas a local frame is updated after convergence of Newton iteration. The use of Löwdin algorithm to calculate orthogonal frames based on the Jacobians of configurations Ωb and Ωa. This corresponds to a definition of orthogonal frames consistent with the polar decomposition. The use of a local constitutive system which is smooth and nonlinear. This is achieved in two steps by replacing the complementarity condition for a non-smooth equation and subsequently by the Chen–Mangasarian replacement function (cf. [13,14]). Local Newton iteration is performed without concerns with the complementarity condition. The GTN and Rousselier constitutive models are explored in detail, with the determination of intersections with axes of pressure–shear diagram (p–q). We include isoerror maps and the explicit determination of the void fraction in both models. The use of finite strain versions of the MINI element (cf. [5]) in triangles and tetrahedra to avoid locking due to quasiincompressibility in the early stages of void fraction evolution. Specializations and the second variation of the variational principle, required to calculate the stiffness matrix. The use of a length-scale in the stress measures and effective plastic strain present in the void growth laws, which was found to be sufficient to produce mesh-insensitive results. This is
[1] P. Areias. Simplas. 〈http://www.simplas-software.com〉. [2] P. Areias, T. Belytschko, Analysis of finite strain anisotropic elastoplastic fracture in thin plates and shells, J. Aerosp. Eng. 19 (4) (2006) 259–270. [3] P. Areias, D. Dias-da-Costa, E.B. Pires, J. Infante Barbosa, A new semi-implicit formulation for multiple-surface flow rules in multiplicative plasticity, Comput. Mech. 49 (2012) 545–564. [4] P. Areias, D. Dias-da Costa, J.M. Sargado, T. Rabczuk, Element-wise algorithm for modeling ductile fracture with the Rousselier yield function, Comput. Mech. 52 (2013) 1429–1443. [5] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, CalcoloXXI (IV) (1984) 337–344. [6] Y. Bai, T. Wierzbicki., A new model of metal plasticity and fracture with pressure and Lode dependence, Int. J. Plast. 24 (2008) 1071–1096. [7] K.-J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Clifls, New Jersey, 1996. [8] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons, Chichester, Wessex, U.K, 2000. [9] P. Betsch, E. Stein, Numerical implementation of multiplicative elastoplasticity into assumed strain elements with application to shells at large strains, Comput. Method Appl. Mech. 179 (1999) 215–245. [10] K.B. Broberg, Cracks and Fracture, Academic Press, San Diego, California, 1999. [11] T.S. Cao, P. Montmitonnet, P.O. Bouchard, A detailed description of the Gurson Tvergaard Needleman model within a mixed velocity pressure finite element formulation, Int. J. Numer. Methods Eng. 96 (9) (2013) 561–583. [12] B. Carlson, Joseph Keller, Orthogonalization procedures and the localization of Wannier functions, Phys. Rev. 105 (January) (1957) 102–103. [13] C. Chen, O.L. Mangasarian, Smoothing methods for convex inequalities and linear complementarity problems, Math. Program. 71 (1) (1995) 51–69. [14] C. Chen, O.L. Mangasarian, A class of smoothing functions for nonlinear and mixed complementarity problems, Comput. Optim. Appl. 5 (1996) 97–138. [15] C.C. Chu, A. Needleman, Void nucleation effects in biaxially stretched sheets, J. Mater. Process. Technol. 102 (1980) 249–256. [16] A. Gurson, Continuum theory of ductile rupture by void nucleation and growth: Part i—yield criteria and flow rules for porous ductile media, J. Eng. Mater. Technol. 99 (1977) 2–15. [17] J. Korelc, Multi-language and multi-environment generation of nonlinear finite element codes, Engineering with Computers, , 18 (4) (2002) 312–327. [18] J. Lemaitre, A Course on Damage Mechanics, second edition, Springer, Berlin Heidelberg, Germany, 1996. [19] P.-O. Löwdin, On the nonorthogonality problem, Adv. Quantum Chem. 5 (1950) 185–199. [20] J. Lubliner, Plasticity Theory, Macmillan, New York, 1990. [21] L. Malcher, F.M. Andrade Pires, J.M.A. César de Sá, An assessment of isotropic constitutive models for ductile fracture under high and low stress triaxiality, Int. J. Plast. 30–31 (2012) 81–115. [22] L. Malcher, F.M. Andrade Pires, J.M.A. César de Sá, An extended GTN model for ductile fracture under high and low stress triaxiality, Int. J. Numer. Methods Eng. 54 (2014) 193–228. [23] B. Moran, M. Ortiz, C.F. Shih, Formulation of implicit finite element methods for multiplicative finite deformation plasticity, Int. J. Numer. Methods Eng. 29 (1990) 483–514. [24] K. Nahshon, J.W. Hutchinson, Modification of the Gurson model for shear failure, Eur. J. Mech. A—Solid 27 (2008) 1–17. [25] K. Nahshon, Z. Xue, A modified Gurson model and its application to punch-out experiments, Eng. Fract. Mech. 76 (2009) 997–1009. [26] D.M. Norris Jr., B. Moran, J.K. Scudder, D.F. Quiñones, A computer simulation of the tension test, J. Mech. Phys. Solids 26 (1978) 1–19. [27] R.W. Ogden, Non-linear Elastic Deformations, Dover Publications, Mineola, New York, 1997. [28] F.J.P. Reis, L. Malcher, F.M. Andrade Pires, J.M.A. César de Sá, A comparison of shear mechanics for the shear mechanisms for the prediction of ductile failure under low stress triaxiality, Int. J. Struct. Integr. 1 (4) (2010) 314–331. [29] Wolfram Research Inc. Mathematica, 2007. [30] F. Roters, P. Eisenlohr, L. Hantcherli, D.D. Tjahjanto, T.R. Bieler, D. Raabe, Overview of constitutive laws, kinematics, homogenization and multiscale
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55 methods in crystal plasticity finite-element modeling: theory, experiments, applications, Acta Mater. 58 (2010) 1152–1211. [31] G. Rousselier, Ductile fracture models and their potential in local approach of fracture, Nucl. Eng. Des. 105 (1987) 97–111. [32] J.C. Simo, Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory, Comput. Method Appl. Mech. 99 (1992) 61–112. [33] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Corrected Second Printing edition, Springer, New York, 2000.
55
[34] V. Tvergaard, Influence of voids on shear band instabilities under plane strain conditions, Int. J. Fract. 17 (4) (1981) 389–407. [35] V. Tvergaard, A. Needleman, Analysis of cup-cone fracture in a round tensile bar, Acta Metall. 32 (1984) 157–169. [36] G. Weber, L. Anand, Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic–viscoplastic solids, Comput. Method Appl. Mech. 79 (1990) 173–202. [37] L. Xue, Constitutive modeling of void shearing effect in ductile fracture of porous materials, Eng. Fract. Mech. 75 (2008) 3343–3366.
Contents lists available at ScienceDirect
Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel
Semi-implicit finite strain constitutive integration of porous plasticity models P. Areias a,d,n,1, T. Rabczuk b, J. César de Sá c a
Department of Physics, University of Évora, Colégio Luís António Verney, Rua Romão Ramalho, 59, 7002-554 Évora, Portugal Institute of Structural Mechanics, Bauhaus-University Weimar, Marienstraße 15, 99423 Weimar, Germany c Mechanical Engineering Department, Faculty of Engineering, University of Porto, Rua Dr. Roberto Frias, s/n, 4200-465 Porto, Portugal d ICIST, Instituto Superior Técnico, Lisboa, Portugal b
art ic l e i nf o
a b s t r a c t
Article history: Received 11 December 2014 Received in revised form 16 May 2015 Accepted 19 May 2015 Available online 12 June 2015
Two porous plasticity models, Rousselier and Gurson–Tvergaard–Needleman (GTN), are integrated with a new semi-implicit integration algorithm for finite strain plasticity. It consists of using relative Green– Lagrange during the iteration process and incremental frame updating corresponding to a polar decomposition. Lowdin's method of orthogonalization is adopted to ensure incremental frameinvariance. In addition, a smooth replacement of the complementarity condition is used. Since porous models are known to be difficult to integrate due to the combined effect of void fraction growth, stress and effective plastic strain evolution, we perform a complete assessment of our semi-implicit algorithm. Semi-implicit algorithms take advantage of different evolution rates to enhance the robustness in difficult to converge problems. A detailed description of the constitutive algorithm is performed, with the key components comprehensively exposed. In addition to the fully detailed constitutive algorithms, we use mixed finite strain elements based on Arnold's MINI formulation. This formulation passes the inf– sup test and allows a direct application with porous models. Isoerror maps for two common initial stress states are shown. In addition, we extensively test the two models with established benchmarks. Specifically, the cylindrical tension test as well as the butterfly shear specimen are adopted for validation. A 3D tension test is used to investigate mesh dependence and the effect of a length scale. Results show remarkable robustness. & 2015 Elsevier B.V. All rights reserved.
Keywords: Constitutive integration Porous plasticity Finite strains Semi-implicit Löwdin's method
1. Introduction For common metals at room temperature, void nucleation and particle debonding, void growth and coalescence are the mechanisms for crack growth from the process region [10]. Detailed multiple-scale analysis of these mechanisms require robust polycrystalline plasticity and representation of dislocation dynamics and grain boundary barriers (Roters et al. [30] present a comprehensive review). Therefore, the so-called phenomenological models (e.g. [18]) are considered a cost-effective tool for predicting damage and fracture in metals when overall quantities, like maximum load and dissipated energy, are of interest. The now well-established GTN model is based on the original yield surface by Gurson [16] and was subsequently modified by Tvergaard and Needleman [34,35]. The void growth part of the GTN model is based on the micromechanics of the ductile process. Phenomenological models formulated based on thermodynamical n
Corresponding author. 1 Researcher ID: A-8849-2013 http://www.researcherid.com/rid/A-8849-2013 Web-of-Science search: areias, p*
http://dx.doi.org/10.1016/j.finel.2015.05.005 0168-874X/& 2015 Elsevier B.V. All rights reserved.
principles have also been proposed such as the Lemaitre model [18] and its subsequent modifications. Among this class of constitutive models, the Rousselier model [31,4] is interesting from a computational perspective as it is simpler and requires the specification of fewer parameters than the GTN model, while possessing the ability to closely match the predictions of the latter by determination of its parameters. It adheres to the idea that the occurrence of localization, crack initiation and propagation emerges as a direct consequence of strain softening due to void growth. Hence, unlike micro-mechanics models of porous metal plasticity, the Rousselier model does not directly model the mechanism of coalescence. However, two properties are worth mentioning: 1. The Rousselier model predicts void fraction growth in pure shear. 2. The Rousselier yield surface has an isochoric vertex causing a non-zero plastic volumetric term component. As a more complete alternative to the Rousselier model, it was recently (see [22]) created an extended version of the GTN model, which includes the shear void fraction as an independent history variable. We here focus on variants of the GTN and Rousselier
42
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
models to assess the isoerror maps and robustness with a new semi-implicit integration algorithm. To compare the GTN and the Rousselier models, we use the following procedure 1. The GTN model uses the same coalescence approach as the classical Tvergaard–Needleman technique (cf. [35]). 2. Both void fraction and damage variables are integrated in closed-form using the Lambert W and the ERF functions. Concerning the semi-implicit integration method with porous plasticity models, some considerations are worth mentioning: 1. Due to their complementarity form (see, e.g. [20]), finite strain plasticity problems often exhibit convergence difficulties for large strain values and severe sensitivity to step size [3]. We have two measures to attenuate these difficulties: the use of a smoothed complementarity condition and the removal of iterative rotation matrices from the constitutive laws in finite strains. 2. It is computationally convenient that both hyperelastic and finite-strain elasto-plastic laws are implemented in a unique algorithm and applicable to any discretization scheme. 3. Classical assumed-strain elements typically do not directly provide the deformation gradient (e.g. [9]), necessary for many constitutive formulations. An estimated deformation gradient can be calculated from the polar decomposition if an approximate rotation matrix is available. Since a constitutive frame is adopted, the rotation matrix is obtained from this frame in two distinct configurations. Compatible frames use Löwdin's method [19]. 4. Since the seminal contributions of Weber and Anand [36] and Simo [32], Kirchhoff stress tensors (i.e. τ ¼ J σ ) are frequently employed in the yield functions. This is a computational convenience, as commonly adopted elasto-plastic and hyperelastic models are often quasi-incompressible. This is not the case of porous plasticity (cf. [4]). We therefore avoid that approach. 5. Semi-implicit formulations, where certain quantities are fixed in the flow vector (but not the flow vector itself), cf. [23,8] can lead to substantial savings in constitutive integration. We further extend the semi-implicit algorithm presented in [3] to achieve very large time steps. Consistent linearization of integrated form of objective rates is intricate and computationally expensive, making it a possible candidate for the explicit integration part of the semi-implicit scheme. However, rigid body motions are exactly accounted for.
t
t ¼σ v u¼u
u
on Γ ta
2.1. Integration algorithm: equilibrium and relative Green–Lagrange strain Cauchy equations of equilibrium for a rotated reference configuration are obtained from the corresponding spatial equilibrium (derivations for the latter are shown in Ogden [27]). Using standard notation (cf. [27]) we write the spatial equilibrium equations as ð1Þ
with the Cauchy tensor components σij (i; j ¼ 1; 2; 3). In (1) i is the direction index and j is the facet index. The components of the body force vector are identified as bi. In (1), coordinates xaj are the spatial, or deformed, coordinates of a given point under consideration. In addition, the following natural and essential boundary conditions
Γta is the ð2Þ
on Γ a
u
ð3Þ
where t is the known stress vector on Γ v is the outer normal and u is the known displacement field on Γua. It is assumed that (1) and (2)– (3) are satisfied for a time parameter t a A ½0; T with T being the total time of observation and for a point with position xa A Ωa belonging to the deformed position domain at the time of analysis. Equilibrium configuration corresponds to the domain Ωa and is identified by the subscript a. In tensor notation, Eq. (1) can be presented as t a,
∇ σT þ b ¼ 0
ð4Þ
with ∇ ¼ ∂=∂xa being the spatial gradient operator. After multiplication _ integration in the deformed configuration Ωa by the velocity field u, and application of integration by parts component-wise, we obtain the _ int is the internal and W _ ext is the external following power form (W power): Z Z Z σ : L dΩa ¼ b u_ dΩa þ t u_ dΓ a ð5Þ Ωa Ωa Γ ta |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl {zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl } _ _ W int W ext _ where L is the velocity gradient: L ¼ ∂u=∂x a . L can be decomposed into a symmetric part (D, called strain rate) and a skew-symmetric part (W , called vorticity), L ¼ D þW . We also introduce the Jacobian J ¼ det F and F is the deformation gradient. The use of equivalence _ int and external W _ ext power for a continuum has between internal W been used to obtain objective rates of the Cauchy stress σ . In the corotational case, we can rotate the Kirchhoff stress and the strain rate to obtain Z _ int ¼ RT0a τR0a : RT0a DR0a dΩ0 ð6Þ W Ω0
From each element nodal arrangement, we use three directions, where e 1a , e 2a and e 3a to define R0a as R0a ¼ ½e 1a ; e 2a ; e 3a
ð7Þ
where each column is a unit vector and e ia e ja ¼ δij . Section 2.3 describes the calculation of R0a . We then use the following convention for writing tensors: the bold symbol indicates tensor components. Therefore, the frame where the tensor components are defined is identified as a superscript. For example, the secondorder tensor S ab with the equilibrium configuration Ωa and reference configuration Ωb has the following components in frame c:S cab . If a global frame 0 is adopted, component transformation between frame 0 and frame c follows: S cab ¼ Rc0 S 0ab RTc0
2. Constitutive integration in finite strains
∂σ ij þ bi ¼ 0 ∂xaj
hold on each part of the boundary Γ a ¼ Γ a [ Γ a where natural boundary and Γua is the essential boundary:
ð8Þ
Note that this approach, although related, does not correspond to the rotated Algorithm of Areias and Belytschko (cf. [2]). Generalizing (8), change of basis2 from d to c results in the following transformation: S cab ¼ Rcd S dab RTcd
ð9Þ
In addition, change in reference configuration reads S 0ac
¼ Rbc S 0ab RTbc
ð10Þ
Using (10) and transformation (9) for coinciding reference configuration and frame, S cac ¼ S bab
2
This corresponds to a change of frame, denoted by a superscript.
ð11Þ
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
From the relation between the time derivative right Cauchy– Green tensor C and the strain rate D, we obtain C_ ¼ 2F T DF
ð12Þ
Simo [33, p. 281, Eq. 8.1.13] derived the following one-step scheme for integrating the strain rate in frame 0, producing the following “strain” from integration of D: h i e0ab;β ffi 12 F Tbβ F Taβ F aβ I F bβ ð13Þ where β is an arbitrary configuration. Using the previous power conjugacy (6), the internal power can be written in the initial configuration Ω0 as Z _ int ¼ S 0a0 : D0a0 dΩ0 W Ω0
where τ and Using a hypoelastic relation, stress updating is performed as follows: S 0a0
¼ RT0a 0a R0a
D0a0
0 ¼ RT0a da R0a .
S 0a0 ¼ S 0b0 þ ΔS a ðRβ0 e0ab;β RTβ0 Þ
S bab ¼ S bbb þ ΔS a ðRβ0 e0ab;β RTβ0 Þ
where Sba0 ¼ Voigt½S ba0 . We note that, omitting indices a, b, and 0, V S ðFÞ can be written as 2
F 211 6 2 6 F 12 6 6 2 6F V S ðFÞ ¼ 6 13 6 F 11 F 12 6 6 4 F 11 F 13 F 12 F 13
F 221
F 231
2F 21 F 11
2F 31 F 11
2F 31 F 21
F 222 F 223
F 232 F 233
2F 22 F 12
2F 32 F 12
2F 32 F 22
2F 23 F 13
2F 33 F 13
2F 33 F 23
F 21 F 22 F 21 F 23
F 31 F 32 F 31 F 33
F 21 F 12 þF 11 F 22 F 21 F 13 þF 11 F 23
F 31 F 12 þ F 11 F 32 F 31 F 13 þ F 11 F 33
F 22 F 23
F 32 F 33
F 22 F 13 þF 12 F 23
F 32 F 13 þ F 12 F 33
S bab ¼ S bbb þ ΔS a ðebab Þ
ð14Þ
with Simo integration scheme reducing to h i ebab ffi 12 Rb0 F Tab F ab I RTb0
ð15Þ
Ωb in (14) as
1 b 1 1 S ¼ S bbb þ ΔS a ðebab Þ J b0 ab J b0 J b0 |fflffl{zfflffl}S ⋆b |fflffl{zfflffl}σ b
ð16Þ
b
Local constitutive quantities, as well as quadratic strain incre Local frames a and b, and therefore the rotation tensor are explicitly integrated.
Given F bab , ebab (both in frame b), R0b and R0a (both in frame 0) and Tb and Ta Recovered from storage F bb0 , σ 0b , ebb0 (σ b is stored in frame 0 for purpose of representation) Represent Cauchy stress in frame b σ b ¼ V RT σ 0 b
Relevant Jacobian determinants
ð17Þ
then we obtain the Cauchy stress in configuration as 1 a 1 J S ¼ S b ¼ b0 S ⋆b J a0 aa J a0 ab J a0 ab
Ωa and frame a ð18Þ
If ΔS a ðebab Þ is homogeneous of degree one, then it follows that b ~b S ⋆b ab ¼ σ b þ ΔS a e ab
0b
b
J b0 ¼ det F bb0 J ab ¼ det F bab Rab ¼ RT0a R0b
Relative rotation and total deformation gradient update
1 b e J b0 ab
ð19Þ
For hyperelastic models, we require the total strain, which is the strain from configuration Ωa with respect to configuration Ω0. Using frame b, we have the following update formula for the Green–Lagrange strain: b b eba0 ¼ ebb0 þ F bT b0 eab F b0
ð20Þ S ⋆b ab
We can observe that is work-conjugate to eba0 in configuration Ωb. Global Cauchy stress can be obtained as σ 0b ¼ R0b σ bb RT0b . In frame a, eaa0 reads eaa0 ¼ RTba eba0 Rba
Thermal effect in the relative strain eb○ ¼ eb þ αðT a T b ÞI ab ab b○ eba0 ¼ ebb0 þ V E F bT b0 eab
Total strain update
b○ Corrected relative strain for UL e~ ab ¼ J1 eb○ b0 ab b ~ b○ and Qab along with (UL) Determine S⋆b ab ¼ σ b þ ΔS a e ab
S a ∂ΔS a 1 ∂Q ab ab , ∂T a , J ~ b○ and ∂Q sensitivity quantities Cab ¼ J1 ∂Δ~ b○ ∂T a b0 ∂e ab b0 ∂e ab (TL) Determine Sba0 eba0 and Qab, along with sensitivity ∂Sb
quantities Ca0 ¼ ∂ea0 b ,
∂Sba0 ∂Q ab ∂T a , ∂eb a0
and
∂Q ab ∂T a
b b 1 S⋆b ab ¼ J b0 V S F b0 Sa0 (TL) Determine relative sensitivities C ¼ 1 V F b C V ðF bT Þ a0 E ab b0 b0 J b0 S b ∂S⋆b ∂Sa0 b 1 ab ∂T a ¼ J b0 V S F b0 ∂T a ∂Q ab ab ¼ ∂Q V E F bT b0 ∂eb ∂eb a0
(TL)Determine relative stresses
where b e~ ab ¼
S
F ba0 ¼ F bab F bb0
1 b S J b0 ab
σ aa ¼ S ⋆a aa ¼
7 7 7 7 7 7 F 31 F 22 þ F 21 F 32 7 7 7 F 31 F 23 þ F 21 F 33 5 F 32 F 23 þ F 22 F 33
and, for the strain quantities, V E ðFÞ ¼ V TS ðFÞ. The full algorithm is given in Algorithm 1. In summary:
Defining S ⋆b ab ¼
3
Algorithm 1. General relative Lagrangian formulation (Voigt notation adopted).
In the case β ¼ b, the result is
ab
b For symmetric tensors, such as S ⋆b ab and eab , it is preferable to use the Voigt notation. Upright bold symbols denote Voigt form of symmetric tensors: 1 b b ⋆b S⋆b V S F b0 Sa0 ð22Þ ab ¼ Voigt½S ab ¼ J b0
ments are implicitly integrated.
where ΔS a is a function of the strain e0ab;β . This term ΔS a is called the constitutive part of the stress. However, since S 0a0 ¼ S bab and S 0b0 ¼ S bbb , we have
We now insert the Jacobian for configuration
43
ð21Þ
a0
ab
α
∂Sab ∂Sab ∂T a ∂T a þ C ab I eaa0 ¼ V E ðRab Þeba0 F aa0 ¼ Rab F ba0 RTab
Update temperature sensitivity Determine strain in frame a
Determine deformation gradient in frame a Determine Cauchy stress in frame 0 σ 0 ¼ 1 V S ðR0a ÞS⋆b ab a J ab
Store F aa0 , σ 0a , eaa0 Return Qab, S⋆b ab , C ab ,
∂S⋆b ab ∂Q ab ∂T a , ∂eb ab
and
∂Q ab ∂T a
44
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
2.2. Local constitutive system
We determine the relevant constitutive quantities by solving a general nonlinear constitutive system. Omitting the frame superscript, it reduces to the root finding for the following nonlinear system: ! zfflfflfflfflffl}|fflfflfflfflffl{Governing φ Sab ; Q ab ; eKab ; eab ; Ta ; χa ¼0 |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}Unknown |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl}Known |{z}Internalvar: ð23Þ where χ a is the internal variable vector. For the present application, it is defined as ( )
χa ¼
εp fa
ð24Þ
where εp is the effective plastic strain and fa is the void fraction variable. Newton iteration for (23) provides, in general, the solution for the constitutive unknowns Sab , Qab, and χ a . For the present application, system (23) has Q ab ¼ 0, T a ¼ 0, eKab ¼ ∅. A smoothed version of (23) is introduced by using the Chen– Mangasarian replacement functions, which depend on a parameter Error such that (cf. [3]) lim φError ¼ φ
Error-0
ð25Þ
Newton iteration on (23), using φError as a replacement provides the following scheme: ( ) h ∂φ i ΔSab ∂φError Error ð26Þ ¼ φError ∂Sab ∂χ a |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}J Δχ a After achieving a solution using the Newton scheme (26), we calculate the sensitivities from the following equation: ( ) dSab ∂φError deab ð27Þ ¼ J 1 dχ a ∂eab Specifically (with further details being provided in [3]), the present elasto-plastic system is therefore given as 8 9 1 ΔS ab nΔγ > > < eab Clinear =
ð28Þ φError ðS ab ; eab ; χ a Þ ¼ μ⋆ Δγ μ⋆ Δγ þ ϕ Error > > : ; Δχ a Δγωðχ a Þ where Δγ is the plastic multiplier increment, μ⋆ is a dimensional parameter described in [3] and the smooth ramp function of Chen and Mangasarian [14] is used for the third equation, which depends on Error. Function ϕ in (28) is the yield function. In addition, the flow vector n is also required for the flow rule, and it is defined as n¼
∂ϕ ∂Sab
ð29Þ
Fig. 1 shows the effect of Error in the satisfaction of the complementarity condition (see also [3]). In (28), function ωðχ a Þ is the internal variable evolution function, such that _ ðSab ; χ a Þ χ_ a ¼ γω
ð30Þ
The motivation for a unique treatment for all yield functions is that any particular yield function is inserted by means of ϕ, n and derivatives of n. This is accomplished with Mathematica and the AceGen add-on [29,17]. Functions hi and hiError are respectively given by hi ¼
þ j j 2
ð31Þ
Fig. 1. Replacement of μ⋆ Δγ μ⋆ Δγ þ ϕ by μ⋆ Δγ 〈μ⋆ Δγ þ ϕ〉Error as a function of a Error parameter (μ⋆ ¼ 1 is depicted).
hiError ¼ þ
Error lnð2Þ ln 1 þ exp lnð2Þ Error
ð32Þ
where hiError is an approximation of hi. 2.3. Löwdin frame With the objective of uniquely defining the element frame compatible with the concept of polar decomposition, use is made of the symmetric orthogonalization formula by Per-Olov Löwdin (cf. [19]). With that purpose, we resort to the singular value decomposition (SVD). Given a real 3 3 matrix J, there exists a decomposition with the form J ¼ U ΣV T
ð33Þ
where U is a 3 3 unit matrix, U U ¼ I, Σ is a 3 3 diagonal matrix and V is a 3 3 unit matrix with V T V ¼ I. Columns of U are normalized eigenvectors of JJ T and columns of V are normalized eigenvectors of J T J. The diagonal matrix Σ contains the singular values of J in its principal diagonal. Singular values are nonnegative reals. Columns of U are the left singular vectors and columns of V are the right singular vectors. The non-zero singular values of J are the square roots of the non-zero eigenvalues of J T J. In the present case columns of J form a basis. Löwdin's formula provides the following orthogonal frame: T
R ¼ UV T
ð34Þ
It was proved by Carlson and Keller [12] that R obtained by this formula (34) is the proper orthogonal matrix whose columns have the minimum Euclidean distance to the corresponding columns of J. A more concise formula for R is R ¼ JV Σ
1
VT
ð35Þ
and, since Σ is a diagonal matrix, we have the following index notation: 1 Rij ¼ pffiffiffiffiJ ik V kl V jl
λl
ð36Þ
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
45
where λl, l ¼ 1; …; n, are the eigenvalues of J T J and columns of V are the corresponding eigenvectors. In our case, and for a tetrahedron, we have the following specialization: 1
R0a ¼ J a V a Σa V Ta
ð37Þ
where J a is the element Jacobian for configuration Ωa, V a are the eigenvectors of J Ta J a . For the configuration Ωb, we have R0b ¼ J b V b σ b 1 V Tb
ð38Þ
The following notation is used for the Jacobian matrices: Ja ¼
∂xa ∂ξ
ð39Þ
Jb ¼
∂xb ∂ξ
ð40Þ
where ξ is the vector of parent-domain coordinates. The motivation for the use of R0a and R0b to obtain a rotation can be summarized as follows: In problems where local frames are required (anisotropic laws, fracture, shells), rotation can be obtained by multiplication Rab ¼ RT0a R0b . For uniformity, we use the Löwdin local frames in all problems. 2.4. The extended GTN model and its integration The version of the Gurson–Tvergaard–Needleman model introduced by Reis et al. [28] is rewritten here in non-dimensional form: " !# 3ST I6 Tdev S q2 ST I3 2 ϕGTN y; S; f ⋆ ¼ 1 þ q3 f ⋆ 2q1 f ⋆ cosh 2y 2y2 ð41Þ The effective void fraction f ⋆ is a function of the void fraction f . We use the nucleation part of the void fraction following the work of Chu and Needleman [15] using the normal distribution: " # f 1 εp εN 2 f_ N ¼ pNffiffiffiffiffiffi exp ε_ p ð42Þ 2 SN SN 2π where fN, SN and εN are constitutive properties for the nucleation (N) part of the void fraction law. The classical void growth term of Gurson (cf. [16]) is adopted: f G ¼ 1 expð εv Þ
ð43Þ
where εv is the trace of deformation. In addition, the shear part of Xue and Nahshon (cf. [37,24,25]) are combined to obtain the following term: 0 1 2 3q6 6 1=3 B det ½S dev C 1=3 2 fS ¼ ð44Þ @1 54 3 Af εpa 2 π ST I6 Tdev S Including the initial void fraction, f0, then the total f follows the following closed-form laws from the above classical relations: f εN ε εp pffiffiffi ERF N pffiffiffi f ¼ f 0 þ N ERF ð45Þ þ 1 expð εv Þ 2 SN 2 SN 2 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} f 0 1 2 3q6 6 1=3 B det ½S dev C 1=3 2 þ @1 54 3 Af εpa 2 π ST I6 Tdev S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}f S
ð46Þ
The additional property q6 controls the void growth in shear. A depiction of (41) with is shown in Fig. 2. In contrast with the Rousselier model the GTN model is symmetric with respect to the
Fig. 2. Extended GTN yield function (y ¼ q1 ¼ q2 ¼ q3 ¼ 1) with p ¼ ST I3 =3 and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 T 2S I6 Tdev S.
q¼
trace of S. The effect of increasing f can be observed as a uniform reduction in radius of the yield surface. The error function ERFðÞ results from closed-form integration of the Gaussian curve and is available in several open-source libraries and is defined as Z 2 2 ERFðÞ ¼ pffiffiffiffi e t dt ð47Þ
π
0
Two distinct types of effective plastic strain are adopted,3 and εp. The latter being given as
ε_ p ¼
γ_ ST n ð1 f Þy ab
εv
ð48Þ
Integration of constitutive quantities follow the unconditionally stable backward-Euler scheme:
εpa ¼ εpb þ
Δγ
ST n ð1 f Þy ab
ð49Þ
In addition, we have the trace of the flow vector, which is used to obtain εv: ε_ v ¼ γ_ nT I3 ð50Þ It is straightforward to obtain the critical ff, corresponding to collapse, which is given by: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q1 þ q21 q3 ff ¼ ð51Þ q3 The effective void fraction f ⋆ is calculated according to the additional property fc, using the Tvergaard–Needleman modification (cf. [35]): 8 f; f of c > > > > < f f f c f f c ; f c rf rf m ð52Þ f⋆ ¼ fcþ f m f c > > > > :ff; f 4f m The effective void fraction (52) is used in the yield function, and is not updated as is the case of f. 3 No power-consistent closed-form solution exists for the effective plastic strain in the classical Gurson model.
46
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
Fig. 3. Rousselier yield function (y ¼ σ 1 ¼ D ¼ 1) with p ¼ ST I3 =3 and q ¼
2.5. The Rousselier model and its integration
Table 1 Constitutive properties for the two porous models (Rousselier and extended GTN).
The Rousselier [31] model for porous metals is characterized by the following yield function in non-dimensional form: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " # T 3 T σ1 f ⋆ ST I3 2S I6 Tdev S þ D exp ϕR y; S; f ⋆ ¼ 1 ð53Þ y 3 1 f ⋆ σ 1 1f ⋆ y in which f ⋆ is the effective void fraction in configuration Ωa and y is the yield stress. In addition to the hardening characteristics reflected in y, the only additional properties are σ1, D and the initial porosity f0. The property σ1 can be estimated from the original yield stress, y0 ¼ yð0Þ and the ultimate tensile strength, ymax ¼ maxεp y as σ estimate ¼ 13 y0 þ ymax ð54Þ 1 Note that ymax must be provided by the actual hardening curve in the porous case. The effective void fraction f ⋆ depends on the corresponding constitutive void fraction f. A depiction of (53) is shown in Fig. 3. Note that the Rousselier model is asymmetric with respect to the trace of S. The effect of increasing f can be observed. Void fraction follows from volume change and is calculated exactly as in the previous section. Fracture value ff coincides with the zero stress condition fulfilling the condition ϕ ¼ 0: ff ¼
y Dσ 1
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 T 2S I6 Tdev S.
ð55Þ
Property
Rousselier
Extended GTN
Elastic properties Hardening law Initial void fraction Void fraction properties
E, ν y yðεp Þ f0 σ1 D
Nucleation
fN SN εN q6 fc fm
E, ν y yðεp Þ f0 q1 q2 q3 fN SN εN q6 fc fm
Shear Coalescence
where Wð) is the Lambert W function. The effective plastic strain rate is given by this power equivalence (see also [4]):
γ_ ε_ p ¼ STab n
ð58Þ
y
and
ε_ v ¼
γ_
ð1 f Þy
STab I3
nT I3
ð59Þ
Integration of constitutive quantities follow the unconditionally stable backward-Euler scheme:
Δγ
It can be observed that D indirectly controls the void fraction rate. We can force the growth near the critical value of void fraction to model coalescence as in the GTN model, but with a slight difference (see [35]): 8 f; f of c > > > > < f f f c f f c ; f c r f rf m ð56Þ f⋆ ¼ fcþ f m f c > > > > :ff; f 4f m
Time stepping is adapted so that Δεp ¼ εpn þ 1 εpn is kept below 5%. This measure is necessary for accuracy reasons in problems with localization.
with fm, fc and ff being additional properties. An mentioned, initial void fraction f0 must also be provided. The equivalent plastic strain increment is computed based on the power equivalence relation: 1 y ε_ p ð57Þ STab ε_ p ¼ ð1 f Þ y 3σ 1 W Df ⋆ exp 3 3σ 1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}y
For Error ¼ 1 10 3 we test the two porous plasticity models, along with the classical von-Mises yield criterion. Properties are shown in Table 1 for the two presented models. Von-Mises is a particularization. Isoerror maps are shown in Fig. 4 for the three yield functions. We use the properties shown in Table 2 for the runs of the isoerror maps. Plane strain condition is assumed (4 direct stress components
εpa ¼ εpb þ
y
STab n
ð60Þ
2.6. Summary of properties for porosity models and isoerror maps
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
47
Fig. 4. Isoerror maps for von-Mises, Rousselier and GTN yield functions (a single step is used).
and 1 shear stress component). For each strain increment a single step is used. The following two conclusions are drawn:
yield function is not smooth for all stress states. It is also acknowledgedly difficult to integrate (cf. [4]).
Very small integration errors are observed for the von-Mises
yield function. This is well known since in a simpler case an exact solution exists. The Rousselier yield function exhibits a larger error than the GTN yield function. This is now relatively well known since this
2.7. Fitting of Rousselier properties to the GTN model To estimate the properties for the Rousselier model from known properties of the GTN model, we can obtain D, σ1 from q1, q2 and q3
48
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
Table 2 Isoerror maps: relevant properties for two prototype materials (consistent units). Property
Rousselier
Property
Extended GTN
E ν y f0
72 400 0.33 352
E ν y f0
72 400 0.33 352
2
2 10 500 106
σ1 D
2 10 1 0.042 0.1 0.2 1 0.06 0.10
fN SN εN q6 fc fm
q1 q2
2 10 1.5 1.0
q3 fN SN εN q6 fc fm
2.25 0.042 0.1 0.2 1 0.06 0.10
[5], see also Bathe [7] and Cao [11], is based on a two-field formulation where:
Pressure is linearly interpolated using the corner nodes. An internal shape function, named a bubble, enriches the velocity or displacement fields.
2
Cauchy stress is calculated from the constitutive stress S⋆b ab as
σ aa ¼
1 ⋆b S J ab ab
ð66Þ
, Cauchy pressure is obtained as T S⋆b I3 ab pa ¼ 3J ab
ð67Þ
We can therefore write S⋆b ab in Voigt form as a sum of deviatoric and pressure terms:
Initial stress, Case I 352 σxx σyy 352 Initial stress, Case II σxx 352 σyy 0
⋆b S⋆b ab ¼ J ab pa I3 þ Tdev Sab
Table 3 Properties for 2024-T351 aluminum (consistent units are adopted). Rousselier
Extended GTN
E ν y
72 400 0.33 352 þ 668ε0:63302 p
E ν y
72 400 0.33 352 þ 668ε0:63302 p
f0 σ1 D fN SN εN q6 fc fm
0.0 174.51 1.07291 0.042 0.1 0.2 1 0.06 0.99
f0 q1 q2 q3 fN SN εN q6 fc fm
0.0 1.5 1.0 2.25 0.042 0.1 0.2 1 0.06 0.99
where Tdev is the 22 13 3 6 1 2 6 3 3 6 6 1 13 6 Tdev ¼ 6 3 60 0 6 6 0 40 0
0
following sparse matrix: 3 13 0 0 0 7 13 0 0 0 7 7 7 2 0 0 07 3 7 0 1 0 07 7 7 0 0 1 05 0 0 0 1
ð68Þ
ð69Þ
In terms of power balance, we use the following relation, where ⋆b ~ It corresponds to a S~ ab depends on the independent pressure p. classical two-field variational principle: 0 T 1 Z Z S⋆b I3 C ab ⋆b T b B _ ext ð70Þ ¼W S~ ab e_ ab dΩb þ @J ab p~ þ Ap~_ dΩb 3 Ωb Ωb |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}W_ int where the relative Jacobian Jab is used to ensure correct volume calculation. We note that the product J ab p~ cannot be used as an unknown field. In (70), we have the following quantities:
at a given void-fraction level fc. Using least-squares, we solve the following least-squares minimization problem Cf. Table 3: 2 2 i 1h ð61Þ qGTN qRousselier þ pGTN pRousselier min σ 1 ;D 2
⋆b ~ 3 þTdev S⋆b S~ ab ¼ J ab pI ab
ð71Þ
Discretization follows the standard MINI formulation: uðξÞ ¼
5 X
N K ξ uK
ð72Þ
K¼1
where
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qGTN ¼ y 1 þf c q3 2q1
ð62Þ
qRousselier ¼ ð1 f c Þ y σ 1 Df c 2
2y 1 þ f c q3 pGTN ¼ arccosh q2 2f c q1 pRousselier ¼ 3ð1 f c Þσ 1 log
ð63Þ
! ð64Þ
y
Df c σ 1
ð65Þ
Independent pressure p~ is interpolated using the corner nodes: ~ ξÞ ¼ pð
4 X
NK ξ p~ K
ð73Þ
K¼1
N1 ðξÞ ¼ 1 ξ1 ξ2 ξ3
ð74aÞ
N 2 ðξ Þ ¼ ξ 2
ð74bÞ
N 3 ðξ Þ ¼ ξ 3
ð74cÞ
N 4 ðξ Þ ¼ ξ 1
ð74dÞ
3. Element technology for quasi-incompressible problems 3.1. Tetrahedron Quasi-incompressible constitutive laws severely reduce the performance of displacement-based elements. Usually, mixed elements are a solution for avoiding locking in quasiincompressible problems. The low-order MINI element by Arnold
The bubble function, for tetrahedra, is given by N5 ðξÞ ¼ ξ1 ξ2 ξ3 1 ξ1 ξ2 ξ3
ð5Þ
For the calculation of the stiffness matrix, the variation of (70) is required. Not all quantities are determined by hand-derivation, and we use Mathematica [29] with the AceGen (cf. [17]) add-on to
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
49
Fig. 5. Axisymmetric representation of tension specimen and relevant contour plots for the GTN and Rousselier models.
calculate some derivatives. Using (70), we obtain Z ⋆b T ⋆b T b _ int ¼ dS~ ab e_ ab þ S~ ab debab dΩb dW Ωb
T 3 dS⋆b I3 7 ab 6 þ 4dJ ab p~ þ J ab dp~ þ 5p~_ dΩb 3 Ωb Z
⋆b
~ 3 ~ 3 J ab dpI dS~ ab ¼ Tdev Cab debab dJ ab pI ð75Þ
2
Specifically, the terms and dJ ab are determined by AceGen. In terms of rotation tensor, we use the following process to obtain e 1 , e 2 and e 3 .
ð76Þ
3.2. Dimensional reduction
ð77Þ
Dimensional reduction for plane-strain requires careful consideration of out-of-plane stress, which is in general non-null. We therefore explicitly enforce the plane-strain condition by setting 13 and 23 components of both stress and strain as zero, but
where the following notation was used: b dS⋆b ab ¼ C ab deab
ð78Þ
debab
50
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
Fig. 6. Void nucleation and growth in the GTN and Rousselier models (10 000 triangles).
h ⋆b i including the out-of-plane stress S~ ab
33
as active. Therefore, for
plane strain, we only omit out-of-plane shear components of stress and strain. In the axisymmetric case, we use reduced integration, which is sufficient to remove locking in the nearincompressible case. 4. Numerical tests All examples are run in the SIMPLAS framework, cf. [1] with the present formulations. MINI elements were implemented to comply with the Löwdin frames and the relative Green–Lagrange strains. Internal degrees-of-freedom of mixed elements are condensed out. 4.1. Assessment of mesh size sensitivity using classical tension specimens An axisymmetric idealization of a simple tension test with smooth geometry (cf. [6] for a complete testing program) is made, with the relevant dimensions shown in Fig. 5. The corresponding properties are shown in Table 1. Also shown in Fig. 5 are the contour plots for the effective plastic strain and void fraction in both models. Rousselier yield function produces earlier failure due to a steeper void fraction growth in the early nucleation stages. We can observe that no hourglass or mesh instabilities are visible and contour plots are relatively smooth. We compare the void fraction evolution for the two models, separating the nucleation and void growth terms, in Fig. 6. It can
Fig. 7. Aluminum specimen, GTN model, mesh size sensitivity: effective plastic strain, void fraction and reaction–displacement results. Step size Δv ¼ 0:09 mm.
be observed that the growth term in Rousselier model4 is much quicker to develop and tends to localize strain in a sharper form (see also Fig. 5). Evolution of void fraction and effective plastic strain for the GTN yield function is shown in Fig. 7. Reaction–displacement results are also shown in the same figure. We obtain two conclusions:
Effective plastic strain and void fraction evolution are close (however not perfectly close in the case of the effective plastic strain) to the ones reported by Malcher et al. [21]. Reaction– displacement results are very similar to the reported experimental curve (cf. [21])
4 Note that properties have been fitted so that intersections in p- and q-axes for both models coincide at f ¼ f c .
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
51
Fig. 9. Butterfly specimen: relevant dimensions, coarse mesh (4483 nodes and 16 473 MINI tetrahedra) and fine mesh (6351 nodes and 24 132 MINI tetrahedra). One symmetry plane (zz) is taken into account.
Fig. 8. Aluminum specimen, Rousselier model, step size sensitivity: effective plastic strain, void fraction and reaction–displacement results. Mesh with 10 000 triangles and 5151 nodes.
Some mesh sensitivity is observed, however it can be considered mild and acceptable. Concerning the Rousselier model, evolution of void fraction and effective plastic strain is shown in Fig. 8, with mild mesh sensitivity. In terms of reaction–displacement results, results are very robust, with near-absence of mesh sensitivity. 4.2. Butterfly test This test, introduced by Bai and Wierzbicki [6] is required to assess the performance in shear of both yield functions. A description of this test is shown in Fig. 9 with the corresponding geometrical data. Both [21,22] perform extensive testing with this test to calibrate the constitutive properties. In this work, and for this test, Aluminum 2024-T351 is used and both mesh size and time step size sensitivities are inspected.
For the coarse mesh (due to size constraints), we show a sequence of relevant contour plots in Fig. 10. Since the MINI mixed element technology is employed, severe deformations are possible without volumetric locking. In addition, we note that a shear band is formed at the specimen weak region. Effect of step size is shown in Fig. 11 for the GTN model and in Fig. 12 for the Rousselier model. Only mild dependence on step size is observed. A complete depiction of the dependence on mesh size is shown in Fig. 13 where again only slight dependence is observed. We note that, if stronger localization is required, then the mesh size has a serious effect in the results and regularization will be required. Mesh size effect on the effective plastic strain and void fraction evolution at the monitored node are shown in Fig. 14. In this case some effect is noted, which is not sufficient to detract from the overall very robust results. 4.3. Mesh dependence and regularization: 3D tension test To inspect the mesh dependence with full localization (until collapse), we use the geometry from Norris tension test [26] (see also [32]). Fig. 15 shows the tension test specimen with effective
52
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
Fig. 10. Butterfly specimen (coarse mesh): sequence of contour plots for the GTN model (effective plastic strain and void fraction).
Fig. 11. GTN model: comparison with experimental results and numerical results by Malcher et al. [22]. Effect of step size.
plastic strain and void fraction contour plots. Two meshes are used: a coarse one with 3175 nodes and a finer one with 5409 nodes (the latter is used in the contour plots). A fixed displacement step of Δw ¼ 0:02 is adopted for this problem. We use the constitutive properties described in Table 4 which produce full localization. Regularization is performed by introducing a lengthscale lc ratio in Eq. (42) with εp -lc εp , and Eq. (43) with εv -lc εv
Fig. 12. Rousselier model: comparison with experimental results. Effect of step size.
for the GTN model and similarly for the Rousselier model. This is sufficient to produce nearly mesh-insensitive results. The lengthscale ratio is defined as follows: lc ¼
min le lf
ð79Þ
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
53
Fig. 13. Effect of mesh-size (von-Mises, GTN and Rousselier).
Fig. 15. 3D tension test: relevant geometrical data and contour plot for the GTN model.
Table 4 Properties for the 3D tension test (consistent units are adopted). Rousselier
Fig. 14. (a) Effective plastic strain and (b) void fraction evolution, comparison with Malcher et al. [21,22] and the model of Bai and Wierzbicki (cf. [6]).
where le is the edge length and lf is the fixed length size, a constitutive property (given in Table 4). Relatively smooth effective plastic strain and void fraction contour plots are obtained. In terms of displacement/reaction results for the GTN problem, we test the unregularized (lc ¼ 1) and regularized cases with the two meshes (see Fig. 16 for the GTN case). The length scale is sufficient to obtain insensitive results. For the Rousselier model, see Fig. 16, the same conclusion is obtained.
Extended GTN
lf E ν y
0.35 72 400 0.33 352 þ 668ε0:63302 p
lf E ν y
0.35 72 400 0.33 352 þ 668ε0:63302 p
f0 σ1 D fN SN εN q6 fc fm
0.002 20.4 9.53 0.08 0.01 0.05 0.5 0.06 0.2
f0 q1 q2 q3 fN SN εN q6 fc fm
0.002 1.5 1.0 2.25 0.08 0.01 0.05 0.5 0.06 0.2
5. Conclusions We presented a general semi-implicit constitutive algorithm and specialized it for the integration of porous plasticity yield
54
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55
shown in the 3D tension test, where both models were tested with two meshes. In conclusion, for the numerical tests, we assessed the mesh-size independence and the choice of yield function in the 2D tension test, butterfly test and 3D tension test. Very robust results are obtained and the algorithm allows severe strains with accompanying mesh distortions in shear without failure of convergence. In particular, for relatively low values of void fraction, results are in general insensitive to mesh and step sizes without regularization. For fully localized problems, a simple length scale is sufficient to remove mesh dependence near the collapse.
References
Fig. 16. Tension test: comparison of two meshes in non-regularized and regularized cases. (a) GTN model (both non-regularized and regularized cases). (b) Rousselier model (regularized cases).
functions. New contributions and distinctive features of the proposed algorithm are:
The combination of relative Green–Lagrange strains with the
exact corotational method. Quadratic strains are used during iteration whereas a local frame is updated after convergence of Newton iteration. The use of Löwdin algorithm to calculate orthogonal frames based on the Jacobians of configurations Ωb and Ωa. This corresponds to a definition of orthogonal frames consistent with the polar decomposition. The use of a local constitutive system which is smooth and nonlinear. This is achieved in two steps by replacing the complementarity condition for a non-smooth equation and subsequently by the Chen–Mangasarian replacement function (cf. [13,14]). Local Newton iteration is performed without concerns with the complementarity condition. The GTN and Rousselier constitutive models are explored in detail, with the determination of intersections with axes of pressure–shear diagram (p–q). We include isoerror maps and the explicit determination of the void fraction in both models. The use of finite strain versions of the MINI element (cf. [5]) in triangles and tetrahedra to avoid locking due to quasiincompressibility in the early stages of void fraction evolution. Specializations and the second variation of the variational principle, required to calculate the stiffness matrix. The use of a length-scale in the stress measures and effective plastic strain present in the void growth laws, which was found to be sufficient to produce mesh-insensitive results. This is
[1] P. Areias. Simplas. 〈http://www.simplas-software.com〉. [2] P. Areias, T. Belytschko, Analysis of finite strain anisotropic elastoplastic fracture in thin plates and shells, J. Aerosp. Eng. 19 (4) (2006) 259–270. [3] P. Areias, D. Dias-da-Costa, E.B. Pires, J. Infante Barbosa, A new semi-implicit formulation for multiple-surface flow rules in multiplicative plasticity, Comput. Mech. 49 (2012) 545–564. [4] P. Areias, D. Dias-da Costa, J.M. Sargado, T. Rabczuk, Element-wise algorithm for modeling ductile fracture with the Rousselier yield function, Comput. Mech. 52 (2013) 1429–1443. [5] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, CalcoloXXI (IV) (1984) 337–344. [6] Y. Bai, T. Wierzbicki., A new model of metal plasticity and fracture with pressure and Lode dependence, Int. J. Plast. 24 (2008) 1071–1096. [7] K.-J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Clifls, New Jersey, 1996. [8] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons, Chichester, Wessex, U.K, 2000. [9] P. Betsch, E. Stein, Numerical implementation of multiplicative elastoplasticity into assumed strain elements with application to shells at large strains, Comput. Method Appl. Mech. 179 (1999) 215–245. [10] K.B. Broberg, Cracks and Fracture, Academic Press, San Diego, California, 1999. [11] T.S. Cao, P. Montmitonnet, P.O. Bouchard, A detailed description of the Gurson Tvergaard Needleman model within a mixed velocity pressure finite element formulation, Int. J. Numer. Methods Eng. 96 (9) (2013) 561–583. [12] B. Carlson, Joseph Keller, Orthogonalization procedures and the localization of Wannier functions, Phys. Rev. 105 (January) (1957) 102–103. [13] C. Chen, O.L. Mangasarian, Smoothing methods for convex inequalities and linear complementarity problems, Math. Program. 71 (1) (1995) 51–69. [14] C. Chen, O.L. Mangasarian, A class of smoothing functions for nonlinear and mixed complementarity problems, Comput. Optim. Appl. 5 (1996) 97–138. [15] C.C. Chu, A. Needleman, Void nucleation effects in biaxially stretched sheets, J. Mater. Process. Technol. 102 (1980) 249–256. [16] A. Gurson, Continuum theory of ductile rupture by void nucleation and growth: Part i—yield criteria and flow rules for porous ductile media, J. Eng. Mater. Technol. 99 (1977) 2–15. [17] J. Korelc, Multi-language and multi-environment generation of nonlinear finite element codes, Engineering with Computers, , 18 (4) (2002) 312–327. [18] J. Lemaitre, A Course on Damage Mechanics, second edition, Springer, Berlin Heidelberg, Germany, 1996. [19] P.-O. Löwdin, On the nonorthogonality problem, Adv. Quantum Chem. 5 (1950) 185–199. [20] J. Lubliner, Plasticity Theory, Macmillan, New York, 1990. [21] L. Malcher, F.M. Andrade Pires, J.M.A. César de Sá, An assessment of isotropic constitutive models for ductile fracture under high and low stress triaxiality, Int. J. Plast. 30–31 (2012) 81–115. [22] L. Malcher, F.M. Andrade Pires, J.M.A. César de Sá, An extended GTN model for ductile fracture under high and low stress triaxiality, Int. J. Numer. Methods Eng. 54 (2014) 193–228. [23] B. Moran, M. Ortiz, C.F. Shih, Formulation of implicit finite element methods for multiplicative finite deformation plasticity, Int. J. Numer. Methods Eng. 29 (1990) 483–514. [24] K. Nahshon, J.W. Hutchinson, Modification of the Gurson model for shear failure, Eur. J. Mech. A—Solid 27 (2008) 1–17. [25] K. Nahshon, Z. Xue, A modified Gurson model and its application to punch-out experiments, Eng. Fract. Mech. 76 (2009) 997–1009. [26] D.M. Norris Jr., B. Moran, J.K. Scudder, D.F. Quiñones, A computer simulation of the tension test, J. Mech. Phys. Solids 26 (1978) 1–19. [27] R.W. Ogden, Non-linear Elastic Deformations, Dover Publications, Mineola, New York, 1997. [28] F.J.P. Reis, L. Malcher, F.M. Andrade Pires, J.M.A. César de Sá, A comparison of shear mechanics for the shear mechanisms for the prediction of ductile failure under low stress triaxiality, Int. J. Struct. Integr. 1 (4) (2010) 314–331. [29] Wolfram Research Inc. Mathematica, 2007. [30] F. Roters, P. Eisenlohr, L. Hantcherli, D.D. Tjahjanto, T.R. Bieler, D. Raabe, Overview of constitutive laws, kinematics, homogenization and multiscale
P. Areias et al. / Finite Elements in Analysis and Design 104 (2015) 41–55 methods in crystal plasticity finite-element modeling: theory, experiments, applications, Acta Mater. 58 (2010) 1152–1211. [31] G. Rousselier, Ductile fracture models and their potential in local approach of fracture, Nucl. Eng. Des. 105 (1987) 97–111. [32] J.C. Simo, Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory, Comput. Method Appl. Mech. 99 (1992) 61–112. [33] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Corrected Second Printing edition, Springer, New York, 2000.
55
[34] V. Tvergaard, Influence of voids on shear band instabilities under plane strain conditions, Int. J. Fract. 17 (4) (1981) 389–407. [35] V. Tvergaard, A. Needleman, Analysis of cup-cone fracture in a round tensile bar, Acta Metall. 32 (1984) 157–169. [36] G. Weber, L. Anand, Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic–viscoplastic solids, Comput. Method Appl. Mech. 79 (1990) 173–202. [37] L. Xue, Constitutive modeling of void shearing effect in ductile fracture of porous materials, Eng. Fract. Mech. 75 (2008) 3343–3366.