This paper investigates both theoretically and using ﬁnite elements the elastoplastic ﬁeld induced by a pres- surized spherical cavity expanding dynamically in an inﬁnite medium modelled using the Gurson-Tvergaard- Needleman porous plasticity approach. The theoretical model, which assumes that the porosity is uniformly distributed in the material and the cavitation ﬁelds are self-similar, incorporates artiﬁcial viscous stresses into the original formulation of Cohen and Durban (2013b) to capture the shock waves that emerge at high cavitation velocities. The ﬁnite element calculations, performed in ABAQUS/Explicit (2013) using the Arbitrary Lagrangian Eulerian adaptive meshing available in the code, simulate the cavity expansion process in materials with uniform and non-uniform distributions of porosity. The ﬁnite element results show that the distribution of porosity has small inﬂuence on the cavitation velocity, as well as on the location of the shock wave, which are primarily deter- mined by the cavity pressure and the average material properties. In contrast, it is shown that the intensity of the shock wave, evaluated based on the maximum value of the plastic strain rate within the shock, depends on the local material porosity. The ability of the theoretical model to reproduce the numerical results obtained for the various distributions of porosity used in this work is exposed and discussed. Keywords: Dynamic spherical cavitation, Shock waves, Gurson materials, Porosity distribution, Finite elements 1. Introduction The work of Hopkins (1960) was among the ﬁrst studies to investigate the problem of a spherical cavity expanding dynamically in an elastoplastic medium. A few years later, the model of Hopkins (1960) was used by Goodier (1965) to analyse the penetration of spherical projectiles into metal targets. While in recent years alternative theoretical approaches to study the problem of a projectile penetrating a metallic target have been developed (e.g. Yarin et al. (1995); Roisman et al. (1997); Rubin et al. (2016)), it is still widely accepted that spherical cavitation models are very useful in understanding the mechanics of penetration (Cohen and Durban, ∗Corresponding author. E-mail address: jarmarti@ing.uc3m.es Preprint submitted to Mechanics of Materials April 29, 2019 2 2013a). An important reason for this is that, as stated by Durban and Masri (2004), spherical cavitation models provide fairly simple, yet accurate, expressions for key parameters like the resisting force and penetration depth . Thus, following the work of Goodier (1965), several authors have contributed over the last decades to the research on dynamic spherical cavitation with application to high velocity penetration problems. For instance, Forrestal et al. (1988) used cavity-expansion approximations to obtain closed-form penetration equations for elastic-perfectly plastic targets (Forrestal and Luk, 1988). The authors obtained predictions for penetration depths that were in reasonable agreement with experimental data obtained from tests in which maraging steel rodswithspherical, ogival, andconicalnoseswereimpactedonto6061-T651aluminumtargetsatvelocitiesranging from 0.4to1.4km/s. Shortly after, Forrestal et al. (1991), following the work of Luk et al. (1991), included the strain hardening of the target material in their analytical model. As in their previous work, the theoretical predictions were compared to experimental data obtained from high-velocity impact tests in which projectiles with diﬀerent geometries were used to penetrate 6061-T651 aluminum targets, and satisfactory agreement between theory and experiments was obtained. In Forrestal et al. (1995), the authors compared their analytical models for incompressible and compressible elastic-perfectly plastic constitutive idealizations to numerical simulations performed by Hallquist (1990), McGlaun et al. (1990) and Chen (1995). The predictions of the analytical models, obtained assuming diﬀerent friction coeﬃcients between projectile and target, were in good agreement with the numerical results. Warren and Forrestal (1998) extended previous penetration models published by Forrestal et al. (1988, 1991, 1995) by including the strain rate sensitivity in the constitutive model of the target material. The theoretical results were compared with experiments previously published by the authors for 6061-T651 aluminum targets, and the inﬂuence of material rate sensitivity in the depth of penetration was brought to light. Some years later, Forrestal and Warren (2008) used the penetration model published in Forrestal et al. (1988, 1992) to identify the range of penetration velocities for which target inertia plays a role in the penetration depth of 6061-T651 aluminum targets impacted by ogive-nose projectiles. In this respect, the recent work of Warren (2016) showed that target inertia eﬀects also depend on projectile geometry and nose shape, projectile density and target material properties. However, in the last ﬁfteen years, the researchers who worked most actively in the development of analytical penetration models based on dynamic spherical cavity expansion equations were Durban and co-workers (Durban and Masri, 2004; Masri and Durban, 2005; Cohen et al., 2010; Cohen and Durban, 2013b). Namely, Durban and Masri (2004) studied dynamic spherical cavitation ﬁelds in pressure sensitive elastoplastic materials modelled with Drucker-Prager plasticity. The results presented by the authors revealed that the cavitation pressure increases when the elastic compressibility decreases and the plastic compressibility increases. Masri and Durban (2005) derived a power expansion solution for the problem of a spherical cavity expanding dynamically in a compressible von Mises solid. Note that this analytical study covered only the case where the cavity expansion velocity was 3 much smaller than the speed of sound in the solid. The authors considered both elastic/perfectly plastic response and strain-hardening behaviour, and applied their analytical model to predict the penetration depths obtained in the experiments of Forrestal et al. (1988, 1991). Good agreement between experiments and analytical predictions was reported. Few years later, Cohen et al. (2010) found out that, for cavity expansion velocities greater than those investigated by Masri and Durban (2005), plastic shock waves emerge in the material. The shocks, which led to discontinuities in the ﬁeld variables, were modelled imposing jump (Hugoniot) conditions regarding mass and momentum conservation at the wave front. It was shown that the level of discontinuities induced by the shock increases with the cavitation velocity. The work of Cohen et al. (2010) was continued by Cohen and Durban (2013b) who investigated hypervelocity spherical cavitation ﬁelds in porous elastoplastic solids described with the Gurson model. Cohen and Durban (2013b) showed that material porosity favours the development of shock waves. The authors also showed that in the spherical cavity expansion problem, for a given cavity expansion velocity, the level of discontinuities induced by the shock increases with the material porosity. A key hypothesis in the analytical works of Durban and co-workers (Durban and Masri, 2004; Masri and Durban, 2005; Cohen et al., 2010; Cohen and Durban, 2013b) is the self-similarity of the cavitation ﬁelds. This hypothesis was veriﬁed by Rodríguez-Martínez et al. (2014) with ﬁnite element simulations performed using ABAQUS/Explicit (2013) for elastoplastic materials modelled with von Mises plasticity and diﬀerent strain hard- ening behaviours. The numerical simulations of Rodríguez-Martínez et al. (2014) also showed that shock waves emerge at high cavitation velocities, in agreement with the analytical predictions of Cohen et al. (2010). However, unlike in the theoretical model, the shock wave in the ﬁnite element calculations did not entail a ﬁnite jump in ﬁeld variables. The authors noticed that, due to the artiﬁcial viscosity included by ABAQUS/Explicit (2013), the shock wave in the ﬁnite element calculations was described as a thin band of steep gradients of stress, strain, density, velocity and increased strain rates. Motivated by these results, in the present paper we incorporate artiﬁcial viscous stresses into the theoretical model of Cohen and Durban (2013b) to capture the shock waves that emerge at high cavitation velocities in porous elastoplastic materials modelled using Gurson plasticity. Our approach, which is based on the works of Wilkins (1980) and Lew et al. (2001), leads to smooth solutions of the set of diﬀerential equations which govern the problem, so that it is not necessary to use the Hugoniot conditions for the calculation of the shock. The results obtained from our theoretical model are validated with ﬁnite element simulations performed in ABAQUS/Explicit (2013). A good agreement between theoretical and numerical predictions is obtained in terms of cavitation velocity, stress and density ﬁelds, and shock wave location. The ﬁnite element calculations also show that the steady-state expansion of the cavity is reached in few microseconds (for the highest cavitation velocity considered), thus verifying the assumption of steady cavitation ﬁelds used in the theoretical model. In addition to the ﬁnite element simulations used to validate the theoretical model (which assumes that the initial porosity is uniformly 4 distributed in the material), we have performed numerical calculations with non-uniform distributions of initial porosity that aim at idealizing the micro-structure of metallic materials with inclusions and second-phase particles randomly distributed in the bulk (Ohno and Hutchinson, 1984; Becker, 1987; Huang, 1993; Bilger et al., 2007). These numerical calculations show that the distribution of porosity has small inﬂuence on the cavitation velocity, as well as on the location of the shock wave, which are primarily determined by the average material properties. In contrast, it is shown that the intensity of the shock wave strongly depends on the local porosity at the shock location. 2. Constitutive framework The mechanical behaviour of the material is described using the Gurson-Tvergaard-Needleman (GTN) consti- tutive model (Gurson, 1977; Tvergaard, 1981, 1982) pre-implemented in ABAQUS/Explicit (2013). For the sake of clarity, the main features of the model are brieﬂy presented in this section. The ﬂow potential has the form: Φ=/parenleftBigg σe σy/parenrightBigg2 + 2q1fcosh/parenleftBigg 3q2σh 2σy/parenrightBigg −1−(q1f)2(1) whereq1andq2are material parameters, and fis the void volume fraction. The eﬀective von Mises stress, σe, and the hydrostatic pressure, σh, are deﬁned by: σe=/radicalbigg3 2s:s;σh=1 3σ:1;s=σ−σh1 (2) whereσis the macroscopic Cauchy stress tensor, sis its deviatoric part, and 1is the unit second order tensor. Moreover,σyis the ﬂow strength of the fully dense matrix material described by the following power-type relation: σy=σ0/parenleftbigg 1 +¯εp ε0/parenrightbiggn (3) being ¯εp=/integraltextt 0˙¯εp(τ)dτthe eﬀective plastic strain, where ˙¯εpis the eﬀective plastic strain rate (see equation (8)). Moreover,σ0is the initial yield stress of the matrix material, nis the strain hardening parameter, and ε0is the reference strain. For the sake of simplicity, we have not considered the temperature and strain rate dependences of the ﬂow strength of the matrix material. This, in turn, facilitates to derive a theoretical model for the spherical cavity expansion problem studied in this work, see Appendix A. 5 The rate of deformation tensor is taken to be the sum of an elastic part, de, and a plastic part, dp, as follows: d=de+dp(4) where the elastic part is related to the rate of the stress by the following hypo-elastic law: Oσ=C:de=C: (d−dp) (5) withOσbeing the objective stress rate (it corresponds to the Green-Naghdi derivative in ABAQUS/Explicit (2013)) andCbeing the tensor of isotropic elastic moduli given by: C=E 1 +νI/prime+E 3(1−2ν)1⊗1 (6) whereEis the Young’s modulus, νis the Poisson’s ratio and I/primeis the unit deviatoric fourth order tensor. The plastic part of the rate of deformation tensor follows the direction normal to the ﬂow potential: dp=˙λ∂Φ ∂σ(7) where ˙λis the non-negative plastic ﬂow proportionality factor. The plastic part of the rate of deformation tensor and the eﬀective plastic strain rate in the matrix material are related by enforcing equality between the rates of macroscopic and matrix plastic work: σ:dp= (1−f)σy˙¯εp(8) Moreover, assuming the incompressibility of the matrix material, the evolution of the void volume fraction is deﬁned as: ˙f= (1−f)dp:1 (9) Note that void nucleation is not considered in the present analysis. Hence, if the initial void volume fraction 6 Symbol Property and units Value ρ0 Initial density (kg/m3) 7600 E Young’s modulus (GPa), Eq. (6) 70 ν Poisson’s ratio, Eq. (6) 0.3 q1 Material parameter, Eq. (1) 1.25 q2 Material parameter, Eq. (1) 1.0 σ0 Reference yield stress (MPa), Eq. (3) 300 n Strain hardening exponent, Eq. (3) 0.1 ε0 Reference strain, Eq. (3) 0.00429 f0 Average value of initial VVF 0.001, 0.01, 0.1 Table 1: Material parameters used in the ﬁnite element calculations (Srivastava et al., 2014; N’souglo and Rodríguez-Martínez, 2018). VVF stands for void volume fraction (porosity). is zero, the macroscopic material is fully dense and follows von Mises plasticity. The ﬂow strength of the matrix material, Eq. (3), has been used along with the GTN model pre-implemented in ABAQUS/Explicit through a user-deﬁned subroutine VUHARD (ABAQUS/Explicit, 2013). The material parameters related to the elastic law, Eqs. (5)-(6), the ﬂow potential, Eq. (1), and the ﬂow strength of the matrix material, Eq. (3), taken from Srivastava et al. (2014) and N’souglo and Rodríguez-Martínez (2018), are given in Table 1. Note that the average initial void volume fraction, from now indistinctly referred to as porosity, is the average value of initial porosity over the whole sample, see section 3.1 for details. 3. Finite element model This section describes the features of the axisymmetric ﬁnite element model developed in ABAQUS/Explicit (2013) to simulate the dynamic, spherical, cavity-expansion problem. The model is identical to the one used by Rodríguez-Martínez et al. (2014), that was based on the previous work of Rosenberg and Dekel (2008). The specimen is a sphere of very large radius Rs= 300mm which has a small cavity of radius Rc= 0.5mm in its center. Due to the symmetry of the model, only the θ≥0half of the specimen has been analyzed (see Fig. 1). The sample is initially at rest and undeformed, while a constant internal pressure pis applied at the cavity wall. The dimensions (radii) of the sphere and the cavity do not inﬂuence the simulation results. We have checked that, for all the numerical simulations, the stress waves generated by the application of the cavity pressure reach the free boundary of the sphere after the steady cavitation ﬁeld has developed. The model has been meshed with 240,000four-node axisymmetric reduced integration elements, CAX4R in ABAQUS notation. The mesh shows radial symmetry, in order to retain the symmetry of the problem, with 200 elements along the circumferential direction and 1200along the radial direction. The elements size is constant along the circumferential direction whereas it decreases along the radial direction as the cavity is approached. 7 Namely, the elements show a bias ratio of 100along the radial direction. The dimensions of the elements located near the cavity are 12µm×4µm. Small elements are needed to capture the large gradients of stress and strain which arise in this region of the specimen. Moreover, these large gradients of stress and strain lead to severe mesh distortion. In order to reduce this problem, the Arbitrary Lagrangian Eulerian (ALE) adaptive meshing available in ABAQUS/Explicit (2013) has been applied to the entire model. ALE adaptive meshing uses a single mesh deﬁnition that is gradually smoothed within analysis steps. The frequency of adaptive meshing is set to 1 and the remeshing sweeps per increment are set to 10. As reported by Rodríguez-Martínez et al. (2014), these values of the controlling parameters of the adaptive meshing ensure a proper aspect ratio of the elements during the simulations. Note that ABAQUS/Explicit (2013) introduces artiﬁcial damping in the calculations in order to attenuate the numerical solution and ensure stability, see Aranda-Iglesias et al. (2017). The code generates bulk viscosity pressures, which are linear and quadratic, respectively, in the volumetric strain rate. This artiﬁcial viscosity is controlled by two parameters /pi1(linear viscosity) and χ(quadratic viscosity). Unless otherwise noted, in the calculations shown in section 5, we have used the default values of the code 0.06and1.2, respectively. Only in few calculations, with large initial porosity and applied pressure, we have increased the value of /pi1, see section 5. Figure 1: Axisymmetric ﬁnite element model used to analyse dynamic spherical cavity expansion. This ﬁgure is adapted from Rodríguez-Martínez et al. (2014). 8 3.1. Porosity distribution The porosity in the GTN model idealizes the void nucleating particles that can be found in many metallic materials. These particles are considered sources of damage and preferential sites for the nucleation of cracks. There are many numerical works in the literature, investigating a wide variety of engineering problems, in which the initial porosity is considered uniformly distributed in the material (e.g. Zhang (1996), Vadillo et al. (2008)). However, for some materials, the void nucleating particles may be non-uniformly distributed in the bulk (e.g. Ohno and Hutchinson (1984), Becker (1987), Huang (1993), Bilger et al. (2007)). Accordingly, in this work we perform numerical simulations with specimens with uniform and non-uniform distributions of initial porosity. Three diﬀerent average values of initial void volume fraction f0are used in the calculations reported in section 5, namely 0.001,0.01and0.1(see Table 1). The specimens with uniform distribution of porosity are modelled assigning the same value of f0to all the nodes of the mesh. The specimens with non-uniform distribution of porosity are modelled using the methodology developed by N’souglo and Rodríguez-Martínez (2018). Namely, the average value of the initial porosity in the samplef0is obtained assigning ﬁnite (positive) values of initial porosity only to some nodes (and zero to the others). The percentage of nodes, over the total number of nodes in the sample, that do have initial porosity is denoted asF. Three values of Fhave been investigated in the simulations: 10%,25%and40%. The speciﬁc nodes with positive initial porosity are selected randomly using the function randperm of MATLAB. Each of these nodes has allotted a random value of initial porosity that was generated using the normal distribution function normrnd of MATLAB. The mean and the standard deviation of this normal distribution are µ=f0/Fandσ=µCv, respectively. The relative standard deviation, Cv, determines the maximum dispersion of the generated random values of porosity. Three diﬀerent values of Cvhave been used in the numerical simulations presented in section 5, namely 5%,15%and25%. Moreover, note that, for any given combination of F,Cvandf0, we could generate an inﬁnite number of random porosity distributions: diﬀerent nodes with positive initial porosity and diﬀerent values of porosity for each one. Nevertheless, in the ﬁnite element simulations reported in section 5 we have used only one random non-uniform distribution of porosity. It is out of the scope of this paper to assess the speciﬁc eﬀect that the randomness of the non-uniform porosity distribution has on the ﬁnite element results. Note that, due to the axisymmetry of the ﬁnite element model, the non-uniform porosity distributions are 2D microstructural idealizations (Tvergaard and Needleman, 1992, 2006). Describing more realistic microstructures require a fully 3D analysis, which is out of the scope of this paper. 4. Main features of the theoretical model In this section we summarize the main features of the theoretical model developed to investigate dynamic spherical cavitation ﬁelds in Gurson materials. The complete formulation can be found in Appendix A. The 9 model used in this paper shows two diﬀerences with respect to the formulation reported in the seminal work of Cohen and Durban (2013b): (1) it includes the parameters q1andq2in the ﬂow potential of the constitutive model, and (2) it uses the artiﬁcial viscosity approach to capture the formation of plastic shocks at high cavitation velocities. Note that solutions of rate-independent nonlinear diﬀerential equations commonly exhibit ﬁrst-order discontinuities associated with shock wave propagation. Due to unbounded derivatives, this singular behaviour causes diﬃculties in the numerical integration of the corresponding diﬀerential equations using typical numerical schemes such as the Runge-Kutta methods. While Cohen et al. (2010) and Cohen and Durban (2013b) solved this problem imposing jump (Hugoniot) conditions regarding mass and momentum conservation at the wave front, in this paper (as previously mentioned) we use the artiﬁcial viscosity approach developed, for instance, by von Neumann and Richtmyer (1950); Landshoﬀ (1955); Wilkins (1980) and Lew et al. (2001). This approach, which stands out for its simplicity, leads to smooth solutions of the set of diﬀerential equations which govern the problem, facilitating the integration process. Recall that the outputs of this theoretical model, which considers a uniform distribution of porosity in the specimen, will be compared in section 5 with the results obtained from the ﬁnite element model presented in section 3, in which uniform and non-uniform distributions of porosity are considered. Following Cohen and Durban (2013b), the problem is modelled as a pressurized spherical cavity of instan- taneous radius aexpanding under self-similar, steady-state conditions in an inﬁnite medium, see Fig. 2. Let (r,θ,φ )stand for the Eulerian coordinates of a material point with reference to a spherical system with the origin located at the center of the cavity denoted by Oin Fig. 2. It is assumed that the only independent variable is the nondimensional radial coordinate ξ=r/a. In the formulation used in this paper, see Appendix A, equations (A.8), (A.9), (A.16), (A.27) and (A.29) provide a system of ﬁve diﬀerential equations with derivatives of ﬁve unknowns: dimensionless radial velocity ( v), dimensionless radial elastoplastic stress ( Σrr), dimensionless circum- ferential elastoplastic stress ( Σθθ), dimensionless yield stress ( Σy) and dimensionless radial viscous stress ( Σv rr). Moreover, porosity ( f), density (ρ) and dimensionless circumferential viscous stress ( Σv θθ) are calculated from the algebraic relations (A.15), (A.21) and (A.28). Integration is performed over the independent variable ξ, from the the elastic wave front ξwup to the cavity wallξ= 1, see Fig. 2. The following boundary conditions have been used for the cavity wall: v= 1, Σrr=−P and for the elastic wave front: v= 0, Σrr=Σθθ= 0 10 Initial distribution of porosityξ and stress-free turbed zone = f0 and ρ = ρ0 Elastic zoneElastoplastic zoneξw ξi ξ=1 Oξp Figure 2: Illustration of the steady cavity expansion ﬁeld in a porous medium. The dimensionless internal pressure P=p/Eis applied at the cavity wall where the dimensionless radial coordinate ξ=r/a= 1. The remote region, at ξ > ξw, is a strain and stress-free undisturbed zone with f=f0andρ=ρ0, while behind the elastic wave front, at ξ=ξw, an elastic range develops. The transition between the elastic and elastoplastic regions is denoted by ξi, and the location of the plastic shock wave by ξp(a plastic shock wave emerges at high cavitation velocities). This ﬁgure is adapted from Durban and Masri (2004). whereP=p/Eis the dimensionless applied pressure. The compatibility of the ﬁeld variables at the interface between the elastic zone and the elastoplastic zone ( ξ=ξiin Fig. 2) has been enforced using the closed-form elastic solution developed by Durban and Masri (2004). As in Cohen et al. (2010) and Cohen and Durban (2013b), the shooting method is used to determine the unknown coeﬃcient that arises from the elastic solution. The reader is referred to the seminal works of Durban and Masri (2004), Masri and Durban (2005), Cohen et al. (2010) and Cohen and Durban (2013b) to obtain further details on the derivation and numerical integration of this theoretical model (with the exception of the artiﬁcial viscosity approach to capture the shocks, which is an original contribution of this paper). 5. Results This section of the paper is split into two parts. In subsection 5.1 we compare theoretical and numerical results for the case in which the porosity is uniformly distributed in the ﬁnite element model, which allows to validate the theoretical model in terms of cavitation velocity, stress and density ﬁelds, etc. In subsection 5.2 the comparison is carried out for the case in which the porosity is non-uniformly distributed in the ﬁnite element model. The goal is to assess the capability of the theoretical model to predict cavitation ﬁelds in metallic materials with a microstructure characterized by the random distribution of void nucleating particles (Tvergaard and Needleman, 11 2006; Srivastava et al., 2014, 2017). As mentioned before, such microstructure is idealized here, using a 2D approach, as initial material porosity (N’souglo and Rodríguez-Martínez, 2018). 5.1. Uniform initial distribution of porosity Fig. 3 shows the variation in dimensionless applied pressure Pwith dimensionless cavitation velocity m, see equation (A.25), for three diﬀerent values of the initial porosity f0:0.001,0.01and0.1. The theoretical and ﬁnite element results are compared within a wide range of cavitation velocities 0.04≤m≤0.35. The agreement between both approaches is excellent, demonstrating the potential of the theoretical model to determine the P−m relation. For the three values of f0investigated, the cavitation velocity increases with the applied pressure, the P−mcurve featuring a concave-upward shape, as in Cohen and Durban (2013b). While the diﬀerence between the results obtained for f0= 0.001andf0= 0.01can hardly be noticed, the P−mcurve obtained for f0= 0.1 shows that, for any given applied pressure within the range of values of Pexplored, the cavitation velocity is greater as the initial porosity increases. The ﬁnite element calculations for f0= 0.1andm> 0.2were performed using a linear viscosity parameter /pi1= 0.6greater than the default value of ABAQUS/Explicit (2013). This was necessary so that the code could complete the calculations, despite the large gradients of stress, strain and strain rate that develop in these computations. The numerical simulations predict that for m= 0.04the time required to reach the steady-state expansion of the cavity is ≈30−35µs, whereas for m= 0.35is≈2−3µs, for the three initial porosities investigated. Loading times of the same order of magnitude were reported by Rodríguez- Martínez et al. (2014) for von Mises materials with diﬀerent strain hardening behaviours. Note that for values of mgreater than 0.35, for the largest value of initial porosity considered f0= 0.1, despite of the artiﬁcial damping and the Arbitrary Lagrangian Eulerian adaptive meshing used in the ﬁnite element simulations, excessive element distortion was found in some calculations at the cavity wall and within the shock wave. Therefore, we selected m= 0.35as the upper bound for the cavitation velocities investigated. Fig. 4 shows the dimensionless von Mises ﬂow stress Σe(also referred to as eﬀective von Mises stress in section 2 and Appendix A), see equations (A.11)-(A.12), versus the dimensionless radial coordinate ξfor calculations conducted with three diﬀerent values of initial porosity, f0= 0.001,0.01and0.1, and two diﬀerent cavitation velocities,m= 0.04and0.35. These are the smallest and greatest cavitation velocities investigated. Theoretical and ﬁnite element results are compared. Note that the ﬂow stress results correspond to the elastoplastic material response, i.e. they are the solution of the constitutive problem, see Appendix A. In each of the six plots of Fig. 4 we present ﬁnite element results for three diﬀerent loading times for which the steady-state cavitation has been reached: 40µs,50µs and 60µs if the cavitation velocity is m= 0.04, and 10µs,20µs and 30µs if m= 0.35. In all cases the agreement between theoretical and ﬁnite element results is good, which shows that the theoretical model, in conjunction with the artiﬁcial viscosity approach, captures the stress ﬁelds that develop in the specimen. Note that the ﬁnite element results obtained for diﬀerent loading times virtually lie on top of each 12 00.050.10.150.2 0 0.1 0.2 0.3 0.4T - f=0.001 N - f=0.001 T - f=0.01 N - f=0.01 T - f=0.1 N - f=0.1Dimensionless applied pressure, P Dimensionless cavitation velocity, mTheoretical model - f0=0.001 Numerical model - f0=0.001 Theoretical model - f0=0.01 Numerical model - f0=0.01 Theoretical model - f0=0.1 Numerical model - f0=0.1 Figure 3: Variation in dimensionless applied pressure Pwith dimensionless cavitation velocity mfor three diﬀerent values of the initial porosityf0:0.001,0.01and0.1. The theoretical and ﬁnite element results are compared. The distribution of initial porosity in the ﬁnite element model is uniform. other, conﬁrming the self-similarity of the cavitation ﬁeld, as it is assumed in the theoretical model, see section 4. Meaningful diﬀerences are only found for the greatest porosity and cavitation velocity analysed, as will be discussed later on. Note that there is an elastic wave –with the front located at ξw(elastic precursor)– induced by the application of the velocity/pressure boundary conditions which compresses the material up to the initial elastic limit which is located at ξi. This elastic wave is followed by a plastic wave that compresses further the material (Czarnota et al., 2017). As the cavitation velocity increases, the elastoplastic interface becomes closer to the elastic precursor, see Cohen and Durban (2015). Furthermore, it is worth mentioning that, as it has been discussed by Cohen et al. (2010), no jump (shock) conditions are expected in the primary ﬁeld variables at the elastic precursor. Accordingly, we have not observed an elastic shock in the calculations performed in this work, for the range of cavitation velocities investigated, and the material behavior considered. Next, we carry out speciﬁc analyses of the 6plots presented in Fig. 4. Figs. 4(a)-(b) show results for f0= 0.001–the lowest porosity considered– and the two cavitation velocities investigated. For m= 0.04, which corresponds to Fig. 4(a), the Σe−ξcurve shows a kinkforξ≈4.77that deﬁnes the interface between the elastic zone (on the right) and the elastoplastic zone (on the left). Recall that in section 4 this interface was denoted by ξi. TheΣe−ξcurve, both to the right and to the left of the kink, features a concave-upward shape. As expected, the ﬂow stress decreases as we move away from the cavity wall (which corresponds to ξ= 1). The values of Σeforξ= 1andξiare≈0.0086and≈0.0043, respectively. Note that, according to the asymptotic analysis developed by Durban and Masri (2004), Σeshould be unbounded at the cavity wall. However, this is not the case for the calculations performed in this work. It seems that, due to 13 the discretization and the numerical solution of the problem, the ﬂow stress has a ﬁnite value at the cavity wall. We have checked that the speciﬁc value of Σeforξ= 1depends on the tolerance of the shooting method in the theoretical model, and the size of the elements at the cavity wall in the ﬁnite element simulations. Moreover, for ξ&15the value of Σeis below 10−5. This radial coordinate deﬁnes the elastic wave front. For m= 0.35, which corresponds to Fig. 4(b), the Σe−ξcurve shows a bump, a brief speed-up in the decrease of Σefollowed by a short plateau, which is located just before the elastic / elastoplastic interface ξi≈2.51. Note that the location of the elastic / elastoplastic interface depends on the cavitation velocity (see Fig. 7). The bump, which is zoomed in the upper right part of the plot, indicates the region of the specimen in which the porosity turns from zero forξ <2.19tof0forξ > ξi. Both to the left of the bump and to the right of ξi, theΣe−ξcurve features a concave-upward shape. The value of Σeat the cavity wall is ≈0.0087, similar to the case of m= 0.04. Note also that the value of Σeis below 10−5forξ&3.32. This radial coordinate deﬁnes the elastic wave front. Figs. 4(c)-(d) show results for the initial porosity f0= 0.01, and the cavitation velocities m= 0.04and m= 0.35. Form= 0.04, which corresponds to Fig. 4(c), the Σe−ξcurve is similar to the one obtained for the smaller porosity f0= 0.001. There is a kink, located at the radial coordinate ≈4.67, which corresponds to the elastic / elastoplastic interface ξi. The branches of the Σe−ξcurve on both sides of the kink also feature a concave-upward shape. In addition, Σe≈0.0086forξ= 1, as in the case of f0= 0.001. Form= 0.35, which corresponds to Fig. 4(d), the bump shown for f0= 0.001has turned into a sudden drop of Σe–which partially bounces back immediately after– followed by a short plateau. While the end of the plateau, as mentioned before, corresponds to the elastic / elastoplastic interface ξi= 2.36, the sudden drop of Σeindicates the emergence of a shock wave ξp≈2.20, as discussed by Cohen and Durban (2013b). Both in the theoretical model and the ﬁnite element simulations, the change in the eﬀective stress induced by the shock, while abrupt, is not a discontinuity, i.e. the shock has a ﬁnite width due to the artiﬁcial viscosity used in both approaches, see sections 3, 4 and Appendix A. In fact, the drop of Σeis less abrupt, i.e. the width of the shock wave is greater, in the ﬁnite element simulations than in the theoretical model (see the zoom in the upper right part of the ﬁgure). We have checked that, in all the ﬁnite element simulations, the element size is smaller than the width of the shock, i.e. there are at least two elements within the shock. The fact that the eﬀective stress is non-monotonic in the shock region (it drops and partially bounces back immediately after) occurs due to the porosity of the material and the small strain hardening considered, see Table 1. We have checked that this eﬀect is reduced increasing the strain hardening, and seems to vanish if the initial porosity is zero (fully dense material). The material does not unload in the shock wave, we have checked in the ﬁnite element simulations that the plastic strain in the material increases monotonically with time. Moreover, we have also checked in our calculations that both circumferential and radial stresses decrease monotonically with ξ. Nevertheless, these results are not shown in this paper for the sake of brevity. 14 Figs. 4(e)-(f) show results for m= 0.04andm= 0.35, and the greatest initial porosity investigated f0= 0.1. Form= 0.04, the kink of the Σe−ξcurve, in comparison with the cases f0= 0.001andf0= 0.01, is shifted leftwardsξi= 3.90(see Fig. 7). To the right of the kink, the Σe−ξcurve features a concave-upward shape, as for the smaller initial porosities considered. However, as we move to the left of the kink, the Σe−ξcurve ﬁrst shows a sort of plateau and after a rapid increase in Σe. In fact, near the cavity, the gradients of ﬂow stress are greater than for f0= 0.001andf0= 0.01. Form= 0.35, which corresponds to Fig. 4(f), the abrupt change in Σethat indicates the emergence of a shock is more pronounced than for f0= 0.01, see Fig. 4(d). The ﬂow stress drops from≈0.0065to≈0.0023forξp≈1.70, immediately bounces back until ≈0.0033, and then increases moderately until reaching a peak which corresponds to the elastic / elastoplastic interface ξi≈1.81. Recall that the ﬁnite element calculation reported in Fig. 4(f) has been performed using a value of the linear viscosity parameter/pi1= 0.6greater than the default value of ABAQUS/Explicit (2013), see section 3. It was necessary to increase the viscosity so that the code could complete the calculation. Moreover, note that the values of ξpandξi are smaller as the cavitation velocity and the material porosity increase, as it will be discussed later on. Fig. 5 shows, for the theoretical model, the elastoplastic ﬂow stress Σe(denoted simply as ﬂow stress in Figs. 4 and 9), the viscous ﬂow stress Σvand the total ﬂow stress ¯Σ=Σe+Σv, see equations (A.22) and (A.29), versus the dimensionless radial coordinate ξform= 0.35andf0= 0.01(as Fig. 4(d)). The Σe−ξand ¯Σ−ξcurves are very similar. The ratio Σv/Σeis≈3 %at the cavity wall and decreases with ξuntil the shock wave location is reached. The emergence of the shock leads to a sharp increase of Σvthat leads to a peak of ¯Σ. For values of the radial coordinate greater than ξpthe value of Σvdrops rapidly and the ratio Σv/Σebecomes smaller than 1 %. Fig. 6 shows the dimensionless material density ρ/ρ0, see equation (A.21), versus the dimensionless radial coordinateξfor the same initial porosities and cavitation velocities that were reported in Fig. 4. The results of the theoretical model and the numerical simulations are compared, and satisfactory agreement is found between both approaches. Meaningful diﬀerences are only found for the greatest porosity and cavitation velocity investigated, as it will be discussed later on. For the lowest porosity considered f0= 0.001, Figs. 6(a)-(b), the density decreases smoothly as we move away from the cavity wall, for both cavitation velocities investigated. The theoretical model predicts that, at the cavity wall, ρ/ρ0≈1.02form= 0.04, and≈1.18form= 0.35. Slightly greater values are predicted by the ﬁnite element model. Note that, according to the asymptotic analysis developed by Durban and Masri (2004), the hydrostatic stress should be unbounded at the cavity wall, and consequently the density should drop to zero. However, this is not the case for the calculations performed in this work, most likely, due to the discretization and numerical solution of the problem (as mentioned before). We have checked that the speciﬁc value of the density for ξ= 1depends on the tolerance of the shooting method in the theoretical model, and 15 00.0020.0040.0060.0080.01 1 3 5 7 9 11 13 15B D F HDimensionless flow stress, Σe Dimensionless radial coordinate, ξf0=0.001 m=0.04Theoretical model Numerical model - t=40 μs Numerical model - t=50 μs Numerical model - t=60 μs ξi=4.77 (a) 00.0020.0040.0060.0080.01 1 2 3 4B D F HDimensionless flow stress, Σ e Dimensionless radial coordinate, ξf0=0.001 m=0.35 Theoretical model Numerical model - t=10 μs Numerical model - t=20 μs Numerical model - t=30 μsξi=2.510.00450.0050.0055 2.2 2.3 2.4zoom ξw=3.32 (b) 00.0020.0040.0060.0080.01 1 3 5 7 9 11 13 15B D F HDimensionless flow stress, Σe Dimensionless radial coordinate, ξf0=0.01 m=0.04Theoretical model Numerical model - t=40 μs Numerical model - t=50 μs Numerical model - t=60 μs ξi=4.67 (c) 00.0020.0040.0060.0080.01 1 2 3 4B D F HDimensionless flow stress, Σ e Dimensionless radial coordinate, ξTheoretical model Numerical model - t=10 μs Numerical model - t=20 μs Numerical model - t=30 μszoomf0=0.01 m=0.35 0.0030.0040.0050.006 2.1 2.2 2.3 ξi=2.36ξw=3.32 (d) 00.0020.0040.0060.0080.01 1 3 5 7 9 11 13 15B D F HDimensionless flow stress, Σe Dimensionless radial coordinate, ξf0=0.1 m=0.04 Theoretical model Numerical model - t=40 μs Numerical model - t=50 μs Numerical model - t=60 μsξi=3.900.00370.00380.00390.0040.0041 2 3 4zoom (e) 00.0020.0040.0060.0080.01 1 2 3 4B D F HDimensionless flow stress, Σ e Dimensionless radial coordinate, ξf0=0.1 m=0.35Theoretical model Numerical model - t=10 μs Numerical model - t=20 μs Numerical model - t=30 μs ξi=1.810.0020.0030.0040.0050.0060.007 1.6 1.7 1.8zoom ξw=3.32 (f) Figure 4: Dimensionless ﬂow stress Σeversus dimensionless radial coordinate ξ. Comparison between the theoretical model and ﬁnite element results obtained for diﬀerent loading times. Three diﬀerent values of the initial porosity f0and two diﬀerent cavitation velocitiesmare considered: (a) f0= 0.001andm= 0.04, (b)f0= 0.001andm= 0.35, (c)f0= 0.01andm= 0.04, (d)f0= 0.01 andm= 0.35, (e)f0= 0.1andm= 0.04, (f)f0= 0.1andm= 0.35. The distribution of initial porosity in the ﬁnite element model is uniform. the size of the elements at the cavity wall in the ﬁnite element simulations. For the intermediate initial porosity f0= 0.01, and the cavitation velocity m= 0.04, the ratioρ/ρ0reaches≈1.03forξ= 1, and it decreases smoothly 16 00.0020.0040.0060.0080.01 1 2 3 4Dimensionless elastoplastic, viscous and total flow stress, Σ e, Σv , Σ f0=0.01 m=0.35 Dimensionless radial coordinate, ξTheoretical modelDimensionless viscous flow stress, Σv Dimensionless total flow stress, ΣDimensionless elastoplastic flow stress, Σe Figure 5: Theoretical model results. Elastoplastic ﬂow stress Σe(denoted simply as ﬂow stress in Figs. 4 and 9), viscous ﬂow stress Σvand total ﬂow stress ¯Σ=Σe+Σvversus dimensionless radial coordinate ξ. The cavitation velocity is m= 0.35and the initial porosity isf0= 0.01. as we move away from the cavity wall, see Fig. 6(c). However, for the same initial porosity and higher cavitation velocitym= 0.35, the ratioρ/ρ0shows an abrupt drop which corresponds to the shock wave. Note that the drop ofρ/ρ0in Fig. 6(d) and the drop of Σein Fig. 4(d) occur for the same value of the radial coordinate. Moreover, the zoom in the upper right part of Fig. 6(d) shows that the drop in density while abrupt, is not a discontinuity. As mentioned before, this is due to the artiﬁcial viscosity included in the theoretical model and the ﬁnite element simulations. The results for the greater initial porosity investigated are shown in Figs. 6(e)-(f). For the cavitation velocitym= 0.04the density decreases smoothly with the radial coordinate. However, for m= 0.35, the ratio ρ/ρ0drops from 1.185to1.007for a value of the radial coordinate ≈1.70due to the emergence of a shock wave (i.e.ξp≈1.70), as discussed in Fig. 4(f). Moreover, note that the ﬁnite element results for diﬀerent times do not lie on top of each other. For ξ<ξp, the numerical results for 10µs are above the numerical results for 20µs and 30µs and these, in turn, are slightly above the theoretical results. Recall that in this ﬁnite element simulation we increased the value of the linear viscosity parameter /pi1from 0.06to0.6. Fig. 7 shows the variation in the locations of the elastic / elastoplastic interface ξiand the shock wave ξpwith the dimensionless cavitation velocity mfor two diﬀerent values of the initial porosity, f0= 0.01andf0= 0.1. The theoretical and ﬁnite element results are compared, and very good agreement is found between the results obtained with both approaches. The results for f0= 0.001are not shown because no shocks were found for this level of porosity for the range of cavitation velocities analysed in this paper. In fact, we have checked that for f0= 0.001shocks emerge for values of mgreater than≈0.43. On the other hand, for 0.01and0.1the minimum values ofmfor which shock waves appear are 0.26and0.14, respectively. The location of the shock wave is taken, 17 11.11.21.31.4 1 3 5 7 9 11 13 15B D F HDimensionless material density, ρ/ρ0 Dimensionless radial coordinate, ξTheoretical model Numerical model - t=40 μs Numerical model - t=50 μs Numerical model - t=60 μsf0=0.001 m=0.04 (a) 11.11.21.31.4 1 2 3 4B D F HDimensionless material density, ρ/ρ0 Dimensionless radial coordinate, ξTheoretical model Numerical model - t=10 μs Numerical model - t=20 μs Numerical model - t=30 μs f0=0.001 m=0.35 (b) 11.11.21.31.4 1 3 5 7 9 11 13 15B D F HDimensionless material density, ρ/ρ0 Dimensionless radial coordinate, ξTheoretical model Numerical model - t=40 μs Numerical model - t=50 μs Numerical model - t=60 μsf0=0.01 m=0.04 (c) 11.11.21.31.4 1 2 3 4B D F HDimensionless material density, ρ/ρ0 Dimensionless radial coordinate, ξTheoretical model Numerical model - t=10 μs Numerical model - t=20 μs Numerical model - t=30 μsf0=0.01 m=0.3511.021.041.06 2.1 2.2 2.3zoom (d) 11.11.21.31.4 1 3 5 7 9 11 13 15B D F HDimensionless material density, ρ/ρ0 Dimensionless radial coordinate, ξTheoretical model Numerical model - t=40 μs Numerical model - t=50 μs Numerical model - t=60 μsf0=0.1 m=0.04 (e) 11.11.21.31.4 1 2 3 4B D F HDimensionless material density, ρ/ρ0 Dimensionless radial coordinate, ξTheoretical model Numerical model - t=10 μs Numerical model - t=20 μs Numerical model - t=30 μsf0=0.1 m=0.350.9511.051.11.151.21.25 1.6 1.7 1.8zoom (f) Figure 6: Dimensionless material density ρ/ρ 0versus dimensionless radial coordinate ξ. Comparison between the theoretical model and ﬁnite element results obtained for diﬀerent loading times. Three diﬀerent values of the initial porosity f0and two diﬀerent cavitation velocitiesmare considered: (a) f0= 0.001andm= 0.04, (b)f0= 0.001andm= 0.35, (c)f0= 0.01andm= 0.04, (d)f0= 0.01 andm= 0.35, (e)f0= 0.1andm= 0.04, (f)f0= 0.1andm= 0.35. The distribution of initial porosity in the ﬁnite element model is uniform. 18 both in the theoretical model and the ﬁnite element simulations, as the radial coordinate for which v/prime6−1.4, wherev/prime=dv/dξ, see equation (A.1). On the other hand, the location of the elastic / elastoplastic interface is taken as the minimum radial coordinate for which the plastic strains are zero. Note that ξiandξpdecrease with mfor both values of f0considered, i.e. the elastic / elastoplastic interface and the shock wave are closer to the cavity wall as the cavitation velocity (i.e. the applied pressure) increases. Note also that, for a given initial porosity, the values ofξiandξpget closer to each other as mincreases. This observation agrees with the theoretical results reported by Cohen et al. (2010) who claimed that increasing mwill eventually lead to ξi=ξp. In addition to favouring the emergence of shocks at lower cavitation velocities, the increase in porosity also brings closer to the cavity wall the elastic / elastoplastic interface and the shock wave. Note that the ξi−mandξp−mcurves for f0= 0.1are below the curves for f0= 0.01, being the gap between curves of diﬀerent porosities rather constant for the whole range of cavitation velocities considered. Recall that for the initial porosity f0= 0.1, the ﬁnite element calculations for m> 0.2are performed using the increased linear viscosity parameter /pi1= 0.6. 12345 12345 0 0.1 0.2 0.3 0.4B D F HJ L N PShock wave location, ξp Dimensionless cavitation velocity, mElastic / elastoplastic interface, ξiξi: Theoretical model - f0=0.01 ξi: Numerical model - f0=0.01 ξi: Theoretical model - f0=0.1 ξi: Numerical model - f0=0.1ξp: Theoretical model - f0=0.01 ξp: Numerical model - f0=0.01 ξp: Theoretical model - f0=0.1 ξp: Numerical model - f0=0.1 Figure 7: Variation in the locations of elastic / elastoplastic interface ξiand shock wave ξpwith dimensionless cavitation velocity m for two diﬀerent values of the initial porosity, f0= 0.01andf0= 0.1. The theoretical and ﬁnite element results are compared. The distribution of porosity in the ﬁnite element calculations is uniform. To be noted that the results presented in Figs. 4, 6 and Fig. 7 conﬁrm the theoretical predictions of Cohen and Durban (2013b): shock waves emerge in dynamic spherical cavitation ﬁelds in Gurson materials, and their intensity increases with the cavitation velocity and the initial material porosity. These theoretical predictions are 19 validated in this paper using ﬁnite element simulations. On the other hand, it has to be noted that Cohen and Durban (2015) showed that for planar shocks, except for small porosities ( f0<0.001), the porosity delays the appearance of plastic shock waves. Therefore, it becomes apparent that the inﬂuence of material porosity on the emergence of shock waves depends on the speciﬁc problem addressed (geometry, boundary conditions, stress state, etc.). Further research is still needed to elucidate the speciﬁc mechanisms which control the role of porosity on the emergence of shock waves in elastoplastic porous media. 5.2. Non-uniform initial distribution of porosity Firstly, we check whether the cavitation velocity in specimens which include non-uniform distributions of porosity matches with the one obtained in samples with the porosity homogeneously distributed in the bulk (those presented in previous section). For that purpose, Fig. 8 shows the radial velocity of a node located at the cavity wall ˙aversus the loading time tfor ﬁnite element simulations with specimens in which the porosity is uniformly and non-uniformly distributed in the material. The node taken for the velocity measurements was identiﬁed in Fig. 1. The results in Fig. 8 correspond to the initial porosity f0= 0.01and the cavitation velocity m= 0.35. Recall that for this combination of f0andm, shock waves are formed in the sample, see section 5.1. Recall also that in the simulations with non-uniform distributions of porosity the value of f0is the average porosity in the specimen, see section 3.1. Results are shown for diﬀerent non-uniform distributions of porosity in which the percentage of nodes that do have initial porosity is F= 10%,25%or40%, and the relative standard deviation of the distribution is Cv= 5%,15%or25%. Note that the ˙a−tcurves obtained from the simulations with uniform and non-uniform distributions of porosity are very similar. The steady-state cavity expansion (i.e. ˙a=constant ) starts att≈2−3µs, as mentioned in section 5.1, and the velocity of the cavity is ˙a≈1060m/s, which corresponds to m≈0.35. On the other hand, note that for the simulations which consider a non-uniform distribution of porosity the radial velocity ˙ashows slight ﬂuctuations which are more accentuated as the value of Cvincreases and the value of Fdecreases, i.e. as the porosity distribution becomes more heterogeneous . Note that for the values of FandCvconsidered, the variation of Fhas greater impact on the ﬂuctuations of the ˙a−tcurves than the variation of Cv. On the other hand, note that we have checked that similar results to those presented in Fig. 8 are obtained for other combinations of f0andm. Secondly, we check whether the cavitation ﬁelds in the samples with non-uniform distribution of porosity are self-similar. For that purpose, Fig. 9 shows the dimensionless ﬂow stress Σeversus the dimensionless radial coordinate ξfor calculations conducted with f0= 0.01, and two cavitation velocities, m= 0.04andm= 0.35. The theoretical results are compared with ﬁnite element results corresponding to diﬀerent loading times. In the ﬁnite element calculations the ﬂow stress is measured along the radial path 1 identiﬁed in Fig. 1. The numerical calculationscorrespondtoaspecimenwithnon-uniformdistributionofporosityforwhich Cv= 15%andF= 25%. For the smallest cavitation velocity, Fig. 9(a), the three Σe−ξcurves obtained from the ﬁnite element simulations 20 50060070080090010001100 0 1 2 3 4 5UDP CV5P10 CV5P25 CV5P40Radial velocity of the cavity wall, a (m/s) Loading time, t ( μs)f0=0.01 m=0.35 Finite element simulations Uniform distribution of porosity Non-uniform distribution of porosity: F=10% Non-uniform distribution of porosity: F=25% Non-uniform distribution of porosity: F=40%Cv=5%. (a) 50060070080090010001100 0 1 2 3 4 5B D F HRadial velocity of the cavity wall, a (m/s). Loading time, t ( μs)f0=0.01 m=0.35 Finite element simulationsCv=15% Uniform distribution of porosity Non-uniform distribution of porosity: F=10% Non-uniform distribution of porosity: F=25% Non-uniform distribution of porosity: F=40% (b) 50060070080090010001100 0 1 2 3 4 5B D F HRadial velocity of the cavity wall, a (m/s). Loading time, t ( μs)f0=0.01 m=0.35 Finite element simulationsCv=25% Uniform distribution of porosity Non-uniform distribution of porosity: F=10% Non-uniform distribution of porosity: F=25% Non-uniform distribution of porosity: F=40% (c) Figure 8: Radial velocity of a node located at the cavity wall ˙aversus the loading time tfor ﬁnite element simulations in which the porosity is uniformly and non-uniformly distributed in the material. The initial porosity and the cavitation velocity are f0= 0.01and m= 0.35. Results are shown for diﬀerent non-uniform distributions of porosity in which the percentage of nodes that do have initial porosity isF= 10%,25%or40%, and the relative standard deviation is: (a) Cv= 5%, (b)Cv= 15%and (c)Cv= 25%. are very similar to the theoretical one. There are some slight ﬂuctuations of the ﬂow stress, i.e. deviations from the theoretical solution, only for values of ξsmaller than ξi(the porosity does not aﬀect the elastic range). For this cavitation velocity, the inﬂuence of the non-uniform distribution of porosity on the ﬂow stress ﬁeld is small. On the other hand, for the greatest cavitation velocity, Fig. 9(b), the non-uniformity in the porosity distribution plays a greater role in the numerical results. Note in the zoom shown in the upper right part of the plot that the gradients of ﬂow stress within the shock wave vary with the loading time. Note also the ﬂuctuations of the ﬂow stress for values of ξslightly smaller than ξp. Nevertheless, we have checked that the theoretical predictions ﬁnd satisfactory agreement with the numerical Σe−ξcurves for any combination of f0,CvandFstudied in this paper. The self-similarity of the cavitation ﬁelds, for the specimens with non-uniform distribution of porosity, is further illustrated in Fig. 10 which shows the dimensionless material density ρ/ρ0versus the dimensionless radial coordinate ξfor the same calculations of Fig. 9. The theoretical model is compared with ﬁnite element results 21 obtained for diﬀerent loading times. In the ﬁnite element calculations, the material density is measured along the radial path 1 identiﬁed in Fig. 1. For m= 0.04, Fig. 10(a), the numerical results show oscillations in the density, which are more accentuated as we approach the cavity wall. The zoom in the upper right part of the ﬁgure shows that these oscillations are dependent on the loading time considered. On the other hand, leaving these oscillations aside, the three ﬁnite element curves are similar to each other, and they are also similar to the theoretical curve. Form= 0.35, Fig. 10(b), the analysis is very much alike. On the one hand, there are ﬂuctuations of the density which are more important close to the cavity wall. On the other hand, the three ﬁnite element curves and the analyticalcurvearesimilar, i.e. theaverageslopeofthe ρ/ρ0−ξcurvebeforetheshockwaveissimilar, thelocation of the shock wave is also the same, etc. Comparable agreement between theoretical and numerical predictions for the density ﬁelds has been found for any other combination of f0,CvandFconsidered in this paper. The results presented in Figs. 9 and 10 show that, while there is an inﬂuence of the porosity distribution on the numerical results, the cavitation ﬁelds in specimens with the non-uniform distributions of porosity used in this paper can be considered self-similar. In addition, it is shown that the theoretical model presented in section 4 provides a good estimation of the cavitation ﬁelds that develop in these specimens. 00.0020.0040.0060.0080.01 1 3 5 7 9 11 13 15B D F HDimensionless flow stress, Σe Dimensionless radial coordinate, ξf0=0.01 m=0.04 Theoretical model Numerical model - t=40 μs Numerical model - t=50 μs Numerical model - t=60 μsF=25% Cv=15% ξi=4.67 (a) 00.0020.0040.0060.0080.01 1 2 3 4B D F HTheoretical model Numerical model - t=10 μs Numerical model - t=20 μs Numerical model - t=30 μsf0=0.01 m=0.35F=25% Cv=15% 0.0030.0040.0050.006 2.1 2.2 2.3zoomDimensionless flow stress, Σe Dimensionless radial coordinate, ξξi=2.36 (b) Figure 9: Dimensionless ﬂow stress Σeversus dimensionless radial coordinate ξforradial path 1 (see Fig. 1). Comparison between the theoretical model and ﬁnite element results obtained for diﬀerent loading times. The numerical calculations correspond to a specimen with non-uniform distribution of porosity for which Cv= 15%andF= 25%. The initial material porosity is f0= 0.01. The cavitation velocity is: (a) m= 0.04and (b)m= 0.35. The ﬁnite element contours of material density ρand eﬀective plastic strain rate ˙¯εpin specimens with uniform and non-uniform distributions of porosity are shown in Figs. 11 and 12, respectively. The cavitation velocity is m= 0.35, the loading time t= 20µs, and the initial porosity f0= 0.01. For the specimen with non-uniform distribution of porosity we have that Cv= 15%andF= 25%, as in Figs. 9 and 10. Note that the contours of material density ρand eﬀective plastic strain rate ˙¯εpare shown in a plane that contains the center of the sphere (this plane has been obtained by sweeping from 0◦to180◦the axisymmetric model built in ABAQUS). 22 11.021.041.061.081.1 1 3 5 7 9 11 13 15B D F HDimensionless material density, ρ/ρ0 Dimensionless radial coordinate, ξTheoretical model Numerical model - t=40 μs Numerical model - t=50 μs Numerical model - t=60 μsf0=0.01 m=0.04F=25% Cv=15% 1.011.021.031.041.05 1 1.1 1.2 1.3 1.4 1.5zoom (a) 11.11.21.31.4 1 2 3 4B D F HTheoretical model Numerical model - t=10 μs Numerical model - t=20 μs Numerical model - t=30 μsf0=0.01 m=0.35F=25% Cv=15%Dimensionless material density, ρ/ρ0 Dimensionless radial coordinate, ξ11.021.041.06 2.1 2.2 2.3zoom (b) Figure 10: Dimensionless material density ρ/ρ 0versus dimensionless radial coordinate ξforradial path 1 (see Fig. 1). Comparison between the theoretical model and ﬁnite element results obtained for diﬀerent loading times. The numerical calculations correspond to a specimen with non-uniform distribution of porosity for which Cv= 15%andF= 25%. The initial material porosity is f0= 0.01. The cavitation velocity is: (a) m= 0.04and (b)m= 0.35. The colour coding in Fig. 11 is such that values of material density ranging from 7600kg/m3(the initial value) to9000kg/m3correlate with a colour scale that goes from blue to red. If the value of the material density is above 9000kg/m3, it remains red. No values of ρbelow 7600kg/m3were obtained. The contours of Fig. 11 show that the material density ﬂuctuates along the circumferential direction in the specimen for which the distribution of porosity is non-uniform (unlike what occurs in the sample with uniform initial porosity). In addition, the non-uniformity of the initial porosity distribution makes density not to decrease monotonically with the radial coordinate (as it was also shown in Fig. 10), i.e. the color scale is not traversed unidirectionally from top to bottom as we move away from the cavity wall. For instance, in the specimen with non-uniform distribution of porosity, along any radial path, there are yellow elements (or green, or light blue) closer to the cavity than red elements (or yellow, or green). Note that this is not the case for the specimen with uniform distribution of porosity. Note also that the transition from light blue to dark blue, which is very well deﬁned in the density contours for both specimens, corresponds to the shock wave. The colour coding in Fig. 12 is such that eﬀective plastic strain rates ranging from 103s−1to105s−1correlate with a colour scale that goes from blue to red. Eﬀective plastic strain rates below 103s−1remain blue and above 105s−1remain red. The contour plots show an outer ring, located at some distance from the cavity wall, in which the plastic strain rate reaches its maximum value. This outer ring, whose width varies between ≈300µm and ≈700µm, is the shock wave. Similar shocks widths were reported in the recent work of Czarnota et al. (2017) for metallic porous materials. Depending on the loading time, there may be from 2to more than 10elements inside the shock (recall that the element size varies with the radial coordinate). The evolution of the eﬀective plastic strain rate with the radial coordinate for the ﬁnite element simulations of Fig. 12 is shown in Fig. 13. 23 The eﬀective plastic strain rate is measured along the radial path 1 . The value of ˙¯εpﬁrst decreases as we move away from the cavity wall, then suddenly increases within the shock wave, and the decreases again eventually dropping to zero at the elastic / elastoplastic interface. The maximum eﬀective plastic strain rate, i.e. the shock wave, is located at ξp= 2.20. The distribution of porosity aﬀects the value of the maximum strain rate attained, which is≈5.2 105s−1for the specimen with uniform distribution of porosity, and ≈2.6 105s−1for the specimen with non-uniform distribution of porosity. Similar values (the same order of magnitude) for the maximum strain rate attained in a shock in a porous material were found by Czarnota et al. (2017). While for the speciﬁc radial path selected to obtain the ﬁnite element results presented in Fig. 13 ( radial path 1 ) the inhomogeneous porosity distribution reduces the peak of strain rate, we have checked that for other radial paths, porosity distributions, and loading times, the opposite eﬀect can be obtained. Note also that, outside the shock wave, the strain rate proﬁles for the specimens with uniform and non-uniform distributions of porosity are very similar. The inﬂuence of the porosity distribution on the strain rate ﬁeld is further illustrated in Fig. 14, which shows ﬁnite element results for specimens with two diﬀerent porosity distributions for which the relative standard deviationCvis15%, and the percentage of nodes with non-zero porosity Fis10%in Fig. 14(a), and 25%in Fig. 14(b). The initial material porosity, cavitation velocity and loading time are f0= 0.01,m= 0.35andt= 20µs, respectively. The two plots of Fig. 14 show the evolution of ˙¯εpwithξfor the two radial paths indicated in Fig. 1:radial path 1 andradial path 2 . The strain rate attained in the shock wave is smaller for the radial path 1 . In addition, the diﬀerence between the strain rate attained in the shock wave for the radial paths 1 and2varies withF. This suggests that the intensity of the shock wave –assuming that the intensity of the shock wave can be deﬁned by the maximum strain rate– depends on both the average porosity of the specimen and the local porosity at the location of the shock. 6. Summary and concluding remarks In this paper we have studied dynamic spherical cavitation ﬁelds in porous metals with uniform and non- uniform distributions of porosity. For that purpose we have developed a theoretical model and performed ﬁnite element simulations in ABAQUS/Explicit (2013). The problem addressed is that of a very large sphere, whose mechanical response is described with the GTN constitutive model, with a small cavity in its center subjected to internal pressure. The theoretical model, based on the seminal work of Cohen and Durban (2013b), has the novelty of using the artiﬁcial viscosity approach developed by Wilkins (1980) and Lew et al. (2001) to capture the shock waves that emerge at high cavitation velocities. This is a simple alternative to the methodology used by Cohen et al. (2010) and Cohen and Durban (2013b) who imposed Hugoniot conditions at the shock wave front to obtain the cavitation ﬁelds. Note that the artiﬁcial viscosity approach used in this paper can also be applied to model 24 Density, ρ(kg/m3) Uniform distribution of porosityContours of material density Initial porosity, f0=0.01 -Cavitation velocity, m=0.35 -Loading time, t=20 μs Rs Rc p pShock wave Shock waveNon-uniform distribution of porosity Cv=15%, P=25% p Figure 11: Finite element results. Contours of material density ρfor specimens with uniform and non-uniform distributions of porosity. For the latter case, the values of the distribution parameters are Cv= 15%andF= 25%. The initial material porosity, cavitation velocity and loading time are f0= 0.01,m= 0.35andt= 20µs, respectively. (For interpretation of the