A combined theoretical–semiempirical penetration model of ballistic penetration of thick section composites

Abstract

A combined theoretical–semiempirical penetration model of ballistic penetration of thick section composites has been developed. The penetration model is developed by breaking down the penetration event into two major phases, that is, (i) the short-time phase I shock compression and (ii) the long-time phase II dynamic penetration. Phase I shock compression considers a one-dimensional stress wave propagation theory and the rate effects on material properties. On the other hand, the phase II dynamic penetration phase incorporates a dynamic increase factor to a semiempirical quasi-static penetration resistance curve. The combined theoretical–semiempirical penetration model is compared with finite element analyses of the same problem and good correlations are obtained.

Introduction

Terminal ballistics^1–6 of finite thickness targets seek theoretical solutions for (i) the ballistic limit velocity ( $V_{\mathrm{B}\mathrm{L}}$ ) as a function of the properties of the projectile–target pair, (ii) a penetration equation correlating the residual velocity of the projectile as a function of the impact velocity, (iii) the penetration resistance (PR) force acting on the projectile as a function of time, and (iv) the time history of projectile erosion. With known PR force and the rate of projectile erosion, if any, the dynamics of the projectile, that is instantaneous rigid body velocity and location of the projectile, can be determined. In Ref. 1, a good summary on terminal ballistic of metals, ceramics, concrete, and soil is presented; however, little is mentioned about composites. In Refs. 2–4, theory of penetration for metals has been outlined. While low velocity impact on composites is a well-studied subject,^7–9 penetration and perforation of thick section composites are not and is the motivation for the present work.

Ballistic limit velocity of a finite thickness target–projectile pair is mostly determined experimentally. In a semiempirical equation in the CGS (centimeter-gram-second) unit formulated by DeMarre and later modified by Lambert for metals, ¹ the parameters should be determined experimentally for each projectile–target pair. A theoretical model was proposed by Forrestal et al. ¹⁰ which can predict the residual velocity of a long rod penetrating the aluminum targets within 15% accuracy. Florence ¹¹ developed a simplified model for ceramic–aluminum targets impacted by a short cylindrical projectile. A more rigorous model is proposed by Walker and Anderson ¹² which require solving complex equations and writing a computer program and later Zaera and Sanchez-Galvez ¹³ used Tate’s equation to develop a penetration model for ceramic–metal targets.

For impact velocities higher than the ballistic limit, Recht and Ipson ¹⁴ and Lambert and Jonas ¹⁵ equations are usually used to predict the residual velocity of the projectile.

Both the Recht and Ipson and Lambert and Jonas models are mostly used in curve fitting the ballistic experiments. Both simple and complex models have been used to predict the ballistic limit for metals, two-component ceramic–composite armor, fabrics and soft laminate^16,17; however, the focus of the present research is to develop penetration models for thick section composites.

While multiscale computational analyses (both Lagrangian and Eulerian) can be used to solve the governing equations of motion and model the physics of penetration under high- and hyper-velocity impacts by modeling the constitutive behavior of the materials, equation of state (EoS), failure mechanisms, and so on^18–32; theoretical solutions are necessary to reduce the computational time and to perform parametric calculations.

For thin section composites, some theoretical models are available^33–40 to predict the ballistic limit and projectile dynamics, but for thick section composites under ballistic impact, we could not find any suitable model that can predict the ballistic limit and projectile dynamics.^41–45 A recent theoretical model^39,40 assumes the penetration of thin composite laminates and uses the conservation of energy and momentum to balance the energy dissipated by different composite damage mechanisms. This model predicts the ballistic limit, and projectile dynamics with some accuracy, however, cannot be used for thick section composites, in which the through-thickness stress wave propagation in the early stage of impact needs to be considered in predicting the initial impact contact force and projectile dynamics. In this regard, we have developed a punch shear–based penetration model ⁴⁴ for ballistic penetration, and a detailed three-dimensional (3D) finite element analysis (FEA) of impact and penetration of thick composites. ⁴⁶ While the punch shear–based penetration model ⁴⁴ can predict the ballistic limit but not the projectile dynamics, the 3D FEA with progressive composite damage material model MAT162 ⁴⁶ can model the correct stress wave propagation and different energy-dissipating damage mechanisms of composites ⁴⁷ such as matrix cracking, delamination and fiber fracture under combined tension shear, compression shear, punch shear, and crush. Our 3D FEA work ⁴⁶ has predicted the projectile dynamics and the PR force as a function of time for a wide range of impact velocities. From these two experimental and FEA work^44,46 and also from other literature cited, ballistic impact and penetration of a composite laminate under a right circular cylinder projectile can be described in two major phases, that is (i) the short-time phase I—shock compression and (ii) the long-time phase II—penetration and perforation. In our experimental work, we have found that the ballistic limit velocity of a 13.2-mm-thick plain weave (PW) S-2 glass/SC15 composite under the impact of a 0.50-cal right circular cylinder of mass 13.8 g is $V_{\mathrm{B}\mathrm{L}}=367\mathrm{m}/\mathrm{s}$ and the ballistic limit energy is $E_{\mathrm{B}\mathrm{L}}=930\mathrm{J}$ . ⁴⁴ We have also developed a quasi-static punch shear test (QS-PST) methodology ⁴⁴ to quantify the energy dissipation capacity under a quasi-static envelope curve (Figure 1). This study on PW S-2 glass/SC15 epoxy concluded that 67% of the ballistic penetration energy (area under the hydrostatic (HS)-envelope curve) could be attributed to material damage. This indicated the relative importance of the dynamic and strain rate behavior on PR of this material during phase I and phase II processes.

Figure 1.

Quasi-Static Penetration Resistance (QS-PR) curve as a function of punch displacement for plain weave (PW) S-2 Glass/SC15 composites. QS-PR curve is obtained from a series of quasi-static punch shear tests (QS-PSTs) on composite specimens of same thickness as a function of circular support span to punch diameter ratio, $SPR=D_S/D_P$ . Adapted from Mandeley Publishers, 2008. ⁴⁴

Through these works, we have identified that there is no theoretical model for composites that can predict the impact force during the impact contact or shock phase, and this is the subject matter of the present article. We have obtained closed form theoretical solution for the shock phase and are also developing theoretical models for the subsequent phases. Since a theoretical model for the later phases are not available yet, we have adopted our quasi-static PR (QS-PR) curve (Figure 1) with a dynamic increase factor (DIF) to develop the present combined theoretical–semi-empirical model of ballistic penetration, which can predict the ballistic limit of a thick section composite plate and can also predict the projectile dynamics at the ballistic limit.

A combined theoretical–semiempirical penetration model of ballistic penetration

Previous finite element modeling ⁴⁶ has identified that the phase I shock stress increases with an increase in impact velocity, and the amplitude of the phase II PR force (with an exponential decay behavior) also increases with an increase in impact velocity. In order to capture the projectile dynamics, stress wave propagation, and PR force as a function of impact velocity; a theoretical one-dimensional (1D) shock compression model for phase I and a semiempirical model for phase II are proposed and are presented in this section.

1D phase I shock compression model of projectile impact on composites

The impact contact force on the projectile as a function of initial impact velocity is the primary unknown quantity during phase I. We will consider stress wave propagation in both the projectile and the composite plate upon impact and also consider a rate-dependent elastic–plastic strain hardening constitutive equation to describe the through-thickness compression behavior of the composite and derive an analytical expression for the PR force. Knowing the PR force as a function of time, the projectile dynamics during penetration will also be evaluated.

Following the rigid projectile models of Galanov et al. ⁴⁸ and Lopatnikov et al. ⁴⁹ for 1D failure/shock wave propagation, we will consider elastic stress wave propagation in a right circular cylinder projectile of mass $m_{\mathrm{P}}$ , cross-sectional area $A_P$ , length $H_P$ , and density ${\mathit{\rho }}_P$ ; impacting on a composite plate of thickness $H_C$ , and density ${\mathit{\rho }}_C$ with zero obliquity and with an impact velocity of $V_I$ (Figure 2(A)). The location and velocity of the projectile composite contact interface (PCCI) are denoted byZ, and $\dot{Z}$ , respectively. The origin of the through-thickness direction, Z (or 3), is located at the PCCI at time $t=0$ . For simplicity, we will assume 1D stress wave propagation throughout the thickness of the composite plate. The stress–strain behavior of the composite under uniaxial strain impact loading will be assumed to be elastic–plastic with strain hardening (Figure 2(B)), where the elastic and plastic modulus are defined as $E_{3e}^{/}=C_{33}$ (Appendix A) and $E_{3p}^{/}$ (will be determined experimentally in the section on Results and discussion ). Upon impact, shock stress will be generated in both the projectile and the composite laminate. If the shock stress in the composite is higher than that in the plain strain yield stress ${\mathit{\sigma }}_{3,C}^{\ast }$ , an elastic and a plastic stress wave front will be generated in the composite laminate and propagate with the wave-front velocities, ${\dot{Z}}_e=c_{3L}^{\ast }$ and ${\dot{Z}}_p=c_{3p}^{\ast \ast }+u_C^{\ast }$ , respectively, with respect to the laboratory frame of coordinate system. Where, $c_{3L}^{\ast }=\sqrt{E_{3e}^{/}/{\mathit{\rho }}_C}$ is the longitudinal wave velocity in the through-thickness direction, and $c_{3p}^{\ast \ast }$ is the unknown plastic-front velocity with respect to the immovable material, which will be calculated using the stress–density relationship of the material and the plastic wave-front conditions. The density and the particle velocity of the composite in the elastic zone is denoted by ${\mathit{\rho }}_C^{\ast }$ and $u_C^{\ast }$ ; and in the plastic zone by ${\mathit{\rho }}_C^{\ast \ast }$ and $u_C^{\ast \ast }$ .

Figure 2.

One-dimensional (1D) shock compression model of projectile impact on a composite plate. (A) 1D shock compression model. (B) Compressive engineering stress-strain model for composite under plain strain impact loading.

Similarly, considering the elastic wave propagation in the projectile, an elastic stress wave front will propagate in the projectile with the wave-front speed $c_{0,P}^{\ast }=\sqrt{E_P/{\mathit{\rho }}_P}$ (bar velocity), in the opposite direction to the projectile velocity with the particle velocity of $u_P^{\ast }$ . The stress in the composite elastic zone, composite plastic zone, and in the projectile elastic zone are denoted by ${\mathit{\sigma }}_{3,C}^{\ast }$ , ${\mathit{\sigma }}_{3,C}^{\ast \ast }$ , and ${\mathit{\sigma }}_{3,P}^{\ast }$ , respectively. If the shock stress in the composite is less than the through-thickness plain strain yield stress ( ${\mathit{\sigma }}_{3,C}^{\ast }$ ) of the composite, only an elastic wave front will be generated and the plastic zone will be absent.

Continuity of the particle velocity at the PCCI suggests that the absolute velocity of the elastic zone in the projectile is equal to the particle velocity in the plastic zone of the composite, ⁵⁰ which is also equal to the velocity of the PCCI, $\dot{H}$

$V_I-u_P^{\ast }=u_C^{\ast \ast }=\dot{Z}$

(1)

Considering linear elastic material behavior for the projectile, the stress in the elastic zone of the projectile can be written as

${\mathit{\sigma }}_{3,P}^{\ast }={\mathit{\rho }}_P{c}_{0,P}^{\ast }{u}_P^{\ast }=z_P^{\ast }{u}_P^{\ast }$

(2)

where $z_P^{\ast }={\mathit{\rho }}_P{c}_{0,P}^{\ast }$ , is the elastic impedance of the projectile. The impedance of the composite in the elastic zone can be defined as $z_C^{\ast }={\mathit{\rho }}_C{c}_{3L}^{\ast }$ .

Conservation of mass, momentum, and front locations for the elastic and plastic fronts in the composite (Figure 2(a)), and the location of the PCCI can be written as ⁴⁹

Elastic front conditions

$\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{s}:{\mathit{\rho }}_C^{\ast }\left(\left(u_C^{\ast }-0\right)-c_{3L}^{\ast }\right)=-{\mathit{\rho }}_C{c}_{3L}^{\ast }$

(3a)

$\mathrm{M}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{u}\mathrm{m}:{\mathit{\sigma }}_{3,C}^{\ast }-0={\mathit{\rho }}_C\left(u_C^{\ast }-0\right)c_{3L}^{\ast }$

(3b)

$\text{Location}:Z_e={\displaystyle \underset{0}{\overset{t}{\int }}{c}_{3L}^{*}}dt$

(3c)

Plastic front conditions

$\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{s}:{\mathit{\rho }}_C^{\ast \ast }\left(\left(u_C^{\ast \ast }-u_C^{\ast }\right)-c_p^{\ast \ast }\right)=-{\mathit{\rho }}_C^{\ast }{c}_{3p}^{\ast \ast }$

(4a)

$\mathrm{M}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{u}\mathrm{m}:{\mathit{\sigma }}_{3,C}^{\ast \ast }-{\mathit{\sigma }}_{3,C}^{\ast }={\mathit{\rho }}_C^{\ast }\left(u_C^{\ast \ast }-u_C^{\ast }\right)c_{3p}^{\ast \ast }$

(4b)

$\text{Location}:Z_p={\displaystyle \underset{0}{\overset{t}{\int }}\left(c_{3p}^{**}+u_C^{*}\right)}dt$

(4c)

PCCI

$\text{Location}:Z={\displaystyle \underset{0}{\overset{t}{\int }}\left(V_I-u_P^{*}\right)}dt={\displaystyle \underset{0}{\overset{t}{\int }}{u}_C^{**}}dt$

(5)

The through-thickness strain–density relationships under plain strain loading for elastic and plastic conditions for the composite material can be expressed as

${\mathit{\rho }}_C^{\ast }=\frac{{\mathit{\rho }}_C}{1-{\mathit{\epsilon }}_e^{\ast }}=\frac{{\mathit{\rho }}_C}{1-\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{E_{3e}^{/}}}=\frac{{\mathit{\rho }}_C}{1-\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{{\mathit{\rho }}_C{\left(c_{3L}^{\ast }\right)}^2}}=\frac{{\mathit{\rho }}_C}{1-\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{z_C^{\ast }{c}_{3L}^{\ast }}}$

(6a)

${\mathit{\rho }}_C^{\ast \ast }=\frac{{\mathit{\rho }}_C}{1-\mathit{\epsilon }}=\frac{{\mathit{\rho }}_C}{1-\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{E_{3e}^{/}}-\frac{{\mathit{\sigma }}_{3,C}^{\ast \ast }-{\mathit{\sigma }}_{3,C}^{\ast }}{E_{3p}^{/}}}$

(6b)

Composite density in the elastic zone given by equation (6a) can also be obtained from the solution of equations (3a) and (3b). However, solutions of equations (4a) and (4b) do not provide equation (6b) for the density of the composite in the plastic zone, rather provides the following solution

${\mathit{\rho }}_C^{\ast \ast }=\frac{{\mathit{\rho }}_C}{\left(1-\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{{\mathit{\rho }}_C{\left(c_{3L}^{\ast }\right)}^2}\right)\left(1-\frac{{\mathit{\sigma }}_{3,C}^{\ast \ast }-{\mathit{\sigma }}_{3,C}^{\ast }}{{\mathit{\rho }}_C^{\ast }{\left(c_{3p}^{\ast \ast }\right)}^2}\right)}$

(6c)

Thus, equating equations (6b) and (6c), one can determine the plastic front velocity with respect to immovable materials as

$c_{3p}^{\ast \ast }=\left(1-\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{E_{3e}^{/}}\right)\sqrt{\frac{E_{3p}^{/}}{{\mathit{\rho }}_C}}$

(7)

Equation (7) is obtained using the mass, momentum, and strain–density relationships, which was not known a priori. The stress in the elastic zone of the composite laminate is given by equation (3b)

${\mathit{\sigma }}_{3,C}^{\ast }={\mathit{\rho }}_C{c}_{3L}^{\ast }{u}_C^{\ast }=z_C^{\ast }{u}_C^{\ast }$

(8)

Using equation (8), the particle velocity in the elastic zone of the composite laminate can be expressed as

$u_C^{\ast }=\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{z_C^{\ast }}$

(9)

Stress in the plastic zone of the composite laminate is given by equation (4b), where the particle velocity, $u_C^{\ast \ast }$ , is unknown

${\mathit{\sigma }}_{3,C}^{\ast \ast }={\mathit{\sigma }}_{3,C}^{\ast }+{\mathit{\rho }}_C^{\ast }{c}_{3p}^{\ast \ast }\left(u_C^{\ast \ast }-u_C^{\ast }\right)$

(10)

In order to determine $u_C^{\ast \ast }$ , we use the continuity of stress at the PCCI

${\mathit{\sigma }}_{3,P}^{\ast }={\mathit{\sigma }}_{3,C}^{\ast \ast }$

(11)

Denoting the impedance of shocked composite by $z_C^{\ast \ast }={\mathit{\rho }}_C^{\ast }{c}_{3p}^{\ast \ast }$ , and using equations (10), (11), and (2), one can write

$z_P^{\ast }{u}_P^{\ast }={\mathit{\sigma }}_{3,C}^{\ast }+z_C^{\ast \ast }\left(u_C^{\ast \ast }-u_C^{\ast }\right),$

(12)

which leads to

$u_C^{\ast \ast }=\frac{z_P^{\ast }}{z_C^{\ast \ast }}{u}_P^{\ast }-\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{z_C^{\ast \ast }}+\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{z_C^{\ast }}$

(13)

Substituting equations (1) and (8) into equation (13) yields

$V_I-u_P^{\ast }=\frac{z_P^{\ast }}{z_C^{\ast \ast }}{u}_P^{\ast }-\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{z_C^{\ast }}\left(\frac{z_C^{\ast }}{z_C^{\ast \ast }}-1\right)$

(14)

which can be simplified to

$u_P^{\ast }=\frac{V_I+\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{z_C^{\ast }}\left(\frac{z_C^{\ast }}{z_C^{\ast \ast }}-1\right)}{\left(1+\frac{z_P^{\ast }}{z_C^{\ast \ast }}\right)}$

(15)

Using equations (1) and (15), one can write

$u_C^{\ast \ast }=\frac{V_I\frac{z_P^{\ast }}{z_C^{\ast \ast }}-\frac{{\mathit{\sigma }}_{3,C}^{\ast }}{z_C^{\ast }}\left(\frac{z_C^{\ast }}{z_C^{\ast \ast }}-1\right)}{\left(1+\frac{z_P^{\ast }}{z_C^{\ast \ast }}\right)}$

(16)

Equation (16) is necessary to calculate the stress of the composite in the plastic zone using equations (10) and (9)

${\mathit{\sigma }}_{3,C}^{\ast \ast }=\frac{z_P^{\ast }{V}_I}{\left(1+\frac{z_P^{\ast }}{z_C^{\ast \ast }}\right)}+\frac{{\mathit{\sigma }}_{3,C}^{\ast }\frac{z_P^{\ast }}{z_C^{\ast }}\left(\frac{z_C^{\ast }}{z_C^{\ast \ast }}-1\right)}{\left(1+\frac{z_P^{\ast }}{z_C^{\ast \ast }}\right)}$

(17)

The phase I dynamic PR force on the projectile can be expressed as

$F_{P,I}^D=A_P{\mathit{\sigma }}_{Z,C}^{\ast \ast }=A_P\left(\frac{z_P^{\ast }{V}_I}{\left(1+\frac{z_P^{\ast }}{z_C^{\ast \ast }}\right)}+\frac{{\mathit{\sigma }}_{3,C}^{\ast }\frac{z_P^{\ast }}{z_C^{\ast }}\left(\frac{z_C^{\ast }}{z_C^{\ast \ast }}-1\right)}{\left(1+\frac{z_P^{\ast }}{z_C^{\ast \ast }}\right)}\right)$

(18)

Equation (18) is the new 1D shock solution for the dynamic PR force acting on the projectile as a function of impact velocity, elastic–plastic composites properties, and the through-thickness crush strength of the composites.

If the shock stress at the PCCI is less than the plain strain yield stress ( ${\mathit{\sigma }}_{3,C}^{\ast }$ ) of the composite laminate, the plastic wave front will be absent. In this case, the elastic stress in the composite and the PR force on the projectile can be determined by setting $z_C^{\ast \ast }=z_C^{\ast }$ in equations (17) and (18)

${\mathit{\sigma }}_{3,C}^{\mathrm{\#}}=\frac{z_P^{\ast }{V}_I}{1+\frac{z_P^{\ast }}{z_C^{\ast }}}$

(19)

$F_{P,I}^{\mathrm{\#}}=A_P{\mathit{\sigma }}_{3,C}^{\mathrm{\#}}=A_P\frac{z_P^{\ast }{V}_I}{1+\frac{z_P^{\ast }}{z_C^{\ast }}}$

(20)

Equations (19) and (20) are the 1D elastic stress wave solutions for 1D impact and are available in literature. ⁵⁰ Derivation of equations (19) and (20) from fundamental relations are presented in Appendix A for reference.

Equations (20) and (18), respectively, provide the elastic and elastic–plastic solution for the shock stress as a function of the impact velocity of the projectile. Figure 3(A) presents the results for a cylindrical steel projectile impacting a PW S-2 glass/SC15 thick section composite plate (properties are listed in the section on Results and discussion). Two values of plain strain yield stress, that is ${\mathit{\sigma }}_{3,C}^{\ast }=$ 850 MPa and 1000 MPa, are used for calculations. Both elastic (equation (19)) and elatic–plastic (equation (17)) solutions are linear in impact velocity, intersect at about V_I = 185 m/s for ${\mathit{\sigma }}_{3,C}^{\ast }=$ 850 MPa, and at about V_I = 220 m/s for ${\mathit{\sigma }}_{3,C}^{\ast }=$ 1000 MPa (Figure 3(A)) and are the transition point above which the elastic–plastic solution is valid. Figure 3(A) can be used to calculate the initial PR force on a known cylindrical projectile at any impact velocity by multiplying the shock stress with the cross-sectional area of the projectile. This 1D solution is based on the elastic properties of the projectile and the elastic–plastic properties of the target composite materials. It is assumed that the projectile remains elastic at all times, which is true in the materials used in this article. However, if the shock stress exceeds the yield strength of the projectile, then our solution will not be valid and one should consider elastic–plastic wave propagation in the projectile. In addition, in the particle velocity in the plastic zone, $u_C^{\ast \ast }$ , exceeds the plastic wave front velocity, ${\dot{Z}}_p=c_{3p}^{\ast \ast }+u_C^{\ast }$ , the present elastic–plastic model will not remain valid. This condition is shown in Figure 3(B), where $u_C^{\ast \ast }=c_{3p}^{\ast \ast }+u_P^{\ast }$ , at V_I = 1,240 m/s.

Figure 3.

Elastic and elastic–plastic solution of shock stress as a function of the impact velocity of the projectile. Plain weave (PW) S-2 Glass/SC15 composite plate, cylindrical steel projectile of diameter, $D_P=$ 12.7 mm. (A) Shock stress vs. impact velocity. (B) Particle velocities vs. impact velocity.

Stepwise solution procedure

The PR force as a function of impact velocity of the projectile given by equation (18) will remain constant until the compressive elastic stress wave front in the projectile reflects back from the free surface of the noncontact end of the projectile as tensile stress wave and reaches the PCCI. The time required for the elastic stress wave front to reach the free end of the projectile is given by

${\mathit{\tau }}_P=\frac{H_P}{c_{0,P}^{\ast }}$

(21)

where ${\mathit{\tau }}_P$ is the characteristic time of the projectile. Once the elastic stress wave front in the projectile reaches the free rear surface of the projectile in time $t={\mathit{\tau }}_P$ , the absolute velocity of the projectile becomes $V_I-u_P^{\ast }$ . When the elastic stress wave front reflects back from the rear surface of the projectile as tensile wave and reaches the PCCI in time $t=2{\mathit{\tau }}_P$ , the absolute velocity of the projectile becomes $V_I-2u_P^{\ast }$ . The particle velocities, stress, and PR force given by equations (15) to (18) remain constant in the time range $0<t<2{\mathit{\tau }}_P$ , and the quantities for this solution step 1 can be denoted as $u_{P(1)}^{\ast }$ , $u_{C(1)}^{\ast }$ , $u_{C(1)}^{\ast \ast }$ , ${\mathit{\sigma }}_{3,C(1)}^{\ast \ast }$ , and $F_{P,I(1)}^D$ . For the subsequent solution steps (i.e. in the time ranges, $2{\mathit{\tau }}_P<t<4{\mathit{\tau }}_P$ , $6{\mathit{\tau }}_P<t<8{\mathit{\tau }}_P$ , . . . , $2\left(\mathrm{k}\mathrm{-}\mathrm{1}\right){\mathit{\tau }}_P<t<2\mathrm{k}{\mathit{\tau }}_P$ , etc.), the following recursive equations can be used

$V_{I(\mathrm{k}\mathrm{+}\mathrm{1})}=V_{I(\mathrm{k}\mathrm{)}}-2u_{P(\mathrm{k})}^{\ast }$

(22)

$\mathrm{\Delta }{u}_{P(\mathrm{k})}^{\ast }=-\frac{2u_{P(\mathrm{k})}^{\ast }}{\left(1+\frac{z_P^{\ast }}{z_C^{\ast \ast }}\right)}$

(23)

$u_{P(\mathrm{k}\mathrm{+}\mathrm{1})}^{\ast }=u_{P(\mathrm{k})}^{\ast }+\mathrm{\Delta }{u}_{P(\mathrm{k})}^{\ast }$

(24)

$u_{C(\mathrm{k}+1)}^{\ast \ast }=u_{C(\mathrm{k})}^{\ast \ast }+\mathrm{\Delta }{u}_{P(\mathrm{k})}^{\ast }\frac{z_P^{\ast }}{z_C^{\ast \ast }}$

(25)

$F_{P,I(\mathrm{k}\mathrm{+}\mathrm{1})}^D=F_{P,I(\mathrm{k})}^D+\mathrm{\Delta }{u}_{P(\mathrm{k})}^{\ast }{A}_P{z}_P^{\ast }$

(26)

where $\mathrm{k}=\mathrm{1}\mathrm{,}\mathrm{2}\mathrm{,}\cdot \cdot \cdot \mathrm{,}\mathrm{n}$ . This procedure can be continued until the elastic stress wave front in the composite reflects back from the rear surface of the composite and reaches the plastic wave front, and the phase I 1D stepwise shock solution is valid up to this phase I characteristics travel time of the composite and can be expressed as

${\mathit{\tau }}_{C,I}=\frac{2H_C}{c_{3p}^{\ast \ast }+c_{3L}^{\ast }+u_C^{\ast }}=\mathrm{T}$

(27)

Discontinuity between solution steps

At the end of solution step 1 and at time $t=2{\mathit{\tau }}_P^{-}$ , the particle velocity in the projectile at PCCI is $u_{P,PCCI}^{\ast -}=V_I-2u_{P(1)}^{\ast }$ . On the other hand, the particle velocity at the PCCI in the shocked plastic zone at the same time is $u_{C,PCCI}^{\ast \ast -}=V_I-u_{P(1)}^{\ast }$ . Since $u_{P,PCCI}^{\ast -}<u_{C,PCCI}^{\ast \ast -}$ , the interface separates, and a release wave propagates in the composite from the PCCI toward the plastic front with a propagation velocity of $c_{3L}^{\ast }$ . One can assume that the particle velocity in the shocked plastic zone undergoes a relaxation process and becomes $u_{C,PCCI}^{\ast \ast +}=u_{C(1)+}^{\ast \ast }=V_I-2u_{P(1)}^{\ast }$ in time $\mathrm{\Delta }t$ , and at the same time, the contact between the projectile and the composite is reestablished. In the meanwhile, the location of the PCCI moves from $Z_{(1)}$ to $Z_{(1)+}$ , the plastic front moves from $Z_{p(1)}$ to $Z_{p(1)+}$ , and the elastic front moves from $Z_{e(1)}$ to $Z_{e(1)+}$ . Conservation of mass and momentum during this time increment, $\mathrm{\Delta }t$ , can be written as

$\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{s}:{\mathit{\rho }}_{C(1)}^{\ast \ast }\left(Z_{p(1)}-Z_{(1)}\right)+{\mathit{\rho }}_{C(1)}^{\ast }\left(Z_{e(1)}-Z_{p(1)}\right)={\mathit{\rho }}_C{Z}_{e(1)}$

(28a)

${\mathit{\rho }}_{C(1)+}^{\ast \ast }\left(Z_{p(1)+}-Z_{(1)+}\right)+{\mathit{\rho }}_{C(1)+}^{\ast }\left(Z_{e(1)+}-Z_{p(1)+}\right)={\mathit{\rho }}_C{Z}_{e(1)+}$

(28b)

Momentum

$\begin{array}{rl}& {\mathit{\rho }}_{C(1)}^{\ast \ast }\left(Z_{p(1)}-Z_{(1)}\right)u_{C(1)}^{\ast \ast }+{\mathit{\rho }}_{C(1)}^{\ast }\left(Z_{e(1)}-Z_{p(1)}\right)u_{C(1)}^{\ast }\\ & ={\mathit{\rho }}_{C(1)+}^{\ast \ast }\left(Z_{p(1)+}-Z_{(1)+}\right)u_{C(1)+}^{\ast \ast }+{\mathit{\rho }}_{C(1)+}^{\ast }\left(Z_{e(1)+}-Z_{p(1)+}\right)u_{C(1)+}^{\ast }\end{array}$

(29)

Setting

$\begin{array}{r}{\mathit{\rho }}_{C(1)}^{\ast }={\mathit{\rho }}_{C(1)+}^{\ast }\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{u}_{C(1)}^{\ast }=u_{C(1)+}^{\ast }\\ u_{C(1)}^{\ast \ast }=V_I-u_{P(1)}^{\ast }\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{u}_{C(1)+}^{\ast \ast }=V_I-2u_{P(1)}^{\ast }\\ Z_{(1)}=\left(V_I-u_{P(1)}^{\ast }\right)2{\mathit{\tau }}_P\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{Z}_{(1)+}=\left(V_I-u_{P(1)}^{\ast }\right)2{\mathit{\tau }}_P+\left(V_I-2u_{P(1)}^{\ast }\right)\mathrm{\Delta }t\\ Z_{p(1)}=\left(c_p^{\ast \ast }+u_{C(1)}^{\ast }\right)2{\mathit{\tau }}_P\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{Z}_{p(1)+}=\left(c_p^{\ast \ast }+u_{C(1)+}^{\ast }\right)\left(2{\mathit{\tau }}_P+\mathrm{\Delta }t\right)\\ Z_{e(1)}=c_{3L}^{\ast }\left(2{\mathit{\tau }}_P\right)\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{Z}_{e(1)+}=c_{3L}^{\ast }\left(2{\mathit{\tau }}_P+\mathrm{\Delta }t\right);\end{array}$

and solving the mass and momentum equations 31 and 32, we can determine the value of $\mathrm{\Delta }t$

$\mathrm{\Delta }t=\frac{2{\mathit{\tau }}_P{u}_{P(1)}^{\ast }}{\left(V_I-2u_{P(1)}^{\ast }\right)-\frac{u_{C(1)}^{\ast }}{c_p^{\ast \ast }}\left(c_p^{\ast \ast }+u_{C(1)}^{\ast }-c_{3L}^{\ast }\right)}$

(30)

after which the second solution phase can begin with continuity of stress and displacement at the PCCI. In our stepwise solution for phase I, we will present this time gap between consecutive solution steps.

Linear approximation of phase I elastic–plastic PR force

If the stepwise solution procedure is continued until time $t=T$ , the number of full solution steps, n, can be written as

$n=T\mathrm{\setminus }{\mathit{\tau }}_P$

(31)

with n additional time discontinuity and one additional partial time step. We can use these step solutions to make a linear approximation of the PR force as a function of projectile displacement

$F_{P,I}^{D,Linear}=F_{P,I}^{D,1}-f_{H}H$

(32)

where the second superscript 1 denotes solution step 1. In order to determine the coefficient, $f_H$ , we will equate the total work done by the stepwise solution to its linear approximation. The work done by the stepwise solution and the linear approximation can be expressed as

$\begin{array}{r}{W}_{T,I}^{Step}=\sum _{i=1}^n{F}_{P,I}^{D,i}\cdot u_{P,i}^{\ast \ast }\cdot {\mathit{\tau }}_P+\frac{1}{2}\sum _{i=1}^n\left(F_{P,I}^{D,i}+F_{P,I}^{D,i+1}\right)\cdot \left(u_{P,i}^{\ast \ast }-u_{P,i}^P\right)\cdot \mathrm{\Delta }t\\ +F_{P,I}^{D,n+1}\cdot u_{P,n+1}^{\ast \ast }\cdot \left(T-n\left({\mathit{\tau }}_P+\mathrm{\Delta }t\right)\right)\end{array}$

(33)

$W_{T,I}^{Linear}={\displaystyle \underset{H=0}{\overset{H_{T,I}}{\int }}\left(F_{P,I}^{D,1}-f_{H}H\right)}dH=F_{P,I}^{D,1}{H}_{T,I}-f_H\frac{{\left(H_{T,I}\right)}^2}{2}$

(34)

where the total projectile displacement of the projectile can be calculated by

$H_{T,I}=\sum _{i=1}^n{u}_{P,i}^{\ast \ast }\cdot {\mathit{\tau }}_P+\sum _{i=1}^n\left(u_{P,i}^{\ast \ast }-u_{P,i}^P\right)\cdot \mathrm{\Delta }t+u_{P,n+1}^{\ast \ast }\cdot \left(T-n\left({\mathit{\tau }}_P+\mathrm{\Delta }t\right)\right)$

(35)

The coefficient $f_H$ , can then be determined by equating equations and

$f_H=\frac{2}{H_{T,I}^2}\left(F_{P,I}^{D,1}{H}_{T,I}-W_{T,I}^{Step}\right)$

(36)

Rate-dependent plain strain yield stress of the composites

In the plain strain elastic–plastic material model for composites presented in Figure 2, the rate dependency of the plain strain yield stress can be modeled using a logarithmic function ⁴⁷

${\mathit{\sigma }}_{Z,{\dot{\mathit{\epsilon }}}_P}^{\ast }={\mathit{\sigma }}_{Z,0}^{\ast }\left(1+C_{FC}ln\left({\dot{\mathit{\epsilon }}}_P/{\dot{\mathit{\epsilon }}}_{P,0}\right)\right)$

(37)

where ${\dot{\mathit{\epsilon }}}_{P,0}$ is a reference plastic strain rate which can be taken as 1/s, and the rate coefficient $C_{\mathrm{F}\mathrm{C}}$ should be experimentally determined. The plastic strain rate under plain strain impact loading can be calculated from the relative locations of the PCCI and the plastic wave front. The location of the PCCI and the location of the plastic wave front at time, ${\mathit{\tau }}_P$ , can be expressed as

$H=u_P^{\ast \ast }{\mathit{\tau }}_P$

(38)

$Z_P=\left(c_P^{\ast \ast }+u_P^{\ast }\right){\mathit{\tau }}_P$

(39)

The average plastic strain over time, ${\mathit{\tau }}_P$ , can be approximated by

${\mathit{\epsilon }}_P=\frac{H}{Z_P}$

(40)

and the plastic strain rate by

${\dot{\mathit{\epsilon }}}_P=\frac{{\mathit{\epsilon }}_P}{{\mathit{\tau }}_P}=\frac{u_P^{\ast \ast }}{\left(c_{3L}^{\ast \ast }+u_P^{\ast }\right){\mathit{\tau }}_P}$

(41)

Since the particle velocities in the elastic and plastic zone of the composite are not known a priori, one can use the quasi-static value of the plain strain yield stress to estimate the particle velocities and then determine the rate-dependent plain strain yield stress using iterative procedures.

Equations (22) to (26) represent the 1D shock compression model of projectile impact on composites during short-time phase I deformation, and equation (32) is a linear approximation of the 1D shock solution. The dynamic effect or the stress wave propagation effect is incorporated in the 1D shock compression model, and the strain rate effect on the through-thickness material properties can be included in the 1D shock compression model using equations (37) and (41). The linear approximation of the phase I 1D shock compression model is combined with the phase II semi empirical quasi-dynamic PR (QD-PR) model to describe the entire penetration event.

Semiempirical (QD-PR) model for phase II ballistic penetration

We assume that at the end of phase I shock compression, the material starts to soften and shows tension shear–dominated PR. Since the behavior of phase II penetration is not known a priori, we propose a semiempirical model to estimate that from QS-PR force–displacement data. The QS-PR curve presented in Figure 1 can be rewritten as a function of the location of the PCCI, H, including a DIF as

$F_{P,II}^{QD}=-\mathrm{D}\mathrm{I}\mathrm{F}\cdot F_{P,II}\left(1-\frac{H_F-H_0}{H-H_0}\right)$

(42)

where $H_0$ is a reference displacement and $H_F$ is the final displacement when the projectile completely stops. The value of $\mathrm{D}\mathrm{I}\mathrm{F}=1.0$ for quasi-static penetration and $\mathrm{D}\mathrm{I}\mathrm{F}>1.0$ for dynamic penetration. The value of DIF is not known a priori, and we will later present a methodology to estimate that from ballistic test results. Equation (42) is the proposed QD-PR model of ballistic penetration, which considers the dynamic effects.

The projectile deceleration as a function of projectile displacement can now be expressed as

${\ddot{H}}_{II}=\mathrm{D}\mathrm{I}\mathrm{F}\cdot \frac{F_{P,II}}{M_P}\left(1-\frac{H_F-H_0}{H-H_0}\right)$

(43)

The linear approximation of phase I 1D shock compression solution (equation (32)) intersects with the phase II QD-PR curve (equation (42)) at displacement, $H=H_{I,II}^D$ . The point of intersection can be determined by equating equations (32) and (42) at $H=H_{I,II}^D$ (Appendix A).

The phase I projectile displacement, velocity, acceleration, and PR force can be determined from 1D shock compression model as a function of time up to a time $t=t_{I,II}^D$ , when $H=H_{I,II}^D$ , and can be used as the initial condition for phase II penetration process. Equation (43) can be used to calculate the acceleration of the projectile ( $\ddot{H}$ ) as a function of projectile displacement (H) with known initial conditions, $H_{I,II}^D,{\dot{H}}_{I,II}^D,{\ddot{H}}_{I,II}^D$ . With an arbitrary small increment of projectile displacement, $\mathrm{\Delta }{H}_i$ , the corresponding increment in time and velocity can be expressed as

$\mathrm{\Delta }{t}_i=\frac{\mathrm{\Delta }{H}_i}{{\dot{H}}_i}$

(44)

$\mathrm{\Delta }{\dot{H}}_i=\frac{{\ddot{H}}_i}{{\dot{H}}_i}\mathrm{\Delta }{H}_i$

(45)

Which then can be used to calculate the time and velocity for the next increment

$t_{i+1}=t_i+\mathrm{\Delta }{t}_i$

(46)

${\dot{H}}_{i+1}={\dot{H}}_i+\mathrm{\Delta }{\dot{H}}_i$

(47)

Using the quasi-dynamic PR force as a function of displacement given by equation (42), the acceleration of the projectile as a function of projectile displacement can be determined by equation (43). The time for the projectile as a function of projectile velocity can be determined using equations (44) and (46), and the projectile velocity as a function of projectile deceleration and velocity can be calculated using equations (45) and (47), which completes the phase II solution for projectile displacement, velocity, acceleration, and PR force as a function of time.

The combined theoretical-semi-empirical penetration model of ballistic penetration

The phase I 1D shock compression model and the phase II QD-PR curve can now be combined to provide the complete solution. The work done by the projectile can be found by integrating the PR force described by equation (32) for phase I, and by equation (42) for phase II with respect to the projectile displacement

$\begin{array}{c}{W}_P={\displaystyle \underset{0}{\overset{H_{I,II}^D}{\int }}{F}_{P,I}^{D,Linear}dH}+{\displaystyle \underset{H_{I,II}^D}{\overset{H_F}{\int }}{F}_{P,II}^{QD}dH}\\ =F_{P,I}^{D,1}{H}_{I,II}^D\left(1-\frac{f_H{H}_{I,II}^D}{2F_{P,I}^{D,1}}\right)-DIF\cdot F_{P,II}\left(\left(H_F-H_{I,II}^D\right)-\left(H_F-H_0\right)\ln \frac{H_F-H_0}{H_{I,II}^D-H_0}\right)\end{array}$

(48)

Equating the work done by the projectile to the initial impact kinetic energy (W_P = M_PV_I* ² /2), the projectile velocity (V_I*) required for complete penetration can be estimated

$V_{I*}=\sqrt{\frac{2F_{P,I}^{D,1}{H}_{I,II}^D}{M_P}\left(1-\frac{f_H{H}_{I,II}^D}{2F_{P,I}^{D,1}}\right)-\frac{2\cdot DIF\cdot F_{P,II}}{M_P}\left(\left(H_F-H_{I,II}^D\right)-\left(H_F-H_0\right)\ln \frac{H_F-H_0}{H_{I,II}^D-H_0}\right)}$

(49)

The combination of 1D shock compression model with QD-PM can now be termed as the new theoretical–semiemperical model of ballistic penetration. In this model, the DIF is unknown and should be estimated following the energy balance procedures. If the experimental ballistic limit velocity is known, the DIF can be estimated by equating the ballistic limit energy to the total energy calculated from the combined penetration model

${\mathrm{E}}_{\mathrm{B}\mathrm{L}}^{\mathrm{B}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}}=\frac{1}{2}{M}_P{V}_{\mathrm{B}\mathrm{L}}^2$

(50a)

$\begin{array}{r}{\mathrm{E}}_{\mathrm{B}\mathrm{L}}^{\mathrm{S}\mathrm{h}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{-}\mathrm{Q}\mathrm{D}\mathrm{-}\mathrm{P}\mathrm{M}}=F_{P,I}^{d,1}{H}_{I,II}^D\left(1-\frac{f_H{H}_{I,II}^D}{2F_{P,I}^{D,1}}\right)\\ -\mathrm{D}\mathrm{I}\mathrm{F}\cdot F_{P,II}\left(\left(H_F-H_{I,II}^D\right)-\left(H_F-H_0\right)ln\frac{H_F-H_0}{H_{I,II}^D-H_0}\right)\end{array}$

(50b)

or

$\mathrm{D}\mathrm{I}\mathrm{F}=\frac{{\mathrm{E}}_{\mathrm{B}\mathrm{L}}^{\mathrm{B}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}}-F_{P,I}^{D,1}{H}_{I,II}^D\left(1-\frac{f_H{H}_{I,II}^D}{2F_{P,I}^{D,1}}\right)}{-F_{P,II}\left(\left(H_F-H_{I,II}^D\right)-\left(H_F-H_0\right)ln\frac{H_F-H_0}{H_{I,II}^D-H_0}\right)}$

(51)

Since $H_{I,II}^D$ is a function of $\mathrm{D}\mathrm{I}\mathrm{F}$ , $H_{I,II}^D$ can be estimated at $\mathrm{D}\mathrm{I}\mathrm{F}=1.0$ to get the first estimate of_DIF, and then using the trial-and-error method of solution the value of both $H_{I,II}^D$ and DIF can be determined. Once a value of DIF is assumed, the theoretical–semiempirical penetration model equations can be solved to predict the time history of projectile velocity, displacement, and PR force by numerical time integration method. Results are presented next.

Results and discussion

The new theoretical–semiempirical penetration model equations are solved using numerical integration procedures in a MATLAB program. Geometric and mechanical properties of the right circular cylinder projectile and the PW S-2 glass/SC15 composites used in the ballistic experiments of Ref. 44 are used to compute the 1D shock compression model parameters and are presented in Table 1. The plastic modulus of the composite under through-thickness plain strain impact loading condition is estimated from QS through-thickness indentation testing of composite plates of thickness H_C = 7.62 mm with right circular cylinder indenters of different diameters (D_P = 6.3 mm, 7.9 mm. 9.5 mm, and 11.1 mm). Figure 4 shows the compressive stress as a function of plastic or nonlinear strain, from where the average plane strain yield stress for this model is calculated to be ${\mathit{\sigma }}_{Z,C}^{\ast }=$ 850 MPa, and the plastic modulus to be $E_{3,P}^{/}=$ 2.5 GPa. The initial softening behavior in the experimental data is due to the finite size indenter which does not represent true plain strain behavior. In order to represent true plain strain behavior, the linear fit to the experimental data, which is represented by the solid line, is shifted to pass through the plain strain yield point and is represented by the dotted line.

Figure 4.

Compressive stress plastic strain behavior of Plain weave (PW) S-2 glass/SC15 composites.

Table 1.

Dimensions and properties of the projectile and the composite plate

Projectile
${\mathit{\rho }}_P^0$ (gm/cm3)	mP (g)	DP (mm)	HP (mm)	EP (GPa)	$c_0^P$ (m/s)	${\mathit{\tau }}_P$ (μs)
7.85	13.8	12.70	13.90	210.00	5172.2	5.37
Composite plate
${\mathit{\rho }}_C^0$ (gm/cm3)	HC (mm)	E1 (GPa)	E2 (GPa)	E3 (GPa)	${\mathit{\nu }}_{12}={\mathit{\nu }}_{21}$	${\mathit{\nu }}_{31}={\mathit{\nu }}_{32}$
1.87	13.20	27.50	27.50	11.80	0.11	0.18
${\mathit{\nu }}_{13}={\mathit{\nu }}_{23}$	Ɗ	C33 (GPa)	$c_{3L}^{\ast }$ (m/s)	$E_{3,P}^{/}$ (GPa)	${\mathit{\sigma }}_{Z,C}^{\ast }$ (MPa)	$c_P^{\ast \ast }$ (m/s)
0.4195	9.192 × 10-32	14.21	2771.6	2.50	850.0	1092.9
$u_P^{\ast }$ (m/s)	T (μs)	$\mathrm{\Delta }t$ (μs)
165.77	6.55	0.3001

Phase I 1D shock-compression model predictions

The particle velocity corresponding to the through-thickness yield stress ( ${\mathit{\sigma }}_{Z,C}^{\ast }$ ) of the composite is calculated to be $u_P^{\ast }=$ 166 m/s, which tells us that impact velocities less than this velocity will only generate elastic stress wave propagation in the composite. Figure 5(a) shows the phase I dynamic PR force on the projectile (equation (18)) as a function of projectile displacement for different impact velocities (200 m/s < V_I < 800 m/s). Equation (18) provides a stepwise solution and the time step duration for the present problem (Table 1) is ${\mathit{\tau }}_P=$ 5.37 μs, the time duration of discontinuity between solution steps is $\mathrm{\Delta }t=$ 0.30 μs, and the total solution time for phase I 1D shock compression model is T = 6.55 μs. Equations (16) and (18) are used to calculate the particle velocity and PR force at the PCCI, which shows that both the velocity and force at PCCI is constant for each time steps ( ${\mathit{\tau }}_P$ ). The constant velocity and time steps are used to calculate the displacement of the PCCI and are presented in Figure 5(A). Phase I theoretical solutions are compared with the QS-PR curve shown in Figure 1. The PR force on the projectile at ballistic limit velocity ( $V_{\mathrm{B}\mathrm{L}}=$ 367 m/s) and at time t = 0 is found to be $F_{P,I}^{D,1}=$ 160 kN. The area under the phase I 1D shock solution at ballistic limit velocity is found to be $E_{\mathrm{B}\mathrm{L}}^{\mathrm{1}\mathrm{D}\mathrm{-}\mathrm{S}\mathrm{h}\mathrm{o}\mathrm{c}\mathrm{k}}=$ 400 J, which is 43% of ballistic experiments ( $E_{\mathrm{B}\mathrm{L}}^{\mathrm{B}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}}=$ 930 J) in this example. This is an important finding that in T = 6.55 μs of time, 43% of the impact energy is transmitted to the composite laminate.

Figure 5.

Phase I dynamic penetration resistance force and plastic strain rate. H_C = 13.2 mm, τ_P = 5.37 us, T = 6.55 us. (A) Penetration resistance force. (B) Plastic strain rate.

By connecting the ends of the theoretical contact forces for different impact velocities, a 1D shock envelope can be plotted to differentiate between phase I and phase II penetration behavior. A linear approximation of the stepwise force–displacement solution is determined using equations (32) to (36) and is presented in Figure 5(A) as solid lines without symbols. The plastic strain rate that develops during the first solution step (equation (41)) is linear with impact velocity and is presented in Figure 5(B). If the through-thickness yield stress and plastic modulus of composites under dynamic plain strain loading can be determined as a function of plastic strain rate from experiments, the present 1D shock model will allow one to consider the effect of strain rate as well.

Combined shock-QS-PR model predictions

The QS-PR curve developed from QS-PST experiments (Figure 1) can be combined with the phase I shock compression model to predict the projectile dynamics and is presented in Figure 6. The shock-QS-PR model can predict the projectile velocity, displacement, and time at and above the ballistic limit velocity of the composite plate. The shock-QS-PR model also provides an analytical expression to calculate the approximate V_BL (equation (49)), value of which is found to be $V_{\mathrm{B}\mathrm{L}}^{\mathrm{S}\mathrm{h}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{-}\mathrm{Q}\mathrm{S}\mathrm{-}\mathrm{P}\mathrm{R}}=$ 356 m/s, which is about 3% less than the ballistic experiment, $V_{\mathrm{B}\mathrm{L}}^{\mathrm{B}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}}=$ 367 m/s. The energy at ballistic limit velocity is calculated to be $E_{\mathrm{B}\mathrm{L}}^{\mathrm{S}\mathrm{h}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{-}\mathrm{Q}\mathrm{S}\mathrm{-}\mathrm{P}\mathrm{R}}=$ 877 J, which is 5.7% less than the ballistic experimental value of $E_{\mathrm{B}\mathrm{L}}^{\mathrm{B}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c}}=$ 930 J. Figure 6(A) shows that the combined shock QS-PR model does not include the complete phase I 1D shock solution/envelope (Figure 5(A)) and thus the prediction of projectile dynamics at higher impact velocities (Figure 6(B)) may not be accurate. The shock QS-PR model considers the phase I dynamic and stress wave propagation effects up to a certain limit but neglects the dynamic effects in phase II penetration and is the subject matter of the shock QD-PM presented in the subsequent section.

Figure 6.

Shock  quasi-static penetration resistance (QS-PR) model input and prediction. (A) Penetration resistance force as a function of projectile displacement. (B) Projectile velocity as a function of time.

The combined theoretical–semiempirical (shock Quasi-Dynamic) penetration model predictions

Dynamic effects in the phase II penetration are considered by introducing the DIF in equation (42), the value of which can be estimated using equation (51). The value of DIF which predicts accurate V_BL based on the energy balance equation (51) is estimated to be 1.10 and its effect on projectile dynamics and on PR is investigated next. The shock QD-PM calculations are compared with the LS-Dyna FEA solutions (validated with experimental results) presented in Ref. 46. Figure 7(A) shows the shock QD-PM input PR force as a function of projectile displacement for DIF = 1.00 and 1.10 at an impact velocity V_I = 367 m/s. The LS-Dyna FEA prediction of PR force ⁴⁶ at V_I = 360 m/s showed rebound of the projectile, while at V_I = 380 m/s showed complete penetration; however, the force–displacement at both the impact velocities appeared to have a small difference up to projectile displacement H = 15 mm. Similarly, the difference in the shock QD-PM force–displacement input for DIF = 1.00 and DIF = 1.10 is also small, however, appears to be the average of the LS-Dyna predictions. A small difference in force–displacement may appear as a large difference in velocity time response as shown in Figure 7(B). The LS-Dyna prediction clearly shows the rebound of the projectile at V_I = 360 m/s with V_R = −26.7 m/s, and complete penetration at V_I = 380 m/s with V_R = 38.3 m/s. Ideally, at ballistic limit velocity V_BL = 367 m/s, the residual velocity should be zero, V_R = 0.00 m/s, however, the shock QD-PM predicted V_R = 0.08 m/s at time t = 871 μs (not shown in Figure 7(B)) for DIF = 1.10, and predicted V_R = 84.6 m/s at time t = 185 μs for DIF = 1.00.

Figure 7.

Validation of shock quasi-dynamic (QD) penetration model with LS-Dyna analyses. (A) Force vs. displacement. (B) Velocity as a function of time.

Clearly, the new theoretical–semi empirical penetration model (shock QD-PM) can predict the ballistic limit velocity, V_BL, and the time history of the projectile velocity accurately with the correct choice of DIF, which is a constant. However, prediction of projectile dynamics above V_BL will require theoretical modeling of phase II penetration and will be addressed elsewhere.

Summary and conclusions

Impact and penetration dynamics of a cylindrical projectile on a thick section composite laminate is presented. A theoretical 1D shock compression model has been developed for short-time phase I shock compression to determine the PR force on the projectile and the projectile dynamics. A rate-sensitive elastic–plastic material behavior is used to model the through-thickness plastic shock and elastic stress wave propagation upon impact. Conservation of mass and momentum at the wave fronts are used to calculate the stress and particle velocities in the shock/plastic and elastic zones to determine the PR force. For the 13.2-mm-thick S-2 Glass/SC15 composite laminate under impact of a 13.8-g cylindrical projectile (D_P = 12.7 mm) at ballistic limit velocity $V_{\mathrm{B}\mathrm{L}}=$ 367 m/s, the 1D shock compression solution predicted that the PR force at time t = 0 is 160 kN, and 43% of the projectile impact energy is transmitted to the composite laminate in this short-time phase I.

The phase I theoretical 1D shock compression solution is combined with the phase II QS-PR curve obtained from the QS-PST experimental methodologies, and a ballistic limit velocity of $V_{\mathrm{B}\mathrm{L}}^{\mathrm{S}\mathrm{h}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{-}\mathrm{Q}\mathrm{S}\mathrm{-}\mathrm{P}\mathrm{R}}=$ 356 m/s (3% less than ballistic experiment) and a ballistic limit energy of $E_{\mathrm{B}\mathrm{L}}^{\mathrm{S}\mathrm{h}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{-}\mathrm{Q}\mathrm{S}\mathrm{-}\mathrm{P}\mathrm{R}}=$ 877 J (5.7% less than ballistic experiment) are predicted. The QS-PR curve is multiplied with a DIF to model the phase II dynamic effects. Energy-based equations are derived to estimate the DIF from QS-PR calculations, and when used, the new combined theoretical–semiempirical penetration model (shock QD-PM) accurately predicted the ballistic limit velocity and energy. The shock-QD-PM prediction of projectile velocity at ballistic limit velocity correlated well with the FEA predictions. However, the shock QD-PM cannot predict the projectile dynamics at higher impact velocities than the ballistic limit and will require further development.

Highlights

• A theoretical formulation based on 1D-shock theory of a RCC projectile on a finite thickness composite plate is presented.

• The plastic wave-front velocity in the composite is determined using mass, momentum, and strain-density relationships (equation (7)). A new finding.

• A new 1D-shock solution for the dynamic penetration resistance force (equation (18)) acting on the projectile is obtained which is a function of impact velocity of the projectile, the elastic-plastic impedances, and the through-thickness crush strength of the composites.

• The new 1D-Shock stepwise solution is discontinuous between solution steps. Duration of this discontinuity is quantified solving the mass and momentum equations (equation (30)) and is necessary for multiple solution steps in very thick composites.

• A dynamic increase factor, DIF, is added to a quasi-static penetration resistance (QS-PR) curve to formulate the Phase II Quasi-Dynamic Penetration Resistance (QD-PR) model. A numerical integration scheme is used to calculate the incremental projectile velocity as a function of time and displacement.

• The theoretical 1D Shock-Compression model is combined with the semi-empirical QD-PR model as the new Theoretical-Semi-Empirical or Shock-Quasi-Dynamic Penetration model to predict the ballistic limit of thick-section composites and the projectile dynamics at the ballistic limit. Good correlations with previous FEA results are obtained.