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 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 [10] which includes specific features to account for the evolution of plastic orthotropy and strength differential effect with accumulated plastic deformation [37]. From a methodological standpoint, the main novelty of this paper with respect to the recent work of N’souglo et al. [32] –which investigated materials with yielding described by the orthotropic criterion of Hill [19]– 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. Keywords: Dynamic formability, Tension-compression asymmetry, Plastic orthotropy, Linear stability analysis, Two-zone model, Finite elements ∗Corresponding author. E-mail address: knsouglo@ing.uc3m.es Preprint submitted to December 2, 2021 2 1. Introduction Formability of sheet metal products generally increases at high strain rates [47, 4, 5, 42, 16]. The stabilizing effect of inertia, which delays necking formation and enhances workpiece ductility [15, 43, 24], 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 [31]. For instance, Balanethiram and Daehn [5] 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 [17] 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 quasi-static loading by forming the aluminum specimens with a hemispherical punch into a round die cavity. Consistent with the earlier work of Balanethiram and Daehn [5], the results of Golovashchenko [17] 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. [13] 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 [5], the strain rate in the workpiece during the tests was estimated to be around 1000 s−1. Dariani et al. [13] 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. [41] 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 [18]. Xu et al. [50, 48] carried out electromagnetic free-bulging forming experiments with AZ31 magnesium sheets 3 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 [49]. The thickness of the samples was 1 mm, as in the experiments of Xu et al. [50, 48]. The tests were performed at different forming speeds, with the maximum strain rate attained in the experiments varying between 1200 s−1and 3500 s−1. As in the free-bulging tests of Xu et al. [50, 48], 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. [49] 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. [27] using a driver plate made of AA5052-O with 2 mm thickness. As in the formability tests performed by Xu et al. [50, 48] 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 quasi-static 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 [26]. The strain rate attained in the experiments was estimated to be around 4000 s−1. The tests of Li et al. [26] 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, Balanethiram and Daehn [4] used electrohydraulic forming to obtain the limit strains of Interstitial Free iron sheets of 0 .84 mm 4 thickness deformed at strain rate of 1000 s−1under nearly plane strain conditions. The comparison with the corresponding quasi-static forming limits showed that the ductility of the material increased by nearly three times at high strain rates. Shortly after, Balanethiram and Daehn [5] 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. [42] 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. [18] 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 [5] that the ben- eficial 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 [26, 50, 48, 49], and also on the mechanical behavior of the material, since the extent of the formability improvement is different for different metal sheets [18, 42]. Nevertheless, as of today, it is not clear which features of the constitutive behavior of the material affect 5 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 [34, 46, 41, 51, 24], which is a crude assumption since most sheet metal products display plastic anisotropy. An exception is the recent work of N’souglo et al. [32], who developed a linear stability analysis and a nonlinear two-zone model, and performed finite element simulations for materials with yielding described by Hill [19] 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 equibiaxial stretching, and to strain rates varying between 100 s−1and 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, two-zone 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−1inertia 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. [32] did not provide a quantitative comparison with experimental data, and also because the Hill yield function [19] fails when it comes to describe the mechanical response of metal sheets which exhibit complex plastic anisotropy [21, 6]. In particular, the theory of Hill [19] 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 [36, 7]. In addition, the formulation of Hill [19] 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 [9, 10, 38]. Therefore, it is evident that the work of N’souglo et al. [32] 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). Hence, in this paper we have developed a linear stability analysis and a nonlinear two-zone model to construct 6 forming limit diagrams for elastic-plastic materials with yielding described by the CPB06 criterion [10]. 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 [5, 27, 50, 48], and reasonable agreement has been obtained for a wide range of loading rates ranging from quasi-static loading to strain rates up to 40000 s−1, and for loading paths ranging from plane strain to 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. 2. 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 [10] which includes specific features to consider the evolution of plastic anisotropy and strength differential effect with accumulated plastic deformation [37]. The formulation of the constitutive model is outlined in Section 2.1, and the values of the constitutive parameters of the materials investigated are provided in Section 2.2. 2.1. Formulation The rate of deformation tensor dis the sum of an elastic part deand a plastic part dp: d=de+dp(1) where the relationship between the elastic part and the rate of the stress is: ▽σ=C:de(2) where▽σis an objective derivative of the Cauchy stress tensor, and Cis the fourth-order isotropic elastic tensor defined as: 7 C= 2GI′+K1⊗1 (3) with1being the unit second-order tensor and I′the unit deviatoric fourth-order tensor. Moreover GandKare the shear and bulk moduli, respectively. Assuming an associated plastic flow rule, the plastic part of the rate of deformation tensor is defined as: dp=˙λ∂¯σ ∂σ(4) where ˙λis the rate of the plastic multiplier and ¯ σis the effective stress [10] given by: ¯σ= ˜mp (|Σ1| −kΣ1)2+ (|Σ2| −kΣ2)2+ (|Σ3| −kΣ3)2 (5) where kis the tension-compression asymmetry parameter, see equation (10a), and Σ i(i= 1, ...,3) are the principal values of the transformed stress tensor Σdefined as: Σ=L:s (6) withLbeing an orthotropic, deviatoric, and symmetric fourth-order tensor, and sthe deviatoric part of the Cauchy stress tensor σ. In the Cartesian coordinate system ( X1, X2, X3) associated to the orthotropy axes of the material, the deviatoric part of the Cauchy stress tensor is represented by a 6-dimensional vector s= (s11, s22, s33, s12, s23, s13) and the fourth-order orthotropic tensor Lis represented by a 6 ×6 matrix: 8 L= L11L12L13 0 0 0 L12L22L23 0 0 0 L13L23L33 0 0 0 0 0 0 L44 0 0 0 0 0 0 L55 0 0 0 0 0 0 L66 (7) with the components of Lbeing the anisotropy coefficients, see equation (10b). Moreover, the parameter ˜ min equation (5) is defined as: ˜m=1p (|Φ1| −kΦ1)2+ (|Φ2| −kΦ2)2+ (|Φ3| −kΦ3)2(8) with Φi(i= 1, ...,3) expressed as follows: Φ1=2 3L11−1 3L12−1 3L13 Φ2=2 3L12−1 3L22−1 3L23 Φ3=2 3L13−1 3L23−1 3L33(9) The evolution of plastic anisotropy and tension-compression asymmetry with accumulated plastic strain is modeled using the methodology developed by Plunkett et al. [37]. Firstly, the values of the parameters kandLij (i, j= 1, ...,6), see equations (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): k(¯εp) =α(¯εp)k(¯εp,n) + (1 −α(¯εp))k ¯εp,n+1(10a) Lij(¯εp) =α(¯εp)Lij(¯εp,n) + (1 −α(¯εp))Lij ¯εp,n+1(10b) 9 where the interpolation parameter αis: α(¯εp) =¯εp,n+1−¯εp ¯εp,n+1−¯εp,n(11) such that α(¯εp,n) = 1 and α ¯εp,n+1 = 0. The yield function is defined as: f= ¯σ−σY (12) where σYis the tensile yield stress in the rolling direction of the material X1, given in the present work by the following relationship: σY=A0−A1exp−A2¯εp(13) with A0,A1andA2being material parameters. Notice that viscous and thermal effects on the yield stress evolution are neglected. Moreover, the effective plastic strain is ¯ εp=Rt 0˙¯εpdτ, where trefers to time, and ˙¯εpis the effective plastic strain rate given by: 10 ˙¯εp=  1 ˜m(1−k)s b2 1+b2 2+3k2−10k+ 3 3k2+ 2k+ 3 b2 3ifb3 |b|≤−(3k2+ 2k+ 3)p 6(k2+ 3)(3 k2+ 1) 1 ˜m(1 +k)s b2 1+b2 2+3k2+ 10k+ 3 3k2−2k+ 3 b2 3ifb3 |b|≥3k2−2k+ 3p 6(k2+ 3)(3 k2+ 1) 1 ˜m(1−k)s3k2−10k+ 3 3k2+ 2k+ 3 b2 1+b2 2+b2 3ifb1 |b|≤−(3k2+ 2k+ 3)p 6(k2+ 3)(3 k2+ 1) 1 ˜m(1 +k)s3k2+ 10k+ 3 3k2−2k+ 3 b2 1+b2 2+b2 3ifb1 |b|≥3k2−2k+ 3p 6(k2+ 3)(3 k2+ 1) 1 ˜m(1−k)s b2 1+3k2−10k+ 3 3k2+ 2k+ 3 b2 2+b2 3ifb2 |b|≤−(3k2+ 2k+ 3)p 6(k2+ 3)(3 k2+ 1) 1 ˜m(1 +k)s b2 1+3k2+ 10k+ 3 3k2−2k+ 3 b2 2+b2 3ifb2 |b|≥3k2−2k+ 3p 6(k2+ 3)(3 k2+ 1)(14) where |b|=p b2 1+b2 2+b2 3, and b1,b2andb3are the principal values of the transformed rate of deformation tensor defined as: b=H:dp(15) where H:G=I′withG=L:I′. Moreover, the work conjugacy relation: σ:dp= ¯σ˙¯εp(16) yields: ˙λ=˙¯εp(17) The formulation of the constitutive behaviour is completed with the Kuhn–Tucker loading–unloading condi- 11 tions: ˙¯εp⩾0, f ⩽0, ˙¯εpf= 0 (18) and the consistency condition during plastic loading: ˙f= 0 (19) 2.2. 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 1). In addition, the three materials display very different yielding behaviour, see Fig. 1, increasing the scope of the analysis carried out in Section 4. Table 1 shows the values of the initial density, the elastic constants and the parameters of the yield stress in the rolling direction [38]. We are aware that using a Voce-type hardening relation (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 [22, 24, 32]. 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 Lij(i, j= 1, ...,6) and the tension-compression asymmetry parameter kfor 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. [38] and Cazacu et al. [11]. Recall from Section 2.1 that the values of Lijand kfor effective plastic strains within the intervals reported in Table 2 are calculated using a linear interpolation function, see equation (10). Moreover, in absence of additional data, we consider that Lijandkdo not evolve for effective plastic strains outside the range of values included in Table 2, so that for larger strains the values of Lij andkfor AZ31-Mg alloy, high purity α-titanium and OFHC copper are the same employed for ¯ εp= 0.1, 0.3 and 12 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. Symbol Property and units AZ31-Mg alloy High purity α-titanium OFHC copper ρ0 Initial density (kg /m3) 1800 4500 8900 G Elastic shear modulus (GPa) 17.31 42.31 46.64 K Bulk modulus (GPa) 37.5 91.67 130.2 A0 Material constant (MPa) 315.4 341 253.8 A1 Material constant (MPa) 140.6 170 224.3 A2 Material constant 16.3 7.03 10.92 Table 1: Numerical values of the initial density, the elastic constants and the parameters of the yield stress in the rolling direction for AZ31-Mg alloy, high purity α-titanium and OFHC copper. Data taken from Revil-Baudard et al. [38]. Anisotropy and tension-compression asymmetry parameters Material ¯εpL11 L22 L33 L12 L13 L23 L44 L55 L66 k 0 1 1.090 3.342 -0.168 0.098 0.243 0.730 7.30 7.74 -0.625 0.05 1 1.090 3.342 -0.168 0.098 0.243 0.730 7.30 7.74 -0.625 AZ31-Mg alloy 0.06 1 1.072 2.905 -0.595 -0.279 -0.096 1.039 10.2 11.02 -0.520 0.08 1 1.099 1.432 -0.817 -0.516 -0.350 1.128 11.21 11.94 -0.215 0.10 1 1.082 0.885 -0.762 -0.657 -0.509 1.058 10.12 11.21 -0.169 0 1 0.850 0.836 0.281 0.381 0.439 0.476 0.527 0.639 -0.024 0.025 1 0.981 1.013 0.847 0.885 0.904 0.128 0.120 0.154 -0.070 High purity α-titanium 0.05 1 0.971 1.013 0.750 0.801 0.824 0.218 0.206 0.259 -0.061 0.075 1 0.917 0.985 0.374 0.478 0.516 0.555 0.528 0.657 -0.084 0.2 1 0.975 0.898 -0.018 0.034 0.052 1.050 0.986 1.285 -0.392 0.3 1 1.032 0.954 -0.024 0.001 0.010 1.161 1.063 1.436 -0.650 0 1 1 1 0 0 0 1 1 1 0 0.03 1 1 1 0 0 0 1 1 1 -0.084 OFHC copper 0.06 1 1 1 0 0 0 1 1 1 -0.127 0.09 1 1 1 0 0 0 1 1 1 -0.137 0.12 1 1 1 0 0 0 1 1 1 -0.137 Table 2: Numerical values of the anisotropy coefficients Lijand the tension-compression asymmetry parameter kfor different values of the effective plastic strain ¯ εpfor AZ31-Mg alloy, high purity α-titanium and OFHC copper. Data taken from Revil-Baudard et al. [38] and Cazacu et al. [11]. Fig. 1 pictures plane stress theoretical yield loci for biaxial loading conditions, σIversus σII, where σIandσ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 13 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 transverse directions is different. Moreover, for OFHC copper, 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 [30] 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. 3. Problem statement The problem addressed is that of a thin plate of initial thickness h0and edges of initial length L0, modeled with the constitutive framework described in Section 2, and subjected to constant and opposed stretching velocities on opposed sides, see Fig. 2. The applied velocities are Vx=±˙ε0 xxL0/2 and Vy=±˙ε0 yyL0/2, with ˙ ε0 xxand ˙ε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 (equibiaxial stretching) in the calculations presented in this paper. Hereinafter, ˙ ε0 xxand ˙ε0 yywill 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 ( X1, X2, X3) is a Lagrangian Cartesian frame defined by an angle ψbetween XandX1, with X1,X2andX3being 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). 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. 14 -500-2500250500 -500 -250 0 250 500ε p = 0.05 ε p = 0.06 ε p = 0.08 ε p = 0.1 (288,0) (-311,0)(297,0) (-341,0)σYMaterial: AZ31-Mg alloyσII (MPa) σI (MPa) (a) -600-400-2000200400600 -600 -400 -200 0 200 400 600σII (MPa) σI (MPa)ε p = 0 ε p = 0.025 ε p = 0.05 ε p = 0.075 ε p = 0.2 ε p = 0.3(323,0) (-433,0)(313,0) (-450,0)σYMaterial: High purity α-titanium (b) -300-200-1000100200300 -300 -200 -100 0 100 200 300ε p = 0 ε p = 0.03 ε p = 0.06 ε p = 0.09 ε p = 0.12(195,0) (-213,0)(195,0) (-213,0)σYMaterial: OFHC copperσII (MPa) σI (MPa) (c) Figure 1: Plane stress theoretical yield loci calculated with CPB06 criterion for biaxial loading conditions, σIversus σII, for different values of the effective plastic strain: (a) AZ31-Mg alloy, (b) high purity α-titanium and (c) OFHC copper. 15 Figure 2: Schematic representation of the geometry and boundary conditions of the problem addressed: a thin plate of initial thickness h0and edges of initial length L0subjected to constant and opposed stretching velocities Vx=±˙ε0 xxL0/2 and Vy=±˙ε0 yyL0/2 on opposed sides. The Lagrangian Cartesian coordinate system associated to the applied velocity field is denoted by ( X, Y, Z ), while ( X1, X2, X3) is a Lagrangian Cartesian frame defined by an angle ψbetween XandX1, with X1,X2andX3being 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. 16 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. [32] –which investigated materials described with the orthotropic yield criterion of Hill [19]– to consider the orthotropic and tension-compression asymmetric yield criterion with distortional hardening described in Section 2. A brief note on the three approaches is provided next, while the reader is referred to the papers of N’souglo et al. [32] and Jacques [24] 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. [40] and N’souglo et al. [32] for a complete description of the finite element model. In addition, in Section 5 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 4. 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 [51, 33, 39] so that the stress multiaxiality effects that develop inside the necked section have been approximated using Bridgman [8] 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. [32] 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 the perturbation wavenumber ξ. The root of the polynomial with physical meaning is the one with the greatest positive real part [14], hereinafter denoted by η+. The wavenumber is related to the normalized perturbation wavelength ¯L=L0/h0as¯L=2πξ−1 h0. Hereinafter, ¯Lwill be also referred to as necking wavelength. On the other hand, the instantaneous growth rate is used to compute the cumulative instability index ¯I=Rt 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 17 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 [24] is an extension of the classical analysis of localized necking proposed by Marciniak and Kuczy´ nski [28] 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 L0(see Fig. 4 in N’souglo et al. [32]). Note that L0also corresponds to the neck spacing. The initial thickness in the imperfection zone is hA,0=h0(1−∆), ∆ being the normalized imperfection amplitude. The ratio between the initial length of the imperfection LA,0 xand the initial unit-cell length is denoted by R=LA,0 x/L0. 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 [24] and was found to remain appropriate for orthotropic materials obeying the Hill [19] yield criterion [32]. The methodology to obtain necking strains with the two-zone model is described in Section 2.5 of Jacques [24]. Computations have been carried out for several normalized cell sizes, 0 .5≤¯L≤6 with ¯L=L0 h0andh0= 2 mm, in order to identify the critical neck spacing and the corresponding necking strains (critical necking strains). The results presented in Section 4 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 4.1) 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=h0−δ 1 + cos2πX L0 , 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 0 ≤X≤L0 2, 0≤Y≤L0 2and 0 ≤Z≤h0 2. 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=L0 h0andh0= 2 mm fixed, and the normalized imperfection amplitude is ∆ =2δ h0= 0.2%. Note that the normalized cell size ¯Lwill 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. [32], which are consistent with the stability analysis and the two-zone model. The unit-cell 18 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 Xaxis. 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µm3(elements are initially shorter along the loading directions). The constitutive model presented in Section 2 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 [20]. The UMAT [2] has been used to obtain the low strain rate finite element results, 0 .0001 s−1and 100 s−1, reported in Figs. 3 and D.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. 4. 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 equibiaxial tension (0 ≤χ≤1, see Section 3) 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. see Dariani et al. [13], Golovashchenko et al. [18] and Li et al. [26], among others), but allows a further exploration of inertia effects. 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 3). 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 3), and the stability analysis is calibrated using finite element simulations different from those presented below (see Appendix B). In Section 4.1, quasi-static forming limit diagrams obtained for the three materials are systematically compared with experimental results available in the literature and with calculations in 19 which the material behavior is modeled with von Mises plasticity (i.e. imposing on the CPB06 criterion k= 1 and Lij=δij, with i, j= 1, ...,6 and δijbeing the Kronecker unit delta) and the same values of initial density, elastic constants and parameters of the yield stress (see Table 1). In Section 4.2, calculations performed for different loading rates are compared to identify the effect of inertia on the necking formability of the materials investigated. 4.1. The influence of material behavior Fig. 3 displays forming limit diagrams (critical major necking strain εc xxversus 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. [32], 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 xxis 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 Lij=δij, with i, j= 1, ...,6 andδijbeing 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εc xx−εc yycurves 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 xxbeing 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. [32] for materials modeled with Hill [19] criterion, the relative order of the εc xx−εc yycurves 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 20 thatr0◦= 1.40 and r90◦= 1.88 for ¯ εp≥0.1, see equation (17) in Cazacu et al. [10]) 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. [12] (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, r0◦= 1.09 and r90◦= 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 [45] 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 [29] 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 kpositive (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. 21 00.10.20.30.40.50.60.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7B D F H J L N P R T VCritical major necking strain, εxx Critical minor necking strain, ε yyccQuasi-static loading FEM: von Mises FEM: CPB06 - ψ=0 FEM: CPB06 - ψ=90 LSA: von Mises LSA: CPB06 - ψ=0 LSA: CPB06 - ψ=90 2ZM: von Mises 2ZM: CPB06 - ψ=0 2ZM: CPB06 - ψ=90 Experimentso o o o o oMaterial: AZ31-Mg alloy ε xx=ε yyc c (a) 00.10.20.30.40.50.60.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7FEM Mises FEM RD FEM TD LSA Mises LSA RD LSA TD 2ZM Mises 2ZM RD 2ZM TD Experiments VCritical major necking strain, εxx Critical minor necking strain, ε yyccQuasi-static loadingMaterial: High purity α-titanium FEM: von Mises FEM: CPB06 - ψ=0 FEM: CPB06 - ψ=90 LSA: von Mises LSA: CPB06 - ψ=0 LSA: CPB06 - ψ=90 2ZM: von Mises 2ZM: CPB06 - ψ=0 2ZM: CPB06 - ψ=90 Experimentso o o o o o ε xx=ε yyc c (b) 00.10.20.30.40.50.60.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7FEM Mises FEM Copper LSA Mises LSA Copper 2ZM Mises 2ZM Copper Experiments PCritical major necking strain, εxx Critical minor necking strain, ε yycc FEM: von Mises FEM: CPB06 LSA: von Mises LSA: CPB06 2ZM: von Mises 2ZM: CPB06 ExperimentsMaterial: OFHC copper Quasi-static loading ε xx=ε yyc c (c) Figure 3: Quasi-static forming limit diagrams obtained with unit-cell calculations (FEM), linear stability analysis (LSA) and two- zone model (2ZM) for: (a) AZ31-Mg alloy with material orientations ψ= 0◦andψ= 90◦, (b) high purity α-titanium with material orientations ψ= 0◦andψ= 90◦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 Lij=δij, with i, j= 1, ...,6 and δijbeing 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. [12], Stachowicz [45] and Melander [29], respectively. (For interpretation of the