1. Introduction
In modern localized warfare, protecting the lives of military and police personnel is of paramount importance. Personal protective equipment such as bulletproof vests plays a crucial role in reducing soldier casualties. Researching how to provide more effective protection when encountering light weapon attacks is of great significance. Bulletproof vests with high ballistic ratings are already capable of effectively stopping small to medium-caliber rifle bullets while ensuring that back deformation does not cause severe injuries to the body. Therefore, studying the interaction between individual protection and projectiles, as well as the features of back deformation, is key to understanding the protective response to impacts.
High impact resistance, high abrasion resistance, high chemical resistance, flexibility, and excellent ballistic properties [
1] make ultra-high molecular weight polyethylene (UHMWPE) fibers ideal for lightweight and high-performance body armor. Bulletproof undershirts are one of the applications of their flexible properties, and typically around 50 sheets of weftless fabric consisting of four layers of prepreg [
2] stacked together can be used as soft body armor to meet the requirements of NIJ Level II ballistic standards [
3]. Hundreds of layers of unidirectional prepregs are stacked together in a two-by-two orthogonal manner and hot-pressed to obtain UHMWPE laminates, which can be used as rigid inserts in composite protection. The laminate could be used as the flexible insert close to the human chest on composite bulletproof vests made of brittle ceramics and flexible laminates. The laminate can also be used as a helmet or to independently withstand the impact of handgun bullets.
For the study of the dynamic mechanical properties of UHMWPE orthotropic lay-up laminates, Shi [
4] investigated the dynamic mechanical properties of UHMWPE orthotropic lay-up laminates with different lay-up angles under high strain rate compression and found that the energy absorption capacity varied with the lay-up angle and that the main damage modes were delamination and compaction in the thickness direction. Asija [
5] utilized the SHPB test to evaluate the performance of the STF-treated UHMWPE composites at high strain rates and showed that STF-treated specimens outperformed untreated specimens in terms of peak stress, strain, strain rate, and impact toughness.
For numerical simulation studies on protective impact response, Long [
6] proposed an improved Taylor approximation method. This new method utilizes the recursive relationship of the inverse Langevin function as an auxiliary tool, enabling a more accurate description of the relationship between stress, strain, and displacement in materials. Compared to the widely-used 5th-order Taylor approximation, this new method offers higher accuracy and stability when dealing with simulations involving large deformations. Nguyen [
7] applied the equivalent laminate discrete method to accurately capture vertical failure. They simulated the impact of different thickness laminates, ranging from 12.7 mm to 102 mm, with calibers of 12.7 mm and 20 mm and velocities ranging from 400 m/s to 2000 m/s. They predicted penetration mechanisms and bulge morphology. Bürger [
8] proposed a strain rate constitutive model to simulate the ballistic impact response of UHMWPE fiber-reinforced composites and predict energy loss. However, the overall modeling failed to simulate laminated damage and bulge height.
In terms of research on theoretical models for composite materials, Long [
9] investigated the influence of temperature on critical energy release rates using damage mechanics material models and element deletion methods. By considering factors such as strain rate, temperature, and stress state, they simulated a three-dimensional fracture specimen using an advanced material model based on test data to find critical energy release rates at different temperatures. The simulation results showed that critical energy release rates increased with temperature. Numerically, the equivalence conditions in brittle and small-scale plastic fracture were proven [
10]. Vaziri [
11] proposed a model to predict the transient response of composite plates under non-penetrating impact. This model used two-dimensional elements to simulate the composite material but could not predict internal phenomena. However, it could predict the structural impact response. Analytical models for studying impacts only provide very limited solutions. This leads to the need for a large number of different analytical models for different scenarios. For example, Ben-dor [
12] compiled around 280 analytical models, with 20 categorized as the most cited models.
The 3D digital image correlation (DIC) measurement and analysis method was first proposed by Luo [
13]. Its basic principle combines the binocular stereo vision principle with digital image correlation matching technology to reconstruct the three-dimensional spatial coordinates of each point on the surface of the object before and after deformation. This enables the determination of surface morphology and three-dimensional deformation information. Compared to other optical measurement methods, 3D-DIC offers advantages such as non-contact, full-field measurement, automation, simplicity of optics, universality, and strong resistance to interference. As a result, it has been widely used in the mechanical property testing of various materials across different fields [
14]. Currently, 3D-DIC technology is relatively mature. Freitas [
15] and Bigger [
16] have used high-speed photography with 3D-DIC to measure the penetration process of bullets of different calibers on various protective materials. They also calculated and analyzed the corresponding deformation. With the growing utilization of high-speed photography in the research of dynamic mechanical properties and impact response, 3D-DIC technology is gradually being applied in the mechanical property and impact testing of UHMWPE laminates.
In the research of impact response in fiber-reinforced composite laminates, the focus is primarily on investigating the influence of laminate thickness, impact velocity, and stacking sequence on penetration characteristics, local damage, and stress wave transmission. While 3D-DIC technology assists in capturing impact features, there is a need for further research on the response of thin laminates and analysis of strain distribution. Studying the strain distribution in UHMWPE laminates contributes to understanding the impact mechanisms and improving numerical models. Numerical simulations primarily aim to predict the dynamic mechanical properties of fiber-reinforced laminates, particularly the consistency of stress–strain curves with experimental data. However, global modeling fails to fully represent delamination phenomena, and the dynamic compressive behavior of interlaminar cohesive strength requires further investigation. Discrete modeling is suitable for laminates of varying sizes and thicknesses. Despite progress made by numerical models in studying back bulging, there is still a need to improve the simulation accuracy of strain distribution and compression wave propagation to deepen the understanding of deformation and damage mechanisms.
To address the limitations in current research, this study investigates the damage mode of UHMWPE laminate under typical pistol bullet impact and quantitatively analyzes the bulge morphology and strain distribution on the back surface using 3D-DIC technology. An optimized numerical simulation model is proposed, employing interlayer contact based on the same traction-separation criterion instead of cohesive units. This is because the large size of the laminate would require a significant amount of time when using cohesive elements. Although, since cohesive effects were not simulated with a mesh, it is difficult to determine the exact moment when interlaminar separation occurs. It has little influence on the focus of this study, which is more on the bulging and the morphological changes of the sub-laminates during the impact process. Therefore, this model is applicable for studying the impact response of laminates on a larger scale while appropriately reducing computational complexity without compromising accuracy. The validity of the numerical model is verified through comparison with experimental results, and further analysis is conducted to explore the interaction characteristics between the bullet and laminate, as well as the variations in energy. The experimental findings provide essential quantitative data for investigating the impact behavior of UHMWPE laminates, while the SPH-FEM-coupled numerical simulation method offers valuable insights into the interaction between pistol bullets and laminates.
4. Numerical Model and Results
4.1. Finite Element Models
To investigate the main features observed in the experiments, numerical simulation models were tested for the impact response before 400 μs. The simulation employed a full-model numerical simulation with symmetric constraints.
Figure 17a shows that the laminate has a denser grid in the impact region, and the overall laminate is simulated using equivalent sublaminates. Along the thickness direction (z-direction), the laminate is divided into 25 equivalent sublaminates, with 2 layers of grid points in each sublaminate. In the x and y-axis directions, there are 51 points on each face, resulting in 2500 elements per face and 5000 elements per equivalent sublaminate. To maintain accuracy and reduce computational cost, the contact relationship between adjacent equivalent sublaminate was handled using an automatic node-to-face contact based on traction-separation criteria, with contact being disconnected when contact option 9 was selected. An interface spacing of 0.001 mm was also set between adjacent equivalent sublaminates to avoid penetration. The bullet was simulated using the Smoothed Particle Hydrodynamics (SPH) method [
20], as it is suitable for large deformations and does not require grid removal, allowing for simulating the bullet dispersion.
Figure 17b shows the full model of the bullet. Contact between the bullet and the laminate was implemented using erosion node-to-face contact, with the soft 1 option opened with a coefficient of 0.2. Additionally, for boundary conditions, displacement in the x, y, and z directions was restrained for a 4 × 4 element grid on each corner of every equivalent sublaminates.
For the material model of the sublaminate, the MAT COMPOSITE FAILURE SOLID MODEL was used. Due to the correlation between the bulge shape and the shear modulus, as well as the correlation between the number of penetrated layers and the shear strength, the shear modulus Gab and shear strength were adjusted based on the height, width, and number of penetrated layers of the bulge observed in the experiments. Please refer to
Table 3 for the specific material parameter values. Due to the laminated structure of the laminate, it can be considered orthotropic [
5]. Therefore, the material properties in the 0°and 90°directions can be assumed to be the same. The commonly used Johnson–Cook material model was used to simulate the metallic materials [
21]. Copper and lead were used for the bullet jacket and bullet core materials, respectively, and the Johnson–Cook material model was applied. Please refer to
Table 4 for the specific material parameters.
To simulate the debonding process of the adhesive layer [
23], the CONTACT_AUTOMATIC_SURFACE_TO_SURFACE_TIEBREAK contact algorithm was used between the equivalent sublaminate. In this contact algorithm, Option 9 was selected as a failure criterion. In this criterion, EN and ET represent the initial slope before normal and tangential separation, respectively, with units of stress/length. T and S represent the ultimate stresses in the normal and tangential directions, respectively. ΔN and ΔS represent the complete separation distances in the normal and tangential directions. Additionally, GIC and GIIC represent the energy release rates in the tangential and normal directions, respectively, with units of stress∙length. To ensure that the separation distance is greater than the unloading distance [
17], the values of EN and ET should be sufficiently large. Please refer to
Table 5 for the specific numerical values [
24,
25,
26].
4.2. Numerical Validation
The temporal evolution of the bulge bottom shape on the back surface, as obtained from numerical simulations, is shown in
Figure 18. By comparing
Figure 7 and
Figure 18, it can be observed that the numerical model results are consistent with the experimental findings. Prior to 220 μs, the bottom shape of the bulge was diamond-shaped in both the simulations and experiments. In the experiments, the bulge bottom gradually transitions into a square shape, starting at 320 μs. The observed trend matches the experimental results.
Figure 19 depicts the evolution of the equivalent stress on the back surface. The overall trend is consistent with
Figure 16. The stress wave is transmitted to the back surface at 40 μs after the initial impact. It can be observed that initially, the maximum equivalent stress is concentrated in the center and has a circular shape. At this stage, only a small region in the center experiences stress. Then, the equivalent stress along the x and y directions appears, gradually expanding towards the boundaries. During this process, the range of the maximum equivalent stress in the center also increases, and its value becomes larger. After 60 μs, the range of the equivalent stress in the center continues to expand, but the numerical values decrease. Additionally, the range of equivalent stress along the x and y directions gradually decreases. During this stage, both the bulge and the fiber direction experience stress. By 120 μs, the equivalent stress along the x and y directions disappears, while the range of equivalent stress in the center continues to expand and maintains a circular shape, with decreasing values. During this stage, the fiber direction no longer experiences stress, and the stress is mainly concentrated in the bulge. Subsequently, the center region of the equivalent stress gradually forms a diamond shape. By 420 μs, the equivalent stress magnitude on the four edges of the diamond shape surpasses that in the center. During this stage, the boundary region of the bulge experiences stress.
Due to the inability to measure stress in experiments, only the stress distribution in the numerical model is analyzed. The analysis of stress results in the numerical model is showed as
Figure 20.
The maximum and minimum values of x-stress are 1.828 × 10
3 MPa and −1.507 × 10
3 MPa, respectively, while for y-stress, they are 1.729 × 10
3 MPa and −1.216 × 10
3 MPa, respectively. It can be observed that the drumhead experiences tension in the middle, while regions near the x and y axes undergo compression at the boundaries. The stress distribution resembles the strain distribution depicted in
Figure 12 and
Figure 13. The maximum and minimum values of z-stress are 9.87 × 10
2 MPa and −9.08 × 10
2 MPa, respectively. At the center of the impact site, the stress is zero, with the drumhead sides experiencing tension and the boundaries experiencing compression. This indicates a certain expansion in the z-direction on the sides of the drumhead and compression at the boundaries in the z-direction. The maximum and minimum values of xy-stress are 6.50×10
2 MPa and −6.72 × 10
2 MPa, respectively. There is no xy-stress at the center of impact and along the x and y axes. The drumhead sides exhibit negative values, suggesting stretching along the 45° direction, while the boundaries show positive values, indicating stretching along the −45° direction. The directions of stress in the second and fourth quadrants are opposite to those in the first and third quadrants. The positive and negative values of strain are depicted in
Figure 21.
The temporal evolution of the bulge width in the x and y directions, as observed in the numerical model and experiments, is shown in
Figure 22. From 0 μs to 150 μs, both the experimental and numerical model results demonstrate a rapid increase in width. During the period from 150 μs to 300 μs, the rate of width increase slows down. However, in the experiments, the width of the bulge increases significantly after 300 μs, gradually transitioning into overall displacement. Correspondingly, in the numerical model, the rate of width increase after 300 μs also becomes larger, although it may not be as pronounced, still exhibiting a trend towards overall displacement.
The specific values are shown in
Table 6. The numerical model and experimental results show minimal differences in the bulge width. The maximum relative error in the x direction is 24.12% and in the y direction is 21.29%. The minimum errors in the x and y directions are 1.49% and 0.85%, respectively. Additionally, both the experimental and numerical model results indicate that the width in the x direction is very close to the width in the y direction, demonstrating the orthotropic consistency of the material in the x and y directions.
4.3. The Morphology of the Bullet Impact Process
The type of erosive damage to the bullet is stubbing. The type of erosive damage to the laminate includes penetration of the first three layers, partial damage of the fourth layer, and bulging expansion of the remaining parts. According to
Figure 23, it can be seen that in the numerical simulation, the front four layers of the equivalent sublaminate are penetrated with a thickness of 1.35 mm, compared to the experimental result of 1.03 mm. The error is 25.18%, indicating that the numerical model is capable of replicating the observed bulge phenomenon in the experiment with reasonable accuracy. In the numerical model, the average diameter of the bullet hole along the x-axis is 7.78 mm, while along the y-axis it is 8.18 mm. Compared to the experimental data, there is an error of 19.23% in the x-axis direction and an error of 24.53% in the y-axis direction for the bullet hole diameter. Detailed data can be found in
Table 7. A noticeable discrepancy between the numerical simulation and the experimental results is that the penetrated equivalent sublaminate does not attach together, as shown in the experiment, but exhibits separation. This separation may be related to the damage characteristics of the interlayer bonding force model used. In the experiment, the interlayer bonding force near the penetrated bullet hole is not completely lost, while in the initial penetration stage of the numerical simulation, the layers have already shown signs of separation, indicating an early failure of the interlayer bonding force. To reduce this interlayer separation, it may be considered to increase the distance parameter for complete separation between the layers, in order to better simulate the actual interlayer bonding behavior.
Figure 24 illustrates the morphological changes of the bullet during the penetration process. Initially, the bullet head is stubbed at the beginning, and as the erosion process progresses, its front end gradually forms a mushroom shape, while the material at the rear end is compressed and concentrated forward. As the penetration continues, the front end of the bullet shell begins to fracture and disperse outward, while the material at the rear end diffuses around. Eventually, as shown in
Figure 24c, the front-end fragments of the bullet core are scattered, but the main body remains concentrated near the impact point. The front end of the bullet shell fractures and disperses under impact, while the rear end forms several larger chunks.
The simulated results are very similar to the observed deformation of the bullet core in the experiment, both exhibiting a mushroom shape. Due to the softness of the lead material, the main body of the bullet does not completely disperse. The bullet shell undergoes an unfolding process in both the simulation and the experiment. However, in the experiment, the bullet shell unfolds and splits into multiple small pieces, scattering and embedding into the laminate. In contrast, in the numerical model, the unfolded bullet shell only forms a few small pieces at the edges and does not completely fracture. This may be due to the bullet’s failure to fully penetrate the fifth layer, leading to partial penetration, and the fifth layer providing support to the bullet shell, restricting further unfolding and dispersal. Since the lateral splitting of the bullet was not restricted by the fifth layer, as shown in
Figure 25, the bullet shell completely fractures into multiple small pieces and disperses in all directions after unfolding. Thus, it can be seen that further research on the interaction between the bullet and the laminates is needed to obtain more accurate simulation results.
Nguyen’s [
17] study pointed out that under the condition of impact velocity less than 500 m/s or laminate thickness less than 10 mm (in this study, the laminate thickness is 11 mm), the damage pattern does not include the shearing and stamping stage but only includes partial deformation penetration and the final bulging but not fully penetrated expansion stage, with partial deformation penetration being considered a transition stage. The specific details of the penetration process are shown in
Figure 26, where the bullet starts to make contact with the laminate at 30 μs when it is a certain distance away from the laminate. During the time period from 40 μs to 100 μs, two significant phenomena were observed: first, a large bulge was formed on the laminate; second, there were noticeable traction marks on the edges of the laminate and along its x and y axes. As shown in
Figure 26a, initially, the bullet compresses the laminate. Once compressed to a certain extent, under the action of compression along the thickness direction and around the edge, three sub-layers of the laminate are penetrated, as shown in
Figure 26c. After penetrating three sub-layers, the bullet continues to move, causing local damage to the first sub-layer but not penetrating it. Subsequently, the penetrated three sub-layers have elastic recovery and separate from the unpenetrated parts. The unpenetrated parts continue to move with the impact of the bullet until the highest point is reached. The final shapes of the impact face and the back face are as shown in
Figure 26d.
At 100 μs, partial fiber layers were observed to undergo penetration, indicating the occurrence of significant deformation penetration process. Subsequently, after 100 μs, the penetrated fiber layers began to partially recover and separate from the non-penetrated parts. Meanwhile, the non-penetrated parts continued to move under the propulsion of the bullet until approximately 400 μs, when the bulge reached its maximum height. This phenomenon corresponds to the characteristic of a protruding bulge that is formed but not penetrated. The experimental results also exhibit two main features: firstly, penetration occurred in approximately 1.35 mm of the fiber layers, with traction observed in both the x and y directions; secondly, the formed bulge ultimately remained unpenetrated.
4.4. Energy Transformation
The energy distribution throughout the entire experimental process, as shown in
Figure 27, reveals an initial kinetic energy of 433.4 J. As shown in
Figure 27b, Before 150 μs, there is a rapid decline in the bullet’s kinetic energy, which is subsequently transferred to the laminate. At this point, the bullet’s kinetic energy is converted into both kinetic energy and internal energy within the laminate, leading to an increase in the laminate’s energy due to energy deposition. After 150 μs, the energy within the laminate marginally exceeds the changes in the bullet’s internal energy, likely attributed to sustained energy absorption and dispersion within the laminate’s larger volume. The variation in the bullet’s energy stabilizes around 200 μs, while the laminate’s energy remains at high levels, indicating a steady state in the late stage of energy transfer and absorption.
At the initial moment, the laminate begins to accumulate energy, possibly due to the formation and propagation of stress waves. By 71 μs, kinetic energy absorption within the laminate peaks at 152.18 J, subsequently trending towards stability, while internal energy absorption within the laminate continues to rise. This suggests potential structural changes or damage within the laminate. At this juncture, the bullet’s kinetic energy stands at 84.4 J, reflecting an 80.64% reduction, with the laminate absorbing 152.18 J of kinetic energy and 133.84 J of internal energy, while the deposited energy within the laminate amounts to 45.9 J. This further confirms substantial energy absorption by the laminate during the early stages of impact, in accordance with the principle of energy conservation.
At approximately 150 μs, the energy within the laminate begins to stabilize, and the energy curves of the system become parallel. This indicates a reduction in the rate of energy transfer, with the laminate likely dispersing energy throughout its entire structure. Subsequently, until 200 μs, the bullet’s kinetic energy remains nearly constant, indicating a stable state of energy transfer and absorption during the late stages of impact, while the energy within the laminate continues to sustain relatively high levels.
The first four equivalent sublaminates absorb a significant amount of energy and exhibit pronounced material erosion under the impact of the bullet. Additionally, compared to other layers, the mesh deformation in these first four equivalent sublaminate layers is more pronounced, possibly due to the propagation of stress waves and energy concentration during the bullet impact.
As shown in
Figure 28a, eroded internal energy is only present in the first four equivalent sublaminates, with no observed erosion in subsequent equivalent sublaminates. This indicates significant mesh deformation in the first four equivalent sublaminates due to erosion. This is because the first four sublaminates have been eroded and damaged enough to become bullet holes. The meshes in these areas become distorted and are deleted. Therefore, the erosion energy is present in these four sublaminates. The energy in the fourth equivalent sublaminate differs from the preceding layers, steadily increasing. At 75 μs, eroded energy stabilizes in the first three layers while indicating more significant variations in the fourth layer, consistent with the trend of the bullet’s kinetic energy dropping below the laminate’s internal energy observed in
Figure 27.
According to
Figure 28b, the internal energy of the first four equivalent sublaminates increases over time, peaking at approximately 200 μs before stabilizing. This suggests significant energy absorption in the initial stages by the first four equivalent sublaminates, consistent with the observation of the bullet’s energy depletion around 200 μs in
Figure 27, indicating the completion of major energy conversion.
In terms of overall energy distribution, the hourglass energy constitutes only a small fraction, indicating that mesh distortion energy in the numerical model is not the primary source of energy dissipation, thereby confirming the accuracy and reliability of the numerical model.