Experimental Analysis and Theoretical Predictions of the Limit Strains of a Hot-dip Galvanized Interstitial-free Steel Sheet

In this work, the formability of a hot-dip galvanized interstitial-free (IF) steel sheet was evaluated by means of uniaxial tensile and Forming Limit Curve (FLC) tests. The FLC was defined using the flat-bottomed Marciniak’s punch technique, where the strain analysis was made using a digital image correlation software. A plastic localization model was also proposed wherein the governing equations are solved with the help of the Newton’s method. The investigated hot-dip galvanized IF steel sheet presented a very good formability level in the deep-drawing range consistent with the measured Lankford values. The predicted limit strains were found to be in good agreement with the experimental data of the hot-dip galvanized IF steel sheet owing to the definition of the localization model geometrical imperfection as a function of the experimental surface roughness evolution and, in particular, to the yield surface flattening near to the plane-strain stress state authorized by the adopted yield criterion.


Introduction
The concept of the Forming Limit Diagram (FLD) was introduced by Keeler 1 in order to assess the formability of sheets under biaxial tension. This earlier study applied for positive values of the smaller principal strain in the plane of the sheet. It was next extended by Goodwin 2 to the full range between uniaxial and equibiaxial tension. Since that time, a very large number of investigations have been devoted to the experimental determination and the theoretical modeling of limit strains in sheets under biaxial tension. The FLD is defined in the axes of minor and major principal strains in the plane of the sheet. The curve obtained by plotting the limit strains obtained for linear strain paths is the Forming Limit Curve (FLC). The FLC is commonly observed to be fairly independent of the test used for its determination, provided the strain paths remain almost linear up to the occurrence of intense flow localization. More generally, the limit strains are assumed by most researchers to represent an intrinsic, although strain-path dependent, property of the material. In other words, structural effects relating to the boundary conditions of the deformation process are implicitly assumed to be without influence on the limit strains.
The prediction of the forming limits set by the occurrence of localized necking in plastically stretched sheets has been the subject of a very large number of theoretical analyses. The most popular approach is based on the model originally proposed by Marciniak and Kuczynski 3 . This approach, hereafter referred to as the M-K model, assumes the existence of an initial thickness imperfection across the sheet. Afterwards, the forming limits result from the process of flow localization in the defective region. Almost all the available descriptions of plastic yielding have been implemented in such a simple plane-stress analysis involving the existence of two individual zones, namely, the homogeneous zone, and the defective, thinner zone. In addition to the beneficial effects of strain-hardening and strain-rate hardening to delay the development of the neck 4,5 , it is now well established that the forming limits strongly depend on the plasticity model adopted, and in particular on the shape of the yield surface [6][7][8] . Thus, considerable advances have been achieved during the last twenty years in the prediction of the sheet metal forming limits by employing yield criteria proposed to provide an improved description of plastic flow under conditions of biaxial loading. For instance, predictions in good agreement with experimental results have been obtained in the case of aluminum alloys 9-10 with the yield criteria proposed by Barlat et al. 11 , Karafillis and Boyce 12 and Barlat et al. 13 .
The effects of the through-thickness normal and shear stress components have also been taken into account within the framework of the M-K model, see the recent works [14][15][16][17][18] . These out-of-plane effects are not negligible in some sheet metal forming processes, such as hydroforming and incremental sheet forming, where significant through-thickness compressive stress components may take place.
On the other hand, the assumption of an ad-hoc groove-like imperfection across the sheet is questionable, especially as the M-K analysis is very sensitive to the size of this imperfection. Several attempts have hence been made to rationalize the physical origin of the M-K imperfection. Its initial interpretation in terms of variations in sheet thickness (geometric imperfection) was amended by considering also possible variations in the strength of the material (material or metallurgical imperfection). This idea was already mentioned in the work of Marciniak et al. 19 , and it was later materialized by assuming spatial variations in the orientation distribution of grains, giving variations in the Taylor factor 20 . Inhomogeneous damage due to the nucleation and growth of microvoids can also occur, for instance by decohesion between matrix and inclusion, or by fragmentation of inclusions, leading to a decrease in the surface area transmitting the load 21 . Indeed, the flow localization process in the presence of an initial thickness imperfection is accelerated by the development of damage, which is more rapid in the thinner, more stretched region 22,23 .
A quasi-random arrangement of hard elastic inclusions was also considered in finite element (FE) calculations to generate inhomogeneous strain fields leading to necking of the representative volume element 24,25 . Wu et al. 26 performed FE computations of the response of a macroscopically infinitely small unit cell made of an assembly of elements, where each element represents a grain following the constitutive laws of crystal plasticity. The orientation distribution of grains is representative of a given texture. Then, the flow localization process can be analyzed at this mesoscopic level under a given average stress-ratio, and compared with the M-K analysis results.
The present work aims at evaluating first the formability of a hot-dip galvanized interstitial-free steel sheet using the Marciniak flat-bottomed punch test technique, where the limit strains were defined with the help of a Digital Image Correlation (DIC) analysis software. Surface roughness measurements were also performed on each FLC blank specimen width in attempt to make a correlation to the level of the accumulated effective plastic-strain. Also, a M-K model is developed within the framework of the flow-theory under the assumption of isotropic work-hardening together with a phenomenological orthotropic plasticity yield criterion. The geometrical imperfection parameter of the M-K model is described as a function of the sheet roughness, mean grain size and the effective plastic-strain. The proposed M-K model is first validated with available data of a low carbon steel sheet, for which FLC tests were performed under linear and bi-linear strain-paths and then, applied to forecast the experimentally obtained limit strains of the investigated hot-dip galvanized interstitial-free steel sheet.

Material
The experimental procedure was performed with a hot-dip galvanized interstitial-free (IF) steel sheet produced at the CSN steel plant in Brazil, which has a nominal thickness of 0.74 mm. Table 1 presents the IF steel sheet  chemical composition obtained by optical emission and  atomic absorption spectrometry (C, S and N). The IF hot-dip galvanized steel microstructure obtained by optical metallography is shown in Figure 1, which is composed by equiaxed ferrite grains with 8.5 ASTM grain size (19 µm).

Uniaxial tensile testing
The uniaxial tensile tests were performed with an Instron 150 kN universal testing machine equipped with a non-contacting video extensometer. The sheet metal specimens were stamped as per NBR ISO 6892:2002 27 with a gauge width of 12.5 mm and 50 mm of gauge length. The hot-dip galvanized IF steel sheet mechanical properties were evaluated from 3 specimens at angular orientations α = 0, 45 and 90 0 with respect to the rolling direction using a constant cross-head rate of 10 mm/min up to the failure. The yield stress was defined by the 0.2% plastic-strain offset method. Also, the true stress (σ) -plastic strain (ε p ) data along the rolling direction was fitted to the modified Swift power law in which the strain-rate effect is considered as 19 : where: K is the strength coefficient (MPa), ε 0 is the pre-strain, N is the strain-hardening exponent, M is the strain-rate sensitivity index, p ε  is the true plastic-strain rate (s -1 ) and 0 ε  is the reference strain-rate (s -1 ). The strain-rate sensitivity was determined using constant cross-head rate of 1 and 10 mm/min, which correspond approximately to nominal strain-rates of 3 × 10 -4 and 3 × 10 -3 s -1 , respectively. The strain-rate sensitivity was defined from the tests at the rolling direction as:   (2) where: subscripts 1 and 2 denote the measures at the lower and higher strain-rates, respectively. From the video measurements of both elongation and width changes together with the cell load recording, the Instron tensile testing machine employs automatic methods to define either the strain-hardening exponent N and the Lankford coefficient or r-value by means of a linear regression fitting of the calculated true stress and true plastic-strain data between 10 and 15% (N 10-15 ) and 2 and 20% (r 2-20 ) of longitudinal plastic-strain, respectively.

Marciniak punch test
Numerous tests have been used to determine the Forming Limit Curve (FLC). Nowadays, there is a trend towards determining the FLC by using a single equipment with samples of different widths. This can be achieved with the Nakazima 28 and the Marciniak 19 tests, where the blank sheet is fixed at its periphery, and deformed by a punch, the shape of which is either spherical (Nakazima) or cylindrical flat-bottomed (Marciniak). The Marciniak test was adopted in this study to obtain the FLC of the hot-dip galvanized IF steel sheet. The tooling geometry is defined by the following parameters: flat-punch diameter, 140 mm, flat-punch nose radius, 10 mm, die-opening diameter, 147 mm, die shoulder radius, 2.5 mm and lockbead diameter, 165 mm, respectively. The Marciniak punch experiments were carried out on an Amsler 1000 kN capacity single-action hydraulic press. The blank length L is equal to 200 mm which has been taken perpendicular with respect to the rolling direction. The blank width W is equal to 80, 90, 100, 120, 130, 150, 160, 170, 180 and 200 mm so as to cover roughly the drawing (ε 1 > 0 and ε 2 < 0) and the stretching (ε 1 > 0 and ε 2 > 0) regions of the FLC. The geometries and dimensions of the blank specimens are shown in Figure 2. Carrier or counter-blanks with the same blank material and geometries were also machined so as to prevent the contact between the blank and the flat-punch surfaces and to assure a homogeneous straining. Central holes with diameters equal to 50 and 60 mm were machined in the carrier blanks for the blank specimen widths W = 200 mm and W = 180, 170 and 160 mm, respectively, whereas counter-blanks were cut in two parts (perpendicular to the blank length) for the blank specimen widths W = 80, 90, 100, 110, 130 and 150 mm. All experiments were conducted with a 0.10 mm Teflon sheet between the punch and the carrier-blank, according to the ISO 12004-2:2008 29 standard. The Marciniak punch tests were stopped immediately after a drop in the punch load in order to avoid the complete separation of the samples.

Strain and roughness measurements
The limit strains are generally determined by carrying out a test up to ductile fracture and, then by analyzing the strain distributions obtained in the vicinity of the fractured zone. The strains are commonly measured using circle or square grids printed or electrochemically etched on the surface of the sheet. The Hecker's method 30 defines the limit strains as the limiting values between principal surface strains, namely, ε 1 and ε 2 , measured on the necked or fractured sites, and on adjacent circles or squares. The accuracy of strain analysis can also be improved with the help of a digital image correlation (DIC) technique to determine the displacements fields. In the present work, the strain analysis was made using electrochemically etched 2.5 mm square grid with the help of the ASAME target model system 31 using both ASTM E 2218 -02 32 and ISO 12004-2:2008 29 standards to define the limit strains. Figure 3 schematically shows the procedure adopted according to the ISO 12004-2:2008 standard to define the limit strains for the specimen blank width W = 120 mm. First, at least three intersections lines nearly perpendicular (within ± 15 0 ) to the fracture must be defined to obtain the true major and minor surface principal strains distributions on either sides of the fracture. For each intersection line, the true strains are then plotted as function of the grid position or node number in order to reject the points subjected to the localized necking. This is achieved by first calculating the local second derivatives from a parabola fitting of each major and minor true strains located at left and right sides of the fracture, wherein a range of 5 points is considered excluding the corresponding fracture strains, as shown in Figure 4a. The remaining nodes including the ones located on either sides of the border necking area, see Figure 4b, are used to fit the 6 th order polynomial. Then, the section limit strains are obtained by replacing the corresponding peak strain node in the fitted 6 th order polynomial equation. Finally, the major and minor limit strains (ε 1 , ε 2 )are defined as the averaged values of the selected sections and the procedure is repeated for each sample width.
Surface roughness measurements were also performed in both as-received and deformed conditions. These measurements were carried out with a portable surface roughness Mitutoyo tester Surftest SJ-301 with a tip radius of 5 µm using a cut-off length of 2.5 µm and the roughness average (R a ) parameter. Three measurements were made in each Marciniak punch specimen width at two equidistant grid points located on either sides of the fractured regions. The surface roughening was then correlated to the level of the surface effective plastic strain, ε , the initial sheet surface roughness, R 0 , and the initial grain size, d 0 , as 33 : In Equation 3, the parameter C is determined from the linear fitting of the plot of R -R 0 versus 0 d ε data. The accumulated effective plastic strain was calculated assuming Figure 3. Schematic of the digital image correlation technique adopted to define the sections measurements from the major and minor in-plane principal strains distributions determined for the blank width W = 120 mm. normal plastic anisotropy along with the corresponding Hill's quadratic 34 conjugated effective plastic strain measure defined under plane-stress conditions as: where: is the average normal plastic anisotropy coefficient while ε 1 and ε 2 are the principal surface plastic strains determined from the ASAME measurements.

Constitutive equations
The constitutive equations are defined within the framework of small strains from the additive decomposition of the total strain-rate tensor, ij ε  , in an elastic part, ij e ε  , and a plastic part, p ij ε  , which can be cast in an incremental form as: In this work, the elastic strains are neglected in the solution of the localization model adopted to predict the plastic limit strains in sheet metal forming processes. Thus, in order to abridge the notation the superscript p is omitted hereafter. The plastic strain components are defined from the associated plastic flow-rule as: Where: dλ is the plastic multiplier, f is the yield function defined under the assumption of isotropic work-hardening as: In Equation 7, F(σ ij ) is a first degree homogeneous stress function which defines the yield surface shape in the stressspace whereas σ is a scalar measure of the effective stress that defines the yield surface size and can be identified by the measures of effective plastic-strain, ε , and plastic-strain rate, ε  , determined from in the uniaxial tensile test. In this description, the plastic multiplier in Equation 6 can be obtained from the plastic loading condition, that is, f = 0, by firstly applying the Euler's identity of homogeneous functions of degree n to the yield function given by Equation 7, that is: n x x (8) and then by calculating the associated plastic-work per unit volume from the equivalent plastic-work principle, namely: Together with the associated flow-rule law, Equation 6, which provides that the plastic multiplier is equal to the effective plastic-strain increment conjugated of the effective stress σ , that is, d d The yield function adopted in this work corresponds to the orthotropic plasticity criterion proposed by Ferron et al. 35 defined for plane-stress conditions as: w h e r e : x 1 = ( σ 1 + σ 2 ) / 2 = ( σ x x + σ y y ) / 2 a n d define the center and the radius of the Mohr's circle, respectively, as a function of the in-plane principal stress components (σ 1 ,σ 2 ) or else the in-plane orthotropic stress (σ xx ,σ yy ,σ xy ). In Equation 10, α = (x,1) = (x,2) is the angle of the orientation between the in-plane principal stress directions (1,2) and the in-plane orthotropic symmetry axes (1,2). For numerical implementation purposes, Ferron's orthotropic plane-stress function is defined as: In Equation 11, the exponents m, n, p and q are previously known while the material parameters A, B, k, a and b can be obtained in two steps. Firstly, the parameters (A, B and k) are calculated from the stress-ratios (σ b /σ 45 ) and (σ b /τ0) between the equibiaxial yield stress, σ b , and the uniaxial yield stresses at 45 0 with respect to the rolling direction, σ 45 , and the shear yield stress, τ0, along the rolling direction, respectively, along with the plastic anisotropy coefficient, r 45 . For steel sheets, the recommended values for the exponents are m = 2, n = 1 or 2, p = 1 or 2 and q = 1 or 2 35 . In particular, a positive value for the parameter k allows a flattening of the yield surface near plane-strain tension/compression and pure shear stress states. Secondly, the parameters a and b describing the initial in-plane sheet metal anisotropy, can be obtained from the plastic anisotropy coefficients (r 0 and r 90 ) or else from the uniaxial yield stresses (σ 0 ,σ 45 ,σ 90 ).
When only uniaxial tension data is available, the parameters defining the yield surface at 45 0 from the rolling direction are determined assuming B = 3A and positive k-value. Then, the parameters defining the initial planar material anisotropy, a and b in Equation 11, can be calculated with the r 0 and r 90 values or else from the stress-ratios (σ 45 /σ 0 ,σ 90 /σ 0 ). Also, Hill's quadratic 34 anisotropic yield criterion can be obtained as a particular case of Equation 11 by setting k = 0, m = 2 and n = p = q =1. The yield loci of Ferron's criterion are plotted in Figure 5 as a function of the angle α in the principal stress space (σ 1 ,σ 2 ) normalized by the equibiaxial yield stress (σ b ). In this type of representation, g(θ,α) stands for the normalized radius defining a point on the yield locus where θ is the corresponding polar angle which defines the current stress state. The function g(θ,α) is defined as an extension of Drucker's 36 isotropic yield criterion as: Where: The principal plastic-strain increments are defined with the help of the parametric description, Equation 12, as 35 : Where: ( , ) . From the metal plasticity incompressibility condition, the plastic anisotropy coefficient defined from an uniaxial tensile test at an angular orientation α with respect to the rolling direction, i.e., r α = dε 2 / dε 3 = -dε 2 /(dε 1 + dε 2 ), can be obtained from Equations 14 and 15 as: where: θ = π/4 defines the uniaxial stress state parallel to the in-plane principal stress σ 1 , as shown by the polar-coordinate description in Figure 5.
On the other hand, the plastic-strain increments in the orthotropic symmetry axes frame (x,y,z) are determined from Equation 6, which depends upon the effective plasticstrain increment and the yield function partial derivatives. The later are obtained by applying the consistency condition, df = 0, to the Ferron's yield function defined in the form of Equations 9 and 10, that is: (17) and then, by defining the terms dx 1 , dx 2 and dα as a function of the variables (x 1 ,x 2 ) along with the help of the relations sin 2 α = σ xy /x 2 and cos 2α = (σ x -σ y )/(2x 2 ), which provides: While the partial derivative corresponding to the sheet thickness F, zz = -(F, x + F, y ) is obtained from the plastic incompressibility condition, that is, dε zz = -(dε x + dε y ).

Marciniak-Kuczynski model
The localization model developed in this section is an extension of the original theory proposed by Marciniak and Kuczynski 3 . This model is based upon the existence of an initial thickness imperfection across the sheet, as depicted in Figure 6, where the superscripts a and b denote the homogeneous and the defective, respectively. Figure 6 indicates the coordinate systems axes adopted in the present M-K model, where a common frame (O, n, t, z) is defined by the normal direction n and tangential t directions to the imperfection and the in-plane normal axis of orthotropy symmetry z. Also, it is assumed that the in-plane orthotropic symmetry axes are coincident to the in-plane principal stress directions in the homogeneous zone a, that is, (x a ≡ 1 a , y a ≡ 2 a ). In the defective zone, the orthotropy directions are not coincident with the principal stress directions. This can be taken into account by defining the orthotropic axes (x b , y b ) in a corotational frame, which rotates as a function of the plastic deformation process 37 , and is defined by the angle ϕ ≡ (x a , x b ) ≡ (y a , y b ) shown in Figure 6. The rotation of the corotational frame axes (x b , y b ) with respect to the principal stress directions in the homogeneous zone a, is defined by 37  The M-K geometrical imperfection has an initial angular orientation defined by the angle: Ψ ≡ (x a , n) ≡ (y a , t), see Figure 6. This angle evolves during the plastic deformation process, which rotation is given by 38 : where: a 1 dε and a 2 dε are the principal strain increments in the homogeneous zone a.
The initial imperfection value is defined by the ratio between the initial thicknesses of the defective zone, 0 b t , and homogeneous zone, 0 a t , as: In this work, the geometrical imperfection is assumed to be dependent on the current sheet roughness R, which is defined by Equation 3, as 33 : In Equation 25, the effective plastic strain measure b ε in the geometrical imperfection zone is used in Equation 3 and, thus, the current geometrical imperfection is defined as: The localization problem within the framework of the M-K analysis is based on the equilibrium of forces between the two zones along the normal (n) and tangential (t) directions of the geometrical imperfection, that is: Which can be written as force per unit width as function of the stress components along the imperfection axes (n,t) and the actual thickness of the two zones, that is: Also, the geometrical compatibility condition imposes that the strain increments along the tangential direction t must be the same between the zones a and b, namely, The solution of the M-K problem is obtained assuming firstly a given strain-ratio defined in the homogeneous zone between the minor and major in-plane principal strain increments, which is expressed by Ferron's yield criterion, see Equations 14 and 15, as: ε θ α π + θ + θ α π + θ ′ ρ = = θ α π + θ -θ α π + θ ′ ε (32) Figure 6. Schematic of the Marciniak-Kuczynski localization model with the homogeneous and geometrical imperfection zones. and, then solving Equation 32 for the angle θ to determine the corresponding stress-state. In the present work, the strain-ratio is varied from the uniaxial tensile stress-state ( / 4 a UT θ = π ) to the equibiaxial straining mode defined by ρ a = 1( a a ES θ = θ ). Secondly, the homogeneous zone a is loaded incrementally for a fixed strain-ratio (or stress-state) with a small effective-strain increment value ( a   4   d 10ε = ). In this way, the stress and strain components in the homogeneous zone a can be straightforwardly determined. With Ferron's orthotropic yield function written in parametric coordinates, see Equation 12, and according to the Figure 6 with the variables 1 In Equations 33 and 34 a σ is the effective stress measure, which in Ferron's yield criterion is identified as the equibiaxial yield stress σ b . Thus, it is more appropriate to define the stress-strain effective measures in the workhardening law given by Equation 1 under the equibiaxial tension stress-state. This is achieved noting from Figure 5 that the ratio between the equibiaxial yield stress and the uniaxial tension yield stress at an angular orientation α with respect to the rolling direction is defined by Equations 33 and 34 as and then by means of the transformations , Equation 1 can be rewritten as: Moreover, the stress and strain increments in the geometrical imperfection axes can be obtained from the following transformation equations: The unknowns variables in the weaker zone b, namely, , are then numerically calculated from the stress and strain states in the homogeneous zone a.
Given that the main governing equations of the M-K problem are defined by two equilibrium forces and one deformation compatibility condition, given by Equations 29-31, an additional relationship is derived from the equivalent principle of plastic-work in the zone b, see Equation 8, so as to obtain a system of four nonlinear equations as proposed by Ganjiani and Assempour 39 : The solution of the system of nonlinear equations . Moreover, a backtracking algorithm is adopted to assure the convergence of the full Newton's method which strongly depends on the initial guess 39,40 . The localized necking criterion is defined when the ratio between the effective plastic-strain increments in the two zones reaches the condition of an unstable flow, which is given by ε , determined as a function of the geometrical imperfection orientation Ψ, which, in turn, is varied between 0 and 90 degrees. Finally, the complete FLC prediction is obtained by varying the strain-ratio in the homogeneous zone ρ a between the uniaxial tensile stressstate and the equibiaxial straining mode.
The effects of non-linear strain-paths are considered assuming bi-linear strain-paths in two stages. In this way, it is possible to analyze the influence of the strain-path upon the forming limits under uniaxial tension, plane-strain tension and equibiaxial stretching pre-straining stages followed by linear proportional strain-paths. The strain components as well the effective plastic strain and the imperfection angular orientation at the beginning of the second stage (2) are initialized with the values determined at the end of the first stage (1):

Plastic properties
The uniaxial tensile tests were carried out at constant cross-head speeds of 1 and 10 mm/min. The average values determined from 3 tests at each angular orientation (α) of the 0.2% off-set yield stress (σ y ), the strain-hardening exponent (N), the ultimate tensile strength (σ u ), the total elongation (∆l) and the Lankford values (r) are listed in Table 2, where the italic values indicate the corresponding standard deviation. The yield stress value increases monotonically between α = 0 and α = 90 degrees, that is, σ 90 > σ 45 > σ 0 . As expected, this stress anisotropy is preserved for the ultimate tensile strength values and is also consistent with the angular evolution of the plastic anisotropy coefficient or Lankford values, which decreases monotonically between α = 0 and α = 90 0 (r 0 > r 45 > r 90 ). The average normal plastic anisotropy coefficient  Table 3 presents the average values of the material parameters fitted according to the work-hardening law defined by Equation 1. The parameters K, ε 0 and N in Equation 1 were obtained using a least-squares method from the uniaxial tensile true stress-true plastic strain data at the rolling direction. The IF steel sheet strain-rate sensitivity parameter M was defined as an average value from the fitted average instantaneous strain-rate sensitivity curve, calculated using Equation 2, as a function of the true plastic-strain level, as shown in Figure 7. The predicted true stress-true plastic strain curves computed with the fitted parameters for the IF hot-dip galvanized steel sheet are in good agreement with the experimental data plotted as a function of the uniaxial tensile test cross-head rate, as shown in Figure 8. Figure 9 shows the photograph of the blanks deformed in the Marciniak punch test, where it is possible to observe that the fractured sites are near to the specimen center. Figure 10 presents the contour values of the major and minor principal strains obtained with the ASAME software from the regions    of interest of each specimen width. The maxima and minima values of the major and minor principal strains (ε 1 ,ε 2 ) varied from (0.87,-0.65) and (0.48, 0.37) for the specimen widths equal to 80 and 200 mm, respectively. It should be noted that none of the probable specimen blank geometries, namely, the widths of 130 and 150 mm, provided the forming limits close to the plane-strain intercept (FLC 0 ). These unexpected results can be ascribed to the irregular lubrication between the two-piece carrier-blank and the flat-bottomed punch, due to the use of a very thin Teflon sheet (0.1 mm thick), and to the insufficient friction condition between the two-piece carrier-blank and the corresponding specimen blank. Figure 10. Major and minor principal strains contour plots determined from the ASAME software analysis as a function of the blank width.

Limit strains and surface roughness
The limit strains determined for the hot-dip galvanized IF steel sheet according to ASTM E 2218-02 32 and ISO 12004 -2:2008 29 standards are shown in Figure 11. Firstly, the limit strains values defined as peer ASTM E 2218-02 32 along with the ASAME DIC 31 software present a very large scattering owing to its strong dependence upon the user's criterion to select the necked and fractured areas at each specimen width. As pointed out, the conducted Marciniak punch tests do not provided points in the FLC near to either sides of the plane-strain intercept due to the uneven friction conditions. The software trend forming limit curves at necking and fracture conditions are also shown in Figure 11a, providing a forecast for the plane-strain intercept FLC 0 . Conversely, the limit strains defined according to ISO 12004 -2:2008 29 , see Figure 11b, clearly show the couple of strains determined for each specimen thanks to the rejection and fitting procedures to delimitate the border of the necking zone from either sides of the peak strain. It should be noted that the fracture limit strains were obtained from the reconstruction of the grid pattern performed with the ASAME software and as the peak strain using the ASTM E 2218-02 32 and ISO 12004 -2:2008 29 standards, respectively.
The IF hot-dip galvanized steel sheet presents a very good formability level in the deep-drawing range thanks to its elevated plastic anisotropy coefficients.
Due to the lack of the experimental limit strains close to the plane-strain intercept, the empirical formula proposed by Keeler and Brazier 41 is adopted here to forecast the engineering FLC 0 value: Where: t is the sheet thickness expressed in mm a n d i s t h e ave r a g e strain-hardening exponent calculated from the uniaxial tensile data listed in Table 2. Henceforth, the corresponding true plane-strain intercept predicted value, , will be adopted together with the limit strains determined according to the ISO 12004 -2:2008 29 for comparison purposes with the theoretical predictions of the present M-K model. Figure 12a correlates the resulting effective plastic-strains calculated with Equation 4 and the measured average roughness (Ra) values as a function of the specimen width. The average initial roughness of the IF hot-dip galvanized  steel sheet is R 0 = 0.93 µm and changed to an average value of 3.55 R = µm due to the in-plane plastic straining. Figure 12b shows the plot of the R -R 0 versus 0 d ε data and the linear fitting according to the Equation 3, from which the resulting C parameter is equal to 0.773 µm 1/2 .

Validation of the M-K model
In order to validate the M-K model developed in the present work, see Section 3.2, we first present the theoretical predictions for linear and bi-linear strain-paths assuming material models as a function of the plastic anisotropy coefficient. Then, the predictions determined with Hill's quadratic 34 and Ferron's 35 yield criteria are compared to the experimental data obtained from linear and bi-linear strain-paths for an AISI-1012 steel sheet with a nominal thickness of 2.5 mm [42] . Table 4 lists the sheet material plastic properties adopted in all the theoretical predictions discussed in the present section. From Equation 25, the initial value of the geometrical imperfection f 0 for the AISI-1012 steel sheet is equal to 0.995. Figure 13a presents linear strain-paths predictions determined as a function of the value of the plastic anisotropy coefficient at 45 0 with respect to the rolling direction. Three material models have been defined assuming r 45 equal to 0.5, 1.0 and 2.5 together with r 0 = r 90 =1.0, respectively. The present material parameter combinations have already been investigated by Barata da Rocha et al. 38 and were adopted here to verify the consistence of the present M-K model predictions. For the planar-anisotropic materials, the numerical simulations were performed with Hill's quadratic yield criterion 34 while the von Mises yield criterion was adopted for the isotropic case. As expected, the same limit strains values are predicted under the uniaxial tensile strain-mode since all the adopted materials models have the same plastic anisotropy coefficient in the rolling direction and, thus, . Conversely, in the biaxial stretching domain the limit strains are strongly affected when the plastic anisotropy coefficient r 45 is > 1.0 whereas the same FLC's are obtained for the isotropic (r 45 = 1.0) and r 45 < 1.0 cases. Moreover, it should be noted that these material parameters provide equal in-plane principal stress-ratio in the equibiaxial stretching stress-state condition, namely, 2 1 1.0 β = σ σ = . These differences are attributed to the geometrical imperfection orientation Ψ which minimizes the computed limit strains, as shown in Figure 13b in the plot of the corresponding critical major principal strains, , 1 a * ε , as a function of the orientation Ψ. The minimum limit strain for the r 45 > 1,0 case is obtained for Ψ ≅ 45 degrees whereas ther 45 < 1.0 material has minima coincident with isotropic reference case for Ψ = 0 and 90 degrees. As pointed out by Barata da Rocha et al. 38 , the original assumption of normality of the geometrical imperfection with respect to major in-plane principal stress 19 is too restrictive since it would provide a unique limit strain in the equibiaxial stretching condition independently of the r 45 value. Figure 14 compares the theoretical predictions for an isotropic (von Mises) material determined from linear strain-paths and bi-linear strain-paths after a pre-straining level of 15% under uniaxial tension (UT), plane-strain tension (PS) and equibiaxial straining (ES) conditions, respectively. In comparison to the linear strain-paths predicted values, it is possible to observe that the limit strains increase in the biaxial strain domain after the same amount of pre-strain under UT and PS conditions. In the deep-drawing range, that is, for negative minor principal strains, the resulting limit strains from the second deformation stage decreased after a UT pre-straining and are improved after a pre-straining under PS condition. Conversely, the influence of the ES pre-straining is more pronounced given that a remarkable decrease of the limit strains in biaxial strain range and an increase in the Table 4. AISI-1012 steel sheet material properties 42 . negative minor principal strain domain. Except for the PS case, the plane-strain intercept of the FLC moved towards the pre-straining mode direction, namely, shifted-up after the UT pre-straining and drastically decreased after the ES pre-straining stage. It should be pointed out that the above bi-linear strain-path calculations were performed without unloading between the first and second stages. Figures 15 and 16 compare the predicted FLCs with the experimental results determined for an AISI-1012 steel sheet 42 from linear and bi-linear strain paths, respectively. The parameters adopted in Ferron's plasticity yield criterion were m = n = 2, p = q = 1 with k = 0.15. The parameters a and b which describe the in-plane sheet metal plastic anisotropy in Ferron's criterion, see Equation 11, were determined from the plastic anisotropic coefficients r 0 and r 90 . The linear strain-path predictions obtained with both Hill's and Ferron's yield criteria are in good agreement with the experimental data near the plane-strain intercept of the AISI-1012 FLC whereas either anisotropic yield criteria underpredict the limits strains for larger values of the inplane negative minor principal strain, as shown in Figure 15. In the biaxial stretching domain, a better approximation is achieved with Ferron's anisotropic yield criterion predictions. This is attributed to the flattening of Ferron's yield surface near to plane-strain tension/compression and pure shear stresses states, thanks to a positive k-value in Equation 11, providing a smaller ratio between the major principal plane-strain yield stress, σ PS1 , and the equibiaxial yield stress, σ b , namely, the material parameter P = σ PS1 /σ b introduced by Barlat 8 , in comparison with Hill's quadratic plane-stress yield locus. Figure 16a compares the theoretical predictions with the experimental AISI-1012 limit strains obtained after 10% of pre-strain in uniaxial tension. As for the linear strain-path, either yield criteria provided a good forecast near to the shifted plane-strain intercept. However, the corresponding theoretical predictions overestimate and underestimate the experimental data for the levels of the minor principal strain -0.20 ≤ ε 2 ≤ -0.10 and -0.40 ≤ ε 2 ≤ -0.30, respectively. Figure 16b shows that either yield criteria predictions underestimate the experimental limit strains after 8% of pre-strain in equibiaxial tension. Recently, Nurcheshmeh and Green 43 developed a M-K model based upon a mixed isotropic-nonlinear kinematic hardening along with Hill's quadratic yield description obtaining improved predictions for the cited AISI-1012 steel sheet limit strains under bi-linear strain-paths.   quadratic and Ferron's criteria, parameters a and b in Equation 11, were determined from the uniaxial tensile data assuming firstly B = 3A and, then, by exchanging the r 0 and r 90 experimental values. Also, the work-hardening law was transformed according to the Equation 36 by setting α = π/2 in Equation 35 since the effective stress measure in Ferron's orthotropic yield criterion is identified by the equibiaxial yield stress (σ 1 = σ 2 = σ b ). Both yield criteria descriptions provided a very good agreement with the experimental limit strains located in the drawing range and a good match with the empirical FLC 0 . This can be attributed to the description of the geometrical imperfection parameter of the M-K model based upon the experimental evolution of the surface roughness as a function of the mean initial grain size and the effective plastic-strain. On the other hand, only Ferron's orthotropic plasticity yield criterion provided a good prediction of the experimental limit strains in the biaxial stretching range. This is accredited to the yield surface shape in the region of interest, namely, the existence of a flattening in Ferron's yield locus between the equibiaxial and plane-strain tension stress-states, which, in turn, provide a decrease in the stress-ratio P = σ PS1 /σ b in comparison to the Hill's quadratic yield surface. Thus, for a given strain-path ρ = ε 2 /ε 1 or else a stress-ratio β = σ 2 /σ 1 , the strain localization in the imperfection zone of the M-K model is favored as the parameter P decreases. This yield surface shape effect is better explained in Figure 18, where the limit strains predictions of the adopted anisotropic yield criteria are plotted in the biaxial stretching range with the computed history strain-paths in both homogeneous (circle-dot lines) and imperfection (dashed lines) zones of the M-K localization model. First, one can observe the strain-path changes towards a plane-strain state in the imperfection zone and, then, the earlier plastic localization with Ferron's yield criterion. Figure 17 compares the experimental limit strains determined for the IF hot-dip galvanized steel sheet with the theoretical predictions of the M-K model based upon the Hill's quadratic 34 and Ferron's 35 yield criteria. From Equation 25 with the initial average roughness value of 0.93 mµ, the initial value of the geometrical imperfection of the M-K model f 0 is equal to 0.997. The parameters adopted in Ferron's orthotropic plasticity yield criterion were m = n = p = q = 2 and k = 0.2. It should be recalled that the experimental limit strains are defined according to the ISO 12004 -2:2008 29 as well as that the plane-strain intercept value has been predicted according to Equation 46. Given that the Marciniak punch tests were performed with the specimen blank length perpendicular to the rolling direction, the material anisotropic parameters of Hill's  the assumption of isotropic work-hardening forecasted only the experimental trends of the low carbon steel FLCs determined from bi-linear strain-paths after uniaxial and biaxial tension pre-straining. The M-K geometrical imperfection parameter was defined more properly from the experimental correlation of the surface roughness changes described as a function of the accumulated effective plastic strain and the mean grain size of the investigated hot-dip galvanized interstitial-free steel sheet. In the deep-drawing range, both Hill's quadratic and Ferron's yield criteria provided a good agreement with the limit strains of the hot-dip galvanized interstitial-free steel sheet. Only Ferron's description forecasted the experimental data in the biaxial stretching domain owing to the effect of the yield surface flattening near to the plane-strain tension stress-state.

Conclusions
In the present work, experimental and theoretical analyses of the limit strains of a hot-dip galvanized interstitial-free steel sheet were performed with the help of the Marciniak punch test technique and based on the M-K localization model, respectively.
The investigated hot-dip galvanized interstitial-free steel sheet presented a very good formability level, which is suitable for exposed parts in which an extra deep-drawing quality is often required. However, the adopted experimental procedure to conduct the FLC tests revealed to be too susceptible to the lubrication conditions between the flatbottomed punch and the carrier-blank. Besides, the limit strains defined as per ISO 12004-2:2008 has proven to be more consistent vis-à-vis the ASTM E 2218 -02 standard procedures. The proposed M-K model based upon both Hill's quadratic and Ferron's anisotropic yield criteria was first validated by means of comparisons with available data of a low carbon steel (AISI 1012) for linear and bi-linear strain-paths. Under linear strain-paths, reasonable agreement was obtained from Ferron's orthotropic yield criterion predictions. Conversely, either adopted yield criteria under