- Email: [email protected]
Fracture in sheet metal forming: Effect of ductile damage evolution
Computers and Structures 85 (2007) 205–212 www.elsevier.com/locate/compstruc
Fracture in sheet metal forming: Effect of ductile damage evolution M. Khelifa
a,*
, M. Oudjene b, A. Khennane
c
a Laboratory LPMM (UMR 7554), ENSAM de Metz, 4 Rue Augustin Fresnel, 57078 Metz Cedex 3, France ERMeP, Insitut Supe´rieur d’Inge´nierie de la Conception, 27 rue d’Hellieule, 88100 Saint-Die´-des-Vosges, France Computational Engineering Research Centre, University of Southern Queensland Toowoomba, QLD 4350, Australia b
c
Received 31 October 2005; accepted 8 August 2006 Available online 25 October 2006
Abstract This work deals with the virtual simulation of the sheet metal stamping process. The main objective is to predict when and where the cracks can appear in the workpiece during the forming operation. A local approach based on the strong coupling between anisotropic elastoplasticity with mixed nonlinear work hardening (isotropic and kinematic) and an isotropic ductile damage is proposed. The theoretical and numerical aspects of the constitutive equations are, first, presented. The resolution of the resulting system of equations is carried out via a Vumat user material, using ABAQUS/Explicit finite element code. The results obtained, in the context of Swift’s benchmark deep-drawing test show the efficiency and the potential interest of the proposed damage model. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Plastic anisotropy; Ductile damage; Large plastic strains; Contact and friction; Metal forming; Deep drawing processes; FEM
1. Introduction Sheet metal forming processes are among the oldest and most widely used industrial manufacturing processes. They allows to produce thin walled parts (also called workpiece) of complicated shape (automotive panels, structural parts in various domains). The process consists, in general, in the plastic deformation of an initial flat blank subjected to the action of the tools (punch, die) while constrained on the periphery by a blank-holder. During the last decade, several commercial codes based on incremental or on inverse approaches are developed and updated for deep-drawing simulation of thin sheets. In particular, the incremental codes are based on implicit or explicit algorithms and allow taking into account large elastoplastic strains and contact conditions with friction between the tools and the sheet. The usefulness of incremental codes (implicit or explicit) has been recognized to evaluate deformation path as well as forming defects
*
Corresponding author. Tel.: +33 3 87 31 53 60. E-mail address: [email protected] (M. Khelifa).
0045-7949/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2006.08.053
such as fracture, wrinkling and springback. The numerical simulations of the forming process have been conducted abundantly and used extensively for the analysis and design of industrial parts to avoid the long and expensive experimental tryout procedures and also to help engineers to better understand metal forming processes [1]. Numerous pertinent references on the subject are available in the literature. It is well known that the metallic material is subjected to large irreversible deformation in a sheet forming process such as stamping or deep-drawing, leading to high strain localization zones and then internal or superficial microdefects (cracks due to the ductile damage) [1]. This damage apparition and its growth cause some quality problems such as necking and fracture which lead to stop the process. These cracks are due to the ductile tearing of the sheet in the zones where it has attained its limit of plastic deformation. This problem has been treated using the Forming Limit Curves (FLC) method, which consists in carrying out simple plastic deformation tests and plotting on a graph the limiting values of the principal strains for which the rupture of the test sample is observed [2]. Thus, one obtains a Forming Limit Curve which characterizes the
206
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
limit of deformation of the sheet. The FLC curve is supposed intrinsic for each sheet and used to predict the risk of appearance of defects (cracks) during any forming process. It allows only to obtain, at each time step, the contours of plastic deformation in the sheet at any process simulation stage. The zones in which the values of the limiting deformations reach the experimental ones, are supposed to be the regions of local rupture [3]. It is clear that this method is not of a predictive nature. It can only be used to compare sheets between them in terms of their forming properties. In particular, it does not make it possible to predict when and where the cracks can appear in the workpiece during the forming operation. Some research works have been published on the damage models based either on continuum damage mechanics [1,4,5], or on Gurson’s theory [6–9], among others. In the present work, we develop and evaluate a damage model, based on the strong coupling between anisotropic elastoplasticity with mixed nonlinear work hardening (isotropic and kinematic) and an isotropic ductile damage. The model is formulated using the theory of engineering plasticity coupled with damage within the thermodynamic framework of irreversible processes [10,11]. Tensile tests carried out on aluminium samples are used to identify the parameters of the model. At present, this paper deals with isotropic damage. Anisotropic damage issue is another difficult subject for future research work. In this paper we, first, recall the main aspects of the formulation of the plasticity-damage model. Then we apply the model in the context of the stamping process via a Vumat user material in the framework of ABAQUS/Explicit finite element code. The data are taken from the Swift’s deep-drawing test. Finally, some discussions and comparisons with experiments are drawn. The obtained results show its ability to predict when and where the ductile damage can appear in the workpiece during the deep-drawing operation. 2. Plasticity-damage model
the isotropic hardening and its associated variable, (D, Y) represents the isotropic ductile damage and its associated variable, K is the fourth-order operator of elastic properties. C and Q are respectively the kinematic and isotropic hardening parameters. The mechanical dissipation is described using the theory, which considers only one potential of dissipation F. It controls both the plastic flow f(r, X, R, D) and the evolution of the damage F(r, X, R, D), with one plastic multiplier k_ [13–15]. kr X k R pffiffiffiffiffiffiffiffiffiffiffiffi ry ¼ 0 ð5Þ 1D 3 a 1 b R2 F ðr; X ; R; DÞ ¼ f þ X :Xþ 4 Cð1 DÞ 2 Qð1 DÞ sþ1 1 S Y Y0 þ ð6Þ b S ð1 DÞ s þ 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with kr X k ¼ ðr X Þ : H : ðr X Þ. ry is the yield
f ðr; X ; R; DÞ ¼
strength. a and b are, respectively, the kinematic and isotropic nonlinear hardening parameters. S, Y0, b and s describe the isotropic damage evolution; while H is the fourth-order operator of Hill [10,16,17] as a function of six constants F, G, H, L, M and N [10,18]. Using the normality rule, complementary relations are derived as follows [10,19]: H : ðr X Þ k_ ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ~n e_ p ¼ k_ pffiffiffiffiffiffiffiffiffiffiffiffi 1D 1 Dkr X k p _ a_ ¼ e_ kaa 1 r_ ¼ k_ pffiffiffiffiffiffiffiffiffiffiffiffi br 1D s k_ Y Y0 _D ¼ b S ð1 DÞ
ð7Þ ð8Þ ð9Þ ð10Þ
where ~n is the normal to the loading surface. The plastic multiplier is determined using the consistency condition f ¼ f_ ¼ 0.
2.1. Theoretical aspects
2l f_ ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ~n : e_ hk_ ¼ 0 1D
The plasticity damage model presented here is based on the thermodynamic approach of irreversible processes with internal variables. The state relations are given hereafter as follows [12]:
where h > 0 is the tangent elastoplastic modulous defined such that: pffiffiffiffiffiffiffiffiffiffiffiffi R h ¼ 2l 1 D þ C ~n : ~n þ Q b pffiffiffiffiffiffiffiffiffiffiffiffi 1D s a 1 Y Y 0 hry i pffiffiffiffiffiffiffiffiffiffiffiffi ~n : X þ ð12Þ S 2 1D ð1 DÞbþ1
r ¼ ð1 DÞK : ee 2 X ¼ Cð1 DÞa 3 R ¼ ð1 DÞQr 1 1 1 Y ¼ ee : K : ee þ Ca : a þ Qr2 2 3 2
ð1Þ ð2Þ ð3Þ ð4Þ
where the couple (ee, r) represents the tensors of elastic strain and Cauchy stress, (a, X) represents the kinematic hardening and its associated variable. (r, R) represents
ð11Þ
The formulation of the above model in large deformations is done using the rotating referential. The rotation Q of the deformed configuration is at each time solution of a kinematic equation [20,21] _ T ; Qðt ¼ 0Þ ¼ I W Q ¼ QQ where I is the identity tensor of order 2.
ð13Þ
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
The elastic strains are assumed infinitesimal in order to ensure the summation of the strain rates D ¼ e_ Je þ Dp
ð14Þ
where e_ Je is the Jaumann’s rotational derivative of the infinitesimal elastic strain. 2.2. Numerical aspects The finite element method is used to discrete the global equilibrium equation expressed as follows:
€ þ fF int g fF ext g ¼ 0 ½M U ð15Þ where [M] is the mass matrix. {Fint} and {Fext} are, respectively, the internal and the external force vectors. The above equation is solved explicitly. The numerical scheme consists in obtaining solutions at each time t+Dt as a function of the known quantities at time t
€ n ¼ ½M1 ½ffext g ffint g U ð16aÞ
Dtnþ1 þ Dtn
€n U_ nþ1=2 ¼ U_ n1=2 þ ð16bÞ U 2
ð16cÞ fU nþ1 g ¼ fU n g þ Dtn þ 1 U_ nþ1=2 The explicit scheme does not require iterating for the determination of the tangent matrix. Its stability and accuracy depend strongly on the size of the time step Dt [22]. The discretization of Eqs. (7)–(10) is associated with an asymptotic integration of Eqs. (8) and (9). Then, the incremental equations become: Dk Dpnþ1 ¼ Dpn þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ nnþ1 1 Dnþ1 ð1 eaDk Þ anþ1 ¼ an eaDk þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ nnþ1 a 1 Dnþ1 ð1 ebDk Þ rnþ1 ¼ rn ebDk þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b 1 Dnþ1 s Y nþ1 Y 0 1 Dnþ1 ¼ Dn þ Dk b S ð1 Dnþ1 Þ
ð17Þ ð18Þ ð19Þ ð20Þ
The constitutive equations need to be integrated over each time to obtain the values of the variables rn+1, Xn+1, Rn+1 and Yn+1. The integration of the constitutive equations allows to reduce the number of equations from (15) to (8) [11,16,23]. The system of equations to solve becomes then fnþ1 ðDk; ~ nnþ1 ; Dnþ1 Þ Rnþ1 krnþ1 X nþ1 k ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ry ¼ 0 1 Dnþ1 1 Dnþ1 hnþ1 ðDk; ~ nnþ1 ; Dnþ1 Þ ¼ H : ðrnþ1 X nþ1 Þ krnþ1 X nþ1 k~ nnþ1 ¼ 0 nnþ1 ; Dnþ1 Þ gnþ1 ðDk; ~ ¼ Dnþ1 Dn Dk
s Y nþ1 Y 0 1 ¼0 S ð1 Dnþ1 Þb
ð21aÞ
ð21bÞ
ð21cÞ
207
with krnþ1 X nþ1 k ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðrnþ1 X naþ1 Þ : H : ðrnþ1 X nþ1 Þ
rnþ1 ¼ rn þ ð1 Dnþ1 Þ½2lDe þ ke trðDeÞ1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2G 1 Dnþ1 Dk~nnþ1 ð1 eaDk Þ aDk X nþ1 ¼ Cð1 Dnþ1 Þ an e þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~nnþ1 a 1 Dnþ1 ð1 ebDk Þ bDk þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Rnþ1 ¼ Qð1 Dnþ1 Þ rn e b 1 Dnþ1
ð22Þ
ð23Þ ð24Þ ð25Þ
where G is the shear modulous; ke and l are the Lame constants. The unknown quantities are Dk; ~nnþ1 and Dnþ1 . The solution of the system of equations is obtained with the Newton–Raphson method [11]. 2.3. Implementation within ABAQUS The determination of the quantities rn+1, Xn+1, Rn+1 and Dn+1 at the time tn+1, is obtained using the radial return algorithm proposed by Simo [23]. (i) Elastic prediction: Given a total strain increment De at the time tn, initially the response is assumed to be purely elastic; that is, the stress is obtained as rnþ1 ¼ rn þ ð1 Dn Þ½2lDe þ ke trðDeÞ1
ð26Þ
Substituting the value of rnþ1 in the expression of the criterion results in r Rn ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ry fnþ1 ¼ pnþ1 ð27Þ 1 Dn with qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ r : H : r nþ1 nþ1 nþ1
ð28Þ
X n ¼ Rn ¼ 0
ð29Þ
if fnþ1 < 0, the mechanical response is purely elastic. The stress at the time tn+1 can be written as follows:
rnþ1 ¼ rnþ1
ð30Þ
Dknþ1 ¼ Dkn ; ~nnþ1 ¼ ~nn ; Dnþ1 ¼ Dn
ð31Þ
Table 1 Material parameters of the used aluminium alloy (AL 5754) Elasticity E (MPa) 70000
m 0.3
Plasticity ry (MPa) 114.1 F 0.89
C (GPa) 2.1 G 0.67
a 21.1 H 0.33
L 1.5
Damage S (MPa) 17
s 13
b 9
Y0 (MPa) 0.01
Q (MPa) 554.2 M 1.5
b 3.67 N 1.65
208
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
(ii) Plastic-damage behaviour (radial correction): However, if fnþ1 > 0, the Newton–Raphson scheme is used to obtain the unknowns parameters Dk; ~ n and D of the system (Eq. (21))
8 > < fnþ1 ðDk; ~nnþ1 ; Dnþ1 Þ ¼ 0 hnþ1 ðDk; ~nnþ1 ; Dnþ1 Þ ¼ 0 > : gnþ1 ðDk; ~nnþ1 ; Dnþ1 Þ ¼ 0
Fig. 1. Equivalent plastic strain and damage distributions for the rolling direction (0°).
ð32Þ
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
Successive corrections are used to tend Eq. (21) to be equal to zero, which results in updated values of the unknowns Dk; ~n and D at the time tn+1. When convergence occurs,
209
the solution is stored and the variables rn+1, Xn+1, Rn+1 and Dn+1 are updated.
3. Numerical applications 3500
3.1. Model parameters identification: tension test
3000
To identify the model parameters, a tensile test was carried out on an aluminium sample. The geometrical characteristics of the tensile specimen are 1 mm thick, 12.52 mm width and 180 mm long. The experiment is also simulated using ABAQUS/Explicit, by considering only the useful part. The finite element mesh consists of 2366 eight-nodes hexahedral elements. The part is meshed as follows: 91 elements along the length, 13 elements in the width direction and 2 elements in the thickness one. The simulation is carried out using a displacement rate of 10 mm/s and by considering a density of 2700 kg/m3. The different material properties and parameters are listed in Table 1. Fig. 1 shows the distributions of effective plastic strain and damage for different displacements in the direction of rolling (0°). For this direction, it appears clearly that there are three stages of damage evolutions. Initially, for a displacement less than u = 16.8 mm, the strain distribution is homogeneous as shown in Fig. 1a. Beyond this limit, the deformation becomes heterogeneous. A highly strained
Load [N]
2500
2000
Experiment Numerical
1500
1000
500
0 0
2
4
6
8
10
12
14
16
18
20
22
Displacement [mm]
(a) : Direction 0° 3000
2500
Load [N]
2000
Experiment Numerical
1500
Table 2 Maximum loads and displacements recorded at failure 0°
1000
Experiment Numerical Error (%)
500
45°
90°
u (mm)
F (daN)
u (mm)
F (daN)
u (mm)
F (daN)
14.9 19.8 32.7
290.5 302.3 4.1
18.4 19.6 6.8
281.6 286.6 1.7
19.1 19.8 3.9
279.8 284.2 1.5
0 0
2
4
6
8
10
12
14
16
18
20
22
Displacement [mm]
(b) : Direction 45° 3000
2500
Load [N]
2000
Experiment Numerical
1500
1000
500
0 0
2
4
6
8
10
12
14
16
18
20
22
Displacement [mm]
(c) : Direction 90° Fig. 2. Load–displacement curves.
Fig. 3. Schematic illustration of the swift process tools.
210
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
zone appears in the middle of the specimen at a displacement of 18.6 mm (Fig. 1b). This is followed by the quick formation of a constriction in the form of two shear bands, before a macroscopic crack starts at the intersection of the two bands at a displacement of 19.8 mm. The value of the damage (SDV21) reaches the maximum value D = 1 in the cracked zones (Fig. 1c). It can be seen that the rupture facies is perpendicular to the longitudinal axis of the sample. The same observations are also made in the 45° and 90° directions with respect to the direction of rolling. The damage evolution also follows three stages.
The obtained results in each direction (0°, 45° and 90°) are also reported in the form of load–displacements curves as shown in Fig. 2. Except for a very slight deviation above a displacement of 6 mm, the obtained results agree very well with their experiment counterparts. The failure loads in the 0° and 45° directions are not as well predicted as for the 90° direction. The reason is that the damage considered here is isotropic. Table 2 shows the recorded values for load and displacement at failure for each direction. Again, it can be seen that the failure load and the corresponding displacements are as well predicted for the 45° and 90° directions.
Fig. 4. Equivalent plastic strain and damage distributions for different punch displacements (travels).
Fig. 5. Predicted sheet shapes for different punch travels (0, 6.63, 9.12, 10.64 and 11.3 mm.
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
After identification of the model parameters, the Swift’s deep-drawing test was used to examine the capability of the developed damage model presented above. This test was also reproduced experimentally at the CETIM of Senlis (France). The schematic deep-drawing test is shown in Fig. 3. The blank is a circular sheet of 73 mm of diameter and 1 mm thick. The blank is meshed using eight-node hexahedral elements with selective reduced integration, while the other tools (punch, die, blank-holder) were regarded as rigid bodies. Coulomb’s friction law was assumed and l = 0.17 was used for the friction coefficient between the sheet and the tools. Fig. 4 shows the distributions of effective plastic strain and damage as a function of the punch displacement. At the beginning (punch travel of u = 7.5 mm), the damage was located along a circumference zone which comes in contact with the punch radius. It can be seen that the damage increases in the circular zone, which decreases in diameter, as the punch displacement increases. A strong localization of damage was induced at a displacement of 10.3 mm, which quickly transformed into macroscopic cracks at a displacement of 11.1 mm.
Fig. 5 presents the predicted shapes of the sheet for different punch travels u = 0, 6.63, 9.12, 10.64 and 11.3 mm. Fig. 6a and b shows a comparison between the experimental rupture and that predicted using the developed damage model. In a deep-drawing test, the risk of rupture is evaluated by using the Forming Limit Curves (FLC) defined in the plane of principal strains as shown in Fig. 7. The strain values at different depths of the punch are measured in order to estimate the likelihood of crack nucleation in the high risk zones. The experimental and numerical points are drawn from the framed zone represented in Fig. 6b. The points located above the FLC curve represent the rupture limit of the sheet. Fig. 8 shows comparison between the experimental and numerical results with regard to the punch force-travel.
0.50
Limit curve Experimental points Numerical points
0.45 0.40
Major strain
3.2. Application to the deep-drawing: Swift’s test
211
0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00
-0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
Minor strain Fig. 7. Forming Limit Curves FLC.
30000
Punch load [N]
25000
Experiment Numerical
20000
15000
10000
5000
0 0
2
4
6
8
10
12
Punch displacement [mm]
Fig. 6. Fracture in the workpiece: comparison between simulation and experiment.
Fig. 8. Punch force–displacement curves: comparison between simulation and experiment.
212
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
The maximum computed punch force is equal to 25.93 kN and it was found in good agreement with the experimental value of 25.92 kN. The rupture was produced at the same time (for punch travel of u = 11.3 mm) in both experimental and simulation procedures. 4. Conclusion A predictive damage model, based on the theory of damage-plasticity, is proposed and successfully implemented in ABAQUS/Explicit finite element code. It allows to predict the onset of macroscopic cracks in thin sheet metals during deep-drawing operations. The experimental FLC method is a post processing at the end of simulation stage. It neglects the interaction between plasticity and damage. The combination of plasticity and damage allows essentially to avoid the over-estimation of the formability of sheet materials. The proposed damage model gives good results in comparison with experiments. The present procedure allows to predict the appearance of ductile rupture and its growth during deep-drawing operations. Acknowledgements CETIM of Senlis and FEDER of ‘‘l’Etat de la region Champagne–Ardenne’’ are gratefully acknowledged for their contributions and financial supports, respectively. References [1] Guo YQ, Li YM, Bogard F, Debrey K. An efficient pseudo-inverse approach for damage modelling in the sheet metal forming process. J Mater Process Technol 2004;151:88–97. [2] Sfar H, et Maillard A. Applications de la simulation de l’emboutissage, CETIM, Rap. No1B8422, De´cembre 2000. [3] Sfar H, et Maillard A. Aides a` la re´alisation de l’emboutissage: simulation et mesure automatique des de´formations, CETIM, Recueil de confe´rences de la journe´e d’information du, 17 Novembre 1998. [4] Lee H, Peng K, Wang J. An anisotropic damage criterion for deformation instability and its application to forming limit analysis of metal plates. Eng Frac Mech 1985;21(5):1031–54.
[5] Zhu YY, Cescotto S, Habraken AM. A fully coupled elastoplastic damage modelling and fracture criteria in metal forming processes. J Mater Process Technol 1995;32:197–204. [6] Gelin JC, Oudin J, Ravalard Y. An imposed finite element method for the analysis of damage and ductile fracture in cold metal forming processes. Ann CIRP 1985;34(1):209–13. [7] Onate E, Kleiber M. Plastic and viscoplastic flow of void containing metal applications to axisymmetric sheet forming problems. Int J Num Meth Eng 1988;25:237–51. [8] Brunet M, Sabourin F, Mghib-Touchal S. The prediction of necking and failure in 3D sheet forming analysis using damage variable. J Phys III 1996(6):473–82. [9] Boudeau N, Gelin JC. Necking in sheet metal forming: influence of macroscopic and microscopic properties of materials. Int J Mech Sci 2000;42:2209–32. [10] Khelifa M, Badreddine H, Belamri N, Gahbiche M-A, Saanouni K, Cherouat A, et al. Int J Form Process 2005;8(N°2). [11] Khelifa M. Simulation nume´rique de l’endommagement en formage des structures minces, The`se de Doctorat, Universite´ de Technologie de Troyes, France, 2004. [12] Saanouni K, Forster C, BenHatira F. On the inelastic flow with damage. Int J Damage Mech 1994;3(2):140–69. [13] Chaboche J-L. Le concept de contrainte effective applique´ a` l’e´lasticite´ et a` la viscoplasticite´ en pre´sence d’un endommagement anisotrope, In: Boehler J-P, editor. Colloque EuroMech 115, Grenoble, Editions du CNRS, 1979, p. 737–60. [14] Lemaitre J. A Continuum damage mechanics model for ductile fracture. J Eng Mat Tech 1985;107:83–9. [15] Lemaitre J. A course on damage mechanics. Springer Verlag; 1996. [16] Chaboche J-L, Cailletaud G. Integration methods for complex plastic constitutive equations. Computer methods in applied mechanics and engineering, vol. 133. Elsevier; 1996. p. 125–55. [17] Hill R. A theory of yielding and plastic flow of anisotropic metals. Royal Soc Lon Proc 1948:281. [18] Lemaitre J, Chaboche JL. Me´canique des mate´riaux solides, Edition Dunod, 1988. [19] Saanouni K, Chaboche J-L. Computational damage mechanics. In: Milne I, Ritchie RO, Karihaloo B, editors. Comprehencive structural integrity. de Borst R, Mang HA, editors. Application to metal forming, numerical and computational methods, vol. 3, 2003, ISBN: 0-08-043749-4 (Chapter 7). [20] Sidoroff F. The geometrical concept of intermediate configuration and elastic–plastic finite strain. Arch Mech 1973;25(N°2):299–308. [21] Dogui A. Plasticite´ anisotrope en grandes de´formations, The`se de doctorat d’e´tat, Universite´ Claude Bernard, lyon I, 1989. [22] ABAQUS, Theory manual, Version 6.4, Hibbit, Karson and Sorensen, Inc., 2004. [23] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998.
Fracture in sheet metal forming: Effect of ductile damage evolution M. Khelifa
a,*
, M. Oudjene b, A. Khennane
c
a Laboratory LPMM (UMR 7554), ENSAM de Metz, 4 Rue Augustin Fresnel, 57078 Metz Cedex 3, France ERMeP, Insitut Supe´rieur d’Inge´nierie de la Conception, 27 rue d’Hellieule, 88100 Saint-Die´-des-Vosges, France Computational Engineering Research Centre, University of Southern Queensland Toowoomba, QLD 4350, Australia b
c
Received 31 October 2005; accepted 8 August 2006 Available online 25 October 2006
Abstract This work deals with the virtual simulation of the sheet metal stamping process. The main objective is to predict when and where the cracks can appear in the workpiece during the forming operation. A local approach based on the strong coupling between anisotropic elastoplasticity with mixed nonlinear work hardening (isotropic and kinematic) and an isotropic ductile damage is proposed. The theoretical and numerical aspects of the constitutive equations are, first, presented. The resolution of the resulting system of equations is carried out via a Vumat user material, using ABAQUS/Explicit finite element code. The results obtained, in the context of Swift’s benchmark deep-drawing test show the efficiency and the potential interest of the proposed damage model. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Plastic anisotropy; Ductile damage; Large plastic strains; Contact and friction; Metal forming; Deep drawing processes; FEM
1. Introduction Sheet metal forming processes are among the oldest and most widely used industrial manufacturing processes. They allows to produce thin walled parts (also called workpiece) of complicated shape (automotive panels, structural parts in various domains). The process consists, in general, in the plastic deformation of an initial flat blank subjected to the action of the tools (punch, die) while constrained on the periphery by a blank-holder. During the last decade, several commercial codes based on incremental or on inverse approaches are developed and updated for deep-drawing simulation of thin sheets. In particular, the incremental codes are based on implicit or explicit algorithms and allow taking into account large elastoplastic strains and contact conditions with friction between the tools and the sheet. The usefulness of incremental codes (implicit or explicit) has been recognized to evaluate deformation path as well as forming defects
*
Corresponding author. Tel.: +33 3 87 31 53 60. E-mail address: [email protected] (M. Khelifa).
0045-7949/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2006.08.053
such as fracture, wrinkling and springback. The numerical simulations of the forming process have been conducted abundantly and used extensively for the analysis and design of industrial parts to avoid the long and expensive experimental tryout procedures and also to help engineers to better understand metal forming processes [1]. Numerous pertinent references on the subject are available in the literature. It is well known that the metallic material is subjected to large irreversible deformation in a sheet forming process such as stamping or deep-drawing, leading to high strain localization zones and then internal or superficial microdefects (cracks due to the ductile damage) [1]. This damage apparition and its growth cause some quality problems such as necking and fracture which lead to stop the process. These cracks are due to the ductile tearing of the sheet in the zones where it has attained its limit of plastic deformation. This problem has been treated using the Forming Limit Curves (FLC) method, which consists in carrying out simple plastic deformation tests and plotting on a graph the limiting values of the principal strains for which the rupture of the test sample is observed [2]. Thus, one obtains a Forming Limit Curve which characterizes the
206
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
limit of deformation of the sheet. The FLC curve is supposed intrinsic for each sheet and used to predict the risk of appearance of defects (cracks) during any forming process. It allows only to obtain, at each time step, the contours of plastic deformation in the sheet at any process simulation stage. The zones in which the values of the limiting deformations reach the experimental ones, are supposed to be the regions of local rupture [3]. It is clear that this method is not of a predictive nature. It can only be used to compare sheets between them in terms of their forming properties. In particular, it does not make it possible to predict when and where the cracks can appear in the workpiece during the forming operation. Some research works have been published on the damage models based either on continuum damage mechanics [1,4,5], or on Gurson’s theory [6–9], among others. In the present work, we develop and evaluate a damage model, based on the strong coupling between anisotropic elastoplasticity with mixed nonlinear work hardening (isotropic and kinematic) and an isotropic ductile damage. The model is formulated using the theory of engineering plasticity coupled with damage within the thermodynamic framework of irreversible processes [10,11]. Tensile tests carried out on aluminium samples are used to identify the parameters of the model. At present, this paper deals with isotropic damage. Anisotropic damage issue is another difficult subject for future research work. In this paper we, first, recall the main aspects of the formulation of the plasticity-damage model. Then we apply the model in the context of the stamping process via a Vumat user material in the framework of ABAQUS/Explicit finite element code. The data are taken from the Swift’s deep-drawing test. Finally, some discussions and comparisons with experiments are drawn. The obtained results show its ability to predict when and where the ductile damage can appear in the workpiece during the deep-drawing operation. 2. Plasticity-damage model
the isotropic hardening and its associated variable, (D, Y) represents the isotropic ductile damage and its associated variable, K is the fourth-order operator of elastic properties. C and Q are respectively the kinematic and isotropic hardening parameters. The mechanical dissipation is described using the theory, which considers only one potential of dissipation F. It controls both the plastic flow f(r, X, R, D) and the evolution of the damage F(r, X, R, D), with one plastic multiplier k_ [13–15]. kr X k R pffiffiffiffiffiffiffiffiffiffiffiffi ry ¼ 0 ð5Þ 1D 3 a 1 b R2 F ðr; X ; R; DÞ ¼ f þ X :Xþ 4 Cð1 DÞ 2 Qð1 DÞ sþ1 1 S Y Y0 þ ð6Þ b S ð1 DÞ s þ 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with kr X k ¼ ðr X Þ : H : ðr X Þ. ry is the yield
f ðr; X ; R; DÞ ¼
strength. a and b are, respectively, the kinematic and isotropic nonlinear hardening parameters. S, Y0, b and s describe the isotropic damage evolution; while H is the fourth-order operator of Hill [10,16,17] as a function of six constants F, G, H, L, M and N [10,18]. Using the normality rule, complementary relations are derived as follows [10,19]: H : ðr X Þ k_ ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ~n e_ p ¼ k_ pffiffiffiffiffiffiffiffiffiffiffiffi 1D 1 Dkr X k p _ a_ ¼ e_ kaa 1 r_ ¼ k_ pffiffiffiffiffiffiffiffiffiffiffiffi br 1D s k_ Y Y0 _D ¼ b S ð1 DÞ
ð7Þ ð8Þ ð9Þ ð10Þ
where ~n is the normal to the loading surface. The plastic multiplier is determined using the consistency condition f ¼ f_ ¼ 0.
2.1. Theoretical aspects
2l f_ ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ~n : e_ hk_ ¼ 0 1D
The plasticity damage model presented here is based on the thermodynamic approach of irreversible processes with internal variables. The state relations are given hereafter as follows [12]:
where h > 0 is the tangent elastoplastic modulous defined such that: pffiffiffiffiffiffiffiffiffiffiffiffi R h ¼ 2l 1 D þ C ~n : ~n þ Q b pffiffiffiffiffiffiffiffiffiffiffiffi 1D s a 1 Y Y 0 hry i pffiffiffiffiffiffiffiffiffiffiffiffi ~n : X þ ð12Þ S 2 1D ð1 DÞbþ1
r ¼ ð1 DÞK : ee 2 X ¼ Cð1 DÞa 3 R ¼ ð1 DÞQr 1 1 1 Y ¼ ee : K : ee þ Ca : a þ Qr2 2 3 2
ð1Þ ð2Þ ð3Þ ð4Þ
where the couple (ee, r) represents the tensors of elastic strain and Cauchy stress, (a, X) represents the kinematic hardening and its associated variable. (r, R) represents
ð11Þ
The formulation of the above model in large deformations is done using the rotating referential. The rotation Q of the deformed configuration is at each time solution of a kinematic equation [20,21] _ T ; Qðt ¼ 0Þ ¼ I W Q ¼ QQ where I is the identity tensor of order 2.
ð13Þ
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
The elastic strains are assumed infinitesimal in order to ensure the summation of the strain rates D ¼ e_ Je þ Dp
ð14Þ
where e_ Je is the Jaumann’s rotational derivative of the infinitesimal elastic strain. 2.2. Numerical aspects The finite element method is used to discrete the global equilibrium equation expressed as follows:
€ þ fF int g fF ext g ¼ 0 ½M U ð15Þ where [M] is the mass matrix. {Fint} and {Fext} are, respectively, the internal and the external force vectors. The above equation is solved explicitly. The numerical scheme consists in obtaining solutions at each time t+Dt as a function of the known quantities at time t
€ n ¼ ½M1 ½ffext g ffint g U ð16aÞ
Dtnþ1 þ Dtn
€n U_ nþ1=2 ¼ U_ n1=2 þ ð16bÞ U 2
ð16cÞ fU nþ1 g ¼ fU n g þ Dtn þ 1 U_ nþ1=2 The explicit scheme does not require iterating for the determination of the tangent matrix. Its stability and accuracy depend strongly on the size of the time step Dt [22]. The discretization of Eqs. (7)–(10) is associated with an asymptotic integration of Eqs. (8) and (9). Then, the incremental equations become: Dk Dpnþ1 ¼ Dpn þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ nnþ1 1 Dnþ1 ð1 eaDk Þ anþ1 ¼ an eaDk þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ nnþ1 a 1 Dnþ1 ð1 ebDk Þ rnþ1 ¼ rn ebDk þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b 1 Dnþ1 s Y nþ1 Y 0 1 Dnþ1 ¼ Dn þ Dk b S ð1 Dnþ1 Þ
ð17Þ ð18Þ ð19Þ ð20Þ
The constitutive equations need to be integrated over each time to obtain the values of the variables rn+1, Xn+1, Rn+1 and Yn+1. The integration of the constitutive equations allows to reduce the number of equations from (15) to (8) [11,16,23]. The system of equations to solve becomes then fnþ1 ðDk; ~ nnþ1 ; Dnþ1 Þ Rnþ1 krnþ1 X nþ1 k ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ry ¼ 0 1 Dnþ1 1 Dnþ1 hnþ1 ðDk; ~ nnþ1 ; Dnþ1 Þ ¼ H : ðrnþ1 X nþ1 Þ krnþ1 X nþ1 k~ nnþ1 ¼ 0 nnþ1 ; Dnþ1 Þ gnþ1 ðDk; ~ ¼ Dnþ1 Dn Dk
s Y nþ1 Y 0 1 ¼0 S ð1 Dnþ1 Þb
ð21aÞ
ð21bÞ
ð21cÞ
207
with krnþ1 X nþ1 k ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðrnþ1 X naþ1 Þ : H : ðrnþ1 X nþ1 Þ
rnþ1 ¼ rn þ ð1 Dnþ1 Þ½2lDe þ ke trðDeÞ1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2G 1 Dnþ1 Dk~nnþ1 ð1 eaDk Þ aDk X nþ1 ¼ Cð1 Dnþ1 Þ an e þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~nnþ1 a 1 Dnþ1 ð1 ebDk Þ bDk þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Rnþ1 ¼ Qð1 Dnþ1 Þ rn e b 1 Dnþ1
ð22Þ
ð23Þ ð24Þ ð25Þ
where G is the shear modulous; ke and l are the Lame constants. The unknown quantities are Dk; ~nnþ1 and Dnþ1 . The solution of the system of equations is obtained with the Newton–Raphson method [11]. 2.3. Implementation within ABAQUS The determination of the quantities rn+1, Xn+1, Rn+1 and Dn+1 at the time tn+1, is obtained using the radial return algorithm proposed by Simo [23]. (i) Elastic prediction: Given a total strain increment De at the time tn, initially the response is assumed to be purely elastic; that is, the stress is obtained as rnþ1 ¼ rn þ ð1 Dn Þ½2lDe þ ke trðDeÞ1
ð26Þ
Substituting the value of rnþ1 in the expression of the criterion results in r Rn ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ry fnþ1 ¼ pnþ1 ð27Þ 1 Dn with qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ r : H : r nþ1 nþ1 nþ1
ð28Þ
X n ¼ Rn ¼ 0
ð29Þ
if fnþ1 < 0, the mechanical response is purely elastic. The stress at the time tn+1 can be written as follows:
rnþ1 ¼ rnþ1
ð30Þ
Dknþ1 ¼ Dkn ; ~nnþ1 ¼ ~nn ; Dnþ1 ¼ Dn
ð31Þ
Table 1 Material parameters of the used aluminium alloy (AL 5754) Elasticity E (MPa) 70000
m 0.3
Plasticity ry (MPa) 114.1 F 0.89
C (GPa) 2.1 G 0.67
a 21.1 H 0.33
L 1.5
Damage S (MPa) 17
s 13
b 9
Y0 (MPa) 0.01
Q (MPa) 554.2 M 1.5
b 3.67 N 1.65
208
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
(ii) Plastic-damage behaviour (radial correction): However, if fnþ1 > 0, the Newton–Raphson scheme is used to obtain the unknowns parameters Dk; ~ n and D of the system (Eq. (21))
8 > < fnþ1 ðDk; ~nnþ1 ; Dnþ1 Þ ¼ 0 hnþ1 ðDk; ~nnþ1 ; Dnþ1 Þ ¼ 0 > : gnþ1 ðDk; ~nnþ1 ; Dnþ1 Þ ¼ 0
Fig. 1. Equivalent plastic strain and damage distributions for the rolling direction (0°).
ð32Þ
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
Successive corrections are used to tend Eq. (21) to be equal to zero, which results in updated values of the unknowns Dk; ~n and D at the time tn+1. When convergence occurs,
209
the solution is stored and the variables rn+1, Xn+1, Rn+1 and Dn+1 are updated.
3. Numerical applications 3500
3.1. Model parameters identification: tension test
3000
To identify the model parameters, a tensile test was carried out on an aluminium sample. The geometrical characteristics of the tensile specimen are 1 mm thick, 12.52 mm width and 180 mm long. The experiment is also simulated using ABAQUS/Explicit, by considering only the useful part. The finite element mesh consists of 2366 eight-nodes hexahedral elements. The part is meshed as follows: 91 elements along the length, 13 elements in the width direction and 2 elements in the thickness one. The simulation is carried out using a displacement rate of 10 mm/s and by considering a density of 2700 kg/m3. The different material properties and parameters are listed in Table 1. Fig. 1 shows the distributions of effective plastic strain and damage for different displacements in the direction of rolling (0°). For this direction, it appears clearly that there are three stages of damage evolutions. Initially, for a displacement less than u = 16.8 mm, the strain distribution is homogeneous as shown in Fig. 1a. Beyond this limit, the deformation becomes heterogeneous. A highly strained
Load [N]
2500
2000
Experiment Numerical
1500
1000
500
0 0
2
4
6
8
10
12
14
16
18
20
22
Displacement [mm]
(a) : Direction 0° 3000
2500
Load [N]
2000
Experiment Numerical
1500
Table 2 Maximum loads and displacements recorded at failure 0°
1000
Experiment Numerical Error (%)
500
45°
90°
u (mm)
F (daN)
u (mm)
F (daN)
u (mm)
F (daN)
14.9 19.8 32.7
290.5 302.3 4.1
18.4 19.6 6.8
281.6 286.6 1.7
19.1 19.8 3.9
279.8 284.2 1.5
0 0
2
4
6
8
10
12
14
16
18
20
22
Displacement [mm]
(b) : Direction 45° 3000
2500
Load [N]
2000
Experiment Numerical
1500
1000
500
0 0
2
4
6
8
10
12
14
16
18
20
22
Displacement [mm]
(c) : Direction 90° Fig. 2. Load–displacement curves.
Fig. 3. Schematic illustration of the swift process tools.
210
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
zone appears in the middle of the specimen at a displacement of 18.6 mm (Fig. 1b). This is followed by the quick formation of a constriction in the form of two shear bands, before a macroscopic crack starts at the intersection of the two bands at a displacement of 19.8 mm. The value of the damage (SDV21) reaches the maximum value D = 1 in the cracked zones (Fig. 1c). It can be seen that the rupture facies is perpendicular to the longitudinal axis of the sample. The same observations are also made in the 45° and 90° directions with respect to the direction of rolling. The damage evolution also follows three stages.
The obtained results in each direction (0°, 45° and 90°) are also reported in the form of load–displacements curves as shown in Fig. 2. Except for a very slight deviation above a displacement of 6 mm, the obtained results agree very well with their experiment counterparts. The failure loads in the 0° and 45° directions are not as well predicted as for the 90° direction. The reason is that the damage considered here is isotropic. Table 2 shows the recorded values for load and displacement at failure for each direction. Again, it can be seen that the failure load and the corresponding displacements are as well predicted for the 45° and 90° directions.
Fig. 4. Equivalent plastic strain and damage distributions for different punch displacements (travels).
Fig. 5. Predicted sheet shapes for different punch travels (0, 6.63, 9.12, 10.64 and 11.3 mm.
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
After identification of the model parameters, the Swift’s deep-drawing test was used to examine the capability of the developed damage model presented above. This test was also reproduced experimentally at the CETIM of Senlis (France). The schematic deep-drawing test is shown in Fig. 3. The blank is a circular sheet of 73 mm of diameter and 1 mm thick. The blank is meshed using eight-node hexahedral elements with selective reduced integration, while the other tools (punch, die, blank-holder) were regarded as rigid bodies. Coulomb’s friction law was assumed and l = 0.17 was used for the friction coefficient between the sheet and the tools. Fig. 4 shows the distributions of effective plastic strain and damage as a function of the punch displacement. At the beginning (punch travel of u = 7.5 mm), the damage was located along a circumference zone which comes in contact with the punch radius. It can be seen that the damage increases in the circular zone, which decreases in diameter, as the punch displacement increases. A strong localization of damage was induced at a displacement of 10.3 mm, which quickly transformed into macroscopic cracks at a displacement of 11.1 mm.
Fig. 5 presents the predicted shapes of the sheet for different punch travels u = 0, 6.63, 9.12, 10.64 and 11.3 mm. Fig. 6a and b shows a comparison between the experimental rupture and that predicted using the developed damage model. In a deep-drawing test, the risk of rupture is evaluated by using the Forming Limit Curves (FLC) defined in the plane of principal strains as shown in Fig. 7. The strain values at different depths of the punch are measured in order to estimate the likelihood of crack nucleation in the high risk zones. The experimental and numerical points are drawn from the framed zone represented in Fig. 6b. The points located above the FLC curve represent the rupture limit of the sheet. Fig. 8 shows comparison between the experimental and numerical results with regard to the punch force-travel.
0.50
Limit curve Experimental points Numerical points
0.45 0.40
Major strain
3.2. Application to the deep-drawing: Swift’s test
211
0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00
-0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
Minor strain Fig. 7. Forming Limit Curves FLC.
30000
Punch load [N]
25000
Experiment Numerical
20000
15000
10000
5000
0 0
2
4
6
8
10
12
Punch displacement [mm]
Fig. 6. Fracture in the workpiece: comparison between simulation and experiment.
Fig. 8. Punch force–displacement curves: comparison between simulation and experiment.
212
M. Khelifa et al. / Computers and Structures 85 (2007) 205–212
The maximum computed punch force is equal to 25.93 kN and it was found in good agreement with the experimental value of 25.92 kN. The rupture was produced at the same time (for punch travel of u = 11.3 mm) in both experimental and simulation procedures. 4. Conclusion A predictive damage model, based on the theory of damage-plasticity, is proposed and successfully implemented in ABAQUS/Explicit finite element code. It allows to predict the onset of macroscopic cracks in thin sheet metals during deep-drawing operations. The experimental FLC method is a post processing at the end of simulation stage. It neglects the interaction between plasticity and damage. The combination of plasticity and damage allows essentially to avoid the over-estimation of the formability of sheet materials. The proposed damage model gives good results in comparison with experiments. The present procedure allows to predict the appearance of ductile rupture and its growth during deep-drawing operations. Acknowledgements CETIM of Senlis and FEDER of ‘‘l’Etat de la region Champagne–Ardenne’’ are gratefully acknowledged for their contributions and financial supports, respectively. References [1] Guo YQ, Li YM, Bogard F, Debrey K. An efficient pseudo-inverse approach for damage modelling in the sheet metal forming process. J Mater Process Technol 2004;151:88–97. [2] Sfar H, et Maillard A. Applications de la simulation de l’emboutissage, CETIM, Rap. No1B8422, De´cembre 2000. [3] Sfar H, et Maillard A. Aides a` la re´alisation de l’emboutissage: simulation et mesure automatique des de´formations, CETIM, Recueil de confe´rences de la journe´e d’information du, 17 Novembre 1998. [4] Lee H, Peng K, Wang J. An anisotropic damage criterion for deformation instability and its application to forming limit analysis of metal plates. Eng Frac Mech 1985;21(5):1031–54.
[5] Zhu YY, Cescotto S, Habraken AM. A fully coupled elastoplastic damage modelling and fracture criteria in metal forming processes. J Mater Process Technol 1995;32:197–204. [6] Gelin JC, Oudin J, Ravalard Y. An imposed finite element method for the analysis of damage and ductile fracture in cold metal forming processes. Ann CIRP 1985;34(1):209–13. [7] Onate E, Kleiber M. Plastic and viscoplastic flow of void containing metal applications to axisymmetric sheet forming problems. Int J Num Meth Eng 1988;25:237–51. [8] Brunet M, Sabourin F, Mghib-Touchal S. The prediction of necking and failure in 3D sheet forming analysis using damage variable. J Phys III 1996(6):473–82. [9] Boudeau N, Gelin JC. Necking in sheet metal forming: influence of macroscopic and microscopic properties of materials. Int J Mech Sci 2000;42:2209–32. [10] Khelifa M, Badreddine H, Belamri N, Gahbiche M-A, Saanouni K, Cherouat A, et al. Int J Form Process 2005;8(N°2). [11] Khelifa M. Simulation nume´rique de l’endommagement en formage des structures minces, The`se de Doctorat, Universite´ de Technologie de Troyes, France, 2004. [12] Saanouni K, Forster C, BenHatira F. On the inelastic flow with damage. Int J Damage Mech 1994;3(2):140–69. [13] Chaboche J-L. Le concept de contrainte effective applique´ a` l’e´lasticite´ et a` la viscoplasticite´ en pre´sence d’un endommagement anisotrope, In: Boehler J-P, editor. Colloque EuroMech 115, Grenoble, Editions du CNRS, 1979, p. 737–60. [14] Lemaitre J. A Continuum damage mechanics model for ductile fracture. J Eng Mat Tech 1985;107:83–9. [15] Lemaitre J. A course on damage mechanics. Springer Verlag; 1996. [16] Chaboche J-L, Cailletaud G. Integration methods for complex plastic constitutive equations. Computer methods in applied mechanics and engineering, vol. 133. Elsevier; 1996. p. 125–55. [17] Hill R. A theory of yielding and plastic flow of anisotropic metals. Royal Soc Lon Proc 1948:281. [18] Lemaitre J, Chaboche JL. Me´canique des mate´riaux solides, Edition Dunod, 1988. [19] Saanouni K, Chaboche J-L. Computational damage mechanics. In: Milne I, Ritchie RO, Karihaloo B, editors. Comprehencive structural integrity. de Borst R, Mang HA, editors. Application to metal forming, numerical and computational methods, vol. 3, 2003, ISBN: 0-08-043749-4 (Chapter 7). [20] Sidoroff F. The geometrical concept of intermediate configuration and elastic–plastic finite strain. Arch Mech 1973;25(N°2):299–308. [21] Dogui A. Plasticite´ anisotrope en grandes de´formations, The`se de doctorat d’e´tat, Universite´ Claude Bernard, lyon I, 1989. [22] ABAQUS, Theory manual, Version 6.4, Hibbit, Karson and Sorensen, Inc., 2004. [23] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998.