Determination of Forming Limit Diagrams Based on Ductile Damage Models and Necking Criteria

Abstract

In this paper, forming limit diagrams (FLDs) for an aluminum alloy are predicted through numerical simulations using various localized necking criteria. A comparative study is conducted for the FLDs determined by using the Lemaitre damage approach and those obtained with the modified Gurson-Tvergaard-Needleman (GTN) damage model. To this end, both damage models coupled with elasto-plasticity and accounting for plastic anisotropy have been implemented into the ABAQUS/Explicit software, through the user-defined subroutine VUMAT, within the framework of large plastic strains and a fully three-dimensional formulation. The resulting constitutive frameworks are then combined with four localized necking criteria to predict the limit strains for an AA6016-T4 aluminum alloy. Three of these necking criteria are based on finite element (FE) simulations of the Nakazima deep drawing test with various specimen geometries, while the fourth criterion is based on bifurcation theory. The simulation results reveal that the limit strains predicted by local criteria, which are based on FE simulations of the Nakazima test, are in good agreement with the experiments for a number of strain paths, while those obtained with the bifurcation analysis provide an upper bound to the experimental FLD.

Keywords:
GTN model; Lemaitre damage approach; plastic anisotropy; forming limit diagram; localized necking criteria; bifurcation theory; Nakazima deep drawing test

1 INTRODUCTION

The good knowledge of the formability of metallic materials is very important for the successful forming of sheet metals. The concept of forming limit diagram (FLD), which was originally introduced by Keeler and Backofen (1963), and later by Goodwin (1968), has been the most widely used tool for the characterization of the formability of sheet metals. This strategy allows delimiting the limit strains for stretched sheet metals that should not be exceeded in order to ensure a good quality of the final product. The FLDs are determined using necking or fracture criteria, which may be based on sound theoretical developments (see, e.g., Hill, 1952; Stören and Rice, 1975; Yamamoto, 1978; Abed-Meraim et al., 2014) or on finite element (FE) simulations (see, e.g., Zhang et al., 2011; Lumelskyy et al., 2012, Martínez-Donaire et al., 2014; Kami et al., 2015). These criteria are generally combined with constitutive models for the prediction of limit strains in sheet metal forming. In order to describe the behavior of sheet metals in a realistic way, advanced elastic-plastic models coupled with damage have been developed in the literature, which can be classified into two main approaches, namely the micromechanics damage modeling (MDM) and the continuum damage mechanics (CDM). The MDM approach has been first developed by Gurson (1977), who considered the initiation of damage as the growth of micro-voids in porous materials surrounded by a rigid-plastic matrix. This preliminary modeling approach has subsequently received a number of extensions to obtain the so-called Gurson-Tvergaard-Needleman (GTN) damage model (see, e.g., Tvergaard, 1982a,1982b; Tvergaard and Needleman, 1984) in order to account for all damage mechanisms (i.e., nucleation, growth and coalescence of voids) as well as the hardening of the dense matrix. Concurrently, the CDM approach has been introduced by the works of Kachanov (1958), and extended later in the framework of irreversible thermodynamics (see, e.g., Lemaitre, 1992; Chaboche, 1999). In the CDM approach, the damage variable represents the surface density of microcracks across a given plane, and may be modeled as an isotropic scalar variable (see, e.g., Lemaitre, 1985, 1992), or a tensor variable for anisotropic damage (see, e.g., Lemaitre et al., 2000; Brünig, 2003).

In this work, both the MDM and the CDM approaches are adopted for the modeling of ductile damage in sheet metals. More specifically, a classical elastic-plastic model with anisotropic plasticity is coupled with the Lemaitre damage theory, while the GTN model is combined with the Hill (1948) anisotropic yield surface to account for plastic anisotropy. The resulting models are implemented into the finite element code ABAQUS/Explicit, through the user-defined subroutine VUMAT, within the framework of large plastic strains and a fully three-dimensional formulation. Each of these models is then used for the FE simulations of the Nakazima deep drawing test with different specimen geometries in order to predict FLDs. The latter are determined using four different criteria for the onset of localized necking. Three of these criteria are based on the FE simulations of the Nakazima deep drawing test, while the fourth one is based on bifurcation theory (see, e.g., Rudnicki and Rice, 1975; Stören and Rice, 1975; Rice, 1976). The numerical FLDs obtained with the current approach are compared with the experimental results taken from Kami et al. (2015).

2 MODELING OF DUCTILE DAMAGE

In this section, the constitutive equations associated with both the GTN damage model and the Lemaitre damage theory are described. Note that both modeling approaches are developed within the framework of large strains and three-dimensional formulation.

2.1 GTN Damage Model

Gurson (1977) proposed a yield condition depending on the void volume fraction, which represents the density of micro-defects within the material. Subsequently, this model has been improved by Tvergaard (1981) and Tvergaard and Needleman (1984) in order to take into account the interaction between voids. The resulting modifications led to the definition of the following plastic yield surface:

${\Phi }_{GTN}={\left(\frac{{\sigma }_{eq}}{Y}\right)}^2+2q_1{f}^{*}cosh\left(\frac{3}{2}{q}_2\frac{{\sigma }_m}{Y}\right)-\left(1+q_3{f}^{*}{}^2\right)\le 0$ (1)

where ${\sigma }_{eq}(\sigma )=\sqrt{3{\sigma }^{\prime }:{\sigma }^{\prime }/2}$ is the macroscopic von Mises equivalent stress, ${\sigma }_m=\frac{1}{3}\sigma :1$ is the hydrostatic stress, _{^{σ’ = σ - σm}} 1 is the deviatoric part of the Cauchy stress σ, with 1 being the second-order identity tensor. The isotropic hardening of the fully dense matrix is described by the variable $$Y\left({\overline{\epsilon }}^{pl}\right)$$, function of the equivalent plastic strain $${\overline{\epsilon }}^{pl}$$. The parameters q ₁, q ₂ and q ₃, introduced by Tvergaard (1981, 1982a), account for void interaction effects. The void coalescence mechanism is considered through the introduction of an effective void volume fraction f ^*(f) (see, e.g., Tvergaard, 1982b; Tvergaard and Needleman, 1984). This function is defined as

$f^{*}(f)=\{\begin{array}{l}f\text{for}f\le f_{cr},\\ f_{cr}+\left(\frac{f_u^{*}-f_{cr}}{f_R-f_{cr}}\right)\left(f-f_{cr}\right)\text{for}{f}_{cr}<f\le f_R,\end{array}$ (2)

where f represents the actual void volume fraction, $$f_u^{*}$$ is the ultimate value of f ^*, while _^fR is the void volume fraction at fracture. The void coalescence phenomenon occurs when the void volume fraction reaches the critical value _^fcr .

In order to account for the plastic anisotropy of the material, the GTN plastic yield surface (see Eq. (1)) is modified by introducing the Hill (1948) equivalent stress instead of the von Mises one (see, e.g., Chen and Butcher, 2013; Kami et al., 2015; Li et al., 2015). The corresponding expression of equivalent stress is given by

${\sigma }_{eq}(\sigma )=\sqrt{{\sigma }^{\prime }:M:{\sigma }^{\prime }}$ (3)

where the fourth-order tensor M contains the six anisotropy coefficients of the Hill (1948) quadratic yield criterion. It is worth noting that the original isotropic GTN model is recovered from the anisotropic one when the Lankford coefficients r ₀, r ₄₅ and r ₉₀ are set to 1.

Based on the principle of equivalence in plastic work rate (Gurson, 1977), the equivalent plastic strain rate $${\dot{\overline{\epsilon }}}^{pl}$$ of the fully dense matrix material is obtained as follows:

${\dot{\overline{\epsilon }}}^{pl}=\frac{\sigma :D^p}{\left(1-f\right)Y}$ (4)

where D^p is the macroscopic plastic strain rate tensor. The latter is defined using the classical normality rule with respect to the yield function, and is expressed as

$D^p=\dot{\lambda }{V}_{GTN}$ (5)

where $\dot{\lambda }$ is the plastic multiplier, and V_GTN = ∂Φ_GTN / ∂σ is the direction of the plastic flow. Isotropic hardening for the dense matrix is assumed in this work, which is defined by the following expression:

$\dot{Y}=h{\dot{\overline{\epsilon }}}^{pl}$ (6)

where $h\left({\overline{\epsilon }}^{pl}\right)$ is the plastic hardening modulus of the fully dense matrix material.

The evolution of the void volume fraction is based on both nucleation of new voids and growth of existing voids (see, e.g., Chu and Needleman, 1980). The associated evolution equation is given by

$\dot{f}={\dot{f}}_{growth}+{\dot{f}}_{nucleation}$ (7)

where

${\dot{f}}_{growth}=\left(1-f\right)D^p:1$ (8)

For the void nucleation, the latter is assumed to be strain controlled in this work. The expression of the nucleation rate is given by

${\dot{f}}_{nucleation}=A_N{\dot{\overline{\epsilon }}}^{pl}$ (9)

where _^AN has been defined in Chu and Needleman (1980) by the following normal distribution law:

$A_N=\frac{f_N}{s_N\sqrt{2\pi }}\exp \left[-\frac{1}{2}{\left(\frac{{\overline{\epsilon }}^{pl}-{\epsilon }_N}{s_N}\right)}^2\right]$ (10)

where _^fN is the volume fraction of the inclusions that are likely to nucleate, _^(N is the mean equivalent plastic strain of nucleation, and _^SN is the corresponding standard deviation.

The consistency condition for the GTN model, which ensures plastic loading, may be written in the following form:

${\dot{\Phi }}_{GTN}=V_{GTN}:\dot{\sigma }+V_Y\dot{Y}+{\delta }_{GTN}{V}_{f^{*}}\dot{f}=0$ (11)

where

$\{\begin{array}{l}{V}_{GTN}=\frac{\partial {\Phi }_{GTN}}{\partial \sigma }=2\frac{M:{\sigma }^{\prime }}{Y^2}+q_1{q}_2\frac{f^{*}}{Y}\sinh \left(\frac{3}{2}{q}_2\frac{{\sigma }_m}{Y}\right)1\\ V_Y=\frac{\partial {\Phi }_{GTN}}{\partial Y}=-2\frac{{\sigma }_{eq}^2}{Y^3}-3q_1{q}_2{f}^{*}\frac{{\sigma }_m}{Y^2}\sinh \left(\frac{3}{2}{q}_2\frac{{\sigma }_m}{Y}\right)\\ V_{f^{*}}=\frac{\partial {\Phi }_{GTN}}{\partial f^{*}}=2q_1\cosh \left(\frac{3}{2}{q}_2\frac{{\sigma }_m}{Y}\right)-2q_3{f}^{*}\end{array}$ (12)

And

${\delta }_{GTN}=\{\begin{array}{l}1forf\le f_{cr}\\ \frac{f_u^{*}-f_{cr}}{f_R-f_{cr}}for{f}_{cr}<f\le f_R\end{array}$ (13)

In the co-rotational frame, which is associated with the Jaumann objective derivative, the Cauchy stress-strain relationship is obtained using the classical hypoelastic law defined by

$\dot{\sigma }=C^e:\left(D-D^p\right)=C_{GTN}^{ep}:D$ (14)

where _^Ce is the fourth-order elasticity tensor, $C_{GTN}^{ep}$ is the elastic-plastic tangent modulus for the GTN model, and D is the strain rate tensor.

By substituting Eqs. (4)-(9) into the consistency condition (11), the expression of the plastic multiplier $\dot{\lambda }$ is obtained as follows:

$\dot{\lambda }=\frac{V_{GTN}:C^e:D}{H_{GTN}}$ (15)

where

$H_{GTN}=V_{GTN}:C^e:V_{GTN}-\frac{1}{\left(1-f\right)}\frac{\sigma :V_{GTN}}{Y}\left[h{V}_Y+A_N{V}_{f^{*}}{\delta }_{GTN}\right]-\left(1-f\right){\delta }_{GTN}{V}_{f^{*}}{V}_{GTN}:1$ (16)

The analytical expression of the elastic-plastic tangent modulus for the GTN model is obtained by substituting the expression of the plastic multiplier (Eq. (15)) into the hypoelastic relationship (14)

$C_{GTN}^{ep}=C^e-\alpha \frac{\left(C^e:V_{GTN}\right)\otimes \left(V_{GTN}:C^e\right)}{H_{GTN}}$ (17)

where α = 0 for elastic loading/unloading, and α = 1 for strict plastic loading.

2.2 Lemaitre Damage Model

The second approach to ductile damage considered in this work is based on the works of Lemaitre (1985, 1992), which was originally introduced by Kachanov (1958). In the literature, this approach is referred to as Continuum Damage Mechanics (CDM) theory, which provides a phenomenological description for ductile damage, in contrast to the micromechanics-based Gurson damage model. In the CDM theory, the damage variable, which may be scalar isotropic or tensor-valued anisotropic, represents the surface density of microdefects. In the current work, the adopted elasto-plastic model coupled with ductile damage takes into account the initial anisotropy of the material, using the Hill’48 quadratic yield function, while hardening is taken to be isotropic.

Based on the strain equivalence principle (Lemaitre and Chaboche, 1978), the material behavior is affected by continuum damage through the introduction of an effective stress tensor $\tilde{\sigma }$ given by

$\tilde{\sigma }=\frac{\sigma }{1-d}$ (18)

where the scalar damage variable d varies between 0 to 1, with d = 0 for an undamaged material, and d = 1 for a fully damaged material.

The plastic yield function F is written in the following form:

$F={\sigma }_{eq}\left(\tilde{\sigma }\right)-Y\le 0$ (19)

where ${\sigma }_{eq}\left(\tilde{\sigma }\right)=\sqrt{{\tilde{\sigma }}^{\prime }:M:{\tilde{\sigma }}^{\prime }}$ is the Hill’48 effective equivalent stress, and ${\tilde{\sigma }}^{\prime }$ is the deviatoric part of the effective stress.

The plastic flow rule is given by the normality law, which defines the plastic strain rate tensor D^P as

$D^p=\dot{\lambda }\frac{\partial F}{\partial \sigma }=\dot{\lambda }{\tilde{V}}_{CDM}=\frac{1}{1-d}\dot{\lambda }{V}_{CDM}$ (20)

where

${\tilde{V}}_{CDM}=\frac{1}{1-d}\frac{M:{\tilde{\sigma }}^{\prime }}{{\sigma }_{eq}}$ (21)

With a special choice of co-rotational frame, which is associated with the Jaumann objective derivative, the constitutive relation is written in the following form:

$\dot{\tilde{\sigma }}=C^e:\left(D-D^p\right)\Rightarrow \dot{\sigma }=\left(1-d\right)C^e:\left(D-D^p\right)-\frac{\dot{d}}{1-d}\sigma$ (22)

The evolution law for the damage variable is expressed by the following equation:

$\dot{d}=H_d\dot{\lambda }=\{\begin{array}{l}\frac{1}{{\left(1-d\right)}^{\beta }}{\left(\frac{Y^e-Y_i^e}{S}\right)}^s\dot{\lambda }if{Y}^e\ge Y_i^e\\ 0otherwise\end{array}$ (23)

where _^Ye is the strain energy density release rate (see, e.g., Lemaitre, 1992; Lemaitre et al., 2000). Its expression is given, in the case of linear isotropic elasticity, by the following relationship:

$Y^e=\frac{J_2{}^2}{2E}\left[\frac{2}{3}\left(1+\nu \right)+3\left(1-2\nu \right){\left(\frac{{\tilde{\sigma }}^H}{J_2}\right)}^2\right]$ (24)

where $J_2=\sqrt{\frac{3}{2}{\tilde{\sigma }}^{\prime }:{\tilde{\sigma }}^{\prime }}$ is the equivalent effective stress in the sense of von Mises, ${\tilde{\sigma }}^H=\frac{1}{3}\tilde{\sigma }:1$ is the hydrostatic effective stress, while E and v denote, respectively, the Young modulus and the Poisson ratio.

It is easy to show that the expression of the Cauchy stress rate tensor given by Eq. (22) can be rewritten in the form

$\dot{\sigma }=C_{CDM}^{ep}:D$ (25)

where $C_{CDM}^{ep}$ is the elastic-plastic tangent modulus for the Lemaitre damage model.

The consistency condition $\dot{F}=0$ allows the determination of the plastic multiplier $\dot{\lambda }$, which writes

$\dot{\lambda }=\frac{V_{CDM}:C^e:D}{H_{CDM}}$ (26)

where

$H_{CDM}=V_{CDM}:C^e:{\tilde{V}}_{CDM}+H_Y$ (27)

where _^HY is the scalar isotropic hardening modulus, which governs the evolution of isotropic hardening (i.e., $\dot{Y}=H_Y\dot{\lambda }$). By substituting the expression of the plastic multiplier (Eq. (26)) into the hypoelastic relationship (22), the elastic-plastic tangent modulus $C_{CDM}^{ep}$ for the Lemaitre damage model is given by

$C_{CDM}^{ep}=\left(1-d\right)C^e-\alpha \left(\frac{\left[C^e:V_{CDM}+H_d\tilde{\sigma }\right]\otimes \left(V_{CDM}:C^e\right)}{H_{CDM}}\right)$ (28)

where α = 1 for strict plastic loading and 0 otherwise.

It is worth noting that the tangent moduli for both the GTN damage model and the Lemaitre damage model (i.e., Eqs. (17) and (28), respectively) are only required for the bifurcation analysis, which will be detailed in Section 5.

2.3 Numerical Implementation of the Constitutive Equations

In this work, both the modified anisotropic GTN model and the Lemaitre damage model are implemented into the commercial finite element code ABAQUS/Explicit via the user-defined material subroutine VUMAT. The same explicit time integration scheme is used for both damage models, which is based on the fourth-order Runge-Kutta method. This algorithm allows updating the stress state and all of the internal state variables at the end of the loading increment starting from a known state at the beginning of the loading increment. This time integration scheme has the advantage of being straightforward and robust, since no iterative procedure is needed, unlike implicit time integration schemes (see, e.g., Aravas, 1987). However, the time increment must be kept small enough to ensure accuracy and stability (see, e.g., Li and Nemat-Nasser, 1993; Kojic, 2002).

It can be shown that the evolution equations for both the GTN damage model and the Lemaitre damage model can be written in the following compact form of general differential equation:

$\dot{u}=h_u\left(u\right)$ (29)

where u encompasses all of the internal variables and stress state, while vector hu(u) includes all evolution laws for each damage model. The above condensed differential equation is then integrated over each loading increment, using the forward fourth-order Runge-Kutta time integration scheme. The resulting algorithms for both damage models are implemented into the finite element code ABAQUS/Explicit, via VUMAT user-defined material subroutines, within the framework of large strains and a fully three-dimensional formulation.

3 DETERMINATION OF MATERIAL PARAMETERS

The material considered in this work is the AA6016-T4 aluminum sheet (see Kami et al., 2015). For this material, the experimental FLD and the material parameters corresponding to the anisotropic GTN damage model have been determined by Kami et al. (2015). In the latter reference, the Swift isotropic hardening law has been considered, which is defined by the following expression:

$Y=k{\left({\epsilon }_0+{\overline{\epsilon }}^{pl}\right)}^n$ (30)

where k, ( ₀ and n are the hardening parameters. The associated elastic-plastic parameters are summarized in Table 1.

Table 1
Elastic properties and Swift’s hardening parameters.

Material	E (MPa)	ν	k (MPa)	ε0	n
AA6016-T4	70,000	0.33	525.77	0.011252	0.2704

As mentioned in Section 2.1, the GTN yield surface has been modified to account for the planar plastic anisotropy. The Hill’48 anisotropy coefficients are determined using the Lankford coefficients r ₀, r ₄₅ and r ₉₀, which were identified by Kami et al. (2015) on the basis of three uniaxial tensile tests performed along three sheet orientations, namely 0°, 45° and 90° with respect to the rolling direction. The corresponding r-values are reported in Table 2. For the GTN damage parameters, the latter were identified in Kami et al. (2015) using an inverse identification procedure that combines the response surface methodology and the simulation of a uniaxial tensile test. The associated parameters are listed in Table 3.

Table 2
r-values for the AA6016-T4 aluminum sheet.

Material	r 0	r 45	r 90
AA6016-T4	0.5529	0.4091	0.5497

Table 3
Damage parameters for the GTN damage model.

Material	f 0	SN	εN	fN	fcr	fR	q 1	q 2	q 3
AA6016-T4	0.00035	0.10	0.30	0.05	0.05	0.15	1.5	1.0	2.25

For the Lemaitre damage model, the modeling of the material hardening and the description of the plastic anisotropy are taken the same as in the case of the GTN damage model. In the current contribution, the Lemaitre damage parameters are calibrated using an inverse identification procedure along with the experimental uniaxial tensile test for the AA6016-T4 aluminum sheet. This inverse identification strategy is based on least-squares minimization of the difference between the experimental and numerical load-displacement response for a uniaxial tensile test. The geometric dimensions and the boundary conditions for the uniaxial tensile specimen are all specified in Figure 1 (see Kami et al., 2014). The specimen is discretized using the eight-node reduced integration solid finite element (C3D8R) from ABAQUS, with an initial mesh size of 0.5 mm. The identified values of the Lemaitre damage model are summarized in Table 4.

Figure 1
Geometry and boundary conditions for the uniaxial tensile specimen.

Table 4
Damage parameters for the Lemaitre damage model.

Material	β	S	s	$Y_i^e$
AA6016-T4	12	4.0	1.3	0.0

In order to better emphasize the identification results as well as the performance of the numerical implementation of both damage models, Figure 2 compares the simulated load-displacement responses, obtained using both damage models, with the experimental counterpart given in Kami et al. (2015). This figure clearly shows that the simulated responses for both damage models are in very good agreement with the experimental curve and, in particular, demonstrates the ability of the implemented models to reproduce the sudden load drop that precedes the final fracture.

Figure 2
Tensile load-displacement responses simulated with the GTN and the Lemaitre damage models, along with the experimental curve taken from Kami et al. (2015).

4 FINITE ELEMENT MODEL

4.1 Nakazima Deep Drawing Test

The implemented GTN and Lemaitre damage models, with their associated material parameters, are used in the simulation of the Nakazima deep drawing test (see Figure 3) in order to determine the FLD of the AA6016-T4 aluminum sheet. The geometric parameters for the Nakazima deep drawing process are summarized in Table 5. According to the standard procedure described in ISO 12004-2 (2008), seven specimens with different geometries are considered in the simulations. Each specimen allows for the reproduction of a particular strain path, which is typically encountered in sheet metal forming processes. The general geometry for a given specimen width is illustrated in Figure 4. The seven specimens are designed by varying the width parameter W from 30 mm to 185 mm, which leads to different strain paths in the central part of the specimens, ranging from uniaxial tension to balanced biaxial tension.

Figure 3
FE representation of the Nakazima deep drawing test with sheet specimen of 70 mm width.

Figure 4
Specimen geometry and dimensions used in the Nakazima deep drawing test.

Table 5
Geometric parameters for the Nakazima deep drawing test.

Punch diameter	100 mm
Die opening diameter	110 mm
Die profile radius	10 mm
Initial sheet thickness	1 mm

During the Nakazima deep drawing simulation, each specimen is clamped all around its circumferential edges so that the material flow is prevented. The Coulomb friction coefficient is taken equal to 0.03 for the contact between the punch and the specimen, while it is taken equal to 0.1 for the contact between the die, the holder and the specimen (see Kami et al., 2015). In addition, a holding force of 100 kN is applied during the forming process. The forming tools are modeled as discrete rigid bodies, while the blank is modeled with the eight-node three-dimensional continuum finite element with reduced integration (C3D8R), which is available in the ABAQUS/Explicit software. Note that particular attention has been paid to optimizing the mesh of the blank, as illustrated in Figure 5. Indeed, the central part of the blank, which is subjected to large plastic strains, is discretized with a fine mesh, while the rest of the blank is discretized with a coarse mesh. Moreover, in order to save computational time, the built-in ABAQUS mass scaling technique is used in what follows, with a target time step of 10^-6 s and time period of 1 s. These numerical parameters, which lead to reasonable computation times, are selected so that the simulation of the Nakazima deep drawing test is achieved under conditions that are quite similar to those of a quasi-static analysis.

Figure 5
FE mesh used for the specimen of 70 mm width.

4.2 Mesh Sensitivity Analysis

As already pointed out by several authors in the literature (see, e.g., Tvergaard and Needleman, 1984; Besson et al., 2001; Peerlings et al., 2001), it is nowadays well known that the mesh size has an important influence on the occurrence of strain localization, and particularly for behavior models exhibiting damage-induced softening. In order to analyze the effect of the mesh size on the numerical results, several meshes are used in the simulation of the Nakazima deep drawing test with the specimen of 30 mm width.

First, the effect of the number of elements in the thickness direction is analyzed by considering, successively, two, three, four and five element layers. Note that for these four different through-thickness FE discretizations, the same in-plane mesh discretization is used, which consists of 0.5×0.5 mm² in the central part of the specimen. Figures 6 and 7 show the effect of the number of elements in the thickness direction on the evolution of the thickness strain and the punch force-displacement response, respectively. These figures reveal that the number of element layers in the thickness direction has a very small effect on the evolution of the thickness strain, while it has a relatively more noticeable effect on the maximum punch force.

Figure 6
Effect of the number of elements in the thickness direction on the evolution of the thickness strain during the Nakazima test with the specimen of 30 mm width.

Figure 7
Effect of the number of elements in the thickness direction on the punch force-displacement response for the Nakazima test with the specimen of 30 mm width.

Then, the impact of the in-plane mesh size on the evolution of the thickness strain and the punch force-displacement response is investigated. To this end, the Nakazima deep drawing test is performed again for the specimen of 30 mm width, using four different in-plane meshes for the central part of the specimen. These in-plane FE discretizations represent coarse, intermediate, fine and very fine meshes, which correspond to mesh sizes of 1×1 mm², 0.75×0.75 mm², 0.5×0.5 mm² and 0.25×0.25 mm², respectively. Note that, for these four different in-plane mesh discretizations, four element layers in the thickness direction are used. Figures 8 and 9 show the simulated results that are obtained with the four in-plane mesh sizes. Similar to the effect of the number of elements through the thickness, the in-plane mesh size has a small effect on the evolution of the thickness strain, while it has a relatively more perceptible effect on the maximum punch force and on the final punch stroke (i.e., after the sudden load drop).

Figure 8
Effect of the in-plane mesh size on the evolution of the thickness strain during the Nakazima test with the specimen of 30 mm width.

Figure 9
Effect of the in-plane mesh size on the punch force-displacement response for the Nakazima test with the specimen of 30 mm width.

In conclusion, the above mesh sensitivity analysis suggests using the fine in-plane mesh (i.e., 0.5×0.5 mm²) with four element layers through the thickness in the remaining simulations of the current study. Indeed, this choice appears as a pragmatic compromise in terms of accurate description of the various nonlinear phenomena and reasonable computational times.

5 LOCALIZED NECKING CRITERIA

In this section, four necking criteria are presented, which will be subsequently used for the prediction of strain localization in sheet metals. Three of these criteria are based on the FE analysis of deep drawing process, while the fourth one is based on bifurcation theory. Note that, none of these criteria requires the introduction of additional user-defined parameters, in contrast to the Marciniak and Kuczynski (1967) criterion.

5.1 Criterion of Maximum Second Time Derivative of Thickness Strain

This criterion is based on the analysis of the evolution of thickness strain during the Nakazima deep drawing test. More specifically, the onset of strain localization is associated with the maximum of the thickness strain acceleration, which is obtained by computing the second time derivative of thickness strain in the localized zone (see, e.g., Situ et al., 2006, 2007, 2011; Zhang et al., 2011; Lumelskyy et al., 2012; Martínez-Donaire et al., 2014). After the maximum in the second time derivative of thickness strain (i.e., thickness strain acceleration) is reached, the localized thinning in the sheet proceeds gradually until the onset of fracture. Based on this numerical criterion, Figure 10 shows an illustration of the onset of localized necking during the simulation of the Nakazima deep drawing test with the specimen of 30 mm width.

Figure 10
Prediction of localized necking based on the maximum of the 2^nd time derivative of thickness strain.

5.2 Criterion Based on the Ratio of Equivalent Plastic Strain Increment

In order to determine the onset of localized necking, a criterion based on the ratio of equivalent plastic strain increment is used as numerical necking criterion (see, e.g., Narasimhan and Wagoner, 1991; Chung et al., 2014). This ratio of equivalent plastic strain increment is associated with two elements: a critical element and its neighboring element. More specifically, the critical element (referred to here as element B) is preliminarily identified during the Nakazima test, which is generally located in the central part of the specimen that is in contact with the punch. Then, the neighboring element is also identified (referred to here as element A), which is located five elements away from the critical element along the rolling direction. With elements A and B thus identified, the onset of localized necking is detected when the ratio of equivalent plastic strain increment in element B to that in element A becomes larger than 10, as illustrated in Figure 11.

Figure 11
Prediction of localized necking based on the ratio of equivalent plastic strain increment.

5.3 Maximum Punch Force Criterion

Several theoretical criteria based on the maximum force principle have been developed in the literature for the prediction of diffuse or localized necking in sheet metals (see, e.g., Swift, 1952; Hora et al., 1996; Mattiasson et al., 2006). Based on these earlier contributions, the maximum in the punch force-displacement response during the simulation of the Nakazima test is used here as numerical criterion for the prediction of localized necking (see the illustration in Figure 12).

Figure 12
Prediction of localized necking based on the maximum in the punch force-displacement response.

5.4 Loss of Ellipticity Criterion

In contrast to the three numerical criteria presented above, a more theoretically-based criterion is proposed here for the prediction of localized necking in sheet metals, which is based on bifurcation theory. This criterion has been established by Rice and co-workers (see, e.g., Rudnicki and Rice, 1975; Stören and Rice, 1975; Rice, 1976) to predict strain localization in the form of an infinite band in a solid otherwise homogeneous. This approach corresponds to the loss of ellipticity (LE) of the partial differential equations governing the associated boundary value problem. The condition of localization, which may be derived from the Hadamard compatibility condition and the static equilibrium equation, is given by the following relation:

$\det \left(A\right)=\det \left(n\cdot L\cdot n\right)=0$ (31)

where A denotes the so-called acoustic tensor, n is the normal to the localization band, while L represents the tangent modulus that relates the nominal stress rate tensor to the velocity gradient. The expression of the latter is given by the following relationship (see, e.g., Haddag et al., 2009; Mansouri et al., 2014):

$L=C^{ep}+L_1-L_2-L_3$ (32)

where _^Cep is the analytical tangent modulus derived from the constitutive equations, which corresponds to $C_{GTN}^{ep}$ for the GTN damage model, and $C_{CDM}^{ep}$

for the Lemaitre damage model (see Eqs. (17) and (28), respectively). The fourth-order tensors L₁, L₂ and L₃, which only depend on Cauchy stress components, result from the large strain framework. Their detailed expressions can be found in Haddag et al. (2009).

The loss of ellipticity condition given by Eq. (31) is numerically assessed by computing the determinant of the acoustic tensor A for each loading increment. The numerical detection of strain localization is achieved when the minimum of the determinant of the acoustic tensor A, over all possible orientations for the normal n to the localization band, becomes non-positive, as illustrated in Figure 13.

Figure 13
Prediction of localized necking based on the loss of ellipticity criterion.

It is worth noting that the LE criterion is based on a three-dimensional bifurcation analysis from a homogeneous pre-localization state. This state of uniform deformation is achieved by considering a single finite element with one integration point, which is subjected to various linear strain paths that are typically applied to sheet metals under in-plane biaxial stretching (i.e., ranging from uniaxial tension to balanced biaxial tension).

6 PREDICTION OF FLDS AND COMPARISON WITH EXPERIMENTS

The necking criteria presented in the previous section are combined here with both the GTN and the Lemaitre damage models to predict the FLDs of the AA6016-T4 aluminum sheet. Figure 14 shows a comparison of the FLDs predicted by the four necking criteria, along with the numerical and experimental FLDs provided by Kami et al. (2015). It is worth noting that the numerical FLD given in Kami et al. (2015) is determined on the basis of another specific procedure, which is very different from the numerical methods and criteria adopted in the current work. Indeed, in Kami et al.’s (2015) numerical FLD determination, which closely mimics their experimental FLD determination, the standard specification of the ISO 12004-2 (2008) is followed, where the numerical strain fields and their experimental counterparts are analyzed by the ARAMIS software to determine the numerical and experimental FLDs.

Figure 14
Prediction of FLDs using the four localized necking criteria along with the FLDs provided by Kami et al. (2015): (a) maximum of the 2^nd time derivative of thickness strain, (b) ratio of equivalent plastic strain increment, (c) maximum of punch force, and (d) loss of ellipticity.

On the whole, the FLDs predicted by the two damage models are close to each other, and are also comparable, in terms of order of magnitude, to the numerical and experimental FLDs provided by Kami et al. (2015). More specifically, the limit strains obtained with the criterion based on the maximum of the 2^nd time derivative of thickness strain are in good agreement with the experimental results in the left-hand side of the FLD (see Figure 14a), while these limit strains are well predicted by the criterion of equivalent plastic strain increment ratio in the neighborhood of the plane-strain tension path (see Figure 14b). However, for the two above-discussed criteria, the predicted limit strains are overall underestimated in the right-hand side of the FLD, which is probably due to the material parameter identification and, particularly, to the identification of damage parameters. Indeed, the latter are identified using only one type of mechanical tests (i.e., a uniaxial tension test) for both the GTN and the Lemaitre damage models, which results in non-negligible error in the right-hand side of the FLD. It is now widely recognized that the accurate calibration of material parameters requires an identification procedure that is based on various types of mechanical tests (i.e., standard uniaxial tension test, plane-strain tension test, Bulge test …), and/or on heterogeneous mechanical tests. Such advanced identification techniques are likely to improve the reliability of the material parameters for various strain paths and, in turn, the accuracy of the corresponding FLD predictions. For the maximum punch force criterion, the predicted FLDs are markedly different from those obtained by the two previous necking criteria, and even the shape of the predicted FLDs does not seem to be usual (see Figure 14c). Indeed, the punch force represents some averaged information during the forming process, and its use to detect local phenomena, such as localized necking, does not seem to be suitable. For the LE criterion, the FLDs predicted by the two damage models are overestimated for almost all strain paths, except in the extreme right-hand side of the FLD, where the limit strains are rather underestimated. It is worth noting that the LE criterion is based on a three-dimensional bifurcation analysis from a homogeneous pre-localization state, with the only consideration of material instability, without taking into account any structural (geometric) effects. Consequently, the FLDs predicted by the LE criterion are expected to set an upper bound to the experimental ones, which is observed here indeed for most strain paths.

7 CONCLUSIONS

In this work, four different necking criteria have been proposed and compared for the prediction of FLDs for the AA6016-T4 aluminum alloy. For the material constitutive modeling, two approaches to ductile damage have also been considered: the Lemaitre continuum damage theory and the GTN damage description, which was extended to the Hill’48 quadratic yield surface to account for the plastic anisotropy of the material. Both damage models have been numerically implemented into the commercial finite element code ABAQUS/Explicit via the user-defined material subroutine VUMAT. The main contributions of the current study and associated conclusions may be summarized as follows:

• The Lemaitre damage parameters have been identified using an inverse identification procedure based on FE fitting of an experimental load-displacement response of a standard tensile test;

• Based on FE simulations of the Nakazima deep drawing test and four different localized necking criteria, numerical FLDs have been determined for the AA6016-T4 aluminum sheet and compared with the experimental FLD. The obtained results suggest that two of the local criteria (i.e., those based on the maximum of the 2^nd time derivative of thickness strain, and the ratio of equivalent plastic strain increment) yield results that are in good agreement with the experiments in the left-hand side of the FLD and in the neighborhood of the plane-strain tension path, while the global criterion based on the maximum punch force does not seem to be suitable to the prediction of localized necking;

• The predictions using the LE criterion provide upper bounds to the classical experimental FLD, which is consistent with the theoretical foundations on which the bifurcation approach is based. On the other hand, this bifurcation approach could be advantageously used to design new materials with improved ductility, by classifying them in terms of formability limits;

• The accuracy of the numerically predicted FLDs with respect to experiments may be improved by considering various mechanical tests in the identification procedure. Indeed, a number of simple and complex strain paths should be included in the identification procedure in order to provide reliable material parameters, thus improving the accuracy of the predicted FLDs.

References

• Abed-Meraim, F., Balan, T., Altmeyer, G. (2014). Investigation and comparative analysis of plastic instability criteria: application to forming limit diagrams. International Journal of Advanced Manufacturing Technology 71:1247-1262.
• Aravas, N. (1987). On the numerical integration of a class of pressure-dependent plasticity models. International Journal for Numerical Methods in Engineering 24:1395-1416.
• Besson, J., Steglich, D., Brocks, W. (2001). Modeling of crack growth in round bars and plane strain specimens. International Journal of Solids and Structures 38:8259-8284.
• Brünig, M. (2003). An anisotropic ductile damage model based on irreversible thermodynamics. International Journal of Plasticity 19:1679-1713.
• Chaboche, J.L. (1999). Thermodynamically founded CDM models for creep and other conditions: CISM courses and lectures, International Centre for Mechanical Sciences. Creep and Damage in Materials and Structures 399:209-283.
• Chen, Z., Butcher, C. (2013). Micromechanics modelling of ductile fracture. Solid Mechanics and Its Applications, Springer, Netherlands.
• Chu, C., Needleman, A. (1980). Void nucleation effects in biaxially stretched sheets. Journal of Engineering Materials and Technology 102:249-256.
• Chun, K., Lee, C., Kim, H. (2014). Forming limit criterion for ductile anisotropic sheets as a material property and its deformation path insensitivity, Part II: Boundary value problems. International Journal of Plasticity 58:35-65.
• Goodwin, G. (1968). Application of strain analysis to sheet metal forming problems in the press shop. SAE Technical Paper No. 680093.
• Gurson, A.L. (1977). Continuum theory of ductile rupture by void nucleation and growth: Part I - yield criteria and flow rules for porous ductile media. ASME Journal of Engineering Materials and Technology 99:2-15.
• Haddag, B., Abed-Meraim, F., Balan, T. (2009). Strain localization analysis using a large deformation anisotropic elastic-plastic model coupled with damage. International Journal of Plasticity 25:1970-1996.
• Hill, R. (1948). A theory of the yielding and plastic flow of anisotropic metals. Proceedings of the Royal Society of London A 193:281-297.
• Hill, R. (1952). On discontinuous plastic states, with special reference to localized necking in thin sheets. Journal of the Mechanics and Physics of Solids 1:19-30.
• Hora, P., Tong, L., Reissner, J. (1996). A prediction method of ductile sheet metal faillure in FE simulation. Numisheet:252-256. Michigan USA.
• ISO 12004-2 (2008). International Organization for Standardization, Metallic Materials - Sheet and Strip - Determination of Forming-Limit Curves, Part 2: Determination of Forming-Limit Curves in the Laboratory, Switzerland.
• Kachanov, L.M. (1958). On Creep Rupture Time. Izv. Acad. Nauk SSSR, Otd. Techn. Nauk 8:26-31.
• Kami, A., Dariani, B.M., Vanini, A.S., Comsa, D.S., Banabic, D. (2014). Application of a GTN damage model to predict the fracture of metallic sheets subjected to deep-drawing. Proceedings of the Romanian Academy, Series A 15:300-309.
• Kami, A., Dariani, B.M., Vanini, A.S., Comsa, D.S., Banabic, D. (2015). Numerical determination of the forming limit curves of anisotropic sheet metals using GTN damage model. Journal of Materials Processing Technology 216:472-483.
• Keeler, S., Backofen, W.A. (1963). Plastic instability and fracture in sheets stretched over rigid punches. ASM Transactions Quarterly 56:25-48.
• Kojic, M. (2002). Stress integration procedures for inelastic material models within finite element method. ASME Applied Mechanics Reviews 55:389-414.
• Lemaitre, J., Chaboche J.L. (1978). Aspect phénoménologique de la rupture par endommagement. Journal de Mécanique Appliquée 2:317-365.
• Lemaitre, J. (1985). A continuous damage mechanics model for ductile fracture. ASME Journal of Engineering Materials and Technology 107:83-89.
• Lemaitre, J. (1992). A course on damage mechanics, Springer, Berlin.
• Lemaitre, J., Desmorat, R., Sauzay, M. (2000). Anisotropic damage law of evolution. European Journal of Mechanics - A/Solids 19:187-208.
• Li, Y.F., Nemat-Nasser, S. (1993). An explicit integration scheme for finite-deformation plasticity in finite-element methods. Finite Elements in Analysis and Design 15:93-102.
• Li, J., Li, S., Xie, Z., Wang, W. (2015). Numerical simulation of incremental sheet forming based on GTN damage model. International Journal of Advanced Manufacturing Technology 81:2053-2065.
• Lumelskyy, D., Rojek, J., Pecherski, R., Grosman, F., Tkocz, M. (2012). Numerical simulation of formability tests of pre-deformed steel blanks. Archives of Civil and Mechanical Engineering 12:133-141.
• Mansouri, L.Z., Chalal, H. and Abed-Meraim, F. (2014). Ductility limit prediction using a GTN damage model coupled with localization bifurcation analysis. Mechanics of Materials 76:64-92.
• Marciniak, Z., Kuczynski, K. (1967). Limit strains in the process of stretch forming sheet metal. International Journal of Mechanical Sciences 9:609-620.
• Martínez-Donaire, A.J., García-Lomas, F.J., Vallellano, C. (2014). New approaches to detect the onset of localised necking in sheets under through-thickness strain gradients. Materials and Design 57:135-145.
• Mattiasson, K., Sigvant, M., Larson, M. (2006). Methods for forming limit prediction in ductile metal sheets. In: Santos, A.D., Barata da Rocha, A. (Eds.), Proc. IDDRG’06:1-9, Porto.
• Narasimhan, K., Wagoner, R.H. (1991). Finite element modeling simulation of in-plane forming limit diagrams of sheets containing finite defects. Metallurgical Transactions A 22:2255-2265.
• Peerlings, R.H.G., Geers, M.G.D., de Borst, R., Brekelmans, W.A.M. (2001). A critical comparison of nonlocal and gradient-enhanced softening continua. International Journal of Solids and Structures 38:7723-7746.
• Rice, J.R. (1976). The localization of plastic deformation. Theorical and Applied Mechanics. Koiter ed.:207-227.
• Rudnicki, J.W., Rice, J.R. (1975). Conditions for the localization of deformation in pressure sensitive dilatant materials. Journal of the Mechanics and Physics of Solids 23:371-394.
• Situ, Q., Jain, M.K., Bruhis, M. (2006). A suitable criterion for precise determination of incipient necking in sheet materials. Material Science Forum 519-521:111-116.
• Situ, Q., Jain, M.K., Bruhis, M. (2007). Further experimental verification of a proposed localized necking criterion. Proceedings of the Numerical Methods in Industrial Forming Processes - 9th Numiform 907-912.
• Situ, Q., Jain, M.K., Bruhis, M., Metzger, D.R. (2011). Determination of forming limit diagrams of sheet materials with a hybrid experimental-numerical approach. International Journal of Mechanical Sciences 53:707-719.
• Stören, S., Rice, J.R. (1975). Localized necking in thin sheets. Journal of the Mechanics and Physics of Solids 23:421-441.
• Swift, H.W. (1952). Plastic instability under plane stress. Journal of the Mechanics and Physics of Solids 1:1-18.
• Tvergaard, V. (1981). Influence of voids on shear band instabilities under plane strain conditions. International Journal of Fracture 17:389-407.
• Tvergaard, V. (1982a). Ductile fracture by cavity nucleation between larger voids. Journal of the Mechanics and Physics of Solids 30:265-286.
• Tvergaard, V. (1982b). Material failure by void coalescence in localized shear bands. International Journal of Solids and Structures 18:659-672.
• Tvergaard, V., Needleman, A. (1984). Analysis of the cup-cone fracture in a round tensile bar. Acta Metallurgica 32:157-169.
• Yamamoto, H. (1978). Conditions for shear localization in the ductile fracture of void-containing materials. International Journal of Fracture 14:347-365.
• Zhang, C., Leotoing, L., Zhao, G., Guines, D., Ragneau, E. (2011). A comparative study of different necking criteria for numerical and experimental prediction of FLCs. Journal of Materials Engineering and Performance 20:1036-1042.

Publication Dates

• Publication in this collection
  Oct 2017

History

• Received
  11 Feb 2017
• Reviewed
  27 June 2017
• Accepted
  07 Aug 2017