Ductile materials subject to impact loading conditions can undergo spall fracture when an incoming compressive stress wave reflects off interfaces and free surfaces as a tensile stress wave. Experimental observations suggest that in ductile materials, spall fracture is driven by the evolution of porosity. However, the presence of initial porosity in ductile materials also introduces plastic compressibility, which can attenuate the incoming compressive stress wave and, as a result, reduce the amplitude of the reflected tensile stress wave. This, in turn, can mitigate or delay spall fracture. In this work, we report on finite deformation finite element calculations that analyze the response of porous ductile materials subjected to impact loading conditions. Two sets of calculations are carried out, in the first set the material contains initial porosity values ranging from 0% to 5% while in the second set the material also undergoes stress-controlled porosity nucleation. Both sets of calculations are carried out for a wide range of imposed impact velocities. Our results show that porosity in ductile materials can, under certain circumstances, mitigate spall fracture by attenuating stress waves. Results correlating the effects of impact velocity, initial porosity, and porosity nucleation on spall fracture are presented and the underlying mechanisms are discussed. Keywords: Porous Material; Plasticity; Stress Wave; Impact Loading; Spall Fracture; Finite Elements 1. Introduction Materials subjected to high strain rate impact loading conditions may undergo spall fracture, due to the interaction of the incoming compressive stress wave with existing inter- faces/free surfaces that result in the reflection of a tensile stress wave (Antoun et al., 2003). Experimentally, spall fracture is analyzed by subjecting one end of a specimen (referred to as the impact face) to flyer plate impact (Bourne et al., 1995; Jones et al., 2016; Chevrier and Klepaczko, 1999), high explosive charge detonation (Rybakov, 2001; Li et al., 2019), high power particle beams (Baumung et al., 1996) or electrical pulse-generated magnetic fields (Davis, 2020), and measuring the velocity profile of the corresponding opposite free end ∗Corresponding author: Email address: ankit.sri@tamu.edu (Ankit Srivastava) Preprint submitted to MECMAT July 30, 2023 (referred to as the response face), for example, by using a Velocity Interferometer System for Any Reflector (VISAR). The response face velocity profile together with knowledge of the material’s initial density and bulk wave velocity can then be used to quantify the spall strength of the material given the material underwent spall fracture (Antoun et al., 2003; Kanel, 2010). In ductile materials, spall fracture in general involves evolution of porosity (Czarnota et al., 2006; Wright and Ramesh, 2008; Me-Bar et al., 1987; Mallick et al., 2020; Thomason, 1999; Chen et al., 2006; Thissell et al., 2000; Cortes and Elices, 1995), which depends on the presence of initial porosity and/or the nucleation of porosity (Johnson, 1981; Tong and Ravichandran, 1995; Nemat-Nasser and Hori, 1987; Flanagan et al., 2022; Cui et al., 2022; Yao et al., 2018; Moore et al., 2018; Jones et al., 2018; Kraus et al., 2009). The evolution of porosity may also depend on the material microstructure (Clayton, 2005; Me-Bar et al., 1987; Remington et al., 2018; Dongare et al., 2010; Cort´ es, 1992; Hahn et al., 2018; Turley et al., 2018; Hashemian et al., 2008). The presence of initial porosity renders a ductile material plastically compressible and can help attenuate the incoming compressive stress wave during impact loading conditions (Herrmann, 1969). Indeed, the attenuation of the response face velocity has been observed during impact experiments on ceramic and metallic porous foams (Jiang et al., 2021; Pramod et al., 2019; Yao et al., 2018; Bonnan et al., 1998; Jones et al., 2020), and composite materials with a layer of porous foam (Boey, 2009; Gama et al., 2001; Flanagan et al., 2022). The impact experiments on porous foam materials have also been complemented with a number of computational studies. These studies have shown that an initial high level of porosity in a material can significantly attenuate the incoming compressive stress wave and/or the response face velocity under impact loading conditions (Boey, 2009; Jiang et al., 2021; Xu et al., 2011; Erhart et al., 2005; Tang et al., 2018). The attenuation of the response face velocity has also been observed during impact experiments on materials with rather low initial porosity levels ( ≤5%) (Lovinger et al., 2020; Wang et al., 2014; Cui et al., 2022; Lauren¸ con et al., 2020; Wang et al., 2007) that also possess adequate tensile strength required for many structural applications. Furthermore, experimental studies of Cui et al. (2022) on a titanium alloy with initial porosity levels 0 .29%, 0.88% and 5 .41% showed that the specimens with two lower initial porosity levels underwent spall fracture when subjected to impact velocities of 500 m/s and 620 m/s. However, the specimens with 5 .41% initial porosity did not undergo spall fracture when subjected to an impact velocity of 500 m/s, although these specimens did undergo spall fracture at an impact velocity of 620 m/s. This suggests that even a low level of initial porosity in materials may be beneficial in mitigating or at least delaying spall fracture. To explore the effect of porosity on spall fracture, we carry out finite deformation finite element calculations to analyze the response of porous ductile materials with low initial porosity levels under impact loading conditions. Specifically, two sets of calculations are carried out. In the first set, the material contains initial porosity values ranging from 0% to 5% while in the second set in addition to any initial porosity, the material undergoes stress-controlled porosity nucleation. Both sets of calculations are carried out for a wide range of imposed impact velocities. The main aims of our calculations are: (i) to investigate 2 the effects of an initial porosity value on the attenuation of the compressive stress wave propagating through the material for a wide range of impact velocities; and (ii) to explore the extent to which the attenuation of the compressive stress wave can mitigate spall fracture due to the evolution of porosity. 2. Problem formulation and numerical method 2.1. Initial/boundary value problem (a) (b) Figure 1: (a) Sketch of a block subject to plane strain boundary conditions on the external faces normal to theX2andX3axes and with a velocity pulse, vL(t), prescribed on the face X1= 0 (labeled as the loading face). (b) The imposed velocity, vL(t), on the loading face with rise time, tRise, hold time, tHold and decay time, tDecay . Finite deformation finite element calculations are carried out to analyze the response of a three dimensional block subjected to uniaxial straining under an imposed one dimensional velocity pulse, as shown in Fig. 1(a). The calculations are based on a convected coordinate Lagrangian formulation of the field equations with the reference configuration taken to be the initial unstressed state. The position of a material point, relative to a fixed Cartesian frame, in the reference configuration is denoted as Xand in the current configuration is denoted as x. The displacement vector uand the deformation gradient Fare written as u=x−X,F=∂x ∂X(1) The dynamic principle of virtual work accounting for the material inertia is written as Z Vs:δFdV=Z ST·δudS−Z Vρ∂2u ∂t2·δudV (2) 3 sis the (unsymmetric) nominal stress tensor, s=F−1·τwithτ= det ( F)σandσbeing the Cauchy stress tensor, T=s·nwithnbeing the surface normal, ρis the mass density in the reference configuration, tis time, and VandSare, respectively, the volume and the surface of the body in the reference configuration. The displacement boundary conditions on X2= 0, w0andX3= 0, h0, where w0andh0 are the dimensions of the block along X2andX3axes, respectively, are u2(X1,0, X3, t) =u2(X1, w0, X3, t) = 0 (3a) u3(X1, X2,0, t) =u3(X1, X2, h0, t) = 0 (3b) and the velocity boundary condition on X1= 0 is ˙u1(0, X2, X3, t) =vL(t) (4) with (˙) denoting ∂( )/∂t. The time variation of the imposed velocity pulse, shown in Fig. 1(b), is specified as vL(t) =  v0t/trise, fort≤trise v0, fortrise< t≤(trise+thold) v0(trise+thold+tdecay−t)/tdecay,for (trise+thold)< t≤(trise+thold+tdecay) 0, fort >(trise+thold+tdecay) (5) The imposed velocity pulse described in Eq. (5) is referred to as flat-top profile and is typical in flyer plate impact experiments (Dandekar and Weisgerber, 1999; Koller et al., 2005) as well as computations (Clayton, 2005) aimed at analyzing spall fracture. 2.2. Constitutive relation The constitutive framework is the modified Gurson (1977) constitutive relation for a vis- coplastic material matrix that has been used in a variety of previous studies, e.g., Tvergaard and Needleman (1992); Srivastava et al. (2014); Osovski et al. (2015); N’souglo et al. (2018); Liu et al. (2019). For completeness, the basic constitutive equations are listed here. The rate of deformation tensor, d= sym( ˙F·F−1), is written as the sum of an isotropic elastic (actually a hypoelastic) part, de, a thermal part, dΘand a viscoplastic part, dp. The elastic strains are assumed to remain small throughout the deformation history and deis given by de=L−1:ˆσ (6) while dΘis taken to be dΘ=α˙ΘI (7) In Eq. (6), ˆσis the Jaumann rate of Cauchy stress and Lis the tensor of isotropic elastic moduli characterized by the Young’s modulus, E, and Poisson’s ratio, ν. In Eq. (7), αis the thermal expansion coefficient and Θis the temperature. 4 Adiabatic conditions are assumed so that ρcp˙Θ =χτ:dp(8) where cpis the specific heat and χis the Taylor-Quinney factor. The plastic part of the rate of deformation tensor, dp, is calculated based on the modified Gurson (1977) flow potential Φ =σ2 e ¯σ2+ 2q1f∗cosh3q2σh 2¯σ −1−(q1f∗)2(9) Here, q1andq2are parameters introduced by Tvergaard (1981, 1982b), ¯ σis the matrix flow strength, σeis the von Mises effective stress, σhis the hydrostatic stress, and following Tvergaard and Needleman (1984), the function f∗is given as f∗=f f < f c fc+ (1/q1−fc)(f−fc)/(ff−fc)f≥fc(10) In Eq. (10), fis the void volume fraction or simply the porosity, fcis the critical porosity to void coalescence and ffis the porosity at failure. The value of fatt= 0 denotes the initial porosity, f0, in the material and the evolution of the porosity follows ˙f=˙fgrowth +˙fnucl (11) where, presuming the elastic volume change to be negligible, the void growth rate or the evolution of porosity is given by the conservation of volume as ˙fgrowth = (1−f)dp:I (12) The second term in Eq. (11) models stress-controlled porosity nucleation following Chu and Needleman (1980), ˙fnucl=  fN sN√ 2πexp" −1 2 σh+¯σ−σN sN2# ˙¯σ ,if (σh+ ¯σ)≥(σh+ ¯σ)maxand ( ˙σh+˙¯σ)>0 0 otherwise(13) where ( ) maxdenotes the maximum over the deformation history. The Eq. (13) is based on the assumption that there is a mean strength, σN, for porosity nucleation and this nucleation strength is distributed in a normal fashion about the mean with standard deviation, sN. The parameter fNin Eq. (13) determines the maximum value of porosity that can nucleate at a material point. The plastic part of the rate of deformation tensor, dp, is given as (Pan et al., 1983) dp=" (1−f)¯σ˙¯ε σ:∂Φ ∂σ# ∂Φ ∂σ(14) 5 with the matrix plastic strain rate, ˙¯ε, taken as ˙¯ε= ˙ε0¯σ g(¯ε,Θ)1/m , g(¯ε,Θ) =σ0G(Θ) [1 + ¯ ε/ε0]N(15) In Eq. (15), ˙ ε0is the reference strain rate, mis the strain-rate sensitivity exponent, σ0is the reference flow strength, Nis the strain-hardening exponent, ¯ ε=R˙¯εdt,ε0=σ0/E, and G(Θ) is taken as G(Θ) = 1 + bGexp(−cG[Θ0−273]) [exp( −cG[Θ−Θ0])−1] (16) where bGandcGare the parameters that control the temperature-dependence of the flow strength, and Θ0is the reference temperature. 2.3. Constitutive parameters The calculations are carried out for a fixed set of values of the constitutive parameters except for the value of the initial porosity, f0, which is taken to vary in the range 0 ≤ f0≤0.05. The fixed set of values of the other constitutive parameters considered herein do not aim at modeling a specific material but follow those used in previous analyses, e.g., Srivastava et al. (2014, 2017); N’souglo et al. (2018); Liu et al. (2019). Also, the constitutive parameters are taken to be spatially uniform so that there is no material related length scale. In the calculations, Young’s modulus Eis taken to be 70 GPa and Poisson’s ratio νis 0.3. In Eqs. (7) and (8), α= 1×10−5K−1,ρ= 7600 kg/m3,cp= 465 J/(kg K), and χ= 0.9. The flow potential parameters in Eq. (9) are specified as q1= 1.25 and q2= 1.0. The values offc= 0.12 and ff= 0.25 are taken for Eq. (10). In Eq. (13), the mean nucleation strength is specified as σN= 1700 MPa, the standard deviation sN= 0.2σ0, and the magnitude fNis either taken to be 0 (when the stress-controlled porosity nucleation is ignored and only the effects of the initial porosity are analyzed) or 0 .04 (when the effects of both the stress-controlled porosity nucleation and initial porosity are analyzed). The material matrix hardening parameters in Eq. (15) are taken to be ˙ ε0= 103s−1,m= 0.01,N= 0.1, and σ0= 300 MPa. The parameters controlling the temperature-dependence of the flow strength in Eq. (16) are taken to be bG= 0.1406, cG= 0.00793 K−1, and Θ0= 293 K. 2.4. Numerical implementation All the calculations are carried out using an in-house three dimensional data parallel finite element code used in the previous studies, e.g., Srivastava et al. (2014); Osovski et al. (2015); Srivastava et al. (2017); N’souglo et al. (2018); Liu et al. (2019). The finite element code employs twenty node brick elements with eight point Gaussian integration for the force vector and twenty seven point Gaussian integration for the diagonal lumped mass matrix. Explicit Newmark β-method time stepping with β= 0 (Belytschko et al., 1976) is used and the constitutive updating is based on the rate tangent modulus method of Peirce et al. (1984). The material failure is implemented via the element vanishing technique proposed by Tvergaard (1982a). In this technique, when the value of the porosity, f, at an integration 6 point reaches 0 .9ff, the value of fis kept fixed thereafter and the material deforms with a very low flow strength. For the initial/boundary conditions described in Section 2.1, the numerical problem is essentially one dimensional. So that the dimensions w0andh0(along X2andX3axes, respectively) of the block shown in Fig. 1(a) have no effect on the predicted response face velocity and the distribution of the stress, strain and porosity field variables. Following this, the dimensions of the block is taken to be l0= 0.15 mm and w0=h0= 0.002 mm, and discretized using a finite element mesh of 100 ×1×1 elements. Limited calculations are also repeated for a block of dimensions l0=w0= 0.15 mm and h0= 0.002 mm discretized using a finite element mesh of 100 ×100×1 elements for the purpose of visualization and presentation of selected results. All the calculations are carried out for the velocity pulse in Eq. (5) specified by trise= 1 ns, thold= 9 ns and tdecay = 1 ns, with the values of v0in the range 125 m /s≤v0≤300 m /s. 3. Uniaxial strain wave propagation The initial/boundary value problem together with the imposed boundary conditions induce a wave that corresponds to a state of uniaxial strain. Although the formulation here is for finite deformation and allows for a complete loss of stress carrying capacity leading to spall fracture, insight into the response is obtained from a one dimensional small deformation linear elastic analysis. With xdenoting the spatial coordinate, σthe stress and vthe velocity, the governing equation for one dimensional wave propagation is ∂σ ∂x=ρ∂v ∂t(17) and with the uniaxial strain denoted by ε,σis given by σ=E(1−ν) (1 +ν)(1−2ν)ε (18) so that the elastic wave velocity, c, is c=s E(1−ν) ρ(1 +ν)(1−2ν)(19) where for the values of E,νandρ,c= 3521 m/s and for an impact velocity v0the stress carried by the wave is σ=ρcv0 (20) Next, given a body with a geometric length that we take here to be l0(see, Fig. 1), we can define a characteristic time tcand a characteristic length lcas tc=l0 c, l c=l0c v0(21) 7 Figure 2: Variation of the hydrostatic stress, σh, with X1(the direction of imposed velocity pulse) for the boundary conditions and material parameters considered. The material is taken to be fully dense (i.e., the porosity, f= 0 throughout the analysis). The impact velocity is v0= 250 m/s, tis time, tcis the characteristic time and l0is the dimension of the block along the X1axis. Thus, although neither the constitutive relation nor the conventional governing equa- tions contain a length scale, there is a characteristic length scale that emerges from the initial/boundary value problem formulation. Figure 2 shows the evolution of the hydrostatic stress, σh, in a fully dense material (i.e., f= 0 throughout the analysis) subjected to an impact velocity, v0= 250 m/s, and the initial/boundary conditions described in Section 2.1. As can be seen in the figure, post- loading using the velocity pulse given in Eq. (5), a compressive stress pulse forms with a very large hydrostatic stress magnitude because the imposed strain state is nearly uniaxial strain. As this pulse propagates through the material, the magnitude of the induced hydrostatic stress decreases. At t/tc= 0.5, the maximum value of the compressive hydrostatic stress is −5.3 GPa and at t/tc= 1.0, the maximum value of the compressive hydrostatic stress has decreased to −4.1 GPa due to plastic dissipation. The maximum value of the compressive hydrostatic stress in the material at any give time ( t <0.5tc) is found to be −6 GPa which is less than the theoretical value of −6.7 GPa estimated using the Eq. (20), also likely due to plastic dissipation. Next, as t→tca low amplitude elastic precursor stress wave first reaches X1=l0and is followed by the stronger compressive hydrostatic stress pulse traveling at a slower velocity. Finally, when the compressive hydrostatic stress pulse reaches the free surface, it is reflected as a tensile hydrostatic stress pulse. The maximum tensile hydrostatic stress in the reflected stress pulse at t/tc= 1.75 is roughly 3 .2 GPa. For a material with a characteristic time, as for the viscoplastic material here, three non-dimensional ratios emerge r0=v0 c, r 1=(v0/l0) ˙ε0, r 2=ρcv0 σ0(22) where r0provides a ratio of the impact velocity to the characteristic wave velocity, r1provides a ratio of the time scale associated with the loading to the material time scale and r2provides a ratio of the stress carried by the loading wave to the reference flow strength. 8 4. Numerical results Two sets of calculations are carried out using values of the parameters specified in Sec- tion 2.3: (i) values of f0in the range 0 ≤f0≤0.05 (with fN= 0) to analyze the effects of initial porosity; (ii) values of f0in the range 0 ≤f0≤0.05 together with stress-controlled porosity nucleation (i.e., fN̸= 0) as given by Eq. (13) to analyze the effects of both initial porosity and stress-controlled porosity nucleation. Both sets of calculations are carried out for the values of the impact velocity, v0, in the range 125 m /s≤v0≤300 m /s. For the range of the values of v0considered, the non-dimensional ratios in Eq. (22) are in the range 0 .036≤r0≤0.085, 833 .33≤r1≤2000 and 11 .15≤r2≤26.76. These values identify the regime considered in the calculations. The value r0<<1 indicates that the impact velocity is a small fraction of the elastic wave speed. The value r1>>1 indicates that the loading rate, v0/l0, is much larger than the material’s characteristic strain rate, ˙ ε0. The value r2>>1 indicates that the stress magnitude carried by the loading wave is significantly greater than the reference strength σ0of the material. Also, in cases where stress-controlled porosity nucleation is considered, the ratio of the stress magnitude associated with the loading wave to the nucleation stress, ρcv0/σN, is a relevant ratio and that ratio is between 1 .97 and 4 .72. The stress, deformation and porosity fields obtained from the finite element calculations are, except for any round off errors, independent of X2andX3axes. When these field quantities are plotted in the X1−X2plane, the values shown are taken as the average of the values at X3= 0 and X3=h0while when such field quantities are shown as functions ofX1, the values shown are the average of those at ( X1, w0/2,0) and ( X1, w0/2, h0). We note that the loading condition prescribed in the calculations is an imposed velocity for all time, so that vL= 0 for t >(trise+thold+tdecay). This differs from impact loading where the impactor may eventually separate from the target material leading to a traction free surface. 4.1. Effect of initial porosity The time evolution of the response face velocity, vR, for materials with initial porosity values of f0= 0, 0 .005, 0 .015 and 0 .03, subjected to impact velocities, v0= 175 m/s and 250 m/s, are shown in Figs. 3(a) and 3(b), respectively. Fig. 3 shows that for all cases the value of vRstarts to increase as t→tci.e., when the elastic loading wave approaches the response face. However, soon after that the value of vRplateaus for a brief period of time before increasing and reaching a maximum value that can be greater than the value of v0. The second increase in the value of vRafter the initial plateau is due to the trailing plastic loading wave arriving at the response face. For a fixed value of v0, the duration of the initial plateau in the value of vRincreases with increasing f0, while the duration of the maximum value of vRdecreases with increasing f0. For a fixed f0, although the duration of the velocity plateau does not strongly depend on the value of v0, the maximum value of vRsignificantly increases with increasing v0. For v0= 175 m/s, Fig. 3(a), the maximum value of vRis 214 m/s, 167 m/s, 78 m/s and 48 m/s, for f0= 0, 0 .005, 0 .015 and 0 .03, respectively, while for v0= 250 m/s, Fig. 3(b), the maximum value of vRis 330 m/s, 300 m/s, 165 m/s and 62 9 (a) (b) Figure 3: Effect of initial porosity, f0, on the variation of the response face velocity, vR, with normalized time, t/tc. (a) Impact velocity, v0= 175 m/s. (b) Impact velocity, v0= 250 m/s. The time evolution of the imposed velocity pulse, vL, is also plotted for reference. In (a) spall fracture occurred for f0= 0.005 and in (b) it occurred for f0= 0.015. When spall fracture occurs, the pullback velocity, ∆ V, is defined as the difference between the maximum value of vRand the subsequent first minimum value of vRi.e., ∆V=vR(tm)−vR(tp), as shown in (a). m/s, for f0= 0, 0 .005, 0 .015 and 0 .03, respectively. After reaching its maximum, the value ofvReither decays or exhibits ringing. The ringing in the value of vRis a consequence of spall fracture, which leads to wave reverberations between the newly formed surface and the response face. For the cases shown in Fig. 3, with f0= 0.005 spall fracture occurred for v0= 175 m/s while with f0= 0.015 spall fracture occurred for v0= 250 m/s. For calculations with sufficiently large values of f0 the post-maximum decay of vRcan lead to negative values of vR. This is a consequence of prescribing vL= 0 for t > t rise+thold+tdecay. Figure 4 shows a detailed view of the evolution of the hydrostatic stress, σh, leading to spall fracture for a calculation with f0= 0.015 and v0= 250 m/s. The results are shown at six values of t/tc. As shown in Fig. 4, at t= 0 the material is initially stress free. A compressive stress pulse then forms and propagates through the material with gradually decreasing amplitude. The maximum compressive hydrostatic stress at t/tc= 0.5 is−4.8 GPa, while at t/tc= 1.0 it is −3.0 GPa. Also, as a consequence of prescribing vL= 0 for t > t rise+thold+tdecay, the material in the wake of the propagating compressive stress pulse is not stress free. This prescribed loading condition subjects the material to a comparatively low magnitude tensile hydrostatic stress field. After reaching the response face the compressive stress pulse is reflected as a tensile stress pulse. At t/tc= 1.5 the maximum tensile hydrostatic stress in the material is 1 .3 GPa. At 1 .5≤t/tc≤1.75 spall fracture occurs at a location slightly away from the response face. Following this the spalled 10 Figure 4: The spatial distribution of the hydrostatic stress, σh, at six values of the normalized time t/tcfor a calculation in which spall fracture occurred. The impact velocity, v0= 250 m/s and the initial porosity, f0= 0.015. material continues to move rather freely along the positive x1axis as seen in the time frames t/tc= 2.0 and 2 .5. The stress wave reverberation within the spalled piece of material can also be seen in these time frames. The evolution of the hydrostatic stress, σh, as a function of normalized position, X1/l0, and normalized time, t/tcfor calculations with initial porosity values of f0= 0.005, 0 .015 and 0 .03 with impact velocities of v0= 175 m/s and 250 m/s are shown in Fig. 5. This figure shows the formation of a compressive stress pulse that propagates through the material with decreasing magnitude. The hydrostatic stress magnitude associated with this compressive stress pulse decreases with increasing f0, for a fixed value of v0while for a fixed value of f0 the hydrostatic stress magnitude increases with increasing v0. For example, for v0= 175 m/s and at t/tc= 0.5 the compressive hydrostatic stress in the pulse is approximately −3.5 GPa, −3.2 GPa and −2.4 GPa for f0= 0.005, 0 .015 and 0 .03, respectively, while at t/tc= 1.0 it is approximately −2.4 GPa, −1.3 GPa and −0.5 GPa for f0= 0.005, 0 .015 and 0 .03, respectively. Similarly, for v0= 250 m/s and at t/tc= 0.5 the magnitude of the hydrostatic 11 Figure 5: The variation of the hydrostatic stress, σh, with normalized position, X1/l0, and time, t/tc, for f0= 0.005,0.015 and 0 .03, and v0= 175 m/s and 250 m/s. The white regions in the plots for v0= 175 m/s, f0= 0.005 and for v0= 250 m/s, f0= 0.015 mark the spall fracture. stress associated with the stress pulse is approximately −5.1 GPa, −4.8 GPa and −4.2 GPa forf0= 0.005, 0 .015 and 0 .03, respectively, while at t/tc= 1.0 it is approximately −3.8 GPa,−3.0 GPa and −1.3 GPa for f0= 0.005, 0 .015 and 0 .03, respectively. Furthermore, the rate of decrease of the hydrostatic stress magnitude associated with this compressive stress pulse also increases with increasing f0. For example, for v0= 250 m/s the magnitude of the hydrostatic stress at X1/l0= 0.5 (which occurs at t/tc≈0.64, 0 .69 and 0 .76) has decreased by approximately 83%, 68% and 46% for f0= 0.005, 0 .015 and 0 .03, respectively. From the contour plots in Fig. 5 the time it takes for a fixed value of the compressive hydrostatic stress to propagate between two positions in the material can be tracked, from which the plastic loading wave velocity can be estimated. Using the time it takes for a value ofσh=−2 GPa to first reach X1/l0= 0.05 and 0 .25, the plastic loading wave velocity forv0= 175 m/s is estimated to be 0 .76c, 0.7cand 0 .6cforf0= 0.005, 0 .015 and 0 .03, respectively, and for v0= 250 m/s it is estimated to be 0 .78c, 0.75cand 0 .68cforf0= 0.005, 12 0.015 and 0 .03, respectively. This shows that the velocity of the front of the plastic loading wave decreases significantly with increasing f0and increases with increasing v0. The contour plots in Fig. 5 also show the evolution of the comparatively low magnitude tensile hydrostatic stress near the impact face that develops in the wake of the compressive stress pulse. This tensile hydrostatic stress is a consequence of the imposed velocity boundary condition as noted previously. The level of this tensile hydrostatic stress increases with increasing porosity f0because of the associated increase in plastic deformation. Also, the magnitude of this tensile hydrostatic stress increases with increasing v0due to the increase in the overall stress magnitude in the material. As for the compressive hydrostatic stress, the magnitude of the reflected tensile hydro- static stress for a fixed value of v0also decreases with increasing initial porosity f0, and for a fixed value of f0increases with increasing velocity v0. For example, with v0= 175 m/s the maximum value of the reflected tensile hydrostatic stress is found to be about 1 .4 GPa and 0 .67 GPa for f0= 0.005 and f0= 0.03, respectively, while for v0= 250 m/s it is approximately 2 .9 GPa and 0 .74 GPa, for f0= 0.005 and f0= 0.03, respectively. The location of the maximum value of the reflected tensile hydrostatic stress is always slightly away from the response face, as can be seen for v0= 250 m/s and f0= 0.005 in Fig. 5 for which the maximum value is at X1/l0≈0.88. This is due to the destructive interference close to the response face between the incoming compressive hydrostatic stress pulse and the reflected tensile hydrostatic stress pulse. The reflected tensile hydrostatic stress can lead to spall fracture. However, the propensity for spall fracture will depend on the magnitude of the reflected tensile hydrostatic stress and the current level of porosity. Both of which depend on the values of v0andf0. For the cases shown in Fig. 5, spall fracture occurs for v0= 175 m/s and f0= 0.005 at t/tc≈1.5 and X1/l0≈0.9, and for v0= 250 m/s and f0= 0.015 at t/tc≈1.52 and X1/l0≈0.92. In Fig. 5, the white regions in the contour plot mark the spall fracture. Figure 6 shows the evolution of porosity, f, as a function X1/l0andt/tcfor calculations with f0= 0.005, 0 .015 and 0 .03, and with v0= 175 m/s and 250 m/s. Although a value of 0.12 is taken for the critical porosity, fc, in Eq. (10), for all calculations shown, the maximum value of fin the contour plots in Fig. 6 is taken to be 0 .04 for clarity of visualization. As shown in the figure and as expected, the value of fat a material point decreases due to the compressive hydrostatic stress associated with the initial pulse. At the early stages i.e., closer to the loading face and loading time, the decrease in the value of f, regardless of the values of v0andf0is large. In all cases, the values of fat these positions and times have decreased below 10−4(but remain positive). Nevertheless, since for a fixed value of v0the rate of decrease of the magnitude of the compressive hydrostatic stress increases with increasing f0, the distance, X1/l0, over which the value of fdecreases significantly also decreases with increasing f0. While for a fixed value of f0the distance, X1/l0, over which the value of fdecreases significantly increases with increasing impact velocity v0. For example, forv0= 175 m/s at t/tc≈1.3 i.e., just before the compressive hydrostatic stress reflects as a tensile hydrostatic stress, f <10−3forX1/l0<0.87,X1/l0<0.58 and X1/l0<0.34 for the cases with f0= 0.005,f0= 0.015 and f0= 0.03, respectively. Similarly, for v0= 250 m/s, f <10−3forX1/l0<0.98,X1/l0<0.89 and X1/l0<0.6 for the case with f0= 0.005, 13 Figure 6: The variation of porosity, f, with normalized position, X1/l0, and time, t/tc, forf0= 0.005,0.015 and 0 .03, and v0= 175 m/s and 250 m/s. The white regions in the plots for v0= 175 m/s, f0= 0.005 and forv0= 250 m/s, f0= 0.015 mark the spall fracture. f0= 0.015 and f0= 0.03, respectively. Due to the tensile hydrostatic stress in the reflected pulse, the porosity fincreases. However, the rate of increase of the value of f, depends on both the current value of f and on the magnitude of the tensile hydrostatic stress in the reflected stress pulse. Both of these in turn depend on the values of initial porosity f0and impact velocity v0. For example, for f0= 0.005 and v0= 175 m/s, the current value of fnear the response face and the magnitude of the tensile stress in the reflected pulse are both sufficiently large for spall fracture to occur. While for f0= 0.005 and v0= 250 m/s, even though the magnitude of the tensile hydrostatic stress is sufficiently high, the current value of fnear the response face is so low that spall fracture does not occur. On the other hand, for f0= 0.03 even though the current value of fis relatively large, the magnitude of the tensile hydrostatic stress in the reflected stress pulse is so low that spall fracture does not occur for both values of impact 14 velocity. 4.2. Effect of stress-controlled porosity nucleation (a) (b) Figure 7: The response face velocity, vRversus normalized time, t/tc, for calculations with stress-controlled porosity nucleation and for various values of initial porosity, f0. (a) Impact velocity, v0= 175 m/s. (b) Impact velocity, v0= 250 m/s. The time evolution of the imposed velocity pulse, vL, is also plotted for reference. In (a) spall fracture occurred for f0= 0 and in (b) it occurred for both f0= 0 and f0= 0.015. Figure 7 shows the time evolution of the response face velocity, vR, for initial porosity values of f0= 0, 0 .015, 0 .024 and 0 .036, subjected to impact velocities of v0= 175 m/s and 250 m/s for calculations with stress-controlled porosity nucleation, Eq. (13). The time evolution of vRuntil vRreaches its maximum value in Fig. 7 follows the same trend as in Fig. 3 where porosity nucleation is not included. However, the time evolution of vRafter vRreaches the maximum value can differ between Figs. 3 and 7 depending on the values off0andv0. For example, with stress-controlled porosity nucleation for both v0= 175 m/s and 250 m/s, the value of vRexhibits ringing due to spall fracture even for f0= 0. Nevertheless, for several cases the time evolution of vRafter vRreaches the maximum value is also the same in Figs. 3 and 7. For example, for f0= 0.015 and v0= 175 m/s or 250 m/s, the time evolution of vRis nearly the same with or without stress-controlled porosity nucleation. Fig. 7(b) includes a case with f0= 0.036 and v0= 250 m/s, for which the value of vRpost-maximum first decays and then tends to increase. This is attributed to the imposed loading condition leading to damage at the impact face as explained later. The evolution of the hydrostatic stress, σh, as a function of normalized position, X1/l0, and normalized time, t/tcfor calculations with stress-controlled porosity nucleation and initial porosity levels, f0= 0, 0 .015 and 0 .036, and impact velocities, v0= 175 m/s and 250 m/s, are shown in Fig. 8. Many features of the results shown here are also similar to those shown in Fig. 5 where only the effects of the initial porosity in the material is 15 Figure 8: The variation of the hydrostatic stress, σh, with normalized position, X1/l0, and time, t/tc, for calculations with stress-controlled porosity nucleation and with f0= 0,0.015 and 0 .036, and v0= 175 m/s and 250 m/s. The white regions in the plots for v0= 175 m/s, f0= 0, for v0= 250 m/s, f0= 0 and for v0= 250 m/s, f0= 0.015 mark the spall fracture. modeled. However, because in these cases voids can nucleate following the stress-controlled nucleation criterion, Eq. (13), spall fracture occurs as long as the magnitude of the tensile hydrostatic stress in the reflected pulse is sufficiently high. Thus, with stress-controlled porosity nucleation, spall fracture can occur even with no initial porosity ( f0= 0) as shown in Fig. 8. The magnitude of the tensile hydrostatic stress in the reflected pulse for a fixed value of impact velocity v0decreases with increasing f0, and for a fixed value of f0increases with increasing v0. Thus, for v0= 175 m/s spall fracture does not occur for sufficiently large values of f0, for example, for f0= 0.015 and 0 .036 in Fig. 8 but for a fixed value of f0spall can occur for a sufficiently large value of v0as is the case for f0= 0.015 with v0= 250 m/s. Figure 9 shows the evolution of porosity, f, as a function X1/l0andt/tcfor calculations with stress-controlled porosity nucleation and with initial porosity values f0= 0, 0 .015 and 0.036, and impact velocities v0= 175 m/s and 250 m/s.The evolution of porosity before the 16 Figure 9: The variation of porosity, f, with normalized position, X1/l0, and time, t/tc, for calculations with stress-controlled porosity nucleation and with f0= 0,0.015 and 0 .036, and v0= 175 m/s and 250 m/s. The white regions in the plots for v0= 175 m/s, f0= 0, for v0= 250 m/s, f0= 0 and for v0= 250 m/s, f0= 0.015 mark the spall fracture. compressive stress pulse reflects as a tensile stress pulse is the same as for the cases without stress-controlled porosity nucleation. However, as shown in the figure with stress-controlled porosity nucleation, the reflected tensile stress pulse can lead to nucleation of voids that can grow and cause spall fracture even for the cases with f0= 0. Also, because the magnitude of the tensile stress in the reflected pulse increases with increasing impact velocity v0, the width over which stress-controlled porosity nucleation occurs also increases with increasing v0. For all cases with f0>0, the magnitude of the tensile stress in the reflected pulse also depends on the value of f0. So that for v0= 175 m/s, and f0= 0.015 and 0 .036 even though the current value of f >0, the stress magnitude in the reflected tensile pulse is not sufficient to cause spall fracture due to porosity nucleation and growth. While for v0= 250 m/s and f0= 0.015 both the value of fnear the response face and the stress magnitude associated with the reflected pulse are sufficiently large for spall fracture to occur. 17 The tensile hydrostatic stress magnitude near the impact face that develops due to the imposed loading conditions can also lead to the nucleation and growth of porosity. Because the magnitude of the tensile hydrostatic stress increases with increasing values of v0andf0, for sufficiently large values of v0andf0, this can lead to impact face damage as shown in Fig. 9 for the case with v0= 250 m/s and f0= 0.036. This impact face damage is the reason that values of vRmay increase following the post-maxima decay for cases with sufficiently large values of v0andf0as in Fig. 7(b). 4.3. Spall fracture and spall strength (a) (b) Figure 10: The dependence of spall fracture (marked by ‘ ×’ symbols) or no spall fracture (marked by ‘o’ symbols) on the initial porosity, f0, and the impact velocity, v0. (a) For calculations with only an initial porosity. (b) For calculations with an initial porosity together with stress-controlled porosity nucleation. The dependence of spall fracture on initial porosity, f0, and impact velocity, v0, is shown in Fig. 10. Fig. 10(a) shows results for calculations with only initial porosity and Fig. 10(b) shows results for initial porosity together with stress-controlled porosity nucleation. Cal- culations for which spall fracture occurred are marked with a red ‘ ×’ and calculations for which spall fracture did not occur are marked with a green ‘o’. Figure 10(a) shows that with only initial porosity, spall fracture does not occur for v0= 125 m/s for the range of values of f0considered. For v0≥150 m/s there is a minimum value of f0needed to trigger spall fracture and this minimum value of f0increases approximately linearly with increasing v0. For v0≥150 m/s there is also a maximum value off0above which spall fracture does not occur. This maximum value of f0also increases roughly linearly with increasing v0. The rate of increase of the maximum value of f0with v0 above which spall does not occur is greater than the rate of increase of the minimum value off0with v0below which spall does not occur. Hence, the range of the values of f0over which spall does occur increases with increasing v0. 18 In Fig. 10(b), for calculations with both initial porosity and stress-controlled porosity nucleation, there is no minimum value of f0below which spall does not occur. Nevertheless, there is a maximum value of f0above which spall fracture does not occur, which also increases nearly linearly with increasing v0. The maximum value of f0above which spall does not occur for a fixed value of v0is approximately the same as for the calculations without stress-controlled porosity nucleation. Figure 11: The variation of the spall strength, σsp, with initial porosity, f0, for three values of the impact velocity, v0, for calculation in which spall fracture occurred. The solid lines and open symbols correspond to calculations with only an initial porosity. The dashed lines and closed symbols correspond to calculations with an initial porosity together with stress-controlled porosity nucleation. Figure 11 shows the variation of the spall strength, σsp, with f0for three values of v0for both the calculations with only initial porosity and those with initial porosity together with stress-controlled porosity nucleation. The value of σspis given by σsp=1 2ρc∆v (23) where ρis the initial density of the material, cis the elastic wave velocity given in Eq. (19) and ∆ Vis the pullback velocity defined in Fig. 3(a). For a fixed value of v0, the value of σspwith only initial porosity initially decreases rapidly with increasing f0and then tends to saturate. For a fixed value of f0, the value of σspincreases with increasing v0. The magnitude of the increase in σspwith increasing v0 depends on the value of f0and increases with decreasing f0. When both initial porosity and stress-controlled porosity nucleation are considered, the value of σspremains nearly constant irrespective of f0andv0until f0exceeds a certain threshold value that increases with increasing v0. For smaller values of f0, i.e., below the threshold value of f0, the values of σspwith stress-controlled porosity nucleation are smaller than the corresponding values without stress-controlled porosity nucleation. However, for 19 f0greater than the v0dependent threshold value of f0,σspdecreases with increasing f0and is identical to the calculations without stress-controlled porosity nucleation. 5. Discussion We have carried out finite deformation finite element calculations to analyze the response of porous ductile materials under impact loading for values of initial porosity, f0, in the range 0≤f0≤0.05 and for values of impact velocity, v0, in the range 125 m /s≤v0≤300 m /s. Our calculations capture the plastic deformation and porosity compaction due to the propagation of the compressive stress pulse. These in turn lead to a decrease in the stress magnitude in the pulse. The results show that an increase in the initial porosity level results in an even greater decrease in the stress magnitude, but subsequently it also decreases the extent of compaction of the porous material. While an increase in impact velocity (for a fixed initial porosity value) leads to an increase in both the magnitude of the compressive stress and the extent of compaction of the porous material. Because our focus is on analyzing the effects of porosity-induced plastic compressibility on the response of porous ductile materials subjected to impact loading, we assumed that the elastic response of the material is linear. However, when a material is subjected to very high impact velocities and/or pressures, the assumption of linear elasticity may no longer be appropriate, and the bulk modulus of the material may vary with pressure. Nevertheless, experiments suggest that nonlinear elasticity only becomes significant at impact velocities greater than several km/s and pressures greater than tens of GPa (Walsh and Christian, 1955; Holmes et al., 1989; Trunin et al., 1989; Reinhart et al., 2001). Meeting either of these conditions requires parameter values an order of magnitude greater than those used in our calculations. In calculations where growth of pre-existing porosity is the only damage mechanism, the propensity for spall fracture depends on the stress magnitude in the reflected pulse and the current level of porosity in the material. Thus, for these materials there is a minimum level of initial porosity below which the current level of the porosity is too small for spall fracture to occur. There is also a maximum level of initial porosity above which the magnitude of the reflected tensile stress is too small for spall fracture to occur. Both the minimum and maximum critical levels of the initial porosity increase with increasing impact velocity, with the increase in the maximum critical porosity level being greater than the increase in the minimum critical porosity level. This results in an increase in the range of values of initial porosity over which spall fracture occurs with increasing impact velocity. When stress-controlled porosity nucleation is included in the calculations, spall fracture can occur due to stress-controlled porosity nucleation, even in the absence of initial poros- ity, as long as the tensile stress magnitude is large enough to nucleate sufficient porosity. Hence, for the values of the parameters characterizing stress-controlled porosity nucleation considered in this work, there is no minimum initial porosity level below which spall frac- ture does not occur. However, there is a maximum initial porosity level above which spall fracture does not occur, and this maximum porosity level increases with increasing impact velocity. Additionally, for the considered parameter values, these maximum initial porosity 20 values are almost the same for calculations, both with and without stress-controlled porosity nucleation. Although not presented here, we also carried out limited parametric analyses to examine the effects of the two parameters characterizing stress-controlled porosity nucle- ation, namely fNandσN. Among the two parameters, fNandσN,σNhad a more significant effect on the propensity of spall fracture to occur because the stress level must reach σN for porosity nucleation. The results of these parametric analyses showed that increasing σN results in a minimum value of the initial porosity (for a fixed value of the imposed impact ve- locity) below which spall fracture does not occur, as in the scenario without stress-controlled porosity nucleation. Also, decreasing σNresults in an increase in the value of the maximum initial porosity level above which spall fracture does not occur. If spall fracture results from the evolution of pre-existing porosity, the spall strength decreases with increasing initial porosity, and it increases with increasing impact velocity. However, when spall fracture is driven by stress-controlled porosity nucleation, over the range investigated, the spall strength is not very sensitive to the value of initial porosity or to the value of impact velocity. Although the values of the constitutive parameters considered in this work do not corre- spond to a specific material but represent a generic porous ductile material, several features of our results are consistent with the limited experimental studies of spall fracture in porous ductile materials. For example, a decrease in the spall strength with increasing initial porosity level has been observed experimentally for Ti-6Al-4V alloy (Cui et al., 2022) and aluminum (Yao et al., 2018), and an increase in the spall strength with increasing impact velocity has also been observed experimentally for Ti-6Al-4V alloy (Cui et al., 2022) and aluminum (Yao et al., 2018; Wang et al., 2007). Furthermore, experiments by Cui et al. (2022) showed that a Ti-6Al-4V alloy with 0 .29% and 0 .88% initial porosity underwent spall fracture when subjected to impact velocities of 500 m/s and 620 m/s, while the same alloy with 5 .41% initial porosity only underwent spall fracture at 620 m/s impact velocity. In our calculations, porosity enters the constitutive relation as a single scalar variable, the void volume fraction. There is no notion of shape, size or spacing of actual discrete voids that exist in real materials. A discrete void in an elastic-viscoplastic matrix, depending on the imposed strain rate and stress state, may collapse into a penny-shape (Nemat-Nasser and Chang, 1990) or fold up and collapse (Swantek and Austin, 2010), both of which will lead to large stress concentrations and can cause early fracture. However, under impact loading conditions similar to those considered here, there is experimental evidence that suggests that discrete voids collapse uniformly without evolving into a shape that induces a large stress concentration (Lovinger et al., 2020; Swantek and Austin, 2010; Lind et al., 2021). Nevertheless, the discreteness of pores may introduce micro-inertia and delay their evolution (Ortiz and Molinari, 1992; Tong and Ravichandran, 1993; Molinari and Mercier, 2001). The stabilizing effect of micro-inertia could also affect the maximum stress level and damage evolution during stress relaxation within the spall plane (Czarnota et al., 2008). Micro-inertia effects are expected to be more important for large discrete voids and for large stress-magnitudes (Sartori et al., 2015; Jacques et al., 2015; Czarnota et al., 2017). Thus, our results that ignore possible micro-inertia effects may be most relevant for materials with a distribution of fine voids that are subject to relatively moderate impact velocities. 21 The impact loading condition prescribed n our calculations is an imposed velocity for all time, such that vL= 0 for t >(trise+thold+tdecay). This differs from the experimen- tally imposed impact loading conditions where the impactor generally separates from the target material, eventually leading to a traction free surface. The imposed velocity loading condition, Eq. (5), can result in damage at the impact face for sufficiently large values of v0andf0for the calculations with stress-controlled porosity nucleation. To confirm that the impact face damage is indeed due to the imposed loading condition, the effects of a finite residual velocity ( vL>0) for t >(trise+thold+tdecay) and of increased decay time, tdecay, were investigated with and without including stress-controlled porosity nucleation. These results (not shown here) indicated that imposing a finite value of the residual velocity or increasing the decay time can alleviate the formation of impact face damage. However, neither the dependence of spall fracture on the initial porosity, f0, or the dependence on impact velocity, v0, were affected. 6. Conclusions Finite deformation finite element calculations were carried out to investigate the effects of the initial porosity on the attenuation of the compressive stress wave propagating through the material for a wide range of impact velocities, and to explore the extent to which the attenuation of the compressive stress wave can mitigate porosity-induced spall fracture. Two sets of calculations are carried out, in the first set the material contains initial porosity values ranging from 0% to 5% while in the second set the material also undergoes stress-controlled porosity nucleation. The key conclusions are: 1. The compressive hydrostatic stress pulse leads to plastic deformation and compaction of the porous ductile material, which in turn leads to a decrease in the stress magnitude in the pulse. 2. An increase in initial porosity results in a decrease in the compressive stress magnitude in the initial wave, but decreases the extent of compaction of the porous material. 3. With no stress-controlled porosity nucleation, there is a minimum level of initial poros- ity below which the remnant porosity level is so low that spall fracture does not occur. There is also a maximum level of initial porosity above which the magnitude of the reflected tensile stress in the pulse is too small for spall fracture to occur. 4. The maximum initial porosity level beyond which the tensile stress magnitude in the reflected pulse is insufficient to cause spall fracture also holds for cases with stress- controlled porosity nucleation. 5. Both with and without stress-controlled porosity nucleation, the maximum level of the initial porosity above which spall fracture does not occur increases with increasing impact velocity. 6. When spall fracture is driven by the evolution of the pre-existing porosity, the spall strength of the material decreases with increasing initial porosity and increases with increasing impact velocity. When spall fracture is driven by stress-controlled porosity nucleation, the spall strength of the material is not very sensitive to the initial porosity level nor to the value of impact velocity, for the range of values considered. 22 Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This work was supported by the U.S. Department of Energy, National Nuclear Secu- rity Administration under Award No. DE-NA0003857. SO and AS also acknowledge the European Union’s Horizon 2020 Programme (Excellent Science, Marie-Sklodowska - Curie Actions, H2020-MSCA-RISE-2017) under REA Grant agreement 777896 (Project QUAN- TIFY). The finite-element calculations were carried out using the high performance research computing resources provided by the Texas A&M University (https://hprc.tamu.edu).