Fethi Abbassi
email: fethi.abbassi@ipeib.rnu.tn

Olivier Pantalé
email: pantale@enit.fr

Sébastien Mistou
email: mistou@enit.fr

Ali Zghal
email: ali.zghal@esstt.rnu.tn

Roger Rakotomalala

Effect of ductile damage evolution in sheet metal forming: experimental and numerical investigations

Keywords: Ductile damage, Fracture, GTN model, Necking, Sheet metal forming, DIC, Numerical simulation

The numerical simulation based on the Finite Element Method (FEM) is widely used in academic institutes and in the industry. It is a useful tool to predict many phenomena present in the classical manufacturing forming processes such as necking, fracture, springback, buckling and wrinkling. But, the results of such numerical model depend strongly on the parameters of the constitutive behavior model. In the first part of this work, we focus on the traditional identification of the constitutive law using oriented tensile tests (0°, 45°, and 90° with respect to the rolling direction). A Digital Image Correlation (DIC) method is used in order to measure the displacements on the surface of the specimen and to analyze the necking evolution and the instability along the shear band. Therefore, bulge tests involving a number of die shapes (circular and elliptic) were developed. In a second step, a mixed numerical-experimental method is used for the identification of the plastic behavior of the stainless steel metal sheet. The initial parameters of the inverse identification were extracted from a uniaxial tensile test. The optimization procedure uses a combination of a Monte-Carlo and a Levenberg-Marquardt algorithm. In the second part of this work, according to some results obtained by SEM (Scaning Electron Microscopy) of the crack zones on the tensile specimens, a Gurson Tvergaard Needleman (GTN) ductile model of damage has been selected for the numerical simulations. This model was introduced in order to give informations concerning crack initiations during hydroforming. At the end of the paper, experimental and numerical comparisons of sheet metal forming applications are presented and validate the proposed approach.

Introduction

The numerical simulation is a widely used tool to optimize the sheet metal forming process and to reduce the development time and therefore the final cost of products in manufacturing such as, automotive and aerospace. An accurate estimation of the material parameters for constitutive models is often a common need for the engineering community [START_REF] Andrade-Campos | On the determination of material parameters for internal variable thermoelastic-viscoplastic constitutive models[END_REF]. With the growth of computing facilities, material models are becoming more and more complex, leading to the need of more complete identification procedures. The identification step precedes the design of the manufacturing process. The results of the identification procedure are used as input data for the numerical simulation of the forming process. Finite Element Modeling (FEM) has been widely introduced into the design of manufacturing products because of its high efficiency in predicting several problems and major defects occurring in sheet metal forming manufacturing process like necking, wrinkling, buckling and surface deflections. In order to predict those defects within a virtual manufacturing system (numerical simulation), an accurate description of the sheet metal mechanical behavior is an essential requirement. In order to increase the predictive capabilities of the virtual forming tools, an accurate modeling of microdefects (i.e. damage) initiation and growth under large transformations have to be used. Concerning, the ductile failure of materials, it is controlled by three stages: micro-void nucleation, growth and coalescence mechanisms (Fig. 1). The damage evolution depends directly on the state of deformation. An explanation of the internal degradation and a review of the mathematical modeling, and the recent development in numerical tool can be found in many publications [START_REF] Chaboche | Continuous Damage Mechanics: a tool to describe phenomena before crack initiation[END_REF][START_REF] Brünig | Nonlocal continuum theory of anisotropically damaged metals[END_REF][START_REF] Rousselier | simplified polycrystalline model for viscoplastic and damage finite element analyses[END_REF][START_REF] Zhiying | The GTN damage model based on Hill'48 anisotropic yield criterion and its application in sheet metal forming[END_REF]. In this investigation, according to the metallographic observations and the spherical inclusion in the used material, a ductile micromechanical damage model proposed by Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth. I. Yield criteria and flow rules for porous ductile media[END_REF][START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] was selected.

Experimental technique and results

Digital Image Correlation in 2D: Full-field optical measuring methods

The image correlation allows the measurement of the displacements fields on a planar surface: a single camera acquires a sequence of images of a planar object under plane strain or stress during the deformation process. The displacements of many points, distributed on the surface of the object are calculated from the grey level analysis of the images. DIC is an application based on the comparison of two images acquired at different states during the deformation. Two subsets are chosen respectively from the reference and the next image used for the computation. The algorithm of correlation used to detect the local displacements of the pixels by comparing the two subsets, is as follows. In practice, a single value is not a unique signature of a point, hence neighboring pixels are used. The matching of images acquired by only one camera, at different time, on an object which becomes deformed is called temporal matching, or tracking. From its principle, the correlation technique can be used correctly only with objects having a surface with a sufficiently random texture.

One CCD (Coupled Charge Device) camera acquires digital images with a 1280 by 1024 pixels definition. These images are then analyzed by the Aramis software developed by the GOM society. This one gives a measure of the displacements and strain fields on the surface of an object within a precision of 200 µm/m. Aramis software is able to match correspondent points of an image by grey level analysis, if the surface of the specimen is covered with a black and white mapping, which forms a random grey-level (see Fig. 2). A gray level coded on 8 or 12 bits corresponds to each pixel of the CCD sensor. A succession of several pixels lying on the same line forms a grey level sequence. The fundamental principle lies in the fact that the distribution does not vary during the deformation of the object. It is therefore sufficient to follow this distribution of grey levels during their displacement to obtain the displacements of the corresponding point. This is made possible by the use of correlation domains which are nn pixel gathering zones.

Fig. 2 The variable grey levels on the surface of a tensile specimen

In the DIC method the relation between the deformed image and the undeformed one is illustrated in Fig. 3.

If we call G the central point of the subset in the initial configuration and G* the corresponding central point of the subset in the deformed configuration, the relationship between the coordinates of these two points can be expressed by:

) y (x, v y y y) (x, u x x * *     (1)
where u(x,y) and v(x,y) represent the homogenous displacement field for one pattern (whole set of pixels). Fig. 3 Relative location of sub-images of the deformed and the undeformed patterns on the surface.

In the Aramis software the displacement field is considered homogeneous and bilinear with respect to x and y:

                   4 3 2 1 4 3 2 1 ) , ( ) , ( b y x b y b x b y x v a y x a y a x a y x u (2)
Where a 4 and b 4 represent the terms of the rigid body motion, a 1 , b 1 , a 2 and b 2 represent the elongation terms and a 3 and b 3 are the shearing terms.

Tensile test

In the first part of this work, a large experimental work has been done. Several low rate tensile tests were carried out in the laboratory on an INSTRON tensile test machine at a strain rate of 10 -3 /s. 12 tensile specimens were obtained using a laser cutting process from the sheet at 0°, 45° and 90° with respect to the rolling direction with standard ASTM-E8. A digital image correlation method was used to measure the displacement fields at the surface of the specimens during the tensile tests. Figure 4 illustrates the specimen's geometry and the equipments used. During the tensile test, many data are collected such as the hardening law, the anisotropic yield criterion and the necking of the material.

Bulge test

The standard tensile test is one of the most important engineering procedures used to characterize the mechanical behavior of materials [START_REF] Celentano | Characterization of the mechanical behaviour of materials in the tensile test: experiments and simulation[END_REF]. As the range of stable uniform strain is restricted to less than half the one sustainable under biaxial stress [START_REF] Michel | Size effects on the constitutive behaviour for brass in sheet metal forming[END_REF], it is obviously desirable to generate the required data directly from biaxial tests. The circular hydraulic bulge test is widely used to determine the biaxial stress-strain curve for sheet metal forming. This test has been used by many authors in order to study the material formability and to illustrate several defects of sheet metal forming process predicted by numerical simulation softwares [START_REF] Rui | One Point quadrature shell elements for sheet metal forming analysis[END_REF].

During this test (see Fig. 5 for geometry and parameters), a circular blank is clamped at its external boundary between a die and blank holder by a drawbeads and a linearly increasing hydraulic pressure is applied on the bottom surface of the blank. Both circular and elliptic dies have been used for those tests. The experimental tests have been performed on a 304L stainless steel. Many experiments have been done in order to study the influence of the rolling direction vs. the orientation of the major axis of the elliptic die (again, at 0°, 45° and 90°). The procedure of a parametric identification problem consists in computing many numerical simulations of a physical process, knowing a set of experimental results, and then to minimize the difference between the set of numerical and experimental responses. The solution of the problem is then composed by the set of parameter for which the minimal difference has been obtained during the process. Concerning the identification procedure, many optimization methods have been developed in order to identify non linear physical laws (behavior of material, damage, necking…). The optimization methods can be classified in three main categories [START_REF] Kleinermann | Parameter identification and shape/process optimization in metal forming simulation[END_REF]: classical structural optimization methods, linearization methods (Gauss-Newton, Levenberg-Marquardt) and convex approximations methods (sub-iterative Globally Convergent Method of Moving Asymptotes (GCMMA)). In this contribution, we used a homemade tool developed in the LGP laboratory called "PILOTE". This one uses a combination of the Monte-Carlo and the Levenberg Marquardt algorithms [START_REF] Abbassi | Parametric sheet metal characterization by using Monte-Carlo and Levenberg-Marquardt: bulge test application[END_REF]. The Monte-Carlo method is a widely used tool based on a random number generation and a set of statistical tools. It has been successfully applied in many fields such as mechanics, electrics or nuclear physics. In this identification procedure, the idea is to use such a factor of chance to solve a physical problem using a set of random parameters as input data for the numerical model. The application of the Monte-Carlo method requires defining a set of initial parameters and the allowed range for each parameter. This method is a good tool for a coarse research. Once the convergence criterion of the Monte-Carlo is satisfied the obtained set is used as input data for the refinement of the solution by the Levenberg-Maraquardt algorithm. The Levenberg-Marquardt optimization procedure is also widely used for nonlinear optimization. It's a stabilized version of the Gauss-Newton method, which works very well in practice and has become the standard for solving nonlinear least squares system of equation.

For the identification, the general procedure is schematically illustrated in Fig. 12. The experimental results obtained by tensile test are injected into the identification algorithm as starting parameters for the identification procedure. The objective function can be defined as in equation 3 to minimize the difference between the numerically computed and experimentally measured displacement of the centre point of the blank. This dome position h can be recorded during the whole bulge test.

            m j FEM r j h j h j h j w m F 1 2 exp exp ] [ ] [ ] [ ] [ 1 (3)
where m is the total number of response, h FEM is the vector of the simulated responses, h exp is the vector of the experimental response and w r is the vector of the responses weights. In this investigation, an isotropic hardening is considered in the numerical model to describe the stress-strain relation of AISI 304 metal. A subroutine (Uhard) has been implemented in the Abaqus. This relationship, a Swift law, is described by the equation:

(4)

Experimental observations

Onset of ductile fracture is initiated by void formation around non-metallic inclusions and second-phase particles in metal matrix that is subjected to plastic strain under influence of external loading. Observations of the damage nucleation mechanisms have been carried out through in situ tensile tests at room temperature. The spherical shape of the second phase particles and non-metallic inclusions present in the stainless steel can be identified as possible spots for the initiation of the ductile fracture (Fig. 7). Several damage models were developed to describe this complex micro-mechanism in sheet metal used in forming process.

Gurson-Tvergaard-Needleman (GTN) Model

The main purpose of this work is to use Continuum Damage Mechanics (CDM) to predict a sheet metal failure. In 1977, Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth. I. Yield criteria and flow rules for porous ductile media[END_REF] developed a constitutive model to describe the micro-growth nucleation mechanical effects of damage in ductile metals, which was extended in 1984 by Tvergaard and Needleman [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] to incorporate some additional parameters (q 1 , q 2 , q 3 ). The yield function describing the plastic constitutive model is represented as follows:

0 1 2 3 cosh 2 2 3 2 1 2 2              ) ) (f q ( σ σ q q f σ σ φ * y m * y (5)
Where q 1 , q 2 and q 3 is a set of fitting parameters used to calibrate the model predicting the periodic array of spherical and cylindrical voids. In this equation, q 3 = (q 1 ) 2 , σ m is the hydrostatic stress, σ y is the yield stress of the material, is the von Mises equivalent stress and f is the void volume fraction which is equal to 1 . The term represent the modified void volume fraction given as follows:

            c c c c * u c c * f f if ) f (f f f f f f f f if f f f (6)
Where is the ultimate value of , is the void volume fraction corresponding to failure, and is the critical volume fraction when rapid coalescence occurs. The instantaneous rate of growth of the void fraction depends both on nucleation of new voids and on the growth of pre-existing voids. It is given by:

          nucleation p growth p kk d a d f df       ) ( ) 1 ( (7)
Where is the plastic hydrostatic strain and "a" is defined as follows:

Spherical inclusions                       2 2 1 exp 2 N N p N N S S f a    (8)
Where is the standard deviation, is the mean effective plastic strain of nucleation and is the nucleation micro-void volume fraction.

Numerical investigations

The objective of the numerical simulations is to predict the damaged zones of the forming parts. The parameters of the GTN model have been taken from literature data [START_REF] Abbassi | Fédération des connaissances de mise en forme dans une plateforme de prototypage virtuel[END_REF][START_REF] Benseddiq | A ductile fracture analysis using a local damage model[END_REF] and the comparison of finite element prediction and experimental results was used to determine the critical damage parameters. Table 2 reports the material parameters used in the all numerical models used hereafter.

Parameters

Element

Bulge test

The GTN damage model was used in the simulation of the elliptical bulge test and compared with the experimental results. Fig. 9 (a) and Fig. 9 (b) present respectively the numerical results and the experimental deformed part after fracture. We observe a good correlation between both results. As the numerical model doesn't include crack propagation, only the initiation of the crack is reported by the model. represent respectively the evolution of the dome position according to the internal pressure of the elliptical bulge test with the major axes oriented 90°,45° and 0° with the rolling direction of sheet. In table 3, a quantitative comparison those results for three different levels of pressure in the bulge test (P=10 MPa, P=20MPa and at fracture) is reported. A quite small difference between those results is observed, mainly linked to the fact that the GTN model is isotropic. The thickness (STH) variation is commonly used as an indicator for the appearance of defects during the forming process. For these reasons, it's recommended to evaluate the thickness variation during the development of a manufacturing process. Gutscher founds that the behavior law parameters has significant influence on the height and thickness in the bulge test [START_REF] Gutscher | Determination of flow stress for sheet metal forming using the viscous pressure bulge (VPB) test[END_REF]. In order to evaluate the identified parameters obtained with the parametric method, we have measured the thickness of the specimen using an optical system based on the use of a fringe pattern projection. The measures have been compared with the numerical results obtained by the numerical modeling of the elliptical bulge test based including the ductile damage.

Results are reported in Fig. 11, and, again, a good correlation was observed between those results.

Erichsen test

Erichsen test, consisting in the deep drawing of a circular sheet using a spherical punch (Fig. 12), provides the so-called Erichsen index (IE) obtained as the punch penetration before rupture [START_REF] François | Techniques de l'Ingénieur -Essais mécaniques des métaux. "Essais d'aptitude à la mise en forme[END_REF]. This instant can be associated to the maximal punch force. The highly localized thickness reduction that takes place beyond this point causes, for slightly higher levels of punch penetration, a sudden circumferential rupture of the specimen. Fig. 13 schematically depicts the problem layout. In the numerical investigation, It is noted that GTN model it is an effective tool to predict the rupture of the specimen.

Conclusions

This paper presents an experimental and a numerical investigation of porous metal plasticity. Physical mechanisms of ductile fracture of stainless steel sheet, nucleation, growth and coalescence of voids have been presented. In this work, the Gurson-Tvergaard-Needleman (GTN) model is used to study the damage evolution and necking in sheet metal forming process. In the experimental section, several tensile tests were presented in order to evaluate the heterogeneity of deformation and the diffuse and local necking processes.

Experimental results have been obtained with optical measurement techniques based on Digital Image Correlation (DIC) and 3D digitalization.

A numerical approach has been proposed based on the use of the GTN model. The main result presented here concerns the experimental measurement of the thickness variation of the whole deformed part during an elliptic bulge test and the comparison with the numerical results. A good correlation between those results was observed, therefore, we can conclude from this that the proposed numerical approach is quite reliable.

  for one pattern a 1 , b 1 , a 2 , b 2 Elongation terms a 3 and b 3 Shearing terms.

Fig. 1 :

 1 Fig. 1: (a) tensile test of bi-entailed plate specimen, (b) evolution of the damage process in a ductile metal related to the macroscopic loading evolution.

Fig. 4

 4 Fig. 4 Experimental setup of uni-axial tensile test

Fig. 5

 5 Fig. 5 Bulge test geometry and parameters

Fig. 6

 6 Fig.6 Schema of parameter optimization by Monte-Carlo and Levenberg-Marquardt algorithm

Fig. 7

 7 Fig. 7 Scanning electron microscope (SEM) views of fracture in tensile specimen (spherical inclusion observations)

Fig. 8a presentsFig. 8 :

 8 Fig.8apresents the distribution of the equivalent plastic strain (PEEQ). These numerical results are obtained with the GTN damage model in a tensile test simulation using the ABAQUS code. A non uniform distribution is clearly obtained due to the localization of the deformation. The load drop due to the presence of a diffuse necking and the localization is clearly captured by the use of the ductile damage model. Fig.8bshows the experimental plastic deformation field just before the rupture of the specimen obtained from the optical measurement system. The good correlation between both results shows the accuracy of the numerical prediction of the shear band.

Fig. 9 (

 9 Fig.9(a) Numerical prediction and (b) experimental specimen after cracking for elliptical bulge test case Fig.10reports the numerical and the experimental evolution of the dome position during the bulge test according to the pressure load. The 3 experimental curves Exp_Ellip_90, Exp_Ellip_45 and Exp_Ellip_0 represent respectively the evolution of the dome position according to the internal pressure of the elliptical bulge test with the major axes oriented 90°,45° and 0° with the rolling direction of sheet. In table 3, a quantitative comparison those results for three different levels of pressure in the bulge test (P=10 MPa, P=20MPa and at fracture) is reported. A quite small difference between those results is observed, mainly linked to the fact that the GTN model is isotropic.

Fig. 10

 10 Fig. 10 Numerical and experimental comparison of dome position evolution according to the pressure load

Fig. 11 .

 11 Fig. 11. Thickness variation in deformed parts (a) numerical results, (b) experimental results

Fig. 12 :

 12 Fig. 12: Erichsen test (a) Geometrical condition [16], (b) numerical model

Fig. 13 :

 13 Fig. 13: Numerical results of Erichsen test (a) Evolution of volume fraction, (b) von Mises stress.

Table 1

 1 reports the results of the parametric identification. Those parameters will be used in the numerical model of the forming process.

		Swift Law coefficients
		0 ε	K (Mpa)	n
	Initial Parameters	0.08	1538	0.6
	Identified Parameters
		0.0368	1580.53	0.53

Table 1 :

 1 Results of parametric identification process

3.1. Continuum Damage Mechanics (CDM)

Table 3 .

 3 Comparison between numerical an experimental results in elliptical bulge test (direction 0°)