The capacity of a structure to absorb the energy generated from an impact or crash is called the energy absorption in crash testing. Several crash tests involve subjecting a material or structure to collision with an object at a predetermined velocity, and the absorbed energy is subsequently determined (Guler et al., 2010; Sun et al., 2017). This aspect is of the utmost importance in designing structures or vehicles to minimize human injuries and vehicle damage. Eq. 1 is the formula for calculating the energy absorption (EA): E A δ = ∫ 0 δ F δ d δ (1) where F is the crushing force, and δ is the displacement. Therefore, the SEA is calculated using Eq. 2: S E A δ = E A δ m = ∫ 0 δ F δ d δ m (2) where m is the mass of the specimen. Crash test data can be used to calculate the MCF. This value is obtained using Eq. 3: M C F = ∑ i = 1 N F i N (3) where N is the number of data points, and i is the index for the i th data point. The CFE is calculated using Eq. 4. C F E = M C F P F (4) The PF is obtained from the force–displacement graph. This represents the maximum force experienced during collision, which typically occurs at the beginning, when the highest impact velocity is reached. The PF measurement plays a crucial role in evaluating the EA capability of a material or structure and in predicting potential damage.

Validation was conducted to ensure the accuracy, reliability, and suitability of various simulation results serving as information for researchers to make decisions. It was achieved by comparing the simulation results with the data obtained from other trusted sources. A validation simulation typically involves a top–down arrangement consisting of an impactor, a thin-walled square tube, and a fixed base. Figure 1 illustrates the meshing and overall arrangement of a rigid shell plane with several 5 mm elements used in the energy-absorbing structure.

The cross-sectional measurements of the square tube were derived from the typical average perimeter of the tubes used in sedans (Tarlochan et al., 2013). The simulated tube had a 300 mm perimeter, 350 mm length, 1.7 kg mass, 75 mm × 75 mm major dimension, and 2 mm thickness. The impact simulation results are depicted by a curve showing the crushing force versus displacement (Figure 2). An insignificant difference exists between the conducted simulation and reference simulation results. In this study, the significance of the differences in the calculated values of the generated curves was determined.

Based on the data processing, it was found that the most significant difference in the crashworthiness indicator values between the simulation results and the reference was 3.9% for the PF, as listed in Table 1. This difference may not be significant; therefore, we may further investigate various optimization processes for tube structures incorporating the crush initiator.

The authors observed the occurrence of the asymmetric mixed collapse mode B in the current problem. This mode was determined using the analytical equations derived by Abramowicz and Jones (1984), which provide the MCF values for a thin-walled square tube and the first-fold wavelength of the basic collapse element, expressed by Eqs 5, 6, respectively: M C F = M 0 1 + 0.41 V c t D 1 / p 46.16 c t h 1 / 3 + 2.14 c t h 2 / 3 + 1.30 (5) 2 H = 1.66 c t 2 / 3 h 1 / 3 (6) where M 0 = σ 0 h 2 / 4 , σ 0 is the flow stress, h is the thickness of the tube, V is the velocity of the impactor, c t is the length of a side of the square tube, D = 6844 s^-1, p = 3.91 , and H is half of the initial distance of the wavelength. Therefore, the simulation results can be cross-referenced with analytical calculations to validate the progressive buckling behavior of a thin-walled square tube obtained via numerical simulations. Eq. 7 provides the strain rate (Abramowicz and Jones, 1984) resulting from the dynamic crushing of a square tube undergoing an asymmetric mixed collapse mode B: ε ˙ = 0.41 V / c t (7) where ε ˙ is the strain rate.

For simplicity, the analytically calculated values of the MCF and first folding wavelength were 113.02 kN and 37.2 mm, respectively. The flow stress was determined using the Johnson–Cook constitutive model for mild steel A36 material (Lacy, 2010), as shown in Figure 3. The material exhibited high sensitivity to strain rate, with increased strain rates resulting in increased flow stress. From Eq. 7, the strain rate generated during the drop impact scenario was approximately 85 s^-1, yielding a corresponding flow stress of 421.24 MPa. The initial distance measured between the upper point ( P u ) and lower point ( P l ) of the first folding wavelength was 40 mm, as shown in Figure 4. The difference in the wavelength of the first folding between the numerical results and the calculated values obtained using Eq. 6 was smaller than the meshing size of 5 mm. Therefore, the numerical results for the impact on the square tube in this study show satisfactory agreement with the analytical predictions.

The experimental design was used to examine the relationship between the experimental parameters and responses. Three parameters were selected as input variables, and their effects on the responses were observed using a face-centered central composite design (CCD). The three factors analyzed were the height, width, and distance between the center point of the crush initiator and the top of the tube, as shown in Figure 5, as well as the detailed geometry illustrated in Figure 6. An empirical correlation was established, based on the results of the experimental design, which were statistically analyzed using an RSM. RSM is a reliable analytical method that employs a second-order polynomial for investigating the effects of these input variables (Sahoo, 2011; Nagaraju et al., 2016; Deshwal et al., 2020). The analysis of variance (ANOVA) was employed to examine the significance of the model and input parameters by calculating the mean squared, degrees of freedom, and sum of the squared deviations for each input (Maqsood et al., 2022). A significance level of 5% was applied to determine the significant input values, and the p-value was used to assess the acceptance or rejection of the null hypothesis by the model. A lower p-value provides more significant evidence to reject the null hypothesis, indicating a more compelling indication that the model is significant. The key aspect of the RSM is the derivation of various mathematical formulas that express the relationship between the response and input variables (Dadrasi et al., 2020). This relationship is expressed as a quadratic equation (Eq. 8): R A , B , C = ¥ 0 + ∑ ¥ i A + ∑ ¥ j B + ∑ ¥ k C + ∑ ¥ i j A B + ∑ ¥ i k A C + ∑ ¥ j k B C + ∑ ¥ i i A 2 + ∑ ¥ j j B 2 + ∑ ¥ k k C 2 (8) where R is the proposed expected response, ¥ 0 is the regression constant, ¥ i , ¥ j , and ¥ k are the linear coefficients, ¥ i j , ¥ i k , ¥ j k are the first-order interaction effects, ¥ i i , ¥ j j , and ¥ k k are the square terms of each factor, and A , B , and C are the input or predictor variables. By generating a three-dimension surface plot based on the derived formula, the optimization process of a square tube structure with an elliptical crush initiator was enabled.

The low and high values of the height were 5 and 15 mm, respectively. The low and high values of the width were 10 and 35 mm, respectively, and those of the distance were 15 and 40 mm, respectively. A CCD with six center points suggests conducting 20 runs with varying input parameters for the geometry and position of the crush initiator in the square-tube structure. The response variables of interest are the PF, MCF, CFE, and SEA, as shown in Figure 5.

After the validation was completed, the square tube structures with the optimized elliptical crush initiators were designed, based on Figure 1. The design of the elliptical crush initiator was based on three geometric parameters: height ( A ), width ( B ), and the center-to-center distance of the cut-out from the top end of the square tube ( C ). Impact tests were conducted at an impactor velocity of 15.6 m/s and a mass of 275 kg to evaluate the crashworthiness of the tube (Witteman, 1999; Tarlochan et al., 2013). Figure 6 depicts the geometry of a square tube with an elliptical cut-out in a line on each side.

In this study, A36 mild steel was used as the tube material. This material has the following parameters (Lacy, 2010): initial yield a = 146.7 MPa, strain hardening b = 896.9 MPa, strain hardening exponent n = 0.320 , strain rate coefficient c = 0.033 , thermal softening M = 0.323 , reference strain rate ε ˙ 0 = 1.0 s^-1, Young’s modulus E = 200 GPa, and density ρ = 7850 kg/m³. A nonlinear finite element code, ABAQUS-Explicit, was used to establish finite element models of square steel tubes. Homogeneous continuum shells with five integration points in the thickness direction of the element were employed to model the square-tube structure. A 5-mm-sized mesh consistent with the square-tube model without a crush initiator, was used. A general contact-component interaction model was employed to simulate the contact between two or more surfaces in contact or adjacent to each other. General contact models are widely used in crash simulations, structural dynamics analyses, and manufacturing process simulations. The contact between the impactor and square tube was modeled as a finite sliding penalty with a Coulomb friction coefficient of 0.2 (Altin et al., 2017). A tie constraint was applied to prevent any relative movements between the square tube and fixed base, serving as an interaction between the two components.

Four different responses were analyzed in this study: PF, MCF, CFE, and SEA. These responses were analyzed by varying the height, width, and center distance of an elliptical crush initiator at the top end of a square steel tube. The experiment was conducted using a face-centered CCD for each input parameter. The results are presented in Table 2. Design-Expert software (Anderson and Whitcomb, 2017) was employed to determine the square regression coefficients ( R 2 ) of the developed model, perform regression fitting, and conduct the response surface study (Kallath et al., 2021). The square regression coefficients ( R 2 ) of the PF, MCF, CFE, and SEA were 0.9308, 0.9939, 0.7200, and 0.9975, respectively. The quadratic-order models developed to represent crashworthiness performance were good, as evidenced by the high R 2 values for each indicator, which were close to 1.

Table 3 lists the ANOVA values for all responses and the suggested model for each response, which involves a quadratic order process to calculate the sum of squares, F-values, and p-values for each input parameter, their squares, and interactions. This table lists data on the sum of the squares of the model, residuals, and lack of fit for each response. Based on the data, the PF model is statistically significant, with an F-value of 14.95 and a corresponding p-value <0.05, indicating its sensitivity to any changes occurring in the input. The PF model is expressed by Eq. 9. P F = 224.15 − 1.11 A − 6.02 B + 8.21 C − 2.81 A B + 3.2 A C + 3.01 B C − 0.5263 A 2 − 5.25 B 2 − 8 C 2 (9)

The statistical analysis of the MCF model reveals an F-value of 182.39 and a corresponding p-value <0.05, indicating that the model is highly sensitive to changes in the input. The MCF model is expressed by Eq. 10. M C F = 119.94 + 0.9632 A − 2.2 B + 4.66 C + 0.5647 A B − 1.49 A C + 2.33 B C − 1.46 A 2 − 1.67 B 2 − 2.76 C 2 (10)

The CFE model had an F-value of 3.54 and a corresponding p-value <0.05, indicating sensitive to changes in the input. The CFE model can be expressed by Eq. 11. C F E = 0.5351 + 0.008 A + 0.0054 B + 0.0001 C + 0.0107 A B − 0.0159 A C − 0.0055 A 2 + 0.0054 B 2 + 0.008 C 2 (11)

Statistical analysis of the SEA model revealed an F-value of 434.93 and a corresponding p-value <0.05, indicating that the model is highly significant and sensitive to changes in the input. The SEA model can be expressed by Eq. 12. S E A = 14.53 + 0.1016 A − 0.2106 B + 0.5344 C + 0.0469 A B − 0.1730 A C + 0.3080 B C − 0.1666 A 2 − 0.1913 B 2 − 0.3166 C 2 (12)

Figure 7 shows the importance of each input parameter, square, and interaction for all the four responses. The results indicate that in the case of the PF, input parameter C 2 was the most significant, whereas A 2 was the least significant ( C 2 > C > B > B 2 > A C > B C > A B > A > A 2 ). In the case of the MCF, input parameter C was the most significant, whereas A B was the least significant ( C > C 2 > B > B C > B 2 > A 2 > A C > A > A B ). In the case of CFE, input parameter A C was the most significant, whereas B C was the least significant ( A C > A = A B = C 2 > B = A 2 = B 2 > C > B C ). In the case of the SEA, input parameter C was the most significant, whereas A B was the least significant ( C > C 2 > B C > B > B 2 > A 2 > A C > A > A B . These findings are consistent with those of Huang et al. (2010), indicating that the position of the elliptical cut-out plays a crucial role in determining the SEA and reducing the PF. Moreover, it is considered optimal when the cut-out is closest to the impact source. Rogala et al. (2021b) reported that when the cut-out was placed at the same location, the cut-out width played a more significant role in decreasing the PF and increasing the SEA than the height. This study confirms the above findings regarding this statement. However, contrasting findings were obtained by Dadrasi et al. (2020), in which the height of the elliptical cutout had a more significant influence than its width. This discrepancy may be attributed to the fact that they only placed a single cutout on one side of the tube, whereas in this study, the cutout was positioned on all four sides of the square tube. Moreover, they adopted a rectangular cross-section, which could potentially lead to distinct crashworthiness behavior. Hence, further analysis of these variations should be conducted in future studies.

The suitability and abnormality of data from the developed model were assessed using several normal plots of the residuals and residual-predicted plots (Anderson and Whitcomb, 2007). A model must have a good fit with the numerical data to be deemed acceptable, and the difference between the predicted and actual values should be minimal and follow a linear trend (Sahoo, 2011; Deshwal et al., 2020; Maqsood et al., 2022), indicating a normal distribution of errors (Chakule et al., 2017). The residual and outlier plots for the PF, MCF, CFE, and SEA models are shown in Figure 8. The data points on the normal plot of the residuals are randomly distributed and do not follow any particular pattern, indicating that the predictions by the proposed model are reliable. All data points in the residual-predicted plot for the entire response model fell within an acceptable range, indicating the absence of outliers, except for two data points related to SEA. The acceptable values for the PF, MCF, and SEA ranged from +4.14579 to −4.14579, whereas those for the CFE ranged from +4.00453 to −4.00453. Despite the presence of a few outliers, the model demonstrated its ability to precisely capture the relationship between the input and response variables with a high R 2 value, even when dealing with noticeably different data points (Hawkins, 2004).

This model has been analyzed and verified, and can help explore the design space. The 3-D surfaces comprehensively represent the responses when the input parameters are altered. Figure 9 shows the results of the numerical simulations and the use of 3-D surfaces and contours to depict the overall response using an RSM. One of the three input variables was kept constant at specific points, because displaying the interaction between more than two variables in a 3-D plot and contour was challenging. The selection of the constant variable was determined based on Figure 7, considering the linear coefficients that exhibited the lowest significance levels. Subsequently, the data points selected for the constant variable generated the highest or lowest values for each response.

Figures 9A,C,G show visual representations of the width ( B ) along the x-axis, distance ( C ) along the y-axis, and responses of the PF, MCF, and SEA along the z-axis. The height variable (A) remained constant because it was considered less significant than the other two input variables in the RSM model. Figure 9A is a 3-D surface graph with the height ( A ) fixed at 15 mm to highlight the range of PF values observed at that point in the data. When the crush initiator was closer to the top of the tube and wider, the PF values decreased. The influence of the height of the crush initiator became more significant as it increased. Figures 9C,G show 3-D surface graphs with the height ( A ) set at 5 mm to demonstrate the range of the MCF and SEA values observed at that specific data point. The MCF values increased when the crush initiator was farther from the top of the tube, and its width was narrower. However, when the MCF values were low at distances very close to the impact source, this can be compensated for by increasing the height of the crush initiator, yielding MCF values even higher than those observed in a tube without a crush initiator. For the SEA response, the values increased with increasing distance of the crush initiator from the top of the tube. When the crush initiator was positioned near the impact source, a reduction in the SEA occurred simultaneously with an increase in the width of the crush initiator. Conversely, when the crush initiator was far from the impact source, the SEA increased with the crush initiator width. Increasing the height of the crush initiator improves the SEA when it is positioned as close as possible to the impact source, but reduces it when the crush initiator is far from the impact source.

Figure 9E shows the height ( A ) along the x-axis, width ( B ) along the y-axis, and the CFE response along the z-axis. The distance ( C ) variable was kept constant because it was less significant than the other two input variables in the RSM model. The distance ( C ) was set at 15 mm to obtain the maximum and minimum values of the CFE at that specific data point using a 3-D surface graph. Longer heights and widths resulted in higher CFE values for crush initiators located close to the impact source. However, for the crush initiators positioned near the impact source, decreasing the width and increasing the height of the crush initiator decreased the CFE. When the crush initiator is far from the impact source and has an extended width, the influence of the height can be disregarded in the changes to the CFE value.

Consequently, designing an elliptical cut-out as a crush initiator in square steel tubes is a complex task because enhancing one crashworthiness performance indicator may compromise another. Hence, the optimization process is crucial for attaining an energy absorber that exhibits the highest performance across all crashworthiness indicators.

Designing an energy-absorbing tube with a square steel shape and crushing it to achieve a high CFE while maintaining an SEA is challenging. The design process must reflect all the crashworthiness performance indices. Typically, the design of an energy absorber aims to minimize the PF and maximize the MCF, CFE, and SEA. These input parameters are set within a specified range to achieve various optimization objectives (Djamaluddin, 2023). The weights and importance factors for all responses were uniform because all crashworthiness performance indicators were regarded as equal in this design. Eq. 13 expresses the expression for multiobjective optimization. Min   P F Max   M C F Max   C F E Max   S E A subject to 5 ≤ A ≤ 15 10 ≤ B ≤ 35 15 ≤ C ≤ 40 (13)

Figure 10 shows the ramp function of the optimization results. By employing an RSM for the optimization, optimal outcomes for the dimensions and positions of the elliptical crush initiator were obtained, with a height of 15 mm, width of 24.784 mm, and distance of 15.08 mm. The optimum values corresponding to these points for the PF, MCF, CFE, and SEA were 200.94 kN, 112.78 kN, 0.564, and 13.70 kJ/kg, respectively. The optimization results revealed an overall desirability ( D ) of 0.635, which served as a metric to gauge how closely a combination of input parameters approached the desired response in the design of an elliptical crush initiator. Desirability scores range from 0 to 1, with a score of 1 indicating that the obtained response matches the desired optimum value (Anderson and Whitcomb, 2017). The optimization results demonstrated that CFE obtained the highest desirability score, whereas PF achieved the lowest score ( d C F E > d M C F > d S E A > d P F ). It should be noted that desirability does not affect the quality of the optimization process, as the main objective is to obtain responses that align with the optimization objectives, rather than to achieve a desirability score of 1.

The dimensions and positions of the optimized elliptical crush initiator were verified through finite element simulations, as shown in Figure 11. Figure 12 shows the deformation patterns that matched the selected points during the impact test. The highest error between the optimization outcome using the RSM and that employing the finite element method was 4.8% for the CFE. This level of disparity is deemed acceptable.

The finite element simulation results for the PF, MCF, CFE, and SEA obtained using the optimized dimensions and position of an elliptical cut-out were 195.06 kN, 115.3 kN, 0.591, and 13.89 kJ/kg, respectively. Based on finite element simulations, it was found that the optimal incorporation of crush initiators determined using the RSM resulted in a 10.12% decrease in PF, a 13.67% increase in CFE, and a 2.23% increase in MCF. However, these outcomes were offset by a slight decrease of 0.82% in the SEA.