#### Abstract

Higher energy absorption efficiency and better crashworthiness performance are always the key objectives for different energy absorbing structures applied in numerous industries including aerospace, rail equipment transportation, and automotive. In this study, a functionally graded thickness (FGT) design method is introduced in the design of a hexagon honeycomb structure to improve energy absorbing efficiency on the basis of a traditional honeycomb with uniform thickness (UT). The validation of a numerical analysis model for a UT honeycomb under axial loading is implemented by a nonlinear finite element code LS-DYNA (V971). Furthermore, the multiobjective crashworthiness optimization of an FGT honeycomb subjected to axial quasi-static compression is conducted to maximize specific energy absorption (SEA) and minimize peak crashing force (PCF). In addition, three surrogate models, including radial basis function (RBF), response surface method (RSM), and kriging (KRG), are compared in the accuracy of predicting SEA and PCF and capacity for optimization design of FGT honeycomb structure; the Nondominated Sorting Genetic Algorithm (NSGA-II) is applied to obtain the Pareto optimal solutions for the maximum thickness, minimum thickness, and thickness variation gradient exponent of a honeycomb wall. The optimal points obtained by different surrogate models subjected to an SEA value of 18.5 kJ/kg, 20 kJ/kg, 22 kJ/kg, and 24 kJ/kg are validated, and corresponding optimal parameters are compared; RBF and RSM are more suitable in crashworthiness optimization design of the FGT honeycomb structure. It is indicated that the FGT honeycomb with optimal geometrical parameters presents remarkable enhancement and energy absorbing potential compared to the traditional honeycomb structure.

#### 1. Introduction

A natural honeycomb consists of numerous hexagon cells, which is prepared for multiple functions of honey storing, living, and reproduction. The section of a cell is constituted as hexagon with minimum material, withstanding external loads and foreign invasion to protect the population safety from destroying. The honeycomb configuration has been widely applied in the energy absorption structure design due to its energy absorbing potential and advantage, which has been validated by some numerical and experimental studies in comparing different section tubes [1, 2]. In general, a commercial metallic honeycomb structure with uniform thickness (UT) is often manufactured by a corrugation process; sheets corrugated are bonded using special binding element to generate many layers. What is interesting is that the actual honeycomb cell shown in Figure 1 is constructed with walls of various thicknesses controlled by a certain graded function, which is different from traditional honeycombs manufactured by corrugation process.

In recent decades, optimization for structural crashworthiness and energy absorption performance has become an important topic of research [3], and there are constant studies of thin-walled structures seeking more efficient utilization of materials to achieve higher energy absorption performance. Aspired by some natural multicell structures, bionic design is widely adopted in energy absorbers [4, 5]. Lightweight and excellent energy sorption indicators have become a common view on energy absorbers design of various industries involving passive protection and safety.

With the development of the method for predicting crush behavior of multicorner prismatic columns under an axial compressive load [6], increasing attention and interests were focused on the energy absorption efficiency and crashworthiness performance for multicell columns, and it is found that they have more energy absorption than a single-cell tube, which is largely due to the generation of more corner elements [7, 8]. With the finding that the tubes with various thicknesses exhibit potential in improving energy absorption performance by some studies [9, 10], the tubes with various thickness distributions have been considered as a new structural configuration with higher energy absorption efficiency of less material utilization. In recent years, functionally graded thickness (FGT) has been proposed and become more and more attractive in the design for energy absorbers. Some interesting works show that an FGT design concept was adopted in different cross-sectional tubes under various loading, including the single-cell tubes with FGT along the longitudinal direction under axial loading [11, 12], circular tubes with symmetrical thickness gradients under a transverse load particularly for three-point and four-point bending [13], and a foam-filled FGT circular tube under oblique loading with different angles [14]. In addition, the multicell square tubes [15, 16] with FGT also exhibit an obvious improvement in crashworthiness performance than traditional UT multicell columns.

The thin-walled structures with functionally graded foam (FGF) and functionally graded strength (FGS) structure also perform well in absorbing crushing energy. FGF filling has been validated as an effective method to improve the deformation stability and energy capacity for the thin-walled structure, and it is found that the energy absorption of multitubes filled with FGF is higher than that of tubes filled with uniform foams [17]. The improvement of energy absorption for tubes using FGF varied with the structural configuration and loading; Koohbor and Kidane [18] studied the crushing response difference of continuously graded and discretely layered foams. Ying et al. [19] proposed an FGS structure, which also brings remarkable promotion in energy absorption capacity.

Some researchers proposed energy structures with functionally graded design in both filler and tubes and generated new structural configurations by assembling. Zhu et al. [20] introduced a double functionally graded (DFG) structure by filling a functionally graded honeycomb in an FGT tube along the axial direction under multiple load cases. Fang et al. [21] proposed an FGT structure filled with FGF along the transverse direction and investigated its crashworthiness performance. All these previous studies indicate that the energy absorption structures considering a functionally graded design method are generally superior to the uniform counterpart in absorbing crushing energy.

It is worthy to explore the energy performance enhancement of a traditional honeycomb with uniform thickness when the FGT was introduced. In this work, the finite element model of the honeycomb structure is established for energy absorption performance analysis, of which the element size is determined by convergence analysis of specific energy absorption (SEA) and mean crushing force (MCF). Then, three surrogate models are established by analyzing samples with different parameters in LS-DYNA, the multiobjective optimization problem of the honeycomb structure is solved using the Nondominated Sorting Genetic Algorithm (NSGA-II) to maximize the SEA and minimize the PCF, and the validation results by partial optimal points indicates that the FGT honeycomb structure has an excellent energy absorption potential than that with uniform thickness (UT).

#### 2. Finite Element Modeling

##### 2.1. Geometrical Parameters

As the honeycomb structure is generally constructed with a repeated hexagon pattern and structural scale is determined by usage requirement, here we established a partial structure constructed with seven cells. The cell wall width is set to 15 mm, and its cross section and functionally graded thickness are described in Figure 2.

**(a)**

**(b)**

Figure 3 presents the function curves of normalized thickness to normalized distance when the gradient exponent varied, its function expression is given in equation (1), and the exponent is the key parameter controlling the thickness distribution of honeycomb wall:

The graded thickness can be computed by equation (1), where is the distance to the middle point of the width and is the half-width of the honeycomb cell (Figure 2). and are the maximum and minimum thicknesses along the half-width of cell wall, respectively. The gradient exponent is defined to change the various trends of thickness.

In this paper, the gradient is set to vary from 0.1 to 10. Figure 3 shows the variation trends of normalized graded thickness distributions along the width under different gradient exponents .

##### 2.2. Crashworthiness Criteria

Generally, multiple indicators are adopted for evaluating crashworthiness of energy absorbing structure, for instance, energy absorption (EA), specific energy absorption (SEA), peak crash force (PCF), mean crash force (MCF), and crash force efficiency (CFE). In particular, several indicators are commonly grouped as objectives in multiobjective optimization problems.

An absorbing structure was expected to absorb more energy with higher energy absorption efficiency as possible in a controllable crushing pattern, and meanwhile, the peak forces level in crash should be reduced as possible to prevent survival zone safety of the passengers from overloading impact. Except peak crush force is directly obtained from force-time curve, other indicators need to be calculated though certain mathematical formations given in equations (2)–(5).

The energy absorption of structure describes the capacity of impact energy absorption, of which mathematically formation can be written aswhere is the impact force as a function of the displacement and is the maximum displacement.

The specific energy absorbed is the key indicator to evaluating the energy absorption by unit mass, and it can be given aswhere is the structural mass participating energy absorption.

From the perspective of energy absorbing stability, the force curve with no obvious fluctuation is generally considered as an ideal energy absorbing structure. Thus, the closeness degree between PCF and MCF has become an important point to indicate higher crush force efficiency. The MCF and CFE can be described as

##### 2.3. Finite Element Model

The finite element constructed for the FGT honeycomb is shown in Figure 4; the cross section consists of seven cells, and the height of structure is 120 mm; the numerical analyses were carried out using the nonlinear explicit analysis code LS-DYNA. The model was established by using Belytschko–Tsay four-node shell elements with five integration points along the thickness direction, the top block with mass was set to compress the honeycomb structure at a constant velocity of 2 m/s, ignoring the strain rate effect, the crushing displacement was set to 80 mm, the bottom block fixed provides vertical displacement constraint for the honeycomb.

In order to deal with the contact condition between the honeycomb structure and the mass block, CONTACT_AUTOMATIC_SURFACE_TO_SURFACE and CONTACT_AUTOMATIC_SINGLE_SURFACE models were employed in dealing with the contact behavior of mass block to honeycomb and self-contact, respectively. The static and dynamic friction coefficients were set to 0.15 for considering the energy dissipation by the sliding surface.

The honeycomb structure was modeled using the material of AA6061O [22], the power law-plastic model with strain hardening (MAT18 in LS-DYNA) was adopted to characterize the mechanical behavior during deformation. The stress-strain curve of the material is given in Figure 5, and the mechanical properties of AA6061O are given as follows: Young’s modulus , initial yield stress , the ultimate stress , Poisson’s ratio , and power law exponent .

##### 2.4. Convergence Analysis

It is found that element size affects the analysis results especially for dynamical solution. The quality of the finite model directly determines the accuracy and reliability of solution, so it is necessary to conduct the model validation under different analytical parameters. Here, we adopted structure with a uniform thickness of 0.2 mm, element size was set to vary from 5 mm to 1.5 mm, and the structure was crashed at a constant velocity of 2 m/s with nonlinear finite analysis code LS-DYNA. The SEA and PCF were chosen as the indicators to analyze the effect of element size on stability and accuracy for the finite model. Figure 6 shows that the SEA and PCF tend to converge when the element size is less than 2 mm. When the element size varied from 2 mm to 1.5 mm, the SEA increased from 9.60 kJ/kg to 9.65 kJ/kg, while the PCF decreased from 7.82 kN to 7.64 kN; the percentage change of SEA and PCF is 0.52% and 2.30%, respectively. It was an acceptable level to ensure the convergence for the finite model analysis; as a result of comprehensively considering the accuracy of FEA results and calculation efficiency, the element size of 1.5 mm was finally determined as the optimal element size in the finite model establishment.

#### 3. Multiobjective Optimization Analysis

##### 3.1. Optimization Definition

The energy absorbing structure is commonly expected to present excellent performance in crushing, for instance, absorbing as much impact energy as possible with less material usage. SEA is a key indicator to be considered and set as one of the optimization objectives. Furthermore, higher PCF should be avoided and constrained under an acceptable level in the crashworthiness design considering the safety of survival space. In this study, the optimization problem was created to obtain the optimal parameters for the FGT honeycomb cell structure with multiobjectives of increasing the special energy absorption and decreasing the peak crush force.

The indicators of SEA, PCF, and CFE were frequently set as objectives in multiobjectives optimization problems [23, 24]. In this study, the optimization problem for the FGT honeycomb can be written as the following multiobjective optimization form considering the two different design criteria:where and are the maximum and minimum thickness along the half-width of the honeycomb cell wall (Figure 2), respectively. denotes the gradient exponent of various thicknesses.

##### 3.2. Surrogate Model

The high computational cost of numerical simulations caused by increasing model scale and solution complexity has limited the solution of optimization design problem, and the surrogate models (also called metamodel) were developed and widely adopted in engineering design as the benefit of well accuracy and high efficiency. Without knowing more information of an optimization design problem, it is indeed necessary to test the prediction accuracy of different surrogate models and select the best performance model. In this study, three typical surrogate models commonly used were adopted in predicting SEA and PCF of the FGT honeycomb structure, including radial basis function (RBF), response surface method (RSM), kriging (KRG) model, of which the quartic polynomial was chosen for RSM more accurately than the lower order polynomial.

The sampling in design space is the basis of establishing a surrogate model, and the flowchart of the implementation of the multiobjective optimization design (MOD) for the FGT honeycomb structure is shown in Figure 7.

The design of experiments (DOE) method provides several means to select the sampling points in the design space efficiently. The different samples method is suitable for certain problem, and the optimal Latin hypercube samples (OLHS) [25] method was proposed by adding a distribution criterion to the traditional Latin hypercube designs (LHDs) [26] to improve the distribution uniformity of sample points in the design space, which remarkably improved the accuracy of the surrogate model as its advantage of well space-filling performance.

In this work, experimental design samples were generated using OLHS in the design spaces, and the samples with fine uniformity are plotted with different perspectives in Figure 8. The numerical analysis of sample points is the most time-consuming section in multiobjective optimization design based on the surrogate model. The FEA process costs about 2 min for solving one experimental design point on the DELL workstation of Intel (R) Xeon (R), CPU E5-2630 @ 2.40 GHz. It costs about 5 h in total for completing all the samples points for multiobjective optimization of the FGT honeycomb structure.

**(a)**

**(b)**

To assess the predicting accuracy of three surrogate models mentioned, here we introduced four commonly used metrics [27] in estimating the fitting error, including *R* square (), relative average absolute error (RAAE), and relative maximum absolute error (RMAE) calculated as follows:where is the number of error validation points; is the exact function value from FEA; is the predicted value by the surrogate model; and is the mean of the exact function value from FEA.

The prediction accuracy of surrogate models is significant to guarantee the solution quality of optimization, and it indicates that the optimal result by surrogate models is remarkably different. As listed in Table 1, we calculated four indexes to assess the accuracy of surrogate models in predicting the response of SEA and PCF. In general, higher and lower and indicate better accuracy of surrogate models; it was found that all the three models have satisfied the accuracy demand (), and the RBF model presents higher accuracy than RSM and KRG by comparing the accuracy criterion of three surrogate models.

##### 3.3. NSGA-II

In our study, the Nondominated Sorting Genetic Algorithm (NSGA-II) [28], the improved version by adding a nondominated sorting principle based on NSGA, is used to search for the Pareto fronts of multiobjectives for the FGT honeycomb. The elitist nondominated sorting coupled with the crowding distance is used in obtaining the nondominated set, and the nondominated front is generated after each generation. Finally, the nondominated set is obtained with convergence of the iteration to form the Pareto fronts. This algorithm has been successfully used in solving crashworthiness optimization design problems of an energy absorbing structure as it has an advantage of obtaining well-distributed optimal points [11, 16, 21, 22, 27], which exhibited better effectiveness and remarkable assistance to more determination in parameters design.

The Pareto set point distribution tends to be apart from the Pareto fronts when adopting insufficient population size and number of generation for genetic algorithm in multiobjective optimization problems. It is necessary to adjust the parameters though the quantity and distribution situation of Pareto is set so as to obtain stable Pareto fronts, and better Pareto set distribution, various population sizes, and number of generation were attempted to be determined for NSGA-II to solve the optimization problem of the FGT honeycomb structure.

The three surrogate models were optimized to generate Pareto fronts by NSGA-II, as shown in Figure 9, and we compared the Pareto fronts distribution by four different groups of algorithm parameters; it is indicated that certain amount of Pareto set points were not close enough to Pareto fronts when the population size was insufficient, and the distribution of the Pareto set tends to be considerably improved when the population size or number of generation increased. While the number of generation was set to 100, the number of the Pareto set points is remarkably more than that corresponding to 80 generation.

**(a)**

**(b)**

**(c)**

To solve the optimization problem defined in equation (6), the parameters settings listed in Table 2 were determined for the NSGA-II implemented in this study, which overall considered the analysis efficiency and sufficient and well-distributed Pareto fronts.

##### 3.4. Results and Discussion

It has been demonstrated that different surrogate models may generate different results in the crashworthiness design, so it is necessary to adopting more than one surrogate model for comparison and validation. For instance, RSM is found to be most suitable for the crashworthiness optimization problems of the FGS thin-walled structures [19], and RBF and KRG with a first-order trend is a most applicative model for CFE and SEA prediction for tapered thin-walled tubes with axisymmetric indentations [23], respectively. RBF and quartic RSM perform well in the foam-filled bitubal circular tube, and RBF model is considered to be more suitable in topological configurations of the multicell thin-walled structures [29]. Based on the comparison results, RBF metamodels of SEA and PCF both perform the best in optimization of foam-filled multicell thin-walled structures. It is indicated that there exists considerable effect on the Pareto fronts when different surrogate models were adopted; the modeling accuracy of different surrogate models may not necessarily guarantee the best optimum. In general, it is recommended to adopt different surrogate models and compare their corresponding Pareto fronts for the multiobjectives in the unknown optimization problem, by which the most suitable surrogate model is determined.

As the single surrogate model is not sufficient to find the optimal Pareto fronts, we carried out the optimization using RBF, RSM, and KRG, respectively, by which the Pareto fronts generated are shown as Figure 10. It is apparent that the range area of optimal points by KRG is short than that of RSM and RBF, meanwhile, and there is an obvious distance between the Pareto fronts by KRG and that by another two models, the RSM and obtains a relative coincidence the RBF on the Pareto fronts points distribution except a smaller number of optimal points.

The half-width of the wall is divided into five discrete thicknesses along the direction from to , named as *T*_{1}, *T*_{2}, *T*_{3}, *T*_{4}, and *T*_{5}. The optimal thickness distributions along the cell wall by three surrogate models under four cases are given in Figure 11. There is a better accordance on the various trend of thickness by three surrogate models for the cases of A, B, and C. However, when the SEA is 24 kJ/kg, the thickness distribution of three surrogate models appears with a noticeable difference; in particular, the thickness trend of RBF exhibits a convex trend while those of another two models still are concave. The optimal parameters of maximum thickness, minimum thickness, and exponent are listed in Table 3, as the increasing of SEA and the maximum and minimum thickness presents a positive rise, and the exponent tends to decrease; the sharp variety of thickness is obviously alleviated.

**(a)**

**(b)**

**(c)**

**(d)**

It is necessary to select some typical optimal design points considering certain requirements for the validation of the optimal Pareto set, four typical SEAs of 18.5 kJ/kg, 20 kJ/kg, 22 kJ/kg, and 24 kJ/kg were determined, and the numerical analyses were carried out for the optimal solution by three different surrogate models, and the comparison results are summarized in Table 4. It is found that the prediction accuracies of SEA and PCF by three surrogate models are almost less than 5%, except the higher error of −35.53% and −28.92% in predicting PCF by KRG under the cases corresponding to SEA of 18.5 kJ/kg and 20 kJ/kg, respectively.

##### 3.5. Validation of Optimal Pareto Point

In order to further investigate the improvement of the FGT honeycomb structure in energy absorption, the numerical analyses were performed for the FGT honeycomb and corresponding UT honeycomb to compare the crashworthiness performance. The boundary conditions are the same to the previous settings in Section 2.3; particularly, the number of cells used is enlarged to 22 for reducing the effect from instable deformation of outer wall. The optimal thickness parameters generated by RBF were selected for validation analysis of FGT honeycombs as better accuracy. Comparison of crushing force and SEA of UT and FGT honeycomb structures with 22 cells is shown in Figure 12, and it can be found there appears an instable zone after 20 ms where the crushing force sharply descends for case A; the similar phenomenon exists for case B, but the descending trend tends to be mild, while those of cases C and D are stable without obvious descending.

**(a)**

**(b)**

**(c)**

**(d)**

The summaries of comparison results are listed in Table 5, which indicates that the FGT honeycomb generates an obvious improvement in SEA without remarkable increase in PCF compared to the corresponding UT honeycomb. What is noticeable is that the enhancement of FGT tends to vary with the different SEA values and the energy absorption efficiencies of cases B and C are higher than those of cases A and D.

The instable zones found in Figure 12 are also demonstrated in the deformation mode comparison plot of FGT and UT honeycombs under four typical cases as shown in Figure 13; unlike the symmetrical ring mode of the UT honeycomb, all the FGT honeycombs under different SEAs deform with the diamond mode and the mixed mode under axial loading. There appears obvious irregular deformation for the FGT honeycomb when the SEA is 18.5 kJ/kg, but it tends to be removed in cases C and D, which is mainly due to the increase of the ratio of thickness and length (*T*/*L*) for the honeycomb cell wall. It is important to notice that the SEA is positively related to the maximum and minimum thickness along the cell wall of honeycomb, and the larger thickness gradient exponent can lead to a sharp thickness change, which tend to bring an instable deformation mode particularly observed in Figure 13.

**(a)**

**(b)**

**(c)**

**(d)**

#### 4. Conclusion

Inspired by the natural honeycomb with various cell wall thicknesses, this study introduced the FGT design concept in the traditional honeycomb structure and numerically investigated the energy absorption performance of the FGT honeycomb under axial loading. Three surrogate models including RBF, RSM, and KRG were established for two indicators of SEA and PCF, and the multiobjectives optimization problem was further carried out in obtaining the optimal parameters to maximize SEA and minimize PCF by using NSGA-II.

By concluding the above results, it was found that the RBF and RSM are relatively more suitable for predicting the crashworthiness response for the FGT honeycomb structure. Compared to the UT honeycomb structure, the FGT honeycomb structure exhibits superior energy absorption capacity and efficiency, the certain enhancement percentage for SEA varied as the different design parameters, the maximum and minimum improvement percentages validated for SEA reach 26.65% and 12.96%, respectively. Furthermore, the deformation mode of the FGT honeycomb under axial loading is apparently different from that of the UT honeycomb, four types of the FGT honeycomb with different parameters present diamond or mixed mode, and the thickness gradient has an essential influence on the stable deformation model. The actual promotion of the FGT honeycomb can be further increased by adopting more honeycomb cells in desirable design size. The FGT honeycomb with optimal design parameters can be an excellent energy absorber according to desirable crashworthiness requirements.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

#### Acknowledgments

The authors gratefully acknowledge the financial support of the National Nature Science Foundation of China (no. 11272070), the Research and Development Program of Science and Technology of China Railway Corporation (no. 2015J007-H), the National Key R&D Program of China (no. 2016YFB1200504).