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

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 (AVMs) resulting from explosive remnants of war (ERW) reached 4,710 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 panels (HSPs) 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 HSP. The introduction of hexagonal HSP significantly decreased 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 (PPR) 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 resistant to indentation [20,21], increased shear modulus [22,23], improvement in fracture toughness [24,25,26], and high energy absorption [20,27]. Due to these superior mechanical properties, auxetic materials have found applications in diverse fields such as medicine, biomedicine, sports, and engineering (for more details, see reference [28]). 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. [29] and Yan et al. [30] studied the responses of ASP and hexagonal HSP to blast loading numerically and experimentally, respectively. Both studies demonstrated that ASP exhibits higher blast resistance than conventional HSP due to increased compressive stiffness as the impact zone of ASP shrinks inward. Lan et al. [31] examined the dynamic response of cylindrical ASP, conventional HSP, and AFSP subjected to blast loading. The results showed that cylindrical ASP outperformed both conventional HSP and AFSP across all design parameter combinations. Qi et al. [32] conducted experiments on close-in blast loading on ASP. The results showed that the introduction of ASP effectively prevented damage to the concrete base under close-in blast loading. Chen et al. [33,34] conducted experimental studies on the blast response of ASP 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. [35] 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 ASP under blast loading indicates that 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. These two objectives often conflict with each other. Multi-objective optimization has been demonstrated as effective in addressing this challenge. Due to the complexity and non-linearity, there is no direct relationship between geometric parameters of auxetic structures and blastworthiness performance. Therefore, surrogate models or metamodels are utilized. 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 non-dominated 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. [36] and Wang et al. [37] conducted optimizations of curved and graded AFSPs, respectively, using NSGA-II with ANN and Kriging models as metamodels. Qi and Jiang et al. [38,39] optimized re-entrant ASP using radial basis function (RBF) and NSGA-II. The optimum design produces better performance by reducing maximum deflection until 33% and increase the areal specific energy absorption until 158.6%. Wang et al. [40] and Lan et al. [41] 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 (GSA) 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 schematic diagram of ASP is shown in Figure 1. The ASP consists of two face sheets with thickness of $t_f$ and auxetic core, which has total dimensions of $L\times W\times H$. The L, W, and H is the length, width, and height of the auxetic core in x, y, and z direction, respectively. The auxetic core consists of repetitive 2D unit cell of auxetic geometry with number of cell of $N_y$ and $N_z$ in y and z direction, respectively. The auxetic core is oriented with its cell’s longitudinal direction parallel to the two face sheets, as research by Walkowiak et al. [42] shows this configuration yields better blastworthiness than when the cell’s longitudinal direction is perpendicular to the face sheets. The size of ASP is determined based on the floor size of the AFV and the availability of additional space as a protective structure. Based on reference [43], typical floor size of small AFV is L = 1400 mm, W = 2500 mm, and H = 120 mm. The stand-off distance (SoD) from blast charge position to the front face is 800 mm. However, to reduce the computational time, the size of ASP is reduced into half of it, where L = 700 mm, W = 1250 mm, and H = 100 m with SoD of 400 mm to compensate the smaller size of ASP.

In this study, four auxetic basic geometries were proposed as the core of 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 cell angle $\theta$ and node radius r for CH, and cell thickness t. The variable 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 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. Machine Learning Model

Machine learning (ML) methods are employed to develop surrogate model capable of predicting outputs based on specific input variables. The predicted output $\widehat{\mathit{Y}}=\{{\widehat{y}}_1,{\widehat{y}}_2,...,{\widehat{y}}_n\}$ is related to the input variables $\mathit{X}=\{x_1,x_2,...,x_m\}$ via a black box function $\widehat{\mathit{f}}\left(\mathit{X}\right)$, where m and n is the dimension 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.2.1. Metamodel

This paper 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 ANN is inspired by the working process of 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 layer 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 l-th layer in 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$$

(2)

where k is the number of layer, $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)}$ represents the input and output layer, 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}$$

(3)

$$\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}$$

(4)

$$\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}$$

(5)

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

2.2.2. Sensitivity Analysis

Sensitivity analysis assesses how variations in input values influence output variables, providing insights into the robustness and reliability of a ML model. SHapley Additive exPlanations (SHAP) [44] is an explainability technique based on the Shapley value concept introduced by Shapley [45] 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 non-random [46]. The Shapley value of 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)$$

(6)

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)$$

(7)

where ${\varphi }_0$ is the output of reference input value or ${\widehat{\mathit{f}}}_{\varnothing }\left({\mathit{X}}_{\varnothing }\right)$. For global sensitivity analysis (GSA), averaged SHAP values are utilized to assess feature sensitivity at a global scale [46]. Averaged SHAP values indicate the average impact of input variables on specific outputs. The averaged SHAP of 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|$$

(8)

2.2.3. Multi-Objective Optimization

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

2.3. Numerical Modeling

In this paper, all numerical simulations were conducted using the non-linear explicit finite element (FE) program LS-DYNA.

2.3.1. Geometry, Boundary Conditions, and Contact Modeling

The ASP was modeled as a quarter-symmetric model, as shown in Figure 3, due to 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 ASP has 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 $y-z$ and $x-z$ plane, respectively, while outer face sheets were fixed on the other two sides. The face sheets and auxetic core were modeled using fully integrated shell with five and two integration points, respectively, to minimize the hourglass energy. For data sampling process, thickness of both face sheets is set to 5 mm. The mesh convergence analysis was performed to obtain the mesh size for auxetic core, which converged at mesh size of 3 mm, while 10 mm was used for the face sheets [43]. AUTOMATIC_SINGLE_SURFACE was applied to all parts of ASP for self-contact, while AUTOMATIC_SURFACE_TO_SURFACE was applied between the auxetic core and both face sheets to prevent the penetration. The upper and lower node of auxetic core was assumed to be fully connected with both face sheets using TIED_NODES_TO_SURFACE. Dynamic and static coefficient of friction for all contact is 0.3 and 0.2, respectively.

2.3.2. Material Modeling

Considering its high strain capability, good manufacturability, and low cost, the ASP was modeled using stainless steel 304 material. Due to strain rate dependency in 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)$$

(9)

The first term represents material plasticity, where ${\epsilon }_{eff}^p$ is the effective plastic strain, A is the yield strength, B and n are strain hardening parameters. The second term accounts for 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 coefficients. 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 [14] and [49,50], as shown in Figure 4. These material parameters were summarized in Table 2.

For the material failure criterion, 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 [51]. Mathematically, the Cockcroft-Latham damage criterion can be written as:

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

(10)

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$ is the strain failure. The value of $D_C$ is 1 [43], 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 armored fighting vehicle sub-system.

2.3.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 [52]. 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)$$

(11)

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 range of 0.147 m/kg^1/3$<Z<$ 40 m/kg^1/3, where $Z=R/m^{1/3}$ is called 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 [53,54]. 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}}$$

(12)

where A, B, $R_1$, $R_2$, and $\omega$ are the constant, 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 was given in Table 3.

3. Results and Discussion

3.1. Validation of FE Model

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). CONWEP and SPH models are detailed in Section 2.3.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 Eulerian and Lagrangian part is 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 SoD of 150, 200, and 250 mm. For the CONWEP model, the 0.15 kg C-4 was equated to 0.18 kg TNT based on mass equivalence [56]. The J-C material model parameters from reference [6] were implemented in the numerical simulations. 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 [55]. 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 summarizes the numerical result compared to the experimental results. MMALE shows the smallest error for SoD of 150 and 200 mm, while CONWEP performs best for 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 optimization 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 paper 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] is relevant as it also investigates the blast response of a sandwich panel, despite using a non-auxetic 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 set-up and schematic diagram of trapezoidal corrugated-core sandwich panel. Their study primarily explores the influence of geometric parameters and SoD on mid-point displacements. The horizontal plane have dimensions of 452 × 440 mm, and the exposed area of the blast or the corrugated core dimensions are 300 × 288 mm.

The numerical model utilizes 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 perfect match indicates strong agreement between numerical predictions and experimental data. The average error across all samples, as detailed in Table 5 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 global deflection in the middle area of the corrugated core and local 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 paper, 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, Design of Experiments (DoE) was generated using numerical methods. To construct the metamodel, 160 samples for each auxetic geometry were generated and simulated using FE solver.

3.2.1. Design Variables, Objective Function, and Constraint

Two conflicting objective variables were optimized in this MOOP, which are permanent displacement ${\delta }_p$ of 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. 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 taking into account the geometric limitations. Based on the design variables and constraints, the design space with the corresponding normalized ${\rho }_r$ for each auxetic model is presented in Figure 10. As shown in Figure 10, higher input variables result in 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 non-linearity, while the output layer used a linear activation function to produce continuous values.

The error parameters of each metamodel are presented in Table 7. 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 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$ have 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 of output value with respect to the reference value. In this case, the reference value is the baseline model. 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 plot 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 non-dominated solution. As expected, the Pareto front illustrates the conflict between the two objectives, positioned at the outer side of the prediction sample marked with blue circles. As ${\delta }_p$ increases, SEA also increases, and vice versa.

Several design points from Figure 15 are listed in Table 8. 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.

A comparison between the Pareto front and selected design points was conducted. First, the baseline model, 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 compromised design, marked with the purple pentagon, is compared with the baseline model. A compromised design is a design point on the Pareto front that has a minimum distance to all existing ideal points, which marked with a green "+". As given in Table 8, although the compromised designs for REH and DAH have larger ${\delta }_p$ than the baseline model 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 model in both ${\delta }_p$ and SEA by 4.63-11.86% and 13.14-49.37%, respectively. These results indicate that the compromised models yield optimal design points that balance multiple objectives and still surpass the baseline model in at least one objective.

The optimized design 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 paper, the optimized design used for comparison is the optimized AFSP conducted by Pratomo et al. [43]. 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 are marked with black "X," while the optimized AFSP is marked with 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. [33], where the DAH structure outperformed the REH structure.

3.3. Blastworthiness Analysis

Figure 18 illustrates the deformation of four models with varying ${\rho }_r$, from low (left) to high (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$ tends 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.

Based on the GSA results, as indicated by the SHAP values (Figure 12, 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 [57] and [58]. 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 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.

As seen in Figure 18, global NPR behavior is not consistently observed in 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 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 due to slender ligament struts makes 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 and SH. In general, auxetic structures configuration enhance blastworthiness as the material flows toward the impact zone (central area), as seen in Figure 19.

3.4. Application of the Optimized ASP for Armored Fighting Vehicle Protection

In this section, the optimized ASP DAH is applied to the sub-system of an AFV to achieve blastworthiness according to NATO STANAG 4569 level 3b standards [59]. 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.

Figure 20 shows the finite element (FE) model of an armored vehicle (AFV) sub-system subjected to an 8 kg TNT blast. This model is adapted from reference [43] 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 TNT blast is modeled using the smoothed-particle hydrodynamics (SPH) method, which has been validated in Section 3.1.2. INTITAL_DETONATION was 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 armor steel 1500 T, while HR 500 is used for the SSP. The material properties for both 1500 T and HR 500 is given in reference [43].

To assess the blastworthiness of ASP, a comparative evaluation was conducted with other structures, which are AFV without protective structures and AFV with 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 optimization process conducted by Pratomo et al. [43], while the configuration of the optimized ASP DAH are presented in Table 8.

Figure 21 shows the simulation results of the AFV sub-system deformation under an 8 kg TNT blast. All deformation and failure modes that occurred in the OSP and SSP are consistent with the failure modes summarized in the research by Jacob et al. [60]. 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, petalling 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 occurred in the OSP. Failure Mode I, i.e. large inelastic deformation, occurred in the OSP for AFV sub-system with 20 mm OSP and protective structure. However, the AFV sub-system with 20 mm OSP suffer very large deformation, compared with AFV sub-system with the protective structures. From Figure 21c, the SSP of AFV sub-system with optimized AFSP experiences failure mode IIc, i.e. completely tearing at the center area and forms ring-shaped or capping. A different failure mode is observed in the SSP of the AFV sub-system 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 sub-system both with and without protective structures, focusing on the mid-point 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 single plates and the optimized AFSP. The optimized ASP DAH structure achieves a maximum displacement reduction of 56.99% and 39.00% compared to the 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 absorbed energy (EA) and specific energy absorption (SEA) for the AFV sub-system with different structural types. It can be seen that in 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 single plate and the optimized AFSP. Notably, the optimized ASP DAH structure achieves SEA improvements of 335.00% and 48.30% compared to the single plate and optimized AFSP, respectively. This highlights the superior effectiveness of auxetic structures in absorbing energy under localized loads compared to other structures. Table 9 summarizes the blastworthiness parameters of AFV sub-systems with and without protective structures. Overall, ASP, particularly ASP DAH, demonstrates significant potential for application as a protective structure in AFVs.

4. Conclusions

The design and multi-objective optimization of four types of auxetic sandwich panels (ASPs) under air blast loading were successfully conducted using a machine learning (ML) approach. The blastworthiness performance of ASPs was evaluated using the finite element (FE) method. The optimization objectives were to reduce permanent displacement (${\delta }_p$) and increase specific energy absorption (SEA). Non-dominated 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 performances of ASP. The optimization process using NSGA-II produces optimal designs efficiently. The optimization results shows that the permanent displacement is conflicting to the SEA.
• Global sensitivity analysis (GSA) using the SHapley Additive exPlanations (SHAP) method indicates that cell thickness in ASPs is the primary factor influencing blast-worthiness 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 SEAs. The compromised 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 AFSP in SEA by 133.33-156.86%. However, only optimized DAH and SH that produces 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. AFV sub-system incorporating optimized ASP DAH achieve a 48.30% increase in SEA over optimized AFSP and a remarkable 335.00% increase over single plates. This indicates 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, as this study does not consider in manufacturing aspect of the ASP and only employs empirical methods for modeling air blast loading. 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, explainable machine learning methods could be utilized in metamodel creation to provide more intuitive insights into the relationship between ASP design variables and blastworthiness performance.