The transverse ballistic impact on a two-dimensional (2D) membrane causes a truncated deformation cone to develop in the wake of tensile implosion waves. Here, the cone wave reflected from the finite boundaries of the elastic membrane has been studied analytically. A first-order linear nonhomogeneous differential equation for the ratio of the reflected cone wave front velocity to the speed of tensile waves is derived, which is further used to calculate the traveling time taken by the reflected cone wave to reach to the projectile surface. Since the reflected wave starts when the membrane is already in a deformed configuration, the speed of the reflected cone wave is a function of radius *r* in the cylindrical coordinates as opposed to almost constant speed of the incoming cone wave studied in the literature. The analytical results are validated with molecular dynamics (MD) simulations of the ballistic impact of projectiles onto a single layer of coarse-grained (CG) graphene. In the second part of the paper, we analyze the membrane impact problem for linear isotropic viscoelastic materials and find that the tensile wave speed for stresses and displacements is the same as that obtained in the case of a linear isotropic elastic material. We also show that only under special conditions, self-similar solutions for the cone wave are possible in viscoelastic materials modeled by Maxwell, Kelvin–Voigt, or a combination of similar models. Our findings lay some grounds on which further studies on the ballistic response of viscoelastic materials can be performed.

## Introduction

The research in effective body armors with low cost, lightweight, and increased ballistic protective performance has a long history [1]. Protective systems against ballistic impact have several other applications too, mainly in the composite parts of aerodynamic vehicles, boat hulls, electronic parts, and coated fabrics. Several analytical, numerical, and experimental studies have been performed to assess ballistic impacts on high performance yarns, fabrics, metals, and other materials [2]. They involve deriving constitutive relations, modeling overall fabric behavior, increasing perforative performance, studying the effect of layers, size, temperature, friction, heat generation, etc. [3].

The analytical literature on transverse impact of two-dimensional (2D) membranes has come a long way from the early works of Grigoryan [4] and Rakhmatulin and Dem'͡$$ianov [5] to the recent works of Walker [6], Phoenix and Porwal [2], Naik et al. [7], and several others. All of these works take into account different stress waves originating in the system after the impact traveling at different speeds. The transverse impact on 2D membrane also causes truncated cone to develop in the wake of tensile implosion waves. Normally, the membrane radius is infinite and reflections of the waves are not accounted for. But Freeston et al. [8] developed a numerical model to describe the strain build-up in yarns due to multiple strain wave reflections from yarn crossover interactions in a woven fabric under the ballistic impact. Ha-Minh et al. [9] have included the reflection of the longitudinal tensile waves from the adjacent yarn in their analytical model for the ballistic impact behavior of 2D plain-weave fabrics. However, the reflection of the cone wave due to finite size of the target has not been considered in the ballistic impact literature. The study of these reflections is increasingly becoming possible after some recent advances in experiments and computer simulations of graphene membranes subject to the ballistic impact. The laser-induced projectile impact tests on finite length (85 *μ*m) multilayer graphene (MLG) have shown that the MLG membranes can be excellent ballistic barrier materials [10]. Inspired by this work, atomistic molecular dynamics (MD) simulations over graphene membranes of nanometer sizes have been carried out which are able to capture stress wave propagation speed, formation of the cone wave, specific penetration energy, etc. [11,12]. In another work [13], with the help of MD simulations, we showed the existence of a critical membrane size below which the cone wave reflections from the boundaries induce perforation resulting into reduction of ballistic proof capability. The present work focuses on a detailed analytical study of the reflection cone wave.

In the first part of the paper, Sec. 2, we study the reflected cone wave occurring due to the finiteness of the specimen size in the analytical framework of the membrane model of Phoenix and Porwal [2]. The model was developed to study the stress and cone waves caused by the ballistic impact on fibrous systems. We derive a first-order linear nonhomogeneous differential equation for the ratio of the reflected cone wave front velocity to the speed of tensile waves. This is used to calculate the time taken by the reflected cone wave from the boundaries of the membrane to reach to the projectile surface. Then, we validate our analytical results with MD simulations of the ballistic impact of projectiles onto a single layer of coarse-grained (CG) graphene. The CG model of graphene maps every four atoms of graphene onto a coarse-grained bead using a mapping that preserves hexagonal symmetry [14]. Our model captures the elasticity, strength, and fracture properties of the monolayer graphene and MLG very well as reported in the previous studies [14] and has the potential to capture nano- and mesoscale failure mechanisms and size effects in the MLG systems. Details regarding the CG model have been provided in the Supplemental material which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection. The projectiles are made of CG beads forming a diamond cubic lattice structure with a lattice constant of 0.72 nm, which is double the actual lattice constant of diamond, in consistent with the CG mapping scheme of graphene. The default projectile bead mass is 96 g/mole, which gives rise to a density of 3.42 g/cm^{3}, close to the actual diamond density of 3.5 g/cm^{3}. More detailed information about the simulation settings can be found in Supplemental material which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection .

In the second part of the paper, Sec. 3, we analyze the 2D membrane impact problem for the linear isotropic viscoelastic materials. The deformation of a thin viscoelastic membrane of infinite size under the ballistic impact with constant projectile velocity is considered. We study the problem with two most general viscoelastic models, Maxwell and Kelvin–Voigt, and find that the tensile wave speed for stresses and displacements is the same as that for the linear isotropic elastic material. We also show that self-similar solutions for the cone wave generated in the wake of the tensile wave are possible only under special conditions. We conclude in Sec. 4 with a summary and suggestions for future work.

## Reflected Cone Wave Analysis

We first consider the deformation of a thin linear isotropic elastic membrane with fixed boundaries due to the ballistic impact by a right circular cylinder (projectile) of radius *r _{p}*. The membrane of radius

*a*is stress-free before the projectile hits at the center of the membrane. Both tensile and cone wave develop in the wake of the impact. Following the incoming cone wave hitting the fixed boundaries of the membrane, the cone wave reflects from the boundary, travels toward the center of the membrane, and starts to deform it further along its traveling path. This has been shown in Fig. 1 with the help of an MD simulation where the membrane radius is

*a*= 1000 Å and the projectile radius is

*r*= 115 Å. The circular membrane dimensions are along the

_{p}*X*- and

*Z*-axes of the MD box and the projectile is traveling along the

*Y*-axis. The front views in Fig. 1 are shown in such a way that the

*Z*-axis goes into the plane and axisymmetric deformation of the membrane becomes visible. Figure 1(a) shows the membrane at the time of impact,

*t*= −

*t*, where

_{p}*t*is the time taken by the cone wave to reach from the projectile surface

_{p}*r*=

*r*to the outermost edge

_{p}*r*=

*a*of the membrane. We take our time origin

*t*= 0 exactly at the point when the incoming wave reaches the membrane boundary and the reflected cone wave starts propagating backward toward the projectile surface as shown in Fig. 1(b). Figure 1(c) shows the propagation of the reflected cone wave at time

*t*=

*t*

_{1}> 0. Here, we locate a membrane element at a distance

*r*from the center under the deformation due to the reflected cone wave. We produce the illustration of the cylindrical coordinate system (

*r*,

*ϕ*,

*y*) and displacements (

*u*, 0,

*v*) of the material element with center coordinates (

*r*,

*ϕ*, 0) in Fig. 1(d) inspired by Fig. 2(b) of Phoenix and Porwal [2], where

*u*and

*v*are displacements along radial

*r*and downward

*y*directions, respectively. Assuming axisymmetric deformation of the membrane, the black curve in Fig. 1(d) shows the curvature of the cone shape at time

*t*= 0. We analyze this reflected cone wave with two basic assumptions: (a) For time

*t*> 0, the projectile has already achieved near zero velocity, i.e.,

*V*(

_{p}*t*= 0) ≈ 0. This makes sure that the deformation caused by the projectile after time

*t*> 0 is negligible. (b) The additional strain caused by the fixed boundary is negligible. Following the MD simulations, we observe that both of these assumptions are justified. The impact velocity of the projectile at time

*t*= −

*t*is

_{p}*V*(

_{p}*t*= −

*t*) = 500 m/s, where

_{p}*t*is approximately equal to 59 ps. The projectile velocity at time

_{p}*t*= 0 is

*V*(

_{p}*t*= 0) ≈ 50 m/s. For our case, the incoming cone wave speed and the tensile wave speed are

*c*≈ 1.5 km/s and

*a*

_{0}≈ 20 km/s, respectively. This means $Vp(t=0)<Vp(t=\u2212tp)<c\u226aa0$. So, the projectile keeps going down with a low velocity after the incoming cone wave reaches the membrane edge but this low velocity does not cause much deformation at the projectile surface. From our simulations, we note that the local in-plane strain defined in Eq. (3) at the projectile surface

*r*=

*r*and time

_{p}*t*= 0 is

*ε*≈ 0.035, whereas at time $t=trp\u2212$, just before the reflected cone wave hits the projectile surface

_{t}*r*=

*r*, we obtain

_{p}*ε*≈ 0.045. We neglect this additional strain during the time period

_{t}*t*= 0 to

*t*=

*t*

_{rp}in which we limit our analysis. This assumption helps simplify our analysis and we leave it upon future study to relax this assumption. Figure 1(e) shows the membrane shape when the reflected cone wave reaches the projectile surface and Fig. 1(f) shows the immediate membrane failure.

*u*,

*v*) caused by the reflected cone wave. Assuming that

*u*always remains negative, the displacements are given by

and *γ* is the angle between the local tangential surface to the membrane element and the ground plane.

*ε*at time

_{t}*t*= 0

^{–}as a function of radial distance

*r*from the center. This can be obtained from Eq. (80) of Phoenix and Porwal [2] as

*r*=

*a*and time

*t*= 0

^{–}, i.e., $\epsilon c0=\epsilon t(r=a;t=0\u2212)$. Further, it can be seen from Eq. (92) of Phoenix and Porwal [2] that the ratio of the speed of the incoming cone wave

*c*to that of the tensile wave

*a*

_{0}is related with $\epsilon c0$ as

Now we intend to verify Eq. (4) with the help of the MD simulation data of the projectile impact on the graphene membrane. A typical case has been chosen with the membrane radius *a* = 1000 Å and the projectile radius *r _{p}* = 115 Å. The plot of

*ε*computed with the help of the MD data at time

_{t}*t*= 0

^{−}for

*r*=

*r*to

_{p}*r*=

*a*is shown in Fig. 2. We try to fit the analytical expression of

*ε*(

_{t}*r*;

*t*= 0

^{–}) given by Eq. (4) over the MD data and obtain an excellent fit. The nonlinear least square solver “lsqcurvefit” of matlab is used for this purpose with the only fitting coefficient $\epsilon c0$. With the help of $\epsilon c0$ obtained from this fitting procedure and Eq. (5), we find that

*c*/

*a*

_{0}= 0.069±0.001 with 95% confidence interval. We also calculate

*c*and

*a*

_{0}independently from the MD data assuming that

*c*remains constant for the incoming cone wave, and as noted above in the beginning of this section, we find that

*c*≈ 1.5 km/s,

*a*

_{0}≈ 20 km/s, and

*c*/

*a*

_{0}≈ 0.075. This is quite close to the analytical fit result, and therefore, validates the analytical relation expressed in Eq. (4).

We seek self-similar solutions for (*u*, *v*) with fixed boundaries at *r* = *a*, i.e., $u(r=a,t)=v(r=a,t)=0$, based on the reflected cone wave front velocity, *c _{r}*, which is a function of time

*t*as the reflected wave propagates. We define two functions

*α*(

_{r}*t*): =

*c*(

_{r}*t*)/

*a*

_{0}, which is the ratio of the reflected cone wave speed to the tensile wave speed, and $z(t):=a\u2212cr(t)t=a\u2212\alpha r(t)a0t$, where $\u222b0t\alpha r(t)a0dt$ is the distance traveled by the reflected cone wave from the boundary

*r*=

*a*in time

*t*. The function

*z*(

*t*) is analogous to the variable

*z*introduced in the analysis of the incoming cone wave (see Eq. (32)) and helps obtain self-similar solutions for the reflected cone wave.

*α*(

_{r}*t*) in the analytical framework of Sec. 3.2 in Phoenix and Porwal [2]. Self-similar solutions mean that functions $u\u2032$ and $v\u2032$ can be defined as $u\u2032(r,z)\u2261u(r,t)$ and $v\u2032(r,z)\u2261v(r,t)$, respectively. Without loss of generality, we drop the primes over the displacement functions $(u\u2032,v\u2032)$ and treat them as (

*u*,

*v*). Equations (1) and (2) in terms of

*z*can be written as

*α*(

_{r}*t*). We note that the forms of Eqs. (8) and (9) are exactly similar to Eqs. (57) and (58) in Phoenix and Porwal [2], except for the fact that

*α*becomes $\alpha r+t\alpha r\u02d9$. Taking $\alpha :=\alpha r+t\alpha r\u02d9$, we recover exactly same equations from Eq. (57) to Eq. (72) of Phoenix and Porwal [2] with only change in

*α*getting replaced by $\alpha r+t\alpha r\u02d9$. Defining the variables

*ξ*=

*r*/

*z*and

*ξ*=

_{a}*ξ*(

*r*=

*a*;

*z*) so that

*ξ*becomes only a function of

*r*for a fixed

*z*, we note that

*ξ*varies from

*ξ*to the outer cone edge

_{a}*ξ*= 1. In this sense,

*α*becomes a parameter for Eq. (76) of Phoenix and Porwal [2] and we further recover Eq. (76) to Eq. (91) of the same paper. Defining

*ε*

_{rc}≡

*ε*(

_{t}*ξ*= 1;

*z*) as the in-plane strain at the outer cone edge of the reflected wave, we rewrite Eq. (92) of Phoenix and Porwal [2] as

*ε*

_{rc}≪ 1 is a function of time

*t*, we can rewrite Eq. (10) as $\alpha r+t\alpha r\u02d9=g(t)$, where $g=\epsilon rc$, and then we solve Eq. (10) as first-order linear nonhomogeneous equation. The complementary homogeneous solution to this ordinary differential equation (ODE) is

*α*=

_{r}*C*/

*t*, where

*C*is a constant. The particular solution to this ODE is $\alpha r=(1/t)\u222b0tg(t)dt$. Therefore, the actual solution is $\alpha r=C/t+(1/t)\u222b0tg(t)dt$. To avoid any singularity at

*t*= 0, we can have

*C*= 0, and therefore, $\alpha r=(1/t)\u222b0tg(t)dt$. This is time-average of

*g*(

*t*), so this should be the same as the average over the distance $|r\u2212a|$ traveled by the reflected cone wave in time

*t*. This means $\alpha r\u0303(r):=\alpha r(t)$ can be expressed as

*ε*

_{rc}can be treated as a function of

*r*. We can safely assume that the in-plane strain at the outer cone edge

*ε*

_{rc}is the sum of two strains: (i) the existing strain

*ε*(

_{t}*r*;

*t*= 0

^{−}) caused by the incoming wave at

*t*= 0

^{−}, and (ii) the strain caused by the reflected wave, which should be approximately the same as the strain $\epsilon c0$ at the outer cone edge of the incoming cone wave. So using Eq. (4),

*ε*

_{rc}(

*r*) can be written as

For our MD simulation case described earlier, we have *r _{p}*/

*a*= 0.115. From Eqs. (15) and (16), we obtain

*t*

_{rp}= 0.575

*a*/

*c*and

*t*= 0.885

_{p}*a*/

*c*. This means

*t*

_{rp}= 0.65

*t*, analytically. We also calculate

_{p}*t*

_{rp}and

*t*from the MD data and find that

_{p}*t*

_{rp}≈ 0.69

*t*. Thus, MD results are very close to the analytical results and the small differences can be attributed to the assumptions made in the analysis. We also plot time ratio $(trp/tp)$ versus

_{p}*ζ*=

*r*/

_{p}*a*in Fig. 3 and find that the time ratio $(trp/tp)$ is a monotonically increasing function of

*ζ*with minimum around 0.63 and maximum around 0.71. This completes our analysis of the reflected wave for the present purpose.

## Ballistic Impact Response for Viscoelastic Materials

### Governing Equations Under Constant Projectile Velocity.

*t*> 0 due to the ballistic impact by a cylindrical projectile of radius

*r*at

_{p}*t*= 0. Working with cylindrical coordinate system (

*r*,

*ϕ*,

*y*), the membrane element with center coordinates (

*r*,

*ϕ*,

*y*= 0) at

*t*= 0 undergoes the displacement (

*u*, 0,

*v*) at time

*t*, where

*u*and

*v*are the in-plane and out-of-plane components of the membrane displacement vector, respectively. It can be shown that these displacements are governed by the following partial differential equations (PDEs) [2,5]:

*ρ*is the membrane initial density,

*σ*and

_{t}*σ*are the local engineering normal stress components along the deformation curve and along the circumferential direction (hoop stress), respectively, and

_{ϕ}*E*and

*ν*are Young's modulus and Poisson's ratio, respectively, and the strains

*ε*and

_{t}*ε*for any linear isotropic material are given by Eq. (3) and

_{ϕ}*f*(

*t*) in time space

*t*to the complex frequency space

*s*such that $L{f(t)}:=f\u0302(s)=\u222b0\u221ef(t)e\u2212stdt$, yield the following results:

*ϕ*and used the zero initial conditions at

*t*= 0, i.e., $u(r,t=0)=(\u2202u/\u2202t)(r,t=0)=0$ and $v(r,t=0)=(\u2202u/\u2202t)(r,t=0)=0$. In order to come up with stress–strain relations analogous to Eqs. (22) and (23) in the complex frequency domain, we need to replace

*E*and

*ν*with their viscoelastic analogies $sEv\u0302$ and $\nu v\u0302$, which are functions of the corresponding transformed operators of time-dependent deviatoric

*P*

_{1},

*Q*

_{1}and dilatational

*P*

_{2},

*Q*

_{2}linear differential operators for the stress–strain constitutive equations [15]

*ν*for elastic membranes, we also assume that the effect of operator $\nu v\u0302$ is negligible in comparison to other operators. This means $Q2\u0302\u2248P2\u0302Q1\u0302/P1\u0302$. We also assume the hydrostatic pressure independence of the elastic response, which yields $P2\u0302=1$, implying $Q2\u0302\u2248Q1\u0302/P1\u0302$. Using Eq. (27), we obtain the transformed modulus $sEv\u0302(s)$ as

*E*and

*η*refer to Young's modulus and viscosity, respectively. Based on Eq. (22) and all the assumptions made earlier, Table 1 also provides the stress–strain relations in the transformed space. Since $\sigma t\u2009cos\u2009\gamma \u302a=sEv\u0302(s)\epsilon t\u2009cos\u2009\gamma \u302a$, Eqs. (25) and (26) become

From these equations, it becomes clear that the stress and discontinuity waves propagate with the speed $a0=(E/\rho )$ for the Maxwell and the K–V model, respectively.^{2} Unlike the pure elastic waves, these stress and discontinuity waves for the Maxwell and the K–V model attenuate and possess a relaxation time *τ* = *η*/*E* so that for time *t* ≪ *τ* the behavior differs little from that of the elastic membrane [17], but for time *t* ≫ *τ*, the stress and discontinuity waves relax and get dissipated.

### Cone Wave Analysis.

*c*with the assumption that

*c*is constant. We let

*α*=

*c*/

*a*

_{0}be the ratio of cone wave speed to tensile wave speed and hope to find a self-similar solution in terms of the variable

*z*

where *t _{p}* =

*r*/

_{p}*c*is the time taken by the cone wave to travel the distance

*r*. We show below that self-similar solutions for the cone wave can only be obtained under special conditions for both Maxwell and K–V models, and therefore, we need to solve full PDEs for complete solutions. For this purpose, we first start with the Maxwell model.

_{p}#### Maxwell Model.

*τ*=

*η*/

*E*is the relaxation time. In terms of

*z*, these PDEs become

*ξ*=

*r*/

*z*, we define two new displacement functions $X(\xi ):=(u+r)/z=(u/r)\xi +\xi $ and $Y(\xi ):=v/z=(v/r)\xi $ in terms of which we have

*ξ*. Substituting these into Eqs. (35) and (36), we obtain

*ξ*. Since

*r*is an independent variable, hence self-similarity transformation is possible only when either (a)

*τ*→

*∞*, i.e., material is viscous, or (b) $X\u2212\xi X\u2032=Y\u2212\xi Y\u2032=0$, which means

*C*

_{1}> 0 and

*C*

_{2}> 0 are constants. The conditions in Eq. (42) are purely geometrical as they can be satisfied for any Maxwell material under proper conditions. Since $X(\xi =r/z)=(u+r)/z$ and $Y(\xi =r/z)=v/z$, the conditions in Eq. (42) also imply $u=(C1\u22121)r$ and

*v*=

*C*

_{2}

*r*, which means that we obtain simple shear motion in the

*r*–

*z*plane. Under these conditions, the PDEs of elastodynamics turn into similar equations as obtained for the linear isotropic elastic membrane for which it has been shown that [2]

#### Kelvin–Voigt Model.

*a*

_{0}and

*τ*have the same meaning as above. In terms of

*z*, these PDEs become

*ε*are given by Eqs. (39)–(41). Here too, we obtain self-similar solutions only for the special cases for which Eqs. (48) and (49) turn into ODEs in

_{t}*ξ*. Since

*z*is an independent variable, this is possible only under the conditions when either (a)

*τ*→ 0, i.e., material is elastic, or (b) $(\xi 2\epsilon t\u2032\u2009cos\u2009\gamma )\u2032=(\xi 2\epsilon t\u2032\u2009sin\u2009\gamma )\u2032=0$, which means

where *C*_{3} and *C*_{4} are constants. We again treat the conditions in Eq. (50) as purely geometrical since any Kelvin–Voigt material can satisfy them under proper conditions. As noted for the Maxwell model, the analysis follows the elastic membrane hereafter for which the cone wave speed is given by the relation in Eq. (43). Using Eq. (41) and the relations in Eq. (50), the constants *C*_{3} and *C*_{4} must satisfy $C32+C42=\epsilon t\u20322|\xi =1$.

When these special conditions are not available, we fail to obtain self-similar solutions and then the governing PDEs for *u* and *v* in Eqs. (35) and (36) for the Maxwell model and Eqs. (46) and (47) for the K–V model need to be solved for the purpose of obtaining the cone wave speed. For any general constitutive model [16] of similar kinds, by extension, it can be said that the self-similar solution for cone wave is possible only for special cases.

## Conclusion

In this paper, we first study the cone wave reflected from the finite boundaries of the elastic membrane in the framework of Phoenix and Porwal [2]. We are able to derive a first-order linear nonhomogeneous differential equation for the ratio of the reflected cone wave front velocity to the speed of tensile waves, which is used to calculate the time taken by the reflected cone wave from the boundaries of the membrane to reach to the projectile surface. The analytical results are verified with the MD simulations of the ballistic impact of projectiles onto a single layer of coarse-grained graphene. Then, we analyze the 2D membrane impact problem for the linear isotropic viscoelastic materials and find that the tensile wave speed for stresses and displacements is the same as that obtained for the linear isotropic elastic material. Then, we also show that self-similar solutions for the cone wave are possible only under special conditions in the absence of which full governing PDEs must be solved. For the present purpose, only the cone wave analysis has been done for the viscoelastic materials. We also plan to analyze the deceleration of projectile due to reactive forces of the viscoelastic membrane. The analysis can also be extended to the impact response of a viscoelastic balanced fabric. Another interesting thing would be to study the widening and thickening of the crazes in polymers under the ballistic impact.

## Acknowledgment

The authors acknowledge the support from the Departments of Civil and Environmental Engineering and Mechanical Engineering at Northwestern University, as well as the Northwestern University High Performance Computing Center for a supercomputing grant.

## Funding Data

Army Research Office (Award No. W911NF-13-1-0241).

A detailed discussion on the wave speeds in the frequency domain can be found in Eringen [16].