5 In this work, we have brought to light the e ect of material orientation on void growth. For that purpose, we have 6 performed nite element calculations using a cubic unit-cell model with a spherical void at its center and subjected to 7 periodic boundary conditions. The behavior of the material is described with an elastic isotropic, plastic orthotropic 8 constitutive model with yielding de ned by Yld2004-18p criterion (Barlat et al., 2005). We have used the multi-point 9 constraint subroutine developed by Dakshinamurthy et al. (2021) to enforce constant values of macroscopic stress 10 triaxialityTand Lode parameter Lin calculations that have been carried out for di erent stress states resulting from 11 the combination of T= 0:33, 1 and 2, with L=1, 0 and 1 (axisymmetric tension, generalized shear and axisymmetric 12 compression, respectively). Firstly, we have performed numerical simulations in which the loading directions are collinear 13 with the orthotropy axes of the material, so that the principal directions of macroscopic stress and strain are parallel. 14 Investigation of the cases for which the minor loading axis coincides either with the rolling, the transverse or the 15 normal direction, has shown that the initially spherical void turns into an ellipsoid whose rate of growth and eccentricity 16 depend on both stress state and material orientation. A key result is that for speci c material orientations the anisotropy 17 switches the e ect of Lode parameter on void growth, reversing the trends obtained for isotropic von Mises materials. 18 Secondly, we have carried out calculations using a novel strategy which consists of including angular misalignments 19 within the range 090, so that one loading direction is parallel to one of the symmetry axes of the material, 20 andis the angle formed between the other two loading directions and the second and third orthotropy axes. In fact, 21 to the authors' knowledge, these are the rst unit-cell calculations ever reported in which the material is modeled using 22 a macroscopic anisotropic yield function with prescribed misalignment between loading and material axes and, at the 23 same time, the macroscopic stress triaxiality and the Lode parameter are controlled to be constant during loading. 24 The nite element calculations have shown that the misalignment between loading and material axes makes the void 25 and the faces of the unit-cell to rotate and twist during loading. Moreover, the main contribution of this work is the 26 identi cation of an intermediate value of the angle for which the growth rate of the void reaches an extreme value 27 (minimum or maximum), so that the numerical results indicate that material orientation and angular misalignment can 28 be strategically exploited to control void growth, and thus promote or delay localization and fracture of anisotropic 29 Corresponding author. Tel. +34 916249588; Fax: +34 916249430. E-mail address:junietof@ing.uc3m.es Preprint submitted to International Journal of Plasticity October 28, 2021 2 metal products. The conclusions of this research have been shown to be valid for three di erent materials (aluminum 30 alloys 2090-T3, 6111-T4 and 6013) and selected comparisons have also been performed using two additional yield criteria 31 (CPB06ex2 and Yld2011-27p). 32 Keywords: 33 Anisotropy, Material orientation, Void growth, Stress triaxiality, Lode parameter, Unit-cell calculations. 34 1. Introduction 35 Most structural metallic alloys display plastic anisotropy. Experimental evidence of the anisotropic behavior of metals 36 was extensively discussed by Benzerga et al. (2004a), Fourmeau et al. (2013), Khadyko et al. (2019) and Benzerga et al. 37 (2019), among others. Generally, the main sources of anisotropic behavior and failure in metals are considered to be: 38 (1) plastic anisotropy resulting from the crystallographic texture, (2) topological anisotropy originated from the spatial 39 distribution of voids and second phase particles and (3) morphological anisotropy, coming from the shape and orientation 40 of voids and particles. Non-random distribution of crystal orientations (texture) is developed, for instance, in extruded 41 pro les and rolled plates. In order to model the plastic behavior of these components, the assumption of isotropy is 42 inadequate and the use of suitable anisotropic constitutive models is essential (Banabic et al., 2000; Abedrabbo et al., 43 2007). At the continuum scale, plastic anisotropy is the result of the distortion of the yield surface shape due to the 44 material micro-structural state. 45 Signi cant e orts have been made over the years to describe the plastic behavior of isotropic and anisotropic metallic 46 materials via di erent kinds of yield functions. For isotropic materials, von Mises (1913, 1928) established a quadratic 47 yield function, which is still widely used today in industry and academia. However, many face-centered-cubic materials, 48 such as some aluminum alloys, while may be isotropic, display a constitutive behavior which is better described by 49 non-quadratic yield surfaces. For instance, Hershey (1954) and Hosford (1972) proposed non-quadratic isotropic yield 50 functions with a variable exponent that considers the material's crystallographic structure. On the other hand, one of 51 the rst anisotropic yield functions for orthotropic materials was proposed by Hill (1948). The criterion of Hill (1948), 52 which is quadratic in stresses and contains six independent parameters to describe the state of anisotropy, was later 53 revisited by himself (Hill, 1979, 1993) to correct what he referred to as inevitable limitations in the range of validity 54 of the 1948 prototype . Nevertheless, the popularity of the Hill48 model persists nowadays due to its simplicity and 55 extensive implementation in commercial nite element codes. 56 Many yield functions have been developed over the years to improve the formulations of Hill (1948, 1979, 1993). 57 For instance, Barlat and Lian (1989) modi ed the isotropic model of Hosford (1972) to include in-plane anisotropy. 58 Later, Barlat et al. (1991, 1997, 2003, 2005) introduced a family of non-quadratic yield functions which use fourth 59 3 order tensors as linear multiplicative operators acting on the stress tensor to introduce material anisotropy. The linear 60 transformation approach provides extended exibility to the yield functions to capture the behaviour of materials with 61 complex anisotropic response by tailoring the number of anisotropy coecients, while ensuring the convexity of the 62 yield surface (Barlat et al., 1991; Kara llis and Boyce, 1993; Barlat et al., 1997, 2003; Bron and Besson, 2004). For 63 instance, using two linear transformations, Barlat et al. (2005) formulated the so-called Yld2004-13p and Yld2004- 64 18p yield criteria. These two anisotropic models were derived from di erent isotropic criteria and contain 13 and 18 65 anisotropic parameters, respectively. In particular, the Yld2004-18p criterion has been widely used to describe the 66 main features of the anisotropic behaviour of di erent grades of aluminum alloys (Yoon et al., 2006; Achani et al., 67 2009; Fourmeau et al., 2011; Achani et al., 2011; Tardif and Kyriakides, 2012). Years later, Aretz and Barlat (2013) 68 derived the so-called Yld2011-18p and Yld2011-27p criteria, which use two and three linear transformations of the stress 69 deviator, respectively, and involve a total of 18 and 27 anisotropic parameters. However, while using multiple linear 70 transformations enables to construct more exible and accurate yield functions, the increase in the number of material 71 parameters requires a larger number of experiments for the calibration of the yield criteria, and the implementation in 72 nite element codes becomes more laborious and generally less ecient. 73 Ductile fracture of metals and alloys has been subject to many studies over the past decades and is known to occur 74 by nucleation, growth or enlargement of voids and nal void coalescence. The voids in the material either preexist 75 (Toda et al., 2004), or nucleate at material second phase particles due to decohesion of the particle-matrix interface, 76 or by particle fracture (Maire et al., 2011). Once the voids nucleate, further plastic deformation enlarges the void size 77 and distorts its shape. During the course of stable deformation, as voids evolve, plastic strain is more or less uniformly 78 distributed in the material. At some point, deformation localizes in a narrow region of the material and adjacent voids 79 ultimately connect or coalesce with each other. 80 The eld of porous plasticity has attracted the attention of many investigators. As recent examples, Torki et al. 81 (2017) developed a micromechanics-based constitutive model to describe plasticity in solids with relatively high levels of 82 porosity. The developed model was obtained by limit analysis considering a cylindrical voided cell loaded under combined 83 tension and shear. Br unig et al. (2018) performed numerical analysis and shear-compression experiments to study the 84 e ect of stress state on damage criteria, with special emphasis on negative stress triaxialities. Torki (2019) proposed 85 a closed-form yield criterion for void growth and coalescence under combined tension and shear stress conditions and 86 validated the analytical results using cell-model nite element calculations. Becker and Callaghan (2020) examined a 87 large population of voids nucleated from randomly distributed spherical particles in order to assess the dependence of 88 load path and mean stress of void growth behavior. Very recently, Shen et al. (2020) developed a new macroscopic yield 89 criterion for a porous Drucker-Prager matrix material and validated the obtained analytical results with numerical nite 90 4 element simulations. 91 Finite element void cell simulations have been widely used for almost the last 50 years (Needleman, 1972; Tvergaard, 92 1982, 1990; Koplik and Needleman, 1988; Worswick and Pick, 1990; Becker and Smelser, 1994; Kuna and Sun, 1994; 93 Pardoen and Hutchinson, 2000) to provide useful information about the e ect of porous microstructure (void shape, 94 void volume fraction, relative void spacing), ow properties of the matrix material, and stress state on the physics of 95 ductile damage. In addition, nite element void cell simulations can o er a basis for assessing the predictive ability of 96 di erent constitutive models before applying them to more complex real-world problems. Most void cell calculations 97 reported in the literature have been performed for isotropic matrix materials and spheroidal voids, where the matrix 98 material is governed by von Mises yield function (e.g., Faleskog et al. (2000), Pardoen and Hutchinson (2003), Kim et al. 99 (2004), Benzerga et al. (2004b), Danas and Ponte Casta~ neda (2012)). Studies dealing with unit-cell simulations involving 100 plastically anisotropic matrix descriptions have been also reported in the literature, but less often than for isotropic 101 ductile materials. A notable exception is the work of Benzerga and Besson (2001), who carried out simulations using 102 axisymmetric unit-cells with Hill (1948) criterion representing the matrix material. On the other hand, three-dimensional 103 unit-cell simulations with a matrix material obeying Hill (1948) plasticity were conducted by Chien et al. (2001) and 104 Wang et al. (2004). Steglich et al. (2010) performed 3D unit-cell calculations for an anisotropic matrix material following 105 the yield criterion developed by Bron and Besson (2004). The in uence of void shape and anisotropy on void growth in 106 di erent materials was analyzed by Keralavarma and Benzerga (2010) and Keralavarma et al. (2011) using axisymmetric 107 unit-cells with the matrix material being de ned by the Hill48 yield criterion. Very recently, Legarth and Tvergaard 108 (2018) also employed the Hill48 criterion in 3D unit-cell simulations that aimed at studying the interaction between 109 initial void spacing, void shape, and plastic anisotropy in void growth evolution in porous media. All these studies 110 noticed the important e ect of the matrix anisotropy in the unit-cell response. Note that the in uence of anisotropy 111 on void growth has been also investigated within the framework of crystal plasticity elsewhere (Wan et al., 2005; Yerra 112 et al., 2010; Lebensohn and Cazacu, 2012; Han et al., 2013; Srivastava and Needleman, 2015; Dakshinamurthy et al., 113 2021). 114 Ductile fracture properties have shown marked dependence on the stress state. The deep in uence of stress triaxiality 115 in the mechanisms of damage and in the failure strain of most metals is known since the pioneering works of Rice and 116 Tracey (1969) and Gurson (1977). In addition, the third invariant of the stress deviator, generally represented by the 117 Lode parameter, has been shown to play an important role on ductile fracture, notably for low values of the stress 118 triaxiality. For large values of the stress triaxiality, the micro-voids contained in the material tend to grow fast and 119 ductile failure occurs by the large increase of void volume fraction which leads to the coalescence of neighboring voids 120 (see the reviews by Garrison and Moody (1987), Tvergaard (1990) and Benzerga and Leblond (2010)). On the other 121 5 hand, for low values of the stress triaxiality, voids tend to atten, and the void volume fraction decreases as the voids 122 become microcracks, which rotate and elongate (with signi cant changes in shape) until the nal interaction between 123 neighboring microcracks leads to material fracture. Moreover, several authors performed unit-cell analyses for shear 124 dominated loadings and showed that the voids quickly change their orientation and shape while they start to deform, 125 which is in uencing the overall behavior of the porous aggregate, as well as the nal failure of the material (McVeigh 126 et al., 2007; Leblond and Mottet, 2008; Scheyvaerts et al., 2011; Nielsen et al., 2012; Tvergaard, 2015). 127 For isotropic materials, the e ect of triaxiality, Lode parameter and void shape can be readily analyzed in 3D cell 128 simulations under prescribed loading conditions where the principal directions of the macro stress are aligned with the 129 axes of the cell and the void. However, for anisotropic matrix materials, the response of the cell also depends on the 130 orientation of the anisotropy axes. Recently, Bryhni Dhli et al. (2017a) performed 3D unit-cell analyses imposing the 131 external stress state of the cubic cell (macroscopic triaxiality and Lode parameter) and considering that the matrix 132 material is governed by the anisotropic yield criterion Yld2004-18p (Barlat et al., 2005). For the sake of simplicity, and 133 to preclude shear e ects in the unit-cell model, the material orthotropy directions were imposed to be collinear with 134 the prescribed loading conditions and with the edges of the cell (i.e., material axes parallel to both the applied stress 135 and cell axes), so that the voids initially spherical turned to be ellipsoidal during loading with semi-axes parallel to 136 the loading directions. The calculations of Bryhni Dhli et al. (2017a) showed that both the exponent of the yield 137 criterion (exponent that governs the curvature of the yield surface) and the directionality of the plastic properties have 138 an important impact on the mechanical response of the unit-cell. In this paper, we extend the analysis of Bryhni Dhli 139 et al. (2017a) performing unit-cell simulations (with three di erent sets of material parameters) in which one loading 140 direction is parallel to one of the symmetry axis of the material while the other two loading directions form an angle 141 with the second and third orthotropy axes, so that the principal directions of macroscopic stress and strain are not 142 aligned. We have investigated angular misalignments within the range 090, and showed that the growth rate 143 of the void reaches and extreme value (minimum or maximum) for an intermediate value of the angle . This conclusion 144 is an original outcome of this paper and it has been substantiated for di erent values of the stress triaxiality T= 0:33, 1 145 and 2, and the Lode parameter L=1, 0 and 1, and for di erent orientations of the unit-cell in which the major or the 146 minor loading direction coincides either with the rolling, the transverse or the normal direction of orthotropy. Note that, 147 while there are not speci c experimental results in the literature to perform a quantitative comparison with the unit-cell 148 nite element calculations presented in this work, some authors have already provided indications on the interplay 149 between plastic anisotropy, material orientation and ductility. The idea is to harness directionality of plastic properties 150 in material design, representing a paradigmatic shift in the way of thinking about the use of plastically anisotropic 151 materials in structural applications. For instance, the in uence of plastic anisotropy, yield strength and work hardening 152 6 on ductile failure was studied by Frodal et al. (2020) using nite element simulations and strain localization analyses of 153 tensile tests in di erent material orientations. Finite element calculations were performed in seven in-plane directions, 154 i.e., 0, 15, 30, 45, 60, 75and 90to the reference direction. The material behavior of the specimen was de ned by 155 the Yld2004-18p yield criterion. Plastic anisotropy was found to have a marked in uence on the tensile ductility, leading 156 to an important variation in the failure strain with the material orientation, in agreement with the trends observed in 157 the experiments of Khadyko et al. (2014, 2019) and Frodal et al. (2019). Moreover, Basu et al. (2017) determined 158 ductility of Mg alloys using smooth and round notched experimental tensile bars. A simple micromechanical model was 159 proposed to rationalize the tests results, showing that plastic anisotropy could be exploited to aid ductility. Nasim et al. 160 (2019) used an invariant parameter called the Anisotropy E ect on Ductility (AED), earlier proposed by Basu et al. 161 (2017), to tailor the formability of Mg alloys. This invariant parameter, which was proven to correctly portray physical 162 formability measurements, was calculated using tension and compression experimental tests of rolled Mg AZ31 alloy. 163 In order to create a particular AED parameter or formability, an automated inverse optimization strategy was used to 164 predict the routes needed to attain target amount of formability. 165 The numerical calculations presented in this paper corroborate the results obtained in the works discussed in the 166 paragraph above, and reinforce the idea that material orientation and angular misalignment can be strategically ex- 167 ploited to control void growth, and thus promote or delay localization and fracture of anisotropic metal products. The 168 manuscript is organized as follows. Section 2 presents the elasto-plastic constitutive model used to describe the mechan- 169 ical behavior of the matrix material. Material yielding is modeled with the Yld2004-18p criterion (Barlat et al., 2005) 170 with parameters corresponding to aluminum alloys 2090-T3, 6111-T4 and 6013. Moreover, Section 3 shows the unit-cell 171 nite element model developed in ABAQUS/Standard (2016), and Section 4 shows calculations in which the principal 172 directions of macroscopic stress and strain are collinear (Section 4.1), and also calculations with prescribed misaligment 173 between the loading and material axes (Section 4.2). A summary of the main ndings of the paper is given in Section 5. 174 2. Constitutive framework 175 We consider an elastic isotropic, plastic orthotropic constitutive model with yielding described by Yld2004-18p 176 criterion (Barlat et al., 2005). This advanced yield criterion is able to successfully describe the main features of the 177 anisotropic behaviour of di erent grades of aluminum alloys (Fourmeau et al., 2011; Achani et al., 2011; Tardif and 178 Kyriakides, 2012). 179 180 The rate of deformation tensor dis assumed to be the sum of an elastic part deand a plastic part dp: 181 7 d=de+dp (1) where the elastic part is related to the rate of the stress by the following linear hypo-elastic law: 182 O=L:de (2) whereOis an objective derivative of the Cauchy stress tensor (which corresponds to the Jaumann derivative in 183 ABAQUS/Standard), and Lis the tensor of isotropic elastic moduli given by: 184 L=E 1 +I0+E 3(12)1 1 (3) withEbeing the Young's modulus, the Poisson's ratio, 1the unit second-order tensor and I0the unit deviatoric 185 fourth-order tensor. 186 187 The yield function is de ned as: 188 =IIIX i=IIIIX j=I ~s0 i~s00 j a(4) where ~s0 iand ~s00 jare the principal values of the second order deviatoric tensors ~s0and~s00, which are de ned by two linear 189 transformations: 190 ~s0=C0s (5a) ~s00=C00s (5b) 191 where sis the deviatoric part of the Cauchy stress tensor , and C0andC00are the matrices which contain the anisotropy 192 coecients: 193 8 C0=2 6666666666666640c0 12c0 13 0 0 0 c0 21 0c0 23 0 0 0 c0 31c0 32 0 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(6a) C00=2 6666666666666640c00 12c00 13 0 0 0 c00 21 0c00 23 0 0 0 c00 31c00 32 0 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(6b) 194 withc0 ijandc00 ijbeing material parameters. Note that in equation (5) we have taken the order of the components 195 of the column vector representing the deviatoric stress tensor to be s=fsxxsyyszzsyzszxsxyg, with the Cartesian 196 coordinate system ( x;y;z ) being associated to the orthotropy axes of the material ( x,yandzaxes correspond to the 197 rolling, transverse and normal directions, respectively). 198 199 The ordered principal values of the tensor ~s0are: 200 ~s0 I= 2q (H0 1)2+H0 2cos 3 +H0 1 (7a) ~s0 II= 2q (H0 1)2+H0 2cos+ 4 3 +H0 1 (7b) ~s0 III= 2q (H0 1)2+H0 2cos+ 2 3 +H0 1 (7c) 201 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, 2ndand 202 3rdinvariant of ~s0(Barlat et al., 2005). The same procedure is followed to obtain the principal values of the tensor ~s00. 203 204 9 The e ective stress associated to this yield criterion is: 205 = 1=a (8) with= 4, so that  reduces to the yield stress in the rolling direction (see equation (9)). 206 207 Moreover, the yield condition is expressed as: 208 f= Y= 0 (9) whereYis the yield stress of the material in the rolling direction (direction x), which is considered to evolve following 209 a power-type relation: 210 Y=0("0+ "p)n(10) where "p=Rt 0_"pd is the e ective plastic strain and _"pis the e ective plastic strain rate. Moreover, 0,"0, andnare 211 material parameters. 212 213 Assuming an associated plastic ow rule, the plastic part of the rate of deformation tensor is: 214 dp=_@ @(11) where _is the rate of plastic multiplier. 215 216 In addition, the work conjugacy relation: 217 :dp= _"p (12) 10 Symbol Property and units Value 0 Initial density (kg =m3) 2700 G Elastic shear modulus (GPa), Eq. (3) 26 :92 K Bulk modulus (GPa), Eq. (3) 58 :33 Table 1: Initial density and elastic constants for the three aluminium alloys considered in this work: 2090-T3, 6111-T4 and 6013. leads to the identity: 218 _"p=_ (13) The formulation of the constitutive model is completed with the Kuhn-Tucker loading-unloading conditions: 219 _>0; f 60; _f= 0 (14) and the consistency condition during plastic loading: 220 _f= 0 (15) 221 The parameters values of the constitutive model corresponding to aluminium alloys 2090-T3, 6111-T4 and 6013 are 222 shown in Tables 1-7. Fig. 1 shows the projections in the deviatoric -plane of the yield loci, where I,IIandIIIare 223 the principal values of the tensor corresponding to the rolling, transverse and normal directions, respectively. This 224 representation of the yield loci makes apparent the anisotropy and the dependence of yielding on material orientation. 225 For instance, note that for aluminium alloy 2090-T3 the minimum and maximum yield stresses correspond to uniaxial 226 tension in the transverse and normal directions, respectively. On the other hand, notice that the Yld2004-18p criterion 227 reduces to isotropic von Mises (1928) plasticity when c0 ij=c00 ij= 1 anda= 2 (Barlat et al., 2005). 228 11 Aluminium alloy 2090-T3 Strain hardening parameters 0Material parameter (MPa), Eq. (10) 646 "0 Material parameter, Eq. (10) 0 :025 n Material parameter, Eq. (10) 0 :227 Table 2: Strain hardening parameters for aluminium alloy 2090-T3. Data after Yoon et al. (2006) and Cvitani c et al. (2008). Aluminium alloy 2090-T3 Anisotropy parameters c0 12c0 13c0 21c0 23c0 31c0 32c0 44c0 55c0 66 0:06989 0:93640 0:07914 1:00306 0:524741 1:36318 1:02377 1:06906 0:95432 c00 12c00 13c00 21c00 23c00 31c00 32c00 44c00 55c00 66 0:98117 0:47674 0:57531 0:86682 1:145010:07929 1:05166 1:14700 1:40462 Degree of homogeneity of the yield function a 8 Table 3: Yld2004-18p yield criterion. Parameters values for aluminium alloy 2090-T3. Data after Barlat et al. (2005). Aluminium alloy 6111-T4 Strain hardening parameters 0Material parameter (MPa), Eq. (10) 503 :7 "0 Material parameter, Eq. (10) 0 n Material parameter, Eq. (10) 0 :233 Table 4: Strain hardening parameters for aluminium alloy 6111-T4. Data after Kim et al. (2010). Aluminium alloy 6111-T4 Anisotropy parameters c0 12c0 13c0 21c0 23c0 31c0 32c0 44c0 55c0 66 0:637630 0:941853 0:693134 1:009620 0:608808 0:730685 1:004578 0:818717 0:616609 c00 12c00 13c00 21c00 23c00 31c00 32c00 44c00 55c00 66 1:281632 1:529747 1:184135 1:4761910:671240:96936 1:005348 1:018549 1:245460 Degree of homogeneity of the yield function a 8 Table 5: Yld2004-18p yield criterion. Parameters values for aluminium alloy 6111-T4. Data after Kim et al. (2010). 12 00.20.40.60.811.2 0306090 120 210 240 270300330Material: AA 2090-T3 σ II/σ Y 0.911.030.91 0.90 0.99 1 1 0.91 1.030.910.900.99σ I/σ Y σ III/σ Y (a) 00.20.40.60.811.2 0306090 120 210 240 270300330σ I/σ Yσ II/σ Y σ III/σ Y10.950.990.92 0.98 0.99 1 0.95 0.99 0.920.980.99Material: AA 6111-T4 (b) 00.20.40.60.811.2 0306090 120 210 240 270300330σ II/σ Y 0.98Material: AA 6013 σ I/σ Y10.93 0.93 10.940.980.9810.94 0.98 1 σ III/σ Y (c) Figure 1: Yld2004-18p yield criterion. Projection in the deviatoric -plane of the yield locus for aluminium alloys: (a) 2090-T3, (b) 6111-T4 and (c) 6013. 13 Aluminium alloy 6013 Strain hardening parameters 0Material parameter (MPa), Eq. (10) 556 :06 "0 Material parameter, Eq. (10) 0 :0062 n Material parameter, Eq. (10) 0 :201 Table 6: Strain hardening parameters for aluminium alloy 6013. Data after Ha et al. (2018). Aluminium alloy 6013 Anisotropy parameters c0 12c0 13c0 21c0 23c0 31c0 32c0 44c0 55c0 66 1:0717 1:1541 1:1746 1:1832 0:8482 1:0672 0:9543 1:0142 1:2007 c00 12c00 13c00 21c00 23c00 31c00 32c00 44c00 55c00 66 0:8119 0:6928 0:7644 0:6693 1:0239 0:9485 1:0467 0:9846 0:7331 Degree of homogeneity of the yield function a 8 Table 7: Yld2004-18p yield criterion. Parameters values for aluminium alloy 6013. Data after Ha et al. (2018). 3. Finite element model 229 We consider a porous material with periodic microstructure modeled by an array of representative volume elements 230 idealized as cubic unit-cells with a spherical void of initial radius R0=a0=b0=c0located at their center, see Fig. 231 2, wherea0,b0andc0denote the initial radius of the void along the axes x1,x2andx3(see Section 4.1). The initial 232 size of the cubic cell is L0L0L0, so that the initial void volume fraction, de ned as the ratio between the initial 233 volume of the void and the initial volume of the cell, is f0=Vvoid 0=Vcell 0= 0:0042. Moreover, the current volumes of 234 the void and the cell are denoted by VvoidandVcell, respectively, so that the current volume of the matrix material is 235 Vmatrix=VcellVvoid. The nite element model has been implemented in ABAQUS/Standard (2016), with the initial 236 con guration of the unit cell being de ned by the domain 0 x1L0, 0x2L0and 0x3L0. The origin of 237 the Lagrangian Cartesian coordinate system ( x1;x2;x3) is located at the bottom left rear corner of the cell, see Fig. 2. 238 The loading directions are determined by the axes x1;x2andx3(see below). The behavior of the matrix material is 239 de ned with the constitutive framework presented in Section 2. 240 We have used the multi-point constraint subroutine developed by Dakshinamurthy et al. (2021) to enforce constant 241 and controlled values of the macroscopic stress triaxiality Tand the macroscopic Lode parameter Lduring the calcu- 242 lations. Note that the approach of imposing TandLthroughout the deformation for the analysis of void growth in 243 14 Figure 2: (a) Unit-cell nite element model and (b) cut-view for x3=L0 2with I, IIand IIIbeing the principal values of the macroscopic stress tensor. The loading directions are aligned with the principal directions of the macroscopic stress tensor, so that the major loading direction corresponds to the principal stress direction associated to  I, and the minor loading direction corresponds to the principal stress direction associated to  III. unit-cell calculations follows from the earlier works of Zhang et al. (2001), Kim et al. (2004), Srivastava and Needleman 244 (2015) and Vadillo et al. (2016), among others. We introduce the macroscopic stress tensor through the following 245 relation with the microscopic (local) Cauchy stress tensor : 246 =1 VcellZ VcelldVcell(16) Moreover, the macroscopic e ective stress , the macroscopic stress triaxiality Tand the macroscopic Lode parameter 247 Lare de ned as: 248  =r 3 20:0; T =h ; L =2IIIIII IIII(17) where: 249 0=h1; h=I+ II+ III 3(18) with I, IIand IIIbeing the principal values of the macroscopic stress tensor with  IIIIII(i.e., I, II 250 and IIIare the major, intermediate and minor principal values of the macroscopic stress tensor, respectively). The 251 15 loading directions are aligned with the principal directions of the macroscopic stress tensor, so that the major loading 252 direction corresponds to the principal stress direction associated to  I(parallel to x1), and the minor loading direction 253 corresponds to the principal stress direction associated to  III(parallel to x3), see Fig. 2. 254 255 The volume average of the microscopic (local) logarithmic strain tensor , expressed as: 256 "=1 VmatrixZ VmatrixdVmatrix(19) leads to the macroscopic e ective strain  ", as de ned by Srivastava and Needleman (2015): 257 "=r 2 3"0:"0;"0=""h1; "h="I+"II+"III 3(20) where"I,"IIand"IIIare the principal values of the tensor ". Expressions (17) and (20) are scalar measures of 258 macroscopic stress and macroscopic strain, respectively, consistent with the de nitions used, for instance, by Srivastava 259 and Needleman (2015), Bryhni Dhli et al. (2017a) and Dakshinamurthy et al. (2021), in unit-cell calculations analyses 260 with di erent anisotropic plasticity theories. Note that expressions (17) and (20) do not determine the speci c e ective 261 stress and strain in a material point for this yield criterion |these are given by equations (8) and (12)|, but they are 262 rather used for the de nition of macroscopic stress scalars (macroscopic Lode and macroscopic triaxiality) and for the 263 graphical representation and interpretation of results. 264 265 In Section 4, we carry out calculations for three di erent values of macroscopic triaxiality T= 0:33, 1 and 2, and 266 macroscopic Lode parameter L=1, 0 and 1, and for macroscopic e ective strains  "up to 0:6. Similar range of values of 267 macroscopic triaxiality and Lode parameter has been studied elsewhere (e.g., Bryhni Dhli et al. (2017a,b)). Moreover, 268 recall thatL=1;0 and 1 represent states of axisymmetric tension, generalized shear and axisymmetric compression, 269 respectively. Note that for L=1 we have that  I>II= III, forL= 0 the principal values of the macroscopic 270 stress tensor are such that  I>II=I+ III 2>III, and forL= 1 we have that  I= II>III. We are 271 aware that the mechanical behavior of anisotropic materials depends on both the principal stresses and the principal 272 stress directions, so that when presenting the results in Section 4 we always indicate the orientation of the loading axes 273 (x1;x2;x3) with respect to the material orthotropy directions ( x;y;z ). Note that providing information about the stress 274 eigenvalues and eigenvectors amounts to giving the six independent stress components in arbitrary basis, namely, the 275 basis related to the material orthotropy axes. In this regard, note that in all the calculations performed in this work the 276 16 major and minor loading directions are parallel to the x1andx3axes, respectively. We will further elaborate on this 277 issue in Section 4. In addition, following the procedure described in Appendix A of Dakshinamurthy et al. (2021), we 278 have applied periodic boundary conditions to the unit-cell, ensuring that the relative displacement between opposing 279 nodes is the same (i.e., the displacement of opposing nodes is coupled, e.g., see Fig. 24 in Section 4.2). Moreover, a 280 mesh convergence study showed that a discretization consisting of 86096 eight-node tri-linear elements with reduced 281 integration and hourglass control (C3D8R in ABAQUS notation) yields numerical results for the evolution of the void 282 growth which are virtually mesh independent (for the range of macroscopic e ective strains and void volume fractions 283 investigated in this paper). This spatial discretization was used in all the unit-cell simulations performed in this work. 284 Moreover, the elastic-plastic constitutive model presented in Section 2 has been implemented in ABAQUS/Standard 285 (2016) through a UMAT subroutine, using the integration algorithm based on numerical approximation of the yield 286 function gradients developed by Hosseini and Rodr guez-Mart nez (2021). 287 4. Results 288 In Section 4.1, we consider that the loading directions are aligned with the orthotropy axes of the material and 289 investigate the cases for which the minor loading axis x3coincides either with the rolling, the transverse or the normal 290 direction (recall that the minor loading direction corresponds to the principal direction associated to  III). The results 291 correspond to AA 2090-T3. Similar calculations were performed by Bryhni Dhli et al. (2017a) using the Yld2004-18p 292 criterion with parameters corresponding to six generic crystallographic textures representative of di erent anisotropic 293 face-centered-cubic metals, including brass and copper. In Section 4.2, we investigate the in uence of the misalignment 294 between loading and orthotropy axes on void growth, exploring di erent material orientations. Results are presented 295 for di erent aluminium alloys: 2090-T3, 6111-T4 and 6013. The main contribution of this paper to the eld of nite 296 element void cell calculations is to consider anisotropic materials modeled with a macroscopic yield criterion in which 297 loading and materials axes are not aligned. In addition, for the sake of completeness, Appendix A presents calculations 298 performed with two additional yield criteria, CPB06ex2 (Plunkett et al., 2008) and Yld2011-27p (Aretz and Barlat, 299 2013), calibrated to model the behavior of AA 2090-T3. 300 301 Recall that the loading directions are xed relative to the Lagrangian coordinate system ( x1;x2;x3) associated to 302 the undeformed cell (see Fig. 2). Moreover, the orthotropy axes ( x;y;z ) are associated to the material (locally), with 303 prescribed directions in the undeformed con guration that can rotate if the material rotates during loading (as in the 304 calculations in which loading and material axes are not aligned). The orientation of the orthotropy axes given in the 305 paper corresponds to the undeformed con guration. Moreover, note that we have not included any failure criterion in 306 17 the simulations, and we have not analyzed localization of plastic deformation leading to coalescence and subsequent 307 material fracture (the reader is referred to the works of Teko glu et al. (2015) and Reboul et al. (2020) where localization 308 and coalescence are analyzed for isotropic materials using unit-cell nite element calculations). The speci c in uence of 309 plastic anisotropy and material orientation on localization of deformation near the void, coalescence and fracture will 310 be investigated in a future work. 311 4.1. Alignment between loading and material axes 312 We perform a parametric analysis on the in uence of macroscopic stress triaxiality (Section 4.1.1), Lode parameter 313 (Section 4.1.2) and material orientation (Section 4.1.3) on void growth, paying speci c attention to the evolution of 314 the volume and the geometry of the void during loading. Note that, following Bryhni Dhli et al. (2017a), we could 315 have exploited the symmetry of the problem and performed the calculations modeling 1 =8 of the unit-cell. However, 316 we have considered the whole unit-cell in order to be consistent with the calculations performed in Section 4.2 in which 317 the misalignment between the principal directions of stress and strain breaks the symmetry of the problem, preventing 318 the use of a simpli ed geometric model. As mentioned before, all the results reported in this section correspond to 319 aluminium alloy 2090-T3. The notation for the material orientations investigated in the calculations performed in this 320 section isxIIIyIzII,yIIIzIxIIandzIIIxIyII, where the rst, second and third letters indicate the orthotropy 321 axes corresponding to the minor, major and intermediate loading directions, respectively, see Fig. 3. 322 Figure 3: Schematic representation for the material orientations investigated in the calculations reported in this section, xIIIyIzII, yIIIzIxIIandzIIIxIyII, where the rst, second and third letters indicate the orthotropy axes corresponding to the minor, major and intermediate loading directions, respectively. 4.1.1. The in uence of stress triaxiality 323 Fig. 4 shows the evolution of the normalized void volume fraction f=f0with the macroscopic e ective strain  "for 324 three values of the macroscopic stress triaxiality T= 0:33, 1 and 2. The Lode parameter is L=1 so that the loading 325 18 in directions x2andx3is the same. The void volume fraction, f, is calculated using the Quickhull algorithm (Barber 326 et al., 1996) to compute the smallest convex set containing the coordinates of the nodes that form the surface of the 327 void. The convex hull is determined at each time increment, thus obtaining the evolution of the void volume fraction 328 during loading. Results are shown for the von Mises material, Fig. 4(a), and for the material modeled with Yld2004-18p 329 yield criterion with orientation xIIIyIzII, Fig. 4(b), so that the minor loading direction x3is parallel to the rolling 330 directionx. Recall from Section 2 that the results for the isotropic von Mises (1928) material are obtained imposing 331 c0 ij=c00 ij= 1 anda= 2 in the Yld2004-18p criterion, see equations (6) and (8). For both yield functions the volume 332 of the void increases with the e ective strain and the triaxiality. Note that, while for T= 1 and 2 the void volume 333 fraction increases slightly faster for the von Mises material (check the intersection of the red and green dashed curves 334 with the axis f=f0= 1:7), forT= 0:33 the increase of porosity is slightly greater for the Yld2004-18p criterion (check the 335 intersection of the solid black curves with the axis  "= 0:4), bringing out that the quantitative e ect of the triaxiality 336 on void growth depends on the plastic anisotropy. Note that we have also checked that the rate of growth of the 337 void depends on the material orientation, i.e., quantitative di erences are obtained for the simulations performed with 338 yIIIzIxIIandzIIIxIyII. Nevertheless, these results are not shown here for the sake of brevity. 339 11.11.21.31.41.51.61.7 0 0.1 0.2 0.3 0.40.33 1 2Normalized void volume fraction, f/f0 Macroscopic effective strain, εT=0.33 T=1 T=2 L=-1von Mises criterion Material: AA 2090-T3 (a) 11.11.21.31.41.51.61.7 0 0.1 0.2 0.3 0.40.33 1 2Normalized void volume fraction, f/f0 Macroscopic effective strain, εT=0.33 T=1 T=2 L=-1Yld2004-18p criterion Material: AA 2090-T3 Material orientation: xΣ yΣ zΣ III I II (b) Figure 4: In uence of the macroscopic stress triaxiality on void growth for AA 2090-T3. Unit-cell nite element calculations for L=1 (axisymmetric tension) and three values of stress triaxiality T= 0:33, 1 and 2. Evolution of the normalized void volume fractionf=f0with the macroscopic e ective strain  "for: (a) von Mises (1928) yield criterion and (b) Yld2004-18p yield criterion (Barlat et al., 2005) with material orientation xIIIyIzII. The dashed yellow line corresponds to the macroscopic e ective strain  "= 0:2. (For interpretation of the