Calibration method of ductile damage model based on hybrid experimental-numerical analysis of uniaxial tensile and hole-expansion tests

Accepted Manuscript Calibration method of ductile damage model based on hybrid experimentalnumerical analysis of uniaxial tensile and hole-expansion tests Hela Soussi, Abdelkader Krichen PII: DOI: Reference:

S0013-7944(18)30180-2 https://doi.org/10.1016/j.engfracmech.2018.07.014 EFM 6080

To appear in:

Engineering Fracture Mechanics

Received Date: Revised Date: Accepted Date:

15 February 2018 29 June 2018 6 July 2018

Please cite this article as: Soussi, H., Krichen, A., Calibration method of ductile damage model based on hybrid experimental-numerical analysis of uniaxial tensile and hole-expansion tests, Engineering Fracture Mechanics (2018), doi: https://doi.org/10.1016/j.engfracmech.2018.07.014

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Calibration method of ductile damage model based on hybrid experimental-numerical analysis of uniaxial tensile and hole-expansion tests Hela Soussi a, Abdelkader Krichen a a

Engineering Production Mechanics and Materials Unit (UGPMM), National Engineering School of Sfax, University of Sfax, B.P. 1173, 3038 Sfax, Tunisia

Abstract The present study details the calibration method of ductile damage model as a useful and sufficiently accurate procedure for modeling ductile failure in sheet metal forming. This calibration method has been performed through a hybrid experimental-numerical analysis of two characterization tests of the material properties namely the uniaxial tensile and the hole-expansion as a newly added feature of this work. Elastic-plastic finite element models with orthotropic anisotropy assumption and experiments of both characterization tests were carried out. The identification of the parameters related to the anisotropic plastic behavior, the strain-hardening after necking and the ductile damage model has been finally achieved by comparing the experimental and numerical results in terms of true tress-true strain and load-displacement relationships for the uniaxial tensile test as well as the state of the final shape resulted from the hole-expansion test. Finally, a flow chart is proposed in this paper as a practical tool which describes properly the calibration procedure. Keywords: ductile damage, calibration method, finite element model, experiment, uniaxial tensile, hole-expansion.

E-mail address: [email protected] (Hela Soussi, corresponding author) [email protected] (Abdelkader Krichen)

1

1.

Introduction

Ductile failure during plastic deformation is a phenomenon of considerable importance in many industrial applications especially, in sheet metal forming. While it may be the condition for the separation that occurs in separating processes like punching, blanking, shearing; it may limit the formability of the metal in other forming processes. In this last case, different kinds of failure regarded as defects such as orange peel aspect, necking, excessive local thinning, crack, tearing, fracture or wrinkling can be reported. From an industrial point of view, the prediction of ductile failure is among the priorities in design and optimization of metal forming processes and products. According to the specification required for the final product and the acceptation of the appearance of such failure kinds, numerous criteria have been proposed to predict and/or modeling the ductile failure. The critical thickness criterion suggested by Huang and Chen [1] has been widely used as failure criterion for predicting the forming limits due to its simplicity. This critical thickness is determined experimentally by considering that it equals to the critical thinning at fracture of a tensile specimen. Based on this experimental measurement to predict failure, most researchers have adopted such criterion for predicting the occurrence of necking [2]; [3] and tearing [1]; [4] in sheet forming processes by controlling thinning of the sheet metal during process. Recently, Soussi et al. [5] have adopted this criterion to predict the crack initiation. In their investigation, the authors have identified explicitly the critical thickness based on the inspection of the homogeneity of the minimum thickness in the necking zone of tensile specimen. They have considered that is equal to the pronounced thickness reduction leading to the crack initiation. Traditionally, the forming limit diagram (FLD) which exhibits the working range where the material can be deformed without failure has been often used as a diagnostic tool for evaluating failure in sheet metal forming. The FLD which was initially introduced by Keeler and Backofen [6] and Godwin [7] is represented in a major-minor strain space. For experimental determination of the limit strains, two main kinds of stretching have been commonly used like Nakazima out-of-plane test [8] and Marciniak in-plane test [9]. It has been established that the FLD is sensitive to the pre-formed conditions of sheet metal, in particular non-linear deformation paths in a forming process. In this context, several authors have proved that the FLD is only applicable for deformation with linear strain ratio [10]; [11]. 2

Nomenclature critical value of the damage indicator variable C1

equivalent plastic strain at the onset of damage in pure shear

C2

equivalent plastic strain at the onset of damage in uniaxial tensile

D

damage variable

Dmax

maximum degradation

e

workpiece thickness

E

Young’s modulus of undamaged material Young’s modulus of damage material

f

void volume fraction

f0

initial void volume fraction

fc

critical value of void volume fraction

f*

effective void volume fraction yield function

F, G, H, L, M, N Hill’s coefficients fracture energy fourth order Hill’s tensor j

clearance between the punch and the die second stress invariant strength coefficient

L

element characteristic length

Max

maximum equivalent plastic strain hardening exponent scalar variable of the isotropic hardening plastic anisotropy coefficients normal anisotropy coefficient

Rij

yield stress ratios deviatoric part of the Cauchy stress tensor

ui

displacement of the onset of fracture

uf

displacement of the complete fracture 3

equivalent plastic displacement equivalent plastic displacement at complete damage weighting function planar anisotropy coefficient elastic strain tensor plastic strain plastic transverse strain plastic thickness strain equivalent plastic strain equivalent plastic strain at the onset of damage equivalent plastic strain at complete damage stress triaxiality mean value of the stress triaxiality in tension ν

Poisson’s ratio Cauchy stress tensor



Cauchy stress effective stress equivalent stress

a.t

initial hole diameter indicator variable of damage

Later, the forming limit stress diagram (FLSD) formulated by the principal stresses (major stress-minor stress) has been proposed by Arrieux et al. [12] to resolve the limitation of the FLD. The concept seemed to be independent of the strain path changes and can be applied for the industrial applications showing complex deformation paths [11]. Another approach called phenomenological damage models has been considered to simulate ductile failure in sheet metal forming. Based on the stress and strain histories during forming process, such approach employs an indicator variable damage initiation when its critical value

to predict the

is reached. This variable is often taken as a

weighted cumulative plastic strain, in which the weighting function of stress state on ductile failure. 4

accounts for the effect

where

denotes the equivalent plastic strain and

represents the equivalent plastic

strain at the damage initiation. According to the form of the function

and the value of

, several phenomenological

damage models also called damage initiation criteria have been increasingly developed in the literature. By considering the function

equal to unity, it can be defined the constant

equivalent strain criterion which states that damage initiates when the equivalent plastic strain

reaches certain predefined value

. The function

can also be dependent on the

maximum stress [13], equivalent stress [14]; [15], hydrostatic pressure [16], the stresstriaxiality defined as the ratio between the hydrostatic stress and the von Mises equivalent stress (Eq. 5) [17]; [18]; [19]; [20] and the lode angle parameter associated to the third invariant of the deviatoric stress tensor [21]; [22]; [23]. These criteria involve a relatively easy calibration and implementation into commercial finite element codes, which makes their suitable for further exploitation in sheet metal forming processes [24]; [25]; [26]; [27]. All the criteria mentioned above are referred to as uncoupled approach which is carried out in post-processing step without participation in the equilibrium of structure i.e. damage does not affect the material behavior. Due to its uncoupled nature, this approach does not consider stress and stiffness softening within the material caused by damage during deformation. Under such conditions, more complex models referred to as coupled approach are currently available in the literature. In such approach, a damage variable that characterizes the damage accumulation is directly incorporated into the constitutive formulation which models softening due to damage. Historically, researches have formulated various models based on micro-mechanics damage modeling (MDM) and continuous damage mechanics (CDM). The micro-mechanics damage modeling takes into account the basic mechanism of ductile damage which is associated of the voids nucleation, growth and coalescence. In 1977, Gurson [28] developed a micro-mechanical model by considering the microscopic idealization of porous metals as aggregates containing spherical voids embedded in a metallic matrix whose behavior is governed by a rigid plastic von Mises yield criterion. This Gurson’s model can be used for metal which has an initial void volume fraction f0 as it can be 5

used for metal initially non porous (f0 = 0) [29].The Gurson’s yield plasticity function is written as dependent on the hydrostatic pressure and the void volume fraction f as an internal variable to represent damage and its softening effects on material strength. Tvergaard [30] modified the Gurson model by taking into account the voids interaction and material strain-hardening to more accurately describe the void growth. Needleman and Tvergaard [31] extended subsequently the formulation of Gurson to include the acceleration in the damage process induced by void coalescence (GTN model) through the introduction of an effective void volume fraction f*. When reaching a critical value fc of void volume fraction, the coalescence of voids begins. The Gurson model was later improved by considering the void shape effect [32], the void size effect [33], the plastic anisotropy [34] and the shear effect [35]. The micro-mechanics damage modeling has been used by several researchers to predict failure in sheet metal forming namely Brunet et al. [36]; Abbassi [37] and Kacem et al. [29]. The MDM is physically more realistic but its calibration is usually more complex. For instance, nine parameters need to be calibrated when GTN model is used for modeling ductile failure. On the other hand, the continuous damage mechanics concept has been developed within the consistent thermodynamic framework, in which the evolution of the phenomenological damage variable D is obtained through a thermodynamic dissipation potential. The CDM was first initiated by Kachanov [38] for creep damage and later was developed by Chaboche [39]; Lemaitre [40], Murakami [41] and others. The effect of damage on the material behavior was accounted for thanks to the definition of effective state variable through the assumption of the strain equivalence, stress equivalence, elastic energy equivalence or total energy equivalence. By considering for example the hypothesis of strain equivalence in the model formulated by Lemaitre [40], it is defined an effective stress

(Eq. 3) which describes the

progressive degradation within the material.

where  is the Cauchy stress and D defines the damage variable as an effective surface density of cracks or cavities in a representative volume element. The evolution law of the damage variable D begins when reaching a critical strain. This is determined experimentally through the stiffness material degradation defined by Eq.4 where E is the Young's modulus of undamaged material. 6

The Lemaitre’s model is mostly used for predicting failure in sheet metal forming as detailed in the works of Teixeira [42]; Bahloul [43] and Maria et al. [44]. However, the use of this model in shear-dominated and complex processes may lead to inaccurate predictions of maximum damage locations. To overcome this shortcoming, Cao et al. [45] have modified the Lemaitre’s model by accounting for the triaxiality and Lode parameter. From a practical point of view, some numerical codes such as Abaqus have been proposed a procedure for modeling failure in ductile metal based on the theory of coupled approach precisely, the continuous damage mechanics concept. This procedure requires the specification of two components. The first one is a damage initiation criterion for the onset of damage that can be one among the phenomenological damage models mentioned previously. Several researchers have adopted this practical procedure of coupling for modeling failure in sheet metal. As damage initiation criterion, most of them need the identification of the equivalent plastic strain at the onset of damage as a constant value [46]; [47] or a function of stress triaxiality [26]; [48]. The second component consists of using damage evolution law including a choice of element removal for modeling the progressive degradation within the material. This damage evolution law is controlled through an overall damage variable D which is activated upon the damage initiation criterion reports the onset of damage (D = 0). By introducing the concept of effective stress

defined by Eq. 3, the

stress and stiffness soften and the elements which failed are finally removed from the finite element model when they satisfy the maximum evolution law for the maximum degradation Dmax≤1. Generally, there has not been a known calibration method to identify the parameters of the damage evolution law. On this basis, various techniques available in the literature have been considered without justification in some investigations reflecting the prediction uncertainty of failure.

In the current study, the coupled approach based on the effective stress concept and offered by Abaqus code was considered along with an accurate calibration methodology of the damage model parameters was developed. This was done by performing numerical simulations and experiments for the uniaxial tensile and hole-expansion tests. A flow chart is finally proposed as a practical tool which describes the calibration procedure. The stress 7

triaxiality-dependent Bao-Wierzbicki (BW) criterion was chosen to appreciate the onset of damage. The softening behavior with stiffness degradation including the choice of element removal was considered for modeling the progressive degradation within the material after the onset of damage. An aluminum alloy 1000 series of a 0.8mm sheet thickness was selected as a study material. 2.

Damage model

2.1.

Damage initiation criterion

The stress-triaxiality dependent Bao-Wierzbicki (BW) criterion was selected as damage initiation criterion. The used phenomenological model in which the material behavior and damage are uncoupled takes into account the dependency of the equivalent plastic strain at the onset of damage

on the stress triaxiality

defined by Eq. 5. It should be noted that

may also depend on the Lode angle as mentioned in [21] and [22], but this dependency has not been considered here.

were

is the Cauchy stress tensor and

the equivalent stress.

This criterion is expressed by Eq. 6 where the indicator variable

indicates the onset of

damage in a given element when it reaches its critical value equal to the unity. It should be noted that

can also be used as a damage indicator to appreciate the damage state.

According to [22], the evolution of

as a function of

can be determined experimentally

using tensile tests with different sample shapes leading to several stress triaxiality [20]; [49]. Because of the difficulty of experiment, the evolution of this work by the relation (Eq. 7) used by [50]; [51] and [26].

8

was approximated analytically in

with

is the mean value of the stress triaxiality in tension, C1 is

C2 is

in uniaxial tensile ( = 1/3). Related to [51], C1 and C2 can be expressed by Eq. 8.

with

the hardening exponent of a power law for isotropic strain-hardening, which can be

written as

where

in pure shear ( = 0) and

is the strength coefficient.

As proposed in [20] and [26], the identification of these parameters (C1, C2,

) was

performed by means of hybrid method which combines both experiment and finite element (FE) model of the tensile test without considering damage model. 2.2.

Damage evolution law

The softening behavior with stiffness degradation based on the effective stress concept was adopted for modeling the damage evolution. Once the damage initiation has been reached (i.e.

= 1), the overall damage variable D which controls the damage evolution law is

activated and the process of softening behavior begins. Because of the strain localization during the damage evolution, the softening response is dependent on the mesh. To alleviate the mesh dependency, Abaqus code adopted a fracture energy method originally proposed by Hillerborg et al. [52] for modeling the quasi-brittle fracture of concrete. With this approach, the softening response after damage initiation is characterized by introducing stress–displacement relation for softening rather than stress–strain relation. The implementation of this stress-displacement concept in a finite element model introduces the definition of an element characteristic length L which depends on the element geometry and formulation. It is defined for example as the square root of the integration point area for 2D elements and the cubic root of the integration point volume for 3D elements [53]. The fracture energy is then given by Eq. 9:

9

This expression introduces the definition of the equivalent plastic displacement conjugate of the yield stress

as the

after the onset of damage. Before the damage initiation,

is 0 that evolves according to

as soon as the damage begins.

In the code Abaqus, the evolution of the damage variable D with the equivalent plastic displacement can be specified in linear, bilinear/piecewise or exponential form. For linear evolution form which was chosen in this study due to its simplicity, the damage variable D increases according to:

This definition ensures that when the equivalent plastic displacement

reaches the value

(Eq. 11), the stress and the stiffness of the material will be fully degraded and the elements which fail will be removed from the model. The linear evolution defines a truly linear stress-strain softening response only if the effective response of the material is perfectly plastic after damage initiation.

In order to account for progressive degradation, the equivalent plastic displacement failure should be identified. By assuming a lower value of

at

equal to 0.1mm, Ribeiro et al.

[48] have adopted a linear damage evolution law in order to predict the fracture behavior for notched and different unnotched specimens subjected to uniaxial tensile test. In the investigation of Thairy and Wang [47], the displacement

was calculated with respect to

Eq. 11. The element characteristic length L was determined based on the mesh size and the value of

was taken from uniaxial tensile test at complete damage. They have considered a

linear evolution law to predict the failure modes of axially compressed steel columns subjected to transverse impact. However, Hogstrom et al. [46] have shown that the equivalent plastic strain at complete damage

taken from uniaxial test is strongly

dependent to the length scale of virtual extensometer over which the strain evaluation has been performed. To take into account this dependence, the authors have considered

as a

function dependent on the length scale by using an expression proposed by Yamada et al. [54] known as Barba’s relation. By comparing between the linear and bilinear damage 10

evolution law, they showed that the bilinear relationship resulted in excellent agreement with the experimental and the numerical true stress–true strain curves. Instead of using analytical methodology, Kacem et al. [26] have identified the displacement

and the

maximum degradation Dmax using inverse method by means of numerical simulations of tensile test. Because of the excessive distortion of the elements generated with the decrease of its stiffness, the authors have used a value of Dmax lower than unity compared to the above-mentioned references which considered the value equal to unity. In the investigation of Kacem et al. [26], the values of

and Dmax leading to the first element deletion

corresponding to the displacement at complete fracture of the tensile specimen taken from experiment, were kept. Accordingly, different ways for identifying

has been considered in the literature due to

the difficulty of its determination directly from an experimental measurement. Despite that it is well known that the fracture point in the uniaxial tensile test seems to be an inaccurate measure, several researches have adopted this fracture point as a measurement tool to determine the value of

. While the authors have noted that the considered value of

and Dmax resulted a good approximation of failure, this does not hide the presence of others values of both parameters leading to an accurate prediction of failure, especially when they have been used without any justification in the literature. In the current study, the choice of the values of

and Dmax was justified related to the

proposed calibration method. Initial simulations were performed with larger values of Dmax close to unity. The results demonstrated that the stiffness of the element decreases leading to its excessive distortion. To solve this problem, the element was not deleted for Dmax close to unity but rather for a critical value Dmax < 1, such as in [26]. Therefore, the description of the total degradation of the material required the identification of both parameters

and

Dmax. This identification was initially implemented by an inverse method using the numerical simulation of the tensile test with damage model. The values of

and Dmax leading to a

more realistic description of the complete fracture as well as a good correlation of the experimental load-displacement curve, are kept. The obtained parameters were subsequently tested using the hole-expansion test for selecting appropriate values of

11

and

Dmax. This newly added feature will be justified later in next section by discussing the limitation observed by using only the uniaxial tensile test. 3.

Material characterization

The 1050-H14 aluminum alloy supplied in 0.8mm thin sheet was chosen as a study material. Its behavior was characterized by performing experiments of the uniaxial tensile and the hole-expansion tests.

4.1.

Tensile test experiments

The experiments of tensile test were carried out on rectangular specimens of gauge dimensions 200mm×20mm using a universal tension-compression machine of maximum load 10kN at a constant tensile speed of 5mm/min. Specimens were initially shear cut along the rolling direction (0°/RD) and then machined to eliminate the hardened area due to cutting. The edges were then locally polished with sandpaper in the center of the sample to reduce the width by some tenths of a millimeter in order to obtain necking in this zone. After clamping, the useful length of the specimen was fixed to be 150mm. The engineering strain was measured by means of a virtual extensometer defined on the specimen, which gauge length L0 is 50mm. The longitudinal true strain is given by:

The true stress also called Cauchy stress

was calculated by the ratio between the load F

and the current section S0 of the specimen before necking. The current section was determined from the initial section of the specimen by assuming an isochoric plastic deformation and a uniform distribution of strain over the gauge length of the specimen.

To evaluate the anisotropic plastic behavior of the material, tensile tests were also performed at 45° and 90° to the rolling direction. Five tests per each configuration were performed in order to ensure the reproducibility of results.

12

The obtained true stress-true strain curves are shown in Fig. 1. Note that, the portion of curve after necking which is initiated when the tensile force reaches its maximum is presented in indicative way. 160 140

Onset of necking

Cauchy stress (MPa)

120

100 80 60

0°/RD 45°/RD 90°/RD

40

20 Uniform strain

0 0

0.01

Localized strain 0.02 0.03 Longitudinal true strain

0.04

Fig. 1 Experimental true stress-true strain curves obtained at 0°, 45° and 90° to the RD

The measurement of a linear interpolation of the evolution of the plastic transverse strain versus the plastic thickness strain (

leads to the plastic anisotropy coefficients

) given in table 1. The normal anisotropy coefficient

significant but the planar anisotropy coefficient

is is rather weak.

The Young’s modulus E of about 68GPa was determined from the experimental curves. The Poisson’s ratio ν was taken equal to 0.28.

Table 1 Anisotropic coefficients of the 1050A-H14

0.37

0.6

0.98

0.64

0.075

The load-displacement curve in the rolling direction is depicted in Fig. 2. The displacement ui corresponding to the onset of fracture can be determined as mentioned in [55] and [26] using either two-or three-dimensional digital image correlation (DIC). Based on the DIC measurements, the instant of onset of fracture (not the location) was defined by the first detectable discontinuity in the measured displacement field of the tensile test. In this work, 13

the displacement ui was identified referred to the work of Soussi et al. [5] using cross section views of necking zone at different positions for three levels of load drop after reaching the maximum load. Based on the inspection of the homogeneity of the minimum thickness in the necking zone, the fracture initiation was identified at the load drop between 15% and 20% of the maximum load. The displacement ui was then considered as the average of both displacements limits at 15% and 20% and it was evaluated equal to 1.95mm as shown in Fig. 2. For the displacement uf corresponding to the complete fracture, it was considered near to the value 2.11mm which corresponds to the point at the end of the experimental loaddisplacement curve (Fig. 2).

2500 2.5

Onset of necking

2 2000

Load (kN)

1.5 1500 1 1000

0.5 500

uf  2.11mm

ui = 1.95mm 00 0

0.4

0.8 1.2 Displacement (mm)

1.6

2

Fig. 2 Experimental load-displacement curve in 0°/RD

4.2.

Experiments of hole-expansion test

The hole-expansion test consists of expanding a hole in a sheet metal with a conical, cylindrical or hemispherical punch into a die and leads to form a flange around the periphery of the hole as shown in Fig. 3. This test is performed by stretching that induces tensile stress component in the circumferential direction at the edge of the hole and caused thinning of

14

the sheet. According to the magnitude of expansion, several defects can be occurred namely necking, crack and tearing [5]; [26].

R6 R4

Displacement

Conical Punch

Blank-holder

Workpiece R0.2

0.8

Blank holding

Flange

R2

Die

R6.8 R15

0.8 30

Fig. 3 Tool geometry parameters of the conventional hole-expansion test In this study, the hole-expansion test was carried out with a dedicated conventional holeexpansion tool including a conical punch, a die and a blank-holder (Fig. 4-a). The punch diameter was 12mm and the half-angle at the top of the cone was 30°. The die diameter of 13.6mm was used to set a clearance j between the punch and the die equal to the workpiece thickness e (j = 0.8mm). The punch profile radius and the die profile radius were taken equal to 4mm and 0.2mm, respectively.

15

(a)

(b) Conical punch

Workpiece Die

Blank-holder

Fig. 4 (a) Experimental tool and (b) finite element model of hole-expansion test (for sake of clarity, the blank-holder is not represented)

The workpiece was mechanically cut from the 1050-H14 aluminum alloy sheet. It was a washer with an external diameter of 30mm and a drilled hole of 4mm diameter in the center. The initial hole diameter a.t was selected in order to generate an excessive expansion in which tearing tends to occur at the flange extremity (Fig. 5) as already detailed in previous work of Soussi et al. [26]. The geometrical characteristics of the tool and the workpiece are illustrated in Fig. 3. To prevent the upward motion of the workpiece during the test, the workpiece was clamped between the die and the blank-holder without applying a large clamping force as shown in the work of Krichen et al. [56]. The experiments were carried out without lubricant using a universal tension-compression testing machine of maximum load 10kN. The punch displacement was kept to be at a constant low rate of 5mm/min.

16

0°

90° Fig. 5 Experimental result of hole-expansion test with a.t = 4mm and j = 0.8mm showing pronounced tearing at the flange extremity

4.

Numerical modeling

5.1.

Uniaxial tensile test

Although the tensile test satisfies the plane stress conditions in the uniform strain range, the localization of plastic strain from necking initiation disrupts the uniformity of the stress field in the gauge length of the specimen. Therefore, because of the development of stress in the thickness direction, the stress field becomes a three-dimensional nature. As for dealing with necking problems, a 3D finite element (FE) model was built. The numerical simulations of the tensile test were conducted with the explicit version of the Abaqus code (v6.10) which uses a dynamic explicit resolution scheme. Similarly to the experimental tensile test, a small notch was executed in the mesh as schematically shown in Fig. 6, the size of which was adjusted for obtaining the onset of the complete fracture located on the edge and in the middle of the specimen. Since Dunand and Mohr [55] have shown that the shell elements are inappropriate for post-necking analysis of local strain fields, the eight-node linear brick reduced-integration solid elements (C3D8R) were used to mesh the specimen. The choice of a reduced integration element was made to speed up the computation as this element uses only one integration point per element instead of four. An automatic mesh generator was used to generate the FE mesh. To determine precisely the equivalent plastic strain at onset of damage

, a finer element size of about 0,15mm3 was

chosen in the useful area of the specimen. The total number of elements and nodes was 17

198090 and 240234, respectively with 5 elements through the sheet thickness. The nodes located on the left side are clamped while those located on the right side are subjected to tensile loading. To account for the anisotropic behavior of the material, a coordinate system (x, y, z) located into the specimen was then added. 150mm

50mm

20mm

Gauge length of the virtual extensometer

y x Local coordinate system of anisotropy 2 1 Global coordinate system

Fig. 6 Finite element model of the tensile test 5.2.

Hole-expansion test

Three-dimensional FE simulations of the hole-expansion test were established by taking into account the anisotropic behavior of the material. The tools (punch, die and blank-holder) were defined as rigid surfaces. The workpiece was modeled as a deformable body. The contact between the tools and the workpiece was assumed with friction governed by the coulomb’s law with constant friction coefficient of 0.2. All simulations were performed by considering the entire geometry of the workpiece. This last was discretized with eight-node linear brick reduced integration elements (C3D8R). To minimize the element size dependency, the adaptive meshing or remeshing techniques were not used but rather the elements in the contact area were refined. In addition, the mesh size (0.15mm3) was chosen similar to the one used in the tensile test simulations (i.e. in the identification of

). The total number of elements in the workpiece was about 124110

elements. The FE model is presented in Fig. 4-b. During the simulation, the die and the blank-holder were fixed and the punch was constrained to move in the axial direction.

18

5.3.

Precautions

It is important to note that the use of explicit version which is mainly dedicated to the simulations of dynamic problems requires some precautions when it is used to solve quasistatic problems. To reduce the simulation time and simultaneously satisfy the quasi-static condition, an appropriate mass scaling factor was used to avoid inertial effects. On this basis, the stability of the numerical solutions obtained in the tensile test or the hole-expansion test were checked by analyzing the ratio between the kinetic energy and the internal energy. Therefore, the mass scaling factor was adjusted to ensure the values of this ratio between 5% and 10% (Fig. 7). In addition, the use of reduced integration elements which are controlled by the “hourglassing control” option requires checking the ratio between the internal energy and the artificial deformation energy. The maximum value of this ratio, which did not exceed 4% in this work, proves that the reduced integration elements did not affect the accuracy of the results.

Ratio between energies

0.1 Energie cinétique / Energie interne Kinetic energy / Internal energy

0.08

Artificial artificielle deformationdeenergy / Internal/energy Energie déformation Energie interne

0.06

0.04

0.02

0 0

1

2

3

4

Displacement (mm)

Fig. 7 Typical evolutions of energies ratio 5.

Material behavior modeling

The material was assumed to have an elastic-plastic behavior with an isotropic hardening. The elastic behavior was assumed to be linear and isotropic with Young’s modulus E and Poisson’s ratio ν. It was described by the Hook’s law in which the Cauchy stress tensor the elastic strain tensor

are related by: 19

and

where

is the identity tensor, the material parameters μ et λ are the classical Lame’s

constants for isotropic elasticity. The material was also assumed to be orthotropic and its plastic behavior was modeled using a standard plasticity theory with an anisotropic quadratic Hill’s 1948 [57] yield criterion and an isotropic hardening model.

6.1.

Hill’s 1948 yield criterion

The Hill’s 1948 criterion was adopted because of its simplicity and its availability in most of the codes of finite element calculations in particular Abaqus. Formally, the yield function is given by:

denotes the second stress invariant that is defined in the Hill’s 1948 yield criterion as follow:

where

is the deviatoric part of the Cauchy stress tensor

and

is the fourth order Hill’s

tensor defined by the six independent coefficients: F, G, H, L, M and N. When the uniaxial tensile along the rolling direction (RD) was considered as a reference state for the strain-hardening and the initial yield stress, the relation G + H = 1 is systematically imposed. In this condition, the Hill’s coefficients are related to the anisotropic coefficients and

,

by the following relations (Eq. 17):

In this paper, the Hill’s coefficients (F, G, H, and N) were determined in the uniform strain range by directly fitting (inverse method) in order to minimize the gap between experimental and simulated stress-strain curves in the three directions (0°, 45° and 90°). This optimization was practiced under the obligation to keep almost constant the value of normal 20

anisotropy coefficient

and by using particularly a volume element (VE) subjected to tensile

loading. It is important to note that the Hill’s coefficients are integrated into the Abaqus code through yield stress ratios Rij called and defined by Eq. 18.

For reasons of simplicity and due to the lack of available data regarding the mechanical behavior in the sheet thickness, the coefficients L and M were kept equal to their isotropic values (i.e., L = M = 1.5). The best fit of Hill’s coefficients and the yield stress ration Rij which lead to the good description of the anisotropic behavior of the 1050A-H14 are given in table 2. The corresponding anisotropic coefficients were calculated from Eq. 17 and they are also reported in the same table. The numerical true stress-true strain curves in the uniform strain range showing a good correlation with experimental results are presented in Fig. 8.

Table 2 Obtained anisotropic coefficients, Hill’s coefficients and yield stress ratios Anisotropic coefficients 0.615

0.64

0.74

Hill’s coefficients

Yield stress ratios

F

G

H

L

M

N

R11

R22

R33

R12

R13

R23

0.515

0.62

0.38

1.5

1.5

1.29

1

1.058

0.94

1.079

1

1

21

160

Cauchy stress (MPa)

140

120 100 Experiment0°/ 0°/RD Expérience DL Experiment45°/ 45°/RD Expérience DL Experiment90°/ 90°/RD Expérience DL

80 60

Model 0°/RD Modèle 0°/ DL Model 45°/RD Modèle 45°/ DL Model 90°/RD Modèle 90°/ DL

40 20 0

0

0.005

0.01

0.015

Longitudinal true strain

Fig. 8 Experimental and numerical true stress-true strain curves in the uniform strain range at 0°, 45° and 90° to the RD

6.2.

Strain-hardening

For the isotropic hardening function described with the scalar variable R (Eq. 15), its identification was performed by distinguishing between the uniform and localized strain ranges which are limited by the occurrence of necking (Fig. 1). In order to establish the evolution of the strain-hardening before necking initiation, the true stress-true plastic strain data measured along the rolling direction in the uniform strain range was used. After necking, this evolution was estimated by fitting a piecewise linear extrapolation. This approach was chosen instead of the classical approach of saturation hardening which consists in extrapolating the hardening curve from the stress value obtained at necking. This choice was adopted because the classical approach leads generally to a poor modeling by presenting a premature drop-off of the load. The piecewise linear extrapolation was identified by performing an adjustment based on a hybrid method involving the experiment and the numerical simulation of tensile test as mentioned in [58] and [26]. This adjustment was practiced by reducing the gap between the experimental and numerical load-displacement curves with respecting initially the tangency condition immediately after necking.

22

The evolution of the strain-hardening after necking estimated by fitting a piecewise linear extrapolation is shown in Fig. 9. The hardening curve obtained from the classical approach of extrapolation is also superimposed in the same figure. 180

necking

160

Cauchy stress (MPa)

140 120 Experiment0°/ DL Expérience

100

Saturation hardening

80

Linear piecewise hardening Ecrouissage extrapolé par morceaux

60 40 20 0 0

0.1

0.2 True plastic strain

0.3

0.4

Fig. 9 Evolution of strain-hardening before and after necking Fig. 10 presents the numerical load-displacement curve of tensile test obtained by applying the piecewise linear extrapolation from necking up to onset of fracture (i.e. ui) showing a good agreement with the experimental curve. The numerical load-displacement curve obtained from saturation hardening is also superimposed in the same figure. This confirmed well that the classical approach cannot describe suitably the load-displacement relationship in the localized strain range. To verify the element size dependency on the strain-hardening evolution, the loaddisplacement curves obtained with a finer mesh size of 0.1mm 3 and a coarse mesh size of 0.3mm3 are superimposed in Fig. 10. The comparison between the obtained results shows that a few dependency of the element size is observed at the end of the curve. This ascertainment is in agreement with the investigation of Dunand and Mohr [55] in which the authors have shown that the element size does not have a large influence on the loaddisplacement curve.

23

2500 2.5

2000 2

Load (kN)

1500 1.5

1000 1

Expérience 0°/ DL Experiment Ecrouissage Saturationsaturant hardening Ecrouissage extrapolé hardening par morceaux(mesh (taille de maillage 0,3mm) Linear piecewise size 0.3mm)

500 0.5

Ecrouissage extrapolé hardening par morceaux(mesh (taille de maillage 0,1mm) Linear piecewise size 0.1mm) Ecrouissage extrapolé hardening par morceaux(mesh (taille de maillage 0,15mm) Linear piecewise size 0.15mm)

00

0

0.5

1 Displacement (mm)

1.5

2

Fig. 10 Experimental and numerical load-displacement curves

6.

Identification of the BW criterion

The equivalent plastic strain at onset of damage C2 in uniaxial tensile was identified according to the methodology used by [55] and [26]. This methodology was based on the numerical simulation of the tensile test without damage and the use of the displacement ui corresponding to the onset of fracture (Fig. 1). The numerical simulation of the tensile test with the best fit of the piecewise linear extrapolation from necking identified in the last section was considered. Post-processing of this simulation gave then access to the evolution of the equivalent plastic strain and the stress triaxiality . Referring to [55], it was assumed that the location of the onset of damage coincides with the location of the maximum equivalent plastic strain (Max

) within the specimen at the instant of onset of fracture. The

corresponding equivalent plastic strain is referred to as the strain of damage initiation

in

tensile test. Fig. 11 displays the numerical repartition of the equivalent plastic strain in the tensile specimen for several displacements (0.5mm, 1mm and 1.95mm = ui). The analysis of the different displacements shows that the maximum of the equivalent plastic strain is located at the geometrical default. The maximum strain C2 is then evaluated to 1.086 for the displacement ui. By using the Eq. 8 and for the hardening exponent

of 0.14, the value of C1

was then calculated to 0.39. The strain C1 can also be determined by performing experiment of pure shear test. 24

The mean stress triaxiality

was evaluated referring to Kacem et al. (2013) [26] from the

evolution of the triaxiality

of the element presenting the maximum equivalent plastic

strain during the tensile test (Fig. 12). These values (C1, C2, as a function of the stress triaxiality

Displacement 0.5mm

0.052 0.048 0.044 0.039 0.035 0.031 0.026 0.022 0.018 0.013 0.009 0.005 0.000

Max

= 0.052

0.126 0.115 0.105 0.094 0.084 0.073 0.063 0.053 0.042 0.032 0.021 0.011 0.000

) lead finally to the evolution of

presented in Fig. 13.

Displacement 1mm

Max

1.086 0.995 0.905 0.814 0.724 0.633 0.543 0.453 0.362 0.272 0.181 0.091 0.000

= 0.125

Displacement 1.95mm = ui

Max

= 1.086

Fig. 11 Repartition of the equivalent plastic strain in the tensile specimen for several displacements

0.8

Stress triaxiality

0.6

η0 = 0.44mm 0.4

0.2

0 0

0.5

1 Displacement (mm)

1.5

2

Fig. 12 Evolution of the stress triaxiality versus the displacement

25

3

𝑝 0

η0 = 0.44mm

2

1

-0.4

-0.2

0

0.2 0.4 Taux de triaxiliaté ( )

Fig. 13 Evolution of

0.6

0.8

as a function of stress triaxiality

7.

Calibration of the softening behavior with stiffness degradation

8.1.

Limitation of the calibration method using only tensile test

The identification of

1

and Dmax was done with the procedure defined previously in section

(2.2). Firstly, the ductile damage model was implemented into the finite element model of the tensile test by introducing the tabular values of as initial values of

and the stress triaxiality (Eq. 7) as well

and Dmax. For this first simulation, the values of

and Dmax equal to

0.5 and 0.55, respectively were given from [26] in which they were identified for the same material of 2mm sheet metal thickness. The load-displacement curve obtained from the numerical simulation is plotted in Fig. 14-a. Different shapes change of the specimen in the necking zone representing the repartition of the equivalent plastic strain from the onset of necking until the complete fracture are also exhibited in Fig. 14-b. They are identified by the values of the corresponding displacement which are marked by six points in blue color in the numerical load-displacement curve (Fig. 14-a). The analysis of the necking zone reveals clearly that the complete fracture (element deletion) of the specimen is generated after an excessive distortion of elements. By considering the experimental displacement of the complete fracture uf of the specimen near to 2.11mm (Fig. 1), it can be noted that this distortion retards relatively the complete

26

fracture for the predicted displacement uf of 5.65mm. This poor description of the complete fracture is proved by comparing both experimental and numerical load-displacement curves. On this basis, several simulations were performed with others values of

and Dmax in a

range which reduces the excessive distortion of elements before removing from model. A typical load-displacement curve for the reduced values (

= 0.2 and Dmax = 0.15) is

superimposed in Fig. 14-a. In addition, various shapes change of the specimen in the necking from the onset of necking until the complete fracture, are also presented in the Fig. 14-c. They are identified by the values of the corresponding displacement which are marked by five points in red color in the numerical curve. The analysis of the obtained results for this specific set of parameters shows that the excessive distortion of elements is significantly reduced. The comparison between the experimental and the numerical curves shows that the considered values of

= 0.2 and Dmax = 0.15 lead to a better description of the

experimental shape. The complete fracture of the specimen is predicted for the displacement uf equal to 2.58mm. As it can be seen in Fig. 15, this result is also obtained for various sets of parameters (

and Dmax) which conclude practically a good correlation of the

experimental load-displacement curve. At this stage, it can be concluded that the calibration methodology based only on the tensile test leads to the multiplicity of the parameters set (

and Dmax) which controls the total

degradation of the material. To consider the above-mentioned limitation, there is the need to the hole-expansion test to be used as a practical method for selecting the appropriate parameters set ( The final calibration using the hole-expansion test is detailed in the next section.

27

and Dmax).

(b) 1

2 0.899 0.824 0.749 0.672 0.599 0.524 0.450 0.375 0.300 0.225 0.150 0.075 0.000 Displacement 1.82mm

0.086 0.078 0.071 0.064 0.057 0.050 0.043 0.036 0.029 0.022 0.015 0.007 0.000 Displacement 0.76mm

(c) 1

0.086 0.078 0.071 0.064 0.057 0.050 0.043 0.036 0.029 0.022 0.015 0.007 Displacement 0.76mm 0.000 0.902 0.826 0.751 0.676 0.601 0.526 0.451 0.376 0.301 0.226 0.151 0.076 0.000

2

Displacement 1.82mm

4

3

2.665 2.443 2.221 1.999 1.777 1.555 1.333 1.111 0.889 0.667 0.445 0.223 Displacement 3.48mm 0.000

Displacement 2.11mm

6 3.349 3.070 2.791 2.512 2.233 1.954 1.675 1.396 1.117 0.838 0.559 0.280 Displacement 5.65mm 0.000

(a) Onset of necking 1 Onset of damage (w = 1) 2

2500 2.5

2 2000

1

3

2

1.5 1500

Experiment 3 Expérience 0.5;;vc Dmax = 0.55 up == 0,5 = 0,56

Load (kN)

5 3.214 2.946 2.678 2.410 2.142 1.875 1.607 1.339 1.071 0.804 0.536 0.268 0.000 Displacement 4.35mm

1.261 1.156 1.051 0.946 0.841 0.736 0.631 0.526 0.421 0.316 0.211 0.106 0.000

1000 1

0.2;;vc Dmax = 0.15 up == 0,2 = 0,15

4

4

500 0.5

5

Complete fracture of 5 specimen

00 0

1

2

3

1.086 0.996 0.905 0.815 0.724 0.634 0.543 0.453 0.362 0.272 0.181 0.091 Displacement 2.11mm 0.000

3 4 Displacement (mm)

4

1.171 1.073 0.976 0.878 0.781 0.683 0.586 0.488 0.391 0.293 0.196 0.098 Displacement 2.35mm 0.000

Complete fracture of specimen 6

5

6

5

1.192 1.093 0.993 0.894 0.795 0.696 0.596 0.497 0.398 0.298 0.199 0.100 Displacement 2.58mm 0.000

Fig. 14 (a) Experimental and numerical load-displacement curves, Shapes change of the specimen showing the repartition of the equivalent plastic strain from the onset of necking until complete fracture for (b) = 0.5, Dmax = 0.55 and (c) = 0.2, Dmax = 0.15

28

2500 2.5

1500 1.5

Expérience 0°/ DL = 0.55 up = 0.5 0,5 ; D vcmax = 0,56 = 0.15 up = 0.35 0,35 ; D vcmax = 0,15 = 0.2 up = 0.21 0,21 ; D vcmax = 0,2 20002 up = 0.2 0,2 ; D vcmax = 0,15 = 0.15

1 1000

1500 1.5

Effort (kN)

Effort (kN)

2000 2

0.5 500

10001 1.5

2 Déplacement (mm)

00 0

1

2

3 4 Déplacement (mm)

5

6

2.5

7

Fig. 15 (a) Experimental and numerical load-displacement curves for different values of and Dmax

8.2.

Final calibration using hole-expansion test

Each parameters set (

and Dmax) able to describe correctly the tensile test in terms of load-

displacement curve and complete fracture was tested by using the hole-expansion test. The parameters set which finally selected was the one whose provide qualitatively a predictive capability of the experiment in terms of flange state. Fig. 16 presents typical final shapes of the flange showing the repartition of the variable D. Contrary to that observed previously for the tensile test, the state of the flange displays a significant sensitivity to the variation of

and Dmax values. By reducing these values, more

pronounced tearing is predicted in the flange extremity. The qualitative comparison of the predicted final shapes with the experimental result (Fig. 5) leads for selecting the values of and Dmax equal to 0.21 and 0.2 respectively, which are able to describe correctly the process of ductile damage of the material subjected to tensile and expansion loading. For more comparison, the numerical simulation using

= 0.5 and Dmax = 0.55 which

predicts a final shape free from tearing and does not describe the experimental result is also presented in Fig. 16.

29

= 0.35 ; Dmax = 0.15

= 0.5 ; Dmax = 0.55 D

0.55 0.51 0.46 0.42 0.37 0.32 0.28 0.23 0.18 0.14 0.09 0.05 0.00

Dmax

D

0.15 0.14 0.13 0.11 0.10 0.09 0.08 0.06 0.05 0.04 0.03 0.01 0.00

Dmax

= 0.21 ; Dmax = 0.2

D

0.20 0.18 0.17 0.15 0.13 0.12 0.10 0.08 0.07 0.05 0.03 0.02 0.00

= 0.2 ; Dmax = 0.15

Dmax

D

0.15 0.14 0.13 0.11 0.10 0.09 0.08 0.06 0.05 0.04 0.03 0.01 0.00

Dmax

Fig. 16 Typical predicted final shapes of the flange for different values of

and Dmax

Finally, in order to achieve the accurate calibration of the ductile damage parameters, there are different steps to be followed. The calibration procedure of the damage model is schematically described by a flow chart proposed in Fig. 17.

30

(1)

Tensile test experiments True stress-True strain curves vs.

Anisotropic coefficients

curves

Initial yield stress ratio Rij

RD = Reference

Numerical modeling of VE subjected to tensile loading

Similarity of true s-s curves up to necking

(2)

Adjustment

Strain-hardening before necking

N

Y Yield stress ratio Rij

Hill’s coefficients

Numerical modeling of tensile test (3D) without damage

Initial strain-hardening after necking Adjustment

(3)

Similarity of L-D in RD up to onset of fracture (ui)

N

Pure shear experiment Hardening exponent

Y Stop at ui

Strain-hardening after necking

Maximal equivalent plastic strain (Max )

C2 2

Displacement ui at onset of fracture Displacement uf at complete fracture

Load-Displacement curves

Locating the element of Max

1 C1

Stress triaxiality vs. displacement BW criterion (

v

Initial values of (

Adjustment

) Numerical modeling of tensile test (3D) with damage

, Dmax)

(4)

Similarity of L-D up to complete fracture

N

Y Values of (

, Dmax)

(6) Numerical modeling of hole-expansion test (3D)

Adjustment

N

Qualitative similarity of the flange state

Boundary conditions (5) Experiment of hole-expansion test in severe condition Flange state

Y Final calibration

Fig. 17 Flow chart of calibration methodology 31

9.

Conclusion

This study describes the proposed methodology for calibrating the ductile damage model proposed by Abaqus which is introduced the effective stress concept for coupling damage and material behavior. This damage model involves the stress-triaxiality Bao-Wierzbicki criterion to appreciate the onset of damage and the softening behavior with stiffness degradation including the choice of element removal for modeling the progressive damage within the material. The calibration methodology enriches the hybrid numericalexperimental analysis of both characterization tests such as the uniaxial tensile and the holeexpansion to identify all the parameters. -The Hill’s coefficients in the uniform strain range was identified in such a way to obtain a good description of the anisotropy behavior in terms of true stress-true strain relationship. This was done through an inverse method using the volume element subjected to tensile loading. - The strain-hardening after necking (localized strain range) was predicted by considering a piecewise linear extrapolation approach which was identified by performing an adjustment of the numerical load-displacement curve resulted from numerical simulation of tensile test compared to the experimental curve. - The identification of Bao-Wierzbicki parameters (C1, C2,

) was performed by means of

hybrid method which uses experimental results and numerical simulation of the tensile test without damage which gave access to the evolution of the equivalent plastic strain and the stress triaxiality. - The identification of the parameters (

and Dmax) which control the total degradation of

the material was carried out through two steps. In the first step, these parameters were initially calibrated by inverse method in order to obtain a good correlation between experiment and numerical results in terms of load-displacement relationship. Under these conditions, it was concluded various parameters sets leading to a good description of the complete fracture of tensile specimen. Because of the multiplicity of the calibration solution resulted from tensile test, these identified parameters are subsequently validated in second step by simulating the hole-expansion test in its severe condition. This was done by 32

comparing qualitatively the experiment and the numerical simulation in terms of the flange state. The used test was finally proven efficient for choosing an adequate solution of the parameters set (

and Dmax).

- The identified damage model can be incorporated into finite element model in order to be used for predicting failure in metal forming process and this will be the future work.

Acknowledgments This work is carried out thanks to the support and funding allocated to the Engineering Production Mechanics and Materials Unit (UGPMM/UR17ES43) by the Tunisian Ministry of Higher Education and Scientific Research.

References [1] Huang YM, Chen JW. Influence of the tool clearance in the cylindrical cup-drawing process. J Mater Process Technol 1996;57:4–13. [2] Chen TC. An analysis of forming limit in the elliptic hole-flanging process of sheet metal. J Mater Process Technol 2007;192–193:373–80. [3] Leu DK, Chen TC, Huang YM. Influence of punch shapes on the collar-drawing process of sheet steel. J Mater Process Technol 1999;88:134–46. [4] Huang YM, Chien KH. Influence of Cone Semi-Angle on the Formability Limitation of the Hole-Flanging Process. Int J Adv Manuf Technol 2002;19:597–606. [5] Soussi H, Masmoudi N, Krichen A. Analysis of geometrical parameters and occurrence of defects in the hole-flanging process on thin sheet metal. J Mater Process Technol 2016;234:228–42. [6] Keeler S, Backofen W. Plastic instability and fracture in sheets stretched over rigid punches. Trans ASM 1963;56:25–48. [7] Goodwin G. Application of strain analysis to sheet metal forming problems in the press shop. SAE Tech Pap 680093 1968. [8] Nakazima K, Kikuma T, Hasuka K. Study on the formability of steel sheets. Yawata Tech Rep 1968;264:8517–30. [9] Marciniak Z, Kuczyński K. Limit strains in the processes of stretch-forming sheet metal. Int J Mech Sci 1967;9:609–20. [10] Butuc MC, Gracio JJ, Barata da Rocha A. An experimental and theoretical analysis on the application of stress-based forming limit criterion. Int J Mech Sci 2006;48:414–29. [11] Uthaisangsuk V, Prahl U, Bleck W. Stress based failure criterion for formability characterisation of metastable steels. Comput Mater Sci 2007;39:43–8. 33

[12] Arrieux C, Bedrin C, Boivin M. Determination of an Intrinsic Forming Limit Stress Diagram for Isotropic Metal Sheets. Proc 12th Bienn Congr IDDRG 1982:61–71. [13] Cockroft M, Latham D. Ductility and the Workability of Metals. J Inst Met 1968;96:33–9. [14] Clift SE, Hartley P, Sturgess CEN, Rowe GW. Fracture prediction in plastic deformation processes. Int J Mech Sci 1990;32:1–17. [15] Oh S, Chen C, Kobayashi S. Ductile fracture in axisymmetric extrusion and drawing: Part 2, workability in extrusion and drawing. J Eng Ind – Trans ASME 1979;101:36–44. [16] Brozzo P, Deluca B, Rendina R. A new method for the prediction of formability limits in metal sheets, Sheet metal forming and formability. Proc Seventh Bienn Conf Int Deep Draw Res Group n.d. [17] Oyane M, Sato T, Okimoto K, Shima S. Criteria for ductile fracture and their applications. J Mech Work Technol 1980;4:65–81. [18] Ayada M, Higashino T, Mori K. Central bursting in extrusion of inhomogeneous materials. Proc First ICTP Adv Technol Plast 1984;1:553–558. [19] Johnson GR, Cook WH. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures. Eng Fract Mech 1985;21:31–48. [20] Bao Y, Wierzbicki T. On fracture locus in the equivalent strain and stress triaxiality space. Int J Mech Sci 2004;46:81–98. [21] Xue L. Damage accumulation and fracture initiation in uncracked ductile solids subject to triaxial loading. Int J Solids Struct 2007;44:5163–81. [22] Bai Y, Wierzbicki T. A new model of metal plasticity and fracture with pressure and Lode dependence. Int J Plast 2008;24:1071–96. [23] Bai Y, Wierzbicki T. Application of extended Mohr–Coulomb criterion to ductile fracture. Int J Fract 2010;161:1. [24] Takuda H, Mori K, Hatta N. The application of some criteria for ductile fracture to the prediction of the forming limit of sheet metals. J Mater Process Technol 1999;95:116– 21. [25] Thipprakmas S, Jin M, Murakawa M. Study on flanged shapes in fineblanked-hole flanging process (FB-hole flanging process) using finite element method (FEM). J Mater Process Technol 2007;192–193:128–33. [26] Kacem A, Krichen A, Manach PY, Thuillier S, Yoon JW. Failure prediction in the holeflanging process of aluminium alloys. Eng Fract Mech 2013;99:251–65. [27] Dewang Y, Hora MS, Panthi SK. Prediction of crack location and propagation in stretch flanging process of aluminum alloy AA-5052 sheet using FEM simulation. Trans Nonferrous Met Soc China 2015;25:2308–20. [28] Gurson AL. 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 1977;99:2– 15. [29] Krichen A, Kacem A, Manach P. Prediction of damage in the hole-flanging process using a physically based approach. Int J Damage Mech 2014;24:840–858. 34

[30] Tvergaard V. Influence of voids on shear band instabilities under plane strain conditions. Int J Fract 1981;17:389–407. [31] Tvergaard V, Needleman A. Analysis of the cup-cone fracture in a round tensile bar. Acta Metall 1984;32:157–69. [32] Gologanu M, Leblond JB, Devaux J. Approximate Models for Ductile Metals Containing Nonspherical Voids—Case of Axisymmetric Oblate Ellipsoidal Cavities. J Eng Mater Technol 1994;116:290–7. [33] Wen J, Huang Y, Hwang KC, Liu C, Li M. The modified Gurson model accounting for the void size effect. Int J Plast 2005;21:381–95. [34] Benzerga AA, Besson J. Plastic potentials for anisotropic porous solids. Eur J Mech ASolids 2001;20:397–434. [35] Nahshon K, Hutchinson JW. Modification of the Gurson Model for shear failure. Eur J Mech - ASolids 2008;27:1–17. [36] Brunet M, Sabourin F, Mguil-Touchal S. The Prediction of Necking and Failure in 3 D. Sheet Forming Analysis Using Damage Variable. J Phys IV Colloq 1996;06:C6-473-C6482. [37] Abbassi F, Pantalé O, Mistou S, Zghal A, Rakotomalala R. Effect of ductile damage evolution in sheet metal forming: experimental and numerical investigations. Key Eng Mater 2010;446:157–69. [38] Kachanov L. Time of the rupture process under creep conditions. Izv Akad Nauk SSSR Otd Teckhnicheskikh Nauk 1958;8:26–31. [39] Chaboche JL. Continuous damage mechanics — A tool to describe phenomena before crack initiation. Nucl Eng Des 1981;64:233–47. [40] Lemaitre J. A continuous damage mechanics model for ductile fracture. J Eng Mater Technol 1985;107:83–89. [41] Murakami S. Mechanical Modeling of Material Damage. J Appl Mech 1988;55:280–6. [42] Teixeira P, Santos AD, Andrade Pires FM, César de Sá JMA. Finite element prediction of ductile fracture in sheet metal forming processes. J Mater Process Technol 2006;177:278–81. [43] Bahloul R. Optimisation of process parameters in flanging operation in order to minimise stresses and Lemaitre’s damage. Mater Des 2011;32:108–20. [44] Maria D, Kerim I, Till C, Sebastian H, Helmut R, Erman TA. Material characterization and validation studies for modeling ductile damage during deep drawing. Procedia Eng 2017;183:77–82. [45] Cao TS, Gaillac A, Montmitonnet P, Bouchard PO. Identification methodology and comparison of phenomenological ductile damage models via hybrid numerical– experimental analysis of fracture experiments conducted on a zirconium alloy. Int J Solids Struct 2013;50:3984–99. [46] Hogström P, Ringsberg JW, Johnson E. An experimental and numerical study of the effects of length scale and strain state on the necking and fracture behaviours in sheet metals. Int J Impact Eng 2009;36:1194–203. 35

[47] Al Thairy H, Wang YC. A numerical study of the behaviour and failure modes of axially compressed steel columns subjected to transverse impact. Int J Impact Eng 2011;38:732–44. [48] Ribeiro J, Santiago A, Rigueiro C. Damage model calibration and application for S355 steel. Procedia Struct Integr 2016;2:656–63. [49] Giglio M, Manes A, Viganò F. Ductile fracture locus of Ti–6Al–4V titanium alloy. Int J Mech Sci 2012;54:121–35. [50] Yu H, Jeong D, Gordon J, Tang Y. Analysis of impact energy to fracture unnotched charpy specimens made from railroad tank car steel. Proc 2007 ASME Rail Transp Div Fall Tech Conf RTDF 2007. [51] Yu HL, Jeong DY. Application of a stress triaxiality dependent fracture criterion in the finite element analysis of unnotched Charpy specimens. Theor Appl Fract Mech 2010;54:54–62. [52] Hillerborg A, Modéer M, Petersson P-E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6:773–81. [53] ABAQUS. User’s Manual (version 6.10) 2010. [54] Yamada Y, Endo H, Pedersen PT. Numerical study on the effect of buffer bow structure in ship-ship collision. Proc Int Offshore Polar Eng Conf 2005;2005:604–11. [55] Dunand M, Mohr D. Hybrid experimental–numerical analysis of basic ductile fracture experiments for sheet metals. Int J Solids Struct 2010;47:1130–43. [56] Kacem A, Krichen A, Manach P-Y. Occurrence and effect of ironing in the hole-flanging process. J Mater Process Technol 2011;211:1606–13. [57] Hill R. A Theory of the Yielding and Plastic Flow of Anisotropic Metals. Proc R Soc Lond Math Phys Eng Sci 1948;193:281–97. [58] Chung K, Ma N, Park T, Kim D, Yoo D, Kim C. A modified damage model for advanced high strength steel sheets. Int J Plast 2011;27:1485–511.

36







The coupled approach offered by Abaqus code and based on the effective stress concept including the choice of element removal for modeling progressive damage is considered. A calibration procedure for identifying the parameters of the damage model is developed through experiments and numerical simulations for uniaxial tensile and hole-expansion tests. A flow chart is proposed in this manuscript as a practical tool to describe the proposed calibration procedure.

37