Design and Multi-Objective Optimization of Auxetic Sandwich Panels for Blastworthy Structures Using Machine Learning Method

Abstract

The design and multi-objective optimization of auxetic sandwich panels (ASPs) are performed to enhance the blastworthiness of armored fighting vehicles (AFVs). Various metastructures in the form of four auxetic geometries are proposed as the sandwich core: re-entrant honeycomb (REH), double-arrow honeycomb (DAH), star honeycomb (SH), and tetra-chiral honeycomb (CH). This paper employs a combination of finite element and machine learning methodologies to evaluate blastworthiness performance. Optimization is carried out using the nondominated sorting genetic algorithm II (NSGA-II) method. The optimization results show significant improvements in blastworthiness performance, with notable reductions in permanent displacement and enhancements in specific energy absorption (SEA). Global sensitivity analysis using SHapley Additive exPlanations (SHAP) reveals that cell thickness is the most critical factor affecting blastworthiness performance, followed by the number of cells and corner angle or radius for CH. The application of optimized ASP on AFVs shows promising results, with no failure occurring in the occupant floor. Furthermore, AFVs equipped with the optimized ASP DAH significantly reduce maximum displacement and acceleration by 39.00% and 43.56%, respectively, and enhance SEA by 48.30% compared to optimized aluminum foam sandwich panels. This study concludes that ASPs have potential applications in broader engineering fields.

1. Introduction

Landmines continue to pose significant threats to both military personnel and civilians. According to casualty statistics from military conflicts, approximately 60% of injuries are attributable to the explosion of mines and improvised explosive devices (IEDs) [1]. The International Campaign to Ban Landmines (ICBL) reported that in 2022, the total number of casualties and specifically anti-vehicle mines resulting from explosive remnants of war reached 4710 and 102, respectively [2]. Consequently, the development of blast-resistant or blastworthy structures for armored fighting vehicles (AFVs) as blast-load mitigation remains critical.

As blast load and IEDs become more powerful, the blastworthy structures of AFVs must be increasingly resilient to withstand such threats. Various methods for blast-loading mitigation have been researched through analytical, numerical, and experimental methods. As the simplest structure, single plates have been widely studied regarding their response to blast loading [3,4,5,6]. To enhance blastworthiness performance, modifications to single plates have been explored, such as adding stiffeners [7,8] and altering the flat plate to a V-shape [9,10]. However, single plates are insufficient for mitigating extreme blast loading. One of the most effective alternative blastworthy structures is the sandwich panel.

Sandwich panels are composite materials in which two or more different materials are combined at a macroscopic scale to produce a new material with superior performance compared to the constituent materials [11]. A sandwich panel consists of two thin, rigid face sheets bonded to a thick, lightweight core. This hybrid design fundamentally increases the moment of inertia of the structure without significantly adding mass. Various core structures have been studied for sandwich panels under blast loading. For example, Dharmasena et al. [12] investigated the blast response of metallic square honeycomb sandwich panel (HSP) structures. The results showed that the deformation of the back face of the sandwich panel was 40–90% smaller than that of an equivalent solid plate. Zhu et al. [13] conducted experimental studies on the blast responses of metallic hexagonal HSPs. The introduction of the hexagonal HSP significantly decreased the structure’s deflection. Zhang et al. [14] experimentally investigated the influence of geometric parameters of metallic trapezoidal corrugated-core sandwich panels subjected to blast loading. The results showed that increasing the sheet and cell thickness, as well as increasing the corrugated angle, enhanced the blast resistance. In addition to honeycomb and corrugated cores, aluminum foam has also proven effective as a core material in sandwich panels for withstanding blast loading. Numerous studies have examined aluminum foam sandwich panels (AFSPs) [15,16,17]; for instance, Hanssen et al. [18] examined the behavior of AFSPs under blast loading experimentally, analytically, and numerically. The results demonstrated that using foam as a sacrificial layer can control the contact stress level, providing local protection to the structure. Liu et al. [19] studied the responses of AFSPs under blast loading and found that the peak load of the sandwich panel was reduced by 61.54–64.69% compared to a single plate.

Conventional honeycomb structures typically exhibit positive Poisson’s ratio characteristics. Over the past three decades, there has been increasing interest in negative Poisson’s ratio (NPR) or auxetic materials due to their unique and contrasting properties compared to conventional materials. These auxetic materials are a prime example of metastructures, which are artificially designed structures with specific geometrical arrangements leading to unusual physical and mechanical properties. Auxetic structures, as a subset of metastructures, expand in all directions when stretched and contract when compressed. This opposite deformation behavior of auxetic materials leads to enhanced mechanical properties, such as more resistance to indentation [20,21], increased shear modulus [22,23], improvement in fracture toughness [24,25,26], and high energy absorption [20,27]. Auxetic materials also outperform conventional materials in other aspects, such as high energy harvesting efficiency for vortex-induced vibration problems [28] and improved acoustic performance [29]. Due to these superior properties, auxetic materials have found applications in diverse fields such as medicine, biomedicine, sports, and engineering (for more details, see reference [30]). Although auxetic materials have superior mechanical properties, the application of ASPs in AFVs remains relatively uncommon and under-researched.

Several studies have investigated auxetic sandwich panels (ASPs) under air-blast loading. For instance, Imbalzano et al. [31] and Yan et al. [32] studied the responses of ASPs and hexagonal HSPs to blast loading numerically and experimentally, respectively. Both studies demonstrated that ASPs exhibit higher blast resistance than conventional HSPs due to increased compressive stiffness as the impact zone of the ASP shrinks inward. Lan et al. [33] examined the dynamic response of cylindrical ASPs, conventional HSPs, and AFSPs subjected to blast loading. The results showed that cylindrical ASPs outperformed both conventional HSPs and AFSPs across all design parameter combinations. Qi et al. [34] conducted experiments on close-in blast loading on ASPs. The results showed that the introduction of the ASP effectively prevented damage to the concrete base under close-in blast loading. Chen et al. [35,36] conducted experimental studies on the blast response of ASPs with re-entrant and double-arrow structures as the core under paper tube-guided air-blast loading. The results showed that the double-arrow structure had 10.9–20.5% less permanent displacement than re-entrant structures. Yan et al. [37] conducted blast experiments on 3D-printed auxetic honeycomb sandwich beams (AHSBs). The results showed that AHSBs possess better blast performance than regular honeycomb. Overall, most research on ASPs under blast loading indicates that the ASP offers superior blastworthiness compared to previous structures.

The primary objective of blastworthiness performance is to protect occupants within a structure. One quantifiable metric is the displacement or acceleration experienced by occupants. Additionally, blastworthy structures must achieve sufficient lightweight construction to facilitate the mobility of AFVs. Given the complex and nonlinear relationships between the geometric parameters of auxetic structures and their blastworthiness performance, machine learning (ML) methods have emerged as valuable tools for modeling these intricate interactions. By utilizing surrogate models or metamodels, researchers can better predict outcomes based on various input variables. Furthermore, ML has been applied to predict specific properties of metamaterials, particularly auxetic structures, such as stiffness [38,39], sound absorption [40], Poisson’s ratio [39,41,42,43], and energy absorption capacity [38,44]. Accurately modeling these properties enhances the design and performance of the blastworthy structures.

As previously mentioned, the objectives of minimizing structural displacement often conflicts with the effort of weight reduction. Multi-objective optimization has been demonstrated as effective in addressing this challenge. Numerous studies have focused on multi-objective optimization of blastworthiness performance; for instance, Cong et al. [10] optimized the design of multi-V hulls for light AFVs using artificial neural networks (ANNs) and the nondominated sorting genetic algorithm II (NSGA-II). They optimized for lower tibia force and hull mass, achieving a 33–42% improvement in blastworthiness and a 43% reduction in mass compared to the baseline design. Qi et al. [45] and Wang et al. [46] conducted optimizations of curved and graded AFSPs, respectively, using NSGA-II with ANN and Kriging models as metamodels. Qi and Jiang et al. [47,48] optimized re-entrant ASPs using the radial basis function and NSGA-II. The optimum design produces better performance by reducing maximum deflection by 33% and increasing the areal specific energy absorption by 158.6%. Wang et al. [49] and Lan et al. [50] further optimized 3D double-arrow sandwich panels.

Previous studies have primarily focused on re-entrant and double-arrow structures as cores for ASPs. However, other fundamental auxetic structures, such as star honeycomb and chiral honeycomb, remain largely unexplored. This study investigates ASPs incorporating four basic auxetic geometries: re-entrant honeycomb (REH), double-arrow honeycomb (DAH), star honeycomb (SH), and chiral honeycomb (CH), specifically for air-blast-loading scenarios. Numerical simulations were employed to evaluate the blastworthiness performance of these ASP configurations. Subsequently, a metamodel was developed to establish relationships between the design variables of ASPs and their blastworthiness performance. Global sensitivity analysis was conducted to highlight the influence of each design variable on blastworthiness. Following this, multi-objective optimization was performed to generate a Pareto-optimal front. Finally, the optimal design from the Pareto front was identified as a protective structure suitable for AFVs subjected to air-blast loading.

2. Materials and Methods

2.1. Geometry Description

The configuration of the auxetic sandwich panel (ASP) is shown in Figure 1. The ASP consists of two face sheets, with a thickness of $t_f$, and an auxetic core, which has the total dimensions of $L\times W\times H$. The L, W, and H are the length, width, and height of the auxetic core in the x, y, and z directions, respectively. The auxetic core consists of repetitive 2D unit cells of auxetic geometry with the number of cells being $N_y$ and $N_z$ in the y and z directions, respectively. The auxetic core is oriented with its cells’ longitudinal direction parallel to the two face sheets. As research by Walkowiak et al. [51] shows, this configuration yields better blastworthiness than when the cells’ longitudinal direction is perpendicular to the face sheets. The size of the ASP is determined based on the floor size of the armored fighting vehicle (AFV) and the availability of additional space as a protective structure. Based on reference [52], the typical floor size of a small AFV is L = 1400 mm, W = 2500 mm, and H = 120 mm. The stand-off distance (SoD) from the blast charge position to the front face is 800 mm. However, to reduce the computational time, the half scale of ASP model was used, where L = 700 mm, W = 1250 mm, and H = 100 m with an SoD of 400 mm to compensate for the smaller size of ASP.

In this study, four auxetic basic geometries were proposed as the core of the ASP, which are re-entrant honeycomb (REH), double-arrow honeycomb (DAH), star honeycomb (SH), and tetra-chiral honeycomb (CH). Figure 2 shows the unit cell configuration of each auxetic geometry. There are four independent geometric parameters which are defined in this study: the cell width w, the cell height h, the corner cell angle $\theta$ and node radius r for CH, and the cell thickness t. The variables w and h are proportional with $N_y$ and $N_z$, respectively, except for the variable h of DAH, which depends on $N_z$ and $\theta$. Other dependent geometric parameters include inner cell spacing d ($d_1$ and $d_2$ for SH), strut length $l_1$ and $l_2$, and inner corner angle ${\theta }^{\prime }$ for DAH.

Table 1. Definition of geometric parameters for each auxetic geometry.

Geometry	Independent Variable	Dependent Variable	Constraint
		$l_1=h/\left(2cos\theta \right)$
REH	$N_y\left(w\right)$, $N_z\left(h\right)$, $\theta$, t	$l_2=\left(w+htan\theta \right)/2$	$\theta <{tan}^{-1}\left(w/h\right)$
		$d=\left(w-htan\theta \right)/2$
		$h=\left[H+(N_z-1)(w/2)tan\theta \right]/N_z$
		$l_1=\left(w/2\right)sec\theta$
DAH	$N_y\left(w\right)$, $N_z$, $\theta$, t	$l_2=\sqrt{h^2+{\left(w/2\right)}^2}$	$\theta <{tan}^{-1}\left(2h/w\right)$
		$d=h-\left(w/2\right)tan\theta$
		${\theta }^{\prime }={tan}^{-1}\left(2h/w\right)-\theta$
SH	$N_y\left(w\right)$, $N_z\left(h\right)$, $\theta$, t	$l_1=w/\left(4cos\theta \right)+htan\theta /\left(4\left(cos\theta +1\right)\right)$	$\theta <\min \left({tan}^{-1}\left(l_1/l_2\right),{tan}^{-1}\left(l_2/l_1\right)\right)$
		$l_2=h/\left(4cos\theta \right)+wtan\theta /\left(4\left(cos\theta +1\right)\right)$
		$d_1=2\left(l_{1}cos\theta -l_{2}cos\theta \right)$
		$d_2=2\left(l_{2}cos\theta -l_{1}cos\theta \right)$
CH	$N_y\left(w\right)$, $N_z\left(h\right)$, r, t	$l_1=\sqrt{w^2-4r^2}/2$	$r<\min \left(w/2,h/2\right)$
		$l_2=\sqrt{h^2-4r^2}/2$

summarizes the geometric parameters and geometric constraints of each auxetic geometry. The relative density is defined as the ratio of the effective density of the structures to the density of the matrix or constitutive material. This parameter indicates the porosity of the cellular structure. Based on those geometric parameters, the relative density (${\rho }_r$) of each auxetic geometry can be expressed as follows:

$$\begin{array}{cc}\hfill & {\rho }_{r,REH}={\displaystyle \frac{t}{w(t+H)}}\left({\displaystyle \frac{2h}{cos\theta }}{N}_z+{\displaystyle \frac{w+htan\theta }{2}}(2N_z+1)\right)\hfill \\ \hfill & {\rho }_{r,DAH}={\displaystyle \frac{2t{N}_z}{wH}}\left({\displaystyle \frac{w}{2cos\theta }}+\sqrt{h^2+{\displaystyle \frac{w^2}{4}}}\right)\hfill \\ \hfill & {\rho }_{r,SH}={\displaystyle \frac{3t}{wh}}{\displaystyle \frac{(w+h)(1+cos\theta +sin\theta )}{2cos\theta (cos\theta +1)}}\hfill \\ \hfill & {\rho }_{r,CH}={\displaystyle \frac{t}{wh}}\left(2\pi r+\sqrt{w^2-4r^2}+\sqrt{h^2-4r^2}\right)\hfill \end{array}$$

(1)

2.2. Numerical Modeling

In this study, numerical analysis was employed for all aspects of the ASP design analysis, including data sampling and the application of the ASP to AFVs. All numerical simulations were conducted using the nonlinear explicit finite element (FE) software LS-DYNA R12.1.0.

2.2.1. Geometry, Boundary Conditions, and Contact Modeling

The ASP was modeled as a quarter-symmetric model, as shown in Figure 3, due to the symmetric nature of the simulation. This approach reduces computational time by up to four times compared to a full model. Therefore, the quarter model of the ASP had dimensions of L = 350 mm, W = 625 mm, and H = 100 mm. Symmetric boundary conditions in x and y directions were applied to the ASP side in the $y-z$ and $x-z$ planes, respectively, while the outer face sheets were fixed on the other two sides. The face sheets and auxetic core were modeled using a fully integrated shell with five and two integration points, respectively, to minimize the hourglass energy. For the data sampling process, the thickness of both face sheets was set to 5 mm. The mesh convergence analysis was performed to obtain the mesh size for the auxetic core, which converged at a mesh size of 3 mm, while 10 mm was used for the face sheets [52]. AUTOMATIC_SINGLE_SURFACE was applied to all parts of the ASP for self-contact, while AUTOMATIC_SURFACE_TO_SURFACE was applied between the auxetic core and both face sheets to prevent penetration. The upper and lower nodes of the auxetic core was assumed to be fully connected with both face sheets using TIED_NODES_TO_SURFACE. The Static and dynamic coefficients of friction for all contact were 0.3 and 0.2, respectively.

2.2.2. Material Modeling

Considering its high strain capability, good manufacturability, and low cost, the ASP was modeled using stainless steel 304 material (SS304). Due to strain rate dependency in the blast simulation, a modified Johnson–Cook (J-C) model was adapted for the simulation. According to this constitutive model, the effective plastic stress ${\sigma }_{eff}^p$ is given by the following equation:

$${\sigma }_{eff}^p=\left(A+B{\left({\epsilon }_{eff}^p\right)}^n\right){\left(1+{\displaystyle \frac{\dot{\epsilon }}{\dot{{\epsilon }_0}}}\right)}^C\left(1-{\left({\displaystyle \frac{T-T_r}{T_m-T_r}}\right)}^m\right)$$

(2)

The first term represents material plasticity, where ${\epsilon }_{eff}^p$ is the effective plastic strain, A is the yield strength, and B and n are strain-hardening parameters. The second term accounts for the strain rate dependency, where $\dot{\epsilon }$ is the strain rate, $\dot{{\epsilon }_0}$ is the reference strain rate, and C is the strain rate sensitivity coefficient. The final term addresses temperature dependency, where T, $T_r$, and $T_m$ are the working temperature, room temperature, and melting temperature, respectively, and m is the thermal softening coefficient.

All material parameters can be determined through curve fitting of experimental data. The modified J-C parameters for SS304 were obtained from curve fitting of the stress–strain curve by references [14,53,54], as shown in Figure 4. These material parameters are summarized in

Table 2. Modified Johnson–Cook parameters of 304 stainless steel.

E	$\mathit{\nu }$	$\mathit{\rho }$	A	B	n	$\dot{{\mathit{\epsilon }}_0}$	C	${\mathit{T}}_{\mathit{r}}$	${\mathit{T}}_{\mathit{m}}$	m	${\mathit{C}}_{\mathit{p}}$	$\mathit{\chi }$	$\mathit{\alpha }$
(GPa)	(-)	(kg/m3)	(MPa)	(MPa)	(-)	(s−1)	(-)	(K)	(K)	(-)	(J/kg·K)	(-)	(K−1)
200	0.3	7900	310	1872	0.96	0.001	0.016	293	1673	1	500	0.9	$1.5\times {10}^{-5}$

.

For the material failure criterion, the Cockcroft–Latham damage evolution rule $\left(D\right)$ was used to simulate the failure of the modified J-C model. This criterion is based on the concept that ductile fracture occurs when the accumulated plastic strain energy in a material reaches a critical value [55]. Mathematically, the Cockcroft–Latham damage criterion can be written as

$$D={\displaystyle \frac{D_C}{W_C}}\max ({\sigma }_1,0){\epsilon }_{eff}^p$$

(3)

where $D_C$ is the critical damage parameter, ${\sigma }_1$ is the maximum principal stress, and $W_C={\int }_0^{{\epsilon }_f}\max ({\sigma }_1,0)d{\epsilon }_{eff}^p$ is the plastic strain energy, with ${\epsilon }_f$ as the strain failure. The value of $D_C$ is 1 [52], while $W_C$ is obtained from the stress–strain curve, with a value of 296 MPa. In this study, failure criteria were not applied during the data sampling process but were only considered for the application of the AFV subsystem.

2.2.3. Blast Load Modeling

Several blast-load models can be implemented in numerical simulations. The simplest model for blast-load simulation is the Conventional Weapons Effects Program (CONWEP), developed by Kingery and Bulmash [56]. This model is selected for its computational efficiency, as it does not require Euler field or particle calculations. The CONWEP model considers two types of pressure from a blast load: reflected pressure and incident pressure, as given in the following equation:

$$P=P_r{cos}^2{\theta }_i+P_i\left(1+{cos}^2{\theta }_i-2cos{\theta }_i\right)$$

(4)

where ${\theta }_i$ is the incidence angle, $P_r$ is the reflected pressure, and $P_i$ is the incident pressure. This empirical model is valid in the range of 0.147 m/kg^1/3$<Z<$ 40 m/kg^1/3, where $Z=R/m^{1/3}$ is called the scaled distance, with R is the stand-off distance (SoD) and m is the mass of TNT. As previously mentioned, the R and m used in this study are 400 mm and 8 kg, respectively, such that Z equal to 0.2 m/kg^1/3, which is still inside the range.

Another method for simulating blast phenomena is smoothed-particle hydrodynamics (SPH), a particle-based computational technique developed by Gingold and Monaghan [57] and Lucy [58]. SPH is a mesh-free Lagrangian method where the system’s state is represented by a finite number of discrete particles whose positions can move according to the velocity field. In this study, the SPH model does not account for the effects of air particles and the casing of the charge. The material model HIGH_EXPLOSIVE_BURN was used for the explosion, characterized by parameters such as density ($\rho$), detonation velocity ($v_D$), and Chapman–Jouget pressure ($P_{CJ}$). The corresponding Jones–Wilkins–Lee (JWL) equation of state (EOS) was employed to compute pressure as a function of initial relative volume and internal energy per unit volume. The EOS of JWL is defined by

$$P=A\left(1-{\displaystyle \frac{\omega }{R_{1}V}}\right)e^{-R_{1}V}+B\left(1-{\displaystyle \frac{\omega }{R_{2}V}}\right)e^{-R_{2}V}+{\displaystyle \frac{\omega E}{V}}$$

(5)

where A, B, $R_1$, $R_2$, and $\omega$ are the constants, E is the internal energy per unit volume, and $V=\rho /{\rho }_0$ is the ratio of the density at the time of the explosion and the initial density of TNT. The material parameters and the corresponding JWL EOS of TNT are given in

Table 3. High explosive burn material and Jones–Wilkins–Lee equation of state parameter of TNT [59].

$\mathit{\rho }$	${\mathit{v}}_{\mathit{D}}$	${\mathit{P}}_{\mathbf{CJ}}$	A	B	${\mathit{R}}_1$	${\mathit{R}}_2$	$\mathit{\omega }$	${\mathit{E}}_0$	${\mathit{V}}_0$
(kg/m3)	(m/s)	(GPa)	(GPa)	(GPa)	(-)	(-)	(-)	(MJ/m3)	(-)
1630	6930	21	371.2	3.231	4.15	0.95	0.3	7000	1

.

2.3. Machine Learning Model

Machine learning (ML) methods are employed to develop a surrogate model capable of predicting outputs based on specific input variables. The predicted output $\widehat{\mathit{Y}}=\{{\widehat{y}}_1,{\widehat{y}}_2,\dots ,{\widehat{y}}_n\}$ is related to the input variables $\mathit{X}=\{x_1,x_2,\dots ,x_m\}$ via a black box function $\widehat{\mathit{f}}\left(\mathit{X}\right)$, where m and n are the dimensions of input and output variables, respectively. The ML model $\widehat{\mathit{f}}\left(\mathit{X}\right)$ serves multiple purposes. Firstly, it can be utilized to explore the design space, providing valuable insights into how each input variable interacts and contributes to the output. Secondly, it is used for optimization, enabling the determination of the Pareto front of the model for specific objective functions.

2.3.1. Metamodel

This study utilizes an artificial neural network (ANN) as the metamodel to predict outputs. A typical ANN works like a black box to predict outputs based on given inputs. The working concept of ANNs is inspired by the working process of a biological neuron system, where many neurons are interconnected to form a complex and nonlinear network. The adaptive nature of input changes allows the ANN method to produce the best output without altering the output criteria. A typical ANN architecture consists of three layers: (1) the input layer, (2) hidden layers, and (3) the output layer. Neurons at the input and output layers represent the number of input and output variables, respectively. Each pair of neurons is connected with the connection weight w, and each neuron possesses a unique bias b. Mathematically, the output of the l-th layer in an ANN algorithm can be expressed as follows:

$${\mathit{o}}^{\left(l\right)}=f_a\left({\mathit{W}}^{\left(l\right)}{\mathit{o}}^{(l–1)}+{\mathit{b}}^l\right),l=1,2,...,k$$

(6)

where k is the number of layers, $f_a$ is the activation function, $W^{\left(l\right)}$ is the weight matrix, and $b^{\left(l\right)}$ is the bias vector. ${\mathit{o}}^{\left(0\right)}$ and ${\mathit{o}}^{\left(k\right)}$ represent the input and output layers, respectively. During the model training process, the weight and bias values undergo meticulous updates through tailored training algorithms and optimizers, aimed at improving prediction accuracy and minimizing the loss incurred during prediction. For regression problems, common loss functions used to evaluate model training include the maximum absolute error (MAX) and mean squared error (MSE). Additionally, R-squared ($R^2$) is employed to quantify how effectively the input variables explain the variability in the output variable. These metrics, MAX, MSE, and $R^2$, are calculated as follows:

$$\begin{array}{cc}\hfill & \mathrm{MAX}=\max \left|y_i-{\widehat{y}}_i\right|\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & \mathrm{MSE}={\displaystyle \frac{1}{n_s}}\sum _{i=1}^{n_s}{\left(y_i-{\widehat{y}}_i\right)}^2\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & R^2=1-{\displaystyle \frac{{\sum }_{i=1}^{n_s}{\left(y_i-{\widehat{y}}_i\right)}^2}{{\sum }_{i=1}^{n_s}{\left(y_i-\overline{y}\right)}^2}}\hfill \end{array}$$

(9)

where $n_s$ is the number of samples, $y_i$ is the actual value, and $\overline{y}$ is the mean value.

2.3.2. Multi-Objective Optimization

A multi-objective optimization problem (MOOP) involves achieving optimal solutions from two or more conflicting objectives. Genetic algorithms (GAs) are evolutionary optimization techniques widely applied to such problems. The GA operates on the principles of natural selection inspired by Charles Darwin’s theory, employing genetic operators including selection, crossover, and mutation [60]. In this study, the second-generation nondominated sorting genetic algorithm or NSGA-II is used as the multi-objective optimization method. NSGA-II represents an advancement over NSGA proposed by Deb [61], recognized as a popular GA for generating the Pareto front in MOOP within engineering. NSGA-II employs two effective sorting principles: elitist nondominated sorting to rank solutions and crowding distance sorting to maintain diversity in the solution set.

2.3.3. Sensitivity Analysis

A sensitivity analysis assesses how variations in input values influence output variables, providing insights into the robustness and reliability of an ML model. SHapley Additive exPlanations (SHAP) [62] comprise an explainability technique based on the Shapley value concept introduced by Shapley [63] in cooperative game theory. Essentially, the Shapley value calculates the average contribution of a feature across all possible subsets F. Unlike variance-based decomposition methods such as Sobol indices, SHAP values offer a more accurate reflection of the input–output relationship in design exploration, as input variables are considered to be nonrandom [64]. The Shapley value of the j-th feature is defined as

$${\varphi }_j\left(\widehat{\mathit{f}},\mathit{X}\right)={\displaystyle \frac{1}{m}}\sum _{S\subseteq F\setminus \left\{j\right\}}{\left({\displaystyle \genfrac{}{}{0pt}{}{m-1}{\left|S\right|}}\right)}^{-1}\left({\widehat{\mathit{f}}}_{S\cup \left\{j\right\}}\left({\mathit{X}}_{S\cup \left\{j\right\}}\right)-{\widehat{\mathit{f}}}_S\left({\mathit{X}}_S\right)\right)$$

(10)

where S is the one subset of F and the right term is the marginal contribution of feature j to a subset S. An essential characteristic of the Shapley value is its efficiency, ensuring that the total sum of Shapley values for all features equals the model output. Therefore, the output value of $\widehat{\mathit{f}}\left(\mathit{X}\right)$ can be expressed as

$$\widehat{\mathit{f}}\left(\mathit{X}\right)={\varphi }_0+\sum _{j=1}^m{\varphi }_j\left(\mathit{X}\right)$$

(11)

where ${\varphi }_0$ is the output of reference input value or ${\widehat{\mathit{f}}}_{\varnothing }\left({\mathit{X}}_{\varnothing }\right)$. For a global sensitivity analysis, averaged SHAP values are utilized to assess feature sensitivity at a global scale [64]. Averaged SHAP values indicate the average impact of input variables on specific outputs. The averaged SHAP of the j-th feature can be expressed as

$$\left|\overline{{\varphi }_j}\right|={\displaystyle \frac{1}{n_s}}\sum _{i=1}^{n_s}\left|{\varphi }_j\left({\mathit{X}}^{\left(i\right)}\right)\right|$$

(12)

3. Results and Discussion

3.1. Numerical Validation

To validate the numerical models employed, the response to blast events was compared with experimental data. This study encompasses two main validation aspects: (1) validating the modeling of blast loads and (2) validating the blast response of sandwich panels.

3.1.1. Validation and Comparison of Air-Blast Model

Firstly, three air-blast models were validated and compared: (1) CONWEP, (2) SPH, and (3) multi-material Arbitrary Lagrangian–Eulerian (MMALE). The CONWEP and SPH models are detailed in Section 2.2.3. MMALE represents an advanced method in which air-blast dynamics are modeled using Eulerian formulations. This approach allows blast waves to propagate through the air domain and interact with Lagrangian structures. In this validation, the contact algorithm used between the Eulerian and Lagrangian parts is a penalty-based method through the keyword CONSTRAINED_LAGRANGE_IN_SOLID.

Experimental data from Børvik et al. [6] used as validation for the air-blast models. A single AL-6XN steel plate specimen with dimension of 406 × 406 × 3.4 mm³ was subjected to a 0.15 kg C-4 explosive charge at an SoD of 150, 200, and 250 mm. For the CONWEP model, 0.15 kg C-4 was equated to 0.18 kg TNT based on mass equivalence [65]. The J-C material model parameters from reference [6] were implemented in the numerical simulations. The SPH and MMALE methods utilized the MAT_HIGH_EXPLOSIVE_BURN model and JWL EOS for modeling the C-4 charge, with material properties and JWL EOS details provided in [59]. The primary comparative metric was the permanent displacement (${\delta }_p$) of the central node of the plate. Figure 5 illustrates the FE models of the three air-blast models.

Table 4. Permanent displacement from various air-blast models and experimental results [6].

SoD (mm)	Exp. (mm)	CONWEP		SPH		MMALE
		Sim. (mm)	Err. (%)	Sim. (mm)	Err. (%)	Sim. (mm)	Err. (%)
150	17	17.81	4.76	17.75	4.41	17.41	2.41
200	12.7	14.77	16.30	14.23	12.05	12.54	−1.26
250	11.3	11.23	−0.62	11.42	1.06	10.27	−9.12

summarizes the numerical results compared to the experimental results. MMALE shows the smallest error for a SoD of 150 and 200 mm, while CONWEP performs best for a SoD of 250 mm. On average, MMALE exhibits the smallest error, followed by SPH and CONWEP. Figure 6 shows the bar plot of ${\delta }_p$ for each air-blast model alongside experimental data. This observation is consistent with the fact that CONWEP relies solely on empirical formulas with inherent limitations. In contrast, SPH and MMALE provide more accurate and detailed simulations by solving the governing equations of fluid and solid dynamics. However, the computational time required by CONWEP is significantly faster, approximately 55 times faster than SPH and 305 times faster than MMALE. Considering computational efficiency, the CONWEP method is chosen for the data sampling process, while the SPH model is used for final evaluation.

3.1.2. Validation of Sandwich Panels

The experiments conducted by Børvik et al. [6] only involved a single plate, whereas this study focuses on the blast response of a sandwich panel. Therefore, it is essential to validate the numerical model of the sandwich panel. Experimental data provided by Zhang et al. [14] are relevant as they also investigated the blast response of a sandwich panel, despite using a nonauxetic core geometry. The experimental setup, objectives, and material properties align closely with those used in this paper. Zhang et al. [14] conducted experiments on metallic trapezoidal corrugated-core sandwich panels under air-blast loads. Figure 7 shows the experimental setup and configuration of a trapezoidal corrugated-core sandwich panel. Their study primarily explored the influence of geometric parameters and SoD on midpoint displacements. The horizontal plane had dimensions of 452 mm × 440 mm, and the exposed area of the blast or the corrugated-core dimensions were 300 mm × 288 mm.

The numerical model utilized a symmetric-quarter model with symmetric boundaries fixed at the outer sides. The material properties of 304 stainless steel were represented using a modified J-C model, with parameters listed in Table 2. A 55 g cylindrical TNT charge, with a radius of 17.5 mm and height of 37.2 mm, was simulated using the SPH method due to its specific shape, while the CONWEP method approximated the charge as a spherical air blast.

Figure 8 shows a correlation plot of ${\delta }_p$, comparing experimental and numerical simulation results. The close alignment of data points around the line of a perfect match indicates strong agreement between numerical predictions and experimental data. The average error across all samples, as detailed in

Table 5. Permanent displacement from numerical and experimental results [14] for trapezoidal corrugated-core sandwich panels.

Specimen	Front Sheet			Back Sheet
	Exp. (mm)	Sim. (mm)	Err. (%)	Exp. (mm)	Sim. (mm)	Err. (%)
TZ-2	28.81	27.97	−2.92	14.14	15.19	7.43
TZ-4	34.33	30.89	−10.02	19.81	18.37	−7.27
TZ-5	25.35	24.00	−5.33	11.47	11.39	−0.70
TZ-6	29.55	27.65	−6.43	16.17	16.46	1.79
TZ-7	26.77	26.86	0.34	13.33	14.02	5.18
TZ-8	31.07	29.49	−5.09	17.15	17.15	−0.06
TZ-9	26.94	26.09	−3.16	13.00	13.81	6.23
TZ-10	29.90	28.13	−5.92	16.49	17.64	6.97
TZ-11	27.37	25.71	−6.07	13.91	14.38	3.38
TZ-12	23.98	23.49	−2.04	16.23	17.65	8.75
TZ-13	35.72	34.53	−3.33	11.26	11.07	−1.69

, is less than 10%. To further validate the numerical model’s reliability, deformation modes of three sandwich panels with different core configurations were compared between numerical simulations and experimental results. As shown in Figure 9 (left), the numerical model accurately reproduces major deformations, including local indentation in the middle area of the corrugated core and plastic buckling of the core web. Figure 9 (right) also shows the cross-sectional displacement profiles, showing close correspondence between numerical and experimental results. Overall, the numerical methods employed in this study demonstrate validity and accuracy suitable for subsequent optimization processes.

3.2. Metamodel and Multi-Objective Optimization

In this study, MOOP was defined and solved to find the optimal solution from two conflicting objectives. As previously mentioned, the metamodel was constructed using an ANN method, while the multi-objective optimization was performed using NSGA-II. Before establishing the ANN algorithm, a design of experiments (DoE) was generated using numerical methods. To construct the metamodel, 160 samples for each auxetic geometry were generated using the Latin Hypercube Sampling (LHS) method and simulated using the explicit FE solver LS-DYNA R12.1.0.

3.2.1. Design Variables, Objective Function, and Constraint

Two conflicting objective variables were optimized in this MOOP, which are the permanent displacement (${\delta }_p$) of the ASP back face and specific energy absorption (SEA) of the auxetic core. The ${\delta }_p$ was obtained by averaging the displacement-time history in the area where the displacement had already converged, while the SEA of the core was calculated by dividing the internal energy E by the core mass m. The objective functions are to minimize the ${\delta }_p$ and −SEA.

The design optimization focused solely on the independent geometric parameters of each auxetic geometry. The two face sheets, material types, and other variables were not optimized since this study focuses only on the influence of auxetic geometry on blast responses. The independent variables w and h were converted into the number of cells $N_y$ and $N_z$, respectively, due to the constraint of fixed ASP dimensions. The design variables for the MOOP for REH, DAH, and SH are specified as follows:

$$\begin{array}{cc}\hfill & \min \left({\delta }_p,-\mathrm{SEA}\right)\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}5\le N_y\le 25\hfill \\ 2\le N_z\le 6\hfill \\ 0<\theta <{\theta }_{max}\hfill \\ 0.5\le t\le 2\hfill \end{array}\right.\hfill \end{array}$$

(13)

and for the CH, the variable $\theta$ is changed with the r:

$$\begin{array}{cc}\hfill & \min \left({\delta }_p,-\mathrm{SEA}\right)\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}5\le N_y\le 25\hfill \\ 2\le N_z\le 6\hfill \\ 2<r<r_{max}\hfill \\ 0.5\le t\le 2\hfill \end{array}\right.\hfill \end{array}$$

(14)

Design constraints for the MOOP are listed in

Table 6. Design constraints for each auxetic model in the multi-objective optimization problem.

Geometry	Design Constraint
REH	$d>4$
DAH	${\theta }^{\prime }>10$
SH	$d_1,d_2>4$
CH	$r<\min \left((w-4)/2,(h-4)/2\right)$

, taking into account the geometric limitations. Based on the design variables and constraints, the design space with the corresponding normalized relative density (${\rho }_r$) for each auxetic model is presented in Figure 10. As shown in Figure 10, higher input variables result in a higher ${\rho }_r$ of the model.

3.2.2. Metamodel Accuracy

In establishing the metamodel, 90% of the DoE samples were used as the train set, while 10% were reserved as the test set. Additionally, 10% of the training set was employed as the validation set during the training process. This ratio was chosen to ensure model accuracy given the limited number of sampling points generated in this study. The Adam optimizer was used due to its efficiency in accelerating the gradient descent algorithm. A learning rate of 0.001 and an early stopping algorithm were also employed to monitor the loss of the validation set, preventing overfitting. All hidden layers used the ReLU (rectified linear unit) activation function to capture model nonlinearity, while the output layer used a linear activation function to produce continuous values.

The error parameters of each metamodel are presented in

Table 7. Error parameters of all auxetic metamodels.

Geometry	${\mathit{\delta }}_{\mathit{p}}$			SEA
	MAX (mm)	MSE (mm2)	${\mathit{R}}^{\mathbf{2}}$ (-)	MAX (J/g)	MSE (J2/g2)	${\mathit{R}}^{\mathbf{2}}$ (-)
REH	3.059	0.713	0.991	0.418	0.009	0.996
DAH	3.307	0.684	0.994	0.302	0.006	0.997
SH	4.524	1.188	0.992	0.250	0.006	0.997
CH	4.309	1.779	0.982	0.498	0.012	0.996

. All MAX values are less than 5 mm for ${\delta }_p$ and 0.5 J/g for the SEA, respectively. Generally, the $R^2$ values of all metamodels are greater than 0.95, indicating that the predicted values from the ANN model are very close to the simulation values. The predicted values of the trained metamodel for each auxetic geometry are shown in Figure 11. The deviations between the predicted values and the simulation sampling points are very small. Based on the three evaluation parameters, all metamodels are confirmed to be accurate and can be used for the subsequent multi-objective optimization process.

3.2.3. Global Sensitivity Analysis

In this section, the sensitivities of each input variable are analyzed using the SHapley Additive exPlanations (SHAP) method on 50,000 generated samples. The average SHAP values for ${\delta }_p$ are presented in Figure 12a. It can be observed that the variables $N_y$ and t exhibit very high average SHAP values, indicating a significant influence on ${\delta }_p$. On average, varying $N_y$ and t results in changes in ${\delta }_p$ of 5.06 mm and 5.68 mm, respectively. In contrast, the variables $N_z$ and $\theta /r$ have relatively low average SHAP values in ${\delta }_p$ of 2.07 mm and 1.42 mm, respectively. The averaged SHAP of SEA is given in Figure 12b. Similar to ${\delta }_p$, the variable t has the highest average SHAP value, contributing an average change of 1.21 J/g, followed by $N_z$ with 0.89 J/g, $N_y$ with 0.77 J/g, and $\theta /r$ with 0.43 J/g. Overall, the variable t has the greatest influence on both blastworthiness performance metrics, while $\theta /r$ has the least influence.

To provide a more intuitive understanding of how each input variable affects the output, summary plots of SHAP for each auxetic geometry are presented in Figure 13. These summary plots display the SHAP values for each variable alongside the input values for each sample. The sign of SHAP values indicates the change in output value with respect to the reference value. In this case, the reference value is the baseline design. The baseline design configuration is determined by selecting the midpoint value of the design variable range, with the thickness set to 1 mm. In general, higher input values result in lower outputs for both ${\delta }_p$ and SEA, as indicated by the gradual change in color. However, the $\theta$ variable shows a nonlinear influence on ${\delta }_p$, as indicated by disordered coloring. Figure 14 shows the scatter plots of the influence of each feature for all auxetic geometries. As illustrated in Figure 14, high input values, especially for variables $N_y$ and t, correspond to the lower-left region of the plot, indicating small ${\delta }_p$ and small SEA, and vice versa. The SHAP values from Figure 13 are consistent with the trends observed in Figure 14, as both depict the same influence of inputs on outputs. This consistency between SHAP values and scatter plots highlights the robustness of the input–output relationship analysis.

3.2.4. Pareto Front and Optimal Solution

Based on NSGA-II optimization, the Pareto front for each auxetic model is shown in Figure 15. Each orange circle in Figure 15 represents a nondominated solution. As expected, the Pareto front illustrates the conflict between the two objectives, positioned at the outer side of the prediction sample marked with light blue circles. As ${\delta }_p$ increases, SEA also increases, and vice versa.

Several design points from Figure 15, including the baseline, ideal, balanced, and optimized designs, are listed in

Table 8. Design variables and blastworthiness performance of selected design point for all auxetic geometries.

Geometry	Type	${\mathit{N}}_{\mathit{y}}$ (-)	${\mathit{N}}_{\mathit{z}}$ (-)	$\mathit{\theta }/\mathit{r}$ (°/mm)	t (mm)	${\mathit{\rho }}_{\mathit{r}}$ (-)	${\mathit{\delta }}_{\mathit{p}}$			SEA
							Pred. (mm)	Sim. (mm)	Err. (%)	Pred. (J/g)	Sim. (J/g)	Err. (%)
REH	Baseline	15	4	30.0	1.00	0.115	27.50	27.64	−0.51	2.13	2.06	3.40
	Ideal min ${\delta }_p$	25	6	45.3	1.99	0.434	10.07	8.32	21.03	0.20	0.22	−9.09
	Ideal max SEA	5	2	0.1	0.54	0.022	49.96	50.89	−1.83	8.52	8.58	−0.70
	Balanced	19	2	26.3	0.53	0.059	31.83	33.99	−6.35	3.98	4.08	−2.45
	Optimized	9	3	61.5	1.63	0.203	18.58	17.92	3.68	1.40	1.19	17.65
DAH	Baseline	15	4	30.0	1.00	0.123	30.39	33.12	−8.24	2.25	2.32	−3.02
	Ideal min ${\delta }_p$	23	6	0.1	2.00	0.310	3.16	2.93	7.85	0.64	0.63	1.59
	Ideal max SEA	6	2	2.6	0.5	0.024	48.88	49.51	−1.27	8.58	8.56	0.23
	Balanced	25	2	0.1	0.68	0.070	31.56	33.90	−6.90	3.66	3.43	6.71
	Optimized	15	5	0.1	1.72	0.205	10.24	9.76	4.92	1.29	1.31	−1.53
SH	Baseline	15	4	17.0	1.00	0.115	34.85	35.41	−1.58	2.44	2.39	2.09
	Ideal min ${\delta }_p$	20	6	16.1	2.00	0.328	11.33	11.78	−3.82	0.72	0.65	10.77
	Ideal max SEA	13	2	1.1	0.51	0.032	49.89	50.25	−0.72	7.42	7.59	−2.24
	Balanced	25	2	0.6	0.89	0.081	32.76	33.77	−2.99	3.69	3.57	3.36
	Optimized	15	5	0.3	1.86	0.207	14.94	14.42	3.61	1.25	1.29	−3.10
CH	Baseline	15	4	5.0	1.00	0.091	38.51	38.95	−1.13	3.17	3.12	1.60
	Ideal min ${\delta }_p$	25	6	6.3	2.00	0.346	14.64	11.81	23.96	0.42	0.53	−20.75
	Ideal max SEA	17	2	2.1	0.62	0.034	49.92	49.91	0.02	6.93	6.92	0.14
	Balanced	18	3	14.5	0.66	0.072	32.57	34.33	−5.13	3.61	3.53	2.27
	Optimized	18	3	14.6	1.88	0.206	18.87	15.44	22.22	1.62	1.28	26.56

. The table shows that the average prediction error from the ML model compared to the simulation is relatively small, with values less than 10%, except for the ideal design with minimum ${\delta }_p$ and the optimized designs for the REH and CH models. However, when examining the absolute value difference between prediction and simulation, it remains below the MAX value from our trained model, as indicated in Table 7. Therefore, considering the small sample size used in this study, we can conclude that the model’s error is within an acceptable range, indicating that the optimization results are reliable.

3.3. Discussions

3.3.1. Comparative Analysis

A comparison between the Pareto front and selected design points was conducted to evaluate improvements achieved through the NSGA-II optimization. First, the baseline design, marked with a red triangle, is located at the lower right of the Pareto front, as shown in Figure 15. The optimization results in an improvement of 12.55–35.60% in ${\delta }_p$ reduction for the same SEA and 26.21–83.68% in SEA enhancement for the same ${\delta }_p$. Therefore, the Pareto front offers superior performance in at least one objective.

Next, the balanced design, marked by the purple pentagon in Figure 15, is compared with the baseline design. A balanced design is defined as a point on the Pareto front that represents the balanced trade-off among conflicting objectives. This point is identified as the one with the minimum distance to the utopia point. The utopia point, an idealized design point represented by dark blue squares in Figure 15, corresponds to the optimal objective values (in this study, the optimal objective values are the lowest ${\delta }_p$ and highest SEA). Therefore, the utopia point is linked to the two ideal design points, marked by green “+” symbols. As given in Table 8, although the balanced designs for REH and DAH have larger ${\delta }_p$ than the baseline design by 22.97% and 2.36%, respectively, they offer superior blastworthiness in SEA, with improvements of 98.06% and 47.84%, respectively. Meanwhile, the SH and CH models outperform the baseline design in both ${\delta }_p$ and SEA by 4.63–11.86% and 13.14–49.37%, respectively. These results indicate that the balanced models yield optimal design points that balance multiple objectives and still surpass the baseline design in at least one objective.

The optimized design of ASP from Table 8 was selected by comparing it with other optimized designs. This is carried out to assess the effectiveness of the optimization results. In this study, the optimized design used for comparison is the optimized aluminum foam sandwich panel (AFSP) researched by Pratomo et al. [52]. The optimized AFSP shares the same dimensions as the ASP with a relative density of 0.6. The material of aluminum foam was modeled as CRUSHABLE_FOAM. The design point from the Pareto front with the same mass as the optimized AFSP was selected as the optimized design for ASP. Figure 16 shows the comparison between the optimized ASP and AFSP. The selected optimized design of each auxetic geometry is marked with a black “X,” while the optimized AFSP is marked with a yellow diamond in Figure 15.

As shown in Figure 15 and Table 8, the optimized ASP for all auxetic geometries outperforms the optimized AFSP in SEA by 133.33–156.86%. However, some auxetic models exhibit worse ${\delta }_p$, specifically REH and CH, by 16.21% and 0.13%, respectively. In contrast, DAH and SH models produce a better ${\delta }_p$ by 36.71% and 6.49%, respectively. Overall, the optimized ASP results from the Pareto front, with the same mass as the optimized AFSP, demonstrate superior blastworthiness. Furthermore, the DAH configuration in Table 8 shows promising performance compared to the optimized AFSP and other auxetic geometries.

Figure 17 shows the comparison of the Pareto front for all auxetic models. For small values of ${\delta }_p$ or SEA, the DAH exhibits the best blastworthiness performance, whereas the REH shows the poorest performance. The SH slightly outperforms the CH in terms of smaller ${\delta }_p$ or SEA. For large vales of ${\delta }_p$ or SEA, the REH and DAH share a similar Pareto front, while the CH performs better than the SH for very large values of ${\delta }_p$ or SEA. Overall, the DAH structure demonstrates the best blastworthiness performance among the structures, followed by the REH, SH, and CH. These results align with the findings from the experiment conducted by Chen et al. [35], where the DAH structure outperformed the REH structure.

3.3.2. Blastworthiness Analysis

Figure 18 illustrates the deformation of four models with varying ${\rho }_r$, from the lowest (left) to the highest (right), at 5 ms. The ideal model with maximum SEA shows significant compression in the central area, whereas the ideal model for minimum ${\delta }_p$ exhibits minimal compression. As seen in Figure 18 and Table 8, the ASP with higher values of ${\rho }_r$ tends to exhibit smaller values of ${\delta }_p$, whereas those with lower ${\rho }_r$ tend to have larger values of ${\delta }_p$. This correlation between ${\rho }_r$ and ${\delta }_p$ is further supported by the SHAP values depicted in Figure 13, where higher values of design variables correspond to smaller ${\delta }_p$. These visual representations illustrate that increasing the design variables or ${\rho }_r$ results in a stiffer structure, thereby enhancing its resistance to blast impulses. However, higher values of ${\rho }_r$ also lead to lower SEA values due to the increased mass of the ASP.

As seen in Figure 18, global negative Poisson’s ratio (NPR) behavior is not consistently observed in the ASP. During high-velocity blast impulses impacting the ASP center, the upper part of the auxetic core collapses rapidly, limiting full shrinkage. NPR behavior is more likely to occur in models with lower ${\rho }_r$, where unit cells are less constrained by neighboring cells, allowing freer deformation. Another influencing factor in NPR behavior of the ASP is the corner angle or node radius in CH models. Increased corner angles tend to promote NPR behavior, while smaller node radii increase the likelihood of NPR behavior in the CH because slender ligament struts make it easier for nodes to rotate. Figure 19 shows the displacement vector of the ASP when subjected to air-blast loading. From Figure 18 and Figure 19, REH and CH show the most dominant NPR behavior, followed by DAH, while SH very rarely exhibits NPR behavior. In general, auxetic structure configurations enhance blastworthiness as the material flows toward the impact zone (central area), as seen in Figure 19.

Based on the global sensitivity analysis results, as indicated by the SHAP values (Figure 12 and Figure 13), the most influential design variable in ASP under air-blast loading is the thickness (t). This finding aligns well with the analytical solutions for dynamic crushing stress or plateau stress of auxetic structures, as discussed in studies by [66,67]. According to the formula for dynamic crushing stress in structures such as REH and SH, the variable t exhibits a quadratic relationship, reflecting its role in the plastic bending moment: $M_p={\sigma }_{y}b{t}^2/4$, where ${\sigma }_y$ is the yield strength and b is the strut length. Therefore, increasing t results in a quadratic increase in crushing stress. Additionally, the variable $N_y$ significantly influences blastworthiness performance. Figure 14 illustrates that design points with higher $N_y$ values tend to cluster near the Pareto front. This is attributed to the increased number of cells stacked at the central blast point, enhancing ASP resistance. In contrast, variables $N_z$ and $\theta /r$ show opposite effects compared to $N_y$, with higher values typically located in the bottom-right region, away from the Pareto front. Although increasing $N_y$ and $\theta /r$ reduces ${\delta }_p$, this reduction is relatively minor and primarily affects SEA reduction, as depicted in Figure 13.

3.4. Application of the Optimized Auxetic Sandwich Panel for Armored Fighting Vehicle Protection

In this section, the optimized ASP DAH is applied to the subsystem of an armored fighting vehicle (AFV) to achieve blastworthiness according to NATO STANAG 4569 level 3b standards [68]. Based on documents from AEP Volume 2, evaluating armored vehicles against blast threats involves three steps: (1) structural integrity test, (2) occupant survivability test, and (3) occupant survivability test overmatch. In the first step, the structure of an armored vehicle is tested for integrity under blast-load conditions. The structure passes this test if there are no holes or cracks that increase the likelihood of blast penetration into the passenger compartment. In this study, the addition of ASP protective structures to armored vehicles was evaluated based on the requirements of the structural integrity test.

3.4.1. Numerical Model of an Armored Fighting Vehicle

Figure 20 shows the FE model of an AFV subsystem subjected to an 8 kg TNT blast. This model is adapted from reference [52] based on NATO STANAG 4569 [68] and simulated to evaluate the blastworthiness of the additional protective structure on the AFV. The structural integrity of an AFV includes the floor or occupant side plate (OSP), struck side plate (SSP), lower hull, and vehicle mass, which are supported by four holders. In this case, the vehicle mass is set to 14,000 kg and modeled as a rigid part. The bottom side of the support and steel pot are fixed. The 8 kg cylindrical TNT with diameter of 270 mm and height of 90 mm is modeled using the SPH method, which was validated in Section 3.1.2. INTITAL_DETONATION is defined at the center of the charge with the SoD of 0.8 m to the floor. The materials chosen for the OSP and lower hull are high-strength martensitic steels of 1500 T, while hot-rolled (HR) steel 500 is used for the SSP. The material properties for both 1500 T and HR 500 are given in reference [52]. The thickness of the OSP and SSP are 10 mm and 6 mm, respectively.

To assess the blastworthiness of the ASP as a protective structure, a comparative evaluation was conducted against other structures, which are an AFV without additional protective structures and an AFV equipped with the optimized AFSP. For the AFV without protective structures, simulations were performed with OSP thicknesses of 10 mm and 20 mm. The optimized AFSP configuration used in this evaluation is the optimum configuration from the optimization process conducted by Pratomo et al. [52], while the configurations of the optimized ASP DAH are presented in Table 8.

3.4.2. Simulation Results

Figure 21 shows the simulation results of the AFV subsystem deformation under an 8 kg TNT blast. All deformation and failure modes that occur in the OSP and SSP are consistent with the failure modes summarized in the research by Jacob et al. [69]. As seen in Figure 21a, the 10 mm OSP without a protective structure fails to meet structural integrity requirements due to severe perforation. According to the Cockcroft–Latham damage failure criterion, petaling damage occurs in the middle of the OSP, allowing blast fragments and debris to enter the occupant compartment. Conversely, other configurations, including the 20 mm OSP without a protective structure, maintain structural integrity. As seen in Figure 21b–d, no perforations occur in the OSP. Failure mode I, i.e., large inelastic deformation, occurs in the OSP for the AFV subsystem with 20 mm OSP and protective structure. However, the 20 mm OSP without a protective structure suffers very large deformation compared with the AFV subsystem with the protective structures. From Figure 21c, the SSP of the AFV subsystem with optimized AFSP experiences failure mode IIc, i.e., completely tearing at the center area and forming a ring shape or capping. A different failure mode is observed in the SSP of the AFV subsystem with optimized ASP DAH, where the SPP experiences failure mode II*c, i.e., partial tearing in the central area of the blast, as seen in Figure 21d.

Figure 22 shows the dynamic responses of the AFV subsystem both with and without protective structures, focusing on the midpoint displacement and acceleration of the OSP. The results demonstrate that applying protective structures to AFVs significantly reduces the displacement and acceleration of the floor or OSP. Furthermore, the optimized ASP DAH structure decreases OSP deformation more effectively than the 20 mm single plate and the optimized AFSP. The optimized ASP DAH structure achieves a maximum displacement reduction of 56.99% and 39.00% compared to the 20 mm single plate and optimized AFSP, respectively. Additionally, it reduces the maximum acceleration by 52.55% and 43.56% relative to the single plate and optimized AFSP, respectively, indicating its superior performance.

Figure 23 presents the energy absorption (EA) and SEA of the AFV subsystem with various structural types. It can be seen that in the AFVs with protective structures, the core structure absorbs the most energy, followed by the SSP and OSP. Additionally, the energy absorbed by the optimized ASP DAH exceeds that of both the 20 mm single plate and the optimized AFSP. Notably, the optimized ASP DAH structure achieves SEA improvements of 335.00% and 48.30% compared to the 20 mm single plate and the optimized AFSP, respectively. This highlights the superior effectiveness of auxetic structures in absorbing energy under localized loads compared to other structures.

Table 9. Blastworthiness parameters of the AFV subsystem without and with various types of protective structures.

Parameter	Single Plate (20 mm)	Optimized AFSP	Optimized ASP DAH
Max. displ. (mm)	156.56	110.40	67.34
			(56.99% */39.00% **)
SEA (J/g)	0.60	1.76	2.61
			(335.00% */48.30% **)
Max. acc. (1000G)	118.01	99.20	55.99
			(52.55% */43.56% **)

* Improvement of optimized ASP DAH with respect to single plate. ** Improvement of optimized ASP DAH with respect to optimized AFSP.

summarizes the blastworthiness parameters of AFV subsystems with and without protective structures.

Overall, the ASP, particularly ASP DAH, demonstrates significant potential for application as a protective structure in AFVs. Compared to conventional structures such as honeycomb and aluminum foam, auxetic structures exhibit superior performance, a finding consistent with other studies [32,33,37]. The enhanced performance is attributed to the unique deformation mechanism of auxetics, specifically their NPR, which causes the structure to densify in regions near the blast source. This densification increases the energy absorption capacity and enhances the blast resistance of the structure. Additionally, the layered arrangement of the auxetic structure in the direction of the blast allows it to progressively absorb impact energy from the outermost to the innermost layers, thereby mitigating the direct effects of blast loads on occupants or OSP. However, a significant challenge in utilizing ASPs as protective structures lies in their manufacturing process. For large-scale production, auxetic cores with basic configurations, such as REH, DAH, and SH, can be manufactured using sheet-folding methods. However, this method requires more effort and time [70] compared to the manufacturing of conventional honeycomb or aluminum foam, leading to higher production costs. Furthermore, ensuring effective bonding between layers presents another challenge in the ASP manufacturing process, as inadequate bonding can compromise the structural integrity and blast-resistance performance of the panels. Future research should focus on addressing the manufacturing and cost-related challenges associated with the large-scale application of ASPs as protective structures.

4. Conclusions

The design and multi-objective optimization of four types of auxetic sandwich panels (ASPs), including re-entrant honeycomb (REH), double-arrow honeycomb (DAH), star honeycomb (SH), and tetra-chiral honeycomb (CH), under air-blast loading were successfully conducted using a machine learning approach. The blastworthiness performance of ASPs was evaluated using the finite element method. The optimization objectives were to reduce permanent displacement (${\delta }_p$) and increase specific energy absorption (SEA). The nondominated sorting genetic algorithm II (NSGA-II) was employed for optimization, utilizing an artificial neural network (ANN) metamodel. The results of the optimization process are summarized as follows:

• The ANN metamodel that was formed was proven to be accurate in predicting the blastworthiness performance of ASPs. The optimization process using NSGA-II produces optimal designs efficiently. The optimization results shows that the permanent displacement conflicts with the SEA.
• Global sensitivity analysis using the SHapley Additive exPlanations (SHAP) method indicates that cell thickness in ASPs is the primary factor influencing blastworthiness performance, significantly affecting stiffness and plastic-bending moments within auxetic cells. Meanwhile, the corner angle and node radius in the CH model are identified as the least influential variables.
• The configuration of auxetic structures enables effective energy absorption, enhancing the blast resistance of sandwich structures, as the material flows toward the impact zone. The occurrence of global negative Poisson’s ratio (NPR) behavior in ASPs under air-blast loading is influenced by rapid and high-velocity blast impulses, leading to localized unit cell collapse near the blast source. Generally, the REH and CH models exhibit more dominant NPR behavior at lower relative densities.
• Multi-objective optimization substantially enhances blastworthiness performance. Compared to baseline models, optimization achieves SEA improvements ranging from 26.21% to 83.68% for equivalent permanent displacements and reduces permanent displacement by 12.55% to 35.60% for equivalent SEA. The balanced REH and DAH models outperform baseline models in SEA only by 98.06% and 47.84%, respectively, while SH and CH models exhibit improvements in both permanent displacement (4.63–11.86%) and SEA (13.14–49.37%). Furthermore, the optimized ASP outperforms the optimized aluminum foam sandwich panel (AFSP) in SEA by 133.33–156.86%. However, only optimized DAH and SH produce better permanent displacement reduction by 36.71% and 6.49%, respectively. Among the four auxetic configurations, the DAH structure offers the best blastworthiness performance.
• ASPs demonstrate promising results in armored fighting vehicle (AFV) applications, meeting structural integrity requirements effectively. In dynamic response scenarios, optimized ASP DAH significantly reduce maximum displacement and acceleration of occupant floors by 39.00–56.99% and 43.56–52.55%, respectively, compared to other structures. The AFV subsystem incorporating optimized ASP DAH achieves a 48.30% increase in SEA over optimized AFSP and a remarkable 335.00% increase over single plates. This indicates that the auxetic core exhibits superior energy absorption capabilities compared to aluminum foam cores and single plates.

The preliminary results of this study indicate that ASPs have significant potential as blastworthy structures for broader applications, including defense in buildings, automotive, and aerospace industries. However, further research is required to comprehensively evaluate the ASP response in practical settings, as this study does not address the manufacturing aspect or experimental testing and relies solely on empirical methods for modeling air-blast loading. Additional investigations, such as nonlinear structural analysis and acoustic analysis due to blast impact, can also further clarify the blast response of the structure. Future research could incorporate additional objectives, such as minimizing acceleration and reducing the injury risk of an occupant, as part of a multi-objective optimization framework. Additionally, using explainable machine learning in metamodels could provide clearer insights into how ASP design variables affect blastworthiness.