Muhammad Waqar Nasir
email: muhammad_waqar.nasir@ensam.eu

Hocine Chalal
email: hocine.chalal@ensam.euc

Farid Abed-Meraim
email: farid.abed-meraim@ensam.eu

Comparison between the Marciniak and Kuczyński imperfection approach and bifurcation theory in the prediction of localized necking for porous ductile materials

Keywords: Marciniak and Kuczyński imperfection approach, Gurson model, Bifurcation theory, Material imperfection, Forming limit diagram, Ductile damage

To prevent the occurrence of localized necking, the concept of forming limit diagram is often used, thus playing an important role in sheet metal forming processes. The aim of the present study is to develop a numerical tool for the theoretical prediction of forming limit diagrams, which would be a cost-effective procedure as compared to experimental measurements. The proposed numerical tool is based on the Marciniak and Kuczyński imperfection approach combined with the Gurson-Tvergaard-Needleman damage model, which is implemented into the MATLAB program within the framework of plane-stress conditions. Forming limit diagrams have been predicted by assuming both geometric (thickness) as well as material initial imperfections in the Marciniak and Kuczyński imperfection approach. These forming limit diagrams, for different sizes of geometric or material imperfections, are also compared with the forming limit diagram obtained by using the bifurcation theory. It is shown that the bifurcationbased forming limit diagram provides an upper bound as compared to the Marciniak and Kuczyński imperfection approach predictions. The results also reveal that, irrespective of the imperfection type considered in the Marciniak and Kuczyński imperfection approach, the corresponding forming limit diagram tends to that predicted by bifurcation theory when the size of initial imperfection tends to zero. Additionally, the predicted ductility limits are lowered as the magnitude of initial imperfection increases; however, the decrease in the ductility limits at balanced biaxial tension is more significant than for the other strain-path ratios. The results for the forming limit diagrams indicate that the predicted ductility limits are more sensitive to the initial imperfection in the thickness and the isotropic hardening coefficient as compared to the other types of material imperfections. Moreover, the initial imperfection in the critical porosity is the most influential one among the Gurson-Tvergaard-Needleman damage parameters.

Declarations

q q q void interaction material parameters Product between a second order tensor and a vector ( i j j = • • ) :

• •

Double contraction between two second-order tensors ( i j i j = • • ) Double contraction between fourth-and second-order tensors ( i j k l k l = • • )

Introduction

Sheet metal forming process has huge applications in various sectors, such as aerospace, automotive, architecture and food industry. The prediction of material formability plays an important role in order to optimize the sheet metal forming processes. Formability of the sheet metal may be limited due to the presence of internal defects, voids, surface roughness, or entirely due to the intrinsic material instability. Better understanding of these material instabilities and the impact of the internal defects ultimately increases the overall productivity of the manufacturing process. In this context, formability of the sheet metal is characterized with the help of a forming limit diagram (FLD), as proposed by Keeler & Backofen [START_REF] Keeler | Plastic instability and fracture in sheet stretched over rigid punches[END_REF] and improved by Goodwin [START_REF] Goodwin | Application of strain analysis to sheet metal forming problems in the press shop[END_REF]. FLD is defined as the plot of in-plane principal strains, namely major and minor strains, produced in the surface of the sheet at the onset of necking under biaxial stretching. Various experimental procedures are available in the literature for the determination of FLDs (see, e.g., [START_REF] Nakazima | Study on the formability of steel sheets[END_REF][START_REF] Hecker | A cup test for assessing stretchability[END_REF]). However, these procedures are very tedious and require advanced and expensive equipment. Therefore, various theoretical and numerical tools have been developed to determine FLDs. These tools require, in particular, the combination of advanced constitutive models and plastic instability criteria.

Considerable efforts have been made in the literature over the last few decades for the accurate modeling of the material response under complex loading conditions. To model the strain hardening phenomenon for simple applications, different isotropic hardening laws, such as Swift's, Ludwig's and Voce's laws, are commonly used. For complex strain paths under cyclic loading, combined isotropic and kinematic hardening laws, e.g., Armstrong-Frederick's and Prager's laws, have been developed in the literature to model some specific phenomena, such as Bauschinger effect, fatigue, ratcheting, etc. To account for the material anisotropy and the realistic modeling of deep drawing process, various anisotropic yield functions have been proposed, which are broadly classified into quadratic and non-quadratic yield functions. Hill [START_REF] Hill | A theory of the yielding and plastic flow of anisotropic metals[END_REF] generalized the isotropic von Mises yield function for anisotropic material with three orthogonal planes of symmetry. Karafillis and Boyce [START_REF] Karafillis | A general anisotropic yield criterion using bounds and a transformation weighting tensor[END_REF] extended their isotropic yield function for anisotropic material by calculating the principal values of transformed deviatoric stresses. As an extension to Karafillis and Boyce [START_REF] Karafillis | A general anisotropic yield criterion using bounds and a transformation weighting tensor[END_REF] yield function, Bron and Besson [START_REF] Bron | A yield function for anisotropic materials application to aluminum alloys[END_REF] proposed a phenomenological non-quadratic yield function, where the anisotropy is represented by 12 parameters in the form of two fourth-order symmetric tensors. Barlat and co-workers came up with several non-quadratic anisotropic yield functions, which involve one [START_REF] Barlat | A six-component yield function for anisotropic materials[END_REF], two [START_REF] Barlat | Linear transformation-based anisotropic yield functions[END_REF] and three [START_REF] Aretz | New convex yield functions for orthotropic metal plasticity[END_REF] linear transformations to the Cauchy stress tensor. The transformation matrices in these yield functions are composed of anisotropy coefficients. Banabic et al. [START_REF] Banabic | A new yield criterion for orthotropic sheet metals under plane-stress conditions[END_REF] extended the Hershey [START_REF] Hershey | The plasticity of an isotropic aggregate of anisotropic face centered cubic crystals[END_REF] isotropic yield function by adding anisotropy weight coefficients. Their yield function was subsequently improved by Banabic and co-workers, and eventually Banabic et al. [START_REF] Banabic | An improved analytical description of orthotropy in metallic sheets[END_REF] presented a yield function with eight material anisotropy parameters. This yield function is commonly known as BBC-2003 yield criterion. For anisotropic material modeling of facecentered polycrystalline metallic materials, Cazacu [START_REF] Cazacu | New yield criteria for isotropic and textured metallic materials[END_REF] extended the isotropic yield criterion using generalized invariants of the Cauchy stress deviator. Her orthotropic yield criterion is general and applicable to three-dimensional stress states. Cazacu [START_REF] Cazacu | New mathematical results and explicit expressions in terms of the stress components of Barlat et al. (1991) orthotropic yield criterion[END_REF] extended the yield functions proposed by Barlat et al. [START_REF] Barlat | A six-component yield function for anisotropic materials[END_REF] for body-centered cubic (BCC) and face-centered cubic (FCC) materials. Her yield functions can be expressed by a simple polynomial in terms of the Cauchy stress components. In addition, non-associative plasticity, in which the yield and potential functions are represented by separate expressions, can also be considered in such advanced constitutive equations (see, e.g., [START_REF] Shen | Forming limit prediction by the Marciniak-Kuczynski model coupled with the evolving non-associated Hill48 plasticity model[END_REF][START_REF] Tang | Study on non-associated plasticity with various forward Euler stress integration algorithms and its prediction of earing in cylindrical cup drawing[END_REF]).

Experimental observations have shown that ductile materials undergo a softening phenomenon at large plastic strains, in which the material gradually loses its load carrying capacity prior to final fracture. Moreover, Rudnicki and Rice [START_REF] Rudnicki | Conditions for the localization of deformation in pressuresensitive dilatant materials[END_REF] and Rice [START_REF] Rice | Localization of plastic deformation[END_REF] have shown that, within the framework of phenomenological constitutive models with smooth yield surface and associative plastic flow rule, the prediction of strain localization requires strongly negative hardening moduli. To model this softening phenomenon, elastic-plastic constitutive equations are generally coupled with ductile damage models. In this regard, two ductile damage theories, i.e., Continuum Damage Mechanics (CDM) and micromechanics of porous media, have been developed and extended over the last four decades. In the CDM approach, cracks are assumed on the surface of a representative volume element (RVE), and the surface density of cracks is treated as damage variable (see, e.g., [START_REF] Lemaitre | A Continuous Damage Mechanics Model for Ductile Fracture[END_REF][START_REF] Bouchard | An enhanced Lemaitre model formulation for materials processing damage computation[END_REF]). Brünig [START_REF] Brünig | Numerical analysis and elastic-plastic deformation behavior of anisotropically damaged solids[END_REF] and Maire and Chaboche [START_REF] Maire | A new formulation of continuum damage mechanics (CDM) for composite materials[END_REF] extended the CDM approach for anisotropic and composite materials, respectively.

During plastic deformations, the damage variable evolves and brings about the softening regime in the material response. The damage variable may be represented by a scalar or a tensor quantity, depending on the damage anisotropy and the softening response. In the second class of damage theories, i.e., micromechanics of porous media, existence of voids inside the RVE is assumed, as a representation of porous materials. The damage variable in this case is the ratio of volume of voids to the volume of RVE. Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF] first developed a ductile damage model based on the micromechanics approach, where a spherical void is assumed inside a spherical RVE. In this model, the damage variable evolves due to the void growth with the plastic deformation, and the effect of hydrostatic stresses is incorporated into the model. In addition to the void growth in the original Gurson model, several extensions have been subsequently introduced, such as void nucleation, void coalescence and void interaction parameters (see, e.g., Chu and Needleman [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF], Tvergaard [START_REF] Tvergaard | Influence of voids on shear band instabilities under plane strain conditions[END_REF][START_REF] Tvergaard | Material failure by void coalescence in localized shear bands[END_REF][START_REF] Tvergaard | On localization in ductile materials containing spherical voids[END_REF], Tvergaard and Needleman [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]). The resulting Gurson-Tvergaard-Needleman (GTN) ductile damage model has been further extended in various domains. For instance, Gologanu et al. [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric prolate ellipsoidal cavities[END_REF][START_REF] Gologanu | Approximate models for ductile metals containing nonspherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF] extended the GTN damage model for oblate and prolate void shapes, while Madou and Leblond [START_REF] Madou | A Gurson-type criterion for porous ductile solids containing arbitrary ellipsoidal voids-I: Limit-analysis of some representative cell[END_REF][START_REF] Madou | A Gurson-type criterion for porous ductile solids containing arbitrary ellipsoidal voids-II: Determination of yield criterion parameters[END_REF] extended it for general (nonspheroidal) voids. Variants of the Gurson model, specifically applicable for low stress triaxialities, have also been developed (see, e.g., Xue [START_REF] Xue | Constitutive modeling of void shearing effect in ductile fracture of porous materials[END_REF], Nahshon and Hutchinson [START_REF] Nahshon | Modification of the Gurson model for shear failure[END_REF], Zhou et al. [START_REF] Zhou | On the extension of the Gurson-type porous plasticity models for prediction of ductile fracture under sheardominated conditions[END_REF] and Daehli et al. [START_REF] Daehli | A Lode-dependent Gurson model motivated by unit cell analyses[END_REF]). Void size effects [START_REF] Dormieux | An extension of Gurson model incorporating interface stresses effects[END_REF], plastic anisotropy [START_REF] Morin | A new technique for finite element limit-analysis of Hill materials, with an application to the assessment of criteria for anisotropic plastic porous solids[END_REF], advanced isotropic and kinematic hardening [40] and coalescence criteria [START_REF] Torki | Theoretical and numerical analysis of void coalescence in porous ductile solids under arbitrary loadings[END_REF][START_REF] Gallican | Anisotropic coalescence criterion for nanoporous materials[END_REF] have also been incorporated into extended versions of the GTN damage model.

In order to predict the onset of plastic instability, various theoretical instability criteria have been proposed in the literature. These instability criteria are mainly classified into three categories. The first category is based on the maximum force concept. Within this approach, Considère [START_REF] Considère | Mémoire sur l'emploi du fer et de l'acier dans les constructions[END_REF] first developed a diffuse necking criterion for uniaxial tension, which has been extended afterwards by Swift [START_REF] Swift | Plastic instability under plane stress[END_REF] for biaxial tension. Later, Hill [START_REF] Hill | On discontinuous plastic states, with special reference to localized necking in thin sheets[END_REF] presented an instability criterion for the prediction of localized necking, commonly known as Hill's zero extension theory. It is noteworthy that the Hill [START_REF] Hill | On discontinuous plastic states, with special reference to localized necking in thin sheets[END_REF] criterion is limited to the range of negative minor strains. Therefore, in practice, Hill [START_REF] Hill | On discontinuous plastic states, with special reference to localized necking in thin sheets[END_REF] criterion is sometimes combined with Swift [START_REF] Swift | Plastic instability under plane stress[END_REF] diffuse necking criterion for the prediction of a complete FLD. Note also that Considère's and Swift's criteria have been extended for the prediction of localized necking on both sides of the FLD by Hora et al. [START_REF] Hora | A prediction method for ductile sheet metal failure in FE-simulation[END_REF]. The second category of plastic instability criteria is based on the bifurcation theory and includes criteria for both diffuse and localized necking. Drucker [START_REF] Drucker | On uniqueness in the theory of plasticity[END_REF] and later Hill [START_REF] Hill | A general theory of uniqueness and stability in elastic-plastic solids[END_REF] proposed the general bifurcation theory for the prediction of diffuse necking. The theory requires the condition of positive definiteness of the tangent modulus for non-bifurcation state. Valanis [START_REF] Valanis | Banding and stability in plastic materials[END_REF] introduced the limit-point bifurcation theory, which predicts diffuse necking corresponding to stationary value of the first Piola-Kirchhoff stress. This criterion is less conservative than the general bifurcation one. In order to predict plastic instability in the form of a planar localization band, Rudnicki and Rice [START_REF] Rudnicki | Conditions for the localization of deformation in pressuresensitive dilatant materials[END_REF], Stören and Rice [START_REF] Stören | Localized necking in thin sheets[END_REF] and Rice [START_REF] Rice | Localization of plastic deformation[END_REF] developed the loss of ellipticity condition for boundary value problems. The criterion predicts localized necking when the acoustic tensor becomes singular. Bigoni and Hueckel [START_REF] Bigoni | Uniqueness and localization-I. Associative and nonassociative elastoplasticity[END_REF] and Neilsen and Schreyer [START_REF] Neilsen | Bifurcations in elastic-plastic materials[END_REF] presented a loss of strong ellipticity criterion for localized necking, which is more conservative than the loss of ellipticity criterion. Among all the above-discussed bifurcation criteria, the general bifurcation condition is the most conservative and provides a lower bound, whereas the loss of ellipticity criterion provides an upper bound. It should be pointed out that the plastic instability criteria based on the bifurcation theory are only applicable for strain-rate independent constitutive equations. [START_REF] Marciniak | Limit strains in the processes of stretch-forming sheet metal[END_REF] laid the foundation of the third category of instability criteria, which is based on the assumption of initial heterogeneity inside the sheet metal. In the original Marciniak and Kuczyński (M-K) model, the initial heterogeneity was assumed in the form of a narrow band of reduced thickness (i.e., geometric imperfection). Later, Marciniak [START_REF] Marciniak | Analysis of necking preceding fracture of sheet metal under tension[END_REF] introduced the band inclination angle, which enables the M-K model to predict realistic limit strains on the left-hand side of the FLD. Moreover, Hutchinson and Neale [START_REF] Hutchinson | Sheet necking-II. Time-independent behavior[END_REF] extended the M-K model by allowing the imperfection band rotation during the plastic deformation. The M-K model has also been extended for visco-plastic materials by Marciniak et al. [START_REF] Marciniak | Influence of the plastic properties of a material on the forming limit diagram for sheet metal in tension[END_REF]. Yamamoto [START_REF] Yamamoto | Conditions for shear localization in the ductile fracture of voidcontaining materials[END_REF] introduced the concept of initial heterogeneity of the mechanical properties (i.e., material imperfection) within a narrow imperfection band. The M-K model has been successfully coupled with various constitutive models, including ductile damage models and has been widely used for the prediction of FLDs. An extensive review on the M-K model and its coupling with various forms of equivalent stresses, polycrystalline yield surfaces, and ductile damage models has recently been presented by Banabic [START_REF] Banabic | A review on recent developments of Marciniak-Kuczynski model[END_REF] and Banabic et al. [START_REF] Banabic | Developments of the Marciniak-Kuczynski model for sheet metal formability: A review[END_REF].

Marciniak and Kuczyński

The coupling of ductile damage models with the M-K imperfection approach could improve the overall predictions of FLDs with respect to experiments [START_REF] Hu | Influence of damage and texture evolution on limit strain in biaxially stretched aluminum alloy sheets[END_REF]. Using an elastic-plastic-damage model combined with the M-K approach, initial imperfection could be introduced in the form of an imperfection band with reduced thickness; increased value of damage; different set of hardening parameters, or any combination of imperfections that results in reduced hardening slope within the imperfection band as compared to the homogeneous zone. In the pioneering work of Yamamoto [START_REF] Yamamoto | Conditions for shear localization in the ductile fracture of voidcontaining materials[END_REF], the original Gurson model was coupled with the M-K model, in which the initial imperfection was introduced in the localization band as a slightly higher value of initial porosity as compared to the homogeneous zone. Ductility limits for two strain paths, i.e. axisymmetric uniaxial tension and plane-strain tension, were analyzed. Needleman and Triantafyllidis [START_REF] Needleman | Void growth and local necking in biaxially stretched sheets[END_REF] extended the approach of Yamamoto [START_REF] Yamamoto | Conditions for shear localization in the ductile fracture of voidcontaining materials[END_REF] by considering more strain paths on the right-hand side of the FLD and allowing for the change in thickness of the imperfection band due to plastic deformation. This approach has been further extended by incorporating the effect of void nucleation (see Chu and Needleman [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF]). A similar approach has been followed by Saje et al. [START_REF] Saje | Void nucleation effects on shear localization in porous plastic solids[END_REF] and Pan et al. [START_REF] Pan | Localization of deformation in rate sensitive porous plastic solids[END_REF] for axisymmetric tension and planestrain tension by incorporating the strain-rate sensitivity and the effect of stress triaxiality during necking. Numerous applications of the approach combining the M-K method with damage models for the prediction of ductility limits have been presented in the literature over the last two decades. For instance, Huang et al. [START_REF] Huang | Failure prediction in anisotropic sheet metals under forming operations with consideration of rotating principal stretch directions[END_REF], Simha et al. [START_REF] Simha | Computational analysis of stress-based forming limit curves[END_REF], Liu et al. [START_REF] Liu | Numerical investigations on the influence of superimposed double-sided pressure on the formability of biaxially stretched AA6111-T4 sheet metal[END_REF], Hosseini et al. [START_REF] Hosseini | Theoretical FLD Prediction Based on MK Model using Gurson's Plastic Potential Function for Steel Sheets[END_REF] and Morin et al. [START_REF] Morin | On the description of ductile fracture in metals by the strain localization theory[END_REF] considered the material imperfection in the form of initial porosity in Gursonbased damage models. Additionally, combination of material and geometric imperfections has also been considered in some investigations. For instance, Regab et al. [START_REF] Ragab | Forming limit diagrams for kinematically hardened voided sheet metals[END_REF][START_REF] Ragab | Effect of void growth on predicting forming limit strains for planar isotropic sheet metals[END_REF] coupled the M-K model with Green's yield function for planar isotropic porous material for the prediction of FLDs. Zadpoor et al. [START_REF] Zadpoor | Formability prediction of high strength aluminum sheets[END_REF] and Hu et al. [START_REF] Hu | A thermal forming limit prediction method considering material damage for 22MnB5 sheet[END_REF] also predicted FLDs by combining the M-K model with a modified version of the Gurson model. Son and Kim [START_REF] Son | Prediction of forming limits for anisotropic sheets containing prolate ellipsoidal voids[END_REF] investigated the effect of void shape on the FLDs by using the M-K model combined with the Gologanu-Leblond-Devaux (GLD) damage model for prolate ellipsoidal voids. Chien et al. [START_REF] Chien | A combined necking and shear localization analysis for aluminum sheets under biaxial stretching conditions[END_REF] implemented the combined material and geometric imperfections for the accurate modeling of shear and necking modes.

In the present contribution, FLDs are predicted by using the M-K imperfection approach combined with the GTN damage model. The effects of both material and geometric imperfections are investigated in the prediction of FLDs. An analogy is established between the M-K imperfection approach and the loss of ellipticity (Rice's bifurcation) criterion for the prediction of localized necking in sheet metals. Indeed, it has been observed by Tvergaard [75] that the M-K imperfection approach reduces to the Stören and Rice [START_REF] Stören | Localized necking in thin sheets[END_REF] bifurcation criterion when the initial imperfection factor is set to zero. Similar analogy has also been pointed out in the literature (see, e.g., [START_REF] Yoshida | Comparison of bifurcation and imperfection analyses of localized necking in rate-independent polycrystalline sheets[END_REF][START_REF] Bettaieb | Investigation of localized necking in substratesupported metal layers: Comparison of bifurcation and imperfection analyses[END_REF]), where J2 deformation theory of plasticity was combined with the M-K method and the Rice bifurcation theory. Abed-Meraim et al. [START_REF] Abed-Meraim | Investigation and comparative analysis of plastic instability criteria: application to forming limit diagrams[END_REF] have also compared the FLDs obtained by using the M-K method and the Rice bifurcation theory, which have been coupled with the CDM. Note that the Rice bifurcation theory assumes homogeneous prelocalization state, according to the pioneering work of Rudnicki and Rice [START_REF] Rudnicki | Conditions for the localization of deformation in pressuresensitive dilatant materials[END_REF]. In other words, structural (geometric) instabilities, such as diffuse necking, and the associated geometric effects, are purposely excluded in the Rice bifurcation theory. The main motivation is to provide an 'intrinsic' strain localization criterion, inherent to the 'material' alone, with no interference with structural (geometric) effects, which eventually leads to higher levels of forming limits especially in the right-hand side of the FLD (see, e.g., [START_REF] Abed-Meraim | Investigation and comparative analysis of plastic instability criteria: application to forming limit diagrams[END_REF]). Consequently, the limit strains predicted with the Rice bifurcation theory are expected to provide an upper bound to those that would be determined experimentally. Therefore, the Rice bifurcation theory could be used as a "qualitative tool" for the prediction of FLDs with respect to experimental data. It can also be used, in the design stage of new materials with enhanced ductility, as an efficient tool for the comparison of formability of several sheets belonging to the same grade of material. However, for the "quantitative comparison" of the predicted limit strains with respect to experimental measurements, it is necessary to reproduce the same operating conditions; namely, the consideration of diffuse necking prior to the occurrence of localized necking, the close correspondence of the boundary conditions between the simulation and the experiments, defects naturally present in the material, such as surface roughness, inhomogeneous distribution of voids inside the material, non-uniform material thickness, etc. Since all of these mechanisms, defects, inhomogeneities and operating conditions are not easy to be incorporated simultaneously during the simulations, therefore the M-K imperfection approach provides a way forward to overcome these difficulties. The concept of initial imperfection factor in the M-K model is considered as a symbolic representation, which encompasses all defects and destabilizing mechanisms typically present inside the sheet metal during forming process. For these reasons, the M-K model provides a very promising numerical tool for the accurate prediction of FLD for any material. With the proposed numerical tool, instead of experimentally determining the FLD for the complete range of loading paths, the limit strains can be measured experimentally for a single strain-path ratio, typically uniaxial tension. Afterwards, the imperfection factor is numerically calibrated to obtain the same limit strains as measured experimentally. In such a way, the proposed numerical tool provides a cost-effective procedure for the determination of FLDs for any given materials. In other words, the calibrated value of the imperfection factor can be regarded as an additional parameter, which could be determined experimentally within an identification procedure including the intrinsic material parameters.

The present research paper is organized as follows. First, the constitutive equations related to the GTN damage model are briefly recalled in Section 2. Then, the M-K imperfection approach and the Rice bifurcation theory are presented in Section 3, followed by the numerical implementation and validations of the developed numerical tool in Section 4. Numerical results in terms of predictions of FLDs are discussed in Section 5. Finally, some concluding remarks are presented in Section 6.

Constitutive equations 2.1 GTN yield function

In this section, elastic-plastic constitutive equations are coupled with porous plasticity modeling for non-viscous material within the framework of associative plasticity. Specifically, Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF] ductile damage model, which is based on the micromechanics of a spherical void present inside a spherical RVE, is considered. In this model, the damage variable, i.e. porosity f , is defined as the ratio of volume of voids to the volume of RVE. In the original Gurson model, the evolution of porosity is attributed to only void growth. Incorporating the effects of void nucleation, void coalescence and void interactions in the original Gurson model, the extended version, i.e. GTN damage model, is expressed by the following yield condition:

2 2 eq * * m Y Σ 3 Σ Φ 2 cosh 1 0 σ 2σ                       2 1 3 q q f q f = + -- ≤ (1)
where ( )

' ' eq Σ 3 : / 2 Σ Σ Σ = represents the macroscopic von Mises equivalent stress, m Σ tr Σ = ( ) / 3 is the hydrostatic stress, ' m Σ Σ Σ I = -
is the deviatoric part of the macroscopic Cauchy stress tensor Σ , and I is the second-order identity tensor. As to σ , it represents the tensile yield stress of the fully dense matrix material, while 1 q , 2 q and 3 q are the void interaction parameters, which have been introduced by Tvergaard [START_REF] Tvergaard | Influence of voids on shear band instabilities under plane strain conditions[END_REF]. Tvergaard and Needleman [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] modeled the void coalescence regime by introducing the concept of effective porosity ( ) * f f . According to their phenomenological coalescence model, when the value of the porosity reaches the critical porosity cr f , void coalescence takes place, which results in the accelerated decay of load carrying capacity of the material. The relationship for the effective porosity is given as: * ( )

cr GTN cr f f δ f f = + - (2)
where GTN δ is the accelerating factor, which is equal to 1 before the onset of coalescence (i.e., cr f f ≤ ), and greater than 1 once the coalescence regime is activated (i.e.,

cr f f > ). Both cr f
and GTN δ are considered as material parameters. Within the framework of associative plasticity, the macroscopic plastic strain rate is assumed normal to the yield function, and is defined by the following classical normal flow rule:

p Y Φ λ λ = Σ D M Σ ∂ = ∂ ɺ ɺ (3)
where p D represents the macroscopic plastic strain rate tensor. The plastic multiplier ɺ λ and Σ M represent the magnitude and direction of the macroscopic plastic strain rate, respectively. The fully dense matrix material is assumed to follow an isotropic strain hardening law defined by the Swift hardening model, i.e., ( )

p p 0 σ( ) n ε K ε ε = + (4)
where K , 0 ε , and n are the hardening parameters and p ε is the equivalent plastic strain of the fully dense matrix material.

Evolution equations for the internal variables

The yield stress defined in Eq. ( 4) can be expressed in a rate form as follows:

p p p p p σ σ( ) ( ) ∂ ∂ ɺ ɺ ɺ ε ε h ε ε ε = = (5)
where p ( ) h ε is the plastic hardening slope of the fully dense matrix material. By adopting the equivalence principle in the plastic work rate, the relationship between the equivalent plastic strain rate p ɺ ε and the macroscopic plastic strain rate tensor p D can be derived as follows:

p p (1 )σ : f ε - = Σ D ɺ (6)
The evolution of the porosity ɺ f can be partitioned as the contribution of void growth ɺ G f and void nucleation ɺ

N f G N f f f = + ɺ ɺ ɺ (7)
Tvergaard [START_REF] Tvergaard | Effect of yield surface curvature and void nucleation on plastic flow localization[END_REF] developed the expression for ɺ G f by using the incompressibility principle of the dense matrix material. The evolution equation for ɺ G f mainly depends on the stress triaxiality and is given by: p

:

G f f = -D I ɺ (8)
In the present contribution, strain-controlled void nucleation is considered, which has been modeled by Chu and Needleman [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF] using the following normal distribution relationship:

p N N f A ε = ɺ ɺ (9) with 2 p 1 exp 2 2π N N N N ψ ε ε A s s - = -                            (10)
where ψ , N ε and N s are the material nucleation parameters.

Determination of the plastic multiplier and the elastic-plastic tangent modulus

The expressions of the plastic multiplier ɺ λ and the elastic-plastic tangent modulus ep C for the GTN damage model are derived in this section. For this purpose, the GTN yield function Y Φ and the plastic multiplier ɺ λ are first expressed in the Kuhn-Tucker form

Y Y Φ 0, 0, Φ 0 λ λ ɺ ɺ ≤ ≥ = (11)
The above relations reveal that when Y Φ 0 < , the plastic flow is zero (i.

* Y Y Y Y * Φ Φ Φ Φ : + σ + 0 σ f f = = Σ Σ       ∂ ∂ ∂                         ∂ ∂ ∂   ɺ ɺ ɺ ɺ (12)
where the different partial derivatives involved in the above equation can be evaluated as:

* ' * Y m 2 * Y m * 2 * eq Y m m σ 3 2 Φ 3 3 Σ + sinh σ σ 2σ Φ 3 Σ M 2 cosh 2 2σ 2Σ Φ 3 Σ 3 Σ M sinh σ σ σ 2σ 1 2 2 2 1 3 f 1 2 2 q q f q q q q f f q q f q = = = = - = = - - Σ Σ M I Σ   ∂         ∂   ∂         ∂   ∂         ∂ (13)
In Eq. ( 12), the Cauchy stress rate tensor Σ ɺ is given by the following hypoelastic law in the corotational material frame: C represents the fourth-order tensor of the elasticity constants, and D is the macroscopic strain rate tensor. By substituting Eqs. ( 5)- [START_REF] Aretz | New convex yield functions for orthotropic metal plasticity[END_REF] and Eq. ( 14) into Eq. ( 12), the final expression of the plastic multiplier ɺ λ is derived as follows:

e : :

H λ λ = Σ M C D ɺ (15) where, ( ) ( ) ( )
* σ e M : :

H : : M + (1 ) : (1 )σ (1 )σ N λ GTN f h A δ f f f = - - - - - Σ Σ Σ Σ Σ Σ M Σ M M C M M I         (16)
Finally, the elastic-plastic tangent modulus ep C is derived by substituting the expression of the plastic multiplier ɺ λ (Eq. ( 15)) into the hypoelastic law (Eq. ( 14)):

( ) ( ) e e ep e
: :

H λ = - Σ Σ C M M C C C ⊗ (17)
3 Localized necking criteria

M-K imperfection approach

The M-K imperfection approach has been developed for the prediction of localized necking in thin sheet metals. In this approach, an initial imperfection is assumed within the sheet in the form of a narrow band. The initial imperfection may be considered as a narrow band of reduced thickness (geometric imperfection), an imperfection band with slightly degraded mechanical properties (material imperfection), or the combination of both geometric and material imperfections. In the original M-K imperfection approach (see Marciniak and Kuczyński [START_REF] Marciniak | Limit strains in the processes of stretch-forming sheet metal[END_REF]), the imperfection band has been assumed perpendicular to the major principal stress axis, which ultimately leads to unrealistically high values of forming limit strains on the left-hand side of the FLD (negative minor strains). Moreover, experimental evidences have shown that the localization band is not perpendicular to the major principal stress axis for the strain paths that cover the left-hand side of the FLD. Therefore, Hutchinson and Neale [START_REF] Hutchinson | Sheet necking-II. Time-independent behavior[END_REF] modified the M-K model, in which the normal to the imperfection band is assumed at an angle 0 θ from the major principal stress axis. The initial orientation 0 θ of the localization band is searched in the plane of the sheet metal in order to minimize the limit strains. Since the imperfection band is considered as weaker than the homogeneous zone, therefore during in-plane biaxial stretching, the strains in the imperfection band always remain larger than those in the homogeneous zone. At a critical loading point, the plastic strains in the imperfection band evolve drastically as compared to those in the homogeneous zone. Localization is assumed to take place when the ratio of equivalent plastic strain rate inside the imperfection band to the homogeneous zone, i.e., p p / ɺ ɺ B H ε ε , exceeds a critical value. The latter is commonly set equal to 10 in the literature. Fig. 1 depicts the schematic illustration of the imperfection band present inside the homogeneous zone. Both the material and geometric imperfections are illustrated in Fig. 1. The material imperfection is represented in this figure by an increased void volume fraction inside the imperfection band, while the geometric imperfection is represented by a reduced thickness of the imperfection band. The initial geometric imperfection factor is defined as:

0 z = B i H i h h ( 18
)
where the superscripts B and H indicate the quantities inside and outside the imperfection band, respectively. The key components of M-K imperfection approach are as follows:

• Definition of the constitutive equations for both zones.

• Kinematic compatibility of strains between the imperfection band and the homogeneous zone.

• The equilibrium of forces across the imperfection band.

Fig. 1 Schematic representation of the initial imperfection band in the M-K imperfection approach

In the present contribution, the constitutive equations defined in Section 2 are assumed to be applicable for both zones H and B. Following the condition of kinematic compatibility, the velocity gradients in the imperfection band G B and in the homogeneous zone G H are related as:

B H = + G G c n ⊗ ɺ ( 19
)
where c ɺ is the jump amplitude vector, and n is the unit vector normal to the imperfection band.

The homogeneous zone is subjected to various linear strain paths under in-plane biaxial stretching for the prediction of FLD. The strain-path ratio is defined as:

22 22 11 11
; 0.5

1 ɺ ɺ H H H H ε G ρ ρ ε G = = - ≤ ≤ (20)
where the particular values 0.5 ρ = - and 1 ρ = represent the uniaxial tension and balanced biaxial tension, respectively. Therefore, the velocity gradient tensor in the homogeneous zone is given as:

              G H H H H ε ρ ε G ɺ ɺ = (21)
By substituting Eq. ( 21) into Eq. ( 19), the velocity gradient tensor inside the imperfection band is calculated as:

11 1 1 1 2 2 1 11 2 2 33 c n c n 0 c n c n 0 0 0               G H B H B ε ρ ε G ɺ ɺ ɺ ɺ ɺ ɺ = + + (22)
The components 33 ( )

z Σ Σ n 0 - ⋅ = B H (23)
where "z" is the value of the geometric imperfection factor in the deformed configuration, which can be calculated as:

0 33 33 33 exp( ) z z exp( ) exp( ) B B B B H i H H H i h ε h ε ε h h ε = = = - (24)
where 0 z is the initial value of the geometric imperfection factor.

Note that the imperfection band rotates during the applied loading. The evolution of the band inclination angle θ is given according to Nanson's formula:

0 11 22
tan( tan( ) exp( )

H H ε ε θ) = θ - (25)
where 0 θ is the initial orientation of the imperfection band with respect to the global axis.

Loss of ellipticity criterion

In this subsection, the Rice localized necking criterion (see Rice [START_REF] Rice | Localization of plastic deformation[END_REF] and Stören and Rice [START_REF] Stören | Localized necking in thin sheets[END_REF]) is derived and compared with the M-K imperfection approach. This strain localization criterion is based on the bifurcation theory and corresponds to the loss of ellipticity of the partial differential equations that govern the associated boundary value problem. By following an updated Lagrangian approach, the equilibrium equation stated in Eq. ( 23) can be rewritten in the following equivalent rate form:

( )

z n N N 0 ⋅ - = B H ɺ ɺ ( 26
)
where N ɺ is the nominal stress rate tensor, which is related to the velocity gradient tensor G as follows:

:

N L G ɺ = ( 27
)
where L is a fourth-order tangent modulus whose expression can be derived in terms of the elastic-plastic tangent modulus ep C as follows (see Mansouri et al. [START_REF] Mansouri | Ductility limit prediction using a GTN damage model coupled with localization bifurcation analysis[END_REF]):

ep 1 2 3 = + -- L C C C C (28)
In the above equation, the analytical elastic-plastic tangent modulus ep C is defined in Eq. ( 17), while 1 C , 2 C and 3 C are composed of convective stress components, which are induced by the large strain framework. Their expressions only depend on Cauchy stress components, and are given by the following expressions:

( ) ( ) 1i jkl ij kl 2ijkl jk il jl ik 3ijkl ik jl il jk C Σ 1 C Σ Σ 2 1 C Σ Σ 2 δ δ δ δ δ = = + = - (29)
By substituting Eq. ( 27) into Eq. ( 26), the equilibrium equations can be written in terms of the velocity gradients:

( ) ( ) ( ) z : : n L G L G 0 ⋅ - = B B H H ( 30
)
Finally, by considering the compatibility condition given by Eq. ( 19), the above equation can be expanded as:

( ) ( ) z : z   ⋅ ⋅ ⋅ ⋅     n L L G n L n c 0 ɺ - + = B H H B ( 31
)
In the case when no material imperfection is introduced, and in the absence of geometric imperfection (i.e.,

0 z z / 1 = = = B H i i h h
), meaning that the material properties within the band remain identical to those in the homogeneous zone (i.e., L L L B H = = ), Eq. ( 31) simplifies to:

( ) ⋅ ⋅ ⋅ n L n c 0 ɺ = (32)
The non-trivial solution for Eq. ( 32) corresponds to the Rice loss of ellipticity criterion. More specifically, the critical condition for bifurcation, which corresponds to the occurrence of localized necking, can be expressed as the singularity of the acoustic tensor

⋅ ⋅ Q n L n = : ( ) ( ) det det 0 ⋅ ⋅ Q n L n = = (33)
In practice, localization bifurcation is predicted numerically during the process of loading when the minimum value of the determinant of the acoustic tensor, for all possible orientations of the normal n to the localization band, becomes non-positive.

It is worth noting that when no geometric and material imperfections are considered, the M-K imperfection approach becomes equivalent to the Rice bifurcation criterion. Note also that, in order to predict localized necking on the right-hand side of FLD using the Rice bifurcation theory, negative hardening moduli are required, which can be achieved by coupling the elasticplastic constitutive equations with ductile damage.

It must be noted that the constitutive equations have been developed in the previous sections by adopting a fully three-dimensional approach. However, the plane-stress conditions are required for the implementation of the M-K imperfection and Rice bifurcation criteria, since the intended application of the present contribution is the prediction of ductility limits of thin sheet metals. Accordingly, the three-dimensional tangent modulus L given by Eq. ( 28) is reformulated for the plane-stress conditions as follows:

33 PS 33 3333 L L L L , with , , , 1, 2 L γδ αβγδ αβγδ αβ α β γ δ = - = (34)
where PS L represents the tangent modulus for plane-stress conditions. Furthermore, the localization band is sought in the plane of the sheet, defined by an in-plane band orientation angle θ (see Fig. 1).

4

Numerical implementation and validations

Numerical implementation of the M-K approach with the GTN model

The GTN damage model discussed in Section 2 is coupled with the M-K imperfection approach, and the resulting numerical tool is implemented into MATLAB user subroutines for the prediction of FLDs. The homogeneous zone of the sheet metal, as illustrated in Fig. 1, is subjected to proportional in-plane biaxial strains, with a strain-path ratio ρ . For each strain-path ratio considered between -0.5 and 1, the initial band orientation angle 0 θ is varied from 0 to 90 . For the given value set of the initial band angle 0 θ and the strain-path ratio ρ , the stress- strain as well as the internal variables are updated for each loading increment, according to an explicit time integration scheme. The simulation is stopped when the ratio of equivalent plastic strain rate inside the imperfection band to the homogeneous zone (i.e. , are the critical strains for the selected value of 0 θ . This procedure is repeated for all possible values of the initial band angle 0 θ , and the smallest values of the critical strains are considered as the forming limit strains for the considered strain-path ratio ρ .

Algorithm and time integration scheme for the M-K approach with the GTN model

The numerical simulations performed in the present contribution are considered to be strain driven. In this work, the equivalent plastic strain rate inside the imperfection band is considered as a prescribed quantity and, throughout the loading, its value is set equal to 1:

p 1 0 ɺ B ε -= (35)
All of the variables at the beginning and at the end of the loading increment are designated with subscript "n" and "n+1", respectively. The evolution of the geometric imperfection factor and the imperfection band angle (see Eqs. ( 24)-( 25)) can be expressed in their incremental forms as follows:

1 33 33 z z exp(Δ ( )) B H n n t G G + = - (36)
and 1 11 tan( tan( )exp(Δ ( 1) )

H n n t ε ρ + θ ) = θ -ɺ (37)
The unit vector 1 n n+ normal to the imperfection band lies in the plane of the sheet, and its components can be updated as:

( ) ( ) T 1 1 1 cos sin 0 n n n n + + + = θ θ     (38)
Applying the plane-stress conditions and exploiting the fact that the third component of the normal vector 1 n n+ is equal to zero, the equilibrium equation stated in Eq. ( 23) can be expanded in the following form: For a given strain-path ratio ρ and an initial imperfection band angle 0 θ , the velocity gradients in the imperfection band and in the homogeneous zone only depend on three unknowns: 11 ɺ H ε , 1 c ɺ and 2 c ɺ (see Eqs. ( 21)-( 22)). By combining Eq. ( 35) with Eq. ( 39), a nonlinear system of three scalar equations is obtained:

( ) ( ) ( ) ( ) 1
( ) ( ) ( ) ( ) 1
( , c , c ) z Σ Σ n Σ Σ n 0 ( , c , c ) Σ Σ n Σ Σ n 0
( , c , c ) 1 0

+ + + + + + + + + + + + + + + + = - + - = = - + - = = -= H B H B H n n n n n n n n H B H B H n n n n n n n n H B ε z ε z z ε ε ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺ R R R (40)
The above system of non-linear equations is solved iteratively at each loading increment to determine the three unknowns 11 ɺ H ε , 1 c ɺ and 2 c ɺ . As to the internal variables involved in the GTN damage model (i.e., porosity f , yield stress σ , equivalent plastic strain p ε , and Cauchy stress tensor Σ ), they can be expressed by a general differential equation of the form:

( ) = ɺ y y g y [START_REF] Torki | Theoretical and numerical analysis of void coalescence in porous ductile solids under arbitrary loadings[END_REF] where the vector y encompasses all of the GTN model variables to be updated at the end of each loading increment. This vector y is updated using an explicit time integration scheme. The complete algorithm for the M-K approach coupled with the GTN model is outlined in Fig. 2.

Fig. 2 Flow diagram for the M-K-GTN algorithm

Numerical implementation of the Rice bifurcation criterion with the GTN model

In the present contribution, the above approach coupling the M-K method with the GTN model is compared with the one combining the Rice bifurcation criterion with the GTN model. To achieve this comparative study, the Rice bifurcation criterion as well as the GTN model are implemented into the finite element code ABAQUS/Standard. It is worth noting that, in contrast to the M-K approach, the Rice bifurcation theory requires a fully homogeneous strain state prior to the occurrence of plastic instability. In other words, only intrinsic strain localization, inherent to the material alone and with no interference with structural (geometric) effects, is considered in the Rice bifurcation approach. To this purpose, a single finite element with one integration point (specifically, C3D8R solid element in ABAQUS) is considered in the simulations (see Fig. 3). This finite element is subjected to various linear strain paths that are those typically applied to sheet metals under in-plane biaxial stretching. The constitutive equations for the GTN damage model are implemented into ABAQUS via a user-defined material (UMAT) subroutine. As to the Rice bifurcation criterion, it is implemented into ABAQUS via a user-defined output variables (UVARM) subroutine. For a given strain-path ratio, the acoustic tensor Q is evaluated for each loading increment for all of the possible orientations of the localization band (i.e., 0° to 90°). The occurrence of strain localization is predicted when the minimum value of the determinant of the acoustic tensor, over all the considered orientations of the localization band, becomes nonpositive.

Fig. 3 Boundary conditions for the single finite element (C3D8R) subjected to in-plane biaxial stretching

Numerical validations of the adopted algorithm

This subsection is devoted to the numerical validations of the numerical implementation of the constitutive equations into the finite element code ABAQUS and into the MATLAB software.

Since the same algorithm and integration scheme are used for the numerical implementation in ABAQUS and in MATLAB, only the numerical results obtained with ABAQUS UMAT subroutine are shown, for conciseness. Some preliminary simulations are first performed, and the results obtained from the present UMAT are compared with the corresponding results obtained from the built-in Gurson-based damage model available in ABAQUS/Standard. It is worth noting that the complete GTN damage model, with growth, nucleation and coalescence mechanisms, is not available in ABAQUS/Standard, which restricts the following validation to only growth and nucleation of voids. The elastic properties, Swift's isotropic hardening parameters and damage parameters used in the simulations are listed in Tables 1 and2. Note that an implicit time integration scheme is considered in the built-in ABAQUS model, where a large loading increment could be used without the expense of accuracy. By contrast, an explicit time integration scheme is adopted in the present work, which requires loading increments to be taken small. Therefore, in order to estimate the size of the loading increment that guarantees a reasonable accuracy when compared with an implicit scheme, Fig. 4 provides the true stressstrain curve under uniaxial tension obtained from the UMAT subroutine for different values of strain increment, along with the stress-strain curve obtained from the built-in ABAQUS model.

Fig. 4 shows that the strain increment must be maintained lower than -5

1×10 to achieve the required accuracy.

(Δ ×10 ) ε = 1
Fig. 4 Comparison of the uniaxial stress-strain curves obtained with the built-in ABAQUS model and the developed UMAT subroutine In order to further evaluate the accuracy of the developed algorithm under in-plane biaxial loading, Fig. 5 shows the true stress-strain curve for balanced biaxial tension, as obtained from the developed UMAT along with the reference solution obtained from the built-in ABAQUS model. This figure shows that the numerical results given by the developed UMAT are in excellent agreement with those obtained with the built-in ABAQUS model, which validates the numerical implementation of the present model.

It is worth noting that, in the case of the M-K imperfection approach, shear stresses and strains may develop inside the imperfection band, although the homogeneous zone is subjected to in-plane biaxial stretching. To validate the accuracy of the developed numerical tool with respect to the shear behavior, Fig. 6 depicts the shear stress-strain curve and the porosity evolution curve for pure shear loading. This figure reveals once again that the numerical results obtained with the developed UMAT are in excellent agreement with the reference results provided by the built-in ABAQUS model.

The final set of simulations are specifically intended to validate the loading and unloading response. For this purpose, the considered material is subjected to cyclic loading under pure shear, and the shear stress-strain response is plotted in Fig. 7. From this figure, it can be observed that the stress-strain curve given by the developed UMAT coincides with that provided by the built-in ABAQUS model.

Prediction of FLDs

In the present contribution, the elastic-plastic-damage (GTN) model, described in Section 2, is combined with two distinct classes of plastic instability criteria: the Rice bifurcation theory and the M-K initial imperfection approach, for the prediction of FLDs. The previously established theoretical link between both plastic instability criteria (see Section 3) is thoroughly investigated in this section through numerical simulations. First, the FLDs are predicted for a fictitious material, which is similar to an aluminum alloy, using the geometric imperfection in M-K approach. Then, by considering various types of material imperfections, a sensitivity study is conducted in order to determine the most influential parameters on strain localization. The FLDs predicted with the M-K approach are systematically compared with the FLD obtained using the Rice bifurcation criterion. Finally, the FLDs for a real AA2024 aluminum material are determined by using different types of initial imperfection in the M-K imperfection approach.

The material parameters for the fictitious material used in the present contribution are listed in Tables 3 and4.

Table 3 Elastic properties and isotropic hardening parameters used in the simulations

E (GPa) ν K (MPa

Geometric imperfection

This subsection outlines the predicted FLDs when the initial imperfection in the M-K approach is assumed in the form of reduced thickness in the localization band, as compared to the thickness in the homogeneous zone. Note that, in the Rice bifurcation approach, no initial imperfection is required to predict localized necking. For different values of the initial thickness imperfection factor, the FLDs predicted by the M-K approach are shown in Fig. 8, along with the FLD obtained by the Rice bifurcation criterion. It can be seen that when the initial geometric imperfection factor increases (i.e., a reduction in the imperfection size), an increase in the limit strains is observed. Furthermore, when the initial geometric imperfection factor becomes equal to 1 (i.e., thickness of the imperfection band becomes equal to that of the homogeneous zone), the FLD obtained with the M-K approach coincides with that given by the Rice bifurcation criterion.

Consequently, the Rice bifurcation criterion provides an upper bound to the limit strains predicted by the M-K approach. This finding is fully consistent with the observations reported in the literature, in which similar analogy between the M-K approach and the Rice bifurcation criterion has been established within the framework of J2 deformation theory of plasticity [START_REF] Yoshida | Comparison of bifurcation and imperfection analyses of localized necking in rate-independent polycrystalline sheets[END_REF][START_REF] Bettaieb | Investigation of localized necking in substratesupported metal layers: Comparison of bifurcation and imperfection analyses[END_REF] and continuum damage mechanics models [START_REF] Abed-Meraim | Investigation and comparative analysis of plastic instability criteria: application to forming limit diagrams[END_REF].

-0.05 0.00 0.05 0.10 0.15 0.20

0.0 0.1 0.2 Rice bifurcation M-K (z 0 =1) M-K (z 0 =0.9999) M-K (z 0 =0.99) M-K (z 0 =0.9)
Major Strain

Minor Strain

Fig. 8 FLDs obtained with the M-K approach (using initial thickness imperfection) and with the Rice bifurcation criterion combined with the GTN damage model

The evolution of the minimum value of the determinant of the acoustic tensor is presented in Fig. 9 along with the evolution of the ratio p p / ɺ ɺ B H ε ε for the initial imperfection factor 0 z = 1.

From this figure, it can be observed that the ratio of the equivalent plastic strain rate, for the M-K approach with no initial imperfection, remains equal to 1 as long as the determinant of the acoustic tensor, for the Rice bifurcation criterion, remains positive. When the condition of loss of ellipticity is met (i.e., det(Q) = 0), the ratio of the equivalent plastic strain rate suddenly becomes very high, as shown in Fig. 9 for the strain path of balanced biaxial tension. Note that similar trends are observed for other strain paths, which are not shown for conciseness. Min det (Q)

Fig. 9 Evolution of the equivalent plastic strain rate ratio for the M-K approach with initial imperfection factor 0 z 1 = , along with the minimum determinant of the acoustic tensor for the Rice bifurcation criterion, under balanced biaxial tension

Material imperfection

In contrast to the original M-K approach, which considers an initial geometric imperfection (i.e., thickness imperfection), Yamamoto [START_REF] Yamamoto | Conditions for shear localization in the ductile fracture of voidcontaining materials[END_REF] has proposed in his pioneering study the introduction of an initial imperfection within the sheet in the form of material imperfection, and more specifically, as an initial porosity imperfection. In this subsection, only material-type imperfections are considered, such as initial porosity imperfection, void nucleation imperfection, void coalescence imperfection and hardening imperfection. The material parameters adopted in the following simulations are the same as those listed in Tables 3 and4 in both the homogeneous zone and the localization band. However, in order to account for a specific initial material imperfection, only one material parameter is varied at a time within the localization band.

Initial porosity imperfection

An imperfection band with slightly increased initial porosity, as compared to the homogeneous zone, is considered here for the M-K approach. The initial material imperfection factor in this case is represented by the difference in initial porosity between the imperfection band and the homogeneous zone, i.e., 0 0 0

B H f f f ∆ = - .
In order to illustrate the effect of the initial porosity imperfection on the material behavior, Fig. 10a depicts the material response under uniaxial tension within the imperfection band for four different values of 0 f ∆ . Higher values of initial porosity are associated with weaker material. As can be seen, higher values of initial porosity imperfection tend to accelerate the porosity evolution, leading to early softening. The corresponding FLDs obtained with the M-K approach are shown in Fig. 10b, along with the FLD predicted with the Rice bifurcation criterion. It should be noted that when 0 0 f ∆ = (i.e., identical initial porosity both in the imperfection band and in the homogeneous zone), the FLD predicted with the M-K approach coincides with that given by the Rice bifurcation criterion, which is consistent with the previous predictions using initial geometric imperfection. On the other hand, when 0 f ∆ is increased, the initial porosity in the imperfection band is greater than that in the homogeneous zone, which ultimately lowers the FLDs. These results are fully consistent with the earlier studies reported in Yamamoto [START_REF] Yamamoto | Conditions for shear localization in the ductile fracture of voidcontaining materials[END_REF] and Needleman and Triantafyllidis [START_REF] Needleman | Void growth and local necking in biaxially stretched sheets[END_REF], in which it has been concluded that the overall shape of the predicted FLD mainly depends on the difference in initial porosity 0 f ∆ , and not on the individual initial porosity in each zone

( ) 0 Δ 0.0015 f = ( ) 0 Δ 0.0005 f = ( ) 0 Δ 0.0002 f = Rice bifurcation M-K M-K M-K M-K Major Strain Minor Strain ( ) 0 Δ 0 f = (a) (b)
Fig. 10 (a) Uniaxial stress-strain curves for different values of initial porosity imperfection, and (b) FLDs predicted with the M-K approach and the Rice bifurcation criterion combined with the GTN damage model

Void nucleation imperfection

Nucleation of voids is considered here as a material imperfection in the M-K approach. To achieve this, the Gaussian amplitude ψ (see Eq. ( 10)), which represents the total volume fraction of cavities that are likely to nucleate, is considered slightly higher inside the imperfection band, as compared to the homogeneous zone. The material imperfection factor in this case is expressed as

B H ψ ψ ψ ∆ = -
, where B ψ and H ψ represent the Gaussian amplitudes inside the imperfection band and the homogeneous zone, respectively. Fig. 11a shows the true stress-strain curve under uniaxial tension within the imperfection band, for four different values of ψ ∆ . Note that the softening regime starts earlier for higher values of ψ ∆ , which is consistent with the physical significance of the nucleation parameter ψ . The corresponding FLDs are presented in Fig. 11b.

It can be observed that the limit strains obtained with the M-K approach are lowered as the considered material imperfection ψ ∆ increases. Moreover, the FLDs predicted with the M-K approach tend to the FLD obtained with the Rice bifurcation criterion when the size of material imperfection ψ ∆ tends to zero, which is consistent with the previous numerical results.

Void coalescence imperfection

Similar to the previous studies, the effect of void coalescence imperfection on the prediction of FLDs is investigated here. Both the critical porosity parameter cr f and the accelerating factor GTN δ , involved in the coalescence mechanism, are considered as separate void coalescence imperfections. Figs. 12a and12c illustrate the effect of the critical porosity and the accelerating factor imperfections, respectively, on the true stress-strain response within the imperfection band under uniaxial tension. As can be seen, significant influence of these coalescence parameters on the softening regime is observed. The corresponding FLDs obtained with the M-K approach are shown in Figs. 12b and12d, for the considered void coalescence imperfections. It should be noted that the overall shapes and trends of the predicted FLDs are similar to what has been obtained with the previous types of material imperfections (see Figs. 10b and11b).

Δ GTN δ = 20
( )

Δ GTN δ = 10
( )

Δ GTN δ = 2.5 Rice bifurcation M-K M-K M-K M-K Major Strain Minor Strain ( ) Δ 0 GTN δ = (c) (d)
Fig. 12 (a) Uniaxial stress-strain curves for different values of critical porosity imperfection, (b) FLDs predicted with the M-K approach and the Rice bifurcation criterion combined with the GTN damage model, (c) Uniaxial stress-strain curves for different values of accelerating factor imperfection, and (d) FLDs predicted with the M-K approach and the Rice bifurcation criterion combined with the GTN damage model

Hardening imperfection

In the earlier work of Needleman and Triantafyllidis [START_REF] Needleman | Void growth and local necking in biaxially stretched sheets[END_REF], it has been shown that qualitatively equivalent FLDs are obtained when either hardening imperfection or initial porosity imperfection is considered in the M-K approach. In order to further investigate this observation, the imperfection band in the M-K approach is assumed with slightly reduced value of Swift's hardening coefficient K (see Eq. ( 4)), as compared to the homogeneous zone. Accordingly, the material imperfection factor in this case is expressed as

H B K K K ∆ = -
, where H K and B K represent the Swift hardening coefficient in the homogeneous zone and inside the imperfection band, respectively. Fig. 13a clearly shows the effect of the hardening imperfection K ∆ on the true stress-strain response within the imperfection band under uniaxial tension. The corresponding FLDs are shown in Fig. 13b. It can be seen that the FLDs are lowered as the initial material imperfection ΔK increases, which is fully consistent with the previous predictions using various types of imperfections. Moreover, the FLD predicted with the M-K approach becomes equivalent to that obtained with the Rice bifurcation criterion when the hardening imperfection vanishes (i.e., 0 K = ∆ ). FLDs predicted with the M-K approach and the Rice bifurcation criterion combined with the GTN damage model

Comparison of the FLDs obtained by using different initial imperfections in the M-K model

From the previous study, it can be clearly concluded that the limit strains are lowered when the size of initial imperfection is increased, whatever the type of initial imperfection. However, in order to determine the most influential initial imperfection factors, and to identify those having a key destabilizing role in triggering plastic flow localization, a comparative study is conducted in Fig. 14. This Fig. 14 illustrates the variation of the normalized major critical strain, predicted with the M-K approach at localization, as a function of the percentage of increase in the initial imperfection factor ξ . Note that, for conciseness, the results are reported for only three particular strain-path ratios, i.e., -0.5, 0 and 1. For the geometric imperfection (i.e., thickness imperfection), the initial imperfection factor ξ is defined in terms of percentage as:

0 z ξ with 0 and     ∆                     H B H B H i i i i i i H H i i h h h h h h h h - = 100 = 100 ≠ ≤ (42)
Similar to Eq. ( 42), the initial imperfection factors associated with the material imperfections are defined as follows:

ξ w       cr B H H B H f H H B H H B H ψ H H H B cr cr cr f H H cr cr f f f f f f f f ψ ψ ψ ψ ψ ψ ψ ψ f f f f f - = 100 = 100; ≠ ≥ - = 100
    ∆       =                   ∆      =              GTN H B H cr cr cr B H H B H GTN GTN GTN δ GTN GTN GTN H H GTN GTN H B H B H K H H f f f δ δ δ δ δ δ δ δ K K K K K K K K ≠ ≤ - = 100 100; ≠ ≥ - = 100 100; ≠ ≤ (43)
From Fig. 14, the following conclusions can be drawn:

•

The critical strains are lowered with the percentage of increase in the initial imperfection factor for all of the considered strain-path ratios. However, the decrease in the critical strains for balanced biaxial tension is comparatively more significant than for the uniaxial tension and plane-strain tension, which is in accordance with the available literature (see, e.g., Hu et al. [START_REF] Hu | Influence of damage and texture evolution on limit strain in biaxially stretched aluminum alloy sheets[END_REF] and Ben Bettaieb and Abed-Meraim [START_REF] Bettaieb | Investigation of localized necking in substratesupported metal layers: Comparison of bifurcation and imperfection analyses[END_REF]).

•

The initial imperfections associated with the thickness and the isotropic hardening have similar effects on the predicted critical strains for all of the strain-path ratios. Additionally, rapid decrease in the critical strains is observed when the initial imperfection ratios 0 z ξ and ξ K are increased, as compared to the damage-related material imperfections. Therefore, it can be concluded that the critical strains, and thus the FLDs, are more sensitive to the initial imperfections associated with the thickness and the isotropic hardening.

•

As to the GTN damage-related material imperfections, the effect of the critical porosity imperfection on the predicted critical strains is the most significant, which reveals that the critical porosity is the most influential parameter on the prediction of FLD among the GTN damage parameters.

FLD prediction for AA2024 aluminum material

In this section, the developed numerical tool based on the M-K imperfection approach and GTN damage model is used to predict FLD for AA2024 aluminum material. The predicted FLDs are also compared with the experimental FLD reported in Janbakhsh et al. [START_REF] Janbakhsh | A comparative study on determination of forming limit diagrams for industrial aluminium sheet alloys considering combined effect of strain path, anisotropy and yield locus[END_REF]. The elastic properties, Swift's isotropic hardening parameters and GTN damage parameters for the studied AA2024 aluminum material have been provided by Wu et al. [START_REF] Wu | An extended GTN model for low stress triaxiality and application in spinning forming[END_REF] and are listed in Tables 5 and6. Figure 15 shows the FLDs predicted using three different types of M-K initial imperfections: thickness, isotropic hardening and critical porosity initial imperfections. These three specific initial imperfections are the most influential parameters on the prediction of FLD, as concluded from the previous investigation (see Section 5.3). Compared to the experimental FLD taken from Janbakhsh et al. [START_REF] Janbakhsh | A comparative study on determination of forming limit diagrams for industrial aluminium sheet alloys considering combined effect of strain path, anisotropy and yield locus[END_REF], Fig. 15a shows that an initial thickness imperfection factor of 0 z 0.995 = leads to excellent agreement with the experimental limit strains. Similar to the consideration of initial thickness imperfection, Fig. 15b shows that the FLD predicted with an initial hardening imperfection of Δ Pa K = 7 Μ is in good agreement with the experimental one.

The effect of critical porosity imperfection on the prediction of FLD is also investigated in Fig. 15c. In contrast with the effect of initial thickness and hardening imperfections on the prediction of FLDs, the critical porosity imperfection has a significant effect on the overall shape of the predicted FLDs (see Fig. 15c). In addition, the predicted limit strains are comparatively closer to the experimental ones for Δ 0.003 cr f = . On the whole, it can be concluded that the above three initial imperfection factors play a crucial role in the prediction of FLD, and their magnitudes should be determined carefully. Also, the experimental identification of these initial imperfections must be combined with that of the material parameters, in order to provide accurate predictions of limit strains.

-0.14 -0.07 0.00 0.07 0.14

Conclusions

In the present contribution, a numerical tool has been developed to predict FLDs by combining the GTN damage model with the M-K imperfection approach and the Rice bifurcation criterion. For the M-K imperfection approach, the effects of thickness-based geometric imperfection as well as various types of material imperfections on the prediction of FLDs have been investigated for a fictitious material. The predicted limit strains have been systematically compared to those obtained with the Rice bifurcation criterion, for which no initial imperfection is required. It has been shown that when the geometric or material imperfection tends to zero (i.e., vanishing initial imperfection), the FLD predicted with the M-K approach tends to that obtained with the Rice bifurcation criterion. Moreover, the limit strains are lowered when the initial imperfection magnitudes are increased for the M-K approach, which is consistent with the basic foundation of this criterion. Accordingly, the Rice bifurcation criterion provides an upper bound for the prediction of strain localization, as compared to the M-K approach. Therefore, the Rice bifurcation criterion could be used as a "qualitative tool" for the prediction of FLDs with respect to experimental data. It can also be used, in the design stage of new materials with enhanced ductility, as an efficient tool for the comparison of formability of several sheets belonging to the same grade of material. As to the M-K approach, the associated numerical tool provides a fairly accurate "quantitative comparison" of the predicted limit strains with respect to experimental measurements. The concept of initial imperfection factor in the M-K approach encompasses all defects and destabilizing mechanisms typically present inside the sheet metal during forming process. In such a way, the proposed numerical tool provides a cost-effective procedure for the determination of FLDs for any given materials.

A parameter sensitivity study has been conducted in order to identify the most influential initial imperfection factors on the prediction of limit strains with the M-K approach. It has been shown that the thickness-based imperfection type as well as the hardening imperfection type have a significant impact on the prediction of limit strains, as compared to damage-based imperfections. As to the latter damage-related material imperfections, the effect of the critical porosity imperfection on the predicted FLDs is the most significant, which highlights the importance of the critical porosity parameter in the interaction between damage evolution and strain localization.

The present approach has additionally been used for the prediction of the FLD of the AA2024 aluminum material. It has been shown that a careful selection of the initial imperfection factor leads to excellent agreement between the predicted limit strains and the experimental ones. These initial imperfection factors, which can be regarded as additional material parameters, could be determined experimentally within an identification procedure including the intrinsic material parameters.

Funding:

  Not applicable Conflicts of interest/Competing interests: The authors declare no conflict of interest. Availability of data and material: Not applicable Code availability: ABAQUS, MATLAB Nomenclature f total void volume fraction (also called porosity) 0 f initial void volume fraction (also called initial porosity) Σ macroscopic Cauchy stress tensor ' Σ deviatoric part of the macroscopic Cauchy stress tensor Σ m Σ hydrostatic part of the macroscopic Cauchy stress tensor Σ

•

  ψ volume fraction of inclusions that are likely to nucleate N ε equivalent plastic strain for which half of inclusions have nucleated N s standard deviation on N ε ep C analytical elastic-plastic tangent modulus defined in the corotational material frame Σ M partial derivative of the yield surface Y Φ with respect to Σ * M f partial derivative of the yield surface Y Φ with respect to * f σ M partial derivative of the yield surface Y Φ with respect to σ e C fourth-order elasticity tensor D total macroscopic strain rate tensor 0 θ initial orientation of the imperfection band with respect to the global axis θ orientation of the imperfection band with respect to the global axis at time t 0 n unit vector normal to the imperfection band at the initial time n unit vector normal to the imperfection band at time t component ij of tensor quantity • i j k l • component ijkl of tensor quantity • n • value of quantity • at the start of loading increment 1 n+ • value of quantity • at the end of loading increment • B value of quantity • inside the imperfection band • H value of quantity • inside the homogeneous zone Mathematical operators Vectorial and tensorial fields are designated by bold letters and symbols Scalar variables and parameters are represented by thin letters and symbols Einstein's convention of summation over repeated indices is adopted •⋅• Product between two second-order tensors ( i j j k = • • )

  is a necessary outcome for strict plastic loading (i.e., 0 ɺ λ > ). The expression Y Φ 0 ɺ = represents the consistency condition and can be used to derive the elastic-plastic tangent modulus ep C as follows:

H G and 33 BG

 33 of the velocity gradients are calculated so that the planestress conditions (i.e.,[START_REF] Madou | A Gurson-type criterion for porous ductile solids containing arbitrary ellipsoidal voids-II: Determination of yield criterion parameters[END_REF] Σ 0 H = and 33 Σ 0 B = ) are satisfied. The equilibrium of forces across the imperfection band can be expressed by the following equation:

11 H

 11 the threshold value of 10. The corresponding values of the principal strains in the homogeneous zone, i.e.

Fig. 5 Fig. 6 Fig. 7

 567 Fig.[START_REF] Hill | A theory of the yielding and plastic flow of anisotropic metals[END_REF] Comparison of stress-strain curves obtained with the built-in ABAQUS model and the developed UMAT subroutine, for balanced biaxial tension

Fig. 11 (

 11 Fig. 11 (a) Uniaxial stress-strain curves for different values of void nucleation imperfection, and (b) FLDs predicted with the M-K approach and the Rice bifurcation criterion combined with the GTN damage model

Fig. 13 (

 13 Fig.[START_REF] Banabic | An improved analytical description of orthotropy in metallic sheets[END_REF] (a) Uniaxial stress-strain curves for different values of hardening imperfection, and (b) FLDs predicted with the M-K approach and the Rice bifurcation criterion combined with the GTN damage model

Fig. 14

 14 Fig.[START_REF] Cazacu | New yield criteria for isotropic and textured metallic materials[END_REF] Normalized major critical strain at localization, as a function of percentage of increase in initial imperfection factor: (a) uniaxial tension (UT), (b) plane-strain tension (PS), and (c) balanced biaxial tension (BBT)

Fig. 15

 15 Fig.[START_REF] Cazacu | New mathematical results and explicit expressions in terms of the stress components of Barlat et al. (1991) orthotropic yield criterion[END_REF] Predictions of the FLD for AA2024 aluminum material using: (a) initial thickness imperfection, (b) hardening imperfection and (c) critical porosity imperfection, in the M-K imperfection approach combined with the GTN damage model

Table 1

 1 Elastic properties and isotropic hardening parameters used in the simulations

	E (GPa)	ν	K (MPa)	n	ε	0
	70	0.33	371.2	0.17	0.00324

Table 2

 2 Growth and nucleation parameters used in the simulations

	f	0	ψ	ε	N	s	N	1 q	q	2	q	3
	0.001 0.27	0.27	0.1	1.5	1	2.15

Table 5

 5 Elastic properties and isotropic hardening parameters for AA2024 aluminum material

			E (GPa)		ν			K (MPa)		n		ε	0
			71.5		0.33			962.14		0.1993	0.0109
	Table 6 GTN damage parameters for AA2024 aluminum material
	f	0	ψ	ε	N	s	N	1 q	q	2	3 q		f	cr	GTN δ
	0.00	0.025 0.122 0.198 1.5	1	2.15	0.014 18.19