Theoretical predictions of dynamic necking formability of ductile metallic sheets with evolving plastic anisotropy and tension-compression asymmetry

In this paper, we have investigated necking formability of anisotropic and tension-compression asymmetric metallic sheets subjected to in-plane loading paths ranging from plane strain tension to near equibiaxial tension. For that purpose, we have used three different approaches: a linear stability analysis, a nonlinear two-zone model and unit-cell finite element calculations. We have considered three materials –AZ31-Mg alloy, high purity α-titanium and OFHC copper– whose mechanical behavior is described with an elastic-plastic constitutive model with yielding defined by the CPB06 criterion (15) which includes specific features to account for the evolution of plastic orthotropy and strength differential effect with accumulated plastic deformation (48). From a methodological standpoint, the main novelty of this paper with respect to the recent work of N’souglo et al. (42) –which investigated materials with yielding described by the orthotropic criterion of Hill (24)– is the extension of both stability analysis and nonlinear two-zone model to consider anisotropic and tension-compression asymmetric materials with distortional hardening. The results obtained with the stability analysis and the nonlinear two-zone model show reasonable qualitative and quantitative agreement with forming limit diagrams calculated with the finite element simulations, for the three materials considered, and for a wide range of loading rates varying from quasi-static loading up to 40000 s− 1, which makes apparent the capacity of the theoretical models to capture the mechanisms which control necking formability of metallic materials with complex plastic behavior. Special mention deserves the nonlinear two-zone model, as it does not need prior calibration –unlike the stability analysis– and it yields accurate predictions that rarely deviate more than 10% from the results obtained with the unit-cell calculations.


Introduction
Formability of sheet metal products generally increases at high strain rates [5,6,21,54,59]. The stabilizing effect of inertia, which delays necking formation and enhances workpiece ductility [20,29,55], is exploited in high velocity metal working operations, such as explosive, electromagnetic and electrohydraulic forming processes, to deform low-ductility lightweight alloys beyond their quasi-static limit [40].
For instance, Balanethiram and Daehn [6] performed electrohydraulic forming experiments with 6061 aluminum sheets at a speed of ≈ 300 m/s. The strain rate attained during the experiments was estimated to be above 1000 s − 1 . The specimens of thickness 1.6 mm were launched onto a conical die with an apex angle of 90 ∘ and a diameter base of 102 mm. The comparison of the electrohydraulic testing results with quasi-static forming limit diagrams taken from the literature [3] showed that the ductility of the metal sheets near plane strain increased by a factor of 4 at high strain rates. Years later, Golovashchenko [22] carried out pulsed electromagnetic forming of AA 6111-T4 and AA 5754 sheets into conical and v-shape dies, and compared the results with forming limit diagrams obtained under quasistatic loading by forming the aluminum specimens with a hemispherical punch into a round die cavity. Consistent with the earlier work of Balanethiram and Daehn [6], the results of Golovashchenko [22] showed that the limit strains for both aluminum alloys under nearly plane strain conditions were up to ≈ 2.5 times greater at high loading rates. Shortly after, Dariani et al. [18] obtained dynamic forming limit diagrams for Al 6061-T6 sheets by explosive forming of specimens of 1 mm thickness into cylindrical and flat type dies. Different specimen designs were used to obtain limit strains corresponding to both the negative and positive minor strain regions of the forming limit diagram. As in the experiments of Balanethiram and Daehn [6], the strain rate in the workpiece during the tests was estimated to be around 1000 s − 1 . Dariani et al. [18] reported an improvement of ≈ 150% in the formability of the material near plane strain in comparison with conventional quasi-static forming, with this increase being more/less important as the loading path moves away from plane strain and approaches equibiaxial/ uniaxial tension. Moreover, Rohatgi et al. [53] investigated the high strain rate formability of 5182-O aluminum sheets using the electrohydraulic forming technique. Free-forming and conical-die forming experiments were performed with specimens of 1 mm thickness tested at strain rates of ≈ 4000 s − 1 . The experimental results showed that the dynamic formability of the aluminum 5182-O sheets, at minor strains of − 0.1 and − 0.05, was ≈ 2.5 and ≈ 6.5 times greater than under quasi-static loading conditions for free-forming and conical-die forming conditions, respectively. The higher increase in ductility for conical-die forming can be explained by the fact that in this case an additional deformation occurs during the impact of the blank against the die [23].
Xu et al. [60,62] carried out electromagnetic free-bulging forming experiments with AZ31 magnesium sheets of 1 mm thickness. The tests were performed with and without an aluminum driver sheet. The driver served to increase the forming velocity, which was limited in the experiments without the driver due to the low electrical conductivity of the magnesium alloy. Specifically, the experiments without driver sheet showed that relative to the formability results obtained under quasi-static loading, the major and minor limit strains increased by 67% and 77%, respectively. On the other hand, the driver sheet was shown to have a positive impact on the formality of the AZ31 specimens, which increased as the thickness of the driver increased. For instance, using a driver sheet of 2 mm thickness led to an increase of the major and minor limit strains, relative to quasi-static loading, of 148% and 184%, respectively. Shortly after, the same authors carried out electromagnetic forming experiments with AZ31 magnesium sheets subjected to uniaxial tension [61]. The thickness of the samples was 1 mm, as in the experiments of Xu et al. [60,62]. The tests were performed at different forming speeds, with the maximum strain rate attained in the experiments varying between 1200 s − 1 and 3500 s − 1 . As in the free-bulging tests of Xu et al. [60,62], in which the samples were subjected to nearly equibiaxial stretching, the tensile formability of the AZ31 sheets in the electromagnetic experiments of Xu et al. [61] was significantly greater than in conventional quasi-static tests, the increase in the major and minor limit strains being 112% and 96%, respectively.
Electromagnetic free-bulging forming experiments on Ti-6Al-4V sheets of 1 mm thickness were performed by Li et al. [36] using a driver plate made of AA5052-O with 2 mm thickness. As in the formability tests performed by Xu et al. [60,62] with AZ31 magnesium specimens, the driver plate was used in order to reach high forming velocities, up to 165 m/s, that otherwise would be difficult to attain due to the low conductivity of the titanium alloy. Postmortem analysis of the specimens showed that the major and minor limit strains were 0.093 and 0.028, respectively. The comparison with quasistatic tests in which the titanium sheets were formed with a hemispherical punch into a round die cavity, for which the major and minor limit strains were 0.079 and 0.029, revealed that the ductility of the titanium sheets increased approximately 20% at high strain rates. Shortly after, the same authors developed a novel experimental design to perform electromagnetic forming experiments on Ti-6Al-4V sheets under uniaxial tension, plain strain and equibiaxial stretching, using different specimen designs [35]. The strain rate attained in the experiments was estimated to be around 4000 s − 1 . The tests of Li et al. [35] showed that, relative to the quasi-static formability of the material, the dynamic ductility of the Ti-6Al-4V sheets increased ≈ 15% for uniaxial tension, ≈ 75% for plain strain and ≈ 25% for equibiaxial stretching.
High velocity forming experiments with metallic materials other than lightweight alloys were also performed by several authors to determine the dynamic ductility of different structural materials. For instance, Balanethiramand Daehn [5] used electrohydraulic forming to obtain the limit strains of Interstitial Free iron sheets of 0.84 mm thickness deformed at strain rate of 1000 s − 1 under nearly plane strain conditions. The comparison with the corresponding quasistatic forming limits showed that the ductility of the material increased by nearly three times at high strain rates. Shortly after, Balanethiram and Daehn [6] carried out additional electrohydraulic forming tests with OFHC copper sheets of 0.34 mm and 0.79 mm thickness. The experimental evidence was that the thicker samples show greater improvement in formability at high strain rates probably due to the increase of inertia effects with the sheet thickness. Moreover, Seth et al. [54] investigated the high velocity formability of several high strength sheet steels with both very low and high quasi-static ductility. The samples were electromagnetically launched at a shaped punch at velocities of 50 − 220 m/s. The thickness of the steel sheets ranged from 0.15 mm to 0.58 mm, and different punch geometries were used to obtain different strain states in the specimens. An aluminum 6111 T4 sheet of the same diameter as the workpiece (80 mm) and 1 mm thick was used as a driver in all the experiments. The electromagnetic forming experiments yielded failure strains for all the steel sheets that were significantly beyond the quasi-static limit, for any strain state between plane strain and equibiaxial stretching, and notably for the materials with lower quasi-static ductility. Moreover, Golovashchenko et al. [23] performed electrohydraulic forming of DP500, DP590, DP780 and DP980 dual phase steel sheets into conical and a v-shape dies so that the samples were subjected to plane strain and equibiaxial stretching, respectively. The thickness of the DP500 specimens was 0.65 mm, while for the samples made of DP590, DP780 and DP980 it was 1 mm. The maximum strain rate attained in the experiments was estimated to be about 20000 s − 1 . The limit strains obtained in the electrohydraulic forming tests were compared with the formability limits resulting from quasi-static limiting dome height tests in which the samples were formed with a hemispherical dome into a round cavity. As in the electrohydraulic forming experiments, different specimen designs were used to obtain plane strain and equibiaxial stretching conditions in the samples during testing. Comparison of the limit strains measured in the electrohydraulic forming experiments and in the limiting dome height tests showed that the ductility of all materials tested increased at high strain rates. For instance, the improvement in plane strain formability varied between 63% and 190%, depending on the grade of dual phase steel.
All these experimental results substantiate the pioneering claim of Balanethiram and Daehn [6] that the beneficial contribution of inertial forces to the formability of metals sheets is both large and general, as it affects most ductile metallic materials. On the other hand, the experimental evidence is that the quantitative effect of inertia on formability improvement depends on the strain state [35,[60][61][62], and also on the mechanical behavior of the material, since the extent of the formability improvement is different for different metal sheets [23,54]. Nevertheless, as of today, it is not clear which features of the constitutive behavior of the material affect the most to the dynamic ductility of the metal sheets. In fact, the number of publications providing insights, either theoretically or using finite element calculations, on the effect of the constitutive behavior of the material on the dynamic formability of metal sheets is relatively small, and generally consider the material to be isotropic [29,44,53,58,63], which is a crude assumption since most sheet metal products display plastic anisotropy.
An exception is the recent work of N'souglo et al. [42], who developed a linear stability analysis and a nonlinear two-zone model, and performed finite element simulations for materials with yielding described by Hill [24] criterion, in order to investigate the role of plastic orthotropy on the dynamic formability of thin sheets of 2 mm thickness subjected to loading paths ranging from plane strain to near equibiaxial stretching, and to strain rates varying between 100 s − 1 and 50000 s − 1 . Forming limit diagrams obtained with the stability analysis and the two-zone model were compared with the finite element results obtained for five different materials, two of them were model materials with mechanical properties specifically tailored to bring to light the interplay between inertia and orthotropy in dynamic formability, while the other three were actual materials, TRIP-780 steel, AA 5182-O and AA 6016-T4. The results obtained with the three approaches -stability analysis, twozone model and finite element calculations-showed qualitative agreement for the five materials considered, and brought out that: (1) the effect of anisotropy on dynamic formability increases as the loading path moves away from plane strain and approaches equibiaxial tension, (2) the ductility improvement due to inertia becomes significant when the applied strain rate is greater than ≈ 2000 s − 1 , and (3) for strain rates greater than 10000 s − 1 inertia seems to be the main factor responsible for neck retardation, over and above the stabilizing effect of strain rate sensitivity. However, the scope of these findings is limited by the fact that N'souglo et al. [42] did not provide a quantitative comparison with experimental data, and also because the Hill yield function [24] fails when it comes to describe the mechanical response of metal sheets which exhibit complex plastic anisotropy [7,26]. In particular, the theory of Hill [24] has been shown not to reproduce the directionality of plastic properties of face-centered cubic metals, such as aluminium alloys, which display greater biaxial flow stress than uniaxial flow stress [8,47]. In addition, the formulation of Hill [24] does not account for the strength-differential effect which characterizes the mechanical response of hexagonal close packed metals, such as titanium and magnesium materials, which generally display both evolving plastic anisotropy and tension-compression asymmetry in yielding [14,15,31,50]. Therefore, it is evident that the work of N'souglo et al. [42] needs to be extended to consider more advanced constitutive models that capture the complex mechanical response generally shown by metal sheets frequently used in high-speed forming operations and needs to be validated with experimental data (see previous paragraphs of this introduction).
Many yield functions have been developed over the years to improve the constitutive model of Hill [24]. For instance, Karafillis and Boyce [30] derived a non-quadratic yield criterion to describe different states of material symmetry, including the case of the fully asymmetric material.
The key feature of the model of Karafillis and Boyce [30] is a fourth order tensor employed as a linear multiplicative operator acting on the stress tensor to introduce material anisotropy. Following the work of Karafillis and Boyce [30], the linear transformation-based anisotropic yield functions became increasingly popular because they allow tuning the number of anisotropy parameters to increase the accuracy of the yield criterion predictions, without this having an impact on the mathematical structure of the constitutive model. For instance, Barlat et al. [10] formulated a plane stress yield function with 9 parameters to describe the mechanical response of aluminum sheets using two linear transformations of the stress tensor. Bron and Besson [12] extended the models of Barlat et al. [10] and Karafillis and Boyce [30] using two linear transformations to obtain a 3D yield function with 16 parameters which provided accurate description of the mechanical behavior of different aluminum sheets. Two linear transformations were also used by Barlat et al. [9] to formulate the so-called Yld2004-13p and Yld2004-18p yield functions, that were derived from different isotropic yield functions, and contain 13 and 18 anisotropic parameters, respectively. Two and three linear transformations were used by Aretz and Barlat [4] to obtain the socalled Yld2011-18p and Yld2011-27p yield functions, which are defined by 18 and 27 parameters, respectively. Moreover, Cazacu and Barlat [14] formulated a yield criterion to describe both the asymmetry and anisotropy in yielding of magnesium and magnesium alloys using the generalized invariants approach developed by Cazacu and Barlat [13]. The yield function of Cazacu and Barlat [14] is a homogeneous function of degree three in stresses and includes 18 material parameters. Additional yield criteria to model the anisotropy and strength differential effect of various metallic materials have been proposed, for instance, by Kondori et al. [32], Lee et al. [34], Nixon et al. [41], Park et al. [45] and Raemy et al. [49].
In this paper we have developed a linear stability analysis and a nonlinear two-zone model to construct forming limit diagrams for elastic-plastic materials with yielding described by the CPB06 criterion [15]. This constitutive model accounts for both plastic anisotropy and yielding asymmetry, and describes the change of the shape of the yield surface during monotonic loading through the evolution with the plastic deformation of the material parameters involved in the expression of the yield function. The predictions obtained with the analytical models have been compared with finite element simulations and with experimental data obtained from the literature for magnesium, titanium and copper [6,36,60,62], and reasonable agreement has been obtained for a wide range of loading rates varying from quasi-static loading to strain rates up to 40000 s − 1 , and for loading paths ranging from plane strain to near equibiaxial stretching. To the authors' knowledge, there are no other publications in the open literature with analytical models predicting the dynamic formability of metal sheets showing evolving plastic anisotropy and tension-compression asymmetry.

Constitutive framework
The mechanical behavior of the materials investigated is described with an elastic-plastic constitutive model with yielding defined by the orthotropic and asymmetric CPB06 criterion [15] which includes specific features to consider the evolution of plastic anisotropy and strength differential effect with accumulated plastic deformation [48]. The formulation of the constitutive model is outlined in Section "Formulation", and the values of the constitutive parameters of the materials investigated are provided in Section "Materials investigated".

Formulation
The rate of deformation tensor d is the sum of an elastic part d e and a plastic part d p : where the relationship between the elastic part and the rate of the stress is: where ▿ is an objective derivative of the Cauchy stress tensor, and C is the fourth-order isotropic elastic tensor defined as: with 1 being the unit second-order tensor and ′ the unit deviatoric fourth-order tensor. Moreover G and K are the shear and bulk moduli, respectively.
Assuming an associated plastic flow rule, the plastic part of the rate of deformation tensor is defined as: where ̇ is the rate of the plastic multiplier and ̄ is the effective stress [15] given by: where k is the tension-compression asymmetry parameter, see Eq. 10a, and Σ i (i = 1,...,3) are the principal values of the transformed stress tensor Σ defined as: with L being an orthotropic, deviatoric, and symmetric fourth-order tensor, and s the deviatoric part of the Cauchy stress tensor σ. In the Cartesian coordinate system (X 1 ,X 2 ,X 3 ) associated to the orthotropy axes of the material, the deviatoric part of the Cauchy stress tensor is represented by a 6-dimensional vector s = (s 11 ,s 22 ,s 33 ,s 12 ,s 23 ,s 13 ) and the fourthorder orthotropic tensor L is represented by a 6 × 6 matrix: with the components of L being the anisotropy coefficients, see Eq. 10b. Moreover, the parameter m in Eq. 5 is defined as: with Φ i (i = 1,...,3) expressed as follows: The evolution of plastic anisotropy and tension-compression asymmetry with accumulated plastic strain is modeled using the methodology developed by Plunkett et al. [48]. Firstly, the values of the parameters k and L ij (i,j = 1,...,6), see Eqs. 5 and 7, are identified for selected levels of effective plastic strain ̄p ,1 <̄p ,2 < ... <̄p ,m . Secondly, a linear interpolation is used to obtain the values of the tension-compression and anisotropy coefficients corresponding to any given level of effective plastic strain ̄p ,n <̄p <̄p ,n+1 (n = 1,...,m − 1): where the interpolation parameter α is: such that (̄p ,n ) = 1 and ̄p ,n+1 = 0.
The yield function is defined as: where σ Y is the tensile yield stress in the rolling direction of the material X 1 , given in the present work by the following relationship: with A 0 , A 1 and A 2 being material parameters. Notice that viscous and thermal effects on the yield stress evolution are neglected. Moreover, the effective plastic strain is ̄p = ∫ t 0̇̄p d , where t refers to time, and ̇̄p is the effective plastic strain rate given by: 2 3 , and b 1 , b 2 and b 3 are the principal values of the transformed rate of deformation tensor defined as: where ∶ = � with = ∶ � . Moreover, the work conjugacy relation: yields: The formulation of the constitutive behaviour is completed with the Kuhn-Tucker loading-unloading conditions: and the consistency condition during plastic loading:

Materials investigated
We investigate 3 materials widely used in different industrial sectors and engineering applications in the form of plates and shells: AZ31-Mg alloy, high purity α-titanium and OFHC copper. Sheet products of these materials are often subjected to dynamic forming operations to produce complex metal parts, see the introduction to this paper (Section "Introduction"). In addition, the three materials display very different yielding behaviour, see Fig. 1, increasing the scope of the analysis carried out in Section "Analysis and results". Table 1 shows the values of the initial density, the elastic constants and the parameters of the yield stress in the rolling direction [50]. We are aware that using a Voce-type hardening relation Eq. 13 to describe the yield stress of the materials investigated may be a crude assumption, particularly at high strain rates, when viscous and thermal effects -due to adiabatic heating-are known to affect dynamic formability limits [27,29,42]. Nevertheless, using a simple constitutive relation reduces the number of material parameters, which helps in the identification of the effects of anisotropy and tension-compression asymmetry on dynamic formability.
The values of the anisotropy coefficients L ij (i,j = 1,...,6) and the tension-compression asymmetry parameter k for the initial yielding and for several individual levels of effective plastic strain are given in Table 2. Data are taken from Revil-Baudard et al. [50] and Cazacu et al. [16]. Recall from Section "Formulation" that the values of L ij and k for effective plastic strains within the intervals reported in Table 2 are calculated using a linear interpolation function, see Eq. 10a. Moreover, in absence of additional data, we consider that L ij and k do not evolve for effective plastic strains outside the range of values included in Table 2, so that for larger strains the values of L ij and k for AZ31-Mg alloy, high purity α-titanium and OFHC copper are the same employed for ̄p = 0.1 , 0.3 and 0.12, respectively. Moreover, notice that, while OFHC copper is considered to be isotropic, high purity α-titanium and AZ31-Mg alloy display both evolving anisotropy and tension-compression asymmetry. Figure 1 pictures plane stress theoretical yield loci for biaxial loading conditions, σ I versus σ II , where σ I and σ II are the components of the Cauchy stress tensor in the rolling and transverse directions, respectively. Results are shown for the three materials investigated in this work and different levels of effective plastic strain. In the case of AZ31-Mg alloy, Fig. 1a, the yield locus evolves from an triangular shape for low values of effective plastic strain (the triangular shape is caused by very pronounced tension-compression asymmetry) to an elliptical shape at large strains. On the other hand, notice that the evolution of the yield locus with the plastic strain is different for high purity α-titanium, see Fig. 1b, such that yield locus is initially ellipsoidal and it becomes triangular at large strains (increasing strength differential effect with plastic straining). Moreover, notice that both AZ31-Mg alloy and high purity α-titanium display plastic anisotropy so that the yield stress in the rolling and the  Table 2 Numerical values of the anisotropy coefficients L ij and the tension-compression asymmetry parameter k for different values of the effective plastic strain ̄p for AZ31-Mg alloy, high purity α-titanium and OFHC copper. Data taken from Revil-Baudard et al. [50] and Cazacu et al. [16] Anisotropy and tension-compression asymmetry parameters  Fig. 1c, the predicted yield loci have an elliptical type shape, with the yield stress in compression being slightly greater than in tension (mild strength differential effect so that the yield loci do not become triangular). Recall that this material is taken to be isotropic and therefore the yield stress in the rolling and the transverse directions is the same. Appendix A provides a comparison with the plane stress yield loci calculated with von Mises [39] criterion, which illustrates the effect of anisotropy and tension-compression asymmetry on the mechanical behavior of AZ31-Mg alloy, high purity α-titanium and OFHC copper.

Problem statement
The problem addressed is that of a thin plate of initial thickness h 0 and edges of initial length L 0 , modeled with the constitutive framework described in Section "Constitutive framework", and subjected to constant and opposed stretching velocities on opposed sides, see Fig. 2. The applied velocities are V x = ±̇0 xx L 0 ∕2 and V y = ±̇0 yy L 0 ∕2 , with ̇0 xx and ̇0 yy being the imposed initial strain rates. The loading condition is determined by the ratio =̇0 yy ∕̇0 xx , which is varied between 0 (plane strain stretching) and ≈ 1 (near equibiaxial stretching) in the calculations presented in this paper. Hereinafter, ̇0 xx and ̇0 yy will be referred to as imposed initial major and minor strain rates, respectively, and χ as loading path. The Lagrangian Cartesian coordinate system associated to the applied velocity field is denoted by (X, Y, Z) , while X 1 , X 2 , X 3 is a Lagrangian Cartesian frame defined by an angle ψ between X and X 1 , with X 1 , X 2 and X 3 being the rolling, the transverse and the through-thickness orthotropy directions of the material, respectively (recall that the OFHC copper is taken to be isotropic so that for this material there is no need to define material axes). Hereinafter the angle ψ is referred to as material orientation (for OFHC copper there is no need to define the angle ψ). We consider the cases for which the orthotropy axes and the loading axes are co-directional ( = 0 • or 90 • ) , so that the principal directions of stress and strain are parallel (this is always the case for OFHC copper). Note that the direction perpendicular to the major stretching direction X is the orientation naturally selected by the material to trigger a neck for the loading paths (0 ≤ χ ≤ 1) and material orientations ( = 0 • or 90 • ) investigated in this work. The origin of both Lagrangian coordinate systems is located at the center of mass of the plate. Note that the out-of-plane direction is assumed of plane stress. The problem is approached using a three-pronged methodology which includes a linear stability analysis, a nonlinear two-zone model and unit-cell finite element calculations. The main novelty is to extend the recent work of N'souglo et al. [42] -which investigated materials described with the orthotropic yield criterion of Hill [24]-to consider the orthotropic and tension-compression asymmetric yield criterion with distortional hardening described in Section "Constitutive framework". A brief note on the three approaches is provided next, while the reader is referred to the papers of N'souglo et al. [42] and Jacques [29] to obtain details on the basic features of the mathematical derivation of stability analysis and non-linear two-zone model, respectively, and to the papers of Rodríguez-Martínez et al. [52] and N'souglo et al. [42] for a complete description of the finite element model (with the difference that in this paper we consider the constitutive framework presented in Section "Constitutive framework"). In addition, in Supplementary material we provide a link to download all the codes developed for the implementation of the linear stability analysis and the nonlinear two-zone model in Wolfram Mathematica with specific instructions for reproducing the theoretical predictions reported in Section "Analysis and results".

Linear stability analysis
This technique is based on the superposition of a small perturbation on the homogeneous solution of the problem, so that if the perturbation grows faster than the background solution, the plastic flow is unstable and a non-homogeneous, neck-like deformation field can develop. We have used the same 2D approach employed in [43,51,63] so that the stress multiaxiality effects that develop inside the necked section have been approximated using Bridgman [11] correction. The fundamental solution of the problem is composed of 16 equations which are nondimensionalized, perturbed using the frozen coefficients method and linearized (note that in the case of N'souglo et al. [42] the number of equations was 17 because the energy balance was taken into account). Imposing that the determinant of the matrix of coefficients of the resulting system of dimensionless linear algebraic equations is equal to zero leads to a third order polynomial for the instantaneous growth rate of the perturbation η which depends on the fundamental solution and on Fig. 2 Schematic representation of the geometry and boundary conditions of the problem addressed: a thin plate of initial thickness h 0 and edges of initial length L 0 subjected to constant and opposed stretching velocities V x = ±̇0 xx L 0 ∕2 and V y = ±̇0 yy L 0 ∕2 on opposed sides. The Lagrangian Cartesian coordinate system associated to the applied velocity field is denoted by (X, Y, Z) , while X 1 , X 2 , X 3 is a Lagrangian Cartesian frame defined by an angle ψ between X and X 1 , with X 1 , X 2 and X 3 being the rolling, the transverse and the through-thickness orthotropy directions of the material, respectively. The origin of both Lagrangian coordinate systems is located at the center of mass of the plate. The out-of-plane direction is assumed of plane stress the perturbation wavenumber ξ. The root of the polynomial with physical meaning is the one with the greatest positive real part [19], hereinafter denoted by η + . The wavenumber is related to the normalized perturbation wavelength L = L 0 ∕h 0 as L = 2 −1 h 0 . Hereinafter, L will be also referred to as necking wavelength. On the other hand, the instantaneous growth rate is used to compute the cumulative instability index Ī = ∫ t 0 + d which tracks the history of the instantaneous grow rate of all the growing modes. Assuming that a perturbation mode turns into a necking mode when the cumulative instability index reaches a critical value which depends on the loading path and the strain rate, the linear stability analysis is used to construct forming limit diagrams which are compared with nonlinear two-zone model predictions and unit-cell finite element simulations. The procedure for calibration of the linear stability analysis is detailed in Appendix B, and the influence of the calibration procedure in the stability analysis predictions is shown in Appendix C.

Nonlinear two-zone model
The non-linear two-zone model [29] is an extension of the classical analysis of localized necking proposed by Marciniak and Kuczyński [37] to account for inertia effects. The model considers an imperfection in the form of a band of reduced thickness, perpendicular to the main straining direction, at the center of a unit-cell of length L 0 (see Fig. 4 in N'souglo et al. [42]). Note that L 0 also corresponds to the neck spacing. The initial thickness in the imperfection zone is h A,0 = h 0 (1 −Δ), Δ being the normalized imperfection amplitude. The ratio between the initial length of the imperfection L A,0 x and the initial unit-cell length is denoted by R = L A,0 x ∕L 0 . In the present work, the value of this parameter was set to R = 0.28. This value was identified in the case of isotropic material [29] and was found to remain appropriate for orthotropic materials obeying the Hill [24] yield criterion [42]. The methodology to obtain necking strains with the two-zone model is described in Section 2.5 of Jacques [29]. Computations have been carried out for several normalized cell sizes, 0.5 ≤L ≤ 6 with L = L 0 h 0 and h 0 = 2 mm, in order to identify the critical neck spacing and the corresponding necking strains (critical necking strains). The results presented in Section "Analysis and results" have been obtained with a normalized imperfection amplitude Δ = 0.2%. This value has been chosen arbitrarily. It would have been probably possible to improve the agreement with the experimental results (see Section "The influence of material behavior") by adjusting the value of Δ for each material investigated. However, in order to illustrate the influence of the material behavior on formability, we have preferred to keep the same value of Δ for all calculations.

Unit-cell finite element calculations
The simulations consider a unit cell with a sinusoidal spatial imperfection defined as h = h 0 − 1 + cos 2 X L 0 , where δ is the amplitude of the imperfection. Exploiting the symmetry of the model, only one eight of the plate has been analyzed, with reference configuration (imperfection-free) defined by the domain Recall that the origin of coordinates of the Lagrangian Cartesian coordinate system (X, Y, Z) is located at the center of mass of the whole cell. As in the two-zone model calculations, we consider several normalized cell sizes 0.5 ≤L ≤ 6 with L = L 0 h 0 and h 0 = 2 mm fixed, and the normalized imperfection amplitude is Δ = 2 h 0 = 0.2% . Note that the normalized cell size L will be also called necking wavelength, to be consistent with the notation employed in the stability analysis and the two-zone model. The finite element calculations are performed with the finite element software ABAQUS under the initial and boundary conditions defined in equations (45) and (46) of N'souglo et al. [42], which are consistent with the stability analysis and the twozone model. The unit-cell calculations have been carried out for selected loading paths χ = 0, 0.125, 0.25, 0.375, 0.5, 0.625 and 0.75. The major stretching direction coincides with the X axis. The finite element model is meshed with eight-node solid elements, with reduced integration and hourglass control (C3D8R in ABAQUS notation), with dimensions ≈ 25 × 25 × 50 μm 3 (elements are initially shorter along the loading directions). The constitutive model presented in Section "Constitutive framework" has been implemented in ABAQUS/Standard [2] and ABAQUS/ Explicit [1] via UMAT and VUMAT user subroutines, respectively, using the stress-update algorithm based on the numerical approximation of the yield function gradients developed by Hosseini and Rodríguez-Martínez [25]. The UMAT [2] has been used to obtain the low strain rate finite element results, 0.0001 s − 1 and 100 s − 1 , reported in Figs. 3 and 13, and the VUMAT [1] has been used to obtain all the other finite element calculations reported in this paper for higher strain rates. Note that both ABAQUS/Standard and Explicit simulations take inertia effets into account. The difference is that the time integration of the equations of motion is based on the central difference scheme in ABAQUS/ Explicit [1] and on the Hilbert-Hughes-Taylor (HHT) method in ABAQUS/Standard [2]. The use of the implicit HHT scheme leads to smaller computation times for low or moderate loading rates, while the explicit scheme is more efficient for highly dynamic loading conditions. Note also that ABAQUS/Explicit uses the Green-Naghdi rate to compute the objective derivative of the stress (see Eq. 2), while the stress rate in ABAQUS/Standard corresponds to the Jaumann derivative. Comparison with calculations in which the material behavior is modeled with von Mises plasticity (i.e. imposing on the CPB06 criterion k = 1 and L ij = δ ij , with i,j = 1,...,6 and δ ij being the Kronecker unit delta) and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1). The experimental data corresponding to AZ31-Mg, high purity α-titanium and OFHC copper are taken from Choi et al. [17], Stachowicz [57] and Melander [38], respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Analysis and results
Calculations have been performed with the three approaches outlined in previous section (stability analysis, two-zone model and unit-cell finite element model) for loading paths varying from plane strain tension to near equibiaxial tension (0 ≤ χ ≤ 1, see Section "Problem statement") and different loading rates varying from quasi-static loading up to ̇0 xx = 40000 s −1 . Note that this value is larger than the strain rate generally attained in dynamic forming processes (e.g. behavior is modeled with von Mises plasticity (i.e. imposing on the CPB06 criterion k = 1 and L ij = δ ij , with i,j = 1,...,6 and δ ij being the Kronecker unit delta) and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1) see Dariani et al. [18], Golovashchenko et al. [23] and Li et al. [35], among others), but allows a further exploration of inertia effects. Recall that formulation of linear stability analysis and nonlinear two-zone model for materials with complex mechanical behavior including evolving anisotropy and tension-compression asymmetry is the main methodological contribution of this work. Results are shown for AZ31-Mg alloy, high purity α-titanium and OFHC copper, with material parameters reported in Tables 1 and 2. For AZ31-Mg alloy and high purity α-titanium, which unlike OFHC copper display plastic anisotropy, calculations are carried out for two material orientations ψ = 0 ∘ and ψ = 90 ∘ (see Section "Problem statement"). Note that all the results obtained in this section using the two-zone model and the stability analysis for the materials modeled with the CPB06 criterion are predictions, as the two-zone model does not need prior calibration (see Section "Problem statement"), and the stability analysis is calibrated using finite element simulations different from those presented below (see Appendix B). In Section "The influence of material behavior", quasi-static forming limit diagrams obtained for the three materials are systematically compared with experimental results available in the literature and with calculations in which the material behavior is modeled with von Mises plasticity (i.e. imposing on the CPB06 criterion k = 1 and L ij = δ ij , with i,j = 1,...,6 and δ ij being the Kronecker unit delta) and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1). In Section "The influence of loading rate", calculations performed for different loading rates are compared to identify the effect of inertia on the necking formability of the materials investigated. Recall that in Supplementary material we provide access to the codes developed for the implementation of the stability analysis and the two-zone model, so that the interested reader can reproduce the results reported in this paper. Figure 3 displays forming limit diagrams (critical major necking strain c xx versus critical minor necking strain c yy ) predicted by the stability analysis, the nonlinear two-zone model and the unit-cell finite element calculations. The critical major necking strain is the logarithmic strain in the major loading direction when necking occurs and it is determined as in N'souglo et al. [42], so that the reader is referred to Section 6.1 therein for the procedure to obtain c xx , and no additional details are provided here for the sake of brevity. The critical minor necking strain c yy is the logarithmic strain in the minor loading direction when c xx is measured. The results correspond to AZ31-Mg alloy with material orientations ψ = 0 ∘ and 90 ∘ , see Fig. 3a, high purity α-titanium with ψ = 0 ∘ and 90 ∘ , see Fig. 3b, and OFHC copper, see Fig. 3c. A comparison is performed with calculations in which the material behavior is modeled with von Mises plasticity (i.e. imposing on the CPB06 criterion k = 1 and L ij = δ ij , with i,j = 1,...,6 and δ ij being the Kronecker unit delta) and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1). The results are obtained for an imposed strain rate of ̇0 xx = 100 s −1 , that for all purposes can be considered quasi-static loading for the calculations performed in this work, as for this loading rate the effect of inertia on the necking formability of the three materials is negligible (for the specimen thickness considered), see Appendix D.

The influence of material behavior
The c xx − c yy curves display a concave-downward shape, with the formability increasing as the loading path moves away from plane strain tension. The results obtained with the stability analysis and the nonlinear two-zone model are in qualitative agreement with the finite element calculations -which are considered as the reference approach to validate the theoretical models-for the three materials investigated. For AZ31-Mg, see Fig. 3a, the formability for the von Mises material (black lines and markers) becomes greater than the results obtained with the CPB06 criterion as the loading path moves away from plane strain; the values of c xx being slightly larger for ψ = 0 ∘ (red lines and markers) than for ψ = 90 ∘ (green lines and markers). By analogy with the results of N'souglo et al. [42] for materials modeled with Hill [24] criterion, the relative order of the c xx − c yy curves obtained with von Mises and CPB06 models depends on the Lankford coefficient when the material is considered anisotropic, the general trend being that if the Lankford coefficient is greater than 1 (as it is the case of AZ31-Mg for ψ = 0 ∘ and 90 ∘ , such that r 0 • = 1.40 and r 90 • = 1.88 for ̄p ≥ 0.1 , see equation (17) in Cazacu et al. [15]) the von Mises material displays greater formability. In addition, Appendix A shows that the CPB06 criterion predicts greater curvature of the yield locus near equibiaxial tension than von Mises, which is assumed to decrease the material formability. The experimental formability results reported by Choi et al. [17] (orange markers) for AZ31B sheets obtained using the limiting dome height test (see Fig. 3b therein) are also included in Fig. 3a. Near plane strain tension, the experimental data are in good quantitative agreement with the finite element predictions, however, as the loading path moves away from plane strain, the agreement gets slightly worse. On the other hand, both stability analysis and two-zone model predictions overestimate the experimental results with the differences increasing as the loading path moves away from plane strain. Fig. 3b pictures the results corresponding to high purity α-titanium. As in the case of AZ31-Mg alloy, the calculations using von Mises criterion yield greater values for the critical major necking strain as the loading path moves away from plane strain tension, and the lower formability corresponds to the CPB06 criterion with orientation ψ = 90 ∘ . As mentioned before, this can be attributed to the fact that the anisotropic material displays a Lankford coefficient greater than 1 for ψ = 0 ∘ and 90 ∘ (namely, r 0 • = 1.09 and r 90 • = 1.11 for ̄p ≥ 0.3 ), and also because the curvature of the yield locus near equibiaxial tension is greater than for the von Mises material, see Appendix A. In Fig. 3b the experimental data reported in Fig. 7 of Stachowicz [57] for high purity α-titanium sheets are included. The finite element predictions are in very good agreement with the experiments. The two-zone model slightly underestimates the experimental data for all the loading paths and the linear stability analysis slightly underestimates/overestimates the experimental data for loading paths near plane strain/equibiaxial tension. We have also included in Fig. 3c the experiments reported in Fig. 10 of Melander [38] for OFHC copper sheets. The numerical simulations and analytical models show good overall qualitative agreement with the experimental results. Moreover, note that for OFHC copper both theoretical models and the finite element calculations show that the predictions using von Mises plasticity are above the results obtained for CPB06 as the loading path moves away from plane strain stretching and approaches equibiaxial tension. We have checked that taking k positive (see Table 2) the trend is reversed (results are not shown for the sake of brevity), and the critical major necking strain is smaller using von Mises plasticity, suggesting that for an imposed value of the initial major strain rate, the formability increases for a material that displays a larger yield stress in uniaxial tension than in uniaxial compression.
These results make apparent the influence of the yield criterion on the materials formability (for the same elastic constants and yield stress parameters), and the capacity of both stability analysis and two-zone model to provide predictions which are in good accord with experiments and consistent with finite element results. Figure 4 shows forming limit diagrams obtained with unitcell calculations, linear stability analysis and two-zone model for AZ31-Mg alloy and different initial major strain rates (5000 s − 1 , 10000 s − 1 , 20000 s − 1 and 40000 s − 1 ) for which the effect of inertia on the formability of the material is noticeable. The agreement between the predictions obtained with the analytical models and the finite element calculations is reasonably good for all loading paths and strain rates. Notice that, as the value of ̇0 xx increases, the c xx − c yy curves are shifted upwards. Moreover, as the loading path moves away from plane strain tension, the von Mises criterion predicts greater formability than the CPB06 model (as in the results shown in Fig. 3a).

The influence of loading rate
The evolution of the critical major necking strain with the imposed initial major strain rate for χ = 0 (plane strain tension) and χ = 0.5 is shown in Fig. 5. The c xx −̇0 xx curves display a concave-upward shape, so that the necking strain is roughly constant for strain rates lower than ≈ 1000 s − 1 (negligible inertia effects, see Appendix D), and then increases nonlinearly for greater strain rates. For plane strain tension, the unit cell finite Mises criterion, and 247% and 239% for the CPB06 criterion with orientations ψ = 0 ∘ and 90 ∘ , respectively. Similar quantitative results are obtained with the stability analysis and the two-zone model. For χ = 0.5, the c xx −̇0 xx curves are shifted upwards since necking is delayed as the loading path moves away from plane strain. The increase of the necking strain with the strain rate is Δ c xx | 40000 s −1 100 s −1 = 88% for the finite element calculations carried out with von Mises criterion, and 214% behavior is modeled with von Mises plasticity (i.e. imposing on the CPB06 criterion k = 1 and L ij = δ ij , with i,j = 1,...,6 and δ ij being the Kronecker unit delta) and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1) and 241% for the calculations performed with CPB06 criterion with ψ = 0 ∘ and 90 ∘ , respectively. This shows that the influence of strain rate on necking strains is strongly dependent on the yield criteria considered. The results presented in Fig. 5 indicate that the stabilizing effect of inertia at high strain rates is enhanced for materials having low formability. The increase of the necking strain with the strain rate computed with the CPB06 model is of the same order as in the experiments of Xu et al. [62] (see Section "Introduction"). Moreover, note that for von Mises the finite element calculations show that the stabilizing effect of inertia strongly decreases as the loading path moves away from plane strain (the same trend was obtained in the unit-cell calculations performed by Rodríguez-Martínez et al. [52] and N'souglo et al. [42] for materials modeled with von Mises [39] and Hill [24] criteria); however, this trend is less evident for the calculations performed with CPB06, showing that the effect of inertia on neck retardation depends on the yield criterion. Figure 6 shows the forming limit diagrams for high purity α-titanium and the same strain rates considered in the calculations performed for AZ31-Mg alloy, i.e. 5000 s − 1 , 10000 s − 1 , 20000 s − 1 and 40000 s − 1 . The relative order of the c xx − c yy curves is independent of the strain rate, with the von Mises material displaying greater formability, being the results obtained with the CPB06 criterion for ψ = 0 ∘ and 90 ∘ very similar (see Section "The influence of material behavior"). While both theoretical models predictions find good correlation with the finite element results, the stability analysis generally overestimates the limit strains obtained from the unit-cell calculations, and the opposite trend is obtained with the two-zone model.
The results for the evolution of the critical major necking strain with the imposed initial major strain rate shown in Fig. 7 for χ = 0 and 0.5 display the ability of both theoretical models to capture the effect of inertia on neck retardation for high purity α-titanium (differences between theoretical model predictions and finite element calculations being less than 10%). For plane strain tension, the increase of the necking strain with the strain rate is Δ c xx | 40000 s −1 100 s −1 = 245% for von Mises criterion, and 212% and 255% for the finite element calculations performed with CPB06 model with ψ = 0 ∘ and 90 ∘ , respectively. For χ = 0.5, the effect of inertia is significantly less, and the increases become 94% for von Mises, and 127% and 153% for CPB06 with ψ = 0 ∘ and 90 ∘ , respectively. These results reinforce the idea that the stabilizing effect of inertia depends on both the loading path and the yield criterion. Figure 8 shows forming limit diagrams for OFHC copper and the same values of the initial major strain rate considered in Figs. 4 and 6, i.e., 5000 s − 1 , 10000 s − 1 , 20000 s − 1 and 40000 s − 1 . Recall from Section "The influence of material behavior" that as the loading path moves away from plane strain the critical major necking strains corresponding to von Mises criterion are above the results obtained with Fig. 7 Critical major necking strain c xx versus imposed initial major strain rate ̇0 xx obtained with unit-cell calculations (FEM), linear stability analysis (LSA) and two-zone model (2ZM) for high purity α-titanium. The loading path is: (a) χ = 0 (plane strain stretching) and (b) χ = 0.5. Comparison with calculations in which the material behavior is modeled with von Mises plasticity (i.e. imposing on the CPB06 criterion k = 1 and L ij = δ ij , with i,j = 1,...,6 and δ ij being the Kronecker unit delta) and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1) CPB06 model because the yield stress in uniaxial compression is greater than in uniaxial tension (i.e. because k < 0, see Table 2). We have checked that, for the whole range of strain rates investigated, considering the same absolute values of the tension-compression parameter k, yet positive, yields critical major necking strains above the von Mises predictions. The reader is referred to the paper of N'souglo et al. [43] for specific analysis and discussion on the effect of the strength differential effect on the dynamic necking formability of ductile sheets subjected to plane strain stretching. The electrohydraulic forming experiments reported by Balanethiram and Daehn [6] for OFHC copper sheets tested behavior is modeled with von Mises plasticity (i.e. imposing on the CPB06 criterion k = 1 and L ij = δ ij , with i,j = 1,...,6 and δ ij being the Kronecker unit delta) and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1). The experimental data are taken from Balanethiram and Daehn [6] at strain rates that were estimated to be ≈ 500 s − 1 -this is a lower bound estimate obtained from a simple analytical model. Subsequently, more refined analyses of the deformation during conical-die electrohydraulic forming based on finite element computations revealed that strain rates larger than 10000 s − 1 are achieved in some parts of the specimen [28]-are compared in Fig. 8c with the results obtained using both analytical approaches and the finite element calculations. The agreement between experimental data, two-zone model and finite element calculations is satisfactory, however, the stability analysis underestimates the formability obtained with the finite elements near plane strain tension for 40000 s − 1 . The evolution of the critical major necking strain with the imposed initial major strain rate for χ = 0 and 0.5 is pictured in Fig. 9a and b, respectively. These graphs make apparent the underestimation by the stability analysis of the stabilizing effect of inertia at high strain rates. A reason for this disagreement is that the stability analysis has been calibrated using the physical constants and yield stress parameters corresponding to AZ31-Mg, which has a mass density ≈ 5 times smaller than copper (see Table 1), so that the dependence of the critical instability index with the strain rate underestimates the effect that inertia has on necking in copper (we have checked that calibrating the set of coefficients of the critical instability index using the initial density, elastic constants and parameters of the yield stress corresponding to OFHC copper, see Table 3, the stability analysis predicts greater stabilizing effect of inertia at high strain rates; see also the results in Fig. 12 where the calibration with OFHC copper yields slightly greater values of critical major necking strain). Moreover, the increase of the necking strain with the strain rate computed with the finite element calculations for plane strain tension is c xx | 40000 s −1 100 s −1 = 352% and 308% for the von Mises criterion and the CPB06 model, respectively. These values are similar to the experimental data reported by Balanethiram and Daehn [6] (see Table 1 therein) and they are greater than for AZ31-Mg alloy and high purity α-titanium (see the paragraphs above), which makes apparent that the contribution of inertia to neck retardation depends on the material. For χ = 0.5, the increase of the necking strain with the strain rate for both von Mises and CPB06 yield criteria is 126% and 191%; i.e., smaller than for plane strain tension.
Recall from Section "Materials investigated" that we have not considered viscous and thermal effects on the behavior of the materials. For specific insights into the role that thermal softening and strain rate sensitivity have on forming limit diagrams of anisotropic materials, the reader is referred to Section 6.4 of N'souglo et al. [42]. The general trend is that, while thermal softening promotes early necking localization, the strain rate sensitivity increases the formability of ductile metals. Furthermore, additional results are presented in Appendix E to shed light on the specific influence of tension-compression asymmetry on the formability. Fig. 9 Critical major necking strain c xx versus imposed initial major strain rate ̇0 xx obtained with unit-cell calculations (FEM), linear stability analysis (LSA) and two-zone model (2ZM) for OFHC copper. The loading path is: (a) χ = 0 (plane strain stretching) and (b) χ = 0.5. Comparison with calculations in which the material behavior is mod-eled with von Mises plasticity (i.e. imposing on the CPB06 criterion k = 1 and L ij = δ ij , with i,j = 1,...,6 and δ ij being the Kronecker unit delta) and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1)

Concluding remarks
In this work, we have used a linear stability analysis, a nonlinear two-zone model and unit-cell finite element calculations to investigate the necking formability of anisotropic and tension-compression asymmetric metallic sheets subjected to in-plane loading for different paths ranging from plane strain tension to near equibiaxial tension, and loading rates varying from quasi-static loading to strain rates up to 40000 s − 1 . We have studied three different materials -AZ31-Mg alloy, high purity α-titanium and OFHC copper-whose mechanical behavior is described with an elastic-plastic constitutive model with yielding defined by the orthotropic and asymmetric CPB06 criterion [15] which includes specific features to consider the evolution of plastic anisotropy and strength differential effect with accumulated plastic deformation [48]. AZ31-Mg alloy and high purity α-titanium display both plastic anisotropy and tension-compression asymmetry while OFHC copper displays only tension-compression asymmetry. The main novelty of this contribution is the extension of the recent work of N'souglo et al. [42] -which investigated materials described with the orthotropic yield criterion of Hill [24]-to consider anisotropic and tension-compression asymmetric materials with distortional hardening. The predictions obtained with the analytical models have been compared with the unit-cell finite element simulations and experimental data obtained from the literature. The linear stability has been calibrated using finite element calculations, assuming that a perturbation mode turns into a necking mode when the cumulative index reaches a critical value that depends on the loading path and the strain rate. As in the work of N'souglo et al. [42], the formability depends on the Lankford coefficient for AZ31-Mg alloy and high purity α-titanium which exhibit anisotropy, the general trend being that greater Lankford coefficient leads to lower formability. For OFHC copper, our results suggest that the formability increases for a material that displays a larger yield stress in uniaxial tension than in uniaxial compression. Furthermore, all the three approaches revealed that the effect of inertia on neck retardation depends on both the loading path and the yield criterion and that the stabilizing effect of inertia at high strain rates is enhanced for materials having low formability. The results obtained with the analytical models are in good accord with the experiments and consistent with the finite element results, for all loading paths and strain rates investigated. The nonlinear two-zone model, unlike the linear stability analysis, does not need prior calibration and yields accurate predictions that rarely deviate more than 10% from the results obtained with the unitcell calculations. All in all, the theoretical tools developed in this paper capture the mechanisms which control necking formability of metallic materials with complex mechanical response, providing guidelines and predictions for the application of low formability metals in high energy forming operations. Figure 10 provides a comparison of plane stress yield loci calculated with von Mises and CPB06 criteria for AZ31-Mg alloy, high purity α-titanium and OFHC copper. The yellow markers indicate the tensile plane strain and equibiaxial tension stress states (the loading paths investigated in this work lie within the yellow markers). For the three materials, the yield locus calculated with CPB06 is exterior to the von Mises ellipsoid, such that the curvature of the yield loci near equibiaxial tension is greater when the anisotropy and tension-compression asymmetry of the material is taken into account. Sowerby and Duncan [56], Parmar and Mellor [46] and N'souglo et al. [42], among others, have reported that for materials modeled with Hill [24] criterion, the formability decreases as the yield locus stretches out (i.e., as the maximum curvature of the yield locus increases compared to a von Mises material), and this trend seems to be consistent with the results obtained in this work with CPB06 criterion, see Section "Analysis and results". The greatest differences between CPB06 and von Mises yield loci correspond to AZ31-Mg alloy and high purity α-titanium, while for OFHC copper the differences are less, as in the forming limit diagrams reported in Section "Analysis and results".

Appendix B: Calibration of the linear stability analysis
The calibration of the linear stability analysis has been performed following the procedure proposed in Section 3.1 of Kumar et al. [33]. As stated in Section "Problem statement", we consider that the stability analysis predicts the formation of a neck when the cumulative instability index reaches a critical value I c which depends on the imposed initial major strain rate and on the loading path (note that in N'souglo et al. [42,43] the critical instability index leading to necking formation was considered to be constant): where a, b, c, d and e are calibration coefficients to be determined. For that purpose, we use a sequential procedure consisting of two-steps in which the results of the stability analysis for the evolution of the major necking strain neck xx with the necking wavelength L are fitted to unit-cell finite element simulations in which the material behaviour is modelled with von Mises plasticity (i.e. imposing on the CPB06 criterion k = 1 and L ij = δ ij , with i,j = 1,..., 6 and δ ij being the Kronecker unit delta) and with values of initial density, elastic constants and parameters of the yield stress corresponding to AZ31-Mg alloy (see Table 1).
1. For plane strain tension χ = 0 and initial major strain rates ̇0 xx = 100 s −1 , 5000 s − 1 and 20000 s − 1 , the best fitting between stability analysis and finite element results is obtained for I c = 2.5, 2.7 and 3.2, see Fig. 11a, leading to the values of a, c and e reported in Table 3. 2. For imposed initial major strain rate ̇0 xx = 5000 s −1 and loading paths χ = 0.5 and 0.75, the best fitting between stability analysis and finite element results is obtained for I c = 4 and 4.4, see Fig. 11b, leading to the values of b and d reported in Table 3.
On the other hand, the value of the critical instability index coefficients depends on the mechanical behavior of the material used for the calibration. Namely, using the initial density, elastic constants and parameters of the yield stress corresponding to high purity α-titanium and OFHC copper, see Table 1, different values of a,...,e are obtained, see Table 3. While all the stability analysis results presented in Section "Analysis and results" are computed with the calibration coefficients obtained with the material parameters corresponding to AZ31-Mg alloy; we have checked that the same qualitative results -same trends and conclusions-are obtained calibrating the stability analysis with the material parameters corresponding to high purity α-titanium and OFHC copper, see Appendix C.
Appendix C: The influence of the stability analysis calibration procedure Figure 12 compares forming limit diagrams, c xx versus c yy , predicted by the linear stability analysis using the three sets of critical instability index coefficients reported in Table 3. Results correspond AZ31-Mg alloy with ψ = 0 ∘ and 90 ∘ , see Fig. 12a, high purity α-titanium with ψ = 0 ∘ and 90 ∘ , see Fig. 12b, and OFHC copper, see Fig. 12c. A comparison is performed with calculations in which the material behavior is modeled with von Mises plasticity and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1). The initial major strain rate is ̇0 xx = 20000 s −1 .
The predictions of the stability analysis are quantitatively different for the three sets of critical instability index coefficients. Namely, the calibration with the material parameters corresponding to OFHC copper yields limit strains slightly larger than the results obtained with the parameters of AZ31-Mg alloy and high purity α-titanium (which are very similar), the differences increasing near equibiaxial tension. On the other hand, note that the relative order of the c xx − c yy curves obtained with CPB06 and von Mises yield criteria is the same for the three sets of coefficients. For instance, in Fig. 12a, for the three sets of critical instability index coefficients, the greater critical major necking strain corresponds to the von Mises material, and the smaller to the CPB06 material with ψ = 90 ∘ . Similar analysis can be performed for the results reported in Fig. 12b and c, which shows that the trends and conclusions obtained in this paper with the stability analysis calculations are independent of the set of critical instability index coefficients employed. Appendix D: The effect of inertia at low strain rates Figure 13 shows forming limit diagrams, c xx versus c yy , calculated with unit-cell finite element simulations for two different values of the initial major strain rate, ̇0 xx = 0.0001 s −1 and 100 s − 1 . The results correspond to AZ31-Mg alloy with ψ = 0 ∘ , see Fig. 13a, high purity α-titanium with ψ = 0 ∘ , see Fig. 13b, and OFHC copper, see Fig. 13c. A comparison is performed with calculations in which the material behavior is modeled with von Mises plasticity and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1). Note that the differences between the results obtained for ̇0 xx = 0.0001 s −1 (black lines) and ̇0 xx = 100 s −1 (red markers) are negligible, which makes apparent that the effect of inertia on necking formability for the materials investigated is very small for ̇0 xx = 100 s −1 (i.e. this strain rate can be considered quasi-static loading for the calculations reported in this paper).  Fig. 12 Comparison between forming limit diagrams obtained calibrating the stability analysis with values of initial density, elastic constants and parameters of the yield stress corresponding to AZ31-Mg alloy, high purity α-titanium and OFHC copper. Results for CPB06 and von Mises criteria: (a) AZ31-Mg alloy, (b) high purity α-titanium and (c) OFHC copper. The initial major strain rate is ̇0 xx = 20000 s −1 Fig. 13 Forming limit diagrams obtained with unit-cell finite element calculations for two different initial major strain rates ̇0 xx = 0.0001 s −1 and ̇0 xx = 100 s −1 . Results corresponding to: (a) AZ31-Mg alloy for ψ = 0 ∘ , (b) high purity α-titanium for ψ = 0 ∘ and (c) OFHC copper. Comparison with calculations in which the material behavior is modeled with von Mises plasticity (i.e. imposing on the CPB06 criterion k = 1 and L ij = δ ij , with i,j = 1,...,6 and δ ij being the Kronecker unit delta) and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Appendix E: Influence of tension-compression asymmetry on the formability Figure 14 shows AZ31-Mg alloy plane stress theoretical yield loci calculated with CPB06 criterion for biaxial loading conditions, σ I versus σ II , for two different values of the effective plastic strain, ̄p = 0.05 and ̄p = 0.1 . The case where both evolving anisotropy and tension-compression asymmetry is considered and the case where only evolving anisotropy is considered (k = 0) are compared. Note that the value of the tension-compression asymmetry parameter affects the shape of the yield locus predicted by the CPB06 criterion in all the four quadrants. The change is shape in the first quadrant is intrinsic to the mathematical construction of the model. The influence of the tension-compression asymmetry parameter in the material formability is illustrated in Fig. 15. The figure shows the forming limit diagrams for AZ31-Mg alloy for the case where both evolving anisotropy and tension-compression asymmetry is considered and for the case where only evolving anisotropy is considered. It is evident that although all the loading paths investigated lie within the first quadrant, the tension-compression asymmetry affects the formability limit.

Conflict of Interests
The authors declare that they have no conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.