Modeling and Simulation of Ballistic Penetration of Ceramic-Polymer-Metal Layered Systems

Numerical simulations and analysis of ballistic impact and penetration by tungsten alloy rods into composite targets consisting of layers of aluminum nitride ceramic tile(s), polymer laminae, and aluminum backing are conducted over a range of impact velocities on the order of 1.0 to 1.2 km/s. Computational results for ballistic efficiency are compared with experimental data from the literature. Simulations and experiments both demonstrate a trend of decreasing ballistic efficiency with increasing impact velocity. Predicted absolute residual penetration depths often exceed corresponding experimental values. The closest agreement between model and experiment is obtained when polymer interfaces are not explicitly represented in the numerical calculations, suggesting that the current model representation of such interfaces may be overly compliant. The present results emphasize the importance of proper resolution of geometry and constitutive properties of thin layers and interfaces between structural constituents for accurate numerical evaluation of performance of modern composite protection systems.


Introduction
Modern protection systems often consist of layers of ceramic, metallic, and/or polymer-based components.Interfaces between layers may strongly influence performance of such systems under ballistic impact.However, the importance of interfacial characteristics-for example, interface thickness, material type, and bonding strength-is not fully understood in many cases.Furthermore, the accuracy of available computational tools to assess such effects has not heretofore been thoroughly quantified.The purpose of this study is the assessment of one computational tool-with typical/default user options enabled-for modeling ballistic impact and penetration of a layered target consisting of one or more ceramic tiles backed by a thick metallic plate, with thin layers of polymer between the tiles in some cases.The current focus is the evaluation of the fidelity of the existing material models (including corresponding property parameters) and related numerical methods; modification of constitutive models or calibration of user-defined parameters to best match experimental ballistic results is beyond the scope of the present study.
As discussed in detail later, the penetrator-target configuration studied in this work duplicates that featured in ballistic experiments of Yadav and Ravichandran [1].Prior to description of the specific problem investigated here and in [1], an overview of literature on the subject is warranted.Other experiments of ballistic impact and penetration of ceramic targets with various interlayers and/or backing materials include those described in [2][3][4][5].Analytical models used to describe and/or predict ballistic penetration and possible perforation of such systems include those presented in [6][7][8].Numerical simulations invoking finite element methods, for example, those of ballistic impact of ceramic systems, are described in [9][10][11].Principles of dimensional analysis applied to relate properties and performance of armor ceramics are developed in [12,13].Comprehensive descriptions of terminal ballistics with applications to brittle solids can be found in several additional lengthy references [14][15][16].Figure 1: Ballistic problem (a) projectile and target (three tiles) [1]; (b) finite element mesh.

Problem Statement
In particular, the penetrator-target configuration simulated in this work depicts that examined in experiments of [1].
As shown in Figure 1(a), a WHA (Tungsten Heavy Alloy) penetrator, cylindrical in shape with flat nose, impacts a target at velocity  ranging from ≈1000 m/s to ≈1200 m/s at null obliquity.The respective length  and diameter  of the penetrator are 50.6 mm and 8.43 mm (/ = 6).
The target consists of one, three, or six tiles of aluminum nitride (AlN), an isotropic polycrystalline armor ceramic.The total thickness of the tile(s) is 38.1 mm in all cases.A thin polyurethane laminate separates neighboring tiles in the experiments when the target contains multiple tiles.Ballistic performance of the ceramic-polymer system (or a single tile in some cases) is quantified by residual penetration depth into a 6061-T6 aluminum (Al) backing block of thickness 76.2 mm, which was sufficient to fully stop the penetrator in all reported experiments [1].
The main result ascertained from the experimental study was that ballistic efficiency was the highest (best) for the three tiles each of thickness 12.7 mm, intermediate for a single tile of thickness 38.1 mm, and the lowest (worst) for six tiles each of thickness 6.35 mm [1].Lateral tile dimensions were 101.6 mm × 101.6 mm.It was speculated that soft polymer layers in the three tile configuration enabled dispersion of the initial, primary compressive shock wave that caused more severe damage in the single tile configuration.On the other hand, bending and tensile failure modes were posited to strongly and negatively influence penetration resistance of the six tile configurations compare to offsetting any benefits obtained by dispersion or attenuation of the initial compressive shock attributed to the presence of compliant polymer layers and weak interfaces.
The computational tool implemented in this study is the EPIC (Elastic Plastic Impact Calculation) finite element code [20] (2013 release).This numerical analysis tool was chosen for two primary reasons: (i) its existing library of material constitutive models and property database are extensive and were thought to be sufficient for the representation of behaviors of each component (i.e., the ceramic, polymer, and metals as listed in Table 1) and (ii) its graphical user interface permits rapid generation of finite element meshes for ballistic penetration simulations of layered targets, as shown, for example, in Figure 1(b).

Mathematical Theory and Numerical Methods
Given initial and boundary conditions, the finite element method for dynamic analysis in a Lagrangian framework seeks solutions of the governing equations of continuum mechanics-conservation of mass, momentum, and energy-written, respectively, in local form as [21]   + ∇ ⋅  = 0, Mass density is ; the particle velocity vector is  = u/; internal energy per unit mass is ; the spatial gradient operator is ∇; partial time derivatives / are taken with respect to fixed material coordinates (i.e., material time derivatives).For a general finite deformation mechanics problem involving an elastic-inelastic solid, the deformation gradient F and volume ratio  obey, with ∇ 0 being the reference gradient operator, 1 being the second-order unit tensor, and u being the particle displacement vector The deformation gradient is decomposed multiplicatively into elastic (superscript ) and inelastic or plastic (superscript ) parts.Assuming small elastic deformation, or assuming an additive decomposition of the rate of stretching into elastic (D  ) and inelastic (D  ) parts in the spatial frame independent from such a multiplicative split, the spatial velocity gradient can be written, with  being skew spin tensor The Cauchy stress tensor  is symmetric and can be split into deviatoric (  ) and hydrostatic parts, letting  denote pressure; also the scalar (Mises) effective stress   is defined in the following equation: ) Denoting by C a state-dependent tangent elastic modulus tensor of order four, the objective Jaumann rate of Cauchy stress obeys a general constitutive equation of the form The state of stress, temperature  (or internal energy ), cumulative inelastic deformation   , and cumulative damage  generally depend on the history of the displacement gradient at each material point via a constitutive model that depends on material type.In practice, for isotropic solids considered herein, ( 5) is replaced by distinct constitutive equations for pressure  (an equation-of-state depending on  and ) and deviatoric stress components.When the material deforms elastically,   < , where  is the effective strength of the solid that generally depends on strain, strain rate, temperature, pressure, and damage.When plastic deformation occurs, the yield condition   =  is enforced numerically via a radial return algorithm.Adiabatic conditions are assumed in (1), a standard practice for dynamic impact problems, leading to the energy balance below, with  being the specific heat per unit mass: Damage variable  is updated via an incremental constitutive equation of the form where  is the instantaneous equivalent strain to fracture that generally depends on state variables.Numerical discretization of global forms of equations in ( 1) is described in [20], for example.Specific materials considered herein are described by particular constitutive equations and associated model parameters, for pressure , strength , and fracture strain .Such equations are listed in what follows next for material classes covered in Table 1.
For the metallic solids (aluminum backing and tungsten rod), the following constitutive equations are used to dictate pressure, strength, and failure behaviors: Here,  * is a dimensionless, normalized total strain rate,  * is homologous temperature, Γ is the Gruneisen coefficient, and , ,   ,   , and   are other material parameters calibrated from experimental data.
For the aluminum nitride ceramic, Here,   is the pressure increment due to bulking,  is a material parameter, as are other terms with  and  subscripts that may take different values for intact and comminuted material [17].The model also accounts for a possible phase change, and corresponding details can be found in [17].
For the polymer, specifically a crushable polyurethane foam, Here,   is the crush pressure beyond which the response is nonlinear with   being the corresponding (elastic) volumetric strain.In (12),   is a material parameter describing volumetric inelastic strain at failure, which also enters a modified form of (7) accounting for volume change as well as cumulative deviatoric plastic strain.Further details can be found in [2,18].As shown in Table 2, cases with and without polymer layers were simulated.In the former, the thickness of polymer layers was restricted by constraints imposed by the mesh generator to a minimum value of 1.054 mm, about four times thicker than the value of 0.254 mm tested experimentally [1].Resolution of the latter very small thickness would require extremely small finite elements, which in turn would drastically increase computational cost through time step reductions imposed by the Courant condition [20], written explicitly later in (14).
Material models were selected from code library options that best matched those of the experiments; details can be found in Table 1.A notable discrepancy is that the density of the polyurethane polymer material used in experiments is somewhat larger (a factor of 3.8) than the most dense polyurethane foam of the available constitutive models.Default options for element failure were imposed in all simulations: tetrahedral elements were eroded [22] when scalar effective strains exceeded a value of 1.5.Nodal masses were conserved upon element erosion, but strength and pressure were zeroed for failed/eroded elements.Frictionless contact between projectile and target was imposed by default along slide-lines.Interfaces were assigned one of two conditions: (i) tied bonding, corresponding to shared nodes and perfect coherence or (ii) free contact, corresponding to duplicate nodes along distinct, interacting frictionless surfaces.In some simulations involving multiple tiles, the polymer layers were excluded.The very thin coating of epoxy used to glue the rearmost tile to the backing block in experiments was not modeled explicitly.Far-field boundary conditions corresponded to free surfaces; that is, the targets were unconfined as in the experiments, though effects of interaction with the mounting apparatus were necessarily excluded in the simulations to maintain a reasonable problem size.
Prior to simulations of the ceramic-polymer-metallic targets, simulations of penetration of the bare backing metal were conducted, similar to those reported experimentally [1].The thickness of the bare metal target was not listed in the experimental study; a value of 6 was used in the simulations, ensuring independence of residual penetration depth  0 from target thickness.A simulation time of 1.0 ms was sufficient for cessation of relative motion of the residual eroded projectile mass to that of the target.Impact velocities of 1030, 1100, and 1160 m/s were considered.
Next, numerical simulations of the layered targets were conducted for the same three impact velocities, as listed in where  is the residual penetration depth into the aluminum backing behind the interface between the backing and rearmost aluminum nitride tile, and  0 is the residual penetration depth into the bare backing at the same impact velocity .
When the projectile completely penetrated the backing metal thickness of 76.2 mm [1], a value of zero was assigned to .In such cases, the residual velocity   of the penetrator at a time of 1.0 ms was recorded (see Table 2) and used as a metric for performance comparisons.Tetrahedral finite element meshes were generated using the EPIC preprocessor, most often with the default fine mesh setting and expanded grid, and the latter feature leading to progressive mesh coarsening with increasing the distance from the penetration zone.This mesh density was found to yield sufficiently mesh size-independent results for residual penetration depths and resulting ballistic efficiency ; in fact, an even coarser medium mesh setting was usually deemed sufficient but was not used.Refer to Table 3 for details comparing fine and medium mesh densities for particular simulations involving one ceramic tile with free bonding, impacted at 1160 m/s.The verification that the results are independent of the time step restriction imposed for numerical integration of the rate (e.g., linear momentum and stress update) equations is also shown in Table 3 for the same target configuration.In particular, Δ max is the maximum allowable time step in a corresponding simulation, which in all cases is smaller than that imposed by the Courant condition necessary for stability of solutions obtained by explicit numerical integration of the equations of motion [20]: with ℎ min and   being the minimum lineal element size in the meshed domain and the effective longitudinal sound speed, respectively.Numerical simulations were executed in parallel mode on 16 processors using the available 2013 version of the EPIC code on the Spirit cluster at the US Air Force Research Laboratory (AFRL).Wall-clock execution times were always less than 24 hours.

Results
Predictions are compared with experiments for the bare backing metal (Al) in Figure 2(a), wherein a linear fit to the data was sufficient to describe computational results for the three impact velocities considered in each case: with particular values of dimensionless constant  0 and constant  1 [s/m] embedded within Figure 2(a).The deformed finite element geometry corresponding to residual penetration at 1.0 ms is shown in Figure 2(b); notice that the damaged zone exceeds the penetration depth of the partially eroded projectile in this case.
Predicted penetration depths significantly exceed experimental values.Reasons for the differences in results cannot be isolated in the present set of complex simulations, but possibilities include the following: the WHA material may be weaker than that depicted by the model, or the Al material may be stronger than that depicted by the model; the erosion criterion invoked in simulations may be too liberal for the Al or too strict for the WHA; omission of friction and commensurate wear between target and eroding projectile may result in larger penetration depths in simulations than those observed in experiments; and/or far-field boundary conditions may artificially affect depth of penetration results at later computation times in finite element simulations.
Representative results from various target configurations and impact velocities are shown in Figure 3, all corresponding to a solution time of 1 ms.In particular, in Figure 3(a), the penetrator barely defeats the single ceramic tile and resides just inside the metal backing plate (/ = 0.136 in Table 2).In Figure 3(b), the entire target-including three ceramic tiles, two layers of polymer, and metal backing plate-has been perforated by the projectile, and all layers of polymer laminate have been highly eroded.The latter result agrees qualitatively with experimental observation of severe damage in polymer layers of recovered targets [1].In Figure 3(c), the initially unbonded six ceramic tiles have been shattered by the projectile which remains lodged at the back-free surface of the aluminum backing at  = 1 ms.
Ballistic efficiencies from simulations and experiments are compared in Figure 4.Note that overlapping data points in Figure 4, for example, those when  → 0 in many instances, can be discerned by examining corresponding numerical values listed in Table 2.In Figure 4(a), simulation results for  for a single ceramic tile exceed those from experiments when the ceramic is perfectly bonded (tied) to the backing plate, while agreement with experiment is closer for free contact between ceramic and backing.For results of the three tile configurations shown in Figure 4(b), experimental values of  exceed simulation predictions regardless of numerical bonding representation or inclusion of polymer layers, though the closest agreement is obtained when polymer layers are omitted in the simulations.In Figure 4(c), the same conclusion is drawn for the six tile configurations; that is, the closest agreement is obtained when the polymer layers are not explicitly represented in the calculations.The simulations do tend to reflect the experimentally observed trend of decreasing ballistic efficiency with the increasing impact velocity.When ranked via descending ballistic penetration resistance, experimental results [1] suggest an ordering of three, one, and then six tiles, while simulation results suggest an ordering of one, three, and then six tiles.
Residual velocities from simulations at 1.0 ms for threeand six-tile target configurations are shown in Figures 5(a) and 5(b), respectively.Recall that complete perforation did not occur in any reported experiment [1].Residual velocities are similar for free interfaces and for tied bonding with polymer, confirming failure and commensurate erosion of the polymer layers, and consistent with efficiency results shown in Figure 4.

Analysis and Discussion
As inferred from examination of solution data in Table 2 and Figures 4 and 5, results suggest that incorporation of compliant polymer layers promotes bending modes and tensile fracture in the ceramic layers, leading to decreased ballistic efficiency relative to simulations wherein polymer is omitted.For example, consider efficiency and residual velocity predictions for the three tile configurations impacted at 1030 m/s and listed in Table 2.When bonding is free, ballistic efficiency  decreases from 0.628 to 0.518 when polymer interlayers are inserted between the ceramic tiles.When bonding is tied, efficiency decreases to zero, and residual velocity becomes nonzero (specifically,   / of 0.115).Similar, but not identical, trends, are evident for the six tile configurations, whereby projectile defeat occurs only when polymer is omitted, with nonzero residual velocities reported whenever polymer layers are included.Furthermore, increasing the number of tiles, while decreasing the individual tiles' thickness, exacerbates this weakness of the target package, especially when more polymer layers are included with an increasing total number of tiles.Stress wave propagation for simulations with three tiles for impact velocities of 1030 m/s is evident in Figure 6, which specifically shows hydrostatic pressure contours (positive in compression) at a time of 0.25 ms after initial impact.The penetration cavity is wider and shallower without polymer (Figure 6(a)), with minor differences in pressure waves emanating from the cavity evident among all cases shown.
As noted already in the context of penetration results for the bare backing metal, the source of discrepancy between model and experiment could not be isolated in these complex multimaterial calculations, but several possibilities can be suggested.Uncertainties in material properties and erosion criteria, omission of contact friction, and possible artifacts of far-field boundary conditions may adversely affect accuracy or precision of results.Another likely source of model discrepancy is the thicker, more compliant polymer representation than that tested experimentally, which would tend to promote target defeat for reasons explained above.
In summary, numerical results listed in Table 2 and shown in Figures 4 and 5 demonstrate how resolution of geometry and behavior of thin interfaces between layers of stiff material in armor systems strongly affects predicted ballistic efficiency.It follows that representation of interfaces should be carefully considered by the numerical analyst when constructing finite element or finite difference models for performance evaluations of such systems.Concurrent experiments and validation simulations on systems of lower complexity are recommended for future work, such that sources of discrepancy between model and experiment can be more precisely identified.Cohesive zone representations of interfacial separation [23][24][25] offer the potential for more realistic modeling of interfacial physics than the fully bonded or free surface interactions prescribed herein among layers.
Constitutive models with a more rigorous basis in finite deformation kinematics [26] and thermodynamics [21] may enable improvements in descriptions of the bulk behavior of metals [23,27] and ceramics [28], albeit at increased model complexity and computational expense.Phase field models [29] of structural transformations (e.g., for high-pressure phase transitions [17] and fracture in AlN) and nonlocal models for inelasticity and damage mechanisms [30] may also offer improvement over usual continuum mechanical treatments available in simulation codes such as EPIC, for example, potential benefits with regard to regularization of numerical solutions.
Experiments in [1] represented by simulations in the present numerical study do not address potential failure mechanisms observed in all possible kinds of ballistic impact problems.For example, adiabatic shear banding, plugging, and/or petal formation in metallic targets (often thin) reported for other armor systems [14,[31][32][33][34][35] are not of primary interest in the present case.Plasticity, fracture, and solid-solid phase transitions are appropriately addressed here for aluminum nitride, but mechanisms prevalent in other brittle targets not relevant here include pore collapse, for example, in concrete targets [34,36] or impacted rocks and minerals [37,38] and stress-induced amorphization, as observed in boron carbide [38,39] and quartz [40].For these different classes of targets not considered herein or in [1], appropriate constitutive models should always be chosen or constructed to represent dominant failure mechanisms observed in corresponding experiments.

Conclusions
Numerical simulations of ballistic impact and penetration of targets consisting of layers of aluminum nitride ceramic tile(s), polymer laminae, and aluminum backing have been conducted over a range of impact velocities on the order of 1.0 to 1.2 km/s.Results for ballistic efficiency have been compared with experimental data.Predicted residual penetration depths often tended to exceed corresponding experimental values, though simulations and experiments both demonstrated a trend of decreasing efficiency with increasing impact velocity.The closest agreement was obtained when polymer interfaces of small but finite thickness were not explicitly resolved, suggesting that the model representation of such interfaces is overly compliant.Results emphasize the importance of proper resolution of geometry and constitutive properties of thin layers and interfaces in numerical evaluation of performance of modern composite protection systems.

Figure 2 :Figure 3 :
Figure 2: Penetration into bare aluminum backing material: (a) depth versus impact velocity for simulation and experiment [1]; (b) simulation result at impact velocity of 1030 m/s.

Table 1 :
Materials and constitutive models.

Table 3 :
Convergence results for mesh size and time integration.