Spherical Cavity Expansion Approach for the Study of Rigid-Penetrator’s Impact Problems

Abstract

In recent years, Spherical Cavity Expansion (SCE) theory has been extensively utilized to model dynamic deformation processes related to indentation and penetration problems in many fields. In this review, the SCE theory is introduced by explaining the different mathematical features of this theory, its solution, and a potential application to model the penetration of a rigid penetrator into a deformable target. First, a chronologically literature review of the most common models used to study this kind of penetration problems is introduced, focusing on the SCE theory. Then, the engineering model of penetration is presented using the SCE approach. The model is then compared and validated with some FE numerical simulations and with previous penetration results. It is concluded that this engineering model based on the SCE theory can be utilized to predict the projectile deceleration and penetration depth into the semi-infinite and finite targets impacted by rigid penetrators.

1. Introduction

Impact-Penetration Dynamics (IPD) is the mechanics of deformation caused by two or more colliding bodies when one of them permanently changes its shape, or its integrity, due to the high forces developed during impact. The study of IPD covers a wide range of problems and applications [1].

To understand the penetration process, many researchers have grouped various IPD problems into a finite number of categories. For example, Jonas and Zukas [2] classified them according to impact velocity, where the response of materials at different velocities is a function of the strain-rates from impact. Hopkins and Kolsky [3] categorized penetration processes by the different physical mechanisms that come into play, which determine the nature of the penetration process. Examples of such mechanisms are plastic, elastic, and thermal behaviors, as well as a decrease in compressibility. Wilbeck [1] suggested that a pair of variables must be chosen to classify impact regimes: (1) the relative magnitudes of the material’s strengths, both penetrator and target strengths (${\sigma }_P$ and ${\sigma }_T$); and (2) the impact pressures P, where pressure is proportional to the material’s density and to the square of the impact’s velocity ($P\sim \rho v^2$). However, adding more variables in the classification makes the analysis more confused. For this reason, Backman and Goldsmith [4] simplified the study of rod penetration according to the impact velocity of the projectile and the respective penetration depth in the impacted target, as shown in Figure 1. This schematic only shows a wide concept for the study of IPD, and it could involve different scenarios, from low velocity impacts to high energy projectiles.

In Figure 1, four different penetration regimes can be analyzed. In the first region, relative velocities are low, and thus, the projectile can remain intact, which produces a deep crater whose diameter is only slightly greater than the projectile’s. In the second region, the projectile deforms while penetrating the target. As the impact velocity increases, a point is reached at which the crater depth increases with the velocity. At this point, the semi-fluid region is reached, where the materials can be modeled as fluids. Over this point, the penetration regime is known as hyper-velocity impact. Finally, the last region involves impacts of relatively ultra-high velocity, the energy of which can generate an explosive reaction in the bodies. Notice that for very low impact velocities (not shown in this figure) impact events are low energy, which could not produce plastic deformation in the bodies. Therefore, this region is not subjected to study in IPD.

Different formulations have been proposed to study each regime in Figure 1; for example, cavity expansion theory [6] for studying Regimes 1 and 2, the Tate formulation [7,8] for studying Regime 2, or fluid dynamics [9] for studying Regime 3. Even though it is possible to find some common concepts among those theories, each penetration regime comes with its own limitations, and for this reason, it is not possible to extrapolate one theory for all penetration regimes. In this review, only the first penetration region is covered, which commonly represent events of a deformable target impacted by rigid rods. The study of rigid rod penetration has different applications. The principal fields of application are:

• The mechanics of indentation and penetration [4,10,11,12,13].
• The soil mechanics for modeling various processes (e.g., penetration) in soils and rocks [14,15,16].
• Geo-mechanical applications: areas of interpretation of in situ tests (both the cone penetrometer and the pressure meter) and also in the prediction of the behavior of piles [17,18,19].
• In mechanical characterization [20,21,22].
• Nuclear structural engineering and related areas [14,23].

The mechanics of the indentation and penetration of a rigid object into thick targets has been the subject of interest for many decades for both experimentalist and numerical analyst [13]. The physical processes happening in the rods are separate from those occurring in the target, but they are constrained by the requirement that stress needs to be continuous across the projectile-target boundary and that displacement needs to be compatible at the boundary. For rigid projectile penetration, satisfying results can be obtained with Cavity Expansion Theory (CET) and internal friction theory [5]. However, in recent years, CET has been more extensively studied than the friction theory [6,24,25,26], and new analytical models have been created around this CET approach. In general, CET consists of two different formulations: spherical or cylindrical cavity expansion methods. The first one is formulated under spherical symmetry conditions, and the second one is formulated using cylindrical coordinates. Both formulations have shown successful results to model different experimental cavity expansion effects [27]. It seems that the spherical formulation is easier to understand and simpler to model because the number of variables can be reduced due to the symmetry of the spherical problem. For this reason, many applications can be recently found using this approach [28,29,30,31].

The main objective of this review is to introduce the Spherical Cavity Expansion (SCE) theory and formulation by explaining the different mathematical features of this approach, its solution, and potential application in penetration problems. For this reason, SCE was utilized to develop an engineering model of penetration for rigid rods penetrating deformable targets. The SCE model is based on Durban et al.’s approach [32]. The paper is outlined as follows: First, a chronological literature review of the most common models used to study this kind of penetration problem is given, focusing on the cavity expansion theory. After that, the SCE theory is reviewed, using Durban et al.’s approach [32]. Next, an engineering model of penetration is proposed. This model will be based on the spherical cavity expansion theory for a strain-rate dependent material. The model is then compared and validated with some FE numerical simulations and with previous penetration results.

2. Literature Review

During the penetration process, the motion of the rigid rod is opposed by a force $F_x$ called penetration resistance. Applying Newton’s law, the motion of the projectile is governed by:

$$m_p\frac{d{V}_x}{dt}=-F_x,$$

(1)

and:

$$V_x=\frac{d{x}_p}{dt};$$

(2)

where $m_p$ and $V_x$ are the mass and instantaneous velocity of the rod and $x_p$ is the coordinate in the direction of the penetration depth. Initial conditions of this problem are: $V_x\left(t=0\right)=V_0$ and $x_p\left(t=0\right)=0$.

Common approaches are well known in penetration mechanics, which are based on making simple assumptions on the force opposing the motion of the projectile; for example, the steady penetration equations developed by Poncelet (Equation (3)) or by Resal (Equation (4)), where the penetration resistance is a function of velocity ($A_P$, $B_R$, and $C_P$ are empirical constants) [33].

$$F_{Poncelet}=A_P+C_P{V}_x^2,$$

(3)

$$F_{Resal}=B_R{V}_x+C_P{V}_x^2.$$

(4)

Although the previous models showed good agreement with most of the penetration impact events for metallic targets, they are fully empirical, and the plastic behavior of the medium was not well understood. More recently, approaches for calculating $F_x$ have assumed that the pressure during the penetration process is constant over all the cavity wall.

Forrestal et al. [34] used this approach to calculate the force for different geometries of rods, as seen in Figure 2. Due to the lack of detailed data of frictional force, it was assumed that the tangential stress on the nose (${\sigma }_t$) was proportional to the normal or radial stress (${\sigma }_r$):

$${\sigma }_r=\mu {\sigma }_t,$$

(5)

where $\mu$ is the sliding-frictional coefficient.

Therefore, assuming a spherical expansion of the cavity wall, the force of the cavity against the movement of the rod, for spherical nose rods, is given as follows:

$$F_x=\pi a^2{\int }_0^{\pi /2}{\sigma }_r\left[sin2\varphi +2\mu si{n}^2\varphi \right]d\varphi ;$$

(6)

for ogival nose rods is given as follows:

$$F_x=2\pi s^2{\int }_{{\theta }_0}^{\pi /2}{\sigma }_r\left[\left(sin\theta -\frac{s-a}{s}\right)\left(cos\theta +\mu sin\theta \right)\right]d\theta ;$$

(7)

and for conical-nose rods is given as follows:

$$F_x=\pi a^2{\sigma }_r\left(V_x\right)\left(1+\mu /tan\phi \right),$$

(8)

where all parameters involved (in the previous three equations) are shown in Figure 2.

Jones et al. [35,36] formulated an equation to calculate $F_x$ force on the wall of a general nose geometry for axis-symmetric rods, as seen in Figure 3. In this figure, $\mathrm{y}=\mathrm{y}\left(x\right)$ is the shape function of an axis-symmetric nose with $0\le x\le a$; this shape function needs to be nose-pointed $\mathrm{y}\left(0\right)=0$, and the shape function must satisfy the condition ${\mathrm{y}}^{\prime }\left(x\right)>0$ and has a base radius of $\mathrm{y}\left(a\right)=a$.

Using those geometrical conditions, F is given by:

$$F_x=2\pi \underset{0}{\overset{a}{\int }}\left(y\left(x\right)y^{\prime }\left(x\right){\sigma }_r+y\left(x\right)F_f\right)dx,$$

(9)

where ${\sigma }_r$ is the radial stress and $F_f$ is the frictional force on the cavity wall. Using Coulomb’s relation:

$${\sigma }_r=\mu F_f,$$

(10)

Equation (9) is reduced to:

$$F_x=2\pi \underset{0}{\overset{a}{\int }}{\sigma }_r\left(y{y}^{\prime }+\mu y\right)dx.$$

(11)

In the previous equations, which are related to calculating $F_x$ (Equations (6)–(8), or Equation (11)), it is necessary to know the radial stress ${\sigma }_r$ on the cavity wall. If it is assumed that the radial stress is constant over all the cavity wall and equal to the pressure on the rod nose ($p_c$), this stress can be calculated using the cavity expansion theory.

Bishop et al. [37] were the pioneers of cavity expansion pressure calculations. Bishop et al. [37] began the study of plasticity in ductile materials for quasi-static conditions. In their study, an approximate formulation was made to calculate the load required to force a cylindrical punch deep into a semi-infinite block of ductile material. Due to the fact that a full theoretical solution of the equations, determining the strains around the head of the punch, is very complex, only two problems were solved:

• Spherical Cavity Expansion (SCE): Starting with a small hollow sphere in the body of an infinite block of ductile material to determine the pressure $p_s$ that will enlarge the spherical hole indefinitely.
• Cylindrical Cavity Expansion (CCE): Starting with a cylindrical hole of infinite length to find the pressure $p_c$ that will enlarge the hole indefinitely.

Bishop et al. [37] also assumed that the true flow stress in compression ${\sigma }_i$, above the yield stress Y, is given by Equation (12), where $ϵ$ is the equivalent strain.

$${\sigma }_i=Y+f\left(ϵ\right).$$

(12)

For this material behavior, pressures $p_s$ and $p_c$ were calculated as:

$$p_s=\frac{2Y}{3}\left(1+ln\left(\frac{E}{\left(1+\nu \right)Y}\right)\right)+2\underset{a}{\overset{c}{\int }}f\left\{2ln\left(\frac{r}{{\left(r^3-a^3\right)}^{\frac{1}{3}}}\right)\right\}\frac{dr}{r},$$

(13)

$$p_c=\frac{Y}{\sqrt{3}}\left(1+ln\left(\frac{\sqrt{3}E}{2\left(\left(1+\nu \right)Y\right)}\right)\right)+\frac{2}{\sqrt{3}}\underset{a}{\overset{c}{\int }}f\left\{\frac{2}{\sqrt{3}}ln\left(\frac{r}{\sqrt{r^2-a^2}}\right)\right\}\frac{dr}{r},$$

(14)

where E and $\nu$ are the elastic modulus and the Poisson’s ratio of the material, a is the radius of the hole, and c is the radius of the plastic region around the hole. If it is assumed that the material is perfectly plastic ($f\left(ϵ\right)=0$ in Equation (12)), the integral terms in Equations (13) and (14) must be omitted.

Hopkins [3] studied the spherically-symmetric cavity formation under conditions of large elastic-plastic deformations for quasi-static and dynamic motion. The most remarkable of Hopkins’ approaches was the calculation of $p_s$ for a perfectly plastic (with a constant dynamic yield stress $Y_d$), incompressible material at dynamic conditions:

$$p_s=\frac{2Y_d}{3}\left(1+ln\left(\frac{2E}{3Y_d}\right)\right)+\rho (a\ddot{a}+\frac{3}{2}{\dot{a}}^2),$$

(15)

where $\dot{\square }$ represents the time derivative; thus, $\dot{a}$ and $\ddot{a}$ are the velocity and acceleration of the spherical radius. Furthermore, Hopkins proposed the differential equation for a perfectly plastic compressible material; however, an explicit relation was not presented.

Knowles and Jakub [38] studied spherically-symmetric motions of an elastic solid containing a spherical cavity and subjected to large strains (e.g., rubber-like materials). A perfectly plastic medium was described by a general strain energy function. Thus, two special cases were solved: Mooney and non-Hookean materials. For those special cases, values of the static stress and dynamic stress at the cavity wall for various values of the dimensionless equilibrium radius were tabulated.

Durban and Baruch [14,39] solved the non-linear problem of a spherical cavity surrounded by an infinite elasto-plastic medium and subjected to uniform and quasi-static radial loads. The governing non-linear equations were solved, in terms of closed integrals, for internal or external pressure conditions.

Using the Hopkins approach, Forrestal and Luk [40] developed an analytical model for the elastic-plastic response of a compressible material from the uniform expansion of a spherically symmetric cavity. The formulation Forrestal and Luk [40] was based on the equations of momentum and mass conservation in Eulerian coordinates. Additionally, the medium was described by a perfectly plastic behavior, and the plastic region was defined with a linear pressure-volumetric strain relation:

$$\mathrm{P}=K\left(1-\frac{{\rho }_0}{\rho }\right)=\frac{1}{3}\left({\sigma }_r+2{\sigma }_{\theta }\right),$$

(16)

and the Tresca yield criterion:

$${\sigma }_r-{\sigma }_{\theta }=Y_d,$$

(17)

where $\mathrm{P}$ is hydrostatic pressure, K is the bulk modulus, ${\rho }_0$ and $\rho$ are the densities in undeformed and deformed configurations, and ${\sigma }_r$ and ${\sigma }_{\theta }$ are the radial and hoop components of Cauchy stress; they are taken as positive in compression.

Forrestal and Luk [40] found an approximate solution of the non-linear problem assuming an incompressible material for the elastic region; although, they also proposed a full non-linear solution using a numerical method to solve the differential equations. At the end, the solution for the spherical cavity pressure can be simplified by Equation (18), where $A_s$ is expressed by Equation (19) and $C_s$ needs to be adjusted to fit cavity expansion results.

$$p_s=A_s+\rho C_s{V}_x^2,$$

(18)

$$A_s=\frac{2Y_d}{3}\left(1+ln\left(\frac{E}{\left(1-\nu \right)Y_d}\right)\right).$$

(19)

Using the same conditions, Forrestal [41] found a solution for the cylindrical cavity pressure, which can be simplified by Equation (20), where $A_c$ is expressed by Equation (21), and $C_c$ needs to be adjusted to fit cavity expansion results.

$$p_c=A_c+\rho C_c{V}_x^2,$$

(20)

$$A_c=\frac{Y_d}{2}\left(1+ln\left(\frac{6E}{\left(5-4\nu \right)Y_d}\right)\right).$$

(21)

Cylindrical and spherical cavity expansion approximations were used to develop a penetration model for a rigid rod with spherical, ogival, and conical noses that penetrated elastic-perfectly plastic targets. The results of the model were compared with penetration data and results from Lagrangian and Eulerian wave codes [34,42]. The predictions were in good agreement with the measured final depths of penetration for aluminum targets at velocities between 400 m/s and 1400 m/s.

Continuing with the Forrestal work, Luk et al. [43,44] developed the dynamic expansion of spherical cavities for elastic-plastic, incompressible and compressible materials with a power-law strain hardening relation. The modified Ludik equation was used to describe the power-law of the plastic behavior:

$$\sigma =\left\{\begin{array}{cc}Eϵ& \sigma \le Y,\\ Y{\left(\frac{Eϵ}{Y}\right)}^n& \sigma >Y,\end{array}\right.$$

(22)

where n is the strain-hardening exponent.

A closed-form solution was found for the incompressible material, and a numerical solution was proposed for the compressible case. Both solutions can be adjusted to Equation (18). However, for the incompressible problem, $A_s$ is given by:

$$A_s=\frac{2Y}{3}\left(1+I{\left(\frac{2E}{3Y}\right)}^n\right),whereI=\underset{0}{\overset{1-\frac{2E}{3Y}}{\int }}\frac{{\left(-lnx\right)}^n}{1-x}dx.$$

(23)

For the compressible case, $A_s$ and $C_s$ must be adjusted to fit cavity expansion results. Thus, a penetration model was numerically programmed, and the results were compared with penetration data. Predictions were in good agreement with the final penetration depths for aluminum targets measured at velocities between 300 m/s and 1000 m/s.

After that, Warren and Forrestal [45] solved the dynamic spherical cavity expansion problem for an elastic-plastic material with a constitutive model that included the effect of strain hardening and strain-rate sensitivity. For a state of uni-axial stress, the stress-strain relation was assumed as:

$$\sigma =\left\{\begin{array}{cc}Eϵ& \sigma \le Y_d,\\ Y{\left(\frac{Eϵ}{Y_d}\right)}^n+\alpha {\left(\frac{\dot{ϵ}}{{\dot{ϵ}}_o}\right)}^m& \sigma >Y_d,\end{array}\right.$$

(24)

$$Y_d=Y+\alpha {\left(\frac{\dot{ϵ}}{{\dot{ϵ}}_o}\right)}^m;$$

(25)

where $Y_d$ is the dynamic yield stress, m is the strain-rate sensitivity exponent, ${\dot{ϵ}}_o$ is a reference strain rate, and $\alpha$ is a curve fitting parameter.

A closed-form solution was found for an incompressible material, which is expressed by:

$$p_s=\frac{2Y}{3}\left(1+{\left(\frac{2E}{3Y}\right)}^n\underset{0}{\overset{1-\frac{2E}{3Y}}{\int }}\frac{{\left(-lnx\right)}^n}{1-x}dx\right)+\rho \left(a\ddot{a}+\frac{3}{2}{\dot{a}}^2\right)+\frac{2\alpha }{3Ym}{\left(\frac{2\dot{a}}{a{\dot{ϵ}}_o}\right)}^m.$$

(26)

The first term of Equation (26) corresponds to the strain hardening component of the solution; the second term corresponds to the inertial part of the solution; and the last term accounts for the strain-rate effects.

For the compressible problem, a numerical solution was proposed. In this case, a spherical expansion solution can be adjusted to:

$$p_s=A_s+\sqrt{\rho Y}{B}_s{V}_x+\rho C_s{V}_x^2;$$

(27)

where $A_s$, $B_s$, and $C_s$ must be adjusted to fit cavity expansion results.

To compute the structural and component responses of a projectile due to three-dimensional penetration events, Warren and Tabara [46,47] implemented the spherical cavity expansion analysis in a transient dynamics finite-element code called PRONTO 3D. The sample demonstrated good agreement between experimental and analytical results. After that, Brown et al. [48] implemented the spherical cavity expansion problem in a numerical solver as a third-party library, which can be used in explicit codes like ABAQUS.

Chen and Li [12,49] used Forrestal et al.’s approach to study the perforation of thick plates by rigid projectiles with various geometrical characteristics. They added plugging formation and movement for blunt-nosed projectiles to their model, obtaining a simple and explicit formula to predict the ballistic limit and residual velocity for the perforation of thick metallic plates, which agreed with available experimental results with satisfactory accuracy.

Teland and Moxnes [50] compared penetration results among analytical results of the cavity expansion theory and numerical simulation results, using two different nose geometries of projectiles (spherical and conical noses) and two different target materials (steel and concrete). Some agreements were achieved among the results, especially in the first stage of penetration (cratering phase). However, in the second phase (tunneling phase), the results suddenly diverge.

More recently, Buchely et al. [51,52,53] developed an engineering model based on the spherical cavity expansion theory for the penetration of a rigid-rod into an elasto-plastic, J2Mises, compressible, Cowper–Symonds strength model. The computations of Buchely and Maranon confirmed previous theoretical considerations and showed that the cavitation fields for rate-dependent materials are functions of time and cannot be modeled with steady-state self-similar approaches. However, during penetration events involving rigid projectiles, the cavity remains constant because the rigid projectile does not undergo a change, and it only moves into the field due to the movement of the projectile. Therefore, a state of stress is placed ahead of the nose of the projectile, and it moves through the material. It was shown that this state of stress is nearly steady during the penetration process at medium and low penetration velocities; therefore, self-similarity transformation can be used in this specific penetration event, even for strain-rate-sensitive materials.

The previous analyses were developed mainly for ductile targets. However, cavity expansion theory has also been extensively used in brittle materials (e.g., ceramics, glasses, and soils).

Forrestaland Lancope [54] developed a model of spherical cavity expansion for an elastic-cracked-plastic medium on the quasi-static conditions. In Forestal and Lancope’s work [54], the plastic region was described by a linear pressure-volumetric strain relation:

$${\sigma }_r-{\sigma }_{\theta }=\lambda \mathrm{P}+\frac{\left(3-\lambda \right)}{3}Y;$$

(28)

where $\lambda$ defines the pressure-dependent shear strength and is less than one ($\lambda <1$) and $\mathrm{P}$ is the hydrostatic pressure defined by Equation (16). Additionally, it was assumed that ${\sigma }_{\theta }=-T$ at the cracked-elastic interface, where T is the tensile strength of the material; and the material in the cracked region was taken as linear elastic with ${\sigma }_{\theta }=0$. Thus, the pressure to open a spherical cavity was calculated as follows:

$$R_t=\frac{Y}{\lambda }\left[\frac{2}{\alpha }{\left(\frac{c}{a}\right)}^{\alpha \lambda }-\frac{\left(3-\lambda \right)}{3}\right];$$

(29)

where the parameters $\alpha$ and $\beta$ are given as follows:

$$\alpha =\frac{6}{3+2\lambda },$$

(30)

$$\beta =\frac{9\left(1-\nu \right)}{2}\left(\frac{Y}{E}\right){\left(\frac{Y}{2T}\right)}^{1/2};$$

(31)

and the relation $\left(c/a\right)$ is given as follows:

$${\left(\frac{a}{c}\right)}^3=\beta +\frac{\beta \left(1-2\nu \right)}{\lambda E}\left[\left(3-\lambda \right)-\left(3+2\lambda \right){\left[\frac{1}{\beta }\right]}^{\frac{\alpha \lambda }{3}}\right].$$

(32)

The solution is simplified when the stress at the cavity surface $R_t<Y$; thus, the material does not reach its shear strength. For this case, the medium has a cracked region bounded by an elastic one, and the pressure to open a spherical cavity is simplified to the equation:

$$R_t=2T{\left(\frac{E}{9\left(1-\nu \right)T}\right)}^{\frac{2}{3}}.$$

(33)

This work was continued by Satapathy and Bless [55,56,57], who extensively studied the spherical cavity expansion of brittle materials. In their analysis, an elastic-cracked-comminuted medium at quasi-static and dynamic conditions was analyzed, using a two-curve pressure–shear behavior of the material (like the Johnson–Holmquist constitutive relation). At the end, closed-form solutions were proposed, and results were compared with experimental penetration data for AD995 alumina targets.

Durban and Fleck [58] also studied the spherical cavity expansion on brittle materials at quasi-static conditions, but using a Ducker–Prager constitutive yield criterion. Some special cases of materials were studied, such as the associative flow rule, elastic-perfectly plastic behavior, and Mohr–Coulomb solids. Their model was numerically solved, and different parameters of the models were studied.

Shaw [13] studied the penetration of rigid objects into semi-infinite compressible solids for penetrators with arbitrary nose geometries. Using smoothed-particle hydrodynamics simulations, he calculated the resistive force of the target through the spherical cavity expansion theory. The effect of friction at the object-target interface was also studied.

Additionally, extra work has been done to model the hyper-velocity expansion cavity, with propagation of shock waves. Cohen and Durban [59] investigated the shock waves in porous materials based on the dynamic cavity expansion theory. This work was continued by Czarnota et al. [60] for ductile metallic materials and by Rodriguez-Martinez et al. [61] for strain-hardening materials.

Summary of the Literature Review

• Penetration of undeformed projectiles into soft targets is valid at low and medium penetration velocities.
• Deceleration of the rigid projectile penetrating targets can be calculated using an analytical expression, e.g., Newton’s second law. Deceleration forces can also be calculated using analytical relations, like the one proposed by Jones et al. [35,36] and given by Equation (11)
• Cavity pressure on the target can be calculated using the spherical cavity expansion theory.

3. Spherical Cavity Expansion Formulation

3.1. Conservation Laws in Spherical Coordinates

In this subsection, the conservation laws for the expansion of a spherical cavity into a deformable target are formulated. This analysis is based on the elasto-plastic, J2 Mises, compressible model developed by Durban and Masri [32,62].

The equations of momentum and mass conservation, in spherical coordinates, are expressed by (ignoring the spin interactions):

$$\frac{\dot{\rho }}{\rho }+{\dot{ϵ}}_r+2{\dot{ϵ}}_{\theta }=0,$$

(34)

$$\frac{\partial {\sigma }_r}{\partial r}+2\frac{{\sigma }_r-{\sigma }_{\theta }}{r}=-\rho \dot{\upsilon };$$

(35)

where $\dot{\square }$ represents the time derivative; $\upsilon$ is the radial particle velocity (positive outward); $\rho$ is the material density; ${ϵ}_r$ and ${ϵ}_{\theta }$ are the radial and hoop strain; and ${\sigma }_r$, ${\sigma }_{\theta }$, and ${\sigma }_{\varphi }$ are the radial, hoop, and meridional components of Cauchy stress taken positive in compression. Additionally, ${ϵ}_{\theta }={ϵ}_{\varphi },{\sigma }_{\theta }={\sigma }_{\varphi }$ from symmetry.

3.2. Material Response

First, the effective stress ${\sigma }_e$ and plastic potential $\varrho$ are defined as follows:

$${\sigma }_e=\tau +\psi {\sigma }_h,$$

(36)

$$\varrho =\tau +\varsigma {\sigma }_h;$$

(37)

where $\tau ={\left(\frac{3}{2}\mathbb{S}··\mathbb{S}\right)}^{1/2}$ is Mises’s effective stress, ${\sigma }_h$ the hydrostatic stress, and parameters $\psi$ and $\varsigma$ represent the plastic pressure sensitivity.

Additionally, the total strain-rate $\mathbf{D}$ in the material can be expressed as follows:

$$\mathbf{D}={\mathbf{D}}^{\mathbf{E}}+{\mathbf{D}}^{\mathbf{P}};$$

(38)

where ${\mathbf{D}}^{\mathbf{E}}$ are the elastic strain-rates and ${\mathbf{D}}^{\mathbf{P}}$ are the elastic strain-rates.

Based on the elasto-plastic flow theory formulation in the hypoelastic form developed by Durban et al. [32,39,58], ${\mathbf{D}}^{\mathbf{E}}$ and ${\mathbf{D}}^{\mathbf{P}}$ are given as follows:

$${\mathbf{D}}^{\mathbf{E}}=\frac{1}{2G}\left(\check{\sigma }-\frac{\nu }{1+\nu }\check{\sigma }··\mathbb{I}\mathbb{I}\right),$$

(39)

$${\mathbf{D}}^{\mathbf{P}}=\frac{{\sigma }_e}{\varrho }{\dot{ϵ}}_{pl}\left[{\left(\frac{3}{2}\mathbb{S}··\mathbb{S}\right)}^{-1/2}\left(\frac{3}{2}\mathbb{S}\right)+\frac{1}{3}\varsigma \mathbb{I}\right];$$

(40)

where $G=E/\left(1+\nu \right)$ is the shear modulus, $\nu$ is Poisson’s ratio, $\mathbb{I}$ is an identity matrix, $\mathbb{S}=\sigma -{\sigma }_h$ are the deviatoric stresses, ${ϵ}_{pl}$ is the effective plastic-strain, and $\check{\sigma }$ is the Jaumann rate of the Cauchy stress defined as follows:

$$\check{\sigma }=\dot{\sigma }+\sigma ·w-w·\sigma ;$$

(41)

where $\sigma$ is Cauchy stress and w is the spin tensor, which can be assumed as zero for solid material flow.

Assuming that the material is hydrostatic pressure independent ($\psi =\varsigma =0$ in Equations (36) and (37)), the total strain-rate is reduced to the next equation:

$$\mathbf{D}=\frac{1+\nu }{E}\check{\sigma }-\frac{3\nu }{E}{\dot{\sigma }}_h\mathbb{I}+{\dot{ϵ}}_{pl}\left(\frac{3}{2}\frac{\mathbb{S}}{\nsim }\right).$$

(42)

Using the next two equations as definitions for $\tau$ and ${\sigma }_h$ in spherical coordinates:

$$\tau ={\sigma }_{\theta }-{\sigma }_r,$$

(43)

$${\sigma }_h=\frac{1}{3}\left({\sigma }_r+2{\sigma }_{\theta }\right);$$

(44)

the first two components of Equation (42) are expressed by the equations:

$${\dot{ϵ}}_r=\frac{1+\nu }{E}{\dot{\sigma }}_r-\frac{3\nu }{E}{\dot{\sigma }}_h+{\dot{ϵ}}_{pl}\left(\frac{3}{2\tau }\left(\sigma -{\sigma }_h\right)\right),$$

(45)

$${\dot{ϵ}}_{\theta }=\frac{1+\nu }{E}{\dot{\sigma }}_{\theta }-\frac{3\nu }{E}{\dot{\sigma }}_h+{\dot{ϵ}}_{pl}\left(\frac{3}{2\tau }\left(\sigma -{\sigma }_h\right)\right).$$

(46)

After some algebraic manipulations, Equations (45) and (46) simplify to the equations:

$${\dot{ϵ}}_r=\frac{1}{E}{\dot{\sigma }}_r-\frac{2\nu }{E}{\dot{\sigma }}_{\theta }-{\dot{ϵ}}_{pl},$$

(47)

$${\dot{ϵ}}_{\theta }=-\frac{\nu }{E}{\dot{\sigma }}_r+\frac{1-\nu }{E}{\dot{\sigma }}_{\theta }+\frac{1}{2}{\dot{ϵ}}_{pl}.$$

(48)

Equations (47) and (48) are called compatibility equations, which are used for modeling the material response.

To complete the formulation, a constitutive relation must be formulated for the material plastic flow. This relation must be a relation between the effective plastic strain and the effective stress, e.g., ${ϵ}_{pl}\left({\sigma }_e\right)$ or ${\sigma }_e\left({ϵ}_{pl}\right)$. Some authors have used simple constitutive relations as the perfectly plastic material ${\sigma }_e=Y_d$, where $Y_d$ is the dynamic yield stress of the material [40,42]; others have used strain-hardening relations ${\sigma }_e=k{ϵ}_{pl}^n$ [43,62], and also, ${\sigma }_e=f\left({ϵ}_{pl},{\dot{ϵ}}_{pl}\right)$ have been used [45].

3.3. Different Cases of the SCE Theory

3.3.1. Compressible-Dynamical Solution

Consider a compressible material under dynamical conditions. The expansion cavity problem can be divided into three regions: undisturbed, elastic, and elasto-plastic, as seen in Figure 4. In this figure, $\xi =r/a$ is the non-dimensional radial coordinate [32]; $\xi =1$ is the cavity wall; ${\xi }_i$ is the unknown elastic-plastic interface where the plastic yields occurs, which depends on the cavity expansion velocity and the material constitutive relation; and ${\xi }_w$ is the rigid-elastic wave front, which can be calculated using the elastic modulus E, Poisson’s ratio $\nu$, the reference material density (${\rho }_0$), and the cavity expansion velocity ($\dot{a}$) [32,40]:

$${\xi }_w=\frac{1}{\dot{a}}\sqrt{\frac{\left(1-\nu \right)E}{\left(1+\nu \right)\left(1-2\nu \right){\rho }_0}}.$$

(49)

In this case, it is assumed that the only independent variable is the non-dimensional radial coordinate $\xi =r/a$. Thus, the time derivative is transformed by the similarity relation [58]:

$$\frac{d}{dt}=\dot{\xi }\frac{d}{d\xi }=\left(\frac{\dot{r}}{a}-\xi \frac{\dot{a}}{a}\right)\frac{d}{d\xi }=\frac{\dot{a}}{a}\left(V-\xi \right)\frac{d}{d\xi },$$

(50)

where $V=\dot{r}/\dot{a}$ is the non-dimensional radial velocity and $\dot{a}$ the cavity expansion velocity.

Using Equation (50)’s definition, Equations (34)–(36), (47), and (48) are transformed to:

$$\left(V-\xi \right)\frac{d}{d\xi }\left[ln\left(\frac{\rho }{{\rho }_0}\right)\right]+\frac{dV}{d\xi }+2\frac{V}{\xi }=0,$$

(51)

$$\frac{d{\mathsf{\Sigma }}_r}{d\xi }-2\frac{\mathsf{\Sigma }}{\xi }={\dot{A}}^2\left(\frac{\rho }{{\rho }_0}\right)\left(V-\xi \right)\frac{dV}{d\xi },$$

(52)

$$\mathsf{\Sigma }={\mathsf{\Sigma }}_{\theta }+{\mathsf{\Sigma }}_r,$$

(53)

$$\frac{dV}{d\xi }=\left(V-\xi \right)\frac{d}{d\xi }\left({\mathsf{\Sigma }}_r-2\nu {\mathsf{\Sigma }}_{\theta }-{ϵ}_{pl}\right),$$

(54)

$$\frac{V}{\xi }=\left(V-\xi \right)\frac{d}{d\xi }\left(-\nu {\mathsf{\Sigma }}_r+\left(1-\nu \right){\mathsf{\Sigma }}_{\theta }+\frac{1}{2}{ϵ}_{pl}\right);$$

(55)

where $\left(\mathsf{\Sigma },{\mathsf{\Sigma }}_r,{\mathsf{\Sigma }}_{\theta }\right)=\left({\sigma }_e,{\sigma }_r,{\sigma }_{\theta }\right)/E$ are the non-dimensional stresses, ${\rho }_0$ is the material density in the undeformed field, and $\dot{A}$ is the non-dimensional cavity expansion velocity defined by:

$$\dot{A}=\frac{\dot{a}}{\sqrt{E/{\rho }_0}}.$$

(56)

The five governing equations (Equations (51)–(55)) determine five of the six unknown parameters of the system ($\mathsf{\Sigma }$, ${\mathsf{\Sigma }}_r$, ${\mathsf{\Sigma }}_{\theta }$, ${ϵ}_{pl}$, V, and $\rho$).

To complete the engineering model of penetration, an additional constitutive relation $\left({\sigma }_e\left({ϵ}_{pl}\right)\right)$ must be formulated for the plastic region of the material. Therefore, the spherical cavity is divided into two deformation regions, according to the non-dimensionalized radial coordinate $\xi$. The elasto-plastic region is bordered by the cavity and elastic plastic interface ${\xi }_i$; the elastic region is bordered by ${\xi }_i$ and ${\xi }_w$.

In the elasto-plastic region $1\le \xi \le {\xi }_i$, the Mises material deforms following the complete strain-displacement equations. First, the non-dimensional radial velocity is found by subtracting Equation (55) from Equation (54), integrating it and using the conditions that V and $\mathsf{\Sigma }$ vanish at $\xi \to \infty$; thus, we obtain:

$$V=\xi \left(1-{\mathrm{e}}^{-\mathsf{\Phi }}\right)\mathrm{with}\mathsf{\Phi }=\left(1+\nu \right)\mathsf{\Sigma }+\frac{3}{2}{ϵ}_{pl}.$$

(57)

Additionally, replacing Equations (54) and (55) in Equation (51), integrating it, and using $\rho ={\rho }_0$ when ${\mathsf{\Sigma }}_r={\mathsf{\Sigma }}_{\theta }=0$ [62], the following is obtained:

$$\rho ={\rho }_0{\mathrm{e}}^{-\mathsf{\Theta }}\mathrm{with}\mathsf{\Theta }=\left(1-2\nu \right)\left(3{\mathsf{\Sigma }}_r+2\mathsf{\Sigma }\right).$$

(58)

Now, replacing Equation (57) in Equation (55):

$$\left(1-2\nu \right){\mathsf{\Sigma }}_r^{{}^{\prime }}+\left(1-\nu \right){\mathsf{\Sigma }}^{{}^{\prime }}+\frac{1}{2}{ϵ}_{pl}^{{}^{\prime }}=\frac{1}{\xi }\left(1-{\mathrm{e}}^{\mathsf{\Phi }}\right);$$

(59)

and using Equations (54), (57), and (58), in the momentum conservation Equation (52), the next equation is obtained:

$${\mathsf{\Sigma }}_r^{{}^{\prime }}-2\frac{\mathsf{\Sigma }}{\xi }={\dot{A}}^2{\xi }^2\left(\left(1-2\nu \right){\mathsf{\Sigma }}_r^{{}^{\prime }}-2\nu {\mathsf{\Sigma }}_{}^{{}^{\prime }}-{ϵ}_{pl}^{{}^{\prime }}\right){\mathrm{e}}^{-\mathsf{\Theta }-2\mathsf{\Phi }};$$

(60)

where $\mathsf{\Phi }$ and $\mathsf{\Theta }$ are given by Equations (57) and (58), and superposed prime denotes differentiation with respect to $\xi$.

Equations (59) and (60) are the governing equations for the stress components ($\mathsf{\Sigma }$ and ${\mathsf{\Sigma }}_r$) in the elastic-plastic region. The boundary conditions for this system are: continuity of the radial stress and plastic yield at the elastic/plastic interface ($\xi ={\xi }_i$) and a continuous velocity on the cavity wall ($V(\xi ={\xi }_i)=1$).

In the elastic region ${\xi }_i\le \xi \le {\xi }_w$, the plastic strain and plastic strain-rate are zero (${ϵ}_{pl}\equiv 0$ and ${ϵ}_{pl}^{{}^{\prime }}\equiv 0$). Therefore, Equations (57), (59) and (60) can be simplified to:

$$V=\xi \left(1-{\mathrm{e}}^{-\mathsf{\Phi }}\right)\mathrm{with}\mathsf{\Phi }=\left(1+\nu \right)\left({\mathsf{\Sigma }}_{\theta }-{\mathsf{\Sigma }}_r\right),$$

(61)

$$-\nu {\mathsf{\Sigma }}_r^{{}^{\prime }}+\left(1-\nu \right){\mathsf{\Sigma }}_{\theta }^{{}^{\prime }}=\frac{1}{\xi }\left(1-{\mathrm{e}}^{\left(1+\nu \right)\left({\mathsf{\Sigma }}_{\theta }-{\mathsf{\Sigma }}_r\right)}\right),$$

(62)

$${\mathsf{\Sigma }}_r^{{}^{\prime }}-2\frac{{\mathsf{\Sigma }}_{\theta }-{\mathsf{\Sigma }}_r}{\xi }={\dot{A}}^2{\xi }^2\left({\mathsf{\Sigma }}_r^{{}^{\prime }}-2\nu {\mathsf{\Sigma }}_{\theta }^{{}^{\prime }}\right){\mathrm{e}}^{\left(1+4\nu \right){\mathsf{\Sigma }}_r-2\left(2-\nu \right){\mathsf{\Sigma }}_{\theta }};$$

(63)

where Equation (53) was used.

One approximation for Equations (62) and Equation (63) is to assume that as $\xi$ increases, ${\xi }^3>>1$ (linearization of the system). Thus, these equations are reduced to:

$$\nu {\mathsf{\Sigma }}_r^{{}^{\prime }}+\left(1-\nu \right){\mathsf{\Sigma }}_{\theta }^{{}^{\prime }}=-\frac{\left(1+\nu \right)}{\xi }\left({\mathsf{\Sigma }}_{\theta }-{\mathsf{\Sigma }}_r\right),$$

(64)

$${\mathsf{\Sigma }}_r^{{}^{\prime }}-2\frac{{\mathsf{\Sigma }}_{\theta }-{\mathsf{\Sigma }}_r}{\xi }={\dot{A}}^2{\xi }^2\left({\mathsf{\Sigma }}_r^{{}^{\prime }}-2\nu {\mathsf{\Sigma }}_{\theta }^{{}^{\prime }}\right).$$

(65)

Equations (62) and (63) (or Equations (64) and (65), if the linear approximation is assumed) are the governing equations for the stress components (${\xi }_r$ and ${\xi }_{\theta }$) in the elastic region. The boundary conditions for this system are: stress-free conditions (${\mathsf{\Sigma }}_r={\mathsf{\Sigma }}_{\theta }=0$) as $\xi \to \infty$ and continuity of the radial stress and plastic yield at the elastic–plastic interface ($\xi ={\xi }_i$).

3.3.2. Incompressible-Dynamical Solution

For an incompressible material, the same conditions of the compressible-dynamical solution can be applied; however, some simplifications could be made to obtain an analytical solution.

First, the mass conservation Equation (51) is reduced to:

$$\frac{dV}{d\xi }+2\frac{V}{\xi }=0.$$

(66)

This equation has the known solution (using the wall condition $V\left(\xi =1\right)=1$):

$$V=\frac{1}{{\xi }^2}.$$

(67)

Using this solution, momentum conservation (Equation (52)) is reduced to:

$$\frac{d{\mathsf{\Sigma }}_r}{d\xi }-2\frac{\mathsf{\Sigma }}{\xi }=2{\dot{A}}^2\left(\frac{{\xi }^3-1}{{\xi }^5}\right),$$

(68)

and the compatibility equations (Equations (54) and (55)) are combined and reduced to (using $\nu =1/2$):

$$\frac{d\mathsf{\Sigma }}{d\xi }+\frac{d{ϵ}_{pl}}{d\xi }=\frac{2}{\xi -{\xi }^4}.$$

(69)

Equations (68) and (69) are the governing equations for the stress components ($\mathsf{\Sigma }$ and ${\mathsf{\Sigma }}_r$) and represent the behavior of the material in the elasto-plastic region ($1\le \xi \le {\xi }_i$), as seen in Figure 4.

In the elastic field ${\xi }_i\le \xi \le \infty$, where plastic strain is zero $\left({ϵ}_{pl}\equiv 0\right)$, the solutions of Equations (68) and (69) using a stress-free condition, as $\xi \to \infty$, are:

$$\mathsf{\Sigma }=\frac{2}{3}ln\left(\frac{{\xi }^3}{{\xi }^3-1}\right),$$

(70)

$${\mathsf{\Sigma }}_r=-\frac{4}{3}\underset{\xi }{\overset{\infty }{\int }}ln\left(\frac{x^3}{x^3-1}\right)\frac{dx}{x}-2{\dot{A}}^2\left(\frac{1}{\xi }-\frac{1}{4{\xi }^4}\right).$$

(71)

One approximation for Equations (70) and (71) is assuming that as $\xi$ increases, ${\xi }^3\gg 1$. Thus, the previous equations are reduced to:

$$\mathsf{\Sigma }=\frac{2}{3{\xi }^3},$$

(72)

$${\mathsf{\Sigma }}_r=-\frac{4}{9{\xi }^3}-2{\dot{A}}^2\left(\frac{1}{\xi }-\frac{1}{4{\xi }^4}\right).$$

(73)

This system can be analytically solved, e.g., assuming an elastic perfectly plastic material ($\mathsf{\Sigma }=Y_d/E={\mathsf{\Sigma }}_0$). Therefore, Equation (68) is transformed to:

$$\frac{d{\mathsf{\Sigma }}_r}{d\xi }-2\frac{{\mathsf{\Sigma }}_0}{\xi }=2{\dot{A}}^2\left(\frac{{\xi }^3-1}{{\xi }^5}\right),$$

(74)

the analytical solution of which, using the boundary condition ${\mathsf{\Sigma }}_r\left(\xi =1\right)=-p_s/E$, is given by:

$${\mathsf{\Sigma }}_r\left(\xi \right)=-\frac{p_s}{E}+2{\mathsf{\Sigma }}_{0}ln\left(\xi \right)+\frac{2{\dot{A}}^2}{\xi }\left(\frac{1}{4{\xi }^3}-1\right)+\frac{3{\dot{A}}^2}{2};$$

(75)

where $p_s$ is the pressure at the cavity wall.

The solution of Equation (75) is obtained using Equation (73) as a boundary condition, and the condition $\mathsf{\Sigma }\left(\xi ={\xi }_i\right)={\mathsf{\Sigma }}_0$:

$$p_s=\frac{2}{3}{Y}_d\left(1+ln\left(\frac{2}{3{\mathsf{\Sigma }}_0}\right)\right)+\frac{3\rho }{2}{\dot{a}}^2.$$

(76)

This is the same result found by Hopkins [3] for the same perfectly plastic solid.

4. Formulation of the Engineering Model of Penetration

In the previous section, the equations of the SCE theory were formulated and analyzed for different scenarios. In this section, an engineering model of penetration is going to be formulated using SCE theory and common plasticity models.

4.1. Projectile Deceleration

Consider the penetration of an axis-symmetrical rigid-rod of mass $m_p$ into a low-strength target, as shown in Figure 5. In this figure, $\mathrm{y}\left(x\right)$ is the shape function of an axi-symmetric nose of the rod; $L_0$ is the cylindrical rod length; a is the rod radius; $V_x$ is the rod velocity; $\beta$ is the meridional angle; and $x_p$ is the rod penetration depth into the target.

According to Newton’s second law, the motion equation of the rigid-rod is given by:

$$m_p\frac{d{V}_x}{dt}=m_p\frac{d^2{x}_p}{d{t}^2}=-F_x;$$

(77)

where $F_x$ is the resisting force produced by the target material onto the projectile, which can be calculated by a balance of forces on the spherical cavity wall in the following way [35,36]:

$$F_x=2\pi {\int }_0^{x_p}\left(\mathrm{y}{\mathrm{y}}^{\prime }{p}_s+\mathrm{y}f\right)dx;$$

(78)

where $p_s$ is the cavity pressure and f is the frictional force. Assuming proportionality between the radial and tangential stresses on the cavity wall, the frictional force can be calculated using Coulomb’s relation [34]:

$$f=\mu p_s;$$

(79)

where $\mu$ is the dynamic sliding-friction coefficient, which was assumed negligible beyond the nose at $x_p\ge a$.

In this analysis, the cavity pressure $p_s$ must be determined to solve the penetration problem. In the next subsection, $p_s$ is predicted using the spherical cavity expansion theory.

4.2. SCE Theory in Penetration Problems

In a previous subsection (see Section 3.3.1), it was shown that the dynamic behavior of a compressible elasto-plastic material can be predicted by solving the non-linear differential equation system comprised by Equations (51)–(55). These five governing equations determine five of the six unknown parameters of the system ($\mathsf{\Sigma }$, ${\mathsf{\Sigma }}_r$, ${\mathsf{\Sigma }}_{\theta }$, V, ${ϵ}_{pl}$, and $\rho$). To complete the formulation, an additional constitutive relation ${\sigma }_e\left({ϵ}_{pl}\right)$ must be formulated for the plastic flow of the material. Therefore, the spherical cavity is divided into two regions (as seen in Figure 4).

In the elastic region ${\xi }_i\le \xi \le {\xi }_w$, the plastic strain and plastic strain-rate are zero (${ϵ}_{pl}\equiv 0$ and ${ϵ}_{pl}^{{}^{\prime }}\equiv 0$). Therefore, the system is simplified into three equations: Equations (61)–(63), or to Equations (64) and (65), if linear approximations are assumed.

In the elasto-plastic region $1\le \xi \le {\xi }_i$, the system can be combined to reduce one variable and to simplify it into four equations: Equations (57)–(60). Combining those equations with the constitutive relation for the plastic flow, the solution of the cavity problem is completed.

4.3. Model Solution

The solution of the penetration model consists of three different sub-processes. The first sub-process is the solution of the cavity expansion problem for the target material (see Figure 6a). The second sub-process is related to expressing the cavity expansion as a function of the expansion velocity (see Figure 6b). Finally, in the third sub-process, the projectile deceleration equation is solved (see Figure 7).

When the cavity expansion problem is solved, the cavity pressure $p_s$ is found for the target material. The inputs for this problem are the material properties ($\rho$, E, $\nu$, and $Y_0$), the constitutive relation ${\sigma }_e({ϵ}_{pl},{\dot{ϵ}}_{pl})$), and the geometrical parameters (a, $\dot{a}$). Notice that the problem is solved for a given cavity radius a and expansion velocity $\dot{a}$. The cavity radius is constant during the event because the projectile is not changing its nose shape (it is rigid), and the expansion velocity is assumed as constant for the solution. First, a guess value ${\xi }_i$ is selected, and the elastic region (${\xi }_i\le \xi \le {\xi }_w$) is solved, using the boundary conditions $\mathsf{\Sigma }\left({\xi }_w\right)={\mathsf{\Sigma }}_r\left({\xi }_w\right)=0$ and the yield criterion $\mathsf{\Sigma }\left({\xi }_i\right)=Y_0$. After the elastic region is solved, ${\mathsf{\Sigma }}_r\left({\xi }_i\right)$ and $V\left({\xi }_i\right)$ are known. These values are taken as boundary conditions for the elasto-plastic region, assuming continuity in the elastic-plastic interface. After that, the elasto-plastic region ($1\le \xi \le {\xi }_i$) is solved, using a constitutive relation for the plastic flow of the material (e.g., Cowper–Symmonds and Johnson–Cook relations, or any other constitutive relation ${ϵ}_{pl}^{{}^{\prime }}\left(\mathsf{\Sigma }\right)$). When the cavity condition at the cavity surface $V\left(\xi =1\right)=1$ is satisfied, the right values of ${\mathsf{\Sigma }}_r\left(\xi =1\right)$ and $\mathsf{\Sigma }\left(\xi =1\right)$ are obtained, corresponding to the chosen value of ${\xi }_i$. Finally, the value of $p_s$ is found using the boundary solution:

$$P_s=p_s/E={\mathsf{\Sigma }}_r\left(\xi =1\right);$$

(80)

where $P_s$ is the non-dimensional cavity pressure.

In the second sub-process, the cavity expansion solution must be solved for a range of non-dimensional expansion velocities $\dot{a}$, as seen in the loop process of diagram in Figure 6b. Usually, this range of $\dot{a}$ is selected specifically for the impact velocities covered during some penetration experiments. After the loop process, and according to some previous available studies that suggest a dependence of $p_s$ on powers of $\dot{a}$ [45,62], the cavity pressure $p_s$ can be approximated using a power function of $\dot{a}$,

$$p_s=a_c+c_c{\dot{a}}^2+d_c{\dot{a}}^3+f_c{\dot{a}}^{n_c},\mathrm{or}$$

(81)

$$P_s=A_c+C_c{\dot{A}}^2+D_c{\dot{A}}^3+F_c{\dot{A}}^{n_c};$$

(82)

where $a_c$ is a coefficient related to the quasi-static solution, because it does not depend on the velocity; $c_c$ is the coefficient related to the inertial movements; and $d_c$, $f_c$, and $n_c$ are the coefficients and exponent related to the strain hardening and the strain-rate contribution. Additionally, $P_s=p_s/E$ is the non-dimensional cavity pressure; therefore, $\left(A_c,C_c,D_c,F_c\right)=\left(a_c,c_c,d_c,f_c\right)/E$ are the non-dimensional coefficients.

In the third sub-process, the deceleration of the projectile and the final penetration depth $d_p$ are predicted. Therefore, Equation (81) is replaced in Equation (78), but first, the expansion velocity $\dot{a}$ and the projectile velocity $V_x$ must be related by:

$$\dot{a}=V_x\cos \varphi ;$$

(83)

where $\varphi$ and $x_p$ are geometrically coupled by:

$$\cos \varphi =1-\frac{x_p}{a}.$$

(84)

The engineering model of penetration is completed with the inputs of the penetration process ($m_p$, a, $V_0$, and $\mu$), and solving Equation (77). Then, the position and deceleration of the projectile can be found during the penetration event (Figure 7).

5. Results and Discussion

In this section, the solution methodology for the engineering model of penetration is compared with previous analytical solutions. In particular, the solution of the SCE problem was compared for two kinds of material constitutive relations: elastic-perfectly plastic, strain-hardening material and power-law relations. Additionally, the engineering model methodology was solved for two different compressible, elasto-plastic aluminum alloys. The particular constitutive relations used and studied were the modified Ludwik relation.

5.1. Comparison to the Forrestal et al. SCE Solution

Forrestal et al. [40] formulated the solution of the SCE field for an elastic-perfectly plastic material, using the constitutive equation:

$$\mathsf{\Sigma }=Y_d/E=\overline{Y};$$

(85)

where $Y_d$ is the dynamic yield stress of the material and E is the material elastic modulus.

Forrestal et al. [40] solved this system for a 6061-T6 aluminum, the parameters of which are shown in

Table 1. Parameters used in the Forrestal Spherical Cavity Expansion (SCE) model [40].

	$\mathit{\rho }$ (kg/m${}^3$)	E (GPa)	$\mathit{\nu }$ (-)	${\mathit{Y}}_{\mathit{d}}$ (MPa)
6061-T6 Al	2700	68.9	0.33	300

, and the pressure at the spherical cavity wall ($P_s=p_s/E$) was plotted versus a non-dimensional cavity expansion velocity, as seen in Figure 8.

As a comparison, the proposed solution methodology was adjusted to these material conditions. First, it is noticed from Equation (85) that ${\mathsf{\Sigma }}^{{}^{\prime }}=0$. Replacing it in Equations (59) and (60), we found that:

$$\left(1-2\nu \right){\mathsf{\Sigma }}_r^{{}^{\prime }}+\frac{1}{2}{ϵ}_{pl}^{{}^{\prime }}=\frac{1}{\xi }\left(1-e^{\mathsf{\Phi }}\right),$$

(86)

$${\mathsf{\Sigma }}_r^{{}^{\prime }}-\frac{2}{\xi }\overline{Y}={\dot{A}}^2{\xi }^2\left[\left(1-2\nu \right){\mathsf{\Sigma }}_r^{{}^{\prime }}+\frac{1}{2}{ϵ}_{pl}^{{}^{\prime }}\right]e^{-\mathsf{\Theta }-2\mathsf{\Phi }};$$

(87)

where $\mathsf{\Phi }=\left(1+\nu \right){\mathsf{\Sigma }}_r+\frac{3}{2}{ϵ}_{pl}$ and $\mathsf{\Theta }=\left(1-2\nu \right)\left(3{\mathsf{\Sigma }}_r+2\overline{Y}\right).$

Solving for ${ϵ}_{pl}^{{}^{\prime }}$ in Equation (86):

$${ϵ}_{pl}^{{}^{\prime }}=-2\left(1-2\nu \right){\mathsf{\Sigma }}_r^{{}^{\prime }}+\frac{2}{\xi }\left(1-e^{\mathsf{\Phi }}\right).$$

(88)

Using Equation (88) as the constitutive equation, the first two processes of the proposed methodology (Section 4.3) were applied to this problem (see Figure 6).

Figure 8 shows the results of $p_s$ from the engineering model (asterisk markers) at different non-dimensional cavity expansion velocities. Notice that the solution of the engineering model agreed well with the fully-nonlinear solution proposed by Forrestal et al. Therefore, it showed that the engineering model could be used to model this kind of behavior (elastic-perfectly plastic materials).

5.2. Comparison to the Masri and Durban SCE Solution

Masri and Durban [62] formulated the solution of SCE for an elastic, strain-hardening material, using the power-law constitutive equation:

$${ϵ}_{pl}=k{\mathsf{\Sigma }}^n.$$

(89)

They solved this system for different kinds of materials, the parameters of which are shown in

Table 2. Parameters used in Masri and Durban’s SCE model [62].

Material	$\mathit{\rho }$ (kg/m${}^3$)	E (GPa)	$\mathit{\nu }$ (-)	k (-)	n (-)
Titanium B120VCA	4400	106	0.333	2.4 × ${10}^{29}$	16.5
Stainless steel	7800	206	0.3	5.78 × ${10}^4$	3
Steel D6AC	7800	213	0.27	2.52 × ${10}^{55}$	28
Aluminum7075-T6	2700	72.4	0.32	3.94 × ${10}^{21}$	10.9

. Furthermore, the pressure at the spherical cavity wall ($P_s$) was plotted versus a non-dimensional cavity expansion velocity, as seen in Figure 9.

As a comparison, the proposed methodology of the solution was adjusted to these material conditions. In this case, the derivative of Equation (89) with respect to $\xi$ was first calculated:

$${ϵ}_{pl}^{{}^{\prime }}=kn{\mathsf{\Sigma }}^{n-1}{\mathsf{\Sigma }}^{{}^{\prime }}.$$

(90)

Using Equation (90) as the constitutive equation, the first two processes of the proposed methodology (Section 4.3) were applied to this problem (see Figure 6).

Figure 9 shows the results of $p_s$ from the engineering model (markers on the plot) for all the materials (Table 2) and at different non-dimensional cavity expansion velocities. Again, notice that the solution of the engineering model agreed well with the solution proposed by Masri and Durban for all the materials. Therefore, it showed that the engineering model could be used to model this kind of behavior (elastic, strain-hardening materials).

5.3. Results of SCE Model for 7075-T6 Aluminum Alloys and Comparison to FE Models

Consider a rigid sphere penetrating a soft target, as seen in Figure 10. In this figure, a is the radius of the projectile with mass $m_p$ and $\beta$ and $x_p$ are the meridional angle and penetration depth target, respectively.

For the hemispherical nose of the sphere, $L=a$ and $\mathrm{y}\left(x\right)=\sqrt{a^2-{\left(a-x\right)}^2}$; thus, ignoring frictional forces, Equation (78) is transformed to:

$$F=2\pi \underset{0}{\overset{x_p}{\int }}{p}_s\left(a-x\right)dx,$$

(91)

Therefore, to solve the penetration problem, the pressure at the cavity wall $p_s$ must be determined, using the spherical cavity expansion theory.

First, it is necessary to formulate a constitutive relation for the target material; for this example, the modified Ludwik relation is used:

$${\sigma }_e=Y+B{ϵ}_{pl}^n.$$

(92)

In this specific case, an aluminum 7075-T6 target material was chosen, the properties of which are shown in

Table 3. Material properties for 7075-T6 aluminum.

$\mathit{\rho }$ (kg/m${}^3$)	E (GPa)	$\mathit{\nu }$ (-)	Y (MPa)	B (MPa)	n (−)
2700	72.4	0.32	601	765	0.09174

.

Solving Equation (92) for ${ϵ}_{pl}$:

$${ϵ}_{pl}={\left(\frac{1}{\overline{B}}\right)}^{1/n}{\left(\mathsf{\Sigma }-\overline{A}\right)}^{1/n};$$

(93)

where $\left(\mathsf{\Sigma },\overline{A},\overline{B}\right)=\left({\sigma }_e,A,B\right)/E$ are the non-dimensional stresses and differentiating Equation (93) with respect to $\xi$:

$${ϵ}_{pl}^{{}^{\prime }}=\frac{1}{n}{\left(\frac{1}{\overline{B}}\right)}^{1/n}{\left(\mathsf{\Sigma }-\overline{A}\right)}^{\left(1-n\right)/n}{\mathsf{\Sigma }}^{{}^{\prime }},$$

(94)

Using Equation (94) as the constitutive equation, the solution of the SCE model was performed, applying the methodology explained in Section 4.3 and structured in Figure 6. After carrying out these two processes, non-dimensional pressure at the cavity wall $P_s$ was calculated as ${\mathsf{\Sigma }}_r\left(\xi =1\right)=P_s=p_s/E$. Figure 11 (circle markers) shows the $P_s$ results at different expansion velocities $\dot{A}$, for 7075-T6 aluminum.

After that, the $P_s$ results were adjusted using Equation (81):

$$P_s=\frac{p_s}{E}=0.0496+1.144{\dot{A}}^2-06202{\dot{A}}^3.$$

(95)

Notice that the last term of Equation (81) is not placed in Equation (95) because the modified Ludwik constitutive relation is not strain-rate dependent. The prediction of Equation (95) is also shown in Figure 11 (solid line). Notice that polynomial approximation agrees well with the exact cavity pressure solution (circle markers). Using Equation (95), the engineering model of penetration was completed replacing Equation (95) into Equation (91), as shown in Figure 7. The results of the engineering model are shown in Figure 12 as lines (solid line at $V_0=250$ m/s and dotted line at $V_0=500$ m/s).

To compare and verify the predictions of the engineering model, a 2D axis-symmetric computational FE model of the penetration problem was also solved using ANSYS/AUTODYN. A Lagrangian mesh was used and optimized in order to minimize the hourglass energy and contact energy among elements to less than $5\%$ of the total energy balance [63]. The target was modeled as a semi-infinite aluminum solid, and the sphere was modeled using a M300 steel alloy with $\varphi =4$ mm. The target was modeled using a rectangular grid of 200 cells × 200 cells. A grade zoning was applied in the contact (front) border, using $0.1$ mm linear increments. On the outer borders, a transmit boundary condition was applied in order to simulate the material as semi-infinite. The sphere was created using a $0.2$ mm grid mesh. Friction between bodies was neglected.

The material parameters used in the FE simulation are shown in Table 3 and

Table 4. Material properties for M300 maraging steel.

$\mathit{\rho }$ (kg/m${}^3$)	E (GPa)	$\mathit{\nu }$ (-)
2700	189	0.3

for the Al target and the steel sphere, respectively. For the target, a linear (elastic) equation of state (EOS) and the strength hardening equation (see Equation (92)) were used. For the sphere, only a linear (elastic) EOS was selected, assuming that the sphere was going to behave as an elastic body without permanent deformation during penetration. Damage models were not included in the simulation

Figure 13 shows a sequence of the FE simulation using a sphere at $V_0=500$ m/s. Notice that the propagation of the stresses showed a spherical pattern, which confirmed the assumption of a spherical symmetry field around the cavity expansion.

A comparison between the engineering model of penetration and FE simulations is shown in Figure 12 at two different initial sphere velocities (250 m/s and 500 m/s). In this figure, markers represent the results of the computational simulations, whereas lines represent the prediction of the engineering model. Good agreement between the engineering model and the computational model, both in penetration profile and final penetration, can be observed, which means that the engineering model could be used to predict the penetration depth in this kind of material.

6. Conclusions

In this review, an engineering model of penetration to predict the final penetration depth when a soft target was impacted by a rigid projectile was formulated. The proposed model was based on the spherical cavity expansion theory using the Durban et al. approach. The input given by the present model was the fact that common constitutive equation models were used, which reduced the experimentation, particularly the dynamic characterization of the materials because the parameters of those models could be found in many references for different kinds of materials.

First, the cavity expansion problem was solved and verified using previous available data from other research:

• Verified using elastic, perfectly plastic, 6061-T6 aluminum.
• Verified using four elastic, strain-dependent materials: titanium B120VCA, stainless steel, D6AC steel, and 7076-T6 aluminum.
• Verified and validated by FE numerical simulations and experimental data using semi-infinite, strain-dependent, 7075-T6 aluminum targets in the range of 250 m/s to 500 m/s.

According to the results obtained, we can conclude that the engineering model can predict the projectile deceleration and penetration depth into the semi-infinite and finite targets.

Some limitations of the model included the fact that the deceleration function of the projectile must be formulated using some simplified models, where the frictional forces were unknown. Additionally, the model did not take into account a failure mechanism into the target. Therefore, some errors could occur, especially when finite targets were modeled, which also limited the use of the model mainly for ductile materials.

Finally, the engineering model of penetration could be used as an inverse problem to characterize the dynamic properties of the materials. Although, only simple relations were used, it was a useful tool due to its simplicity and faster solution of the inverse problem.