This paper investigates the steady-state dynamic radial expansion of a pressurized circular cylindrical cavity in an inﬁnite porous medium modeled with the constitutive framework developed by Monchiet et al. (2008), which considers the material to display a periodic porous microstructure with spheroidal voids and matrix described by the orthotropic yield criterion of Hill (1948). For that purpose, we have extended the formulation of dos Santos et al. (2019b) to consider oblate and prolate voids, which allows to assess the role of the initial voids shape on the elastoplastic-anisotropic ﬁelds that develop near the cavity. The theoretical development follows the cavity expansion formalism of Cohen and Durban (2013) and employs the artiﬁcial viscosity approach of Lew et al. (2001) to avoid singularities in the ﬁeld variables due to the formation of plastic shock waves. The main outcome of this work is a relationship between the critical cavity expansion velocity for which plastic shocks emerge and the initial aspect ratio of the spheroidal voids. The results show that the formation of shocks is delayed for oblate voids, in comparison with spherical and prolate voids. These ﬁndings have been substantiated for diﬀerent anisotropic behaviors and initial void volume fractions. Keywords: Cylindrical cavity expansion, Matrix anisotropy, Ductile materials, Spheroidal voids, Plastic shocks 1. Introduction 1 Additive manufacturing (3D printing) has been recognized as a promising process to aid both design and 2 production activities, such that, in addition to rapid prototyping, it has also been employed to produce functional 3 products, including metallic parts (Thompson et al., 2016). It is well known that printed materials have intrinsic 4 features, such as porosity and anisotropy, which are expected to inﬂuence the mechanical behavior of additive 5 manufactured metallic components (de Terris et al., 2019; Sola and Nouri, 2019). Driven by the growing employ- 6 ment of this emerging technology, large eﬀorts are being made in order to understand the mechanical behavior of 7 additive manufactured parts, which is signiﬁcantly less studied in comparison with the mechanical performace of 8 components and structures fabricated using conventional processes, such as forging and/or machining (Gray et al., 9 ∗Corresponding author. E-mail address: tiago.santos@ufsm.br Preprint submitted to International Journal of Mechanical Sciences 9th October 2021 2 2017; du Plessis et al., 2020). In particular, when it comes to impact applications, the interplay between material 10 porosity, plastic anisotropy and inertia on the perforation/penetration mechanics of printed structures needs to 11 be assessed. In the words of Costas et al. (2021), additive manufacturing has a large unrealised potential in the 12 ﬁeld of ballistic protection, with the largest drawback being the uncertainty surrounding the mechanical properties 13 of additive manufactured materials under impulsive loading . 14 Currently, there are only a few experimental studies concerning the perforation/penetration of printed struc- 15 tures. For instance, Kristoﬀersen et al. (2020) investigated the ballistic perforation resistance of additive manu- 16 factured aluminum plates employing both experimental and numerical analysis. The authors observed negligible 17 diﬀerence in ballistic performance between the additive manufactured and die-cast materials. As an important 18 outcome, they concluded that a printed material may have similar or even better ballistic properties than a 19 traditionally-manufactured equivalent material. Moreover, Costas et al. (2021) addressed the ballistic properties 20 of a selective laser melting additive manufactured maraging high-strength steel and demonstrated the potential of 21 this material in protective applications, which could be even better had the presence of large voids and particles 22 been reduced. 23 In absence of experimental data, theoretical and numerical analysis are essential to advance in the knowledge 24 of the mechanisms that control the perforation/penetration of printed materials. In this context, the cavity 25 expansion theory has been accepted as a simple and fairly accurate approach to understand the underlying 26 mechanics of penetration/perforation problems (Forrestal et al., 1990, 1991; Durban and Masri, 2004; Masri and 27 Durban, 2005). The key advantage of the cavity expansion theories is to provide simple expressions to predict the 28 pressure-velocity relationship during penetration and also to estimate the mechanical ﬁelds close to the cavity. 29 For example, the cylindrical cavity analogy was used by Forrestal et al. (1990) and Forrestal and Warren (2009) 30 to obtain closed-form perforation equations for the ballistic limit and residual velocities for conical and ogival 31 nose rigid projectiles penetrating aluminum targets. Masri and Durban (2009) studied deep penetration problems 32 using the dynamic cylindrical cavitation model, and developed a practical procedure for estimating penetration 33 the depth of rigid axisymmetric projectiles, using diﬀerent plasticity theories to model the mechanical behavior of 34 the targets, including von Mises and Tresca. Following the formalism of Cohen and Durban (2013), the model of 35 Masri and Durban (2009) was extended by dos Santos et al. (2019b) in order to study cylindrical cavity expansion 36 in orthotropic porous materials modeled with the constitutive framework of Benzerga and Besson (2001). The 37 authors investigated the inﬂuence of key material features, such as plastic anisotropy, initial porosity and strain 38 hardening, on the plastic shock wave emergence and on the elastoplastic ﬁelds that develop near the cavity. Among 39 other conclusions, it was found that the formation and location of shock waves is practically independent of the 40 anisotropic nature of the matrix material and the initial porosity was shown to promote the emergence of shocks 41 for lower cavitation velocities. 42 Although the framework developed by dos Santos et al. (2019b) accounted for plastic anisotropy, the adopted 43 3 porous plasticity model (Benzerga and Besson, 2001) was based on the assumption of persisting spherical voids. 44 However, additive manufactured materials also contain non-spherical pores (de Terris et al., 2019; Sola and Nouri, 45 2019), and the shape of the voids aﬀects the response of porous ductile materials (Gologanu et al., 1993, 1994, 46 1997), specially when the matrix material is anisotropic (Monchiet et al., 2008; Keralavarma and Benzerga, 2010; 47 Morin et al., 2015; Tvergaard and Nyvang Legarth, 2019). Speciﬁc micromechanics-based models have been 48 proposed in order to consider the interplay between anisotropy and voids shape in the mechanical response of 49 the material. For example, Monchiet et al. (2008) extended the Gologanu-Leblond-Deveaux’s model (Gologanu 50 et al., 1993, 1994), assuming spheroidal (prolate/oblate) voids embedded in a Hill (1948) matrix material, so that 51 the authors obtained a closed-form macroscopic yield function depending on the void volume fraction, the void 52 shape, and the directionality of the plastic properties. Moreover, Keralavarma and Benzerga (2010) extended the 53 Gologanu-Leblond-Perrin-Deveaux’s model (Gologanu et al., 1997), which employed diﬀerent kinematic conditions 54 in comparison to Gologanu et al. (1993, 1994), and derived plastic constitutive relations for anisotropic porous 55 materials containing spheroidal voids in an plastically anisotropic matrix (Hill, 1948). Unlike Monchiet et al. 56 (2008), the authors also proposed evolution equations for the void volume fraction, the void shape, and the 57 orientation. The model of Keralavarma and Benzerga (2010) was assessed in a subsequent work using unit cell 58 ﬁnite element simulations (Keralavarma and Benzerga, 2015). This constitutive model was shown to provide a 59 robust framework to be used in large-scale simulations of ductile fracture accounting for the combined eﬀects of 60 stress triaxiality, void shape, void orientation, matrix orthotropy, and loading path (see Kweon et al. (2016)). 61 Morin (2015) compared the models of Monchiet et al. (2008) and Keralavarma and Benzerga (2010) against unit 62 cell ﬁnite element calculations. The results demonstrated that, when spheroidal voids are considered, both models 63 provide more precise predictions for the eﬀective yield loci, in comparison with porous plasticity models based on 64 the assumption of persisting spherical voids (Gurson, 1977; Benzerga and Besson, 2001). Moreover, Morin et al. 65 (2015) extended Madou-Leblond’s proposal (Madou and Leblond, 2012a,b) and developed a criterion for general 66 ellipsoidal cavities embedded in a plastically anisotropic material also obeying Hill’s criterion. In contrast to the 67 works of Monchiet et al. (2008) and Keralavarma and Benzerga (2010), the study of Morin et al. (2015) was not 68 restricted to spheroidal voids, so that ellipsoidal voids with arbitrary shapes were considered at the expense of a 69 more complex formulation. 70 In this work, we extend the dynamic cavity expansion model of dos Santos et al. (2019b) in order to assess 71 the coupled eﬀects of voids shape and plastic anisotropy on the elastoplastic ﬁelds that develop near the cavity, 72 and on the formation of plastic shocks at high cavitation velocities. The theoretical development follows the 73 cavity expansion framework of Cohen and Durban (2013) using the artiﬁcial viscosity method (Lew et al., 2001) 74 in order to avoid singularities related to the emergence plastic shock waves. The main novelty of this work is to 75 model the material behavior with the yield criterion of Monchiet et al. (2008), which has a simpler formulation 76 compared with the model of Keralavarma and Benzerga (2010), being especially suited for its implementation in 77 4 analytical models, with the material microstructure, void volume fraction and void shape, evolving according to 78 the constitutive equations of Gologanu et al. (1993, 1994). As a main result, this work shows the inﬂuence of voids 79 initial aspect ratio and matrix anisotropy on the critical cavity expansion velocity which leads to the formation 80 of plastic shocks. 81 2. Constitutive framework 82 The mechanical behavior of the material is described using the constitutive model for anisotropic porous metals 83 developed by Monchiet et al. (2008), who extended the formulation of Gologanu et al. (1993, 1994) to consider 84 materials with matrix described by the orthotropic yield criterion of Hill (1948). 85 86 The model of Monchiet et al. (2008) assumes that the material displays a periodic porous microstructure 87 that can be approximated by an array of representative volume elements idealized as spheroidal unit-cells which 88 contain a confocal spheroidal void, see Fig. 1. The semi-axes of the spheroidal void and the spheroidal unit-cell 89 are(a1,b1)and(a2,b2), respectively, and the corresponding eccentricities are deﬁned as: 90 e1=c a1ande2=c a2(prolate) (1) e1=c b1ande2=c b2(oblate) (2) wherec=/radicalBig/vextendsingle/vextendsinglea2 1−b2 1/vextendsingle/vextendsingle=/radicalBig/vextendsingle/vextendsinglea2 2−b2 2/vextendsingle/vextendsingleis the focal distance. 91 92 Moreover, the void volume fraction, f, and the void shape parameter, S, can be expressed as: 93 f=a1b2 1 a2b2 2(3) S= ln/parenleftbigga1 b1/parenrightbigg (4) so thatS >0stands for a prolate void and S <0for an oblate void. Note that given the void volume fraction f 94 and the void shape parameter S, the eccentricities can also be calculated using the following expressions (Gologanu 95 et al., 1997): 96 e2 1=  1−1 w2(prolate) 1−w2(oblate)ande3 2/parenleftbig1−e2 2/parenrightbign=fe3 1/parenleftbig1−e2 1/parenrightbignwithn=  1 (prolate) 1/2 (oblate)(5) wherew=a1 b1. 97 5 x2x3 x1 x1x3 x2 a1a2b2 b1 a2b2 b1 a1 Figure 1: Spheroidal unit-cell containing a confocal spheroidal void referred to a Cartesian coordinate system (x1,x2,x3), withx3as the symmetry axis: (a) prolate void and (b) oblate void. 2.1. General equations for elastic/plastic materials 98 The rate of deformation tensor dis additively decomposed into an elastic part, de, and a plastic part, dp, as 99 follows: 100 d=de+dp(6) The elastic rate of deformation tensor is related to the Cauchy stress tensor, σ, by the following hypo-elastic 101 law: 102 Oσ=C:de(7) withO ()being an objective derivative and Cthe isotropic elastic tensor: 103 C=E 1 +νI/prime+E 3 (1−2ν)1⊗1 (8) whereEis the Young’s modulus, νis the Poisson’s ratio, 1is the unit second-order tensor and I/primeis the unit 104 deviatoric fourth-order tensor. 105 106 The plastic rate of deformation tensor is related to the Cauchy stress tensor by an associated ﬂow rule: 107 dp=˙λ∂Φ ∂σ(9) where Φis the ﬂow potential, see Eq. (12), and ˙λis a non-negative rate of plastic multiplier satisfying the 108 loading/unloading Kuhn-Tucker conditions: 109 6 ˙λ≥0; Φ≤0; ˙λΦ = 0 (10) and the consistency condition: 110 ˙λ˙Φ = 0 (11) where the superimposed dot denotes diﬀerentiation with respect to time. 111 2.2. Monchiet et al. (2008) yield criterion 112 Using a limit analysis approach, Monchiet et al. (2008) derived the following ﬂow potential: 113 Φ =/parenleftBigg ¯σe σy/parenrightBigg2 + 2 (1 +g) (f+g) cosh/parenleftBigg κσp σy/parenrightBigg −(1 +g)2−(f+g)2(12) wheregis a void shape parameter, see Eq. (23), κis a function that depends on the void shape, the void volume 114 fraction and the anisotropy of the material, see Eq. (20), σpis a scalar stress measure, see Eq. (16), ¯σeis the 115 generalized eﬀective stress, see Eq. (15), and σyis the yield strength of the matrix material deﬁned as follows: 116 σy=σy0/parenleftbigg 1 +¯εp ε0/parenrightbiggn (13) whereσy0is the initial yield strength, nis the strain hardening exponent and ε0is the reference strain. Moreover, 117 ¯εp=/integraltextt 0˙¯εp(τ)dτis the eﬀective plastic strain of the matrix material and ˙¯εpis the eﬀective plastic strain rate, the 118 latter being the work conjugate of the yield strength of the matrix material. The eﬀective plastic strain rate is 119 related to the plastic rate of deformation tensor assuming that the rate of macroscopic plastic work is equal to 120 the rate of eﬀective plastic work in the matrix material: 121 σ:dp= (1−f)σy˙¯εp(14) The generalized eﬀective stress ¯σeis given by: 122 ¯σe=/radicalBigg σ2e+1 1−h2ζ/parenleftBig h2η2σ2p+ζσ2q+ 2ησpσq/parenrightBig (15) with: 123 σe=/radicalbigg3 2σ:M:σ;σp=1 3σ:X;σq=σ:Q (16) whereQis a deviatoric second-order tensor, Xis a second-order tensor that depends on the shape of the voids 124 through a parameter α2, see Eq. (24), and Mis the anisotropy tensor proposed by Hill (1948). 125 7 Q= 1 20 0 01 20 0 0−1 (17) X= 3 2(1−α2) 0 0 03 2(1−α2) 0 0 0 3 α2 (18) M= G+H−H−G 0 0 0 −H F +H−F 0 0 0 −G−F F +G0 0 0 0 0 0 2 L 0 0 0 0 0 0 2 M 0 0 0 0 0 0 2 N (19) For the sake of simplicity, we consider the axes of the spheroidal voids to be aligned with the orthotropy 126 directions of the material (see Fig. 1 and Section 3), so that the tensor Mis given in terms of six anisotropy 127 coeﬃcients: F,G,H,L,MandN. See in Eq. (19) the representation based on Voigt’s notation. 128 129 The functions κ,ζandη, see Eqs. (12) and (15), are: 130 κ=/radicalBigg 6f(1−f) (1 +g) (f+g) [(β1−fβ2) (h1+ 3h2−4h3) + 6 [α1(1−α1)−fα2(1−α2)]h2+ [1−α1−f(1−α2)] (h1−3h2+ 4h3)] (20) 131 ζ=κ2(1 +g) (f+g) (α2−α1)2 (1−f)2(21) η=κ2(1 +g) (f+g) (α2−α1) 1−f(22) and the void shape parameters g,αiandβi, withi= 1,2, see Eqs. (12), (18) and (20)-(22), are: 132 g=  0 (prolate) h1+3h2+4h3 2h1+h2+2h315π 8(1+α2 1)e3 2√ 1−e2 2(oblate)(23) 8 αi=  1−e2 i e3 iarctanh (ei)−1−e2 i e2 i(prolate) −√ 1−e2 i e3 iarctan/parenleftbigg ei√ 1−e2 i/parenrightbigg +1 e2 i(oblate)(24) βi=  (1−3αi)1 e2 i(prolate) −(1−3αi)1−e2 i e2 i(oblate)(25) Notice that αias deﬁned in equation (24) relates with αGologanu iintroduced by Gologanu et al. (1993, 1994) 133 asαGologanu i =1−αi 2. 134 135 Moreover, the parameters hi, withi= 1,2,3, see Eqs. (15), (20) and (23), are given by: 136 h1=1 4(2H11+ 2H22−H33+ 2H66) ;h2=3 2H33;h3=1 2(H44+H55) where the constants H11,H22,H33,H44,H55, andH66are components of the tensor Hwhich fulﬁlls the condition 137 H:M=K, whereK=I−1 3(1⊗1), withIbeing the unit fourth-order tensor, so that the following expressions 138 are obtained: 139 H11=1 94F+G+H FG+FH+GH;H22=1 9F+ 4G+H FG+FH+GH;H33=1 9F+G+ 4H FG+FH+GH H44=1 L;H55=1 M;H66=1 N 140 Assuming the incompressibility of the matrix material, the evolution of the void volume fraction is given as: 141 ˙f= (1−f)dp:1 (26) Moreover, following Gologanu et al. (1993, 1994) and Keralavarma and Benzerga (2010), the void shape 142 parameterSevolves according to: 143 ˙S=3 2dp/prime 33+1 2/parenleftbigg3α1−1 f+ 1−3α2/parenrightbigg dp:1 (27) wheredp/prime 33is the deviatoric part of the plastic strain rate along the symmetry axis of the voids x3, see Fig. 1, with 144 dp/prime=dp−1 3(dp:1)1. 145 9 Symbol Property and units Value ρ0Initial density/parenleftbigkg/m3/parenrightbig, Eq. (58) 7600 EYoung’s modulus (GPa), Eq. (8) 70 νPoisson’s ratio, Eq. (8) 0.3 σy0Reference yield stress, Eq. (13) 300 nStrain hardening sensitivity, Eq. (13) 0.1 ε0Reference strain, Eq. (13) 0.00429 f0Initial void volume fraction 0.001,0.01,0.1 KAnisotropy coeﬃcient 0.75,1,1.5 S0Initial void shape parameter 1.098(prolate), 0(spherical) and−1.098(oblate) Table 1: Material parameters. Data after Srivastava et al. (2014) and dos Santos et al. (2019b). 2.3. Specialization of the constitutive equations: transversely isotropic matrix material 146 The calculations presented in Section 4 consider the matrix material to be transversely isotropic, so that the 147 anisotropy coeﬃcients of the tensor M, see Eq. (19), fulﬁll the following relations (e.g., see Section 2.1. of dos 148 Santos et al. (2019b) for details): 149 F=G;N=G+ 2H;L=M= 1 (28) leading to: 150 F=G=1 3K2;H=2K2−1 3K2;N=4K2−1 3K2(29) whereK=/radicalBig G+H 2Gis the sole independent anisotropy coeﬃcient. According to Hill (1948), K≥/radicalBig 1 2(see also 151 Section 2.1 in dos Santos et al. (2019b)), being the isotropic material for K= 1. 152 153 The assumption of a transversely isotropic material represents many practical applications (Hill, 1948), and 154 it is particularly convenient from a modeling standpoint because it allows to derive an axisymmetric cylindrical 155 cavitation model when the material is isotropic in the r−θplane, see Section 3. For the speciﬁc formulation 156 developed in Section 3, if K > 1the yield stress of the material is greater parallel to the axis of the cylindrical 157 cavity, and if K < 1the yield stress is greater in the r−θplane. 158 2.4. Material parameters 159 Table 1 shows the values of the initial density and the elastic constants, and the values of the parameters 160 corresponding to the ﬂow strength of the matrix material, which have been taken from Srivastava et al. (2014) 161 and dos Santos et al. (2019b). In Section 4, we carry out a parametric analysis to investigate the eﬀect of plastic 162 anisotropy, initial void volume fraction and initial void shape in the cavitation ﬁelds. In the parametric analysis 163 diﬀerent values of K,f0andS0are used, see Table 1. The values of K,f0are the same considered by dos Santos 164 et al. (2019b), and the values for the initial aspect ratio of the voids are taken from Monchiet et al. (2008). 165 10 166 3. Dynamic cylindrical cavity expansion model 167 The problem is stated as that of a circular cylindrical cavity subjected to internal pressure expanding dynam- 168 ically under plane strain conditions in an inﬁnite porous medium described with the constitutive model presented 169 in Section 2. The novelty lies in the extension of the earlier model of dos Santos et al. (2019b) to consider oblate 170 and prolate voids. As stated in Section 2.3, the matrix material is considered to be transversely isotropic, so that 171 the symmetry axis of both material and voids is parallel to the longitudinal direction of the cylindrical cavity, and 172 therefore the principal directions of stress and strain coincide. Moreover, the artiﬁcial viscosity approach derived 173 by Lew et al. (2001) is used to capture the plastic shock waves that emerge at high cavity expansion velocities, 174 preventing the formation of discontinuities in the elastoplastic ﬁelds that develop near the cavity, and facilitating 175 the numerical integration of the set of diﬀerential equations which govern the problem (dos Santos et al., 2019a,b). 176 177 The cavity expansion model is formulated relative to an Eulerian cylindrical coordinate system (r,θ,z )with 178 origin located at the center of the cavity, see Fig. 2, using notations customarily considered (e.g. see Masri 179 and Durban (2009) and Cohen and Durban (2013), among others). Following the works of Durban and Fleck 180 (1997) and Durban and Masri (2004), the cavity, of instantaneous radius a, is assumed to expand at constant 181 velocity under self-similar, steady-state conditions, so that the only independent variable of the problem is the 182 dimensionless radial coordinate ξ=r/a. Accordingly, the derivative with respect to the radial coordinate becomes: 183 d() dr=d() dξdξ dr= ()/prime1 a(30) where the prime superscript denotes diﬀerentiation with respect to ξ. In addition, the time derivative can be 184 transformed by the following self-similarity relation (Durban and Fleck, 1997; Masri and Durban, 2005): 185 ˙() =d() dξ˙ξ= ()/prime˙a a(v−ξ) (31) wherev=˙r/˙ais the dimensionless radial velocity, with ˙a=constant (constant velocity expansion), so that the 186 dimensionless acceleration reads ˙v=¨r/˙a. 187 188 The principal values of the rate of deformation tensor are: 189 drr=d˙r dr=˙a av/prime(32) dθθ=˙r r=˙a av ξ(33) 11 Uniform porosity distributionz ξξ=1ξpξiξwStrain and stress-free undisturbed zoneElastoplastic zoneElastic zone x1x3p Oξ ξw ξi ξ=1 Oξpx1x2Uniform porosity distribution Figure 2: Schematic of self-similar ﬁeld developed during dynamic expansion of a circular cylindrical cavity in an elastoplastic porous medium with spheroidal voids. The internal pressure pis applied at the cavity wall where ξ= 1. The remote ﬁeld, at ξ > ξw, is undisturbed while behind the elastic wave front, located at ξ=ξw, an elastic range develops. The transition between elastic and elastoplastic regions is denoted by ξi, and the location of the plastic shock wave that emerges at high cavity expansion velocities by ξp. 12 dzz= 0 (34) where equations (30)-(31) have been used to obtain the expressions for the radial and circumferential strain rates. 190 191 Moreover, the hypo-elastic law introduced in equation (7) can be rewritten as: 192 de rr=˙σrr E−ν E( ˙σθθ+ ˙σzz) (35) de θθ=˙σθθ E−ν E( ˙σrr+ ˙σzz) (36) de zz=˙σzz E−ν E( ˙σrr+ ˙σθθ) (37) where, in the absence of material spin, the objective derivative in (7) has been replaced by a time derivative, with 193 σrr,σθθandσzzbeing the principal values of the Cauchy stress tensor. Moreover, the ﬂow rule introduced in 194 equation (9), using the ﬂow potential (12), the work conjugacy relation (14), the generalized eﬀective stress (15), 195 and relationships (16), can be rewritten as: 196 dp rr= (1−f)˙¯εpNrr (38) dp θθ= (1−f)˙¯εpNθθ (39) dp zz= (1−f)˙¯εpNzz (40) where 197 Nrr=1 23(G(Σrr−Σzz)+H(Σrr−Σθθ)) Σy+φΣp Σy+χΣq Σy+κ(1−α2) (g+ 1) (f+g) sinh/parenleftBigκΣp Σy/parenrightBig /parenleftBig¯Σe Σy/parenrightBig2+ (g+ 1) (f+g) sinh/parenleftBigκΣp Σy/parenrightBigκΣp Σy(41) Nθθ=1 23(F(Σθθ−Σzz)+H(Σθθ−Σrr)) Σy+φΣp Σy+χΣq Σy+κ(1−α2) (g+ 1) (f+g) sinh/parenleftBigκΣp Σy/parenrightBig /parenleftBig¯Σe Σy/parenrightBig2+ (g+ 1) (f+g) sinh/parenleftBigκΣp Σy/parenrightBigκΣp Σy(42) Nzz=1 23(F(Σzz−Σθθ)+G(Σzz−Σrr)) Σy+ 2ϕΣp Σy+ 2ψΣq Σy+ 2α2κ(g+ 1) (f+g) sinh/parenleftBigκΣp Σy/parenrightBig /parenleftBig¯Σe Σy/parenrightBig2+ (g+ 1) (f+g) sinh/parenleftBigκΣp Σy/parenrightBigκΣp Σy(43) with 198 13 φ=(1−α2)h2η2+η 1−h2ζ;χ=(1−α2)η+ζ 1−h2ζ;ϕ=α2h2η2−η 1−h2ζ;ψ=α2η−ζ 1−h2ζ(44) where Σrr=σrr/E,Σθθ=σθθ/EandΣzz=σzz/Eare the dimensionless radial, circumferential, and axial stresses, 199 respectively, ¯Σe=¯σe/Eis the dimensionless generalized eﬀective stress and Σy=σy/Eis the dimensionless yield 200 strength of the matrix material. Moreover, Σp=σp/EandΣq=σq/Eare dimensionless scalar stress measures. 201 202 The additive decomposition of the rate of deformation tensor, see equation (6), can be rewritten using expres- 203 sions (32)-(34), (35)-(37) and (38)-(40), as follows: 204 v/prime= (v−ξ)/parenleftBigg Σ/prime rr−ν/parenleftbigΣ/prime θθ+ Σ/prime zz/parenrightbig+ (1−f)∂¯εp ∂ΣyΣ/prime yNrr/parenrightBigg (45) v ξ= (v−ξ)/parenleftBigg Σ/prime θθ−ν/parenleftbigΣ/prime zz+ Σ/prime rr/parenrightbig+ (1−f)∂¯εp ∂ΣyΣ/prime yNθθ/parenrightBigg (46) 0 = (v−ξ)/parenleftBigg Σ/prime zz−ν/parenleftbigΣ/prime rr+ Σ/prime θθ/parenrightbig+ (1−f)∂¯εp ∂ΣyΣ/prime yNzz/parenrightBigg (47) where the relation 3Σh= Σrr+ Σθθ+ Σzzand the time derivative (31) have been used. The derivative ∂¯εp/∂Σy, 205 which is obtained from equation (13), is calculated as: 206 ∂¯εp ∂Σy=ε0 nΣy0/parenleftBigg Σy Σy0/parenrightBigg1 n−1 (48) with Σy0=σy0/Ebeing the dimensionless initial yield strength of the matrix material. 207 208 Moreover, using equations (38)-(40), the expression for the evolution of the void volume fraction given in (26) 209 can be written as follows: 210 f/prime= (1−f)2(¯εp)/prime(Nrr+Nθθ+Nzz) (49) where the time derivative (31) has been used. 211 Similarly, using equations (31) and (38)-(40), the evolution of the void shape parameter given in (27) can be 212 expressed as: 213 S/prime= (1−f) (¯εp)/prime/parenleftbigg (Nzz−Nrr) +1 2/parenleftbigg3α1−1 f+ 1−3α2/parenrightbigg (Nrr+Nθθ+Nzz)/parenrightbigg (50) 214 The balance of mass takes the form: 215 14 ˙ρ+ρtr (d) = 0 (51) which can be rewritten as: 216 (v−ξ)ρ/prime ρ+v/prime+v ξ= 0 (52) using equations (31)-(34). Moreover, inserting expressions (45) and (46) into equation (52) leads to: 217 ρ/prime ρ=−/parenleftbigg 3(1−2ν)Σ/prime h+f/prime (1−f)/parenrightbigg (53) where expression (49) has been used. Following Cohen and Durban (2013), equation (53) is integrated from a 218 given coordinate ξto the elastic wave front ξw, so that the following algebraic relation for the density ratio is 219 obtained: 220 ρ ρ0=/parenleftbigg1−f 1−f0/parenrightbigg exp (−3(1−2ν)Σh) (54) 221 Moreover, the balance of linear momentum along the radial direction takes the form: 222 ∂σt rr ∂r+σt rr−σt θθ r=ρ¨r (55) whereσt rr=σrr+σv rrandσt θθ=σθθ+σv θθ, withσv rrandσv θθbeing the artiﬁcial viscous radial and circumferential 223 stresses deﬁned as (Wilkins, 1980; Lew et al., 2001): 224 σv rr=  c2 0ρ(∆r)2|drr|drr+cLρ∆rVsdrrfordrr<0 0 for drr≥0(56) σv θθ=  c2 0ρ(∆r)2|dθθ|dθθ+cLρ∆rVsdθθfordθθ<0 0 for dθθ≥0(57) wherec0≈2andcL≈1areconstantparameters, ∆risagivenincrementalongtheradialdirectionand Vs=/radicalbig E/ρ 225 (see dos Santos et al. (2019a) for details). 226 227 Using expressions (30) and (31), equation (55) can be rewritten as: 228 /parenleftBig Σt rr/parenrightBig/prime+1 ξ/parenleftBig Σt rr−Σt θθ/parenrightBig =m2ρ ρ0(v−ξ)v/prime(58) 15 where Σt rr= Σrr+ Σv rrandΣt θθ= Σθθ+ Σv θθ, with Σv rr=σv rr/EandΣv θθ=σv θθ/E. Notice that the artiﬁcial viscosity 229 approach is only used in the deﬁnition of the stresses entering into the balance of linear momentum (see dos 230 Santos et al. (2019a,b) for details). Moreover, m=˙a√ E/ρ0is the dimensionless cavitation velocity. Section 4 shows 231 results for 0.05≤m≤0.4, which correspond to dimensional cavitation velocities ranging between 0.1574Km/s 232 and1.2139Km/s for the material parameters given in Table 1. Similar range of cavity expansion velocities was 233 investigated by Cohen and Durban (2013) and dos Santos et al. (2019a). 234 235 The problem is completed by the following boundary conditions: 236 v= 1and Σrr=−Patξ= 1 (59) v= 0and Σrr= Σθθ= 0 atξ=ξw (60) whereP=p Eis the dimensionless applied pressure. An iterative shooting method is employed to enforce the 237 velocity boundary condition, while the stress boundary condition is used to calculate the corresponding applied 238 pressure. The initial conditions, with an unknown constant, and the location of the elastic wave front are deter- 239 mined adapting the elastic solution derived by Masri and Durban (2009) to the constitutive framework used in 240 this paper (see also dos Santos et al. (2019b)). 241 3.1. Plastic shock condition 242 The condition for the plastic shock formation is derived adapting the approach of Cohen and Durban (2013) 243 to the constitutive model presented in Section 2. 244 245 Firstly, the ﬂow potential, Eq. (12), is linearized as follows: 246 Φ/prime≈∂Φ ∂ΣrrΣ/prime rr+∂Φ ∂ΣθθΣ/prime θθ+∂Φ ∂ΣzzΣ/prime zz+/parenleftBigg ∂Φ ∂Σy+∂Φ ∂f(1−f)2(Nrr+Nθθ+Nzz)∂¯εp ∂Σy/parenrightBigg Σ/prime y (61) +∂Φ ∂S(1−f)/parenleftbigg Nzz−1 2(Nrr+Nθθ) +1 2/parenleftbigg3α1−1 f+ 1−3α2/parenrightbigg (Nrr+Nθθ+Nzz)/parenrightbigg∂¯εp ∂ΣyΣ/prime y≈0 where equations (49) and (50) have been used, and the shape functions η,ζ,κandghave been considered to be 247 constant within a small variation of plastic deformation. The derivative∂Φ ∂Sis calculated numerically, while the 248 expressions for∂Φ ∂Σrr,∂Φ ∂Σθθ,∂Φ ∂Σzz,∂Φ ∂Σyand∂Φ ∂fare given by: 249 ∂Φ ∂Σrr=3 (G(Σrr−Σzz) +H(Σrr−Σθθ)) Σ2y+(φΣp+χΣq) Σ2y+κ Σy(1−α2) (1 +g) (f+g) sinh/parenleftBigg κΣp Σy/parenrightBigg (62) 16 ∂Φ ∂Σθθ=3 (H(Σθθ−Σrr) +F(Σθθ−Σzz)) Σ2y+(φΣp+χΣq) Σ2y+κ Σy(1−α2) (1 +g) (f+g) sinh/parenleftBigg κΣp Σy/parenrightBigg (63) ∂Φ ∂Σzz=3 (F(Σzz−Σθθ) +G(Σzz−Σrr)) Σ2y+2 (ϕΣp+ψΣq) Σ2y+ 2α2κ Σy(1 +g) (f+g) sinh/parenleftBigg κΣp Σy/parenrightBigg (64) ∂Φ ∂Σy=−2/parenleftBigg (1 +g) (f+g) sinh/parenleftBigg κΣp Σy/parenrightBigg κΣp Σ2y+¯Σ2 e Σ3y/parenrightBigg (65) ∂Φ ∂f= 2/parenleftBigg (1 +g) cosh/parenleftBigg κΣp Σy/parenrightBigg −(f+g)/parenrightBigg (66) 250 Secondly, the determinant of the matrix of coeﬃcients of the derivatives with respect to ξin Eqs. (45)-(47), 251 (58), and (61) is imposed to be zero, leading to the following singularity condition: 252 ∆ =m2ρ ρ0(v−ξ)2−(1−f) (1+ν)∂¯εp ∂Σy/parenleftBig (Nθθ+νNzz)∂Φ ∂Σθθ+ (Nzz+νNθθ)∂Φ ∂Σzz/parenrightBig −(1−ν)/parenleftBig ∂Φ ∂Σy+ωf∂Φ ∂f+ωS∂Φ ∂S/parenrightBig (1−f)∂¯εp ∂Σy/parenleftBig ωrr∂Φ ∂Σrr+ωθθ∂Φ ∂Σθθ+ωzz∂Φ ∂Σzz/parenrightBig −(1−2ν) (1 +ν)/parenleftBig ∂Φ ∂Σy+ωf∂Φ ∂f+ωS∂Φ ∂S/parenrightBig= 0 (67) with: 253 ωrr= (1−ν)Nrr+ν(Nθθ+Nzz) (68) ωθθ= (1−ν)Nθθ+ν(Nrr+Nzz) (69) ωzz= (1−ν)Nzz+ν(Nrr+Nθθ) (70) ωf= (1−f)2∂¯εp ∂Σy(Nrr+Nθθ+Nzz) (71) ωS= (1−f)/parenleftbigg Nzz−1 2(Nrr+Nθθ) +1 2/parenleftbigg3α1−1 f+ 1−3α2/parenrightbigg (Nrr+Nθθ+Nzz)/parenrightbigg∂¯εp ∂Σy(72) Hereinafter, ∆will be called the shock wave function, and the dimensionless cavity expansion velocity for 254 which the plastic shock wave emerges will be referred to as critical dimensionless cavitation velocity mc. The 255 17 location of the plastic shock wave corresponds to the radial coordinate for which the condition ∆ = 0is met 256 within the elastoplastic zone 1≤ξ≤ξi, see Fig. 2. 257 4. Results 258 This section of the paper is split into two parts. The eﬀect of initial void shape in the elastoplastic ﬁelds that 259 develop near the cavity is investigated in Section 4.1, and the eﬀect of anisotropy in Section 4.2. 260 4.1. The eﬀect of initial void shape 261 Fig. 3 displays the variation in dimensionless applied pressure Pwith dimensionless cavitation velocity mfor 262 three diﬀerent values of the initial void shape parameter: S0= 1.098(prolate witha1 b1= 3),S0= 0(spherical 263 witha1 b1= 1) andS0=−1.098(oblate witha1 b1= 0.33), see equation (4) and Fig. 1. All the results presented 264 in this section correspond to initial void volume fraction f0= 0.01and anisotropy parameter K= 1(isotropic 265 plasticity). The P−mcurves feature a concave-upward shape, consistent with the cylindrical cavity expansion 266 results presented by Masri and Durban (2009) and dos Santos et al. (2019b) for fully dense materials and for 267 porous materials with spherical voids, respectively. The novelty brought to light in this paper is that the eﬀect of 268 the initial shape of the voids in the P−mcurves is negligible. 269 00.050.10.150.20.25 0 0.1 0.2 0.3 0.4 0.5F D BDimensionless applied pressure, P Dimensionless cavitation velocity, mK=1 (Isotropic material) f0=0.01S0=1.098 S0=0 S0=-1.098 Figure 3: Variation in dimensionless applied pressure Pwith dimensionless cavitation velocity mfor three diﬀerent values of the initial void shape parameter: S0= 1.098(prolate witha1 b1= 3),S0= 0(spherical witha1 b1= 1) andS0=−1.098(oblate witha1 b1= 0.33), see equation (4) and Fig. 1. The initial void volume fraction is f0= 0.01, and the anisotropy parameter is K= 1(isotropic plasticity). 270 In contrast, the initial shape of the voids aﬀects the elastoplastic ﬁelds that develop near the cavity. Fig. 271 4 illustrates the variation in (a) dimensionless radial stress Σrr, (b) dimensionless circumferential stress Σθθ, (c) 272 dimensionless axial stress Σzz, (d) dimensionless generalized eﬀective stress ¯Σeand (e) void volume fraction fwith 273 the dimensionless radial coordinate ξ. The dimensionless cavitation velocity is m= 0.2. As in Fig. 3, calculations 274 18 have been performed for S0= 1.098,S0= 0andS0=−1.098. Notice that, while the cavitation ﬁelds obtained for 275 prolate and initially spherical voids are very similar, signiﬁcant diﬀerences appear for oblate voids at intermediate 276 values of the dimensionless radial coordinate. 277 Fig. 4(a) shows that the dimensionless radial stress decreases monotonically with the radial coordinate. At 278 ξ= 1, the absolute value of the radial stress matches with the cavity pressure. The elastic/elastoplastic interface 279 is located at ξi≈4.4. The location of the elastic/elastoplastic interface is taken as the smaller radial coordinate 280 for which the plastic strains are zero. The calculations performed with initially spherical and prolate voids show a 281 sudden drop in the radial stress for an intermediate value of the radial coordinate which indicates the emergence 282 of a plastic shock wave at ξp≈3.75. Recall from Section 3.1 that the location of the shock wave corresponds 283 to the radial coordinate for which the condition ∆ = 0is met within the elastoplastic zone 1≤ξ≤ξi. Note 284 that the change in the radial stress induced by the shock, while abrupt, is not a discontinuity. The artiﬁcial 285 viscosity included in the model provides the shock wave with a ﬁnite width (dos Santos et al., 2019b). On the 286 other hand, the calculations performed with oblate voids predict a more gentle decrease of the radial stress, with 287 no indications of the formation of plastic shocks (the shock wave function does not meet the condition ∆ = 0 288 within the elastoplastic zone). 289 Figs. 4(b) and 4(c) show the circumferential and axial stress ﬁelds. The Σθθ−ξandΣzz−ξcurves are 290 qualitatively and quantitatively similar to the radial stress ﬁelds (note that Σh≈−0.08at the cavity wall). The 291 calculations performed with initially spherical and prolate voids display a sudden drop in the circumferential and 292 axial stress which reveals the formation of a shock wave, while the rate of decrease of ΣθθandΣzzwithin the 293 elastoplastic range is largely constant for the calculations performed with oblate voids. 294 Fig. 4(d) displays the evolution of the dimensionless generalized eﬀective stress with the radial coordinate. 295 Near the cavity wall, ¯Σeis one order of magnitude less than the radial, circumferential and axial stresses –compare 296 Figs. 4(a), 4(b), 4(c) and 4(d)–, i.e., the hydrostatic stress is one order of magnitude greater than the eﬀective 297 stress. The shock wave formation in the calculations for S0= 1.098andS0= 0leads to a sudden drop of 298 eﬀective stress, which partially bounces back immediately after, leading to a short plateau that extends up to 299 the elastic/elastoplastic interface. In the calculation performed with S0=−1.098, the drop of dimensionless 300 generalized eﬀective stress that occurs at ξ≈2.9does not indicate a shock wave –the condition ∆ = 0is not met 301 and the artiﬁcial viscosity is not required to integrate the set of diﬀerential equations which govern the cavity 302 expansion model–, but the total collapse of the oblate voids (see the discussion below). 303 Fig. 4(e) illustrates the evolution of the void volume fraction with the radial coordinate. For the initially 304 prolate and spherical voids, the void volume fraction is zero between the cavity wall and the plastic shock location 305 1≤ξ≤ξp. Within the shock, the value of fsuddenly increases, reaching the initial void volume fraction 306 f0= 0.01at the elastic/elastoplastic interface. The model predicts that for S0= 1.098andS0= 0the voids 307 become increasingly prolate during the compaction process and eventually collapse with Sgrowing unbounded 308 19 (the evolution of Swithξis not shown for the sake of brevity). The calculation corresponding to S0=−1.098 309 reveals that the total collapse of the initially oblate voids occurs at ξ≈2.9, i.e., at shorter distance from the 310 cavity wall. In addition, the closure of the initially oblate voids is not abrupt but gradual and continuous until the 311 elastic/elastoplastic interface is reached. The model predicts that for S0=−1.098the voids become increasingly 312 oblate during the closure process and eventually collapse with Sdecreasing unbounded. 313 314 The inﬂuence of the initial voids shape in the elastoplastic ﬁelds that develop near the cavity is signiﬁcantly 315 reduced for lower and higher cavitation velocities. For instance, Fig. 5 pictures the dimensionless generalized 316 eﬀectivestress ¯Σeandthevoidvolumefraction fversusthedimensionlessradialcoordinate ξforsmallercavitation 317 velocitym= 0.05. The ¯Σe−ξcurves obtained for prolate, spherical and oblate voids virtually overlap each other, 318 see Fig. 5(a). Moreover, while the closure of the voids starts slightly nearer the cavity for oblate voids, the 319 diﬀerences between the f−ξcurves obtained for the three values of S0considered are small, see Fig. 5(b). 320 Fig. 6 shows the evolution of ¯Σeandfversus the dimensionless radial coordinate ξfor greater cavity expansion 321 velocitym= 0.4. The results for the dimensionless generalized eﬀective stress pictured in Fig. 6(a) show that 322 a shock wave is formed at ξp≈2.25for the three initial void shapes considered, including for oblate voids 323 S0=−1.098, for which no plastic shock was observed at m= 0.2. This is an original result of this investigation 324 which indicates that the initial shape of the voids aﬀects the cavity expansion velocity for which plastic shocks 325 ﬁrst emerge. We will further elaborate on this issue in the following paragraphs. The f−ξcurves reported in Fig. 326 6(b) show that the void volume fraction increases rapidly from zero to the initial void volume fraction within the 327 shock wave for prolate, spherical and oblate voids. The diﬀerences in the porosity evolution for the three values 328 ofS0considered are smaller than for lower cavity expansion velocities, see Figs. 4(e) and 5(b), which suggests 329 that the inﬂuence of the initial shape of the voids in the porosity ﬁeld becomes less when plastic shocks emerge. 330 331 Fig. 7 shows the variation in the location of the elastic/elastoplastic interface ξiand the shock wave ξpwith 332 the dimensionless cavitation velocity m. Consistent with the results presented by dos Santos et al. (2019b), both 333 elastic/elastoplastic interface and plastic shock are located closer to the cavity wall as the cavitation velocity 334 increases. The evolution of ξiandξpwithmdisplays a concave-upward shape, so that as the cavity expansion 335 velocity increases, the location of the plastic shock is nearer to the elastic/elastoplastic interface (this result is 336 consistent with the theoretical predictions of Cohen et al. (2010)). Note that the eﬀect of the initial voids shape 337 S0on the location of the elastic/elastoplastic interface is very small. However, as anticipated in the previous 338 paragraph, the minimum cavity expansion velocity at which the shocks form, i.e. the critical cavitation velocity 339 mc, depends on the initial shape of the voids. Namely, Fig. 7 shows that mcis0.255,0.160and0.150, for 340 S0=−1.098,S0= 0andS0= 1.098, respectively. The values of mcare calculated with the condition ∆ = 0, see 341 equation (67), and coincide with the minimum cavity expansion velocity for which the artiﬁcial viscosity approach, 342 20 -0.1 -0.08 -0.06 -0.04 -0.02 0 1 2 3 4 5 6F D BDimensionless radial stress, Σrr Dimensionless radial coordinate, ξm=0.2 K=1 (isotropic material) f0=0.01 shock wave, ξpelastic / elastoplastic interface, ξizoom S0=1.098 S0=0 S0=-1.098-0.04 -0.02 0 3.5 4 4.5 (a) -0.1 -0.08 -0.06 -0.04 -0.02 0 1 2 3 4 5 6F D BDimensionless circumferential stress, Σθθ Dimensionless radial coordinate, ξshock wave, ξpelastic / elastoplastic interface, ξizoom m=0.2K=1 (isotropic material) f0=0.01 S0=1.098 S0=0 S0=-1.098-0.04 -0.02 0 3.5 4 4.5 (b) -0.1 -0.08 -0.06 -0.04 -0.02 0 1 2 3 4 5 6F D BDimensionless axial stress, Σzz Dimensionless radial coordinate, ξm=0.2 shock wave, ξpelastic / elastoplastic interface, ξiK=1 (isotropic material) f0=0.01 S0=1.098 S0=0 S0=-1.098zoom-0.04 -0.02 0 3.5 4 4.5 (c) 00.0020.0040.0060.0080.01 1 2 3 4 5 6F D B Dimensionless radial coordinate, ξDimensionless generalized effective stress, Σe elastic / elastoplastic interface, ξi shock wave, ξpm=0.2K=1 (isotropic material) f0=0.01 S0=1.098 S0=0 S0=-1.098total collapse of oblate voids (d) 00.0020.0040.0060.0080.010.012 1 2 3 4 5 6F D B Dimensionless radial coordinate, ξVoid volume fraction, fK=1 (isotropic material) f0=0.01m=0.2 S0=1.098 S0=0 S0=-1.098shock wave, ξp total collapse of oblate voids (e) Figure 4: Results corresponding to three diﬀerent values of the initial void shape parameter: S0= 1.098(prolate witha1 b1= 3),S0= 0 (spherical witha1 b1= 1) andS0=−1.098(oblate witha1 b1= 0.33), see equation (4) and Fig. 1. The initial void volume fraction is f0= 0.01, and the anisotropy parameter is K= 1(isotropic plasticity). The dimensionless cavitation velocity is m= 0.2. Variation in (a) dimensionless radial stress Σrr, (b) dimensionless circumferential stress Σθθ, (c) dimensionless axial stress Σzz, (d) dimensionless generalized eﬀective stress ¯Σeand (e) void volume fraction fwith dimensionless radial coordinate ξ. 21 00.0020.0040.0060.0080.01 1 3 5 7 9 11 13 15 17 19 21 23 25F D B Dimensionless radial coordinate, ξDimensionless generalized effective stress, Σe K=1 (isotropic material) f0=0.01m=0.05S0=1.098 S0=0 S0=-1.098 (a) 00.0020.0040.0060.0080.010.012 1 3 5 7 9 11 13 15 17 19 21 23 25F D B Dimensionless radial coordinate, ξm=0.05Void volume fraction, fK=1 (isotropic material) f0=0.01 zoomS0=1.098 S0=0 S0=-1.098 00.0020.0040.0060.0080.01 1 2 3 4 (b) Figure 5: Results corresponding to three diﬀerent values of the initial void shape parameter: S0= 1.098(prolate witha1 b1= 3),S0= 0 (spherical witha1 b1= 1) andS0=−1.098(oblate witha1 b1= 0.33), see equation (4) and Fig. 1. The initial void volume fraction is f0= 0.01, and the anisotropy parameter is K= 1(isotropic plasticity). The dimensionless cavitation velocity is m= 0.05. Variation in (a) dimensionless generalized eﬀective stress ¯Σeand (b) void volume fraction fwith dimensionless radial coordinate ξ. 00.0020.0040.0060.0080.01 1 2 3F D B Dimensionless radial coordinate, ξDimensionless generalized effective stress, Σe K=1 (isotropic material) f0=0.01m=0.4 shock wave, ξpelastic / elastoplastic interface, ξi zoomS0=1.098 S0=0 S0=-1.098 0.0020.0030.0040.0050.0060.007 2.2 2.25 2.3 2.35 2.4ξp (a) 00.0020.0040.0060.0080.010.012 1 2 3F D B Dimensionless radial coordinate, ξVoid volume fraction, fK=1 (isotropic material) f0=0.01m=0.4 zoomS0=1.098 S0=0 S0=-1.098 00.0020.0040.0060.0080.010.012 2.2 2.3 2.4 (b) Figure 6: Results corresponding to three diﬀerent values of the initial void shape parameter: S0= 1.098(prolate witha1 b1= 3),S0= 0 (spherical witha1 b1= 1) andS0=−1.098(oblate witha1 b1= 0.33), see equation (4) and Fig. 1. The initial void volume fraction is f0= 0.01, and the anisotropy parameter is K= 1(isotropic plasticity). The dimensionless cavitation velocity is m= 0.4. Variation in (a) dimensionless generalized eﬀective stress ¯Σeand (b) void volume fraction fwith dimensionless radial coordinate ξ. 22 see equations (56)-(57), is required to solve the set of diﬀerential equations which govern the cavity expansion 343 model. 344 1357911 1357911 0 0.1 0.2 0.3 0.4 0.5F D B L J Hξi ξpShock wave location, ξpElastic / elastoplastic interface, ξiK=1 (Isotropic material) f0=0.01 oblate void Dimensionless cavitation velocity, mS0=1.098 S0=0 S0=-1.098 S0=1.098 S0=0 S0=-1.098prolate void spherical void Figure 7: Results corresponding to three diﬀerent values of the initial void shape parameter: S0= 1.098(prolate witha1 b1= 3),S0= 0 (spherical witha1 b1= 1) andS0=−1.098(oblate witha1 b1= 0.33), see equation (4) and Fig. 1. Variation in the location of the elastic/elastoplastic interface ξiand the plastic shock wave ξpwith the dimensionless cavitation velocity m. The initial void volume fraction isf0= 0.01, and the anisotropy parameter is K= 1(isotropic plasticity). 345 The evolution of mcwithS0is further investigated in Fig. 8 for values of the void shape parameter within 346 the range−1.60≤S0≤1.60(which corresponds to 0.2≤a1 b1≤5). The cavity expansion model predicts that for 347 initially oblate voids the relationship between mcandS0is nonlinear, featuring a concave-downward shape with 348 a maximum value of mcforS0≈−1. For initially prolate voids, the value of mcdisplays a mild and roughly 349 linear decrease with the increase of the void shape parameter S0. Note that the case of initially spherical voids 350 S0= 0corresponds to a kink in the mc−S0curve. These results reinforce the idea that the emergence of shock 351 waves in porous materials is modulated by the shape of the voids. Nevertheless, these predictions are speciﬁc for 352 the constitutive model employed in this work, see Section 2, and further research is required to obtain deﬁnite 353 conclusions on the eﬀect of voids shape on shocks formation, see the discussion and perspectives for future work 354 reported in Section 5. 355 4.2. The eﬀect of anisotropy 356 Fig. 9 pictures the variation in dimensionless generalized eﬀective stress ¯Σeand void volume fraction fwith 357 dimensionless radial coordinate ξfor calculations performed with initially prolate voids S0= 1.098. The initial 358 void volume fraction is f0= 0.01and the dimensionless cavitation velocity is m= 0.2. A comparison is performed 359 between calculations corresponding to two diﬀerent values of the anisotropy parameter: K= 0.75andK= 1.5. 360 Recall that if K > 1the yield stress of the material is greater along the axis of the cylindrical cavity, and if K < 1 361 23 0.10.150.20.250.3 -2 -1 0 1 2Critical dimensionless cavitation velocity, mc Initial void shape parameter, S0K=1 (isotropic material) f0=0.01prolate voids oblate voids Figure 8: Variation in the critical dimensionless cavitation velocity mcwith the initial void shape parameter S0, see equation (4) and Fig. 1. The initial void volume fraction is f0= 0.01, and the anisotropy parameter is K= 1(isotropic plasticity). it is greater in the r−θplane. The anisotropy has little eﬀect on the eﬀective stress and porosity ﬁelds. The 362 only noticeable diﬀerence is that the location of the shock wave and the elastic/elastoplastic interface are slightly 363 further from the cavity wall in the case of K= 1.5. A shock is formed at ξp≈3.75for both anisotropy parameter 364 values considered. 365 00.0020.0040.0060.0080.01 1 2 3 4 5 6B D Dimensionless radial coordinate, ξDimensionless generalized effective stress, Σe m=0.2 S0=1.098 f0=0.01K=0.75 K=1.5 elastic / elastoplastic interface, ξi shock wave, ξp zoom 0.0030.0040.0050.006 3.5 4 4.5 (a) 00.0020.0040.0060.0080.010.012 1 2 3 4 5 6B Dm=0.2 S0=1.098 f0=0.01 Dimensionless radial coordinate, ξK=0.75 K=1.5Void volume fraction, fzoom 00.0020.0040.0060.0080.010.012 3.6 3.7 3.8 3.9 4 (b) Figure 9: Results corresponding to prolate voids S0= 1.098(a1 b1= 3) for two diﬀerent values of the anisotropy parameter: K= 0.75 andK= 1.5. The initial void volume fraction is f0= 0.01, and the dimensionless cavitation velocity is m= 0.2. Variation in (a) dimensionless generalized eﬀective stress ¯Σeand (b) void volume fraction fwith dimensionless radial coordinate ξ. 366 These results stand in contrast to the calculations performed for initially oblate voids S0=−1.098for the 367 same cavity expansion velocity m= 0.2, which show a noticeable inﬂuence of the anisotropy parameter on the 368 24 eﬀective stress and porosity ﬁelds, see Fig. 10. Speciﬁcally, the drop in the eﬀective stress for intermediate values 369 of the radial coordinate, which corresponds to the total collapse of the voids, occurs further from the cavity for 370 K= 0.75, see Fig. 10(a). In addition, the closure of the voids with the radial coordinate is more gradual for 371 K= 1.5, see Fig. 10(b), with the elastic/elastoplastic interface being located at a slightly greater value of the 372 radial coordinate. 373 00.0020.0040.0060.0080.01 1 2 3 4 5 6B D Dimensionless radial coordinate, ξDimensionless generalized effective stress, Σe m=0.2 S0=-1.098 f0=0.01 K=0.75 K=1.5elastic / elastoplastic interface, ξi total collapse of oblate voids (a) 00.0020.0040.0060.0080.010.012 1 2 3 4 5 6B D Dimensionless radial coordinate, ξm=0.2 S0=-1.098 f0=0.01Void volume fraction, f K=0.75 K=1.5total collapse of oblate voids (b) Figure 10: Results corresponding to oblate voids S0=−1.098(a1 b1= 0.33) for two diﬀerent values of the anisotropy parameter: K= 0.75andK= 1.5. The initial void volume fraction is f0= 0.01, and the dimensionless cavitation velocity is m= 0.2. Variation in (a) dimensionless generalized eﬀective stress ¯Σeand (b) void volume fraction fwith dimensionless radial coordinate ξ. 374 Nevertheless, the inﬂuence of the parameter Kin the elastoplastic cavitation ﬁelds depends on the cavity 375 expansion velocity, so that the eﬀect of the anisotropy is generally less for other values of m. 376 Results for initially oblate voids S0=−1.098and lower cavity expansion velocity m= 0.05are shown in Fig. 377 11. The ¯Σe−ξcurves obtained for both values of Kconsidered are very similar, see Fig. 11(a), the diﬀerence only 378 appears in the location of the elastic/elastoplastic interface, which is slightly further from the cavity for K= 1.5. 379 Moreover, the f−ξcurves pictured in Fig. 11(b) show that the collapse of the voids occurs slightly closer to the 380 cavity forK= 1.5. 381 Results for initially oblate voids S0=−1.098and greater cavity expansion velocity m= 0.4are shown in 382 Fig. 12. A shock wave is formed at ξp≈2.25for both values of Kconsidered. Small diﬀerences in the ¯Σe−ξ 383 curves appear within the elastoplastic zone at intermediate values of ξso that the eﬀective stress for K= 0.75is 384 slightly below the predictions obtained for K= 1.5, see Fig. 12(a). Moreover, the porosity ﬁelds for K= 0.75 385 andK= 1.5are very similar, with a sudden variation in the void volume fraction within the shock, see Fig. 386 12(b). Consistent with the results presented in Fig. 9(b) for prolate voids, the anisotropy shows little eﬀect on 387 the porosity ﬁelds for cavity expansion velocities leading to shock formation. 388 389 Incontrast, anisotropydoesplayasigniﬁcantroleonthecriticalcavitationvelocity. Fig. 13showsthevariation 390 25 00.0020.0040.0060.0080.01 1 3 5 7 9 11 13 15 17 19 21 23 25B D Dimensionless radial coordinate, ξDimensionless generalized effective stress, Σe S0=-1.098 f0=0.01K=0.75 K=1.5m=0.05 elastic / elastoplastic interface, ξi (a) 00.0020.0040.0060.0080.010.012 1 3 5 7 9 11 13 15 17 19 21 23 25B D Dimensionless radial coordinate, ξVoid volume fraction, fm=0.05 K=0.75 K=1.5 00.0020.0040.0060.0080.01 1 2 3 4zoomS0=-1.098 f0=0.01 (b) Figure 11: Results corresponding to oblate voids S0=−1.098(a1 b1= 0.33) for two diﬀerent values of the anisotropy parameter: K= 0.75andK= 1.5. The initial void volume fraction is f0= 0.01, and the dimensionless cavitation velocity is m= 0.05. Variation in (a) dimensionless generalized eﬀective stress ¯Σeand (b) void volume fraction fwith dimensionless radial coordinate ξ. 00.0020.0040.0060.0080.01 1 2 3B DS0=-1.098 f0=0.01K=0.75 K=1.5m=0.4 Dimensionless radial coordinate, ξDimensionless generalized effective stress, Σe zoomshock wave, ξpelastic / elastoplastic interface, ξi 0.0020.0030.0040.0050.0060.007 2.2 2.25 2.3 2.35 2.4ξp (a) 00.0020.0040.0060.0080.010.012 1 2 3B Dm=0.4 K=0.75 K=1.5 Dimensionless radial coordinate, ξVoid volume fraction, f 00.0020.0040.0060.0080.010.012 2.2 2.3 2.4zoomS0=-1.098 f0=0.01 (b) Figure 12: Results corresponding to oblate voids S0=−1.098(a1 b1= 0.33) for two diﬀerent values of the anisotropy parameter K= 0.75 andK= 1.5. Initial void volume fraction f0= 0.01and dimensionless cavitation velocity m= 0.4. Variation in (a) dimensionless generalized eﬀective stress ¯Σeand (b) void volume fraction fwith dimensionless radial coordinate ξ. 26 inmcwith the initial void shape parameter S0forK= 0.75andK= 1.5. While the mc−S0curves feature 391 similar shape for both anisotropy parameter values considered, there are quantitative diﬀerences in the predictions 392 obtained from the cavity expansion model. Speciﬁcally, for initially oblate voids S0<0, the formation of shocks 393 occurs for greater value of the cavity expansion velocity for K= 1.5, and the opposite trend is obtained for initially 394 prolate voids S0>0, with the absolute value of the slope of the S0>0branch increasing as Kincreases. Note that 395 for initially spherical voids S0= 0, the cavity expansion velocity for which shocks ﬁrst emerge is roughly the same 396 for both values of K. These ﬁndings provide new insights into the interplay between voids shape, anisotropy and 397 shock wave formation in porous materials. However, as mentioned before, the predictions of the cavity expansion 398 model are speciﬁc for the constitutive framework presented in Section 2, so that further research is needed to 399 broaden the generality of our conclusions. 400 0.10.150.20.250.3 -2 -1 0 1 2B D prolate voids oblate voids f0=0.01K=0.75 K=1.5Critical dimensionless cavitation velocity, mc Initial void shape parameter, S0 Figure 13: Results corresponding to two diﬀerent values of the anisotropy parameter: K= 0.75andK= 1.5. Initial void volume fractionf0= 0.01. Variation in the critical dimensionless cavitation velocity mcwith the initial void shape parameter S0, see equation (4) and Fig. 1. 401 The eﬀect of initial porosity on the critical cavitation velocity mcis investigated in Fig. 14 for calculations 402 performed with anisotropy parameter values K= 0.75,K= 1(isotropic plasticity) and K= 1.5, and two diﬀerent 403 initial void volume fractions, f0= 0.001andf0= 0.1in Figs. 14(a) and 14(b), respectively. There is qualitative 404 agreement between the mc−S0curves obtained for both values of f0, and the general trends obtained for the 405 eﬀect ofKon the evolution of the critical cavitation velocity with the initial shape of the voids are virtually 406 independent of the initial void volume fraction (see also the results for f0= 0.01in Fig. 13). On the other hand, 407 the comparison between the results presented for f0= 0.001andf0= 0.1shows that decreasing/increasing the 408 initial porosity delays/promotes the formation of shocks waves (Cohen and Durban, 2013), so that the mc−S0 409 27 curves are shifted upwards/downwards for the whole range of values of S0considered. Note that the eﬀect of initial 410 void volume fraction seems to be particularly important for oblate voids, so that as f0increases, the maximum of 411 themc−S0curves becomes closer to S0= 0(spherical voids). 412 0.20.250.30.350.4 -2 -1 0 1 2B D FCritical dimensionless cavitation velocity, mc Initial void shape parameter, S0f0=0.001prolate voids oblate voids K=0.75 K=1 K=1.5 (a) 0.050.0750.10.130.15 -2 -1 0 1 2B D FCritical dimensionless cavitation velocity, mc Initial void shape parameter, S0f0=0.1K=0.75 K=1 K=1.5prolate voids oblate voids (b) Figure 14: Results corresponding to three diﬀerent values of the anisotropy parameter: K= 0.75,K= 1(isotropic plasticity) and K= 1.5. Variation in the critical dimensionless cavitation velocity mcwith the initial void shape parameter S0for two diﬀerent values of the initial void volume fraction: (a) f0= 0.001and (b)f0= 0.1. 5. Summary and concluding remarks 413 In this work, we have obtained a self-similar solution for the steady-state expansion of a pressurized cylindrical 414 cavity in an inﬁnite porous medium described with the constitutive framework developed by Monchiet et al. 415 (2008). The theoretical model is based on previous developments of Durban and co-workers (Durban and Masri, 416 2004; Masri and Durban, 2005; Cohen et al., 2010; Cohen and Durban, 2013) and includes the artiﬁcial viscosity 417 approach of Lew et al. (2001) to capture the formation of plastic shocks at high cavity expansion velocities. We 418 have carried out a comprehensive parametric analysis which encompasses calculations with diﬀerent initial void 419 shapes−1.6≤S0≤1.6, a wide range of cavity expansion velocities 0.05≤m≤0.4, and diﬀerent anisotropy 420 parameter values varying between 0.75≤K≤1.5. The theoretical model predicts that the relationship between 421 applied pressure and cavity expansion velocity is virtually independent of the shape of the voids included in the 422 porous medium. In contrast, considering that the inﬁnite porous medium contains initially spherical and prolate 423 voids leads to important diﬀerences in the elastoplastic ﬁelds that develop near the cavity when compared to 424 initially oblate voids. Speciﬁcally, spherical and prolate voids favor early formation of shock waves, leading to 425 strong gradients in the radial proﬁles of stress and porosity. In addition, we have derived an analytical expression 426 for the critical cavity expansion velocity leading to shock waves formation, mc, which shows that for initially oblate 427 voids the relationship between mcandS0is nonlinear, featuring a concave-downward shape with a maximum value 428 ofmcfor an intermediate value of S0, while for initially prolate voids the value of mcdecreases roughly linearly 429 28 with the increase of S0. Furthermore, we have shown that the eﬀect of the anisotropy parameter Kon the value 430 ofmcis such that increasing Kdelays the formation of plastic shocks for initially oblate voids, while the opposite 431 trend is generally obtained for initially prolate voids. These results are substantiated for diﬀerent initial void 432 volume fractions within the range 0.001≤f0≤0.1. 433 All the theoretical results presented in this work have been compared with ﬁnite element simulations performed 434 in Abaqus/Explicit using a VUMAT subroutine, in which the model of Monchiet et al. (2008) was implemented 435 following the integration scheme developed by Hosseini and Rodríguez-Martínez (2021). The agreement between 436 numerical simulations and theoretical model predictions for all the cavitation velocities and material features 437 investigated in this work provides additional veriﬁcation of the cavity expansion solution, supporting the results 438 and conclusions presented in this paper (the comparison is not shown for the sake of brevity). However, validation 439 of the theoretical predictions shown herein will be pending until relevant experimental results become available. 440 Note also that the results presented in this paper are speciﬁc for the constitutive framework of Section 2 (which 441 is based on speciﬁc evolution laws for the shape of the voids during deformation), in which the material porosity 442 stands for the single damage parameter with no contribution of the voids size, so that the eﬀect of microinertia 443 –local acceleration ﬁelds that develop near the pores during dynamic compaction– on the mechanical response of 444 the material is not accounted for (Czarnota et al., 2017, 2020). This is a potential limitation of the self-similar 445 cavity expansion model that shall be further assessed using ﬁnite element simulations including the actual porous 446 microstructure of the material (see the recent work of Marvi-Mashhadi et al. (2021)). Such calculations are part of 447 awork in progress that will help to elucidate the speciﬁc role of voids size distribution on the structure of plastic 448 shock waves in porous metallic materials. 449 Funding 450 The research leading to these results has received funding from the European Research Council (ERC) under 451 the European Union’s Horizon 2020 research and innovation programme. Project PURPOSE, grant agreement 452 758056. 453 454 TdS wishes to acknowledge the support of FAPERGS, Fundação de Amparo à Pesquisa do Estado do Rio 455 Grande do Sul, grant agreements 19/2551−0001054−0and21/2551−0000744−3. 456 Conﬂict of interest 457 The authors declare that they have no conﬂict of interest. 458 29 Supplementary material 459 In the following link we provide access to the MATLAB code with the implementation of the dynamic cavity 460 expansion model: https://doi.org/10.21950/AZ6YF7. 461 Author contributions 462 T. dos Santos : Conceptualization; Formal analysis; Funding acquisition; Investigation; Software; Methodol- 463 ogy; Project administration; Resources; Supervision; Validation; Visualization; Writing - original draft; Writing 464 - review & editing. J. C. Nieto-Fuentes : Formal analysis; Investigation; Validation; Writing - original draft; 465 Writing - review & editing. N. Hosseini : Formal analysis; Investigation; Software; Validation; Writing - review 466 & editing. J. A. Rodríguez-Martínez : Conceptualization; Formal analysis; Funding acquisition; Investigation; 467 Methodology; Project administration; Resources; Supervision; Validation; Visualization; Writing - original draft; 468 Writing - review & editing. 469