Simulation-augmented learning for ballistic performance prediction of ultra-high-strength steel

Raw materials

The G33 ultra-high-strength steel is a low-cost, high-performance steel that can develop either martensitic or bainitic microstructures upon exposure to heat. Its chemical composition is presented in Table 1.

The measured phase transition temperatures are as follows: AC3 at 800 °C, AC1 at 644 °C, Ms at 250 °C, and Bs at 360 °C. These phase transition points were used as the basis for designing a “three upsetting and three drawing” process, which was implemented for the steel ingots. The specific steps are as follows:
(1) The initial forging temperature was set at 1,150 °C. The ingot was upset to half of its original height and then drawn back to its original height.
(2) The ingot was reheated to a temperature of 1,080 °C, upset again to half its height, and then drawn back to the original height.
(3) The ingot was reheated to 1,050 °C, upset to half its height for the last time, and then forged into rods measuring Ф15 × 1,000 mm.

Finally, the rods underwent annealing treatment at 680 °C to effectively relieve the internal stresses induced by the forging process.

A comprehensive range of G33 microstructures was captured for the training dataset by designing heat-treatment schedules that encompass the major bainitic and martensitic transformation pathways. Consequently, we selected specific combinations of the temperature and holding time. Adjustment of the cooling rate and tempering temperature enabled systematic control of the microstructural morphology and distribution, along with the size and dispersion of the precipitated carbides. These microstructural changes strongly influence the toughness, hardness, corrosion resistance, and overall mechanical response of the steel.

Bainitic ultra-high-strength steel was produced via a newly designed two-step heat treatment route consisting of austenitization and isothermal quenching, the latter of which was performed in a salt bath containing 50 wt.% of NaNO₃ and KNO₃. Austenitization was accomplished by increasing the temperature in the range from 900 to 1,100 °C in increments of 100 °C. Because the temperature and holding time of isothermal quenching strongly affect the resulting microstructure, upper bainite was targeted at 345 and 355 °C with a holding time of 30 h, whereas lower bainite was targeted at 245 and 255 °C with holding times of 10, 20, and 30 h. The steel specimen was cooled in oil in the final stage to reduce quenching stress and distortion. This bainitic route involved six key process parameters and ten procedural variants.

Martensitic ultra-high-strength steel was produced by adopting a heat-treatment route with up to four stages - austenitization, quenching, cryogenic treatment, and tempering. Austenitization was conducted from 700 to 1,100 °C in increments of 50 °C to systematically examine its effect on the microstructural evolution. Water, oil, and air-cooling were used to compare the different quenching conditions. Some schedules also included cryogenic treatment at -73 °C for 1 h to evaluate the effect thereof on the phase morphology and distribution. Tempering was performed over low- (150-200 °C), medium- (200-400 °C), and high-temperature (400-700 °C) ranges with holding times of 2, 4, and 5 h, respectively. Some routes included double tempering to broaden the range of tempering strategies and resulting microstructures. Overall, the martensitic route involved as many as 12 key process parameters and 115 different processing schedules. In addition, nine conventional annealing treatments commonly used for other metallic materials were included for comparison with the pearlitic structures obtained under these conditions.

In total, 134 heat treatment schedules were designed and implemented, which yielded nearly 3 t of test material. Details of the heat-treatment design and the expected primary microstructures are provided in Supplementary Table 1.

Mechanical property dataset derived from experiments

The mechanical properties that were used to construct the training dataset were acquired directly through experimental testing and comprise 13 basic mechanical properties. The specific parameters are detailed in Supplementary Table 2.

(1) Quasi-static compression and tensile tests were performed at room temperature
Quasi-static compression and tensile tests were performed on the 134 steel variants using an INSTRON 5985 universal testing machine. Testing adhered to GB/T 7314-2017 and GB/T 228.1-2010 standards at a strain rate of 1 × 10^-3 s^-1. Compression specimens measured Φ6 mm × 9 mm, and tensile specimens had a gauge section of Φ5 mm × 30 mm. Each test was repeated at least three times to obtain quasi-static stress-strain curves at room temperature.

(2) Dynamic compression and tensile tests at room temperature
Dynamic compression and tensile tests were conducted using split Hopkinson pressure and tension bar systems in compliance with the GJB 8374-2015 and GB/T 34108-2017 standards. The dynamic compression tests were conducted with strain rates of 3,000 and 5,000 s^-1, using samples sized Φ5 mm × 5 mm. The dynamic tensile tests were carried out at a strain rate of 2,000 s^-1, with sample gauge lengths of Φ2 mm × 5 mm. Each test was repeated at least three times to obtain the stress-strain curves for room-temperature dynamic compression and tension.

(3) Vickers hardness tests
The Vickers hardness values of the 134 steel variants were measured using a Wolpert 450SVD tester in the same laboratory. These tests were conducted under a 5 kg load, following the standard GB/T 4340.1-2009. Each material was tested 10 times, and the average Vickers hardness value was used as the final result.

In total, we standardized and collected 1,072 stress–strain curves across various loading conditions. The tests produced 1,742 measurements of the mechanical properties of G33 steel under different heat treatment conditions. These measurements were then visualized, analyzed, and used to construct the dataset. As shown in Figure 1, the resulting dataset spans a broad numerical range, shows a well-distributed variation, and has a high overall quality. These mechanical property measurements formed the input feature space for subsequent service-performance prediction.

Figure 1. G33 dataset overview: radial map of dynamic yield strength (3,000 s^-1) and SEMs of the lowest and highest specimens. SEMs: Scanning electron microscope images.

Service performance dataset generated via high-efficiency simulations

After establishing the mechanical property dataset, the next step was to obtain the corresponding service-performance labels. Because direct penetration testing for all 134 steel variants would be prohibitively expensive, a combined experimental–simulation strategy was adopted. G33 ultra-high-strength steel is primarily used in supersonic kinetic-energy penetrators operating at Mach 2-3; therefore, its service-performance evaluation is centered on penetration efficiency. For deep-penetration applications, two criteria are particularly important: whether the projectile maintains structural integrity after impact and whether it can retain sufficient penetration capability to reach the specified target depth. The essential service-performance parameters used for these evaluations are listed in Supplementary Table 3.

Scaled penetration experiment using simulated projectiles

The penetration performance of G33 steel was evaluated cost-effectively by performing scaled projectile penetration tests, considering that scaled models can reproduce the key penetration dynamics of actual projectiles at significantly lower costs. These tests also allow repeated testing under controlled conditions, which accelerates service data collection and supports the rapid optimization of both the structure and material.

Among the 134 steel variants, 30 representative samples were selected to cover a broad range of heat-treatment routes and microstructural states. For each category, at least one subgroup was selected to highlight steels with distinctly different performance characteristics to provide a comprehensive evaluation of the behavior of G33 steel under high-speed impact conditions. High-speed photographs of the projectile penetration process are presented in Figure 2A, where each sample underwent a minimum of three tests at different penetration speeds. A total of 95 penetration tests were completed. Figure 2B shows a selection of recovered targets and scaled simulated projectiles to illustrate the experimental outcomes.

Figure 2. (A) High-speed camera images; (B) Target plates and scaled projectile samples.

(1) Scaled simulated projectile and target
An ogive-nosed scaled projectile with dimensions of 200 mm × 200 mm × 20 mm was used to impact the naval steel target plates.

(2) Sub-caliber launch platform
A 30-37 caliber smoothbore cannon served as the launch platform for the scaled simulated projectiles, as illustrated in Supplementary Figure 1A.

(3) Velocity measuring device
In this experiment, the firing distance was set to 3 m and the firing angle was set to 0° to the normal. The initial velocity was determined using a foil interrupt target positioned 1 m from the gun muzzle, as depicted in Supplementary Figure 1B. Additionally, a high-speed camera was employed to validate the launch velocity of the projectile and document the interaction between the projectile and the target, including the post-penetration posture and velocity of the projectile, as shown in Supplementary Figure 1C.

(4) Target frame and projectile recovery device
The steel plate was secured using F-type clamps to ensure that the clamping width around the perimeter did not exceed 20 mm. Rubber pads were applied to the front and back faces of the clamped sections for cushioning, as shown in Supplementary Figure 1D. To retrieve the projectile and any resulting post-penetration fragments, a recovery mat was strategically placed on the after-effect target stand, situated 50 mm behind the steel target, as shown in Supplementary Figure 1E.

A schematic diagram of the scaled projectile penetration test is shown in Supplementary Figure 2.

Efficient acquisition method for penetration simulation parameters

The scaled penetration experiments provided reliable reference data but were insufficient in number for robust machine-learning training. Therefore, finite element simulations were introduced to expand the service-performance dataset, provided that the simulation parameters could be accurately identified. Numerical simulations are an effective tool for augmenting the service performance data. Although the scaled penetration tests provided reliable experimental data, most of the service-performance data used in this study were generated by simulations. Because the validity of these results depends on the accuracy of the simulation parameters, an efficient method is required to accurately compute these parameters for the different G33 microstructures. The finite element method (FEM) is widely used in engineering to solve boundary value problems and analyze the response of complex physical systems under service conditions. LS-DYNA is a widely used finite element package that includes models of approximately 300 materials and 16 equations of state. Once reliable simulation parameters are obtained, high-performance computing can be used to lower the experimental costs, shorten the computation time, and support the selection of materials.

Figure 3 details the parameter acquisition process for penetration simulations. By incorporating these parameters into the models and utilizing large-scale parallel computing, this approach shifts the acquisition of service performance data from costly experimental testing to efficient numerical simulation. This innovative approach, which requires data from only a limited set of physical tests to rapidly compute the service performance parameters for G33 ultra-high-strength steel, substantially augments the sample size of the service performance dataset. Such enhancement significantly bolsters the accuracy of predictions made by machine learning models.

Figure 3. Flowchart of the parameter acquisition and optimization method for numerical simulation of penetration tests.

The JC constitutive and failure model is widely used in impact and dynamic loading analyses because this model accounts for strain hardening, strain rate sensitivity, and thermal softening. Accurately optimized parameters are necessary to improve the fidelity of the penetration simulations.

To ensure that the simulations accurately captured the behavior of the scaled projectiles penetrating the steel targets, we focused on the deflection of the projectile. A full-scale simulation model was built to match the experimental dimensions by defining the initial velocity using the *INITIAL_VELOCITY_GENERATION parameters. Both the projectile and the steel target were modeled using the JC constitutive and failure equations (*MAT_JOHNSON_COOK). The target JC parameters derived from experimental fitting are detailed in Table 2. An option to control hourglassing (*CONTROL_HOURGLASS) was included to prevent the zero-energy modes associated with single-point integration. The contact between the projectile and steel target was defined as the eroding surface-to-surface contact (*CONTACT_ERODING_SURFACE_TO_SURFACE). Fixed constraints (*BOUNDARY-SPC_SET) were applied to the contact area between the steel target and target frame to simulate the support provided by the frame.

Table 2

JC constitutive and failure equation parameters for marine steel

Constitutive parameters	A	B	n	C	m
	0.0068	0.0055	0.2000	0.0200	1.0000
Failure parameters	D 1	D 2	D 3	D 4	D 5
	0.1580	0.2100	3.9700	0.0010	1.0000

JC: Johnson-Cook.

In addition to the constitutive and failure parameters, the penetration simulation requires auxiliary parameters. Notably, the PC value (keyword: *MAT_JOHNSON_COOK_PC: Failure stress or pressure cutoff) is an integral failure criterion. Without an explicit physical interpretation, this parameter cannot be calibrated directly through experimentation. In the simulation model, once the first principal stress within an element reaches a predefined threshold, it is eliminated automatically. This mechanism, in concert with the JC failure parameters, effectively dictates the failure dynamics of the material.

Overall, the penetration simulation model required 11 parameters, as summarized in Table 3.

Table 3

Parameters of the numerical simulation model for penetration tests

Category	Required parameters
JC constitutive equation	A	B	n	C	m
JC failure equation	D 1	D 2	D 3	D 4	D 5
Other supporting parameters	Failure criterion PC value

JC: Johnson-Cook.

The efficient determination of the parameters listed in Table 3 is described below.

(1) The JC constitutive equation is^[23]:

(1) $$ \sigma =(A+B\varepsilon^n)(1+C\mathrm{ln}\dot{\varepsilon}^n)(1-(T^*)^m), $$

where ε denotes the equivalent plastic strain, $$ \dot{\varepsilon}^*=\dot{\varepsilon}/\dot{\varepsilon}_0 $$ is the dimensionless plastic strain rate, and $$ \dot{\varepsilon} $$ represents the strain rate. The reference strain rate is set as $$ \dot{\varepsilon}_0 $$ = 0.001/s in this study. T^* = (T - T_r)/(T_m - T_r) signifies the dimensionless temperature in Kelvin (K), T is the absolute temperature, T_r is the room temperature at 298 K, and T_m is the melting point of the ultra-high-strength steel. The three terms in the equation correspond to strain hardening, strain-rate hardening, and thermal softening, respectively. Parameters A, B, n, C, and m are the material constitutive parameters that require fitting. These five parameters were determined using a clustering global optimization algorithm previously developed by our research group^[24].

(2) Understanding the behavior of a material under complex stress states is fundamental to the design of new structural materials. Macroscopic failure results from damage accumulation until a critical threshold is reached. This damage is associated with the microstructural changes induced by temperature, loading, and other service conditions, which initiate microdefects, cause them to grow and coalesce, and eventually form fractures. Because ductile fracture in metallic materials is strongly influenced by the stress triaxiality, strain rate, and temperature, an appropriate failure equation is required to describe the damage behavior of ultra-high-strength steel during loading.

The form of the JC failure equation is^[23]:

(2) $$ \varepsilon _f=[D_1+D_2\mathrm{exp}(D_3\sigma ^*)(1+D_4ln\dot{\varepsilon}^*)(1+D_5T^*) ] $$

In the equation, ε_f represents the failure strain, σ^* denotes the stress triaxiality, $$ \dot{\varepsilon}^*=\dot{\varepsilon}/\dot{\varepsilon}_0 $$ is the dimensionless plastic strain rate, and $$ \dot{\varepsilon} $$ represents the experimental strain rate. In this study, the reference strain rate is set at $$ \dot{\varepsilon}_0 $$ = 0.001/s. T^* = (T - T_f)/(T_m - T_f) indicates the dimensionless temperature in Kelvin (K), T is the temperature, T_f is the reference temperature, taken as room temperature (298 K) in this study, and T_m is the melting point of the ultra-high-strength steel. Parameters D₁-D₃ were derived by fitting the diameters of the fracture surfaces of the tensile test specimens with varying notch radii. Parameters D₄ and D₅ were derived by fitting the fracture surface diameters from the dynamic tensile specimens across different strain rates and from the high-temperature quasi-static specimens across various temperatures.

Traditional approaches for determining JC failure parameters require extensive empirical testing and often ignore parameter interdependencies. This disjointed fitting process is inefficient and rarely yields globally optimal solutions. By leveraging insights into the parameter influence patterns and utilizing global optimization algorithms, a novel methodology was developed to obtain the JC failure parameters from quasi-static and dynamic tensile tests. This approach enables accurate and relatively high-throughput parameter identification to build the service performance dataset.

The chaos optimization algorithm leverages the ergodicity and the sensitivity of chaotic systems to the initial conditions to enhance the search diversity. Unlike conventional algorithms, this approach effectively prevents premature convergence to local optima (a common complication in conventional optimization algorithms) making it highly suitable for complex multidimensional problems. Because it can search for multiple variables simultaneously, it is well-suited to the coupled multiparameter inversion problem addressed here.

Unlike traditional metaheuristic methods such as genetic algorithms (GA) and particle swarm optimization (PSO), which may converge prematurely in highly nonconvex multidimensional parameter spaces, the chaos optimization algorithm exploits the ergodicity and pseudo-randomness of chaotic maps. This feature renders it suitable for the inversion of JC parameters in the present study.

According to chaos theory, seemingly random behaviors can arise in deterministic systems. Although chaotic trajectories appear disordered, they follow deterministic rules and can provide broad, nonrepetitive coverage of a bounded interval after sufficient iteration. Chaotic variables can be generated using several mappings, including logistic, cubic, and infinite folding maps. Among these, the infinite folding map provides a relatively balanced distribution with fewer iterations, which helps accelerate the convergence in the present optimization procedure. Therefore, it was selected for the present study.

The logistic equation is shown in Equation (3), with the initial value μ in the range (0, 1). When ν = 4, the system is in a state of chaos^[25].

(3) $$ \mu _{i+1}=v\mu_i(1-\mu_i),\mu\in(0,1),\mu \ne 0.25,\ 0.5,\ 0.75 $$

The cubic mapping equation is given in Equation (4). Within the interval (-1, 1), chaos occurs when the initial value is not zero. After multiple iterations, the generated points were randomly distributed, and even a slight difference in the initial value could cause significant changes in the output^[26].

(4) $$ \mu _{i+1}=4(\mu_i)^3-3\mu_i,\mu\in(-1,1),\mu \ne 0 $$

The infinite folding mapping equation is presented in Equation (5). Within the interval (-1, 1), chaos occurs when the initial value is not zero. After multiple iterations, the generated points are randomly distributed, thereby facilitating traversal of the entire system space^[27].

(5) $$ \mu _{i+1}=\mathrm{sin}(\frac{2}{\mu_i}), \mu\in(-1,1),\mu \ne 0 $$

Chaos optimization linearly maps chaotic variables onto the search space of decision variables, using the resulting trajectories for global optimization. The algorithm explores candidate solutions through a chaotic search and progressively narrows the search interval to accelerate convergence. Because chaotic motion combines randomness with ergodicity, this procedure improves the chance of identifying a near-global or global optimum.

The simulations were performed using the actual dimensions of the experimental setup. The details are presented in Supplementary Figure 3.

In this study, specimen #121 of the G33 ultra-high-strength steel was examined with designated ranges for the JC failure parameters D₁-D₅. Numerical simulations of the quasi-static and dynamic tensile tests were performed to ascertain the impact of these parameters on the simulation results. Detailed analyses of the fracture elongation rates for both testing modalities were conducted, and the results are shown in Figure 4. The results indicated that increments in the JC parameters D₁ and D₂, along with a decrease in D₃, significantly affected the fracture elongation rates observed in the numerical simulations, leading to continuous linear variations in the outcomes. This reveals the critical influence of these parameters on the material performance. In contrast, the modulation of parameters D₄ and D₅ has minimal impact on the experimental results, highlighting their limited effect on the material properties.

Figure 4. Influence of JC failure parameters D₁-D₅ on numerical simulation results.

The fracture failure behavior of G33 ultra-high-strength steel under quasi-static tensile testing was systematically analyzed to determine the optimal ranges for the JC failure parameters D₁-D₃, as well as specific values for D₄ and D₅. The fracture elongation rates were thoroughly examined, and the results are shown in Figure 5. The dataset was divided into five intervals, and one representative material was selected from each interval for fitting, resulting in five samples. Materials whose fracture elongation was close to the median value of each interval were selected, while differences in microstructure and heat-treatment history were also considered to ensure representativeness. The fitted JC failure parameters are listed in Table 4.

Figure 5. Statistical chart of fracture elongation in quasi-static tensile tests of ultra-high-strength steel.

Table 4

Statistical table of JC failure equation parameters for 5 types of ultra-high-strength steel

ID	D 1	D 2	D 3	D 4	D 5
#25	0.0422	0.0736	-2.1799	-0.0185	0.5472
#80	0.0357	2.4998	-4.3184	-0.0068	-0.7419
#93	0.2073	1.6284	-2.8741	-0.0182	0.1258
#109	0.1139	3.1154	-5.6164	-0.0158	0.6378
#121	0.1569	2.6815	-5.4423	0.0062	-1.3978

JC: Johnson-Cook.

Within each interval, the fitted values for parameters D₁-D₃ were expanded both upward and downward to create numerical ranges that served as input domains for the chaotic optimization algorithm. Parameters D₄ and D₅ were assigned constant values common to all ultra-high-strength steel samples within these ranges. To precisely optimize parameters D₁-D₃ , simulation models were constructed for both the split Hopkinson tension bar (SHTB) test and the quasi-static tensile test. An infinite folding chaotic map is then used to generate candidate parameter sets within the prescribed ranges for global optimization. The optimization procedure is illustrated in Figure 6 and described below.

Figure 6. Chaotic optimization path for JC failure equation parameters. JC: Johnson-Cook; SHTB: split Hopkinson tension bar.

The model optimization procedure is summarized as follows:

① Initialization of Chaotic Variables and Setting Optimization Goals: Experimental results are used as target values for simulation calculations. Three distinct nonzero numbers, μ₁, μ₂, and μ₃, within the range of (-1, 1) are selected as initial values for the chaotic variables.

② Determination of Optimization Intervals for Parameters D₁-D₃: Incorporate the optimization interval (-1, 1) and adjust its size according to practical considerations to enhance the computational efficiency.

③ Mapping Chaotic Variables to Parameter Values: Using Equation (6), the chaotic variables are mapped into the optimization intervals to generate the parameters D₁-D₃, which are then input into the simulation model for computation^[28].

(6) $$ x=\frac{b+a}{2}+\frac{b-a}{2}\times \mu $$

④ Extraction of Simulation Results and Calculation of Simulation Errors^[29]: The error in the quasi-static tensile simulation is calculated as shown in Equation (7), where A_ZJTLS represents the fracture elongation rate from the quasi-static tensile simulation, and $$ true_{A_{ZJTLS}} $$ is the fracture elongation rate obtained experimentally. The error in the SHTB simulation is outlined in Equation (8), where A_SHTB denotes the fracture elongation rate from the SHTB simulation, and $$ true_{A_{SHTB}} $$ represents the experimentally obtained SHTB fracture elongation rate.

(7) $$ Error_{ZJTLS}=\left | \frac{A_{ZJTLS}-true_{A_{ZJTLS}}}{true_{A_{ZJTLS}}} \right | $$

(8) $$ Error_{SHTB}=\left | \frac{A_{SHTB}-true_{A_{SHTB}}}{true_{A_{SHTB}}} \right | $$

⑤ Calculation of Errors for Parameters D₁-D₃^[29]: Given the significant impact of parameters D₁-D₃ on both models, their error is computed as the sum of errors from the two simulation models, as outlined in Equation (9).

(9) $$ Error_{D_{1-3}}=Error_{ZJTLS}+Error_{SHTB} $$

⑥ Establishment of the Initial Optimal Solution Set: The number of solutions in the initial optimal set is defined as N, with the first N points generated by the chaotic variables serving as the initial optimal solutions. Let D_1,i denote the values of D₁ after the ith iteration. If the error for D_1,i ($$ Error_{D_1} $$) is less than the maximum error in the optimal set, f(D_1,j), where D_1,j corresponds to the point with the highest error, then D_1,j is replaced by D_1,i in the optimal set. Similarly, the optimal solution sets for D₂ and D₃ are updated according to the same principle.

⑦ Iteration and Convergence Criteria: A parameter K was specified. If the optimal solution set for a given parameter remains unchanged for K consecutive iterations, the bounds of that set are adopted as the new optimization interval for the next cycle, and the procedure returns to Step 2.

⑧ Termination of Optimization: The optimization process concludes when the errors for all parameters simultaneously meet the predefined target errors, that is, when $$ Error_{D_1} $$, $$ Error_{D_2} $$, and $$ Error_{D_3} $$ all satisfy the error criteria, and the optimized results are output. If the criteria were still not satisfied after several iterations, the parameter set with the smallest error in the optimal solution set was considered as the final result.

By combining experimental testing, numerical simulations, and parameter optimization, the JC constitutive and failure parameters were efficiently fitted for all 134 G33 steel variants.

(3) Validation of numerical simulation parameter accuracy in penetration testing. To verify the accuracy of the service performance data derived from the simulation, one ultra-high-strength steel sample was randomly selected from each of the three groups: samples lying on the fitted surface and the outliers above and below it. The samples were subjected to scaled projectile penetration simulations. As shown in Table 5, all the simulation errors were below 10%, thereby satisfying the dataset accuracy requirement. Supplementary Figure 4 compares the recovered projectile morphology after the penetration test at 727 m/s with the corresponding simulation results for G33 sample #128. The projectile shapes were consistent, and the residual velocity error after target penetration was only 3.6%. These results support the use of high-throughput simulations to generate service performance data for subsequent machine learning analysis.

Table 5

Comparison of target test results and numerical simulation results for the critical projectile fracture penetration speed of G33 ultra-high-strength steel

ID	Heat treatment process	Target test results (m/s)	Numerical simulation results (m/s)	Error rate (%)
#40	1,100 °C/1 h + OQ + -73 °C/1 h + 560 °C/5 h + AC	728	680	6.6
#83	1,050 °C/1 h + OQ + -73 °C/1 h + 200 °C/4 h + AC	830	915	9.3
#128	1,000 °C/1 h + WQ + 250 °C/5 h + AC	735	680	7.5

OQ: Oil quenching; AC: air cooling; WQ: water quenching.

This study established an efficient method for determining the numerical simulation parameters required for penetration analysis. Through a combined experimental validation and high-fidelity simulation, a paired dataset containing 483 sets of mechanical and service performance data was constructed (95 experimental service performance samples and 388 additional high-fidelity simulation samples). The dataset includes five target metrics: scaled projectile penetration velocity, penetration occurrence, post-penetration projectile integrity, residual velocity after perforation, and critical velocity for projectile fragmentation. Our strategy substantially reduced the experimental cost while improving the data acquisition efficiency, thereby providing a solid basis for the predictive modeling of the relationship between the mechanical properties and service performance.

Overview of the methodology

Based on the paired dataset generated in Section “Materials and Methods”, a two-stage machine-learning framework was constructed to map conventional mechanical properties to ballistic service responses. The overall workflow of this framework is illustrated in Figure 7. The framework was designed according to the physical sequence of the penetration process. In the first stage, a random forest classifier determined whether penetration occurred. For penetrated cases, the samples were passed to the second stage, in which a MLP regression was used to predict the residual velocity and critical fragmentation velocity, while a separate random forest classifier evaluated post-penetration projectile integrity. This staged design prevented the regression model from assigning physically unrealistic residual velocities to non-penetrated cases and improved the consistency between the model outputs and the actual penetration process.

Figure 7. Flowchart of the two-stage model development for predicting operational performance parameters. MLP: Multilayer perceptron.

Operational performance classification model based on random forest

Random forest classification models were developed to predict two service outcomes: whether penetration occurred and whether the projectile remained intact after penetration. Random forest is an ensemble learning method based on decision trees and is well known for its robustness and predictive accuracy on tabular data. It also provides estimates regarding the importance of features that are valuable for physical interpretations in material engineering. A decision tree consists of root, internal, and leaf nodes. The root node contains all the training data, each internal node splits the data according to a selected feature, and the leaf nodes represent the final categories.

The hyperparameters of the random forest models were optimized using GridSearchCV, which exhaustively searches for a predefined parameter grid and evaluates each combination through cross-validation. The parameter grid used in this study (“param_grid”) included the major random forest hyperparameters and their candidate values.

The specific grid parameters of the developed random forest classification model are listed in Supplementary Table 4.

Operational performance regression model based on MLP

The MLP is a feedforward artificial neural network composed of an input layer, one or more hidden layers, and an output layer. As each layer is fully connected to the next layer, the MLP can learn nonlinear relationships and is suitable for complex regression tasks. In this study, a relatively shallow MLP was used to capture the nonlinear mapping between the mechanical properties and ballistic performance while limiting the risk of overfitting.

To predict the continuous service performance metrics, namely the post-penetration residual velocity and critical fragmentation velocity, we developed a three-layer MLP regression model. The input layer contained 128 neurons, the hidden layer contained 64 neurons, and the output layer contained a single neuron that produced a predicted value. The detailed model configurations are listed in Supplementary Table 5.

The model was trained with the Adam optimizer using an initial learning rate of 10^-3 and a decay of 10^-5. This setting enabled stable convergence during training. All calculations were performed using Python 3.

Comparison with baseline models and ablation analysis

The benefits of the proposed architecture were assessed by comparing its performance with that of a conventional single-stage baseline model (for example, unified MLP). As depicted in Supplementary Figures 5 and 6, the single-stage model proved inherently deficient in that it was unable to map the abrupt discontinuity at the ballistic limit, and often erroneously predicted positive residual velocities for unpenetrated targets. In contrast, the two-stage RF + MLP framework imposed a physically meaningful classification boundary before regression, eliminated these zero-boundary anomalies, and reduced the mean absolute error (MAE) from 33.7 to 30.4 m/s.

Ablation analysis further demonstrated that each module contributed to the final performance. Removing the feature selection step introduced redundant dimensionality and weakened generalization, whereas omitting the chaos-based optimization compromised the quality of the simulation-generated data because the parameter search became more prone to converging to local optima. These results indicate that an integrated framework is necessary for accurate and physically consistent ballistic prediction.

Impact of data augmentation and generalization analysis

A key limitation of data-driven terminal ballistics research is the scarcity of experimental service performance data. In this study, the reliance on only 95 physical firing tests creates a serious risk of overfitting, in which the model learns the training samples but fails to generalize them to new material states.

The role of chaos-optimized finite element data augmentation was therefore evaluated by constructing a learning curve [Supplementary Figure 7] to track the training and cross-validation performances as the dataset became larger. The addition of increasing numbers of high-fidelity simulation samples steadily increased the cross-validation accuracy. The curve indicated that the model generalization stabilized only after the dataset exceeded 125 samples.

The evaluation with the learning curve indicates that the high predictive accuracy of the proposed framework is not simply a consequence of fitting a small dataset. Rather, it arises from the combination of physics-informed high-throughput simulation augmentation, which expands the dataset to 483 reliable samples, and a two-stage learning architecture that simplifies the prediction task and improves generalization.

Error analysis, validation, and model uncertainty

Stratified 10-fold cross-validation enabled the generalization to be rigorously assessed. Stratification preserved the class imbalance between penetration and non-penetration cases across all folds, thus reducing the evaluation bias toward the majority class.

The predictive reliability was further examined using a classification confusion matrix [Supplementary Figure 8] and a regression residual plot [Supplementary Figure 9]. The confusion matrix showed that the Stage-1 random forest accurately identifies the ballistic limit boundary. The residual plot for the Stage-2 MLP revealed that the residual velocity errors were approximately symmetric about zero without any obvious systematic bias.

Despite the high overall accuracy, analysis of the worst predicted cases provided important insights into the model uncertainty. The highest prediction variances were clustered exactly at the threshold of the ballistic limit (V₅₀). From a materials science perspective, the terminal ballistic response at this critical boundary is highly stochastic. The transition between partial embedment and complete perforation is highly sensitive to random microstructural inhomogeneities, such as the localized distribution of microvoids or inclusions. These microscopic defects introduced inherent aleatoric uncertainties into the dynamic fracturing process. Considering that the current model accepts deterministic macroscopic mechanical properties as its input, the prediction uncertainty increases near this transition point. Recognizing this stochastic boundary is important for the practical use of data-driven models in armor design.

Operational performance prediction based on the entire dataset

The input and output data were standardized, and a 10-fold cross-validation was used to divide the dataset into training and testing subsets. This was accomplished by randomly splitting the full dataset into ten subsets, with each iteration using one subset as the test set and the remaining nine subsets for training. This procedure yielded ten models, and the performance was evaluated based on the average accuracy across ten folds. The results reported in this section correspond to the full feature model before the feature selection.

The predictive accuracies of the various models on their test sets are detailed in Figure 8. The classification model for projectile penetration achieved a high accuracy of 98.8% with minimal variability (4.2%), thereby demonstrating robust generalization and stability across data partitions. In contrast, the classification model assessing the integrity of the projectile post-penetration achieved an accuracy of 72.9%, with a substantial variability of 18%, indicating potential dependency and sensitivity to specific datasets, which could constrain the predictive performance. The regressor for the residual velocity after penetration achieved an accuracy of 92.7%, with a variability of 5.6%, indicating good generalization and stability. Here, variability refers to the range of results obtained from 10-fold cross-validation.

Figure 8. Accuracy of the test set for the mechanical and operational performance prediction model (full-feature model): (A) Classification model for projectile penetration; (B) Classification model for post-penetration projectile integrity; (C) Regression prediction model for residual velocity after projectile penetration.

For the regression tasks, the accuracy was defined as (1 - MAPE) × 100%, where mean absolute percentage error (MAPE) is given by:

Here, n represents the total number of samples, y_i is the observed value of the ith sample, and $$ \hat{y}_i $$ is the predicted value of the ith sample.

Operational performance prediction based on parameter selection

The preceding analysis utilized a dataset containing 1,742 values relating to the mechanical properties, and predicted the service performance metrics with high accuracy. However, collecting a large number of mechanical parameters is time-consuming, expensive, and sensitive to experimental variations. In addition, many parameters were correlated with each other. To reduce redundancy and streamline data collection, the Pearson correlation coefficients were calculated between all feature pairs and used for feature selection.

The equation for calculating the Pearson correlation coefficient is presented in Equation (11)^[30]:

where x_i and y_i denote individual parameter values, whereas $$ \bar{x} $$ and $$ \bar{y} $$ are the means of the respective parameter sets, and S_x and S_y are their standard deviations. The sign of the Pearson correlation coefficient indicates whether the relationship is positive or negative and its magnitude indicates the strength of the correlation. Values closer to 1 or -1 indicate stronger correlations, with positive and negative signs representing direct and inverse relationships, respectively.

Figure 9 presents the Pearson correlation coefficient matrix for the 13 sets of parameters. In this research, a threshold of 0.9 was established. Any pair of parameters with an absolute correlation coefficient exceeding this value was considered highly redundant, which necessitated the retention of only one parameter from each correlated pair. Based on the analysis depicted in the figure, five parameter sets were removed: the quasi-static tensile strength, quasi-static compressive yield strength, Vickers hardness, dynamic compressive strength at 3,000 s^-1, and dynamic compressive strength at 5,000 s^-1.

Figure 9. Pearson correlation coefficient matrix.

After feature selection, eight parameters relating to the mechanical properties together with the penetration velocity were used as inputs to the two-stage model, and the service performance predictions were recalculated using 10-fold cross-validation. The results of the final optimized feature selection model are shown in Figure 10. The classification model that determined projectile penetration maintained a high accuracy of 98.3% with a variability of 4.2%, demonstrating robust predictive precision. This indicates that the elimination of the five redundant parameters had a minimal impact on the predictive performance of the model, whereas the nine retained parameters sufficiently encapsulated the critical determinants of the model output. The classification model assessing the integrity of the projectile post-penetration improved the output accuracy of the model to 74.8% with a reduced variability of 18%, thereby significantly lowering its sensitivity to specific data subsets. Meanwhile, the regression model predicting the residual velocity of the projectile after penetration enhanced the accuracy of the model output to 92.9% while maintaining a variability of 5.6%. These improvements further bolstered the generalization capabilities of the model.

Figure 10. Predictive accuracy of the test set (final optimized feature-selected model) for the mechanical and operational performance prediction model after feature optimization: (A) Classification model for projectile penetration; (B) Classification model for post-penetration projectile integrity; (C) Regression prediction model for residual velocity after projectile penetration.

By removing redundant parameters, the input dimensionality was reduced, training became more efficient, and model generalization was improved. This refinement allowed the model to capture the key factors governing the output with fewer parameters, thereby reducing the data collection effort and shortening the prediction process.

However, even after feature reduction, the predicted service performance variables remained relatively scattered, which made a direct comparison among materials less straightforward. Additional statistical processing is still required to efficiently rank the materials. Therefore, we selected the critical projectile fracture penetration velocity, which is defined as the highest penetration velocity at which the projectile remains intact, as a single comparative indicator that directly links the mechanical properties to service performance.

Figure 11 illustrates the predictive accuracy distribution of the regression model for assessing the critical projectile fracture penetration speed of G33 ultra-high-strength steel. The model achieved an average predictive accuracy of 89.8% with a variability of only 8.9%, demonstrating both high precision and robust stability in predicting critical speeds. The selected experimental data and the corresponding model predictions are presented in Table 6. This model enables the direct estimation of the critical projectile fracture penetration velocity from mechanical property data and thus provides a practical basis for rapid service performance assessment and material optimization.

Figure 11. Distribution of predictive accuracy for the regression model test set of critical projectile fracture penetration velocity (final optimized model) for G33 ultra-high-strength steel.

Feature interpretability and physical mechanisms

In materials informatics, model interpretability is essential for guiding material design. Therefore, SHapley Additive exPlanations (SHAP) analysis was performed on the nine key input features retained after feature selection to examine whether the proposed two-stage framework captured meaningful structure–property relationships. As shown in Supplementary Figure 10, the SHAP results indicate that the penetration response is mainly controlled by high-strain-rate mechanical properties, while quasi-static properties and impact velocity play secondary roles. Among all features, Dynamic Tensile Elongation at 2,000 s^-1 shows the highest importance, suggesting that dynamic ductility strongly affects plastic energy dissipation and resistance to back-face tearing or spalling under impact loading. Dynamic Compressive Yield Strengths at 5,000 s^-1 and 3,000 s^-1 are the next most influential features, reflecting their role in resisting localized plastic deformation and delaying adiabatic shear band formation during high-rate compression. In contrast, the lower importance of Impact Velocity and Quasi-Static Tensile Yield Strength further suggests that terminal ballistic performance is governed primarily by dynamic, rather than quasi-static, mechanical responses. These results confirm that the model predictions are physically consistent with established mechanisms of dynamic impact and fracture.