In this paper, we have investigated the e ect of material orthotropy on the formability of metallic sheets sub- jected to dynamic biaxial stretching. For that purpose, we have devised an original three-pronged methodology which includes a linear stability analysis, a nonlinear two-zone model and nite element calculations. We have studied 5 di erent materials whose mechanical behavior is described with an elastic isotropic, plastic anisotropic constitutive model with yielding based on Hill (1948) criterion. The linear stability analysis and the nonlinear two-zone model are extensions of the formulations developed by Zaera et al. (2015) and Jacques (2020), respec- tively, to consider Hill (1948) plasticity. The nite element calculations are performed with ABAQUS/Explicit (2016) using the unit-cell model developed by Rodr guez-Mart nez et al. (2017), which includes a sinusoidal spatial imperfection to favor necking localization. The predictions of the stability analysis and the two-zone model are systematically compared against the nite element results {which are considered as the reference approach to validate the theoretical models{ for loading paths ranging from plane strain stretching to equibiaxial stretching, and for di erent strain rates ranging from 100 s1to 50000 s1. The stability analysis and the two-zone model yield the same overall trends obtained with the nite element simulations for the 5 materials investigated, and for most of the strain rates and loading paths the agreement for the necking strains is also quantitative. Notably, the di erences between the nite element results and the two-zone model rarely go beyond 5%. Altogether, the results presented in this work provide new insights into the mechanisms which control dynamic formability of anisotropic metallic sheets. Keywords: Dynamic formability, Necking, Plastic orthotropy, Linear stability analysis, Two-zone model, Finite elements Corresponding author. Tel. +33 2 98 34 89 36; fax. +33 2 98 34 87 30. E-mail address: nicolas.jacques@ensta-bretagne.fr Preprint submitted to Journal of the Mechanics and Physics of Solids September 14, 2020 2 1. Introduction Interest of aeronautical and automotive industries in high energy rate forming (HERF) processes, such as explosive, electromagnetic and electrohydraulic forming, has been growing in recent years due to their potential to form aluminium, titanium, high-strength steel and other low formability metallic materials (Song et al., 2004). The high strain rates achieved in HERF processes in comparison to conventional quasi-static metal forming operations lead to signi cant changes in the constitutive behavior of the materials and give rise to inertial e ects which contribute to delay necking localization and ductile fracture (Mamalis et al., 2004; Gayakwad et al., 2014). The rst experimental evidence that inertia e ects improve the formability of metal sheets was most likely provided by Wood (1967), who showed that the ductility of several materials (stainless steel, vascojet tool steel, Ren e 41 superalloy) is signi cantly increased during high-explosive forming operations, such as tube and dome bulging. More recently, Balanethiram and Daehn (1992) obtained the dynamic formability diagrams of Interstitial Free iron sheets with 0 :84 mm thickness tested at strain rates of 103s1using electrohydraulic forming (the sheets were formed inside a conical die with 90apex angle). The comparison with the corresponding quasi-static forming limits showed that the ductility of the specimens near plane strain increased by three times at high strain rates. Shortly after, Balanethiram and Daehn (1994) presented additional electrohydraulic experiments showing that the formability of AA 6061 and OFHC copper sheets also increases at high strain rates, suggesting that the bene cial contribution of inertial forces to the ductility of metals sheets is both large and general (Balanethiram and Daehn, 1994). The experiments on OFHC copper were performed using sheets of 0 :34 and 0:79 mm thickness, with the thicker samples showing much greater improvement in formability at high strain rates due to the increase of inertia e ects with the specimen thickness. Hence, Daehn et al. (1995) concluded that when inertial forces are at work in stabilizing necking, size e ects become important. Years later, Seth et al. (2005) investigated experimentally the dynamic formability of sheet specimens made of ve di erent cold-rolled steels with thicknesses varying from 0 :15 to 0:38 mm. The tests consisted of launching electromagnetically the steel sheets at a shaped punch at velocities of 50{220 m =s. Similarly to the experimental evidence presented by Balanethiram and Daehn (1992, 1994) for Interstitial Free iron, AA 6061 and OFHC copper, the ductility of these ve steels was also dramatically increased at high strain rates as compared to the quasi-static forming limit diagrams. The experimental results reported by Golovashchenko (2007) for the electromagnetic forming of AA 6111-T4 and AA 5754 sheets into conical and v-shape dies also showed that the ductility of the electromagnetic formed samples can be 2:5 times greater than under quasi-static loading. Dariani et al. (2009) 3 obtained the forming limit diagrams of AA 6061-T6 and 1045 steel under quasi-static and explosive loading (in the latter case the strain rate reaches 103s1). These experiments revealed that the ductility of the explosively formed specimens increases by a factor of 2 :5 for AA6061-T6 and by factor of 1 :6 for 1045 steel. Golovashchenko et al. (2013) investigated the dynamic formability of four di erent dual phase steels (DP500, DP590, DP780 and DP980) at strain rates up to 2104s1using electrohydraulic forming. As compared to the quasi-static case, the relative improvement in plane strain formability in electrohydraulic forming conditions was between 63% and 190%, depending on the grade of dual phase steel. Moreover, Rohatgi et al. (2014) performed electrohydraulic forming experiments which showed that the dynamic formability of AA 5182-O sheets of 1 mm thickness, relative to its corresponding quasi-static formability, increases by 2:5 times and6:5 times under free-forming and when forming inside a conical die, respectively. Theoretical investigations to rationalize the e ect of inertia on neck retardation have been conducted using both linear and nonlinear analyses. Most of the linear analyses are based on the linear perturbation theory, which evaluates the growth of an in nitesimal perturbation superimposed on the homogenous background solution of the problem at hand. Fressen- geas and Molinari (1985) showed that if the perturbation grows faster than the background solution, the plastic ow is unstable and a neck-like deformation eld can develop. Fressengeas and Molinari (1994) extended the plane strain linear stability analysis of Hutchinson et al. (1978b) to include inertia e ects, and they showed that inertia decreases the growth rate of perturbations, slowing down necking development and enhancing material ductility. In addition, it was demonstrated that at high strain rates, unlike what happens under quasi-static loadings, there is a perturbation of nite wavelength (so-called the critical wavelength) for which the growth rate of the perturbation is maximum. This critical wavelength was used by Fressengeas and Molinari (1994) to determine the average spacing between necks in the multiple localization patterns that emerge in shaped-charge jets, and satisfactory agreement was found between theory and experimental data available in the literature. Mercier et al. (2010) extended the plane strain linear stability analysis developed by Fressengeas and Molinari (1994) to thermoviscoplastic materials to predict the critical conditions for the nucleation of necks and the spacing between necks in hemispherical shells made of copper and tantalum, and subjected to dynamic expansion at strain rates of 104s1. The stability anal- ysis predictions were in reasonable agreement with the number of necks obtained in the experiments and showed that inertia delays the formation of necks (which indeed nucleate for strains greater than the Consid ere (1885) strain, see also Vaz-Romero et al. (2017)). Moreover, Zaera et al. (2015) extended the linear stability analysis 4 developed by Dudzinski and Molinari (1991) to consider inertia e ects and study the formation of multiple necks in thin metal sheets modeled with von Mises plasticity and subjected to biaxial loading at high strain rates. Using the concept of e ective instability introduced by Dudzinski and Molinari (1991), which states that a perturbation mode turns into a necking mode when the grow rate of the perturbation reaches a critical value, Zaera et al. (2015) used the results of the stability analysis to construct dynamic forming limit diagrams. The predictions of the stability analysis were in qualitative agreement with nite element results which showed that the average spacing between necks increases as the loading path moves away from uniaxial tension and approaches equibiaxial tension (see also Rodr guez-Mart nez et al. (2017)). The stability analysis also predicted in agreement with the nite ele- ment calculations the orientation of the necking bands that form the localization pattern, which are aligned with the direction of zero stretch for loading paths between uniaxial tension and plane strain (Hill, 1952), and they are perpendicular to the principal loading direction when the loading path lies between plane strain and equibiaxial tension. In addition, the nite element calculations indicated that the formability of metal sheets increases as the loading rate increases, in agreement with the experimental evidence discussed above (see also Rodr guez-Mart nez et al. (2017)). Moreover, the nonlinear analyses are generally based on the hypothesis of Marciniak and Kuczy nski (1967) that necking is initiated by an ab initio geometric imperfection of the specimen in the form of a thickness variation, and they are usually known as two-zone models or long-wavelength analyses (Hutchinson et al., 1978b). For instance, Fressengeas and Molinari (1985) developed a one-dimensional model based on the long-wavelength approximation including material inertia to investigate the e ect of strain rate on necking localization under uniaxial tension. Despite the limitations of the model, that did not consider the multiaxial stress state that develops in a necked section, Fressengeas and Molinari (1985) succeeded to provide one of the rst theoretical veri cations that neck formation is delayed at high strain rates. Years later, Xue et al. (2008) formulated a dynamic two-zone analysis to investigate necking in thin plates subjected to plane strain extension. The theoretical predictions were compared against nite element simulations in which the plate was modeled as an array of unit-cells with sinusoidal periodic geometric imperfections. The calculations were performed for strain rates ranging from 100 s1to 2000 s1(but for a plate with large thickness of 20 mm, which promotes inertia e ects). The two-zone analysis did a good job capturing neck retardation due to inertia e ects, however, because triaxiality e ects were not included in the formulation of the model (as in the one-dimensional model of Fressengeas and Molinari (1985)), short necking wavelengths were not suppressed, and the two-zone analysis failed to predict the emergence of a critical necking wavelength associated to the minimum necking strain (which determines the average spacing between necks in 5 multiple localization problems, as mentioned before). Very recently, Jacques (2020) developed a two-zone model that extends the work of Xue et al. (2008) to consider any arbitrary in-plane loading and includes stress triaxiality e ects using Bridgman (1952) approximation to account for the hydrostatic stresses which develop inside a necked region. The predictions of the model of Jacques (2020) were validated against the dynamic forming experiments of Golovashchenko et al. (2013) and the unit-cell nite element calculations of Rodr guez-Mart nez et al. (2017), who extended the computational model of Xue et al. (2008) to consider loading paths ranging from plane strain to equibiaxial tension. Excellent agreement was found between theoretical model predictions and nite element calculations, providing a de nite theoretical con rmation of the stabilizing e ect of inertia on neck development for di erent strain rates and loading paths, ranging from uniaxial tension to equibiaxial tension. The two-zone model of Jacques (2020) also succeeded to capture the decrease of the critical necking wavelength with the strain rate, and its increase as the loading path moves away from plane strain to equibiaxial tension. The increase in ductility at high strain rates predicted by the two-zone model was also found in agreement with the experimental data of Golovashchenko et al. (2013). However, most of the theoretical works cited in this introduction were speci cally formulated for isotropic materials. To the authors' knowledge, neither stability analysis nor two-zone models have been yet developed to study the in uence of inertia in the dynamic formability of anisotropic metal sheets. This is an important gap in the literature, as it is known that most sheet metal products are anisotropic. Thus, in this paper, we extend the stability analysis of Zaera et al. (2015) and the two-zone model of Jacques (2020) to consider orthotropic Hill (1948) plasticity. The necking strains calculated with both theoretical models are compared to nite element calculations performed with the unit-cell model used by Rodr guez-Mart nez et al. (2017), for a wide range of strain rates ranging from 100 s1to 50000 s1, and for loading paths varying from plane strain to equibiaxial tension. We consider that the loading axes are parallel to the orthotropy axes of the material. Therefore, for this range of loading paths, the orientation of the necking band is perpendicular to the principal loading direction, so that the orientation of the geometric imperfection to be included in the unit-cell nite element model is known. The stability analysis and the two-zone model predict the same overall trends obtained with the nite element calculations, and for most of the strain rates and loading paths investigated, the agreement for the necking strains is also quantitative, notably the di erences between the nite element results and the two-zone model rarely go beyond 5%. Altogether, this paper presents several methodologies to predict dynamic formability diagrams for anisotropic materials, and provides additional insights into the stabilizing role of inertia on neck development for 6 di erent loading paths. 2. Constitutive framework The mechanical behavior of the materials studied in this paper is described using an elastic isotropic, plastic anisotropic constitutive model with yielding based on Hill's criterion (Hill, 1948). The formulation of the con- stitutive model is outlined in Section 2.1, and the values of the parameters corresponding to the ve materials investigated are given in Section 2.2. 2.1. Formulation We consider the additive decomposition of the total rate of deformation tensor dinto an elastic part deand a plastic part dpas follows: d=de+dp(1) The elastic part of the rate of deformation tensor is related to the rate of the stress by the following linear hypo-elastic law: O=C:de(2) whereOis an objective derivative of the Cauchy stress tensor, which corresponds to the Green-Naghdi derivative in ABAQUS/Explicit (2016) (when a VUMAT user-subroutine is used), and Cis the tensor of isotropic elastic moduli given by: C= 2GI0+K1 1 (3) with1being the unit second-order tensor and I0being the unit deviatoric fourth-order tensor. Moreover Gand Kare the shear and bulk moduli, respectively. The yield condition is de ned as: 7 f= Y= 0 (4) where is the e ective stress associated to the Hill (1948) criterion: = (G+H)2 11+ (F+H)2 22+ (F+G)2 33 2H11222F22332G1133+ 2N2 12+ 2L2 23+ 2M2 131=2(5) withij(i;j= 1;:::;3) being the components of the Cauchy stress tensor in the Cartesian coordinate system (X1;X2;X3) associated to the orthotropy axes (see Figs. 3 and 4). Moreover F,G,H,L,MandNare the anisotropy parameters (see Section 2.2). Note that the von Mises (1928) yield criterion is recovered if F=G= H= 0:5 andL=M=N= 1:5. In equation (4), Yis the yield stress of the material in uniaxial tension in the rolling direction, which is considered to evolve with the e ective plastic strain  "p, the e ective plastic strain rate _"pand the temperature T, according to the following power-law type relation: Y=0 "0+ "pn_"p _"refmT Tref (6) where0,"0,n, _"ref,m,Trefandare material parameters (see Section 2.2). Moreover, most of the results presented in this paper have been obtained with a rate-independent behaviour ( m= 0) under isothermal conditions of deformation ( = 0, see equation (10)). This facilitates the identi cation of the speci c role played by material orthotropy in the dynamic formability of metallic materials (see Section 6). The in uence of strain-rate sensitivity and thermal softening is speci cally discussed in Section 6.4. The e ective plastic strain is related to the e ective plastic strain rate as follows: "p=Zt 0_"p( )d (7) 8 Assuming an associated plastic ow rule, the plastic part of the rate of deformation tensor is: dp=_"p@ @(8) where the e ective plastic strain rate _"pis de ned as: _"p=F(dp 11)2+G(dp 22)2+H(dp 33)2 FG+FH+GH+2(dp 12)2 N+2(dp 23)2 L+2(dp 13)2 M1=2 (9) withdp ij(i;j= 1;:::;3) being the components of the plastic part of the rate of deformation tensor in the Cartesian coordinate system ( X1;X2;X3) associated to the orthotropy axes (see Figs. 3 and 4). Moreover, assuming adiabatic conditions of deformation (no heat ux) and considering that plastic work is the only source of heating, the evolution of temperature is given by: _T= _"p Cp(10) where is the Taylor-Quinney coecient, is the current material density and Cpis the speci c heat. The formulation of the constitutive behaviour is completed with the Kuhn{Tucker loading{unloading condi- tions: _"p>0; f 60; _"pf= 0 (11) and the consistency condition during plastic loading: _f= 0 (12) 9 Hill's model (Hill, 1948) highlights due to its simplicity, versatility and easy calibration procedure; so that 70 years after its publication, the criterion of Hill (1948) still is widely used to describe the directionality of plastic properties of metallic materials in scienti c and technological applications (Padmanabhan et al., 2009; Koubaa et al., 2017; Korkolis et al., 2018; Greco et al., 2018; Baral et al., 2018). On the other hand, we are aware that Hill's model, which is quadratic in stresses and contains only six independent parameters to describe the state of anisotropy (see equation (5)), cannot capture complex anisotropic behaviors observed for some materials (Hu, 2007; Banabic, 2010). Thus, in future works, we will extend the three approaches developed in this paper to consider more advanced yield criteria speci cally developed to model the mechanical response of highly textured metallic materials which exhibit complex plastic anisotropy (e.g. the models of Kara llis and Boyce (1993), Bron and Besson (2004), Cazacu et al. (2006) and Aretz and Barlat (2013)). 2.2. Materials investigated We investigate 5 di erent materials denoted as Material 1 ;:::;Material 5. Materials 1 and 2 are model materials with elastic properties and initial density representative of steel, and anisotropy parameters and parameters of the yield stress in the rolling direction that have been speci cally tailored to bring to light the roles of inertia and orthotropy in dynamic formability. Materials 3, 4 and 5 are actual materials, namely, TRIP-780 steel, aluminium alloy 5182-O and aluminium alloy 6016-T4. These three materials are widely used in di erent industrial sectors and engineering applications, and their sheet products are often used in forming operations (Butuc et al., 2002; Li et al., 2010; Ju et al., 2015; Zhou et al., 2018). The values of the parameters corresponding to Materials 1 ;:::;5 are shown in Table 1. Note that all ve materials show full in-plane anisotropy and through thickness isotropy. We are aware that using the simple power-law type relation given in equation (6) to describe the yield stress of the materials investigated in this work is a strong assumption (see the works of Zerilli and Armstrong (1987), Nemat-Nasser and Li (1998) and Rusinek and Klepaczko (2001) where more advanced constitutive relations are developed), especially for TRIP-780 steel, for which the martensitic transformation that occurs at large strains may lead to stress-strain curves with sigmoidal shape (see the works of Tomita and Iwamoto (1995), Iwamoto et al. (1998) and Choi et al. (2010)), which cannot be modeled by equation (6). Nevertheless, using a simple power-law type relation limits the number of material parameters, which facilitates to perform the parametric analysis carried out in Section 6, and it helps in the interpretation of results. Fig. 1 shows the normalized plane stress theoretical yield loci for biaxial loading conditions, I=Yversus 10Symbol Property and units Material 1 Material 2Material 3 TRIP-780 steelMaterial 4 AA 5182-OMaterial 5 AA 6016-T4 0Initial density (kg =m3) 7850 7850 7850 2700 2700 G Elastic shear modulus (GPa), Eq. (3) 75.19 75.19 75.19 25.94 25.94 K Bulk modulus (GPa), Eq. (3) 196.1 196.1 196.1 67.65 67.65 Cp Speci c heat (J =Kg K), Eq. (10) 460 460 460 910 910 Taylor-Quinney coecient, Eq. (10) 0, 0.9 0 0 0 0 0Yield stress parameter (MPa), Eq. (6) 300 300 1351 585.2 417.85 "0Yield stress parameter, Eq. (6) 0.023 0.023 0.007 0.005 0.01 n Yield stress parameter, Eq. (6) 0.12 0.12 0.169 0.3232 0.245 m Yield stress parameter, Eq. (6) 0, 0.01, 0.02 0 0 0 0 _"ref Yield stress parameter (1 =s), Eq. (6) 1000 1000 1000 1000 1000  Yield stress parameter, Eq. (6) 0.5 0 0 0 0 Tref Yield stress parameter (K), Eq. (6) 300 300 300 300 300 F Orthotropy parameter, Eqs. (5) and (9) 0.276 0.208 0.391 0.472 0.729 G Orthotropy parameter, Eqs. (5) and (9) 0.332 0.625 0.532 0.609 0.556 H Orthotropy parameter, Eqs. (5) and (9) 0.668 0.375 0.468 0.391 0.444 N Orthotropy parameter, Eqs. (5) and (9) 1.229 1.250 1.386 1.663 1.194 L Orthotropy parameter, Eqs. (5) and (9) 1.5 1.5 1.5 1.5 1.5 M Orthotropy parameter, Eqs. (5) and (9) 1.5 1.5 1.5 1.5 1.5 Table 1: Values of the parameters for Material 1 ;:::;Material 5. Numerical values of initial density, elastic constants, physical constants, parameters of the yield stress in the rolling direction and plastic orthotropy parameters. Initial density, elastic properties and physical constants of Materials 1 and 2 ( model materials ) are representative of steel and the parameters of the yield stress in the rolling direction and the plastic orthotropy parameters have been speci cally tailored to bring out the roles of inertia and orthotropy in dynamic formability (see Section 6). The parameters of Materials 3, 4 and 5 which correspond to TRIP-780 steel, AA 5182-O and AA 6016-T4, respectively, are taken from Panich et al. (2013), Banabic et al. (2004) and Butuc et al. (2003). 11 II=Y, whereIandIIare the components of the Cauchy stress tensor in the rolling and transverse directions, respectively, for the ve materials considered in this paper. We only show the tension-tension quadrant because all the calculations presented in Section 6 correspond to I>0 andII>0. Speci cally, we investigate the loading cases that lie within the yellow markers that correspond to plane strain tension. Fig. 1a compares the results for Material 1 and Material 2 with the von Mises material. For Material 1 the yield locus is exterior to the von Mises ellipsoid, showing large di erences near equibiaxial tension which are gradually reduced as the I=Y= 0 axis is approached. The curvature of the yield locus of Material 1 near equibiaxial tension is signi cantly greater than for the von Mises material. For Material 2, the yield locus intersects to the von Mises ellipsoid, so that the locus of Material 2 turns from being slightly interior to the von Mises ellipsoid to being exterior and showing important di erences in the sector II> I. Note also that the yield loci of Materials 1 and 2 intersect at I=Y0:87, such that for I=Y&0:87 the yield locus of Material 1 is exterior to Material 2, and vice versa for I=Y.0:87. Fig. 1b compares the results for TRIP-780 steel, AA 5182-O and AA 6016-T4 with the von Mises material. The yield loci of TRIP-780 steel and AA 6016-T4 are exterior and interior to the von Mises ellipsoid, respectively, with increasing di erences near I=Y= 0. Moreover, the yield loci of von Mises material and AA 5182-O intersect at I=Y0:72, so that for greater/smaller values of I=Ythe yield locus of AA 5182-O is interior/exterior to the von Mises ellipsoid. Note also that the results obtained for TRIP-780 steel and AA 5182-O coincide for I=Y= 0 andII=Y= 0, yet the yield locus of TRIP-780 steel is exterior for any other stress state. Fig. 2 shows the evolution of the Lankford coecient (ratio between in-plane transverse and through-thickness plastic strains under uniaxial tension) versus loading orientation for the ve materials investigated in this paper. Note that for isotropic materials the Lankford coecient is constant r= 1. For materials described with the Hill (1948) yield criterion, the value of ris a nonlinear function of which shows maxima/minima for = 0and 90. In addition, for speci c combinations of the anisotropy parameters, materials modeled with Hill (1948) yield criterion show an additional maximum/minimum for an intermediate value of the loading orientation (Dasappa et al., 2012; Huh et al., 2013; Shen et al., 2019; Cazacu and Rodr guez-Mart nez, 2019). Namely, for Material 1 the Lankford coecient is greater than 1, with absolute minimum and maximum for 42and 90, respectively. For Material 2 and TRIP-780 steel the Lankford coecient increases with the loading orientation, with rsmaller than 1 for <45and greater than 1 for  >45. Note that Material 2 is comparatively more anisotropic than TRIP-780 steel because the variation of the Lankford coecient with is more important. For AA 5182-O the 12 00.20.40.60.811.21.41.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6MISES MATERIAL1 MATERIAL 2 PσII / σY σ I / σ Yσ I=σ II von Mises Material 1 Material 2Plane strain tension (a) 00.20.40.60.811.21.41.6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6MISES TRIP-780 AA5182-O AL6016-T4 PσII / σY σ I=σ II σ I / σ YPlane strain tension von Mises Material 3: TRIP-780 steel Material 4: AA 5182-O Material 5: AA 6016-T4 (b) Figure 1: Normalized plane stress theoretical yield loci for biaxial loading conditions, II=YversusI=Y, according to Hill (1948) criterion. (a) Comparison between Material 1, Material 2 and von Mises material. (b) Comparison between Material 3: TRIP-780 steel, Material 4: AA 5182-O, Material 5: AA 6016-T4 and von Mises material. In Section 6 we investigate the loading cases that lie within the yellow markers that correspond to plane strain tension. value ofrincreases, reaches a maximum for = 50, and then decreases. The Lankford coecient is slightly greater than 1 within the range 4062, and smaller than 1 for any other loading orientation. The value ofrfor= 90is greater than for = 0. For AA 6016-T4 the value of ris smaller than 1 for any , reaching the minimum for = 47, with the Lankford coecient for = 0being greater than for = 90. Sowerby and Duncan (1971) and Parmar and Mellor (1978) among others noticed that for Hill (1948) materials there is a strong coupling between the Lankford coecient and the material formability. The general trend being that for r>1 the elliptical plane stress yield locus stretches out (increasing the maximum curvature of the yield surface, see Material 2 in Fig. 1a) and the material formability decreases (with respect to the isotropic von Mises (1928) material), and viceversa for r<1 (see Section 6). The di erences in the plane stress yield loci and the Lankford coecient shown by the 5 materials studied will be used to explain the results presented in Section 6, and to determine the role of material orthotropy in dynamic formability. 13 00.511.522.53 0 15 30 45 60 75 90MATERIAL1 MATERIAL 2 TRIP-780 AA5182-O AL6016-T4 Loading direction, θ ο)Lankford coefficient, rθTransverse direction Rolling direction Material 1 Material 2 Material 3: TRIP-780 steel Material 4: AA 5182-O Material 5: AA 6016-T4 Figure 2: Lankford coecient rversus loading orientation . Comparison between Material 1, Material 2, Material 3: TRIP-780 steel, Material 4: AA 5182-O and Material 5: AA 6016-T4. 3. Linear stability model The linear stability model is an extension of the 2D approach developed by Zaera et al. (2015) to orthotropic materials with yielding based on Hill (1948) criterion. We consider a thin plate of initial thickness h0and edges of initial length L0 xandL0 y, see Fig. 3. The plate is subjected to constant and opposed stretching velocities Vx=_"0 xxL0 x=2 andVy=_"0 yyL0 y=2 on opposed sides, where _ "0 xxand _"0 yyare 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. The Lagrangian Cartesian coordinate system associated to the applied velocity eld is denoted by ( X;Y;Z ) while the corresponding Eulerian counterpart is denoted by ( x;y;z ). Note that the direction Zof the plate is assumed of plane stress. The plate displays orthotropic symmetry with the Cartesian frame of orthotropy ( X1;X2;X3) de ned by an angle betweenXand X1(see Fig. 3), with X1,X2andX3being the rolling, transverse and through-thickness directions, respectively. Our attention is limited to the scenarios where the orthotropy axes and the loading axes are co-directional ( = 0 or 90), i.e. the principal directions of stress and strain coincide. In this paper the angle is also referred to as material orientation. Notice that varying the material orientation from = 0to = 90amounts to exchange the value of the anisotropy coecients FandG, see equations (5) and (9). The elastic deformations are neglected and hereafter  "= "pwill be named e ective strain and _"=_"pe ective strain rate. 14 Figure 3: Linear stability analysis. Schematic representation of the geometry and boundary conditions of the problem addressed: a plate of initial thickness h0and edges of initial length L0 xandL0 ysubjected to constant and opposed stretching velocities Vx=_"0 xxL0 x=2 andVy=_"0 yyL0 y=2 on opposed sides, where _ "0 xxand _"0 yyare the imposed initial strain rates. The Lagrangian Cartesian coordinate system associated to the applied velocity eld is denoted by ( X;Y;Z ), while the corresponding Eulerian counterpart is denoted by (x;y;z ). Note that the direction Zof the plate is assumed of plane stress. The plate displays orthotropic symmetry with the Cartesian frame of orthotropy ( X1;X2;X3) de ned by an angle betweenXandX1. The unknowns of the problem are: velocities ( Vx;Vy), strain rates ( dxx;dyy), deformation gradients ( Fxx;Fyy;Fzz), thickness (h), e ective strain ( "), e ective strain rate ( _"), temperature ( T), uniform stresses ( xx;yy), average stress in the major stretching direction ( avg xx), e ective stress ( ), yield stress in the rolling direction ( Y) and Bridgman's correction factor argument ( '). The fundamental equations governing the loading process are presented below: Kinematic relations: dxx=@Vx @x dyy=@Vy @y(13) 15 Fxx=_Fxx dxx Fyy=_Fyy dyy(14) Incompressibility condition: FxxFyyFzz= 1 (15) Current thickness: h=h0Fzz (16) Momentum balance: Fyy@(havg xx) @X=h0@Vx @t Fxx@(hyy) @Y=h0@Vy @t(17) whereis the material density (assumed constant, i.e. =0). The momentum balance equation in the out of plane direction has been disregarded, see Zaera et al. (2015). Average stress in the major stretching direction: avg xx=B(')xx (18) where 16 B(') =p 1 +'1ln 1 + 2'+ 2p '(1 +') 1 (19) is the correction factor introduced by Bridgman (1952) to model the hydrostatic stresses that developed in a necked section (see Zhou et al., 2006; Zaera et al., 2015). Notice that, while in this paper we have used the classical Bridgman (1952) approximation, in future works we shall consider the one-dimensional model based on an asymptotic analysis recently developed by Audoly and Hutchinson (2016, 2019) to determine the stress state inside the neck. Argument of Bridgman's correction factor: '=h 8@2h @x2(20) E ective stress: = (G+H)2 11+ (H+F)2 222H1122+ 2N2 121=2(21) Previous expression is a specialization of equation (5) for plane stress conditions, where 11,22and12are given by: 11=xxcos2 +yysin2 22=xxsin2 +yycos2 12= (xxyy) sin cos (22) Flow rule: 17 dxx=_"@ @xx dyy=_"@ @yy(23) where _"=Fd2 11+Gd2 22+H(d11+d22)2 FG+FH+GH+2d2 12 N1=2 (24) Previous expression is a particularization of equation (9) for plane stress conditions, where d11,d22andd12 are given by: d11=dxxcos2 +dyysin2 d22=dxxsin2 +dyycos2 d12= (dxxdyy) sin cos (25) The governing equations are completed with the expressions of the yield condition, the yield stress in the rolling direction, the e ective strain and the temperature evolution equation, given by the relations (4), (6), (7) and (10), respectively. The solution of the problem at any given loading time, usually referred to as homogeneous or fundamental solution: S1= (V1 x;V1 y;d1 xx;d1 yy;F1 xx;F1 yy;F1 zz;h1;_"1;"1;T1;1 xx;1 yy;avg;1 xx;1;1 Y;'1)T(26) can be readily obtained by integration of the governing equations with the following initial conditions: 18 Vx(X;Y; 0) = _"0 xxX;Vy(X;Y; 0) = _"0 yyY;ij(X;Y; 0) = 0 "(X;Y; 0) = 0; h(X;Y; 0) =h0;T(X;Y; 0) =T0=Tref(27) and boundary conditions: Vx L0 x=2;Y;t =Vx L0 x=2;Y;t = _"0 xxL0 x=2 Vy X;L0 y=2;t =Vy X;L0 y=2;t = _"0 yyL0 y=2(28) where, as already mentioned, _ "0 xxand _"0 yyde ne the initial strain rates in XandYdirections of the plate, respectively. The applied loading is selected through input values of _ "0 xxand = _"0 yy=_"0 xx. Note that _ "0 xxwill be also referred to as imposed initial major strain rate. The details on the kinematics of the fundamental solution are given in Appendix A. The stability of the homogeneous solution is tested by introducing at some time t1a small perturbation Sof the form: S(X;t)t1=S1eiX+(tt1)(29) where S1= (Vx;Vy;dxx;dyy;Fxx;Fyy;Fzz;h; _";";T;xx;yy;avg xx;;Y;')T(30) is the perturbation amplitude, is the wavenumber (also called perturbation mode) and is the growth rate of the perturbation at time t1. The perturbation is imposed on lines X=constant which is the orientation naturally selected by the material to trigger a neck for the loading paths (0  1) and the orientations of the orthotropy axes ( = 0or 90) investigated in this work (see Zaera et al., 2015). Notice that, while in this paper we have assumed the classical frozen coecients theory and the exponential form of the time dependency of the perturbation for the sake of simplicity, Xavier et al. (2020) have recently developed a promising method to relax these assumptions. 19 The perturbed solution of the problem is: S=S1+S (31) withjSj S1 . A system of linear algebraic equations is obtained by substituting equation (31) into the governing equations and keeping only the rst-order terms in the increments S1. This system of equations, which is not explicitly shown in this paper for the sake of brevity, has been adimensionalized by introducing the following dimensionless groups: ^Vx=Vx h0_"1; ^Vy=Vy h0_"1; ^dxx=dxx _"1; ^dxx=dxx _"1; ^Fxx=Fxx; ^Fyy=Fyy ^Fzz=Fzz; ^h=h h0;^_"=_" _"1; ^"= "; ^T=T T0; ^xx=xx 0 ^yy=yy 0; ^avg xx=avg xx 0; ^= 0; ^Y=Y 0; ^'=' ^= _"1; ^=h0(32) A non trivial solution for the adimensionalized perturbation amplitude ^S1is obtained only if the determinant of the matrix of coecients of the system of adimensionalized linear algebraic equations is equal to zero. Application of this condition leads to a third order polynomial in ^ : B3(^S1)^3+B2(^S1;^)^2+B1(^S1;^)^+B0(^S1;^) = 0 (33) with 20 B3(^S1) = 8~c(^F1 xx)3 ^M1(^d1 xx)22^R2^d1 xx^d1 yy+^R1(^d1 yy)2+^Q2 ^M1^R1^R2 2 B2(^S1;^) =8(^F1 xx)2 ^F1 xx ^Q3 1 ^R2 2^M1^R1 +^R1(^d1 yy)2 + ~c^Q1 ^R2 2^M1^R1 + ^M1^Q3(^d1 xx)22 ^Q3^R2^d1 xx^d1 yy ^2~c~H2^h1^F1 yy (^d1 yy)2+^M1^Q2 B1(^S1;^) =2 3^2~H2^h1^F1 yy 12(^F1 xx)2 ^Q3 (^d1 yy)2^M11 + ~c ^d1 xx (^d1 yy)2+^M1^Q2 ^M1^Q1  ~c^1 xx 12(^F1 xx)2^2(^h1)2 ^M1(^d1 xx)22^R2(^d1 xx)2(^d1 yy)2+^R1(^d1 yy)2+^Q2 ^M1^R1^R2 2 B0(^S1;^) =2 3^2~H2^h1^F1 yy 12^d1 xx(^F1 xx)2 ~c^M1^Q1 ^Q3 (^d1 yy)2^M11 ^1 xx 12(^F1 xx)2^2(^h1)2  ^Q3 1 ^R2 2^M1^R1 +^R1(^d1 yy)2 + ~c^Q1 ^R2 2^M1^R1 + ^M1^Q3(^d1 xx)22 ^Q3^R2^d1 xx^d1 yy(34) where ^R1=0@2 @2xx t1,^R2=0@2 @xx@yy t1and ^M1=0@2 @2yy t1are the non-dimensional second-derivatives of the e ective stress with respect to the Cauchy stress tensor components, ^Q1=1 0@Y @" t1,^Q2=_"1 0@Y @_" t1 and ^Q3=T1 0@Y @T t1are the non-dimensional strain, strain-rate and temperature sensitivities of the material, respectively, ~ c=CpT0 0is the non-dimensional speci c heat, and1 ~H2= h0_"12 0represents the inertial resistance to motion. Equation (33) has three roots in ^ which depend on the adimensionalized fundamental solution of the problem ^S1and the adimensionalized wavenumber ^. Note also that the adimensionalized wavenumber is related to the normalized perturbation wavelength L0=h0as:L0=h0=2 ^. The requisite for unstable growth of the perturbation is given by<(^)>0. Thus, the root with the greatest positive real part denoted by ^ +represents the physical solution of the problem (Dudzinski and Molinari, 1991). Note that ^ +is usually referred to as instantaneous instability index (El Ma  et al., 2014). Following Fressengeas and Molinari (1994), we introduce a cumulative instability index de ned as I=Rt 0^+_"1d which tracks the history of the instantaneous grow rate of all the growing modes. To calculate I, we introduce the perturbation at di erent times and sum the instantaneous growth rate obtained for each loading time. The stabilizing e ect of stress multiaxiality and inertia on small and large perturbation wavelengths, respectively, boosts the growth of a nite number of intermediate perturbation modes (Fressengeas and Molinari, 1994). The wavelength of the mode with the greatest value of the cumulative instability index Iis called the critical cumulative wavelength, and the corresponding value of Iis called the critical 21 cumulative instability index. The critical cumulative instability index, which evolves with the loading time, and thus with the strain, will be used in Section 6 to determine when a perturbation mode turns into a necking mode, allowing to perform quantitative comparisons between the results obtained with the linear stability analysis set out in this section, and the nonlinear two-zone model and the nite element model presented in Sections 4 and 5, respectively. Notice that Dudzinski and Molinari (1991) and Zaera et al. (2015) relied on the instantaneous instability index ^ +to determine when a neck is formed, in order to construct forming limit diagrams. However, in this paper we choose Ibecause, as discussed by El Ma  et al. (2014) and Vaz-Romero et al. (2017), it provides predictions for necking formation which are generally closer to nite element results. 4. Nonlinear two-zone model The nonlinear two-zone model is an extension to orthotropic materials with yielding based on Hill (1948) crite- rion of the formulation recently developed by Jacques (2020) to include inertia e ects in the classical localization approach of Marciniak and Kuczy nski (1967). The model consists of a plate idealized as an array of unit-cells with spatial imperfections in the form of bands with reduced initial thickness. Fig. 4 shows the schematic representation of the unit-cell. Following the notation used in Section 3, the Lagrangian and Eulerian Cartesian coordinate sys- tems associated to the applied loading conditions are ( X;Y;Z ) and (x;y;z ), respectively. Similarly, ( X1;X2;X3) is the Cartesian coordinate system associated to the orthotropy axes de ned by an angle betweenXandX1. The initial lengths of the cell are L0 xandL0 y, and the zones inside and outside the band are denoted by AandB, respectively. The initial thicknesses of zones AandBare: hA;0=h0(1) hB;0=h0(35) where  is the normalized amplitude of the imperfection. Note that superscripts ()iwithi=A;B are used to denote quantities in the zones AandB, respectively. The initial and boundary conditions are consistent with those used in the linear stability model (see Fig. 3, and equations (27) and (28)). Note that, in the present work, the imperfection (band with reduced thickness) is assumed perpendicular to the major stretching direction X. Indeed, as earlier stated in Section 3, this is the orientation naturally selected by the material to trigger a neck for the loading paths (0  1) and material 22 orientations ( = 0or 90) investigated in this paper. However, this assumption can be easily relaxed to consider other material orientations. Note also that, similarly to the stability model (see Section 3), the direction Zof the cell is assumed of plane stress and elastic deformations are neglected. Figure 4: Two-zone model. Schematic representation of the unit-cell with two zones of uniform thickness AandB. The initial lengths of the cell are L0 xandL0 y, and the initial thicknesses of zones AandBarehA;0andhB;0, respectively, with hB;0>hA;0. The initial and boundary conditions are consistent with those used in the linear stability model, see Fig. 3. The zonesAandBare assumed to be regions of homogeneous deformation and the unknowns of the problem in these two zones are: strain rates ( di xx;di yy), deformation gradients ( Fi xx;Fi yy;Fi zz), thicknesses ( hi), e ective strains ("i), e ective strain rates ( _"i), temperature ( Ti), uniform stresses ( i xx;i yy), average stress in zone Ain the major stretching direction ( A;avg xx ), yield stresses in rolling direction ( i Y) and Bridgman's correction factor argument in zoneA('A). The equations governing the loading process are presented below: Kinematic relations: The components of the overall deformation gradient in the whole unit-cell ( Fxx;Fyy;Fzz) are given by equation (A.7) and the continuity of the displacement between the zones AandBleads to the following relations 23 linking the deformation gradient components in the zones AandB, and the overall deformation gradient components: FA xxR+FB xx(1R) =Fxx FA yy=FB yy=Fyy(36) where R=LA;0 x L0x(37) is the ratio between the initial length of the imperfection, LA;0 x, and the initial length of the cell in the major stretching direction, L0 x. The strain rate components in the zones AandBare given by: di xx=_Fi xx Fixx di yy=_Fi yy Fiyy(38) Virtual work principle: Jacques (2020) stated that including inertia e ects in the classical localization approach of Marciniak and Kuczy nski (1967) implies that the momentum balance equation cannot be ful lled at any point of the cell, and he circumvented this diculty using the virtual principle work (i.e. the weak form of the momentum balance equation), which leads to the following di erential equation to describe the evolution of FA xx: FA xx=1 Cd hB;0B xx FBxxhA;0A;avg xx FAxx! (39) 24 with Cd=(L0 x)2 12 R2hA;0+R(1R)hB;0 (40) Note that all the unknowns of the problem can be readily computed once FA xxis obtained. The required steps to derive equation (39) are given in Jacques (2020). Using the plane stress condition, i zz= 0, the Cauchy stress tensor components in zones AandBare computed as: i xx=si xxsi zz i yy=si yysi zz(41) where si xx;si yy;si zz are the components of the stress deviator referred to the coordinate system associated to the loading axes. Moreover, with respect to the orthotropy axes, they can be computed as: si 11=i Y _"iF (FH+FG+HG)di 11 si 22=i Y _"iF (FH+FG+HG)di 22 si 33=i Y _"iF (FH+FG+HG)di 33 si 12=i Y _"idi 12 N(42) where the relationship between si xx;si yy;si zz and si 11;si 22;si 33;si 12 is calculated using the corresponding rotation matrix, as in equation (22). These expressions have been obtained using equations (8) and (21). Note that the expression of the shear component si 12is given for completeness but, for the material orientations considered in the present study ( = 0and = 90),si 12is equal to zero. The governing equations are completed with the expressions of the yield stress in the rolling direction, the 25 e ective strain, the temperature evolution equation, the incompressibility condition, the current thicknesses, the average stress in the major stretching direction in zone Aand the e ective strain rate. These expressions are given by equations (6), (7), (10), (15), (16), (18) and (24), respectively. Note that the argument of the Bridgman correction factor is calculated using the following expression (see Jacques (2020) for more details): 'A=2hA hBhA 16(LA;0 xFAxx)2(43) The equations governing the loading process have been integrated using the explicit Euler method. In the calculations we consider that L0 x=L0 y=L0. The results of the model are obtained for several normalized cell sizes 0:25L0=h010 andh0 xed (h0= 2 mm). The normalized cell size in the two-zone model is the counterpart of the normalized perturbation wavelength in the linear stability analysis, see Section 3, and thus it will be denoted the same L0=h0. The normalized amplitude of the imperfection is  = 0 :2% in most two-zone calculations presented in this work (unless otherwise noted). This value is arbitrary and we take it because the same imperfection amplitude was used in the nite element calculations performed by N'souglo et al. (2020). The speci c in uence of the imperfection amplitude in the two-zone model results is investigated in Section 6.5. Moreover, the two-zone model results are dependent on the value of R(see equation (37)). In the case of isotropic materials, Jacques (2020) found that the value R= 0:28 gives the best agreement between the nite element simulations and the dynamic two-zone model, and therefore the same value is used in the present study. It will be seen in Section 6 that this value also yields an excellent agreement between the two-zone model and the nite element computations for the di erent anisotropic materials and within the whole range of loading conditions considered. 5. Finite element model The nite element model is based on the unit-cell model developed by Xue et al. (2008), and later adopted by Rodr guez-Mart nez et al. (2017) and N'souglo et al. (2020), to study the in uence of constitutive behavior and inertia on neck formation in thin metal plates subjected to dynamic loading. The model consists of a plate modeled as an array of unit-cells with sinusoidal spatial imperfections de ned by the following expression: h=h0 2 2 1 + cos2X L0x (44) 26 whereis the amplitude of the imperfection. Due to the symmetry of the model, only one eight of the cell has been analyzed (see Fig. 5), with reference con guration (imperfection-free) given by the domain 0 XL0 x=2, 0YL0 y=2 and 0Zh0=2, withL0 x=L0 y=L0. The origin of coordinates of the Lagrangian Cartesian coordinate system ( X;Y;Z ) is located at the center of mass of the whole cell. Similarly to the dynamic two- zone model (see Section 4), we have investigated several normalized cell sizes 0 :25L0=h010 withh0 xed (h0= 2 mm). As in the two-zone model, the normalized imperfection amplitude is  =2 h0= 0:2% in most of the nite element calculations (unless otherwise noted). Nevertheless, the role of  in the nite element results is discussed in Section 6.5. The nite element calculations are performed with ABAQUS/Explicit (2016) under the following imposed initial conditions: Vx(X;Y;Z; 0) = _"0 xxX Vy(X;Y;Z; 0) = _"0 yyY Vz(X;Y;Z; 0) = _"0 xx+ _"0 yy Z T(X;Y;Z; 0) =T0=Tref(45) and boundary conditions: Ux(0;Y;Z;t ) = 0; Uy(X;0;Z;t) = 0; Uz(X;Y; 0;t) = 0 Vx L0 x=2;Y;Z;t = _"0 xxL0 x=2;Vy X;L0 y=2;Z;t = _"0 yyL0 y=2(46) consistent with those used in the theoretical approaches (see Sections 3 and 4) and the symmetry imposed to the nite element model (note that Ux,UyandUzdenote displacements along the X,YandZdirections, respectively). As shown by Rodr guez-Mart nez et al. (2017), these initial conditions allow to minimize the propagation of stress waves within the unit-cell during the loading process. Similarly to the theoretical models (see Sections 3 and 4), the nite element calculations have been carried out for loading paths ranging between plane strain stretching and equibiaxial stretching (0  1), namely, we have performed calculations for = 0, 0:125, 0:25, 0:375, 0:5, 0:625, 0:75 and 1. Each marker in Figs. 7, 8, 9, 10, 12, 13, 14 and C.16 corresponds to a value of . Notice that unit-cell calculations with ab initio geometric imperfections do not allow for exploring loading paths such that <0. The problem is that the orientation that has to be assigned to the imperfection such that it is aligned with the direction of zero stretch when necking occurs is not known in advance (see Hill (1952)). Moreover, the major stretching 27 direction of the cell, which coincides with the Xaxis (see Fig. 5), has been aligned either with the rolling direction ( = 0in the theoretical models) or the transverse direction ( = 90in the theoretical models). Figure 5: Finite element model. Mesh and boundary conditions for a normalized cell size of L0=h0= 1. A large imperfection amplitude has been shown in the gure for better illustration of the geometric perturbation. For interpretation of the