In-Silico Model Validation of Impact on a Composite Helmet Shell

High-velocity impacts to the helmeted head remain a significant source of injury and death despite advances in materials and design [1]. Composite helmets have improved protection by reducing the risk of projectile penetration, but have introduced a different injury mode in the process. Helmets used by military or law enforcement, e.g., prevent ballistic penetration by absorbing the projectile energy through deformation, absorbing the energy of the projectile over a larger distance [2–4]. If the projectile energy is high enough, the helmet will deform inward to strike the head, transferring the remaining energy onto the head and causing a class of injury known as behind helmet blunt trauma (BHBT) [5–7]. BHBT encompasses a large group of head injuries, including skull fracture, hematoma, and traumatic brain injury [8] caused by blunt force trauma due to large deformation in the helmet shell. Accurate quantification of the deformation in such helmet systems is critical to preventing BHBT. Finite element (FE) simulations are particularly useful for identifying how the helmet and the head interact during impact, since the use of experimental approaches is limited by safety challenges and speed of data acquisition. However, FE simulations of high-velocity impact require careful validation against existing experimental data. To accurately capture the underlying mechanics that cause BHBT, a critical step in simulation is to examine whether the predicted deformation of FE helmet models closely matches that of their physical counterparts.

Modern composite helmets generally consist of a layered shell of woven Kevlar fibers in a thermoset resin matrix [2,3,9]. During a projectile impact, fiber composite helmets deform to a failure threshold at which the Kevlar fibers begin to break. As the helmet further deforms and an increasing number of fibers break, the strength of the composite around the impact site decreases [3]. Strain-softening due to fiber breakage becomes an important consideration during the numerical modeling of composite materials since the strength of the material depends on the local strain. In high-velocity impact problems where there are large and highly localized strains, the performance of the numerical composite model is highly dependent on the ability to accurately capture such softening behavior. Because the calculated local strain depends on the element size in FE models, the material failure of the composite model depends on the discretization of the domain. As a result, a simple convergence of the results with increasing mesh density, seen in simpler material models, does not exist for composites undergoing failure. Instead, the mesh size must be tuned until the simulation performance matches the experimentally observed results using the same material [10].

Previous numerical studies of impact mechanics on composite helmets have been performed with various justifications for their mesh selection. Most commonly, variable mesh sizes provide detail at the impact site while saving computational costs elsewhere on the helmet [11–16]. However, due to the strain-softening or damage models incorporated into many such studies, a variable mesh size will result in areas of the helmet with a coarser mesh exhibiting a stiffer material response relative to areas with a finer mesh. Other researchers have used a consistent mesh size, and selected the size based on balancing resolution with computational efficiency [17,18]. Tan et al. selected a mesh size of 4 mm based on such criteria and their results showed differences of up to 60% between the deformation results of their experiments and their simulations. In contrast, Rodríguez-Millán et al. performed a calibration of mesh and material properties by comparing the results of their simulations to experimental tests of ballistic impacts on a composite helmet [10]. As a result of the calibration performed, they achieved a much closer agreement between their numerical and experimental results. Jazi et al. used simplified elasticity models to represent helmet behavior [19]. Such simplified models eliminate the convergence and stability challenges but introduce their own limitations associated with neglecting the underlying damage behavior of the composite.

Due to the curvature of the helmet, particularly across the front, the incident angles between the projectile and the helmet depend on the location of the impact points. The incident angle affects the ratio of energy absorbed by the helmet versus the energy retained by the projectile through residual velocity or plastic deformation [12]. The ballistic impact simulations in the studies discussed were setup to mimic the experimental test standard for composite helmets used by military and law enforcement. For example, in the United States, the experimental test standard for projectile impact on composite helmets is defined by the NIJ-STD-0106.01 standard [20]. In experimental testing, precise control of the impact point of the projectile is limited. As such, the test standard defines a large region in which an impact to the helmet is considered valid. Projectile hits are considered valid for the test if they impact the front of the helmet less than 90 mm above the basic plane, defined as passing through the center of the ear openings to the lower edge of the eye socket on the head, and not more than 50 mm from the mid-sagittal plane. Similar profiles are also defined for side and rear impacts. In contrast to experimental testing, the exact impact point in the FE simulations can be precisely controlled, but the point of impact is often not explained in previously published studies.

Such studies present a wide range of reported deformation of the helmet system. It is not clear how their different selections of mesh size or impact position may contribute to the different responses observed since the selection criteria for both factors are not fully explained. Such limitations make it difficult to compare the results of existing studies with each other or to use the results they generate as validation of future research efforts. The goal of our work presented herein was to clarify a validation process for a composite helmet model. Subsequently, we evaluated the effects of impact position on the results of helmet impact simulations.

Due to the strain-softening nature of the composite materials used in helmets, the behavior of the helmet shell does not converge with mesh size as seen with most simple linear elastic models. Instead, the mesh size must be calibrated along with the material properties used against experimental data. A simplified case of a 9-mm full metal jacket (FMJ) projectile striking a flat composite plate was used to determine the effects of mesh size on the deformation of the composite. The back face deformation (BFD), defined as the deformation of the face of the plate opposite the impact, was used as the measure of composite performance since it is the metric used to qualify armors for protective use in the NIJ test standard [20]. The BFD calculated in these simulations was compared with experimental data to determine an appropriate mesh size that would align the predicted values of the simulated composite with the experimental counterpart.

The procedure for the plate-impact simulations was based on previous work performed by Rodríguez-Millán et al., who performed experimental testing of ballistic impacts on a flat composite plate [10]. In their study, a 400 mm square plate constrained at the edges was impacted with a 9-mm FMJ projectile at approximately $425±15m/s$⁠. The BFD of the plate was measured at approximately 18 mm for the first impact, which was used as the experimental performance baseline for the following analysis.

Three meshing parameters were selected for investigation: the element size across the surface of the plate, the number of layers of elements through the thickness of the plate, and the element size of the projectile. Each parameter was varied separately to determine its effect on the BFD of the composite. A summary of the parametric test conditions can be found in Table 3. The FE meshes for the projectile and plate, and later the helmet, were generated using eight-node hexahedral elements. One-point constant stress Gaussian integration (ELFORM 1) was selected to reduce computational cost and improve performance at large deformations. The use of under-integrated elements introduces the risk of non-physical shear locking, or hour-glassing [27]. Tests were performed to evaluate the effectiveness of the different hourglass controls available in LS-DYNA. The Flanagan–Belytschko viscous hourglass control (IHQ 2) [28,29] was found to create the greatest reduction in hourglass energy without loss of simulation stability. IHQ 2 was applied to the plate, projectile, and helmet for all subsequent analyses.

The contact relationship between the projectile and the plate was controlled by surface-to-surface contact. Excessively deformed elements that failed according to the above criteria were removed from the analysis. The contact algorithm used deleted the excessively distorted elements near the impact point and adjusted the contact surface as the elements at the crash-front were eroded to improve the simulation stability. All simulations were performed by transient analysis using explicit time integration. The time-step size was selected automatically and adjusted on each time-step based on the size of the smallest element in the simulation [29].

The results of the parametric analyses performed are shown in Fig. 2. In Fig. 2(a), as the mesh density increased, the resultant BFD of the plate increased as well. Below an element size of 1.5 mm, the weakening effect of the mesh refinement resulted in full penetration of the composite plate. The strain-softening effect caused excessive element distortion, which resulted in transient simulations that did not converge. Figure 2(b) shows the effect of the number of elements through the thickness of the shell on BFD. The BFD increased initially when changing from two to six layers of elements. Subsequently, an increased number of element layers resulted in a decrease in BFD. Similarly, the maximum BFD for varying mesh sizes on the projectile is shown in Fig. 2(c). At a mesh size of less than 1.25 mm, the projectile mesh had a negligible effect on the performance of the composite.

The results of the parametric analysis we performed were used to select the appropriate parameters for the helmet model. Based on the results produced by Rodríguez-Millán et al. [10], the smallest stable mesh size of 1.5mm which would approximate the BFD observed in the plate impact experiments they performed. As such, a mesh size of 1.5mm was used to generate the subsequent helmet model.

We employed the meshing parameters established in Sec. 2 and prepared a finite element model of a composite helmet shell [30], as shown in Fig. 3. The helmet mesh was generated symmetrically using a semi-automated meshing process. First, a shell mesh of one half of the helmet was generated. The shell model was then extruded into three dimensions and mirrored to create the full helmet model. The material properties of the helmet were the same as those of the composite plate listed in Table 1. Since the material model used for the helmet (MAT-054) is orthotropic, it was necessary to define the primary direction of fiber orientation. Due to the complex curvature of the helmet, it was not feasible to individually define the primary and secondary fiber directions for each element of the mesh. Instead, the simplified method developed by Lee and Gong was used [13]. The helmet was divided into five subsections, as shown in Fig. 3. For each section, the fiber direction was defined by a constant vector relative to a global coordinate system. The helmet decomposition was selected to minimize the angle between the curvature of the helmet and the global vector for that section, minimizing any weakening effect that might have occurred because of simplifying the fiber orientation to a uniform vector.

The effects of impact position on the helmet were tested using two concentric nine-point grids, as illustrated in Fig. 4. The location of the center of each grid was selected so that all points on the grid were considered valid hits according to the experimental test standard [20]. The center of the front impact grid was centered left-to-right on the helmet and located 30 mm up from the front rim of the helmet. The center of the side grid was centered front-to-back on the helmet and 95 mm from the bottom rim. The back grid was centered left-to-right and located 60 mm from the bottom, and the top grid was centered front-to-back and left-to-right on the helmet. The projectile was aligned with the global coordinate axes for each impact direction. For front, back, and side impacts, the projectile was oriented in horizontal plane and parallel to the x- or y-axis. For top impacts, the projectile was oriented parallel to the z-axis. The angle of impact was not changed for individual points on each grid to mimic the experimental test standard, where all projectiles are fired from the same position and orientation. The projectile velocity was set to 360 m/s to conform to the NIJ test standard [20].

The time progression of helmet deformation shortly after impact for a representative case of a frontal impact is shown in Fig. 5. The initial impact of the projectile created a large deformation around the impact site. Subsequently, as the elastic energy around the impact site rebounded the helmet shell, the rim of the helmet deflected inward and eventually overshot the deformation at the impact site. As a result, the point of maximum helmet deformation was found on the rim of the helmet, rather than at the impact site. Such an effect becomes more pronounced as the impact location moves closer to the helmet rim.

The resultant maximum BFD for back, front, side, and top impacts on the helmet can be seen in Fig. 6. Note that the asymmetry observed in the results is attributed to the numerical round-off and small variations in the contact algorithm in this highly non-linear problem. Otherwise, it was seen that small differences in impact position resulted in a significant change in the helmet deformation. The positional change in deformation was greatest on the front of the helmet, where the BFD increased by 35% when the impact point was moved down the front of the helmet by 20 mm. The effect was less pronounced on the side and back of the helmet and negligible on the top. Figure 7 shows a comparison of the mean and spread of observed BFD for different impact directions between the present study and experimental results obtained from impact tests performed on 48 composite helmets [31]. The numerical model developed herein displayed a lower spread in observed BFD for all impact directions. The numerically calculated BFD was lower for front and back impacts than the experimental results, but side and top impacts show similar behavior between the experimental and numerical results.

The simulations of projectile impact on plate show that the mesh has a significant effect on the prediction of the performance of the composite due to the failure related material softening. Smaller element size resulted in greater softening of the material around the impact location, as these high-velocity impacts exhibit large deformation focused on a small area. Below a certain element size, the projectile fully penetrated the plate, which was not observed in the experimental analysis [10]. Moreover, below a critical element size, the simulation was not able to achieve a stable solution. To select the appropriate mesh size for the helmet, the smallest mesh size where full penetration was not observed in a stable manner was used.

Using the smallest stable mesh size that did not experience full penetration, we still observed that the results of the plate-impact simulations underestimated the experimentally measured BFD [31]. Other researchers have previously investigated adjusting the properties of the material, rather than the mesh, to align the performance of their models with experiments [17,32]. In both of their studies, the through-thickness shear moduli were adjusted to calibrate the model performance. However, it is not clear whether the impact points were matched between the experimental and numerical results. The present study focuses only on how meshing decisions affect the performance of numerical helmet models.

In the initial plate-impact studies, it was found that the composite responded in an artificially stiff manner with a low number of elements through the thickness of the composite. As the number of element layers increased, there was an increase in the BFD followed by a consistent decrease at higher layer counts. The inclusion of more element layers helped to converge toward accurately representing the number of composite laminae in the helmet, whereas thicker elements resulted in an overly weak material response by acting more like a bulk material. The number of element layers had a smaller effect on the helmet deformation compared with changing the mesh density across the plate surface. Regardless, we conclude that the two mesh density criteria of layers should be adjusted concordantly to achieve the desired composite performance.

In the helmet impact simulations, the maximum BFD of the helmet varied largely depending on the impact direction and location. The BFD observed for side and top impacts was closely in agreement with previous experimental data [31], while the front and back impacts underestimated the expected BFD. The spread of the simulated deformation was less than the reported experimental data for front, back, and side impacts as seen in Fig. 7. The larger spread in experimental data may be due to less control over the impact conditions in an experimental setup when compared with the present simulations.

For front, back, and side impacts on the helmet, the BFD increased as the impact point moved lower on the helmet. The increase was more apparent at the front of the helmet than on the sides or back, likely due to the sharper curvature and the reduced total area available for energy absorption. Since an increase in BFD is observed as the impact location moves closer to the rim of the helmet, part of the increase in BFD may also be attributed to the impact point being closer to that unsupported boundary. Such a conclusion is supported by observing the progression of helmet deformation, shown in Fig. 5.

The large spread in the BFD calculated at different impact points shown in Fig. 6 demonstrates two critical limitations to compare results across different impact studies: small experimental sample size and lack of impact position information. In previous experimental studies impacts on a small number of helmet samples were examined [10,17,33]. Such studies have not provided sufficient description of the impact point to enable accurate simulated reconstruction. We have shown through the analysis in Fig. 6 that changes in impact position as small as 5 mm can result in significant alterations of the BFD. As such, generating an accurately simulated recreation of a single impact test from a limited description of the impact point without precise measurements of its position is not feasible. Our results show that a wide range of different helmet deformations could be mimicked by making small adjustments to the impact point until the BFD matches an experimentally derived result. Consequently, model validation by reconstruction of single experimental impact cases is unreliable.

The consequences of mesh size discussed in the present study are applicable to any impact modeling of composite materials using similar strain-softening models. Additionally, the effects of impact position should be considered regardless of helmet model when comparing numerical performance against experimental results for similar impact problems.

In this paper we have described, for the first time, a set of best practices for calibrating a model of high-velocity impact on a composite helmet, considering the role of both the physically relevant material properties and the numerically based mesh characteristics. The performance of a numerical model of a composite system depends on the selection of both mesh and material properties. The material properties can be calibrated independently of the system geometry using well-established material testing methods. However, the present study shows that the mesh size cannot be picked arbitrarily or solved for by a traditional FE mesh convergence study. The stiffness of the composite varies significantly with the mesh size chosen, and stable solutions cannot be achieved below a certain element size. A mesh calibration must be performed to align the performance of the numerical model with its physical counterpart in addition to an accurate assignment of material models.

The exact position of high-velocity impact has a significant effect on the calculated BFD of the helmet. The position of projectile impact is rarely specified in the literature. Since the impact position so drastically affects the BFD, without exactly matching the impact point between two studies a comparison of the results cannot be drawn. As a result, calibration against small sample size experimental studies is insufficient to ensure proper correlation of numerical model performance since individual results can be replicated through small changes in impact position. Instead, calibration must be performed against a statistical distribution of deformation from a large sample size of experimental tests.

Funding for this project was provided by the United States Army Research Lab Soldier Protection Program Cooperative Agreement #W911NF2120208. The views presented in this paper are those of the authors and do not necessarily represent the views of the sponsor.

There are no conflicts of interest. This article does not include research in which human participants were involved. Informed consent not applicable. This article does not include any research in which animal participants were involved.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

$m$ =

Johnson–Cook thermal softening parameter

$n$ =

Johnson–Cook hardening exponent

$A$ =

Johnson–Cook yield stress

$B$ =

Johnson–Cook hardening parameter

$C$ =

Johnson–Cook strain rate sensitivity parameter

$I$ =

Cockcroft–Latham damage parameter

$T$ =

temperature

$ec$ =

Chang–Chang compressive fiber damage parameter

$ed$ =

Chang–Chang compressive matrix damage parameter

$ef$ =

Chang–Chang tensile fiber damage parameter

$em$ =

Chang–Chang tensile matrix damage parameter

$Sc$ =

shear strength

$Tr$ =

reference room temperature

$Tm$ =

melting temperature

$Wcr$ =

Cockcroft–Latham parameter

$Xc$ =

longitudinal compressive strength

$Xt$ =

longitudinal tensile strength

$Yc$ =

transverse compressive strength

$Yt$ =

transverse tensile strength

$a,b,c$ =

primary, secondary, and tertiary composite orientations

$β$ =

shear strength weighting factor

$ε$ =

strain

$ε˙$ =

strain rate

$ε0˙$ =

reference strain rate

$εf¯$ =

failure strain

$dε¯$ =

strain increment

$σij$ =

stress component in the $ith$ plane and $jth$ direction

$σmax$ =

maximum principal stress

$σY$ =

Johnson–Cook yield surface

BFD =

back face deformation

BHBT =

behind helmet blunt trauma

FE =

finite element

FMJ =

full metal jacket

J–C =

Johnson–Cook material model