Dynamic cylindrical cavity expansion in orthotropic porous ductile materials 1

Abstract


Introduction
Spherical and cylindrical cavity expansion theories have been widely used over the last 50 years to model high velocity penetration and perforation problems in which the target exhibits ductile hole growth.The cavity expansion models use the solution of the pressure for steady expansion of a cavity in an infinite medium as the contact pressure experienced by a projectile that penetrates a target.While the assumption of spherical or cylindrical flow fields in the target material during the penetration may be considered a rough approximation to the actual flow field (Yarin et al., 1995;Roisman et al., 1997;Rubin et al., 2016), cavity expansion theories have been widely accepted as a basic physical model underlying perforation and penetration mechanics in metallic materials (Masri, 2016).According to Durban and Masri (2004) and Cohen and Durban (2013a), key to the success of the cavity expansion theories is that they provide simple expressions for critical parameters like the resisting force, the ballistic limit and the penetration depth (e.g.Forrestal andWarren (2008, 2009) and Cohen et al. (2010a)).
Spherical cavity expansion approximations have been primarily used to model deep penetration problems (Forrestal and Warren, 2009).Goodier (1965) employed the dynamic solution developed by Hopkins (1960) to describe the penetration of spherical impactors into metal targets modelled as elastoplastic, incompressible solids.Forrestal et al. (1988) developed a spherical cavity expansion model for rate-independent, elastic-perfectly plastic materials that was validated against penetration depths obtained from tests in which projectiles with different geometries were impacted onto aluminum targets at velocities which ranged from 0.4 − 1.4 km/s.Forrestal et al. (1991) included the strain hardening of the material in the penetration equations developed by Forrestal et al. (1988), and showed the influence of material compressibility in the predictions of the cavity expansion model.Warren and Forrestal (1998) incorporated the strain rate sensitivity into the model of Forrestal et al. (1991), and concluded that it was necessary to take into account the compressibility, the strain hardening, and the strain-rate sensitivity of the target material to obtain good agreement with previously published depth of penetration data for striking velocities between 300 and 1200 m/s.Durban and Masri (2004) developed a dynamic spherical cavity expansion model for Drucker-Prager solids with a non-associated flow rule.The model, which assumed a self-similar expansion field, showed that plastic pressure sensitivity causes an increase in the cavitation pressure.The self-similar spherical cavity expansion approach was also used by Masri and Durban (2005) to develop penetration equations for compressible von Mises solids.The penetration depth predictions obtained with the cavity expansion model were in reasonable agreement with the experimental data reported by Forrestal et al. (1988Forrestal et al. ( , 1991) ) for different impact velocities and projectiles nose shapes.Cohen et al. (2010b) investigated the elastoplastic fields induced by a pressurized spherical cavity expanding in von Mises solids at higher velocities than those investigated by Masri and Durban (2005).They showed that when the speed of the cavity is ≈ 50% of the elastic wave velocity in a long rod, shock waves emerge in the material.Cohen and Durban (2013b) found similar results for porous materials modelled with Gurson plasticity, and showed that the increase of the material porosity favours the development of shock waves at lower cavitation velocities.
On the other hand, cylindrical cavity expansion approximations have been primarily used to model perforation problems (Forrestal and Warren, 2009).Forrestal et al. (1990) obtained perforation equations based on the plane strain expansion of a cylindrical cavity which provided residual velocity and ballistic limit predictions in agreement with experimental data obtained from tests in which conical nose projectiles were impacted onto 5083-H131 aluminium armour plates of various thickness.The analytical model considered linear elastic, power law strain hardening plastic solids.The effect of target inertia, that was neglected in the perforation equations developed by Forrestal et al. (1990), was investigated by Forrestal and Warren (2009) which showed that the effect of inertia in the perforation process increases with the target thickness.
A dynamic cylindrical cavitation model valid for moderate expansion velocities was developed by Masri and Durban (2006) for incompressible elastoplastic materials with a definite yield point.Shortly after, Masri and Durban (2009) developed solutions for orthotropic materials, including the isotropic von Mises and Tresca solids as special cases.The analysis was valid for any hardening response, including materials with no definite yield point, and accounted for elastic compressibility in the linear elastic zone.Masri (2014) also included the effect of temperature in the cylindrical cavity expansion model, and showed that material resistance to cavitation is weakened with adiabatic thermal softening.Very recently, Johnsen et al. (2017) extended the cylindrical cavity expansion theory to include the Voce strain hardening rule.
Overall, many numerical and theoretical studies addressing the perforation of isotropic metallic materials have been conducted over the last decades.However, as discussed by N'souglo et al. (2019), only few studies considering anisotropy effects on penetration problems have been reported in the literature (Grytten et al., 2009;Masri and Durban, 2009;Seidt et al., 2013;Mars et al., 2015;Koubaa et al., 2017).In this paper we adapt the dynamic spherical cavity expansion model developed by Cohen et al. (2010b) and Cohen and Durban (2013b) for isotropic solids to the case of a cylindrical cavity expanding under plane strain conditions in anisotropic porous materials.Further, following the recent work of dos Santos et al. (2019), we incorporate into the cylindrical cavity expansion model the artificial viscosity approach developed by Lew et al. (2001) to capture the shock waves that emerge at high cavitation velocities.The respective influence of anisotropy and porosity in the cavitation pressure and in the shock waves formation is exposed and discussed.In addition, new results are presented regarding the stress and density fields that evolve near the cavity for porous materials with varying degrees of mechanical anisotropy.Finally, the cylindrical cavity expansion model is used in conjunction with the Recht and Ipson (1963) formulas to assess the ballistic performance of porous anisotropic targets against high velocity perforation.The outline of the paper is as follows.Section 2 presents the constitutive model developed by Benzerga and Besson (2001).
The cylindrical cavity expansion model for porous anisotropic solids is developed in Section 3. Theoretical results for materials with different degrees of mechanical anisotropy and initial void volume fractions are presented in Section 4. We conclude in Section 5 with a summary of the main findings of this work.

Constitutive framework
The mechanical behaviour of the material is described using the elastic isotropic, plastic orthotropic constitutive model developed by Benzerga and Besson (2001) for porous ductile metals.The model con-sists of extending the constitutive framework developed by Gurson (1977) and Tvergaard (1981Tvergaard ( , 1982) ) for spherical voids to the case of an orthotropic material.For that purpose, the mechanical behaviour of the matrix material is described using Hill's model (Hill, 1948).Only the main features of the Benzerga-Besson model are presented in this section.For additional details, the reader is referred to the work of Benzerga and Besson (2001) and the recent book of Cazacu et al. (2019).
In cylindrical coordinates (r, θ, z), the flow potential reads: where q 1 and q 2 are material parameters, f is the void volume fraction, and σ yr is the matrix yield strength in the radial direction (see equation ( 4)).Moreover, σ h = 1 3 σ : 1 is the mean stress, with σ being the macroscopic Cauchy stress tensor and 1 the unit second order tensor.The constant ψ is defined as: where F0 , G0 , H0 , L0 , M0 , and Ñ0 are dimensionless anisotropy coefficients.Moreover, σe is the equivalent stress expressed as: where σ rr , σ θθ and σ zz are the radial, circumferential and axial normal stresses, and σ rθ , σ θz and σ zr are the corresponding shear stresses.Benzerga and Besson (2001) showed that if F0 = G0 = H0 = 1 2 and L0 = M0 = Ñ0 = 3 2 then ψ = 2, and equation ( 1) reduces to the flow potential derived by Gurson (1977) and Tvergaard (1981Tvergaard ( , 1982)).Following the work of Benzerga and Besson (2001), it can be shown that equation ( 1) reduces to the flow potential derived by Hill (1948) provided that f = 0 and the following relations are fulfilled: where F0 , G0 , H0 , L0 , M0 , and Ñ0 represent the definitions of the dimensionless anisotropy coefficients in equation ( 2).Moreover, F 0 , G 0 , H 0 , L 0 , M 0 , N 0 are the anisotropy coefficients introduced by Hill (1948), which are defined as: being σ 0r , σ 0θ and σ 0z the initial yield stresses of the matrix material along the radial, circumferential and axial directions, respectively, and σ 0rθ , σ 0θz and σ 0zr the corresponding initial yield stresses in shear.
Based on previous relationships, both the anisotropy parameter ψ, equation ( 2), and the equivalent stress σe , equation (3), can be given in terms of the initial yield stresses of the matrix material.Consequently, it is noticed that isotropy is recovered when The matrix yield strength along the radial direction σ yr , which enters into equation (1), is described by the following power-type relation: being n the strain hardening parameter and ε 0 a reference strain.As in Hill (1950), we assume that the strain hardening evolves proportionally for any material direction, i.e. the ratio between the current yield stresses, e.g.σyz σyr , remains constant through the deformation process.Recall that σ yz and σ yr are the matrix yield strengths along the axial and radial directions, respectively.Moreover, εp = t 0 εp (τ ) dτ is the equivalent plastic strain, where εp is the equivalent plastic strain rate in the matrix material.According to Hill's rationale (Hill, 1950), the equivalent plastic rate εp is the work conjugate of σy (see equation ( 10)), which is defined as: (5) Note that, for the sake of simplicity, we have not considered the temperature and strain rate dependences of the yield strength of the matrix material.This, in turn, facilitates to derive the theoretical model for the cylindrical cavity expansion problem presented in section 3.
The rate of deformation tensor is taken to be the sum of an elastic part, d e , and a plastic part, d p , as follows: where the elastic part is related to the rate of the stress by the following hypo-elastic law: with σ being an objective stress rate and C being the isotropic elastic tensor given by: where E is the Young's modulus, ν is the Poisson's ratio, and I is the unit deviatoric fourth order tensor.
The plastic part of the rate of deformation tensor follows the direction normal to the flow potential: where λ is the non-negative plastic flow proportionality factor.
The plastic part of the rate of deformation tensor and the equivalent plastic strain rate in the matrix material are related assuming that the rate of macroscopic plastic work is equal to the rate of equivalent plastic work in the matrix material: Moreover, assuming the incompressibility of the matrix material, the evolution of the void volume fraction is defined as: Note that void nucleation is not considered in the present analysis.Hence, if the initial void volume fraction is null, the macroscopic material is fully dense.
The formulation of the constitutive model is completed with the loading/unloading Kuhn-Tucker conditions: and the consistency condition:

Specialization of the constitutive equations
The constitutive equations presented above are specialized to describe four particular cases that will be analysed in section 4, namely: isotropic fully dense material, transversely isotropic fully dense material, isotropic porous material and transversely isotropic porous material.

• Isotropic fully dense material
The constitutive framework presented in section 2 reduces to standard von Mises plasticity.The initial material porosity f 0 is zero.Moreover, the initial yield stress is the same for any direction, which leads to:

• Transversely isotropic fully dense material
The constitutive framework presented in section 2 reduces to transversely isotropic Hill plasticity, with mechanical properties that are symmetric about the z − axis.The initial material porosity f 0 is zero.Moreover, the initial yield stress is the same for any direction in the r − θ plane, which leads to: The transverse isotropy about the z −axis leads to a constant Lankford coefficient for uniaxial traction along any direction in the r−θ plane.This constant coefficient may be calculated considering a uniaxial traction along the r-direction: The assumption of a transversely isotropic material, in addition to represent many practical applications (Hill, 1950), is particularly convenient from a modeling standpoint, since an axisymmetric cylindrical cavitation model can be adopted when the material is isotropic in the r − θ plane.
Moreover, since σ 0r = σ 0θ , the ratio between the initial yield stresses along the axial and radial directions is: The combination of equations ( 14) and ( 15) provide the following relations between R and K: According to Hill (1948), R is expected to be always positive, and therefore K ≥ 1 2 ≈ 0.7071.
Previous relations allow to express F0 , G0 , H0 and Ñ0 as a function of a single scalar parameter, R or K, as follows: Moreover, for sake of simplicity, we assume that L0 = M0 = 3 2 (the isotropic value, see discussion after equation ( 3)).This simplification has also been assumed by Cazacu et al. (2019) in the absence of sufficient experimental data.At this point, it is readily noticed that, for a transversely isotropic material, the main features associated to its plastic anisotropy are represented by single scalar parameter, R or K.Note that if K = 1 (or R = 1), then F0 = G0 = H0 = 1 2 and Ñ0 = 3 2 , and the material follows von Mises plasticity.

• Isotropic porous material
The constitutive framework presented in section 2 reduces to standard Gurson plasticity.The initial material porosity f 0 is different from zero.The isotropy of the material leads to: • Transversely isotropic porous material The constitutive framework presented in section 2 reduces to transversely isotropic Benzerga-Besson plasticity, with mechanical properties that are symmetric about the z − axis.The initial material porosity f 0 is not zero.As in the case of the transversely isotropic fully dense material, the plastic anisotropy parameters, expressed as a function of K, are: and the parameter ψ, expressed as a function of K or R, becomes: where the value assigned to the parameters L0 and M0 was 3 2 (isotropic value).If K = 1 (or R = 1), then F0 = G0 = H0 = 1 2 , Ñ0 = 3 2 and ψ = 2, and the material follows Gurson plasticity.

Material parameters
Table 1 shows the values of the initial density and the elastic constants, and the values of the parameters corresponding to the flow potential and the flow strength of the matrix material, which have been taken from Srivastava et al. (2014) andN'souglo andRodríguez-Martínez (2018).In section 4 we carry out a parametric analysis to investigate the effect of initial porosity, strain hardening and plastic anisotropy in the cavitation fields.In the parametric analysis different values of f 0 , n and K are used, see Table 1.Note that, while the initial radial yield stress is constant in the analysis developed in section 4, the initial axial yield stress σ 0z changes with K, see equation ( 15).

Theoretical model
The theoretical model developed in this paper is based on the seminal contribution of Cohen and Durban (2013b) who studied the problem of a spherical cavity expanding dynamically in an isotropic porous medium.
standard numerical methods.While Cohen and Durban (2013b) solved this problem imposing Hugoniot conditions regarding mass and momentum conservation at the shock front, in this paper we use the artificial viscosity approach developed, for instance, by Wilkins (1980) and Lew et al. (2001).This approach, which leads to smooth solutions of the set of differential equations which govern the problem and stands out for its simplicity, was applied -and validated with finite element simulations-in the recent paper of dos Santos et al. (2019) to capture the emergence of shock waves in isotropic porous materials.We are aware that there are more physically sound approaches to model shocks in porous materials; e.g.considering microinertia effects related to the rapid collapse of the voids (Czarnota et al., 2017).Nevertheless, analysing the role of microinertia in the cavitation fields of elastoplastic porous materials is out of the scope of this research, and it will be addressed in a future work.

Steady self-similar dynamic cylindrical cavity expansion
The problem is modelled as a pressurized cylindrical cavity of instantaneous radius a expanding under self-similar, steady-state conditions in an infinite porous medium, see Fig. 1.Let (r, θ, z) stand for the Eulerian coordinates of a material point with reference to a cylindrical system with the origin located at the center of the cavity denoted by O in Fig. 1.The axial direction z is under plane strain constraint.
Following Durban and Fleck (1997) and Durban and Masri (2004), we take the nondimensional radial coordinate ξ = r a as the only independent variable of the problem.Accordingly, the derivative with respect to the radial coordinate is: The internal pressure p is 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 = f 0 and ρ = ρ 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 figure is adapted from Durban and Masri (2004).
where the prime superscript denotes differentiation with respect to ξ.Moreover, as shown by Durban and Fleck (1997), the time derivative is: where v = ṙ ȧ is the dimensionless radial velocity.Since ȧ is constant (steady-state expansion), the dimensionless acceleration reads v = r ȧ .
We assume that the axes of orthotropy of the material are aligned with the radial, circumferential and axial directions, i.e. the principal directions of stress and strain coincide.Therefore, under cylindrical expansion, the only active components of the stress tensor are σ rr , σ θθ and σ zz .Similarly, the only active components of the rate of deformation tensor d are: with due to the plane strain condition imposed in the z direction.Note that, to obtain expressions ( 18) and ( 19), equations ( 16) and ( 17) have been employed.For the problem at hand, since the principal axes are fixed, i.e. no rotation is included, the rate of deformation tensor coincides with the strain rate tensor.Moreover, the hypo-elastic law introduced in equation ( 7) can be rewritten as: ) Note that, in absence of material spin, the objective derivative () in equation ( 7) has been replaced by a time derivative ().
Using the expressions for the strain rates given in ( 18)-( 20), and the additive decomposition of the rate of deformation tensor (6), the following three equations are obtained for the radial, circumferential and axial directions: where the elastic deformation rates given in equations ( 21)-( 23), and the plastic deformation rates given in equations ( 24)-( 26), have been employed.To obtain expressions ( 30)-( 32), the time derivative ( 17) and the relation 3Σ h = Σ rr + Σ θθ + Σ zz have been used.The derivative ∂ εp ∂Σyr , which is obtained from equation ( 4), is calculated as: being Σ 0r = σ 0r E the dimensionless initial yield stress of the matrix material along the radial direction.
Moreover, using equations ( 24)-( 26), the expression for the evolution of the void volume fraction given in ( 11) can be written as follows: Using the time derivative (17) and equations ( 27)-( 29), previous expression is rewritten as: This equation shows, as noted by Benzerga and Besson (2001), that the porosity rate depends on the anisotropy of the material through the constant ψ.If ψ < 2 the porosity rate is greater than for an isotropic material.On the other hand, if ψ > 2 the porosity rate decreases in comparison with an isotropic material.Benzerga and Besson (2001) also noted that, since ψ appears in the hyperbolic sine, any small variation of ψ has significant influence in the porosity rate.
Moreover, the ratio between current material density ρ and initial material density ρ 0 can be obtained from the balance of mass which, in Eulerian description, takes the form: where l is the spatial velocity gradient.In absence of material spin l = d, and, using equations ( 17)-( 20), previous expression becomes: Inserting equations ( 30)-(32) into equation ( 37), the following expression is obtained: where the definition for f given in ( 35) has been used.Following Cohen and Durban (2013b), integrating equation ( 38) from a given coordinate ξ to the elastic wave front ξ w , where f = f 0 , ρ = ρ 0 , and Σ h = 0 (see Fig. 1), a closed-form relation for the density ratio is obtained: Moreover, the balance of linear momentum along the radial direction, in Eulerian description, takes the form: where σ t rr and σ t θθ are the total radial and circumferential stresses given by the sum of the elastoplastic stresses and the viscous stresses, i.e. σ t rr = σ rr + σ v rr and σ t θθ = σ θθ + σ v θθ .On the one hand, σ rr and σ θθ are related to the elastic and plastic deformation rates by the constitutive equations presented in section 2. On the other hand, σ v rr and σ v θθ are related to the total deformation rates by means of an artificial viscosity (Lew et al., 2001) to be presented in the sequel.The artificial viscosity approach is only used in the definition of the stresses that enter into the balance of linear momentum.Employing equations ( 16)-( 17), the linear momentum equation ( 40) can be rewritten as: where is the dimensionless cavitation velocity defined as the ratio between the cavity expansion velocity ȧ and the elastic wave velocity in a long rod E /ρ 0 .
Following the development presented in the recent work of dos Santos et al. ( 2019), which is based on previous works of Wilkins (1980) and Lew et al. (2001), the dimensionless viscous stresses, Σ v rr = σ v rr E and , that enter into equation ( 41), are calculated as: and being the radial and circumferential artificial viscosity parameters, respectively.In expressions ( 44) and ( 45 2019), and satisfactory agreement between both approaches was found for a wide range of cavitation velocities.
The integration is performed over the independent variable ξ, from the elastic wave front ξ w up to the cavity wall ξ = 1.An iterative shooting method is used to enforce the velocity boundary condition at the cavity wall.The initial conditions, with an unknown constant, and the location of the elastic wave front are determined adapting the elastic solution derived by Masri and Durban (2009) to the constitutive framework used in this paper.

Elastic zone solution
The solution of Masri and Durban (2009) is based on the fact that, within the elastic zone, where the plastic strain rates are zero, the dimensionless stresses are expected to be small Thus, equations ( 30), ( 31), ( 32), ( 37) and ( 41) can be solved analytically to obtain the following expressions for the radial, circumferential and axial elastic stresses: and where B and C are integration constants to be determined.The parameter M = ȧ/C E can be regarded as the cavity expansion Mach number (Masri and Durban, 2009), being Moreover, the corresponding material density ρ and velocity v are calculated by the following expressions: and where the approximation exp has been used.

Determination of the integration constants
The elastic wave front is located at ξ = ξ w , where v = 0. Thus, equation ( 52) leads to: Knowing that at the elastic wave front Σ rr = 0, previous expression is inserted into equation ( 48) to obtain: Moreover, the constant C is determined imposing field continuity conditions at the (unknown) elastic / elastoplastic interface ξ = ξ i , where equation ( 1) is rewritten as: Using equations ( 48)-( 50), the mean stress Σ h and the equivalent stress Σe are calculated as: where Then, equations ( 56) and ( 57) are inserted into (55), and the resulting expression is solved in terms of the constant C. Within the shooting method, the elastic / elastoplastic interface location is iterated until the velocity boundary condition (v = 1) is reached at the cavity wall.

Results
This section of the paper is split into two parts.Subsection 4.1 shows results for fully dense materials, for a wide range of dimensionless cavitation velocities 0.04 ≤ m ≤ 0.6.For the material parameters of The results obtained for porous and fully dense materials are compared.Fig. 3(a) shows that Σ rr is negative at any value of the radial coordinate ξ.The radial stress decreases monotonically with ξ.The kink shown by the Σ rr − ξ curves at ξ ≈ 8 corresponds to the interface between the elastic zone (on the right) and the elastoplastic zone (on the left).Note that, for this cavitation velocity, the specific location of the elastic / elastoplastic interface, from now on denoted by ξ i (see Fig. 1), slightly depends on K, as it will be discussed later on (see Fig. 5).The influence of the anisotropy parameter K in the radial stress is greater near the cavity, where Σ rr increases as K decreases, in agreement with the results presented in Fig. 2 (Σ rr = −P at the cavity wall).As we move away from the cavity, the Σ rr − ξ curves get closer to each other until they intersect at ξ ≈ 6.For 6 ξ ≤ ξ i the radial stress is slightly greater as the anisotropy parameter K increases.In the elastic range, ξ i ≤ ξ ≤ ξ w , the influence of K in the radial stress is negligible.Moreover, the location of the elastic wave front ξ w = 11.6 is independent of the plastic properties of the material.Fig. 3(c) shows that the axial stress Σ zz is also negative for any value of the radial coordinate ξ.

Fully dense materials
Similarly to the circumferential stress, the axial stress first increases rapidly, reaches a maximum, and then decreases more slowly.Near the cavity the axial stress is slightly greater as K decreases.As we move away from the cavity, the Σ zz − ξ curves get gradually closer and eventually intersect at ξ ≈ 7. From ξ ≈ 7 until the elastic / elastoplastic interface, the axial stress is slightly greater as K increases.In the elastic zone, Fig. 3(d) shows that the equivalent stress Σ e decreases monotonically with the radial coordinate ξ.Note that the decrease of Σ e is particularly fast near the cavity.In fact, according to the asymptotic analysis developed by Durban and Masri (2004), Σ e should be unbounded at the cavity wall.However, this is not the case for the calculations performed in this work.It seems that, due to the numerical solution, the equivalent stress has a finite value at ξ = 1 (dos Santos et al., 2019).We have checked that the specific value of Σ e for ξ = 1 slightly depends on the tolerance of the shooting method used to solve the problem (see section 3).Near the cavity, the equivalent stress is slightly greater as K decreases.The difference between the curves obtained for different values of K is gradually reduced with ξ.The four Σ e − ξ curves intersect at ξ ≈ 6.5 and their order is reversed until the the elastic / elastoplastic interface is reached.In the elastic zone the influence of K in the Σ e − ξ curves is small.is slightly greater than the axial stress, which is in turn slightly greater than the circumferential stress.Fig.

4(d)
shows that the equivalent stress is a decreasing function of the radial coordinate.Note that, unlike the radial, circumferential and axial stress profiles, the Σ e − ξ curves display a concave-upward shape in the elastoplastic zone.Note also that in the elastoplastic zone the values of Σ e are one order of magnitude smaller that the values Σ rr , Σ θθ and Σ zz , i.e. the mean stress is one order of magnitude greater than the equivalent stress.In addition, the influence of K in the equivalent stress is more important than in the radial, circumferential and axial stress profiles, particularly, near the cavity wall.
Note that all the stress fields shown in the paper correspond to the elastoplastic stresses.The viscous stresses lack of physical meaning and their role in the theoretical model is limited to facilitate the integration of the governing equations at high cavitation velocities.The reader is referred to Appendix A and the work of dos Santos et al. ( 2019) to obtain additional insights into the respective contribution of elastoplastic and viscous stresses into the total stresses that enter in equation ( 41). for m ≥ 0.37.As the cavitation velocity increases ξ p decreases, i.e. the shock wave becomes closer to the cavity.Note that the difference between ξ i and ξ p is small for any value of m considered (differences in the second decimal place), and becomes smaller as m increases.This observation agrees with the theoretical and numerical results reported by Cohen et al. (2010b) who showed for porous isotropic materials that increasing m will eventually lead to ξ i = ξ p .

Porous materials
Fig. 6 shows the variation in dimensionless applied pressure P with dimensionless cavitation velocity m for four different values of the anisotropy parameter K: 0.75, 1 (Gurson plasticity), 1.5 and 2. Different values of initial void volume fraction f 0 and strain hardening coefficient n are considered.The results for the fully dense material and K = 2 are also included.
-    coefficient has been increased to 0.5.Note also that the scale of the y-axis has been increased one order of magnitude.The increase of the strain hardening boosts the equivalent stress at the cavity wall.Note also that, in the elastoplastic zone, the variation of K influences the equivalent stress such that, as K decreases, the equivalent stress is slightly greater.However, in the elastic zone, the impact of the anisotropy parameter is imperceptible.cavitation velocity has been increased to 0.6.The influence of the anisotropy parameter on the equivalent stress in the elastoplastic zone is important.The shock wave, which is located at similar radial coordinate than the elastic / elastoplastic interface (see Fig. 13(c)), is illustrated by the sudden drop of Σ e located at ξ ≈ 1.68.In the elastic zone, the influence of K is negligible.
Note that in the four plots of Fig. 7 the curves of f 0 = 0 and f 0 = 0.001 for K = 2 practically lie on top of each other, which suggests that 0.1% porosity does not affect significantly the cavitation fields, in agreement with the results presented in Figs.6(a) and 6(c).Next, we analyse the cavitation fields for a greater value of initial porosity.The comparison between the fully dense material and the porous ones reveals that considering 10% of initial porosity in the calculations has an important impact on the predicted stress fields.For instance, Fig. 8(a) shows that the radial stress for the fully dense material is greater than for the porous materials for any value of ξ.Note also that for the fully dense material the radial stress decreases smoothly with ξ.
However, the Σ rr − ξ curves for the porous materials show a steep decline of the radial stress at ξ ≈ 2.8.
For ξ < 2.8 the radial stress is greater as K decreases, and vice versa for ξ > 2.8.Fig. 8(b) illustrates that the circumferential stress is also significantly greater for the fully dense material.Unlike in Fig. 8(a), the Σ θθ − ξ curve for K = 2 is above the curves for K = 1.5, 1 and 0.75 for any radial coordinate.Fig. 8(c) shows the profiles of axial stress, which are similar, both quantitatively and qualitatively, to the profiles of radial stress.For the porous materials the Σ zz − ξ curves also show a steep decline of Σ zz at ξ ≈ 2.8 such that for ξ < 2.8 the axial stress is slightly greater as K decreases, and vice versa for ξ > 2.8.The equivalent stress profiles are shown in Fig. 8(d).Note that, for the sake of clarity, the scale of the y-axis has been reduced in comparison with Figs.8(a)-(b)-(c).Near the cavity, the equivalent stress is 2 to 4 times smaller than the radial, circumferential and axial stresses, i.e. the mean stress is greater than the equivalent stress.
For small values of the radial coordinate, the equivalent stress is smaller for the fully dense material than for the porous materials.The Σ e − ξ curve of the fully dense material intersects the Σ e − ξ curves of the porous materials at ξ ≈ 2. The specific location of the intersection depends on the parameter K.For greater values of ξ, the equivalent stress for the fully dense material is significantly greater.Note also that, while the equivalent stress decreases monotonically with the radial coordinate for the fully dense material, this is not the case if f 0 = 0.1.For the porous materials the equivalent stress first decreases, reaches a minimum at ξ ≈ 2.8 (the radial coordinate corresponding to the steep decline of the stress in Figs.8(a) and 8(c)), and then slightly increases until a kink appears in the Σ e − ξ curves.The minimum indicates the incipient formation of a shock wave (a small increase in m will lead to the emergence of a shock wave, see Fig. 13(b)) and its location is practically independent of K.The kink corresponds to the elastic / elastoplastic interface, and its location depends on the value of K, as further discussed in Fig. 13(b).To the left of the kink, the equivalent stress is slightly greater as K decreases.On the contrary, in the elastic zone, to the right of the kink, the equivalent stress in greater as K increases.0.75, 1 (Gurson plasticity), 1.5 and 2. The corresponding values of ψ are: 2.052, 2, 2.18 and 2.466.Recall that the porosity rate depends on the anisotropy through the constant ψ.As in Fig. 8 the strain hardening exponent n, the cavitation velocity m, and the initial void volume fraction are 0.1.The results for the fully dense material and K = 2 are also included in the graph which shows the density profiles.Fig. 9(a) shows that the f − ξ curves feature a sigmoidal shape.Near the cavity wall, where large compressive stresses appear, the porosity is zero.The porosity takes positive values for a value of the radial coordinate that depends on K such that the greater the value of ψ, the smaller the value of ξ for which f > 0. Note that the f − ξ curves intersect at ξ ≈ 2.86, which is the radial coordinate for which the shape of the curves turns from concave-upward to concave-downward.For ξ < 2.86 the porosity is greater as ψ increases, and vice versa for ξ > 2.86.The porosity reaches the initial value f 0 = 0.1 at the elastic / elastoplastic interface.Fig. 9(b) shows that the dimensionless material density ρ/ρ 0 decreases with the radial coordinate.For the porous materials, the ρ/ρ 0 − ξ curves feature a concave-downward shape for ξ < 2.86, and a concave-upward shape for ξ > 2.86.Note that for ξ < 2.86 the density is greater as K decreases, and vice versa for ξ > 2.86.Unlike in the case of the porous materials, for the fully dense material and K = 2, the ρ/ρ 0 − ξ curve features a concave-upward shape from the cavity wall until the elastic wave front.Near the cavity, for ξ < 2.86, the density in the fully dense material is smaller than in the porous materials.The plastic compressibility of the porous materials is responsible for the increased density close to the cavity.However, at ξ ≈ 2.86 the ρ/ρ 0 − ξ curves for the porous and the fully dense material intersect such that for ξ > 2.86 the density in the fully dense material is greater.These results reveal that considering 10% of initial porosity in the calculations has an important impact on the predicted density fields.We are aware that most metals and alloys display values of initial porosity smaller than 10%.However, considering a wide range of values of f 0 helps to determine the effect of porosity in the elastoplastic fields and in the formation of shock waves.1.5 and 2. The strain hardening exponent is n = 0.1.The cavitation velocity is m = 0.1.In the plot which shows the dimensionless material density ρ/ρ 0 , results for the fully dense material and K = 2 are also included.
the stresses decrease smoothly, featuring a concave-downward shape.At ξ ≈ 1.54 there is a sudden drop of the stresses which indicates the emergence of a shock wave.The location of the shock wave practically coincides with the location of the elastic / elastoplastic interface, i.e. ξ p ≈ ξ i , as it will shown later in Fig. 13(b).The comparison of the results for the porous materials with the fully dense material show that, in the latter case, the radial, circumferential and axial stresses are significantly higher.Note also that for the fully dense material the drop of stress that indicates the shock wave occurs at a greater value of ξ.Moreover, Fig. 10(d) shows that the evolution of the equivalent stress for the porous materials also depends very little on the parameter K.Note that the equivalent stress is ≈ 50 times smaller than the radial, circumferential and axial stresses, i.e. the mean stress is ≈ 50 times greater than the equivalent stress.Moreover, unlike in the case of the fully dense material, the equivalent stress for the porous materials is not a decreasing function of the radial coordinate.The development of the shock wave leads to a sudden decrease of Σ e that partially bounces back until reaching a local maximum that coincides with the elastic / elastoplastic interface (i.e.ξ p ≈ ξ i as mentioned before).Note that the equivalent stress for the fully dense material is below the equivalent stress for the porous materials until, in the latter case, the shock wave emerges.The strain hardening exponent is n = 0.1.The cavitation velocity is m = 0.6.The results for the fully dense material and K = 2 are also included.Fig. 11(a) shows that from the cavity wall, until the shock wave location, the porosity is zero.Within the shock the porosity increases rapidly and reaches the initial value f 0 = 0.1 at the elastic / elastoplastic interface.Within the elastic range, for ξ > ξ i , the value of f is constant.Fig. 11(b) shows that, for the porous materials, the influence of K on the dimensionless density is negligible.The density decreases smoothly as we move away from the cavity wall featuring a concave-downward shape until the shock wave location is reached.The shock leads to a sudden drop of the density.Moreover, the comparison between the results for f 0 = 0.1 and f 0 = 0 shows that, for the latter case, the density is 10% smaller near the cavity and the shock wave is located at a greater value of ξ (as mentioned before).the Σ e − ξ curves for the porous materials and the fully dense material are very similar.In addition, the influence of plastic anisotropy in the Σ e − ξ curves for the porous materials is small.In contrast, Fig. 12(b) shows that the effect of porosity and plastic anisotropy in the equivalent stress is important for m = 0.6.
For this cavitation velocity, the Σ e − ξ curve for the fully dense material is below the curves for the porous materials until the shock wave emerges for the latter at ξ p ≈ 1.55.Note that, unlike in Fig. 10 Fig. 12(c) shows that, for m = 0.1, the porosity decreases nonlinearly as the cavity is approached.Note that, at the cavity wall, while the porosity is small, it is not zero.Near the cavity the porosity is slightly smaller as the parameter ψ decreases.The porosity profiles depicted in Fig. 12(d) for m = 0.6 show that f = 0 if ξ < ξ p and f = 0.1 if ξ > ξ p .Within the shock wave there is an abrupt variation of the porosity.Fig. 13 shows that, irrespective of the values of f 0 and n, the elastic / elastoplastic interface and the shock wave are located closer to the cavity as m increases.In addition, we have that ξ i ≈ ξ p for the highest cavitation velocities investigated.Moreover, note that the anisotropy parameter K hardly affects the shock wave location.On the other hand, the elastic / elastoplastic interface is located further from the cavity as K increases.The influence of K on ξ i is more important as m decreases.These results find qualitative agreement with the data reported in Fig. 5 for fully dense materials.The comparison between Figs. 13(a) and 13(b), and between Figs. 13(c) and 13(d), shows that the increase of the initial porosity brings the elastic / elastoplastic interface closer to the cavity, boosts the influence of K on ξ i , and favours the emergence of shock waves at lower cavitation velocities.Note that the increase of the strain hardening exponent n has the opposite effect.aluminum targets tested in Rosenberg and Forrestal (1988), the theoretical results for n = 0.1 are in good agreement with the experiments.Rosenberg and Forrestal (1988).

Summary and concluding remarks
In this paper we have developed a theoretical model to study the self-similar elastoplastic fields induced by a cylindrical cavity expanding dynamically in a rate-independent porous medium.The main novelties of the model with respect to the work of Cohen and Durban (2013a) are that it takes into account the plastic anisotropy of the material, and it uses viscous damping to capture the shock waves that emerge at high cavitation velocities.The viscous damping is introduced using the artificial viscosity approach proposed by Lew et al. (2001).The plastic anisotropy is included using the constitutive framework developed by Benzerga and Besson (2001) to model the porous medium.We have investigated materials with several degrees of anisotropy, K = 0.75, 1.5 and 2, and different values of initial porosity, f 0 = 0.001, 0.01 and 0.1.The special cases for which the material is fully dense, f 0 = 0, and/or isotropic, K = 1, have also been considered.Different strain hardening exponents have been used in the calculations: n = 0.1, 0.3 and 0.5.In addition, a wide range of cavitation velocities has been investigated 0.04 ≤ m ≤ 0.6, and the influence of material porosity and anisotropy on the emergence and development of shock waves has been analysed.Finally, the cylindrical cavity expansion model has been used to calculate the Recht-Ipson curves and assess the ballistic performance of porous anisotropic targets against high velocity perforation.The main outcomes of the research are summarized below: • The cavitation pressure P increases as the parameter K decreases, i.e. as the axial yield stress decreases.The effect of anisotropy on the cavitation pressure is greater for low cavitation velocities and it becomes more important as the material strain hardening n increases.The strain hardening also boosts the cavitation pressure.
• The formation and location of shock waves is practically independent of the anisotropy parameter K.In contrast, the initial porosity is shown to promote the emergence of shocks for lower cavitation velocities.On the contrary, the increase of the material strain hardening delays the shocks formation.
Irrespective of the values of K, f 0 and n, as the cavitation velocity increases, the shock forms closer to the cavity.
• The elastic / elastoplastic interface is closer to the cavity wall as the cavitation velocity m increases.
At high cavitation velocities, the interface and the shock wave are located at similar radial coordinate.
The decrease of the anisotropy parameter K and the increase of the initial porosity f 0 also move the elastic / elastoplastic interface closer to the cavity.The effect of anisotropy is particularly important at low cavitation velocities.On the contrary, the increase of the strain hardening parameter n reduces the effect of K on the interface location, which is situated further from the cavity as n increases.
• The radial, circumferential and axial stresses, Σ rr , Σ θθ and Σ zz , are of the same order of magnitude than the equivalent stress Σ e for low cavitation velocities.However, the mean stress Σ h increases significantly with m such that, for the greatest cavitation velocity tested m = 0.6, within the elastoplastic zone, Σ h can be up to fifty times greater than Σ e for low values of the strain hardening exponent.
While Σ rr , Σ θθ and Σ zz are greater for the fully dense material than for the porous materials, the opposite behaviour is observed for the equivalent stress Σ e near the cavity wall.
• The increase of the cavitation velocity m decreases the influence of the plastic anisotropy in the radial, circumferential and axial stresses.On the contrary, the anisotropy parameter K plays a meaningful role in the equivalent stress for any value of m considered.The specific influence of K on the equivalent stress depends on the radial coordinate ξ, the strain hardening n and the initial porosity f 0 .
• The porosity f decreases as the cavity is approached.For low cavitation velocities the evolution of f with the radial coordinate is smooth, generally being the porosity in the material greater as ψ increases.The increase of the strain hardening slows down the porosity decrease.For high cavitation velocities, the variation of porosity is localized within the shock wave such that for ξ < ξ p the porosity is zero and for ξ > ξ p the porosity is f 0 .
• The material density, near the cavity wall, is greater for the porous materials than for the fully dense.This behavior is attributed to the plastic compressibility of the porous materials.Far from the cavity, the opposite behaviour is observed.For the porous materials, the influence of the anisotropy parameter K in the density fields decreases as the cavitation velocity increases.For low values of m, near the cavity, the density is smaller as K increases.Far from the cavity, the trend is reversed.
• The Recht-Ipson curves show that the residual velocity is greater for K = 2 than for K = 0.75, which revels the influence of anisotropy on the ballistic performance of metallic targets.For n = 0.1, decreasing K from 2 to 0.75 increases the ballistic limit ≈ 8%.For n = 0.5, the increase of the ballistic limit velocity is significantly greater ≈ 23%.Moreover, increasing the initial porosity f 0 from 0 (fully dense material) to 0.1 increases the residual velocity, and thus decreases the energy absorption capacity of the target.The increase of the strain hardening exponent n from 0.1 to 0.5 has the opposite effect.
Appendix A. The respective contribution of elastoplastic and viscous stresses to total stresses and n = 0.5.Note that the elastoplastic radial stress Σ rr was denoted simply as radial stress in Figs. 3, 4, 8 and 10.For ξ < ξ p the contribution of the viscous radial stress to the total radial stress is very small and the ratio Σ v rr/Σ rr is less than 0.001.The emergence of the shock leads to a sharp increase of Σ v rr , such that the ratio Σ v rr/Σ rr reaches values greater than 1.Note that the viscous stress increases the width of the shock wave.The viscous stress quickly drops to zero for values of the radial coordinate greater than ξ p .
The observed behavior is readily understood relying on equation ( 42) which shows that gradients of velocity (which is the case within the shock wave) lead to large values of Σ v rr . -

Figure 1 :
Figure1: Scheme of self-similar field in dynamic expansion of a circular cylindrical cavity in an elastoplastic porous medium.The internal pressure p is 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 = f 0 and ρ = ρ 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 figure is adapted fromDurban and Masri (2004).
), c 0 ≈ 2 and c L ≈ 1 are constant parameters, and ∆ξ is a given increment along the dimensionless coordinate ξ.The parameters c 0 and c L play an important role in the artificial viscosity method.If the values of c 0 and c L are too small, the emergence of shock waves leads to discontinuities in the elastoplastic fields, and no solution is reached.If the values of c 0 and c L are too large, the method provides unrealistic results, overestimating the actual stresses in the material.Nevertheless, the values used in this paper, c 0 ≈ 2 and c L ≈ 1, which are commonly used in the literature, provide reasonable results.Note that theoretical results using the artificial viscosity method were compared with finite element simulations in dos Santos et al. (

Fig. 2 Figure 2 :
Fig.2shows the variation in dimensionless applied pressure P with dimensionless cavitation velocity m for four different values of the anisotropy parameter K: 0.75, 1 (von Mises plasticity), 1.5 and 2. Recall that for all the calculations reported in this subsection f 0 = 0 and n = 0.3.For the four values of K investigated the cavitation velocity increases with the applied pressure, the P − m curve featuring a concave-upward shape.Fig.2reveals that the increase of the anisotropy parameter K shifts downwards the P −m curve.As K increases, less pressure is required to reach a given cavitation velocity.In other words, for a fixed value of the initial radial yield stress, the cavitation pressure decreases as the initial axial yield stress increases,

Fig. 3
Fig. 3(b)  shows that the circumferential stress Σ θθ is compressive at any value of the radial coordinate ξ.As we move away from the cavity, the circumferential stress first increases rapidly, reaches a maximum,

Fig. 4 Figure 3 :
Fig.4illustrates the radial, circumferential, axial and equivalent stresses (Σ rr , Σ θθ , Σ zz and Σ e ), versus the dimensionless radial coordinate ξ.The difference with Fig.3is that the cavitation velocity has been increased to m = 0.6.Fig. 4(a) shows that the radial stress is negative and decreases with ξ.Unlike in the case of m = 0.1 illustrated in Fig. 3(a), for m = 0.6 the Σ rr − ξ curves feature a concave-downward shape.The sudden drop of the radial stress indicates the emergence of a shock wave at ξ p ≈ 1.65.For this value of m the locations of the shock wave and the elastic / elastoplastic interface are very close to each other ξ i ≈ ξ p , as it will be shown in Fig. 5.Note that the change in the radial stress induced by the shock, while abrupt, is not a discontinuity, i.e. the shock has a finite width due to the artificial viscosity included in the model.The influence of the anisotropy parameter K in the radial stress profiles is very small.Similar qualitative results are shown in Figs.4(b) and 4(c) for the circumferential and axial stresses.Quantitatively, the radial stress

Fig. 5
Fig.5shows the variation in the locations of the elastic / elastoplastic interface ξ i and the shock wave ξ p with the dimensionless cavitation velocity m.The location of the shock wave is taken as the radial coordinate for which v ≤ −1.4.The condition v ≤ −1.4 is chosen arbitrarily.The value −1.4 is chosen because it was not possible to solve the problem without artificial viscosity for v ≤ −1.4.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.The elastic / elastoplastic interface gets closer to the cavity as the cavitation velocity increases.Note that the influence of the plastic anisotropy parameter K in the location of the elastic / elastoplastic interface decreases as m increases.The location of the shock wave is virtually independent of the anisotropy parameter K.The shocks emerge, for the four values of K investigated,

Fig. 6
Fig. 6(a)  shows results for f 0 = 0.001 and n = 0.1.The P − m curves obtained for different values of K practically lie on top of each other.Moreover, attending to the results for the fully dense material and K = 2, it is inferred that the influence of this initial porosity in the cavitation pressure is negligible.In Fig.6(b) the initial porosity has been increased to 0.1.The results for different values of K also show no meaningful differences.On the other hand, this value of initial porosity leads to a drop in the cavitation pressure which becomes increasingly important as m increases.Note that this large value of initial porosity f 0 = 0.1 leads to a meaningful decrease of the material flow stress.Fig.6(c) shows results for f 0 = 0.001 and n = 0.5.The difference with Fig.6(a) is that the strain hardening exponent n has been increased to 0.5.The increase of n shifts upwards the P − m curves because the strain hardening exponent increases the material flow stress.In addition, unlike in the case of n = 0.1, the parameter K does influence the P − m curves, especially for the lowest values of m considered.For any given cavitation velocity, the cavitation pressure increases as K decreases.Moreover, for K = 2, the difference between the results for f 0 = 0 and f 0 = 0.001 is imperceptible.In Fig.6(d) the initial porosity has been increased to 0.1, keeping the strain hardening exponent n equal to 0.5.The increase of porosity from f 0 = 0.001 to f 0 = 0.1 has shifted downwards the P − m curves.

Figure 6 :
Figure 6: Porous materials.Variation in dimensionless applied pressure P with dimensionless cavitation velocity m for four different values of the anisotropy parameter K: 0.75, 1 (Gurson plasticity), 1.5 and 2. Different values of initial void volume fraction f 0 and strain hardening coefficient n are considered: (a) f 0 = 0.001 and n = 0.1, (b) f 0 = 0.1 and n = 0.1, (c) f 0 = 0.001 and n = 0.5 and (d) f 0 = 0.1 and n = 0.5.The results for the fully dense material and K = 2 are also included.

Figure 7 :
Figure 7: Porous materials: f 0 = 0.001.Dimensionless equivalent stress Σ e versus dimensionless radial coordinate ξ for four different values of the anisotropy parameter K: 0.75, 1 (Gurson plasticity), 1.5 and 2. Different values of cavitation velocity m and strain hardening coefficient n are considered: (a) m = 0.1 and n = 0.1, (b) m = 0.6 and n = 0.1, (c) m = 0.1 and n = 0.5 and (d) m = 0.6 and n = 0.5 .The results for the fully dense material and K = 2 are also included.

Fig. 8
Fig. 8 shows the radial, circumferential, axial and equivalent stresses versus the dimensionless radial coordinate.Calculations with four different values of K are presented.Both the strain hardening exponent

Fig. 9 Figure 8 :
Fig. 9 illustrates the variation of the void volume fraction f and the dimensionless material density ρ/ρ 0 with the dimensionless radial coordinate ξ.Calculations with four different values of K are presented:

Fig. 10 Figure 9 :
Fig. 10 illustrates the variation of the radial, circumferential, axial and equivalent stresses with the dimensionless radial coordinate.Calculations with four different values of K are presented.The strain hardening exponent is n = 0.1, and the initial void volume fraction is 10%.The results for the fully dense material and K = 2 are also included.The difference with the results presented in Fig. 8 is that the cavitation velocity m has been increased to 0.6.Figs.10(a)-(b)-(c) show that the profiles of radial, circumferential and axial stress for the porous materials are very similar, qualitatively and quantitatively, and practically independent of the anisotropy parameter K.The value of Σ rr , Σ θθ and Σ zz at the cavity wall is ≈ −0.42.As we move away from the cavity,

Fig. 11
Fig.11shows the void volume fraction f and the dimensionless density ρ/ρ 0 versus the dimensionless radial coordinate ξ.Calculations with four different values of K are presented.The dimensionless density ρ/ρ 0 for the porous materials is compared with the fully dense material and K = 2.As in Fig.10, the strain hardening exponent n and the initial void volume fraction are 0.1.The difference with the results presented in Fig.9is that the cavitation velocity m has been increased to 0.6.

Figure 11 :
Figure 11: Porous materials: f 0 = 0.1.Void volume fraction f and dimensionless material density ρ/ρ 0 versus dimensionless radial coordinate ξ. Results are shown for four different values of the anisotropy parameter K: 0.75, 1 (Gurson plasticity),1.5 and 2. The strain hardening exponent is n = 0.1.The cavitation velocity is m = 0.6.In the plot which shows the dimensionless material density ρ/ρ 0 , results for the fully dense material and K = 2 are also included.

Fig. 12
Fig.12shows, for porous materials, the dimensionless equivalent stress Σ e , void volume fraction f and dimensionless density ρ/ρ 0 versus the dimensionless radial coordinate ξ.The initial void volume fraction is f 0 = 0.1.As in previous graphs, calculations for K = 0.75, 1, 1.5 and 2 are presented.Recall that the corresponding values of ψ are: 2.052, 2, 2.18 and 2.466.In Figs.12(a)-(c)-(e) the cavitation velocity m is 0.1 and in (b)-(d)-(f) is 0.6.The difference with the results presented in Figs. 8, 9, 10 and 11 is that the strain hardening has been increased to n = 0.5.

Fig. 12
Fig. 12(a)  shows that the equivalent stress for m = 0.1 decreases monotonically with ξ.Note that this was not the case in Fig.8(d)where the strain hardening exponent n was 0.1.Moreover, unlike in Fig.8(d), (d), the equivalent stress is monotonic in the shock region.Fig.12(b) also shows that, for the porous materials and ξ < ξ p , the equivalent stress is smaller as K increases.The influence of the anisotropy parameter K in the Σ e − ξ curves is greater in Fig.12(b) than in Fig.10(d) where the strain hardening exponent was 0.1.

Fig. 12
Fig. 12(e)  shows the density profiles for m = 0.1.For the porous materials, the ratio ρ/ρ 0 first increases with ξ near the cavity wall, reaches a maximum and then decreases smoothly.The influence of K in the density profiles is small.For ξ < 3.90 the density of the porous materials is greater than the density of the fully dense material.The opposite behaviour occurs for ξ > 3.90.Fig.12(f) shows that the density profiles for the porous materials and m = 0.6 are hardly dependent on the anisotropy parameter K.The density for the fully dense material is ≈ 10% smaller than the density for the porous materials until the shock wave emerges for the latter.Note that similar qualitative and quantitative results were shown for n = 0.1 in Fig.11(b).

Figure 14 :
Figure 14: Porous materials: f 0 = 0.1.Residual velocity v r versus impact velocity v 0 (Recht-Ipson curves).Results are shown for two different values of the anisotropy parameter K: 0.75 and 2. Two values of the strain hardening exponent are considered: (a) n = 0.1 and (b) n = 0.5.Results for the fully dense material and K = 2 are also included.The theoretical results are compared with the experimental data reported byRosenberg and Forrestal (1988).

Table 1 ,
this range of values of m corresponds to dimensional cavitation velocities ranging between 0.1214 and 1.8209 Km/s.Four different values of the anisotropy coefficient, K = 0.75, 1, 1.5 and 2, are used in the calculations.Recall from equation (15) that the minimum possible value of K is 1 /2.The strain hardening exponent is n = 0.3.Subsection 4.2 shows results for porous materials with two different initial void volume fractions f 0 = 0.001 and 0.1.The cavitation velocities and anisotropy parameters investigated are the same as in subsection 4.1.Two different values of the strain hardening exponent n are considered: 0.1 and 0.5.

Table 2 :
Ballistic limit velocities v b for different values of the anisotropy coefficient K, the initial porosity f 0 and the strain hardening exponent n.These results are obtained using equation (58), in conjunction with the cylindrical cavity expansion model.The impact velocity is calculated using the cavity expansion velocity in the theoretical model v 0 = ȧ/tan ϕ.The effective resisting stress σ rs is taken as the quasi-static cavitation pressure p s in the theoetical model for m = 0.04.