6 In this paper, we have modi ed the stress integration scheme proposed by Choi and Yoon (2019), which is 7 based on the numerical approximation of the yield function gradients, to implement in the nite element code 8 ABAQUS three elastic isotropic, plastic anisotropic constitutive models with yielding described by Yld2004- 9 18p (Barlat et al., 2005), CPB06ex2 (Plunkett et al., 2008) and Yld2011-27p (Aretz and Barlat, 2013) criteria, 10 respectively. We have developed both VUMAT and UMAT subroutines for the three constitutive models, and 11 have carried out cylindrical cup deep drawing test simulations and calculations of dynamic necking localization 12 under plane strain tension, using explicit and implicit analyses. An original feature of this paper is that these 13 nite element simulations are systematically compared with additional calculations performed using (i) the 14 numerical approximation scheme developed by Choi and Yoon (2019), and (ii) the analytical computation of the 15 rst and second order yield functions gradients. This comparison has shown that the numerical approximation 16 of the yield function gradients proposed in this paper facilitates the implementation of the constitutive models, 17 and in the case of the implicit analyses, it leads to a signi cant decrease of the computational time without 18 impairing the accuracy of the nite element results. In addition, we have demonstrated that there is a critical 19 loading rate below which the dynamic implicit analyses are computationally more ecient than the explicit 20 calculations. 21 Keywords: 22 Advanced yield criteria, Anisotropy, Stress integration, Finite elements, Metal forming 23 1. Introduction 24 The pioneering work of Hill (1948) was a turning point in the mathematical theory of anisotropic plasticity. 25 By analogy with the von Mises (1928) yield criterion for isotropic materials, Hill (1948) proposed a quadratic 26 yield function to model the macroscopic behaviour of orthotropic metals which depends only on the deviatoric 27 stresses and is pressure independent. The model highlights due to its simplicity, versatility and easy calibration 28 Corresponding author. Tel. +34 916249904; Fax: +33 916249430. E-mail address:jarmarti@ing.uc3m.es Preprint submitted to Finite Elements in Analysis and Design February 23, 2021 2 procedure; so that 70 years after its publication, the criterion of Hill (1948) still is the anisotropic model 29 most widely used to describe the directionality of plastic properties of metallic materials used in scienti c 30 and technological applications (Padmanabhan et al., 2009; Koubaa et al., 2017; Korkolis et al., 2018; Greco 31 et al., 2018; Baral et al., 2018; Cazacu and Rodr guez-Mart nez, 2019). However, due to the fact that Hill's 32 yield function is quadratic in stresses and contains only six independent parameters to describe the state of 33 anisotropy, it cannot capture with the same accuracy both the yield stress and Lankford coecients (Hu, 34 2007; Banabic, 2010). Hill (1979) acknowledged what he referred to as inevitable limitations in the range of 35 validity of the 1948 prototype and proposed an improved yield criterion based on arbitrary linear combinations 36 of equal non-integer powers of the principal values of the deviatoric stress, together with the same powers of 37 their di erences. This criterion, speci cally developed to describe the mechanical response of textured sheets 38 for situations in which the orthotropy and the loading are co-directional (Hill, 1990; Tong, 2016), included 7 39 parameters to describe material anisotropy, and provided extended exibility without sacri cing too much of 40 the simplicity that was a main attraction of the 1948 proposal. Nevertheless, nor the original model, neither the 41 1979 improvement can describe the mechanical response of materials which exhibit complex plastic anisotropy. 42 Many yield functions have been developed over the years to improve the formulations of Hill (1948, 1979). 43 For instance, Barlat and Lian (1989) extended to orthotropy the 2D model of Barlat and Richmond (1987) using 44 a linear transformation of the components of the stress tensor. The anisotropic criterion of Barlat and Lian 45 (1989) captured the behavior of orthotropic sheet metals under a full plane stress state including 4 independent 46 material parameters, and showed remarkable capacity to describe the yield surfaces computed with the Bishop 47 and Hill (1951) polycrystal model. Shortly after, Barlat et al. (1991) extended to orthotropic anisotropy the 48 isotropic criterion proposed by Hosford (1972) in order to generalize the 2D model of Barlat and Lian (1989) 49 to any complex three-dimensional stress state. However, the criterion of Barlat et al. (1991) cannot be used 50 for materials with symmetries other than the orthotropic symmetry. This motivated to Kara llis and Boyce 51 (1993) to propose a non-quadratic yield criterion capable of describing di erent states of material symmetry, 52 including the most generalized case of the totally asymmetric material. Following previous ideas of Barlat and 53 Lian (1989) and Barlat et al. (1991), among others, Kara llis and Boyce (1993) used a fourth order tensor as 54 a linear multiplicative operator acting on the stress tensor to introduce material anisotropy. The predictions 55 of the model of Kara llis and Boyce (1993) were found to be in good agreement with experimental results 56 for a variety of materials. The fourth order tensorial operator introduced by Kara llis and Boyce (1993) was 57 further used by Barlat et al. (1997) to improve the model developed in Barlat et al. (1991), introducing three 58 additional material parameters which provided increased exibility to describe the anisotropic response of 59 3 di erent materials. The model of Barlat et al. (1997) was demonstrated to be particularly valid for materials 60 which exhibit large pure shear yield stress. 61 Following the works of Kara llis and Boyce (1993) and Barlat et al. (1997), the linear transformation-based 62 anisotropic yield functions became increasingly popular because, among other reasons, they allow tailoring the 63 number of anisotropy parameters without signi cant impact on the mathematical structure of the constitutive 64 model. Barlat et al. (2003) developed a plane stress yield function with 9 parameters to describe the mechanical 65 response of aluminum sheets using two linear transformations of the Cauchy stress tensor. It was shown that 66 the accuracy in the description of the ow stress anisotropy and Lankford coecients is increased as the number 67 of linear transformations, and thus the number of anisotropy parameters, increases. Bron and Besson (2004) 68 extended the models of Barlat et al. (1991) and Kara llis and Boyce (1993) using two linear transformations. 69 The resulting 3D yield function, that was demonstrated to be convex, contained 16 parameters and provided 70 accurate description of the plastic anisotropy of various aluminum sheet samples. Two linear transformations 71 were also used by Barlat et al. (2005) to formulate the so-called Yld2004-13p and Yld2004-18p yield functions. 72 These two anisotropic models were derived from di erent isotropic criteria, so that they contain 13 and 18 73 anisotropic parameters, respectively. The values of these parameters were computed for 2090-T3 and 6111- 74 T4 aluminum alloys. The comparison between yield functions predictions and experiments showed that both 75 models capture the main features of the anisotropic behaviour of 2090-T3 and 6111-T4. Notably, the Yld2004- 76 18p model described very accurately the ow stresses in uniaxial tension and the Lankford coecients of both 77 sheet materials. Two and three linear transformations were used by Aretz and Barlat (2013) to obtain the 78 so-called Yld2011-18p and Yld2011-27p yield functions, which contain 18 and 27 parameters, respectively. 79 While the Yld2011-18p model can be regarded as complementary to the yield function Yld2004-18p proposed 80 few years earlier (same number of parameters but di erent behavior during calibration of the experimental 81 data), the Yld2011-27p provided extended exibility and improved accuracy to describe the complex plastic 82 anisotropy of AA 2090-T3 and AA 3104-H19. 83 Most of the anisotropic yield criteria cited so far were speci cally developed to model metals with cubic 84 structure. Hence, they are usually represented by even functions of stresses and do not capture the ten- 85 sion/compression asymmetry in yielding characteristic of hexagonal closed packed metals. Cazacu and Barlat 86 (2004) proposed a yield criterion to describe both the asymmetry and anisotropy in yielding of magnesium 87 and magnesium alloys using the generalized invariants approach developed by Cazacu and Barlat (2001). The 88 yield function of Cazacu and Barlat (2004) is an homogeneous function of degree three in stresses and contains 89 18 material parameters. Shortly after, Cazacu et al. (2006) developed the so-called CPB06 yield criterion, 90 4 in which the degree of homogeneity is a material parameter to be determined, to capture simultaneously the 91 anisotropy and tension/compression asymmetry of hexagonal closed packed metals. Similarly to Kara llis and 92 Boyce (1993) and Barlat et al. (1997), the anisotropy was introduced using a linear transformation of the stress 93 deviator. For 3D stress conditions, this yield criterion involves 9 independent anisotropy coecients, 1 param- 94 eter to describe the strength di erential e ect and 1 more to determine the degree of homogeneity of the yield 95 function. Plunkett et al. (2008) extended the CPB06 model using two linear transformations and derived the 96 so-called CPB06ex2 yield criterion which contains a total of 21 parameters, and represents with great accuracy 97 both the tensile and compressive anisotropy in yield stresses and Lankford coecients of materials with hcp 98 crystal structure and of metal sheets with cubic crystal structure. During the last few years, additional yield 99 functions to model the anisotropy and tension/compression asymmetry in yielding of various metallic materials 100 have been proposed, for instance, by Nixon et al. (2010), Khan et al. (2012), Raemy et al. (2017), Kondori 101 et al. (2019), Park et al. (2019), Sanz et al. (2019) and Lee et al. (2020). 102 As a general rule, the ability of the yield functions to describe complex anisotropic and tension/compression 103 asymmetric behaviors increases with the number of material parameters. However, accuracy comes with a price. 104 As the number of parameters increases: (i) larger number of tests is needed to determine the material parame- 105 ters (if experimental data are not available, modelers also rely on numerical data obtained from polycrystalline 106 calculations, see Yoon et al. (2006) and Grytten et al. (2008)), (ii) the implementation of the yield functions 107 in nite element codes becomes an increasingly laborious task and (iii) the computational eciency decreases. 108 In particular, a key step in the nite element implementation process is the computation of the yield function 109 gradients (i.e. the derivatives of the e ective stress with respect to the stress components) which are required 110 to calculate the e ective plastic strain increment and the tangent sti ness tensor, the latter being necessary in 111 implicit nite element solutions (e.g. Yoon et al. (2004) and Vadillo et al. (2016)). For advanced yield criteria, 112 in which anisotropy is included using more than one linear transformation (see previous paragraphs of this in- 113 troductory section), the analytical forms of the rst and second order yield function gradients are lengthy and 114 complex expressions, which increases the diculties associated to the implementation of constitutive models in 115 nite element codes. In order to alleviate this problem, Aretz (2007) proposed to use numerical approximations 116 of the yield function gradients, which is also an e ective way to get rid of the singularities that appear in their 117 analytical counterparts (Barlat et al., 2005; Aretz and Barlat, 2013). Very recently, Choi and Yoon (2019) 118 modi ed the stress integration algorithm developed by Aretz (2007) to implement in ABAQUS/Standard vari- 119 ous anisotropic constitutive models with distortional hardening under associated and non-associated ow rules. 120 In addition, the simple mathematical form of the numerical approximations suggests that their application in 121 5 the nite element implementation process could lead to signi cant computational time reduction, facilitating 122 the utilization of advanced yield functions in large scale numerical simulations. Nevertheless, to the authors' 123 knowledge, the anticipated computational eciency of using numerical approximations of the yield function 124 gradients has not been thoroughly investigated so far. In this regard, note that Choi and Yoon (2019) only 125 performed relatively small simulations: single element calculations to evaluate the accuracy of the numerical 126 approximations of the yield function gradients, and deep drawing calculations where the blank is modeled 127 with a maximum of 2000 elements to check the number of ears predicted by the various constitutive models 128 investigated. 129 In this paper, we modify the numerical di erentiation schemes proposed by Aretz (2007) and Choi and 130 Yoon (2019) to implement in the nite element code ABAQUS three elastic isotropic, plastic anisotropic 131 constitutive models with yielding described by Yld2004-18p, CPB06ex2 and Yld2011-27p criteria, respectively. 132 We have developed VUMAT and UMAT subroutines for the three constitutive models, using both the analytical 133 expressions and the numerical approximations of the rst and second order yield function gradients. We have 134 carried out cylindrical cup deep drawing test simulations using the con guration reported in Yoon et al. 135 (2006), and calculations of dynamic necking localization under plane strain tension with the nite element 136 model developed by N'souglo et al. (2020), using both implicit and explicit analyses. We have identi ed that 137 there is a critical loading rate below which the dynamic implicit analyses are computationally more ecient 138 than the explicit calculations. Moreover, the systematic comparison between the nite element results obtained 139 with the analytical expressions and the numerical approximations of the yield functions gradients has shown 140 that, while in the case of the explicit analyses using a VUMAT the procedure used to calculate the stress 141 derivatives hardly a ects the computational time, the utilization of the numerical approximations in implicit 142 analyses using a UMAT indeed leads to an important decrease of the computational cost. This decrease reaches 143 up to 70% for some of the numerical simulations presented in this paper, and it increases as the number of 144 elements in the nite element mesh increases. We have also shown that the computational time required for the 145 dynamic necking calculations is between 5% and 45% smaller for the numerical di erentiation scheme proposed 146 in this work in comparison with the formulation developed by Choi and Yoon (2019). 147 2. Constitutive framework 148 We consider three elastic isotropic, plastic orthotropic constitutive models with yielding described by 149 Yld2004-18p (Barlat et al., 2005), CPB06ex2 (Plunkett et al., 2008) and Yld2011-27p (Aretz and Barlat, 150 2013) criteria, respectively. These advanced yield criteria were speci cally developed to capture the mechanical 151 6 response of highly textured metallic materials which exhibit complex plastic anisotropy. Section 2.1 sets out 152 the general formulation of the three constitutive models, while the speci c features of each yield criterion are 153 shown in Sections 2.2, 2.3 and 2.4, respectively. 154 155 2.1. Basic features of the constitutive models 156 Elasto-plastic constitutive equations in nite element codes for large deformation solid mechanical applica- 157 tions are mostly based on an ad hoc extension of classical small-strain elasto-plasticity to the nite deformation 158 rang. We assume the additive decomposition of the total rate of deformation tensor dinto an elastic part de 159 and a plastic part dp: 160 d=de+dp (1) where the elastic part of the rate of deformation tensor is related to the rate of the stress by the following 161 linear elastic law: 162 O=L:de (2) whereOis an objective derivative of the Cauchy stress tensor, and Lis the tensor of isotropic elastic moduli 163 given by: 164 L= 2GI0+K1 1 (3) with 1andI0being the unit second-order tensor and the unit deviatoric fourth-order tensor, respectively. 165 MoreoverGandKare the shear modulus and the bulk modulus (they are assumed constant). 166 167 In order to achieve incremental objectivity of the constitutive equations, ABAQUS considers the objective 168 stress derivative to be de ned as: 169 7 O=_+!! (4) where!is a skew-symmetric tensor so thatOcoincides with the Green-Naghdi objective derivative used by 170 ABAQUS/Explicit if != , where =_RRTwithRbeing the polar rotation tensor, and with the Jaumann 171 objective derivative used by ABAQUS/Standard if !=w, wherewis the spin tensor (i.e. the skew-symmetric 172 part of the velocity gradient tensor). Note that_() denotes di erentiation with respect to time. 173 174 The general expression for the yield criteria considered in this paper is: 175 =  ~S0 i;~S00 i;~S000 i (5) where ~S0 i;~S00 i;~S000 i , withi= 1;2;3 (see equation (8)), are the principal values of the second order deviatoric 176 tensors ~s0,~s00and~s000, which are de ned by three linear transformations: 177 ~s0=C0s (6a) ~s00=C00s (6b) ~s000=C000s (6c) 178 where sis the deviatoric part of the Cauchy stress tensor , and C0,C00andC000are the matrices which contain 179 the anisotropy coecients: 180 8 C0=2 666666666666664c0 11c0 12c0 13 0 0 0 c0 21c0 22c0 23 0 0 0 c0 31c0 32c0 33 0 0 0 0 0 0 c0 44 0 0 0 0 0 0 c0 55 0 0 0 0 0 0 c0 663 777777777777775(7a) C00=2 666666666666664c00 11c00 12c00 13 0 0 0 c00 21c00 22c00 23 0 0 0 c00 31c00 32c00 33 0 0 0 0 0 0 c00 44 0 0 0 0 0 0 c00 55 0 0 0 0 0 0 c00 663 777777777777775(7b) C000=2 666666666666664c000 11c000 12c000 13 0 0 0 c000 21c000 22c000 23 0 0 0 c000 31c000 32c000 33 0 0 0 0 0 0 c000 44 0 0 0 0 0 0 c000 55 0 0 0 0 0 0 c000 663 777777777777775(7c) 181 withc0 ij,c00 ijandc000 ijbeing material parameters. Note that for Yld2004-18p and Yld2011-27p yield criteria c0 ii= 182 c00 ii=c000 ii= 0 withi= 1;2;3 (see Sections 2.2 and 2.4), i.e. for these two constitutive models the matrices which 183 contain the anisotropy coecients only have 9 non-zero coecients. In equation (6), we have taken the order of 184 the components of the column vector representing the deviatoric stress tensor to be s=fs11s22s33s23s31s12g, 185 with the Cartesian coordinate system (1 ;2;3) being associated to the orthotropy axes of the material (1, 2 and 186 3 axes correspond to the rolling, transverse and normal directions, respectively). 187 188 The ordered principal values of the tensor ~s0are: 189 9 ~S0 1= 2q (H0 1)2+H0 2cos 3 +H0 1 (8a) ~S0 2= 2q (H0 1)2+H0 2cos+ 4 3 +H0 1 (8b) ~S0 3= 2q (H0 1)2+H0 2cos+ 2 3 +H0 1 (8c) 190 with= arccosq p3=2 ,p= (H0 1)2+H0 2andq=2 (H0 1)3+ 3H0 1H0 2+ 2H0 3 2, whereH0 1,H0 2andH0 3are the 1st, 191 2ndand 3rdinvariants of ~s0. Analogous procedure is followed to obtain the principal values of the tensors ~s00192 and~s000. 193 194 Moreover, the yield condition is expressed as: 195 f= Y= 0 (9) where = ( ) is the e ective stress (the speci c form of  for the three yield criteria considered is given 196 in equations (19), (21) and (25)) and Yis the tensile yield stress in the rolling direction (see also equations 197 (19), (21) and (25)). The hardening of the material, which is considered to be isotropic and governed by the 198 e ective plastic strain  "p, is assumed to be of the form: 199 Y=0 "0+ "pn(10) where0,"0, andnare material parameters, see Table 1. 200 201 The e ective plastic strain is de ned as: 202 "p=Zt 0_"pd (11) 10 with _"pbeing the work-conjugate of the speci c e ective stress associated to the yield criterion (see equations 203 (19), (21) and (25)). Recall that a superposed dot denotes di erentiation with respect to time. 204 205 Moreover, assuming an associated plastic ow rule, the plastic part of the rate of deformation tensor is: 206 dp=_@ @(12) where _is the rate of plastic multiplier. 207 208 Since the e ective stress is a rst order homogeneous function in stress (see equations (19), (21) and (25)), 209 we have from the Euler's homogeneous function theorem that: 210 =:@ @(13) Therefore, the work conjugacy relation: 211 :dp= _"p (14) leads to the identity: 212 _"p=_ (15) The formulation of the constitutive model is completed with the Kuhn-Tucker loading-unloading conditions: 213 _>0; f 60; _f= 0 (16) and the consistency condition during plastic loading: 214 11 Symbol Property and units Value 0Initial density (kg =m3) 2700 G Elastic shear modulus (GPa), Eq. (4) 26 :92 K Bulk modulus (GPa), Eq. (4) 58 :33 0Material parameter (MPa), Eq. (10) 646 "0Material parameter, Eq. (10) 0 :025 n Material parameter, Eq. (10) 0 :227 Table 1: Initial density, elastic constants and strain hardening parameters for AA 2090-T3. Data after Yoon et al. (2006) and Cvitani c et al. (2008). _f= 0 (17) 215 Table 1 shows the initial density, the elastic constants, and the hardening parameters used in the nite 216 element simulations shown in Section 7. The parameters values correspond to AA 2090-T3 and they are taken 217 from Yoon et al. (2006) and Cvitani c et al. (2008). 218 2.2. Yld2004-18p yield criterion 219 The anisotropy is introduced using 2 linear transformations of the stress deviator, equations (6a) and 220 (6b), which provide 18 anisotropy coecients, equations (7a) and (7b), whose values must be determined from 221 experiments (recall that for Yld2004-18p yield criterion c0 ii=c00 ii= 0, withi= 1;2;3, and therefore the matrices 222 C0andC00only have 9 non-zero coecients). Note that, if C0=C00, i.e. if the formulation accounts for only 223 one linear transformation, the Yld2004-18p criterion reduces to the Yld91 model (Barlat et al., 1991), provided 224 thatc0 12=c0 21,c0 13=c0 31andc0 32=c0 23. 225 226 The yield function is: 227 Yld04 =3X i=13X j=1 ~S0 i~S00 j a (18) and the associated e ective stress is: 228 12 Anisotropy parameters c0 11c0 12c0 13c0 21c0 22c0 23c0 31c0 32c0 33c0 44c0 55c0 66 00:06989 0:93640 0:07914 01:00306 0:524741 1:36318 01:02377 1:06906 0:95432 c00 11c00 12c00 13c00 21c00 22c00 23c00 31c00 32c00 33c00 44c00 55c00 66 0 0:98117 0:47674 0:57531 00:86682 1:145010:07929 01:05166 1:14700 1:40462 Degree of homogeneity of the yield function a 8 Table 2: YLD2004-18p yield criterion. Parameters values for AA 2090-T3. Data after Barlat et al. (2005). = Yld04 1=a (19) with= 4, so that  reduces to the yield stress in the rolling direction (see equation (9)). 229 230 The parameters values corresponding to AA 2090-T3 are given in Table 2. 231 2.3. CPB06ex2 yield criterion 232 The model accounts for both the anisotropy and tension-compression asymmetry of the material, which 233 are introduced, as in the case of the Yld2004-18p model (Barlat et al., 2005), using 2 linear transformations 234 of the stress deviator (equations (6a) and (6b)). The CPB06ex2 yield criterion contains 21 parameters: 18 to 235 describe the anisotropy ( c0 ijandc00 ijin equations (7a) and (7b)), 2 for the tension-compression asymmetry ( k0236 andk00in equation (20)), and 1 to determine the degree of homogeneity of the yield function. Note that when 237 C0=C00andk0=k00, the CPB06ex2 model reduces to the yield criterion developed by Cazacu et al. (2006). 238 239 The yield function is: 240 CPB06 =3X i=1(j~S0 ijk0~S0 i)a+3X i=1(j~S00 ijk00~S00 i)a(20) and the associated e ective stress is: 241 13 Anisotropy parameters c0 11c0 12=c0 21c0 13=c0 31c0 22c0 23=c0 32c0 33 0:453 0:841 1:248 1:058 2:284 3:201 c0 44c0 55c0 66c00 11c00 12=c00 21c00 13=c00 31 111:0260:453 0:7051:148 c00 22c00 23=c00 32c00 33c00 44c00 55c00 66 0:139 0:5190:878111:978 Degree of homogeneity of the yield function a 12 Tension-compression asymmetry parameters k0k00 0:054 0:027 Table 3: CPB06ex2 yield criterion. Parameters values for AA 2090-T3. Data after Plunkett et al. (2008). = CPB06 1=a (21) with 242 =3X i=1(j'0 ijk0'0 i)a+3X i=1(j'00 ijk00'00 i)a(22) where 243 '0 1=(2c0 11c0 12c0 13) 3; '0 2=(2c0 12c0 22c0 23) 3; '0 3=(2c0 13c0 23c0 33) 3 '00 1=(2c00 11c00 12c00 13) 3; '00 2=(2c00 12c00 22c00 23) 3; '00 3=(2c00 13c00 23c00 33) 3(23) Notice that the function is such that  reduces to the tensile yield stress in the rolling direction (see 244 equation (9)). 245 246 The parameters values corresponding to AA 2090-T3 are given in Table 3. 247 14 2.4. Yld2011-27p yield criterion 248 The anisotropy is introduced using 3 linear transformations of the stress deviator, equations (6a)-(6c), 249 leading to 27 anisotropy coecients, equations (7a)-(7c), thereby providing extended exibility to capture the 250 mechanical response of materials with complex anisotropic behaviour. Recall that for the Yld2011-27p model 251 c0 ii=c00 ii=c000 ii= 0, withi= 1;2;3. 252 253 The yield function is: 254 Yld11 =3X i=13X j=1 ~S0 i+~S00 j a +3X k=1 ~S000 k a (24) and the associated e ective stress is: 255 = Yld11 1=a (25) with 256 =4 3a + 52 3a + 61 3a (26) Similarly to the Yld2004-18p and CPB06ex2 criteria, the function is chosen so that  reduces to the yield 257 stress in the rolling direction (see equation (9)). 258 259 The parameters values corresponding to AA 2090-T3 are given in Table 4. 260 3. Integration of the constitutive equations 261 The numerical integration is a strain-driven process where the increment of the total strain tensor 4" 262 is given at a quadrature point, and the stress tensor should be updated for the corresponding time incre- 263 ment. The incremental integration is local in space and occurs at each quadrature point inside the nite 264 elements. The three constitutive models presented in Section 2 have been implemented in ABAQUS/Explicit 265 and ABAQUS/Standard through user subroutines VUMAT and UMAT, respectively, using the classical return 266 mapping scheme (Sim o and Taylor, 1985; Sim o and Hughes, 1998). The speci c algorithms used in the VUMAT 267 15 Anisotropy parameters c0 11c0 12c0 13c0 21c0 22c0 23c0 31c0 32c0 33c0 44c0 55c0 66 0 0:12667 1:29525 0:28407 0 1:42858 0:604524 1:36079 0 1 1 0:8364 c00 11c00 12c00 13c00 21c00 22c00 23c00 31c00 32c00 33c00 44c00 55c00 66 0 1:44953 0:99770 1:34779 0 1:00332 1:786760:08678 0 1 1 1 c000 11c000 12c000 13c000 21c000 22c000 23c000 31c000 32c000 33c000 44c000 55c000 66 00:61319 0:392401:31670 02:31010 0:72917 0:75833 0 1 12:78319 Degree of homogeneity of the yield function a 12 Table 4: Yld2011-27p yield criterion. Parameters values for AA 2090-T3. Data after Aretz and Barlat (2013). and UMAT subroutines are shown in Sections 3.1 and 3.2, respectively. The goal is to calculate the increment 268 of the e ective plastic strain which is required to update the stress state. In what follows, superscripts nand 269 n+ 1 denote the previous and current time steps, respectively. The return is performed at time n+ 1. 270 3.1. Convex cutting-plane algorithm 271 We rewrite the yield condition as: 272 fn+1= n+1n+1 Y= 0 (27) The preceding relation is a non-linear algebraic equation in 4"pwhich can be solved using an iterative 273 Newton-Raphson procedure. For that purpose, the previous expression is linearized as follows: 274 fn+1 (k+1)fn+1 (k)+@f @n+1@n+1 @n+1@n+1 @4"p (k)"p(k)+ @f @n+1 Y@n+1 Y @"n+1p@"n+1 p @4"p! (k)"p(k)= 0 (28) wherekstands for the iterative index and: 275 @f @n+1= 1;@n+1 @4"p=L:@n+1 @n+1 @f @n+1 Y=1;@"n+1 p @4"p= 1 (29) 16 From equation (28), "p(k)is calculated as: 276 "p(k)=f @ @:L:@ @+@Y @"p n+1 (k)(30) wheren+1 (k)denotes that "p(k)is computed at the current time step n+ 1 for a given iteration k. 277 278 The increment of the e ective plastic strain and the updated stress state are calculated after every iteration 279 as follows: 280 4"p(k+1)=4"p(k)+"p(k) (31) n+1 (k+1)=n+1 (k)L:@n+1 @n+1 (k)"p(k) (32) The iterative procedure is performed until fn+1 (k+1)is lower than a user-de ned tolerance. Then, the cor- 281 responding stress state is accepted as the current stress state. Note that the convex cutting-plane algorithm 282 requires to obtain the rst order yield function gradients (i.e. the rst derivative of the e ective stress with 283 respect to the components of the stress tensor) to determine the increment of the e ective plastic strain. 284 3.2. Implicit backward-Euler algorithm 285 We rewrite the yield condition and the ow rule as: 286 fn+1= n+1n+1 Y= 0 (27 revisited) gn+1=4"p4"p@n+1 @n+1= 0 (33) These expressions are non-linear algebraic equations in 4"pwhich can be solved using an iterative Newton- 287 17 Raphson procedure. For that task, the previous expressions are linearized as follows: 288 fn+1 (k+1)fn+1 (k)+@f @n+1@n+1 @n+1@n+1 @4"p (k)"p(k)+ @f @n+1 Y@n+1 Y @"n+1p@"n+1 p @4"p! (k)"p(k)= 0 (28 revisited) gn+1 (k+1)gn+1 (k)+ @g @4"p@4"p @"n+1p@"n+1 p @4"p! (k)"p(k)+@g @n+1@n+1 @4"p (k)"p(k)+@g @4"p (k)"p(k)= 0 (34) where 289 @f @n+1= 1 @f @n+1 Y=1;@"n+1 p @4"p= 1(29 revisited) and 290 @g @4"p=I;@4"p @"n+1p=I;@"n+1 p @4"p=L1:@n+1 @4"p @g @n+1=4"p@2n+1 @n+1@n+1 @g @4"p=@n+1 @n+1(35) where Iis the unit fourth-order tensor. 291 292 From equation (34), we obtain: 293 @n+1 @4"p (k)"p(k)=En+1 (k) gn+1 (k)+@n+1 @n+1 (k)"p(k)! (36) where 294 En+1 (k)= L1+4"p(k)@2n+1 @n+1@n+1 (k)!1 (37) 295 18 Inserting equation (36) into equation (28), we obtain: 296 "p(k)=f@ @:E:g @ @:E:@ @+@Y @"p n+1 (k)(38) 297 The increment of the e ective plastic strain and the updated stress state are calculated after every iteration 298 as follows: 299 4"p(k+1)=4"p(k)+"p(k) (31 revisited) n+1 (k+1)=n+1 (k)En+1 (k) gn+1 (k)+@n+1 @n+1 (k)"p(k)! (39) The iterative procedure is performed until fn+1 (k+1)andgn+1 (k+1)are lower than a user-de ned tolerance. Then, 300 the corresponding stress state is accepted as the current stress state. Note that the implicit backward-Euler 301 algorithm requires to obtain the rst and second order yield function gradients (i.e. the rst and second 302 derivatives of the e ective stress with respect to the components of the stress tensor) to determine the increment 303 of the e ective plastic strain. 304 4. Algorithmic tangent modulus 305 In this section, we derive the general expression of the algorithmic tangent modulus, a critical factor in the 306 UMAT implementation of the constitutive models to preserve the quadratic rate of asymptotic convergence of 307 the iterative Newton-Raphson procedure (see equations (27) and (33)). See, for instance, Yoon et al. (1999) 308 and Yoon et al. (2004) for additional details. 309 310 The stress at time n+ 1 is: 311 n+1=n+L:4"4"pL:@n+1 @n+1(40) 19 Di erentiation of previous expression leads to: 312 =L "@n+1 @n+1"p4"p@2n+1 @n+1@n+1 (41) which can be recast into: 313 =L "@n+1 @n+1"p (42) where: 314 L= L1+4"p@2n+1 @n+1@n+11 (43) Di erentiation of the consistency condition, equation (27), leads to: 315 f=@f @n+1@n+1 @n+1 + @f @n+1 Y@n+1 Y @"n+1p@"n+1 p @4"p! "p= 0 (44) Substituting equation (42) into equation (44) gives the following relation: 316 "p=@n+1 @n+1:L:" @n+1 @n+1:L:@n+1 @n+1+@n+1 Y @"n+1p(45) Inserting equation (45) into equation (42) yields: 317 =Lep:" (46) where: 318 Lep=0 BBB@LL:@n+1 @n+1 L:@n+1 @n+1 @n+1 @n+1:L:@n+1 @n+1+@n+1 Y @"n+1p1 CCCA(47) is the algorithmic tangent modulus. Note that we need to obtain the rst and second order yield function 319 gradients to calculate the algorithmic tangent modulus. 320 321 The idea of giving details on the formulation of the constitutive models in Section 2, and on the integration 322 20 algorithms in Sections 3 and 4, stems from our goal of providing a self-contained manuscript, that it is easy to 323 follow and fully reproducible. 324 5. Numerical approximations for the stress derivatives 325 This section presents a numerical approximation of the rst and second derivatives of the e ective stress 326 with respect to the components of the stress tensor. This is a simple alternative to the analytical derivatives, 327 which are given in Appendix A, that reduces the computational time without impairing the accuracy of 328 the nite element calculations (see Section 7). As mentioned before, these derivatives are required for the 329 integration of the constitutive equations, namely, for the calculation of the e ective plastic strain increment 330 in the convex cutting-plane and the implicit backward-Euler algorithms, equations (30) and (38), and for the 331 calculation of the algorithmic tangent modulus, see equation (47). 332 5.1. First derivative approximation 333 The interpolation polynomial of Lagrange, for a set of k+1 data points ( x0;h(x0)),...,(xj;h(xj)),...,(xk;h(xk)), 334 where no two xjare the same, is: 335 h(x)kX j=0h(xj)Lj(x) (48) where 336 Lj(x) =kY m=0;m6=jxxm xjxm=xx0 xjx0:::xxj1 xjxj1xxj+1 xjxj+1:::xxk xjxk(49) For a set of 3 equally-spaced data points ( x0;h(x0));(x1;h(x1));(x2;h(x2)) which ful ll the conditions 337 x0=x1aandx2=x1+a, the derivative of equation (48) leads to: 338 @h(x) @x x1h(x2)h(x0) 2a(50) Particularization of previous expression for the derivative of the e ective stress with respect to the stress 339 components leads to: 340 @ @ij(ij+ ij)(ijij) 2ij(51) 21 Notice that previous expression is equivalent to the central di erence scheme used by Choi and Yoon (2019), 341 see Appendix B, and di erent from the forward di erence scheme reported in Appendix A.2 of Banabic et al. 342 (2005), Appendix C of Barlat et al. (2005) and Aretz (2007). 343 5.2. Second derivatives approximation 344 The rst three terms of the Taylor series expansion for approximating h(x1+a) andh(x1a) aboutx1 345 are: 346 h(x)jx1+ah(x)jx1+@h(x) @x x1a+1 2@2h(x) @x@x x1a2(52) and 347 h(x)jx1ah(x)jx1@h(x) @x x1a+1 2@2h(x) @x@x x1a2(53) Adding equations (52) and (53) we obtain: 348 @2h(x) @x@x x1h(x1a)2h(x1) +h(x1+a) a2(54) Particularization of previous expression for the second derivative of  with respect to the stress components 349 whenij=klleads to: 350 @2 @ij@ij(ijij)2+ (ij+ ij) 2 ij(55) Notice that previous expression is the same used by Choi and Yoon (2019), see Appendix B. 351 352 Moreover, the rst three terms of the Taylor series expansion for approximating h(x1+a;y1+a) and 353 h(x1a;y1a) about (x1;y1) are: 354 h(x;y)j(x1+a;y1+a)h(x;y)j(x1;y1)+@h(x;y) @x (x1;y1)a+@h(x;y) @y (x1;y1)a +1 2@h2(x;y) @x@x (x1;y1)a2+@h2(x;y) @x@y (x1;y1)a2+1 2@h2(x;y) @y@y (x1;y1)a2(56) 22 and 355 h(x;y)j(x1a;y1a)h(x;y)j(x1;y1)@h(x;y) @x (x1;y1)a@h(x;y) @y (x1;y1)a +1 2@h2(x;y) @x@x (x1;y1)a2+@h2(x;y) @x@y (x1;y1)a2+1 2@h2(x;y) @y@y (x1;y1)a2(57) Adding equations (56) and (57) leads to: 356 @h2(x;y) @x@y (x1;y1)h(x1a;y1a)2h(x1;y1) +h(x1+a;y1+h) 2a21 2@h2(x;y) @x@x (x1;y1)1 2@h2(x;y) @y@y (x1;y1)(58) Particularizing previous expression for the second derivative of  with respect to the stress components 357 whenij6=klleads to: 358 @2 @ij@kl(ijij;klkl)2+ (ij+ ij;kl+ kl) 2ijkl1 2@2 @ij@ij1 2@2 @kl@kl(59) where ij= kl. Notice that previous expression is di erent from the formulation used by Choi and 359 Yoon (2019), see Appendix B, and thus an original contribution of this paper. All the calculations shown in 360 Section 7 are performed with  ij= kl= 0:01 MPa. Nevertheless, we also carried out simulations with 361 ij= kl= 0:001 and 0:0001 MPa and the di erences in the numerical results (stress and strain elds in 362 the specimen) were negligible. Moreover, we have checked that the relative error of the numerical derivatives 363 with respect to their analytical counterparts, for the three stress steps  ij= kl= 0:1, 0:01 and 0:001 MPa, 364 is less than 108for all the calculations performed in this work, thus having negligible in uence in the nite 365 element results (see Section 7.1). 366 367 Notice that the simple form of equations (51), (55) and (59), in comparison with the analytical derivatives 368 shown in Appendix A, greatly facilitates the implementation of the constitutive models in nite element 369 codes. It was also pointed out by Barlat et al. (2005) that the singularities of the analytical derivatives shown 370 in Appendix A.4 can be avoided by using numerical approximation of the yield function gradients. As it will be 371 shown in Section 7, using the numerical approximation of the second derivative leads to a signi cant decrease 372 of the computational time in simulations performed with UMAT subroutines in ABAQUS/Standard implicit 373 23 calculations. In addition, the computational cost of the numerical approximations proposed in this paper 374 is smaller than that of the numerical di erentiation scheme used by Choi and Yoon (2019), see Section 7.2. 375 The comparison performed in Section 7 between the computational eciency of the numerical di erentiation 376 schemed used in this paper and the approach of Choi and Yoon (2019) is a novel aspect of this manuscript. In 377 this regard, note that Choi and Yoon (2019) only performed small scale simulations: single element calculations 378 to evaluate the accuracy of the numerical derivatives, and deep drawing simulations where the blank is modeled 379 with a single layer of shell elements to check the number of ears predicted by the various constitutive models 380 investigated. The larger scale calculations presented in Section 7, with nite element models including tens of 381 thousands of elements, are an original contribution of this paper. 382 6. Finite element models 383 This section presents the nite element models developed in ABAQUS (2016) to carry out cylindrical cup 384 deep drawing test calculations (Section 6.1) and simulations of dynamic necking localization under plane strain 385 tension (Section 6.2) with AA 2090-T3 specimens modeled with the three yield criteria presented in Section 2. 386 6.1. Cup deep drawing 387 Fig. 1 shows a schematic representation of the cylindrical cup deep drawing process. The speci c tools and 388 blank dimensions are given in Table 5. 389 Figure 1: Tools and blank geometries for cylindrical cup deep drawing. This gure is adapted from Yoon et al. (1999, 2004). Taking advantage of the problem symmetry and the orthotropy of the material, only a quarter section of 390 the cup is analyzed. The loading condition is a prescribed constant velocity imposed to the punch of 0 :4 m=s. 391 24 Dp(mm)Dd(mm)Db(mm)rp(mm)rd(mm)e(mm)t(mm) 97.46 101.48 158.76 12.7 12.7 2.7 1.6 Table 5: Tools and blank dimensions for the cup deep drawing test. Data after Yoon et al. (2006). Similarly to Aretz and Barlat (2013), a blank-holder force of 5 :5 kN for the quarter model (corresponding to 392 approximately one percent of the initial yield stress of the material in the rolling direction) has been applied 393 during the simulations. The friction is modeled with Coulomb's law, with a friction coecient equal to 0 :1 for 394 all contact surfaces (Yoon et al., 2006). 395 396 Punch, holder and die are modeled as rigid solids. Moreover, the blank has been meshed with tri-linear 397 solid elements with 8-nodes (C3D8 in ABAQUS notation). We have used 5 di erent meshes (mesh 1,...,mesh 398 5) in which the number of elements trough the thickness increases from 1 to 5. The number of elements in 399 each mesh is: 5915, 11830, 17745, 23660 and 29575. The goal is to illustrate the computational time saving 400 obtained with the numerical approximation of the stress derivatives as the number of elements of the model 401 increases. The mesh design of the blank with 5 elements through the thickness is shown in Fig. 2. 402 Figure 2: Blank mesh used in the nite element simulations. Five elements through the thickness of the sheet. Total number of elements: 29575. 403 To the authors' knowledge, there are just few relevant papers using several solid elements through the 404 thickness of the blank to simulate deep drawing of circular cups (Rousselier et al., 2009). Deep drawing 405 calculations are usually performed with a single solid element through the thickness of the specimen, notably 406 25 when advanced constitutive models are used (see Yoon et al. (2006); Vladimirov et al. (2010); Shutov et al. 407 (2015)), often with the aim of limiting the computational time. Shell type elements have also been used in 408 several works (Yoon et al., 2004, 2010; Aretz and Barlat, 2013; Choi and Yoon, 2019). As mentioned before, the 409 cylindrical cup deep drawing test calculations performed in this paper with nite element models that include 410 tens of thousands of elements, are an original contribution of this work. 411 6.2. Plane strain tension 412 This section presents the 2D nite element model used to simulate dynamic necking localization under 413 plane strain tension. Material points are referred to using a Cartesian coordinate system with positions in the 414 reference con guration denoted as ( X;Y ). The specimen is under plane strain constraint in the out-of-plane 415 direction. The origin of coordinates is located in the center of mass of the specimen, see Fig. 3. We use the 416 same nite element model as N'souglo et al. (2020) who, following the earlier work of Xue et al. (2008), modeled 417 a plate subjected to plane strain tension, and with geometric periodic perturbations, as an array of unit cells 418 with sinusoidal spatial imperfections given by the following expression: 419 h=h0 1 + cos2X L0 (60) 420 wherehis the perturbed thickness of the unit cell and is the amplitude of the imperfection. Moreover, L0421 andh0are the initial length and the unperturbed thickness of the unit cell, respectively (see Fig. 3). Due to 422 the symmetry of the problem, only a quarter of the unit cell is analyzed, see Fig. 3. For all the calculations 423 presented in Section 7.2, the rolling direction is aligned with the Xaxis. Finite element simulations will be 424 conducted for several normalized cell lengths varying within the range 0 :5<L0 h06, with  =2 h0= 0:2% and 425 h0= 2 mm. Hereinafter, the normalized cell length will be also referred to as cell size. Note that the spatial 426 imperfection is imposed on lines X=constant which is the orientation naturally selected by the material to 427 trigger a neck under plane strain tension when the orthotropy and the loading are codirectional. The initial 428 and boundary conditions used in the nite element model are: 429 Vx(X;Y; 0) = _"0 xxX Vy(X;Y; 0) =_"0 xxY(61) 26 and 430 Ux(0;Y;t) = 0; U y(X;0;t) = 0 Vx L0=2;Y;t = _"0 xxL0=2(62) where _"0 xxis the initial strain rate imposed on the loading direction. 431 Figure 3: Finite element model. Mesh and boundary conditions. A large imperfection amplitude has been shown in the gure for better illustration of the geometric perturbation. This gure is taken from N'souglo et al. (2020). The specimen is under plane strain constraint in the out-of-plane direction. The nite element model is meshed with four-node bilinear plane strain elements with reduced integration 432 and hourglass control (CPE4R in ABAQUS notation), with initial dimensions 2525m2. 433 7. Sample results 434 In this section, we present sample nite elements results for deep drawing (Section 7.1) and dynamic necking 435 localization under plane strain tension (Section 7.2). We have performed implicit and explicit nite element 436 simulations using UMAT and VUMAT subroutines, respectively, with the analytical derivatives of Appendix 437 A, the numerical approximations of Section 5, and the numerical di erentiation scheme developed by Choi and 438 Yoon (2019), which is shown in Appendix B. Note that we provide access to the UMAT and VUMAT codes 439 developed in this paper. Recall that the convex cutting-plane algorithm used to code the VUMAT involves only 440 the rst derivative of the yield function gradients, which takes the same form in the numerical di erentiation 441 scheme developed in Section 5 and in the formulation proposed by Choi and Yoon (2019), see Appendix B. All 442 27 the calculations have been performed with a personal computer using a single CPU with the following features: 443 Intel(R) Core(TM) i78550U CPU @ 1 :80 GHz 1:99 GHZ. 444 7.1. Cup deep drawing 445 This section presents the results corresponding to the cup deep drawing nite element calculations. The 446 simulations have been performed using a dynamic implicit analysis in ABAQUS/Standard and UMAT subrou- 447 tines for the yield functions implementation. The calculations are terminated when the displacement of the 448 punch reaches 60 mm, so that we make sure that the deep drawing operation has been completed (additional 449 punch displacement does not lead to further specimen deformation, see Figs. 6 and 7). Table 6 shows the CPU 450 times corresponding to the simulations performed with the three yield criteria, using ve di erent meshes for 451 the blank in which the number of elements trough the thickness varies from 1 to 5 (see Section 6.1). The CPU 452 times obtained with the analytical derivatives of Appendix A (denoted with -A) and the numerical approxima- 453 tions of Section 5 (denoted with -N) are compared. Using the numerical approximations leads to a signi cant 454 decrease of the computational time for the three yield criteria, and this decrease is greater as the number of ele- 455 ments in the mesh increases. This is an original result of this paper. Namely, for the Yld2004-18p criterion, the 456 CPU times obtained using the numerical approximations of the yield function gradients are 47% and32% 457 of those obtained using the analytical derivatives for the meshes with 1 and 5 elements trough the thickness, 458 respectively. For the CPB06ex2 model these percentages are 57% and46% (i.e. smaller computational 459 time decrease than in the case of Yld2004-18p criterion), and for the Yld2011-27p model they are 47% and 460 35%, respectively. It becomes apparent that the decrease of the computational time using the numerical 461 approximations enables to perform larger scale ABAQUS/Standard calculations. For instance, the CPU time 462 for Yld2004-18p-N with 5 elements through the thickness is only 1:35 times the CPU time corresponding to 463 Yld2004-18p-A with 1 element trough the thickness. 464 28CPU time (units in kiloseconds) Yld2004-18p-A Yld2004-18p-N CPB06ex2-A CPB06ex2-N Yld2011-27p-A Yld2011-27p-N 1 Element 128:752 60 :520 107 :991 61 :940 132 :977 62 :957 2 Elements 211:607 80 :443 184 :322 87 :444 229 :964 90 :574 3 Elements 319:536 109 :008 265 :738 118 :439 327 :210 120 :531 4 Elements 411:566 132 :924 355 :614 159 :813 435 :562 154 :776 5 Elements 544:324 174 :803 448 :508 205 :446 567 :211 199 :222 Table 6: CPU time for cup deep drawing nite element simulations performed with the three yield criteria presented in Section 2. The results obtained with the analytical derivatives of Appendix A (denoted with -A) and the numerical approximations of Section 5 (denoted with -N) are compared. The simulations have been performed using a dynamic implicit analysis in ABAQUS/Standard and UMAT subroutines for the yield functions implementation. 29 465 A key point is that, despite the signi cant reduction in computational time, the results are practically the 466 same using the analytical derivatives and the numerical approximations of the yield functions gradients. Fig. 4 467 shows nite element results for the cup height uas a function of the angle with respect to the rolling direction 468 . The cup height is measured when the displacement of the punch reaches 60 mm (when the simulation 469 is terminated). The calculations are performed with 1 element trough the thickness. The nite elements 470 results correspond to: (a) Yld2004-18p-A and Yld2004-18p-N, (b) CPB06ex2-A and CPB06ex2-N and (c) 471 Yld2011-27p-A and Yld2011-27p-N (same notation used in Table 6). The ucurves obtained with the 472 analytical derivatives and the numerical approximations are virtually on top of each other (we have checked 473 that this agreement occurs for the 5 meshes investigated). Moreover, note that we also compare the nite 474 element results with the experimental data reported by Yoon et al. (2006) for AA 2090-T3. The qualitative 475 agreement between the predictions of the three yield criteria and the experiments is satisfactory (note that 476 Yoon et al. (2006) already compared these experiments with the predictions of the Yld2004-18p-A criterion), 477 with a maximum of the cup height for an intermediate value of , and minima near the rolling and transverse 478 directions, with the cup height being greater for = 0than for= 90. However, the three models show 479 quantitative di erences with the experiments, especially in the prediction of the cup height for = 0, and in 480 the maximum value of ufor50. According to Yoon et al. (2006), the discrepancies between the numerical 481 predictions and the experiments could be partially attributed to the fact that the center of the blank was not 482 aligned properly with the centers of the die and the punch during the drawing test. 483 484 Fig. 5 shows the comparison of the cup height results obtained with the numerical approximation of the 485 yield function gradients of Section 5 for the meshes with 1 and 5 elements trough the thickness of the sheet, for 486 the three yield criteria: (a) Yld2004-18p-N, (b) CPB06ex2-N and (c) Yld2011-27p-N. Note that the results with 487 1 element trough the thickness were already presented in Fig. 4. The increase in the number of elements shifts 488 upwards the ucurves, such that the cup height, for any angle , is1 mm greater when 5 elements trough 489 the thickness are used, leading to slightly better quantitative agreement with the experiments of Yoon et al. 490 (2006). Moreover, increasing the number of elements generally helps to describe plastic localization phenomena, 491 reducing the pathological mesh sensitivity obtained with rate-independent materials in problems where inertia 492 forces are not relevant (Needleman, 1988). The contours of e ective plastic strain presented in Fig. 6 for the 493 Yld2004-18p criterion show a comparison between the results obtained with the analytical derivatives and the 494 numerical approximation of the yield function gradients proposed in Section 5, for simulations performed with 495 30 40424446485052 0 15 30 45 60 75 90B D FCup height, u (mm) Angle with respect to the rolling direction, θ (o)FEM: Yld2004-18p-A FEM: Yld2004-18p-N Experiments: Yoon et al. ( 2006 )1 Element (a) 40424446485052 0 15 30 45 60 75 90B D F Angle with respect to the rolling direction, θ (o)Cup height, u (mm)FEM: CPB06ex2-A FEM: CPB06ex2-N Experiments: Yoon et al. ( 2006 )1 Element (b) 40424446485052 0 15 30 45 60 75 90B D F Angle with respect to the rolling direction, θ (o)Cup height, u (mm)1 Element FEM: Yld2011-27p-A FEM: Yld2011-27p-N Experiments: Yoon et al. ( 2006 ) (c) Figure 4: Comparison between nite element simulations (FEM) and experimental results. Cup height uas a function of the angle with respect to the rolling direction . The experimental data are taken from Yoon et al. (2006). The nite elements results are obtained with the analytical derivatives of Appendix A (denoted with -A) and the numerical approximations of Section 5 (denoted with -N). The calculations are performed with 1 element trough the thickness of the sheet. (a) Yld2004-18p-A and Yld2004-18p-N. (b) CPB06ex2-A and CPB06ex2-N. (c) Yld2011-27p-A and Yld2011-27p-N. The simulations have been performed using a dynamic implicit analysis in ABAQUS/Standard and UMAT subroutines for the yield functions implementation. 31 1 and 5 elements trough the thickness. The contours are taken when the displacement of the punch reaches 496 60 mm (when the simulations are terminated). All the isocontours have the same color coding, so that e ective 497 plastic strains ranging from 0 to 0 :7 correlate with a colour scale that goes from blue to red. E ective plastic 498 strains greater than 0 :7 remain red, so that the red zones in the cups correspond to areas with highly localized 499 plastic deformation and important thinning of the thickness of the cup. These red zones, which are located 500 between the 8 ears that develop in the cups (which are numbered in Fig. 6(a)) have greater extension for 501 the calculations performed with 5 elements trough the thickness of the specimen, illustrating the in uence of 502 the mesh in simulating plastic localization processes (when rate-independent materials are used). As stated 503 by Yoon et al. (2006), the Yld2004-18p criterion, unlike many other phenomenological yield functions, shows 504 the capability of predicting more than four ears in cup drawing. Moreover, note that the contours of plastic 505 strain for the calculations performed with the analytical derivatives and the numerical approximation of the 506 yield function gradients proposed in Section 5 are virtually identical (compare Figs. 6(a) and 6(c) with Figs. 507 6(b) and 6(d), respectively). In future works, it is recommended to include strain rate and temperature e ects 508 in the evolution of the yield stress of the material, in order to develop a more comprehensive analysis on the 509 material aspects that control plastic localization in deep drawing operations. 510 511 512 Fig. 7 shows contours of e ective plastic strain  "pcorresponding to calculations performed with 5 elements 513 trough the thickness of the sheet for CPB06ex2 and Yld2011-27p criteria. The color coding is the same used in 514 Fig. 6. The results obtained with the analytical derivatives of Appendix A and the numerical approximation 515 of the yield function gradients of Section 5 are compared: (a) CPB06ex2-A, (b) CPB06ex2-N, (c) Yld2011- 516 27p-A and (d) Yld2011-27p-N. As for the Yld2004-18p criterion, the procedure used to calculate the stress 517 derivatives has no in uence on the plastic strain contours (compare Figs. 7(a) and 7(c) with Figs. 7(b) and 518 7(d), respectively). Eight ears are formed in the four calculations shown in Fig. 7, being the height of numbers 519 2, 4, 6 and 8 greater than the height of the other four. On the other other hand, the yield criterion plays a 520 role in the localization zones between the ears of the cup, which have slightly greater extension in the case 521 of CPB06ex2 model (compare Figs. 7(a) and 7(b) with Figs. 7(c) and 7(d), respectively). Nevertheless, the 522 distribution of plastic strains in the specimen is similar for both yield criteria. 523 32 40424446485052 0 15 30 45 60 75 90B D F Angle with respect to the rolling direction, θ (o)Cup height, u (mm)FEM: 1 Element FEM: 5 Elements Experiments: Yoon et al. ( 2006 )Yld2004-18p-N (a) 40424446485052 0 15 30 45 60 75 90B D F Angle with respect to the rolling direction, θ (o)Cup height, u (mm)FEM: 1 Element FEM: 5 Elements Experiments: Yoon et al. ( 2006 )CPB06ex2-N (b) 40424446485052 0 15 30 45 60 75 90B D F Angle with respect to the rolling direction, θ (o)Cup height, u (mm)FEM: 1 Element FEM: 5 Elements Experiments: Yoon et al. ( 2006 )Yld2011-27p-N (c) Figure 5: Comparison between nite element simulations (FEM) and experimental results. Cup height uas a function of the angle with respect to the rolling direction . The experimental data are taken from Yoon et al. (2006). The nite elements results are obtained with the numerical approximations of Section 5 (denoted with -N). The calculations are performed with 1 element and 5 elements through the thickness of the sheet. (a) Yld2004-18p-N. (b) CPB06ex2-N. (c) Yld2011-27p-N. The simulations have been performed using a dynamic implicit analysis in ABAQUS/Standard and UMAT subroutines for the yield functions implementation. 33 Figure 6: Finite element simulations. Cup deep drawing calculations. Contours of e ective plastic strain  "pcorresponding to calculations performed with Yld2004-18p yield criterion. Comparison between results obtained with the analytical derivatives of Appendix A (denoted with -A) and the numerical approximations of Section 5 (denoted with -N). Calculations performed with 1 and 5 elements trough the thickness. (a) Yld2004-18p-A and 1 element trough the thickness. (b) Yld2004-18p-N and 1 element trough the thickness. (c) Yld2004-18p-A and 5 elements trough the thickness. (d) Yld2004-18p-N and 5 elements trough the thickness. The simulations have been performed using a dynamic implicit analysis in ABAQUS/Standard and UMAT subroutines for the yield function implementation. 34 Figure 7: Finite element simulations. Cup deep drawing calculations. Contours of e ective plastic strain  "pcorresponding to calculations performed with 5 elements trough the thickness. Comparison between results obtained with the analytical derivatives of Appendix A (denoted with -A) and the numerical approximations of Section 5 (denoted with -N). (a) CPB06ex2-A. (b) CPB06ex2-N. (c) Yld2011-27p-A. (d) Yld2011-27p-N. The simulations have been performed using a dynamic implicit analysis in ABAQUS/Standard and UMAT subroutines for the yield functions implementation. 35 7.2. Plane strain tension 524 This section presents the results corresponding to the dynamic plane strain tension calculations. We have 525 performed simulations with both UMAT and VUMAT subroutines, using the numerical approximation of the 526 yield function gradients of Section 5 and the numerical di erentiation scheme developed by Choi and Yoon 527 (2019), which is shown in Appendix B. While we do not show the results for the sake of brevity, we have 528 checked that the explicit calculations with the analytical expressions of the yield function gradients are not 529 computationally more expensive, which brings out that the numerical approximation of the yield function 530 gradients is computationally advantageous only if the second order yield function gradients are involved in 531 the integration scheme (i.e. for the UMAT, as shown in Section 7.1). The calculations are carried out for 532 two imposed initial strain rates _ "0 xx= 100 s1and 10000 s1(see Section 6.2). This range of strain rates is 533 usually attained in high energy rate forming processes, such as electromagnetic and electrohydraulic forming, 534 in which the formation of dynamic necks imposes limitations to the material formability (Balanethiram and 535 Daehn, 1992, 1994). 536 537 Fig. 8 shows the evolution of the necking strain  "neckwith the cell size L0=h0, for dynamic implicit 538 and explicit simulations performed for _ "0 xx= 100 s1with the numerical approximation of the yield function 539 gradients and the three yield criteria. The implicit calculations are performed with the numerical approximation 540 of the stress derivatives of Section 5 and the numerical di erentiation scheme developed by Choi and Yoon 541 (2019), see Appendix B. Recall that explicit calculations only require the rst derivative of the yield function 542 gradients, which is the same in the scheme developed in Section 5 and in the formulation proposed by Choi and 543 Yoon (2019). The necking strain corresponds to the e ective plastic strain measured in the nite element located 544 at the upper right corner of the thickest section of the cell (see the yellow marker in Figure 3) when the e ective 545 plastic strain rate equals 103s1in that speci c element (similar procedure was applied in Rodr guez-Mart nez 546 et al. (2017) and N'souglo et al. (2020, 2021)). This condition indicates that plastic strain has been localized 547 in the center of the specimen and the ends of the cell are unloaded. The necking strain decreases nonlinearly 548 with the cell size, displaying a concave-upwards shape, with a large slope for small values of L0=h0, which is 549 gradually reduced with the increase of the cell length. The large values of the necking strain for short cells 550 are due to the stabilizing e ects of stress multiaxiality, which delay necking formation (Rodr guez-Mart nez 551 et al., 2017; N'souglo et al., 2020, 2021). The results obtained for the three yield criteria are very similar, 552 qualitatively and quantitatively, with slightly larger necking strains for the CPB06ex2 model. However, these 553 di erences are hardly noticeable unless the  "neckL0=h0curves are plotted in the same graph. Still, we do not 554 36 choose such representation of results because it would make the comparison between the UMAT and VUMAT 555 results less clear. Notice that the results presented in Fig. 8 bring out that the necking strains predicted by 556 the implicit and the explicit calculations are virtually the same. The key outcome is that the implicit analyses 557 performed with the UMAT, using both the numerical di erentiation scheme developed in this paper and the 558 formulation of Choi and Yoon (2019), are signi cantly cheaper in terms of computational cost, see Table 7 (the 559 CPU time is taken approximately when the necking condition is met). The maximum computational saving 560 corresponds to L0=h0= 0:5, for which the CPU time of the implicit calculations is 50 times less, and the 561 minimum computational saving corresponds to the longest cell length considered L0=h0= 6, for which the 562 implicit calculation is completed in a CPU time which is 30 times less. Moreover, note that the CPU time 563 for the calculations with the UMAT increases with L0=h0due to the increase of the number of elements with 564 the cell size. In contrast, for the explicit calculations, the computational time decreases, reaches a minimum 565 for an intermediate cell size, and then increases. In this case, the CPU time seems to be a balance between 566 the increase of the number of elements with L0=h0, and the decrease of the simulation/necking time with the 567 increase of the cell size (due to the decrease of the necking strain). The small computational time of these 568 simulations with UMAT will allow to perform ecient large scale calculations of dynamic formability problems 569 (and other localization related problems), within a range of strain rates <1000 s1which is less studied in 570 the literature, partly due to the elevated computational cost of performing explicit analyses for intermediate 571 strain rates (e.g. in our previous works, Zaera et al. (2015), N'souglo et al. (2020) and Rodr guez-Mart nez 572 et al. (2017), the minimum loading rate studied is 500 s1, 4000 s1and 5000 s1, respectively). This research 573 is left as future work. In addition, notice that for the implicit calculations the computational time obtained 574 with the numerical di erentiation scheme developed in this paper is generally between 5% and 25% lower than 575 in the case of the formulation of Choi and Yoon (2019). This is an original result of this paper which brings 576 out the e ect that the speci c form of the numerical approximation used for the second derivatives of the yield 577 function gradients has on the computational cost of these simulations. 578 37 00.20.40.60.81 0 1 2 3 4 5 6B D F Cell length, L0/h0ε0 =100 s-1 xx.Yld2004-18pNecking strain, ε neckUMAT: dynamic implicit VUMAT: dynamic explicit UMAT: dynamic implicit - Choi and Yoon (2019) (a) 00.20.40.60.81 0 1 2 3 4 5 6B D F Cell length, L0/h0ε0 =100 s-1 xx.CPB06ex2Necking strain, ε neckUMAT: dynamic implicit VUMAT: dynamic explicit UMAT: dynamic implicit - Choi and Yoon (2019) (b) 00.20.40.60.81 0 1 2 3 4 5 6B D YOON Cell length, L0/h0ε0 =100 s-1 xx.Necking strain, ε neck Yld2011-27pUMAT: dynamic implicit VUMAT: dynamic explicit UMAT: dynamic implicit - Choi and Yoon (2019) (c) Figure 8: Finite element results. Evolution of the necking strain  "neckwith the cell length L0=h0. Comparison between dynamic implicit and explicit simulations performed with UMAT and VUMAT subroutines, respectively. The imposed initial strain rate is _"0 xx= 100 s1. The results are obtained with the numerical approximation of the yield function gradients of Section 5 and the numerical di erentiation scheme proposed by Choi and Yoon (2019), which is shown in Appendix B. The calculations are performed with three yield criteria: (a) Yld2004-18p-N, (b) CPB06ex2-N and (c) Yld2011-27p-N. 38CPU time (units in kiloseconds) Yld2004-18p CPB06ex2 Yld2011-27p Cell length UMAT UMAT - Choi and Yoon (2019) VUMAT UMAT UMAT - Choi and Yoon (2019) VUMAT UMAT UMAT - Choi and Yoon (2019) VUMAT 0:5 0:186 0 :252 52 :877 0:251 0 :216 54 :118 0:265 0 :273 61 :560 1 0:260 0 :368 23 :880 0:288 0 :346 24 :600 0:321 0 :479 28 :500 2 0:409 0 :525 19 :500 0:400 0 :522 20 :820 0:519 0 :728 26 :580 3 0:581 0 :762 22 :440 0:568 0 :763 23 :760 0:728 0 :994 31 :020 4 0:753 0 :991 26 :700 0:757 0 :977 28 :980 0:918 1 :233 36 :600 5 0:940 1 :254 32 :580 0:945 1 :213 35 :160 1:143 1 :529 39 :660 6 1:107 1 :445 37 :980 1:093 1 :392 37 :800 1:657 1 :796 52 :260 Table 7: CPU time for plane strain nite element simulations performed with the three yield criteria presented in Section 2. The CPU time is taken approximately when the necking condition is met. The imposed initial strain rate is _ "0 xx= 100 s1. Comparison between dynamic implicit and explicit simulations performed with UMAT and VUMAT subroutines, respectively. The results are obtained with the numerical approximation of the yield function gradients of Section 5 and the numerical di erentiation scheme proposed by Choi and Yoon (2019), which is shown in Appendix B. 39 579 Fig. 9 shows the evolution of the necking strain  "neckwith the cell length L0=h0, for greater strain rate 580 _"0 xx= 10000 s1. Unlike in the case of 100 s1, the value of  "neckreaches a minimum for an intermediate value of 581 the cell size3, and increases for large values of L0=h0due to the in uence of inertia, which is more important 582 as the strain rate increases (Fressengeas and Molinari, 1994; Mercier and Molinari, 2003; Rodr guez-Mart nez 583 et al., 2017). The necking strains calculated with the implicit and explicit analyses are very similar, as in 584 the case of the results for the lower strain rate presented in Fig. 8. In contrast, unlike for 100 s1, the CPU 585 time is greater for the implicit calculations performed with the UMAT, see Table 8. For L0=h0= 0:5, using 586 the VUMAT leads to a computational saving which depends on the constitutive model and ranges between 587 5% and 20% with respect to the implicit calculations performed with the numerical di erentiation scheme of 588 Section 5, and for L0=h0= 6 the computational saving increases up to 250% for the Yld2011-27p yield criterion. 589 These results bring out that there is a critical loading rate from which implicit analyses are comparatively less 590 ecient than the explicit calculations. The key point is that, in comparison with the simulations performed 591 for 100 s1, the CPU time decreases two orders of magnitude for the explicit calculations, whereas it increases 592 for the implicit computations, using both the numerical di erentiation scheme developed in this paper and the 593 formulation of Choi and Yoon (2019), for most of the cell lengths investigated (it is multiplied by 2 for the 594 greater values of L0=h0). Moreover, the computational cost of both implicit and explicit calculations increases 595 with the increase of the cell size. Notice that this was not the case for the explicit calculations for the lower 596 strain rate 100 s1, for which the CPU time was minimum for an intermediate value of L0=h0. The di erence 597 is the increase of the necking strain for long wavelengths for 10000 s1, which is added to the fact that the 598 number of elements increases with L0=h0, resulting in a longer CPU time as the cell size increases. Note also 599 that for the implicit calculations, the computational time obtained with the numerical di erentiation scheme 600 developed in this paper is between 15% and 40% lower than that obtained with the formulation of Choi and 601 Yoon (2019). This comparative analysis is an original contribution of this paper. 602 40 00.20.40.60.81 0 1 2 3 4 5 6B D F Cell length, L0/h0Yld2004-18p ε0 =10000 s-1 xx.Necking strain, ε neckUMAT: dynamic implicit VUMAT: dynamic explicit UMAT: dynamic implicit - Choi and Yoon (2019) (a) 00.20.40.60.81 0 1 2 3 4 5 6B D F Cell length, L0/h0ε0 =10000 s-1 xx. CPB06ex2Necking strain, ε neckUMAT: dynamic implicit VUMAT: dynamic explicit UMAT: dynamic implicit - Choi and Yoon (2019) (b) 00.20.40.60.81 0 1 2 3 4 5 6B D F Cell length, L0/h0ε0 =10000 s-1 xx. Yld2011-27pNecking strain, ε neckUMAT: dynamic implicit VUMAT: dynamic explicit UMAT: dynamic implicit - Choi and Yoon (2019) (c) Figure 9: Finite element results. Evolution of the necking strain  "neckwith the cell length L0=h0. Comparison between dynamic implicit and explicit simulations performed with UMAT and VUMAT subroutines, respectively. The imposed initial strain rate is _"0 xx= 10000 s1. The results are obtained with the numerical approximation of the yield function gradients of Section 5 and the numerical di erentiation scheme proposed by Choi and Yoon (2019), which is shown in Appendix B. The calculations are performed with three yield criteria: (a) Yld2004-18p-N, (b) CPB06ex2-N and (c) Yld2011-27p-N. 41CPU time (units in kiloseconds) Yld2004-18p-N CPB06ex2-N Yld2011-27p-N Cell length UMAT UMAT - Choi and Yoon (2019) VUMAT UMAT UMAT - Choi and Yoon (2019) VUMAT UMAT UMAT - Choi and Yoon (2019) VUMAT 0:5 0:223 0 :444 0 :196 0:224 0 :266 0 :210 0:277 0 :344 0 :227 1 0:411 0 :755 0 :228 0:412 0 :492 0 :237 0:502 0 :650 0 :259 2 0:752 0 :983 0 :272 0:755 0 :912 0 :294 0:916 1 :533 0 :339 3 1:126 1 :336 0 :370 1:110 1 :353 0 :391 1:351 1 :853 0 :475 4 1:509 1 :890 0 :551 1:582 1 :817 0 :560 1:826 2 :368 0 :653 5 1:934 2 :436 0 :701 1:977 2 :327 0 :743 2:304 2 :928 0 :884 6 2:322 3 :556 0 :996 2:310 5 :379 1 :017 2:839 3 :576 1 :150 Table 8: CPU time for plane strain nite element simulations performed with the three yield criteria presented in Section 2. The CPU time is taken approximately when the necking condition is met. The imposed initial strain rate is _ "0 xx= 10000 s1. Comparison between dynamic implicit and explicit simulations performed with UMAT and VUMAT subroutines, respectively. The results are obtained with the numerical approximation of the yield function gradients of Section 5 (denoted with -N) and the numerical di erentiation scheme proposed by Choi and Yoon (2019), which is shown in Appendix B. 42 8. Concluding remarks 603 In this paper we have implemented in ABAQUS three advanced yield criteria using both analytical deriva- 604 tives and numerical approximations of the yield functions gradients. Speci cally, we have developed UMAT 605 and VUMAT subroutines for the Yld2004-18p, CPB06ex2 and Yld2011-27p models, with material parameters 606 values corresponding to aluminium alloy 2090-T3. We have shown that using the numerical approximation 607 of the yield functions gradients greatly simpli es the implementation of the constitutive models, which is a 608 particularly laborious task for the yield criteria considered in this paper, if the analytical expressions of the 609 stress derivatives are used. Moreover, we have performed cylindrical cup deep drawing test simulations and 610 calculations of dynamic necking localization under plane strain tension using implicit and explicit analyses in 611 ABAQUS/Standard and ABAQUS/Explicit, respectively. The nite element results have shown that, while the 612 explicit analyses performed with the VUMAT are virtually insensitive to the procedure used for the calculation 613 of the yield function gradients, the implicit analyses with the UMAT, which involve the computation of second 614 order yield function gradients, greatly bene t from using the numerical approximation for the stress derivatives, 615 leading to a decrease in the computational time up to 70% for some of the deep drawing calculations performed. 616 In addition, the dynamic necking simulations have brought out that there is a critical value of loading rate 617 below/above which the implicit/explicit analysis is computationally more ecient. For the lower strain rate 618 investigated 100 s1, the implicit analyses are between 30 and 50 times computationally cheaper than the 619 explicit simulations. In contrast, for the greater strain rate investigated 10000 s1, the implicit calculations are 620 up to 2:5 times more expensive. Note also that the implicit dynamic necking calculations have been performed 621 with both the numerical di erentiation scheme developed in this paper and the formulation proposed by Choi 622 and Yoon (2019), the latter being between 5% and 45% computationally more expensive for the nite element 623 simulations carried out in this work. 624 Acknowledgements 625 The research leading to these results has received funding from the European Research Council (ERC) 626 under the European Union's Horizon 2020 research and innovation programme. Project PURPOSE, grant 627 agreement 758056. 628 629 The authors express sincere gratitude to Dr. Oana Cazacu (University of Florida) and Dr. Guadalupe 630 Vadillo (University Carlos III of Madrid) for helpful discussions on plastic anisotropy and implementation of 631 43 the constitutive models in ABAQUS. 632 Appendix A. Analytical expressions for the stress derivatives 633 This Appendix shows the analytical expressions for the rst and second derivatives of the e ective stress 634 with respect to the components of the stress tensor for the three constitutive models presented in Section 2. 635 Appendix A.1. Yld2004-18p yield criterion 636 The rst derivative of the e ective stress with respect to the components of the stress tensor is: 637 @ @ij=@ @ Yld04@ Yld04 @ij(A.1) with: 638 @ Yld04 @ij=3X a=13X b=16X c=1@ Yld04 @~S0a@~S0 a @H0 b@H0 b @~s0c@~s0 c @ij+@ Yld04 @~S00a@~S00 a @H00 b@H00 b @~s00c@~s00 c @ij(A.2) where ~S0 a,H0 band ~s0 care the principal values, the invariants and the independent components, respectively, of 639 the tensor s0(see equation (6)). Similarly, ~S00 a,H00 band ~s00 cdenote the same quantities for the tensor s00. 640 641 The second derivative is: 642 @2 @ij@kl=@2 @ Yld04@ Yld04@ Yld04 @ij@ Yld04 @kl+@ @ Yld04@2 Yld04 @ij@kl(A.3) with: 643 44 @2 Yld04 @ij@kl=3X a=13X b=13X c=13X d=16X e=16X f=1@2 Yld04 @~S0a@~S0 b @~S0 a @H0c@H0 c @~s0e@~s0 e @ij! @~S0 b @H0 d@H0 d @~s0 f@~s0 f @kl! +@2 Yld04 @~S00a@~S00 b @~S00 a @H00c@H00 c @~s00e@~s00 e @ij! @~S00 b @H00 d@H00 d @~s00 f@~s00 f @kl! +@2 Yld04 @~S0a@~S00 b @~S0 a @H0c@H0 c @~s0e@~s0 e @ij! @~S00 b @H00 d@H00 d @~s00 f@~s00 f @kl! +@2 Yld04 @~S00a@~S0 b @~S00 a @H00c@H00 c @~s00e@~s00 e @ij! @~S0 b @H0 d@H0 d @~s0 f@~s0 f @kl! +3X a=13X b=13X c=16X d=16X e=1@ Yld04 @~S0a@2~S0 a @H0 b@H0c@H0 b @~s0 d@~s0 d @ij@H0 c @~s0e@~s0 e @kl +@ Yld04 @~S00a@2~S00 a @H00 b@H00c@H00 b @~s00 d@~s00 d @ij@H00 c @~s00e@~s00 e @kl +3X a=13X b=16X c=16X d=1@ Yld04 @~S0a@~S0 a @H0 b@2H0 b @~s0c@~s0 d@~s0 c @ij@~s0 d @kl +@ Yld04 @~S00a@~S00 a @H00 b@2H00 b @~s00c@~s00 d@~s00 c @ij@~s00 d @kl +3X a=13X b=16X c=1@ Yld04 @~S0a@~S0 a @H0 b@H0 b @~s0c@2~s0 c @ij@kl+@ Yld04 @~S00a@~S00 a @H00 b@H00 b @~s00c@2~s00 c @ij@kl(A.4) Appendix A.2. CPB06ex2 yield criterion 644 The rst derivative of the e ective stress with respect to the components of the stress tensor is: 645 @ @ij=@ @ CPB06@ CPB06 @ij(A.5) where: 646 @ CPB06 @ij=3X a=13X b=16X c=1@ CPB06 @~S0a@~S0 a @H0 b@H0 b @~s0c@~s0 c @ij+@ CPB06 @~S00a@~S00 a @H00 b@H00 b @~s00c@~s00 c @ij(A.6) and the second derivative is: 647 @2 @ij@kl=@2 @ CPB06@ CPB06@ CPB06 @ij@ CPB06 @kl+@ @ CPB06@2 CPB06 @ij@kl(A.7) with: 648 45 @2 CPB06 @ij@kl=3X a=13X b=13X c=13X d=16X e=16X f=1@2 CPB06 @~S0a@~S0 b @~S0 a @H0c@H0 c @~s0e@~s0 e @ij! @~S0 b @H0 d@H0 d @~s0 f@~s0 f @kl! +@2 CPB06 @~S00a@~S00 b @~S00 a @H00c@H00 c @~s00e@~s00 e @ij! @~S00 b @H00 d@H00 d @~s00 f@~s00 f @kl! +3X a=13X b=13X c=16X d=16X e=1@ CPB06 @~S0a@2~S0 a @H0 b@H0c@H0 b @~s0 d@~s0 d @ij@H0 c @~s0e@~s0 e @kl +@ CPB06 @~S00a@2~S00 a @H00 b@H00c@H00 b @~s00 d@~s00 d @ij@H00 c @~s00e@~s00 e @kl +3X a=13X b=16X c=16X d=1@ CPB06 @~S0a@~S0 a @H0 b@2H0 b @~s0c@~s0 d@~s0 c @ij@~s0 d @kl +@ CPB06 @~S00a@~S00 a @H00 b@2H00 b @~s00c@~s00 d@~s00 c @ij@~s00 d @kl +3X a=13X b=16X c=1@ CPB06 @~S0a@~S0 a @H0 b@H0 b @~s0c@2~s0 c @ij@kl+@ CPB06 @~S00a@~S00 a @H00 b@H00 b @~s00c@2~s00 c @ij@kl(A.8) Appendix A.3. Yld2011-27p yield criterion 649 The rst derivative of the e ective stress is: 650 @ @ij=@ @ Yld11@ Yld11 @ij(A.9) with: 651 @ Yld11 @ij=3X a=13X b=16X c=1@ Yld11 @~S0a@~S0 a @H0 b@H0 b @~s0c@~s0 c @ij+@ Yld11 @~S00a@~S00 a @H00 b@H00 b @~s00c@~s00 c @ij+@ Yld11 @~S000a@~S000 a @H000 b@H000 b @~s000c@~s000 c @ij(A.10) The second derivative reads: 652 @2 @ij@kl=@2 @ Yld11@ Yld11@ Yld11 @ij@ Yld11 @kl+@ @ Yld11@2 Yld11 @ij@kl(A.11) with 653 46 @2 Yld11 @ij@kl=3X a=13X b=13X c=13X d=16X e=16X f=1@2 Yld11 @~S0a@~S0 b @~S0 a @H0c@H0 c @~s0e@~s0 e @ij! @~S0 b @H0 d@H0 d @~s0 f@~s0 f @kl! +@2 Yld11 @~S00a@~S00 b @~S00 a @H00c@H00 c @~s00e@~s00 e @ij! @~S00 b @H00 d@H00 d @~s00 f@~s00 f @kl! +@2 Yld11 @~S000a@~S000 b @~S000 a @H000c@H000 c @~s000e@~s000 e @ij! @~S000 b @H000 d@H000 d @~s000 f@~s000 f @kl! +@2 Yld11 @~S0a@~S00 b @~S0 a @H0c@H0 c @~s0e@~s0 e @ij! @~S00 b @H00 d@H00 d @~s00 f@~s00 f @kl! +@2 Yld11 @~S00a@~S0 b @~S00 a @H00c@H00 c @~s00e@~s00 e @ij! @~S0 b @H0 d@H0 d @~s0 f@~s0 f @kl! +3X a=13X b=13X c=16X d=16X e=1@ Yld11 @~S0a@2~S0 a @H0 b@H0c@H0 b @~s0 d@~s0 d @ij@H0 c @~s0e@~s0 e @kl +@ Yld11 @~S00a@2~S00 a @H00 b@H00c@H00 b @~s00 d@~s00 d @ij@H00 c @~s00e@~s00 e @kl +@ Yld11 @~S000a@2~S000 a @H000 b@H000c@H000 b @~s000 d@~s000 d @ij@H000 c @~s000e@~s000 e @kl +3X a=13X b=16X c=16X d=1@ Yld11 @~S0a@~S0 a @H0 b@2H0 b @~s0c@~s0 d@~s0 c @ij@~s0 d @kl +@ Yld11 @~S00a@~S00 a @H00 b@2H00 b @~s00c@~s00 d@~s00 c @ij@~s00 d @kl +@ Yld11 @~S000a@~S000 a @H000 b@2H000 b @~s000c@~s000 d@~s000 c @ij@~s000 d @kl +3X a=13X b=16X c=1@ Yld11 @~S0a@~S0 a @H0 b@H0 b @~s0c@2~s0 c @ij@kl+@ Yld11 @~S00a@~S00 a @H00 b@H00 b @~s00c@2~s00 c @ij@kl+@ Yld11 @~S000a@~S000 a @H000 b@H000 b @~s000c@2~s000 c @ij@kl(A.12) Appendix A.4. Singular cases 654 The derivatives @~S0 a=@H0 b, witha= 1;2;3 andb= 1;2;3, are: 655 @~S0 a @H0 1=(~S0 a)2 (~S0a)22H0 1~S0aH0 2(A.13a) @~S0 a @H0 2=~S0 a (~S0a)22H0 1~S0aH0 2(A.13b) @~S0 a @H0 3=2 3 (~S0a)22H0 1~S0aH0 2 (A.13c) 47 which are not de ned when: 656 (~S0 a)22H0 1~S0 aH0 2= 0 (A.14) which leads to: 657 ~S0 a=H0 1q (H0 1)2+H0 2 (A.15) The comparison of previous expression with equation (8) shows that if = 0 the derivatives @~S0 2=@H0 band 658 @~S0 3=@H0 bare not de ned; and if =the derivatives which are not de ned are @~S0 1=@H0 band@~S0 2=@H0 b(see 659 Barlat et al. (2005)). This problem can be circumvented by direct calculation of the derivatives @ =@H0 b, which 660 requires to express as a function of H0 b. Nevertheless, such procedure introduces additional complications to 661 the calculation of the e ective stress derivatives shown in Appendix A.1, Appendix A.2 and Appendix A.3. 662 Obviously, the same singular cases are encountered in the derivatives @~S00 a=@H00 band@~S000 a=@H000 b. 663 48Appendix B. Numerical approximations for the stress derivatives proposed by Choi and Yoon (2019) 664 This Appendix shows the numerical approximation of the rst and second derivatives of the e ective stress with respect to the components 665 of the stress tensor proposed by Choi and Yoon (2019) to implement in ABAQUS/Standard various anisotropic constitutive models with 666 distortional hardening under associated and non-associated ow rules. Note that the di erence with respect to the numerical di erentiation 667 scheme developed in this paper resides in equation (B.3), which in the case of Choi and Yoon (2019), who took the expression from Aretz 668 (2007), involves the calculation of additional terms, increasing the computational time. 669 Appendix B.1. First derivative approximation 670 @ @ij(ij+ ij)(ijij) 2ij(B.1) Appendix B.2. Second derivatives approximation 671 @2 @ij@ij(ijij)2+ (ij+ ij) 2 ij(B.2) @2 @ij@kl(ijij;klkl)(ijij;kl+ kl)(ij+ ij;klkl) + (ij+ ij;kl+ kl) 4ijkl(B.3) 49