Multi-objective robust optimization of foam-filled double-hexagonal crash box using Taguchi-grey relational analysis

In this paper, a novel thin-walled double-hexagonal crash box is first proposed and then multi-objective robust optimized for better overall crashworthiness under multi-angle impact loading, using a proposed hybrid method combining aluminum foam-filling and Taguchi-grey relational analysis (GRA). Specifically, the finite element (FE) models of the regularly-shaped double-hexagonal column (DHC) extracted from original irregularly-shaped crash box under multi-angle impact loading, including hollow (H-DHC) and foam-filled (F-DHC), are first built and validated by experiments. On this basis, a comprehensive crashworthiness comparison is conducted to explore relative merits of F-DHC over original H-DHC under multi-angle impact loading. After that, the F-DHC is multi-objective robust optimized for maximizing overall specific energy absorption (SEAθ) and minimizing overall initial peak crushing force (IPCF0) simultaneously under multi-angle impact loading, using a hybrid method of Taguchi-GRA. At last, a bumper-crash box integrated crashworthiness analysis under multi-angle impact loading is executed to further verify the optimization. The optimal F-DHC and the optimized crash box within the optimal F-DHC demonstrate evident improvement of crashworthiness compared to their respective initial designs, indicating aluminum foam-filling combined with Taguchi-GRA could be an effective approach for multi-objective robust optimization of the novel crash box and other similar vehicle structures.


Introduction
2][3] In the process of a car frontal collision, the collision energy dissipation mainly depends on the body frontend structures composed of a bumper, two crash boxes, and two front rails at the two ends of the bumper (Figure 1).The crash box, which is a classical thinwalled structure (TWS) connecting the bumper and the front rail, is classified as the main crash energyabsorbing member in the front-end structures and is expected to be capable of absorbing kinetic energy through collapse deformation during a frontal impact, as well as maintaining vehicle deceleration within a safe limit. 4,5Hence, crashworthiness investigation and optimization design of thin-walled crash boxes have always been research hotspots.
To improve crashworthiness of TWSs without adding too much weight, scholars worldwide have proposed many effective methods during past decades, such as multi-cell design, 6,7 functionally-graded design, [8][9][10] filling structures (or materials), [11][12][13] employment of lightweight yet high strength composite materials, [14][15][16] hybrid materials, and structures, 17 and application of advanced manufacturing processes, 18,19 etc. Considering existing car bodies are generally composed of a large number of hollow TWSs possessing specific demands or constraints on strength, lightweight, spatial layout, or cost, aluminum foam-filling is considered a relatively more efficient and practical approach to improve the crashworthiness of TWSs of a car body due to its excellent compressive energyabsorption characteristics and light weight. 202][23][24][25][26][27][28] For instance, Zhang et al. 29 investigated dynamic response and energy absorption performance of aluminum foam-filled sandwich circular tubes under internal blast loading.Wang et al. 30 studied the energy absorption behavior of an aluminum foam-filled circular-triangular nested tube energy absorber under impact loading; Sun et al. 31 conducted experimental and numerical investigation into the crashworthiness of aluminum-foam-composite hybrid structures; Qi et al. 32 executed comparative study on empty and foam-filled hybrid material double-hat beams under lateral impact.Xu et al. 33 explored mechanical properties of aluminum foam filled re-entrant honeycomb with uniform and gradient designs, etc. Besides, there also have been a few studies on crashworthiness characteristics of foam-filled thinwalled crash boxes.Santosa and Wierzbicki 34 studied the effect of low density filler material, such as aluminum honeycomb, or foam, on the crushing resistance of a square box column under axial quasi-static loading.Toksoy and Gu¨den 35 investigated the crushing behavior of partially Al closed-cell foam filled commercial 1050H14 Al crash boxes under quasi-static and dynamic deformation velocities.Wang et al. 4 carried out crashworthiness study on a novel aluminum foamfilled crash box with partial filling and combined trigger design under axial impact loading, etc.
In addition to the proposed effective methods mentioned above, the optimization techniques, including experimental design, 36,37 surrogate modelling, 38,39 and optimization algorithms, 40,41 etc., which are extensively employed to handle all kinds of optimization issues, [42][43][44][45][46][47][48] have also been widely employed to further enhance the crashworthiness of TWSs or thin-walled crash boxes.For instance, Ebrahimi et al. 40 executed crashworthiness efficiency optimization for two-directional functionally graded foam-filled tubes under axial crushing impacts using multi-objective particle swarm optimizations (MOPSOs).Meric x and Gedikli 49 performed multiobjective optimization of energy absorbing behavior of foam-filled hybrid composite tubes based on feedforward artificial neural networks (FFNN).Sun et al. 41 conducted multi-objective optimization for thin-walled structures with functionally graded thickness by combining response surface method (RSM) and Non-dominated sorting genetic algorithm-II (NSGA-II).Yin et al. 24 performed multi-objective robust optimization of foam-filled bionic thin-walled structures through a hybrid method of ensemble meta-model, NSGA-II, 3-sigma robust design, and Monte Carlo simulation (MCS).Zhou et al. 37 carried out multi-objective design optimization of a novel NPR crash box based on optimal Latin square design method, RSM and NSGA-II.Hou et al. 50carried out a multi-objective optimization of a novel crash box composed of Al shell and 3Dprinted lattice core for maximum specific energy absorption (SEA) and minimum mass using optimal Latin hypercube sampling (OLHS), RSM and NSGA-II.Tan et al. 36 conducted a multi-objective optimization of a novel auxetic hierarchical honeycomb crash box by combining RSM with NSGA-II and archive-based micro genetic algorithm (AMGA).Wang et al. 51 conducted a many-objective optimization design of a novel hexahedral pyramid (H-P) crash box composed of H-P negative Poisson's ratio (auxetic) inner core structure and outer shell layer by integrating the RSM method and the multi-objective evolutionary algorithm based on decomposition with detecting and escaping strategy (MOEA/D-DAE).Xie et al. 38 carried out a multiobjective crashworthiness optimization of energyabsorbing box with gradient lattice structure by combining OLHS, RSM-RBF neural network (RBFNN)-Kriging hybrid meta-models, NSGA-II and entropygrey relational analysis (E-GRA), etc.
Through extensive literature review and research, it has been found that the crashworthiness of TWSs or thin-walled crash boxes with conventional crosssectional shapes (circular, triangular, square, hexagonal, polygon, elliptical, hat-shaped, etc.) and different modes (hollow, filled) under uni-directional (axial, transverse, etc.) crushing loading have been extensively studied, while few studies have been reported concerning crashworthiness of novel double-hexagonal thinwalled crash boxes under multi-angle impact loading, which is more realistic considering collision angles may vary in actual car collisions and thus affect the energy absorption characteristics.3][54][55][56] Engineering designs of crash box obtained from such deterministic optimization would have a risk of violating performance requirements and place occupants at potential risks of injuries during a car collision. 3In addition, conventional optimization approaches based on experimental design, surrogate modeling (RSM, RBF, Kriging, etc.) and optimization algorithms (NSGA-II, MOPSO, MOEA, etc.) are extensively employed to handle such engineering optimization issues, while they generally rely on a large number of calculating samples to obtain sufficiently accurate surrogate models, which usually leads to a significant increase in computational costs.
Hence in the present study, a novel thin-walled double-hexagonal crash box is first proposed (as shown in Figure 1) and then multi-objective robust optimized for better overall crashworthiness under multi-angle impact loading, using a proposed hybrid method combining aluminum foam-filling and Taguchi-grey relational analysis (GRA).Considering the proposed crash box is irregularly-shaped longitudinally with the upper end face inclined to match the bow-shaped bumper, a substructure with regular double-hexagonal cross section and flat upper and lower end faces at the lower end of the crash box, herein named as hollow doublehexagonal column (H-DHC), is extracted and taken as object for aluminum foam-filling, crashworthiness investigation, and multi-objective robust optimization.Specifically, the finite element (FE) models of the DHCs under multi-angle impact loading, including hollow (H-DHC) and foam-filled (F-DHC), are first built and validated by experiments.On this basis, a comprehensive crashworthiness comparison is conducted to explore relative merits of F-DHC over original H-DHC under multi-angle impact loading.After that, the F-DHC is multi-objective robust optimized for maximizing overall specific energy absorption (SEA u ) and minimizing overall initial peak crushing force (IPCF 0 ) simultaneously under multiple-angle impact loading, using a hybrid method of Taguchi-GRA.At last, a bumper-crash box integrated crashworthiness analysis under multi-angle impact loading is executed to further verify the optimization.Figure 2 outlines the flowchart of the present study.
The remaining part of the present study is organized in this manner: section 2 introduces the finite element (FE) modeling and verification of DHCs (H-DHC, F-DHC), section 3 elaborates the crashworthiness indicators and comparison of H-DHC and F-DHC under multi-angle impact loading, followed by section 4 with the multi-objective robust optimization of F-DHC using Taguchi-GRA, and section 5 concludes the present study at last.

Geometrical configuration
In the present study, a novel thin-walled crash box with double-hexagonal cross section is first proposed, whose upper end face is inclined designed to match the bowshaped bumper, as shown in Figures 1 and 3.For convenience of research and without losing generality, a substructure with regular double-hexagonal cross section and flat upper and lower end faces at the lower end of the crash box, herein named as hollow doublehexagonal column (H-DHC), is extracted and taken as object for aluminum foam-filling, crashworthiness investigation, and multi-objective robust optimization, as shown in Figure 3.The H-DHC is composed of inner plate and outer plate.To improve energy absorption without adding too much weight, the space between inner plate and outer plate of H-DHC is fully filled with lightweight aluminum foam to constitute the foamfilled double-hexagonal column (F-DHC).The geometric parameters of the initial DHCs (H-DHC and F-DHC) are as follows: column length L = 140 mm, hexagonal side length a = 50 mm, inner plate thickness t 1 = 2 mm and outer plate thickness t 2 = 2.2 mm.To monitor the load and boundary conditions of a thinwalled crash box in real vehicle crash scenarios, the bottom ends of the DHCs are placed on a rigid base fully fixed, while the top ends are subjected to a rigid wall at an initial velocity of 15 m/s.Impact angles ranging from 0°to 30°at an interval of 5°(namely 0°, 5°, 10°, 15°, 20°, 25°, and 30°) are considered to take into account the uncertainty of impact direction.

Material Behavior
The column walls of DHCs, namely the inner and outer plates, are both made of aluminum alloy with material parameters as follows: density r = 0.27 g/cm 3 , Young's  modulus E = 68.2GPa and Poisson's ratio v = 0.3.The stress-strain curve of the aluminum alloy is shown in Figure 4, 57 and the elasto-plastic model with modified piecewise linear plasticity is employed to constitute the material behaviors.In contrast, A principal structure model proposed by Deshpande and Fleck 58 is used to constitute the behavior of the aluminum foam, whose yield criterion could be described as: where s y denotes the yield stress and s e represents the equivalent stress.The s e could be further described as: where s v represents the von Mises effective stress and s m denotes the mean stress.Besides, a is a parameter about yield surface shape control, which could be expressed by Poisson's ratio v p with a function relationship.
where the Poisson's ratio v p equals 0 for the aluminum foam. 59According to equation (3), a ' 2.12 can be obtained.The strain hardening rule could be summarized as: where e e represents the equivalent effect strain.s p , a 2 , g, b, E P , and e d denote the foam material parameters, which are related to the foam density by the formula as follows: where r f 0 and r f represents density of the base material and foam, respectively.[62] Finite element (FE) modeling In the present study, the explicit non-linear FE code LS-DYNA is employed to perform numerical investigation on the crashworthiness characteristics of the DHCs (H-DHC, F-DHC) under multi-angle impact loading.
Consistent with the geometric model, the bottom ends   of the DHCs are placed on a rigid base fully fixed, while the top ends are subjected to a rigid wall at an initial velocity of 15 m/s from a set of various angles (0°, 5°, 10°, 15°, 20°, 25°, and 30°).The inner plate and outer plate of DHCs are modeled using Belytschko-Tsay quadrilateral shell elements with five integration points across the thickness and one in the element plane.In contrast, the aluminum foam is modeled using eightnode hexahedron elements with one-point reduced integration mechanism coupled with stiffness based hourglass control. 63The inner plate, outer plate and the foam filler are all meshed with an element size of 4 mm to reduce the computational cost while maintaining simulation accuracy. 59To avoid interpenetration, the ''automatic single surface'' contact is defined for both the inter plate and the outer plate, the ''automatic surface to surface'' contact is defined among the inter plate, out plate and the foam filler, and the ''automatic node to surface'' contact is defined among the upper rigid wall, middle DHCs, and the lower rigid base.For all the contacts defined above, the static and dynamic coefficients are uniformly set to 0.3 and 0.2, respectively. 63Figure 5 demonstrates the built FE models of DHCs, including H-DHC and F-DHC, under multiangle impact loading.

Verification of FE models
To verify the accuracy of the developed FE model of H-DHC, the original hollow and irregularly-shaped crash box is manufactured and selected to perform experimental verification, as shown in Figure 6.Herein the material of the crash box is aluminum alloy AA6063, whose stress-strain curve is also shown in Figure 6.The bottom end of the crash box is placed on a rigid base fully fixed, while the top end is compressed by the machine platform at a velocity of 10 mm/min.
To avoid test error and for better comparison, 6 specimens of crash box with identical structural parameters are tested.Figure 7 compares the crushing force versus displacement curves for the FE simulation and experimental tests, from which it can be seen that the six groups of test results are basically the same, and for each specimen, the simulation result can well captures the deformation progress in the test.In addition, the corresponding final deformation modes of the six specimens and that of a simulation with six views are compared as shown in Figure 8. From which, it is easily seen that the final deformation modes of experimental results agree fairly well with that of simulation.Without losing generality, the foam-filled square column (F-SC) is employed to indirectly verify the accuracy of the developed FE model of F-DHC, considering there are few relevant reports concerning crushing experiment of foam-filled double-hexagonal column as far as the authors know.Figure 9 compares the crushing force versus displacement curve and the final deformation mode between the FE simulation and the experiment performed by Li et al., 64 while Table 2 further compares the extracted initial peak crushing force (IPCF) and mean crushing force (MCF) on the whole stroke bet-ween the experiment and the simulation.From which it also can be seen that there is reasonable agreement between the simulation and experimental results.
Considering the effect of strain rate of the aluminum alloys employed in the present study are insensitive and can be ignored, the FE models verified under the quasistatic crushing can be extended to dynamic scenarios.In a summary, the developed FE models of DHCs, including H-DHC and F-DHC, are of acceptable accuracy and provide considerable confidence for us to explore the crashworthiness characteristics and optimization approach under multi-angle impact loading.

Crashworthiness indicators
To quantitatively evaluate the crashworthiness of the DHCs under dynamic impact loading, it is essential to predefine crashworthiness indicators.There are many crashworthiness indicators available, among which Energy absorption (EA), specific energy absorption (SEA), peak crushing force (PCF), mean crushing force (MCF), and crushing force efficiency (CFE) have been widely utilized to estimate the crashworthiness of thinwalled structures. 60Specifically, EA denotes the absorbed energy of structures via their plastic deformation, which demonstrates the total capacity for energy absorption.SEA, which denotes the energy absorption per unit mass, is often considered a more objective and effective indicator to measure the energy absorption efficiency of material.PCF and MCF respectively  denote the maximum and average value of instantaneous crushing force on the whole stroke, while CFE demonstrates the ratio of MCF relative to PCF, which measures the load consistency and uniformity.The above crashworthiness indicators, as shown in Figure 10, are normally calculated as follows: Nevertheless, under multi-angle impact loading scenarios, the varied impact angles could largely affect the value of energy absorption.In this regard, an ideal energy absorber is expected to demonstrate superior overall energy absorption characteristics under multiangle impact loading.Thus, taking into account the influence of different impact angles, a novel overall crashworthiness indicator SEA u is employed and defined as follows: where SEA u k represents the SEA under the specific impact angle u k (k = 1, 2, ., N u ), N u denotes the total number of impact angles considered.v uk stands for the   relative weight (importance or occurrence probability) of u k .and can be described as equation (12).In addition, in the automotive engineering, the initial peak crushing force (IPCF), which denotes the first peak of the crushing force and normally measures the load severity when initial plastic folding happens, is often considered critical to survival rate of occupants when vehicle collision occurs.More specifically, high IPCF may lead to severe injury or even death of occupants.Further, under multi-angle impact loading scenarios, the IPCF under impact angle of 0°(marked as IPCF 0 ) is generally the maximum among all impact angles, which should be primarily controlled.Hence, for DHCs under multi-angle impact loading in the present study, IPCF 0 is adopted as another key indicator to evaluate overall crashworthiness.Thus, SEA u and IPCF 0 are finally adopted to evaluate overall crashworthiness under multi-angle impact loading.Obviously, for DHCs in the present study, higher values of SEA or SEA u indicating a higher capacity for energy absorption are preferred.In contrast, PCF, PCF 0 , IPCF, or IPCF 0 are all ''smaller-the-better'' to reduce the load severity and risk of occupant injuries.

Crashworthiness comparison
To explore the potential relative merits of F-DHC compared with H-DHC under multi-angle impact loading, a comprehensive crashworthiness comparison is first conducted.Except for the aluminum foam filler, other structural parameters of H-DHC, and F-DHC are kept the same.Figure 11 compares the crashworthiness indicators, including EA, SEA, IPCF, and MCF, between H-DHC and F-DHC under impact angles ranging from 0°to 30°with an interval of 5°, while Table 3 further compares the corresponding final deformation modes.Obviously, for both H-DHC and F-DHC, the values of all above crashworthiness indicators decrease gradually along with increasing impact angles, which well confirms the importance and necessity of considering multi-angle impact loading scenarios.Besides, it can be seen that the EA and MCF of F-DHC under each impact angle outweigh that of H-DHC due to the introduction of aluminum foam, and the folding lobes of F-DHC is relatively more regular than that of H-DHC especially under larger impact angles, which both reveal the effect of aluminum foam filler on improving total energy absorption capacity.Nevertheless, it can also be seen that the SEA of F-DHC under each impact angle is smaller than that of H-DHC, which discloses the fact that the aluminum foam filler has no relative advantages over column walls, namely inner plate and outer plate, from the point of material energy absorption efficiency.Besides, note that the IPCF of F-DHC under each impact angle is relatively larger than that of H-DHC, which is due to the enhancement of stiffness caused by introduction of aluminum foam filler.To sum up, the F-DHC outperforms the H-DHC concerning energy absorption, while it is inferior to the latter concerning load severity which is closely related to occupant injuries in a collision.Considering the above two aspects normally are competing and even conflicting, hence multi-objective optimization approach for F-DHC needs further investigation to obtain better overall crashworthiness under multi-angle impact loading.

Multi-objective robust optimization of F-DHC
For F-DHC, as mentioned earlier, direct employment of deterministic optimization results would encounter large potential risks in actual vehicle collision scenarios.Hence, multi-objective robust optimization for F-DHC is more preferable and adopted in the present study.Compared with other robust optimization methods, the Taguchi robust method, which adopts Taguchi orthogonal array (OA) including design and noise factors and handles the experimental results based on signal-to-noise ratio (SNR) that to be maximized, demonstrates notable merits on evaluating robustness of experiments, influence of factors on response and exploring optimal combination of factor levels with high efficiency and dramatic simplicity.Nevertheless, the conventional Taguchi robust method is generally limited to optimizing a single response at a time, which is not suitable to handle the multi-objective robust optimization of F-DHC herein.However, this could be well overcome by coupling Taguchi robust method with gray relational analysis (GRA), namely Taguchi-GRA.Through GRA, which is acknowledged as a relatively accurate method for multiple criteria decision making (MCDM), multiple responses could be converted to a single response called gray relational grade (GRG), which stands for the relative closeness of a candidate to the ideal referential alternative and normally is ''larger-the-better'' (LTB).Therefore, Taguchi-GRA is employed herein to perform multi-objective robust optimization of F-DHC for better overall crashworthiness under multi-angle impact loading.

Taguchi-GRA methodology
The main procedure of Taguchi-GRA are summarized as follows: Step 1: Taguchi experimental design.In general, the main steps of Taguchi experimental design include: (1) Define the design variables and corresponding levels, including control factors and noise factors.(2) Define the responses (or objectives) and corresponding quality characteristics.(3) Build the Taguchi orthogonal array (OA), with control factors constituting the inner table and noise factors constituting the outer Table 4 Execute the experimental samples one by one and count the response data.
For multi-objective robust optimization of F-DHC in the present study, structural parameters including foam density (r), inner plate thickness (t 1 ) and outer plate thickness (t 2 ) are selected as three control factors and for each factor five levels are assigned, thus a L 25 (53) orthogonal array is built to constitute the inner table.Meanwhile, considering the inherent uncertainty in the material properties of the plates of F-DHC, the density (D), elastic modulus (E), Poisson's ratio (P), and yield strength (Y) of plate material are chosen as four noise factors, and for each factor three levels are

Impact angle 0°5°10°15°20°25°30°H
-DHC F-DHC assigned, thus a L 9 (34) orthogonal array is built to constitute the outer table.Considering multi-angle impact loading scenarios, the overall SEA (SEA u ) is taken as one response (or objective), while IPCF under 0°impact loading (IPCF 0 ) is taken as an another based on the fact that IPCF of F-DHC obtains maximum under pure axial impact loading, as disclosed by Figure 11.Table 4 demonstrates design variables and corresponding levels, while Table 5 demonstrates the constituted OA table for Taguchi experimental design, which is a specially designed array that permits the investigation of main effects and interaction effects of design variables on responses.
Step 2: Single-to-noise ratio (SNR).After the Taguchi experimental design has been completed and the response data for each factor level has been counted, the signal-to-noise ratio (SNR), which measures the magnitude of the mean of a process compared to its variation, is employed to evaluate the response performance according to the corresponding characteristics. 67pecifically, the SNR of a response with ''larger-thebetter (LTB)'' characteristic could be calculated using equation (13), while the SNR of a response with ''smaller-the-better (LTB)'' characteristic could be calculated using equation ( 14): where SNR ij represents the SNR for the jth response of ith experiment, y ijk denotes the experimental result for jth response of ith experiment at kth measurement, n denotes the total number of measurements per experiment.
For multi-objective robust optimization of F-DHC in the present study, two conflicting objectives namely SEA u and IPCF 0 , are selected to measure and robustly optimize the overall crashworthiness of F-DHC under multi-angle impact loading.Obviously, SEA u is largerthe-better (LTB) to improve overall energy absorption, while IPCF 0 is smaller-the-better (STB) to control load severity.Hence, equations ( 13) and ( 14) are respectively adopted for SEA u and IPCF 0 to calculate corresponding SNR.
Step 3: Gray relational analysis (GRA).In general, for a single-response (or objective) optimization, based on the obtained SNR of response, the Taguchi method is capable of evaluating the influence of each factor on the response and further determining the optimal combination of factor levels via main effect analysis.However, as aforementioned earlier, the traditional Taguchi method is limited to handling single-objective optimization, which therefore is not suitable to address multi-objective robust optimization of F-DHC herein.To eliminate this obstacle, GRA is incorporated to Taguchi method to convert the multiple responses to a single gray relational grade (GRG) that is to be maximized.Specific procedures of GRA are as follows: Gray relational generation.First, the experimental data sequence (from Taguchi OA) are transformed to a uniform range of [0,1] through different normalization methods according to response characteristics, which is generally named as gray relational generation.Specifically, for original sequence with ''larger-the-better'' (LTB) response 68,69 : Else, for original sequence with ''smaller-the-better'' (STB) response: Otherwise, for original sequence with ''nominal-thebest'' (NTB) response: where x o i (k) and x Ã i (k) denote the original and normalized value of kth response at ith experiment, respectively.max x o i (k) and min x o i (k) denote the maximum and minimum value of x o i (k) for the kth response among i experiments, respectively.x o i is a target value, n and m denote the total number of experiments and responses, respectively.
Gray relational coefficient (GRC).Second, the GRC j i (k) is calculated from the normalized experimental data to measure the correlation between the current experiment and the desired solution as: where x Ã r (k) denotes the ideal normalized value, that is, the maximum normalized SNR (since SNR is LTB), for the kth response among all experiments,D ri (k) denotes the difference between x Ã r (k) and x Ã i (k), d is the identification coefficient, d 2 (0,1) and d = 0.5 is generally recommended and adopted in this paper. 68,69ay relational grade (GRG).Finally, the gray relational grade (GRG) which measures the correlation of current experiment and the ideal solution can be calculated by averaging the GRC for multiple responses as: Nevertheless, taking into account the difference of response preference or significance, different weights are assigned to different responses to calculate the weighted GRG as follows: where GRG i and GRG iv denote the GRG and weighted GRG for ith experiment, v k denotes the weight of kth response, n denotes the total number of responses.For multi-objective robust optimization of F-DHC in the present study, n = 2 as SEA u and IPCF 0 are selected as two responses.The GRG could be treated as an overall response of the crashworthiness characteristics of F-DHC instead of multiple responses -SEA u and IPCF 0 .Obviously, larger value of GRG i or GRG iv signifies that the corresponding combination of factor levels is closer to the optimal condition.
Step 4: Main effect and interaction effect analyses.Through GRA, the multi-objective robust optimization has now been converted to a single-objective robust optimization of GRG, which then could be handled by Taguchi analysis to explore the optimal combination of factor levels.Specifically, the average response value of GRG of each factor level is first calculated, then the main effect value of each factor, which is the difference between the maximum and minimum average GRG of factor levels, could be obtained.The process could be described as follows: where f u is the mean GRG value for the uth level of factor f, N is the number of uth level experiments of factor f, and G uv is the GRG value of vth experiment with the uth level of factor f, M represents the total number of levels of factor f while Df is the main effects.
In general, the average response value of each factor level, the main effect value of each factor on response and the interaction effect value of any two factors on response could be obtained and demonstrated in a factor response table, a main effect plot and an interaction effect plot, respectively.Based on which, the main effect of any one factor and the interaction effect of any two factors on the response could be figured out.More significantly, the optimal combination of factor levels for optimal response could be determined.Normally, the lower the main effect value of a factor, the smaller the influence of the factor on the response, and vice versa.In contrast, the interaction plots measure the influence of factor interaction on response.Parallel lines in an interaction plot indicate no interaction, and the greater the difference in slope between the lines, the higher the degree of interaction.

Results and discussion
The experimental results based on the Taguchi OA are demonstrated in Table 5, based on which the SNR of SEA u (marked as SNR1) and IPCF 0 (marked as SNR2) for each experiment are respectively calculated according to equations ( 13) and ( 14), and the corresponding results are counted and displayed in Table 6.
To explore the influence of control factors (r, t 1 and t 2 ) on responses, the response tables for SNR1 of SEA u and SNR2 of IPCF 0 are calculated and demonstrated in Tables 7 and 8, respectively.For a certain factor, the mean value at a certain level represents the global average value of SNR for all experiments with this specific level, and the Delta denotes the maximum difference of SNR mean among different levels.Generally, a factor with a higher Delta among the factors indicates a relatively larger influence on the response compared with others.Therefore, the most influential factor on SEA u is outer plate thickness (t 2 ), followed by inner plate thickness (t 1 ), and finally foam density (r) according to Table 7.In contrast, the most influential factor on IPCF 0 is also the outer plate thickness (t 2 ), followed by foam density (r), and finally inner plate thickness (t 1 ) according to Table 8.Therefore, both SEA u and IPCF 0 receive the largest influence from outer plate thickness (t 2 ).Besides, the main effect plots of control factors on the SNR1 of SEA u and SNR2 of IPCF 0 are demonstrated in Figure 12, from which the influence of each factor on the SNR of responses can also be figured out by comparing the slopes of the lines.Nevertheless, a more significant function of main effect plots is to determine the optimal combination of factor levels for each response.To be more specific, the optimal combination of factor levels for maximizing SNR1 (namely maximizing SEA u ) is F5-I5-O5 (namely r = 0.5 g/cm 3 , t 1 = 2.8 mm and t 2 = 3.0 mm.), as boldly marked in Table 7.In contrast, the optimal combination of factor levels for maximizing SNR2 (namely minimizing IPCF 0 ) is F1-I1-O1 (namely r = 0.1 g/cm 3 , t 1 = 1.2 mm, and t 2 = 1.4 mm), as boldly marked in Table 8.
In addition, to explore the influence of interactions among control factors on a single response, the interaction effect plots of control factors on SNR1 of SEA u and SNR2 of IPCF 0 are demonstrated in Figures 13  and 14, respectively.An interaction plot is a plot of means for each level of a factor with the level of a second factor kept constant, which could be employed to evaluate the two-way interactions among factors on the responses.To be more specific, there is no interaction if the lines are parallel, otherwise there is, and the greater the lines deviate from being parallel, the larger the interaction.It can be seen that for SNR1 (Figure 13), the interaction effect is observed between either two of the three factors (foam density, inner plate thickness, and outer plate thickness), among which the interaction between inner plate thickness and outer plate thickness appears relatively larger, indicating the relationship between inner plate thickness, and SEA u , depends on outer plate thickness and vice versa.Likewise, for SNR2 (Figure 14), the interaction effect is also observed between either two of the three factors, while the corresponding interactions are restively smaller than that of SNR1.
As just revealed, the optimal combination of factor levels for maximizing SEA u (F5-I5-O5) are entirely different with that for minimizing IPCF 0 (F1-I1-O1).Actually, the optimal combination of factor levels that make one goal the best will make the other goal the worst.However, the automobile crash box prefers the F-DHC be manufactured at the same level of design to obtain optimal SEA u and IPCF 0 simultaneously, which    is not suitable to be handled with traditional Taguchi method aimed to optimize a single response at a time.Therefore, to achieve multi-objective robust optimization of F-DHC, the two conflicting responses (SEA u and IPCF 0 ) are converted into a single response, that is, gray relational grade (GRG) that to be maximized (LTB), through gray relational analysis (GRA).For this purpose, the gray relational generation based on equation (15) for SEA u and equation ( 16) for IPCF 0 , the gray relational coefficient based on equation (18)  and the gray relational grade (GRG) based on equation (19) or equation ( 20) are obtained and demonstrated in Table 9.Note that herein to explore the effect of different weights on the optimization results, SNR1 is given three different weights (v 1 = 0.4,0.5,0.6) and the weights of SNR2 (v 2 = 0.6,0.5,0.4) are changed accordingly, resulting in three different groups of GRG as shown in Table 9.
The response table for GRG under three different weights are calculated and demonstrated in Tables 10 to 12, respectively.Obviously, the influence of control factors (r, t 1 and t 2 ) on GRG differs along with different weights.Specifically, the most influential factor on GRG when v 1 = 0.4 is foam density (r), followed by outer plate thickness (t 2 ), and finally the inner plate thickness (t 1 ).In contrast, the most influential factor on GRG when v 1 = 0.5 is foam density (r), followed by inner plate thickness (t 1 ), and finally outer plate thickness (t 2 ), which is identical to the results when v 1 = 0.6.Considering the main effect and interaction effect plots for GRG under such different weights are of little difference, the main effect and interaction effect plots for GRG when v 1 = 0.4 are selected as a representative to be shown in Figures 15 and 16, respectively.According to Figure 15, the optimal combination of factor levels for GRG when v 1 = 0.4 is F1-I5-O1 (namely r = 0.1 g/cm 3 , t 1 = 2.8 mm, and t 2 = 1.6 mm), as boldly marked in Table 10.In contrast, the optimal combination of factor levels for both GRG when v 1 = 0.5 and GRG when v 1 = 0.6 is F1-I5-O5, as boldly marked in Tables 11 and 12, respectively.In addition, according to Figure 16, relatively larger interaction effect on GRG is observed among either two of the three factors compared with that on SEA u or IPCF 0 solely just discussed, indicating that the three factors demonstrate large interaction effects on overall crashworthiness (SEA u and IPCF 0 ) under multi-angle impact loading.

Verification of optimal F-DHC
Now that the optimal combination of factor levels, as well as corresponding structural parameters, for better overall crashworthiness (SEA u and IPCF 0 ) of F-DHC is obtained, the subsequent step is to verify the effectiveness of the optimization results.Considering the optimal combination of factor levels under v 1 = 0.5 and v 1 = 0.6 is identical, hence two case are designed to compare the verification: (1) Case 1 (v 1 = 0.4): F1-I5-O1, namely r = 0.1 g/ cm 3 , t 1 = 2.8 mm and t 2 = 1.6 mm; (2) Case 2 (v 1 = 0.5 or 0.6): F1-I5-O5, namely r = 0.1 g/cm 3 , t 1 = 2.8 mm and t 2 = 3.0 mm.
A detailed crashworthiness comparison among the optimal design in case 1, the optimal design in case 2 and the initial design (F3-I3-O3, namely r = 0.3 g/cm 3 , t 1 = 2.0 mm and t 2 = 2.2 mm) of F-DHC under multiangle impact loading is conducted.Table 13 presents the comparison of overall crashworthiness indicators of F-DHC under multi-angle impact loading.It can be seen that: (1) In case 1, the optimal design obtains 24.34% higher of SEA u and 27.08% lower of IPCF 0 ; (2) In case 2, the optimal design 49.59% higher of SEA u and 27.61% higher of IPCF 0 .Obviously, compared with that in case 1, the optimal design in case 2 gains more increase of SEA u (49.59%vs 24.34%) due to larger weight (v 1 = 0.5 or 0.6 vs v 1 = 0.4), while it obtains less decrease of IPCF 0 (227.61%vs 27.08%) due to smaller weight (v 2 = 0.5 or 0.4 vs v 2 = 0.6), indicating the relative weight directly affects the optimization result.More specifically, for a certain response (or objective), the larger the weight, the better the optimization for it.Nevertheless, considering the application of F-DHC in actual crash box, an ideal optimal design should not pursue maximum SEA u at the expense of larger IPCF 0 that represents worse load severity, therefore the optimal design in case 1 is considered a relatively better choice for optimization of F-DHC concerning overall crashworthiness.
Besides, for better comparison, the crashworthiness indicators (SEA, IPCF) of optimal F-DHC in case 1 (abbreviated as Case 1 in the table) and case 2 (abbreviated as Case 2 in the table) under single-angle impact loading (0°, 5°,10°, 15°, 20°, 25°, and 30°) are also compared with that of the initial design, as shown in     decrease, the IPCF under most impact angles (0°, 5°,10°, 15°, 20°, and 25°) are dramatically increased.Like discussed earlier, the optimal F-DHC in case 2 gains large increase of SEA at the expense of elevating IPCF dramatically under almost all impact angles, which is not an ideal solution in actual application on crash box.In contrast, the optimal F-DHC in case 1 gains moderate increase of SEA under small impact angles and slight decrease of SEA under large impact angles, meanwhile reducing IPCF significantly under all impact angles, which is considered a more ideal solution compared to that of case 2, considering actual vehicle crashes are more likely to take place under small impact angles, and during which the load severity (IPCF) should be controlled of more priority to guarantee occupant safety.In addition, the corresponding final deformation modes of the optimal design in case 1 and the initial design under single-angle impact loading are also compared as shown in Table 15, in which only subtle difference are observed, indicating the overall crashworthiness of F-DHC under multi-angle impact loading have been improved under the condition of slight changes in deformation modes.To sum up, the multi-objective robust optimization of F-DHC is reasonable and effective.

Verification of optimization for crash box
As explained at the beginning of this study, the regularly-shaped H-DHC is extracted from an irregularly-shaped automobile crash box to perform foam-filling and multi-objective robust optimization, whether the obtained optimal F-DHC improves the crashworthiness of the crash box or not needs further investigation.To this end, the FE model of an automobile body with the crash box studied under multi-angle impact loading is built, which is further simplified to a bumper-crash box integrated model under multi-angle impact loading as shown in Figure 19.Specifically, the below rigid bases of the crash box are fully constrained while the bumper is impacted by a rigid wall at an initial velocity of 15 m/s from multiple angles.On this  basis, the crashworthiness of the bumper-crash box assembly with the optimal F-DHC (abbreviated as optimized design) under multi-angle impact loading are compared with that of the initial design.Note that for bumper-crash box assembly, the output forcedisplacement curve is complicated by the influence of the bumper, therefore IPCF is no longer a suitable indicator to evaluate load severity of the while assembly.In contrast, PCF, which indicates the global maximum of crushing force on the whole stroke, is relatively more suitable and thus adopted instead.Likewise, PCF 0 , which denotes the PCF under 0°impact angle and generally the maximum among all angles, is employed to evaluate the load severity under multi-angle impact loading.
Consequently, the crashworthiness of the optimized design and the initial design are compared, among which Table 16 presents the overall crashworthiness indictors of bumper-crash box assembly under multi-angle impact loading, Table 17 and Figure 20 demonstrates the crashworthiness indictors under individual single-angle impact loading and Figure 21 shows the corresponding final deformation modes.According to Table 16, the optimized design gains 18.31% higher of SEA u and 20.76% lower of PCF 0 simultaneously compared with the initial design.In addition, compared with the initial design, the optimized design gains different percent of increase of SEA and obtains different percent of decrease of PCF under all single-angle impact loading simultaneously, as shown in Table 17 and Figure 20.Meanwhile, no evident change is observed in final deformation modes before and after optimization, as shown in Figure 21, which is desired especially for component-level optimization.Hence, the crash box is effectively multi-objective optimized with the optimal F-DHC.To sum up, the multi-objective robust optimization of the crash box is reasonable and effective.

Conclusion
The present study first proposes a novel doublehexagonal crash box, and then puts forward a multiobjective robust optimization method combining aluminum foam-filling and Taguchi-gray relational analysis (GRA) to optimize the crash box for better overall crashworthiness under multi-angle impact loading.The main findings, limitations, and recommendations of this study are as follows: (1) The constructed finite element (FE) models of DHCs, including hollow (H-DHC) and foamfilled (F-DHC), are of good accuracy with simulation results basically consistent with that of experiments.Therefore, the subsequent crashworthiness comparison and multiobjective robust optimization based on FE models are reasonable and credible.(2) The foam-filled DHC (F-DHC) outperforms the hollow DHC (H-DHC) concerning total energy absorption, while it is inferior to the latter concerning load severity.Hence, multi-objective optimization for F-DHC is desired to seek optimal design for better overall crashworthiness under multi-angle impact loading.(3) According to the main effect analysis, the sorting of factor influence differs for single response alone (SEA u , IPCF 0 ) and multiresponse together (GRG), which is affected by the weights assigned for each response.In addition, obvious interaction effects are observed among foam density (r), inner plate thickness (t 1 ) and outer plate thickness (t 2 ) for SEA u , IPCF 0 , and GRG.(4) Compared with the initial design, the optimal F-DHC obtained through Taguchi-GRA gains 24.34% higher of SEA u and 27.08% lower of IPCF 0 simultaneously, which further improves the bumper-crash box assembly with 18.31% higher of SEA u and 20.76% lower of PCF 0 simultaneously, indicating the multi-objective robust optimization for the crash box is reasonable and effective, thus aluminum foam-filling combined with Taguchi-GRA could be an effective approach for multi-objective robust

Figure 1 .
Figure 1.Novel double-hexagonal crash box under multi-angle impact loading.

Figure 2 .
Figure 2. The flowchart of the present study.

Figure 5 .
Figure 5. FE models of DHCs under multi-angle impact loading.

Figure 8 .
Figure 8.Comparison of the final deformation modes: (a) Experiments and (b) simulation.

Figure 9 .
Figure 9.Comparison of the experiment 64 and simulation for F-SC: (a) Force-displacement curve and (b) Final deformation modes.

Figure 12 .
Figure 12.Main effect plots of factors on SNR1 of SEA u and SNR2 of IPCF 0 : (a) Mean of SNR1 and (b) Mean of SNR2.

Figure 13 .
Figure 13.Interaction effect plots of factors on SNR1.

Figure 14 .
Figure 14.Interaction effect plots of factors on SNR2.

Figure 15 .
Figure 15.Main effect plot of factors for GRG when v 1 = 0.4.

Figure 16 .
Figure 16.Interaction effect plots of factors for GRG when v 1 = 0.4.

Figure 19 .
Figure 19.FE model of automobile body and bumper-crash box assembly under impact loading.

Figure 21 .
Figure 21.Comparison of final deformation modes under single-angle impact loading.

Table 3 .
Final deformation modes of H-DHC and F-DHC.

Table 4 .
Control factors, noise factors, and corresponding levels.

Table 5 .
Taguchi experimental layout and results.

Table 6 .
The calculation results of SNR.

Table 7 .
Response table for SNR1 of SEA u .

Table 14 .
Comparison of crashworthiness indicators of F-DHC under single-angle impact loading.

Table 16 .
Comparison of overall crashworthiness indicators of bumper-crash box assembly under multi-angle impact loading.

Table 15 .
Comparison of final deformation modes under single-angle impact loading.

Table 17 .
Comparison of crashworthiness indicators under single-angle impact loading.