A Study of Multi-Objective Crashworthiness Optimization of the Thin-Walled Composite Tube under Axial Load

In recent decades, thin-walled composite components have been widely used in the automotive industry due to their high specific energy absorption. A large number of experimental and numerical studies have been conducted to characterize the energy absorption mechanism and failure criteria for different composite tubes. Their results indicate that the energy absorption characteristics depend highly on the failure modes that occur during the impact. And failure mechanism is dependent on fiber material, matrix material, fiber angle, the layout of the fibers, as well as the geometry of structure and load condition. In this paper, first, the finite element (FE) model of the CFRP tube was developed using the Tsai-Wu failure criterion to model the crush characteristics. The FE results were validated using the published experimental. Then, a series of FE simulations were conducted considering different fiber directions and the number of layers to generate enough data for constructing the GMDH-type neural network. The polynomial expression of the three outputs (energy absorption, maximum force, and critical buckling force) was extracted using the GMDH algorithm and was used to perform the Pareto-based multi-objective optimizations. Finally, the failure mechanism of the optimum design point was simulated in LS-DYNA. The main contribution of this study was to successfully model the CFRP tube and damage mechanism using appropriate material constitutive model’s parameters and present the multi-objective method to find the optimum crashworthy design of the CFRP tube.


Introduction
Having better crashworthy parameters for the structure include energy absorption characteristics has been widely considered a key factor to improve the safety of the occupants during a crash in the vehicle or aerospace industry [1][2][3][4]. Some of the main typical energy-absorbing structures are thin-walled structures that have proven fairly effective and been widely studied and employed in vehicles [5][6][7]. To enhance the energy absorption capacity without adding much weight, composite materials have become very effective and popular in many engineering fields [8][9][10][11]. Among composite materials, carbon-reinforced composite tubes are considered as substantially efficient energy-absorbing components that have widely applied in the automotive and aerospace industry. Their crashworthiness behavior has been the subject of several studies in recent years [12][13][14]. A large number of experimental and numerical studies have been conducted to gain a better understanding of the failure mechanism of composite tubes [15][16][17]. Their results indicate that the mechanisms of energy absorption and dissipation for composite tubes are far more complex than those observed in conventional materials. The composite tubes absorb the impact energy through a combination of complex damage mechanisms including fiber breakage, matrix cracking, delamination, and friction [18]. Therefore, developing the numerical models that can replicate the impact response of these materials can provide a better understanding of energy absorption characteristics and crush failure mode for the composite tube.
Noticeably, FE simulation is an efficient and economic method to model the crushing mechanism of the composite tube [19]. Early simulation works were limited by computational technology [20], but recently, as computational power has improved significantly, extensive work has been done on numerical simulations using explicit FE codes such as LS-DYNA and ABAQUS [21][22][23]. These studies have been carried out for special geometric shapes such as square or circular tubes, frusta, honeycomb, etc [24]. Even though crash simulations using the FE method have become routine work in the automotive industry, their development for composite structures and failure prediction have yet to be explored [25,26] because crush simulation and failure prediction, even under quasi-static loading, is a difficult task. Recently, this topic has been the subject of several studies [19,21,[27][28][29][30]. They tried to characterize the failure mechanism from different aspects, including but not limited to the micro-mechanics model, fiber types, geometry of the tube, and how they change the crushing mechanism. Although these published studies have made significant contributions to our current understanding of composite crash behavior, only a few attempts have been made to optimize those behaviors. In addition, the failure mechanism and damage sequences of CFRP tubes are still not well understood [15]. This study focuses on characteristics of the failure mechanism for CFRP composite tube under axial load using appropriate material constitutive model and then optimize its crashworthy behavior.
On the other hand, design optimization of vehicle structures considering the maximum capacity of energy absorption with minimizing the mass and cost of the vehicle is one of the key areas in the vehicle industry [14]. Over a decade, there has been extensive research on genetic algorithms to design artificial neural networks. Also, many research studies have been conducted to implement evolutionary methods in identifying the system parameters [31,32]. Among them, the group method of data handling (GMDH) algorithm is a self-organizing approach that characterizes the complicated models based on their performance as a set of multiple inputs and single output pairs [33,34]. This method was proposed by Ivakhnenko [35] and can model complex systems without having prior knowledge about that system. In addition, Because GAs is useful for dealing with complex problems that have a large number of variables, their application in the GMDH-type neural network was investigated [36].
In this paper, first, the FE model of the CFRP tube was developed and using an appropriate failure criterion that can represent the crush characteristics. The results were then validated using the experimental data reported in [37]. The fiber direction and thickness of the tube were considered as input variables and several simulations were conducted to generate enough data to construct the GMDH-type neural network. The polynomial expression of the three outputs was extracted using GMDH models which were then used to perform the Pareto-based multi-objective optimizations. These results were then used for the Pareto-based multi-objective optimization which provides some valuable insights into the optimum crashworthy design of the CFRP tube. The results can be applied in any structural design that requires crash safety such as the automotive and aerospace industry. The main contribution of this study was to successfully model the CFRP tube and damage mechanism using appropriate material constitutive model's parameters and present the multi-objective method to optimize the crashworthy behavior of the CFRP composite tube.

Finite Element Model
A FE model of an unidirectional carbon/epoxy composite tube was developed using LS-DYNA software (Livemore Software Technology, Livermore, CA, USA) and it was validated against the experimental data published by Kim, Yoon [37]. Figure 1 shows the dimensions of the [90/0]7 tube that was utilized in that experimental study. The cylinder had a 30 mm internal diameter and a height of 100 mm with seven layers of carbon fibers with a thickness of 0.13 mm for each layer. The total mass of the tube was 28.3 g. Table 1 lists the mechanical properties of the composite tube and the resin epoxy that was extracted from the experimental tests [37]. The FE model has three main components: a cylinder using 19900 shell elements with a size of 0.5 mm, a loading plate using solid elements, and a stationary plate using the Rigid_Wall option. The composite tube was modeled using 14 integration points. The 60 × 60 × 10 mm 3 rigid plate was modeled to apply the quasi-static load. The details of the geometry and FE model components are illustrated in Figure 1. The details of the trigger mechanism are illustrated in this figure. The three layers of shell elements were used to obtain the closest possible angle in experimental model. The failure mechanism of the composite tube was modeled using the Tsai-Wu theory [38]. The MAT_ENHANCED_COMPOSITE_DAMAGE (MAT 54/55 keyword in LS-DYNA) was used to model the damage in the composite tube. Based on the Tsai-Wu theory the failure criterion for the compressive fiber mode, tensile matrix mode, and the combination of tensile and compressive matrix mode are given as following equations respectively: If it fails, then Ea = Eb = Gab = vba = vab = 0 For the compressive longitudinal direction: If failed, then Ea = vba = vab = 0; for the compressive transverse direction: If failed, then Ea = Eb = Gab = vba = vab = 0; where the σ ij is stress in the (ij) direction, Ea is Young's modulus in the (i) direction, Gij is in-plane shear modulus in the (ij) direction, vij is Poisson's ratio in the (ij) direction, S is shear strength and X and Y is strength, respectively.
Other parameters that have a significant influence on failure occurs during the simulation cannot be extracted through the experimental tests. Therefore, their values were obtained through the trial and error process which are shown in Table 2. These parameters are the maximum strain for the fibers in tensile (DFAILT), compression (DFAILC) mode, and the maximum strain for the matrix for the tensile/compression (DFAILM) and shear (DFAILS) modes. The load was applied as a quasi-static with a maximum velocity of 0.1 mm/s. Figure 2 shows the velocity profile implemented in our simulation over time. Three different contacts were used for this simulation include NODE_TO_SURFACE between the loading plate and tube, AUTOMATIC_SINGLE_SURFACE between all components, and ERODING_SINGLE_SURFACE to prevent the penetration of the composite layers. Because the explicit solver was selected in our analysis, the mass scaling option was also used to improve the simulation time.   The strength reduction parameters were used to degrade the pristine fiber strength of a ply under a compressive load. The SOFT parameter, which cannot be measured experimentally, was selected by trial and error during several crash simulations. It should be noted that studies have shown that while the material strength parameters affect the stress-strain curve significantly, changing the material strength has a very small effect on the total energy of simulation [27]. In addition, the failure strain was identified as the most critical parameter for the MAT 54 that affects the mechanical responses. Since for the entire numerical simulations in this study only one type of composite tube and load condition was selected to investigate the crashworthy behavior, the values of parameters may not be applicable when those parameters change. Figure 3 shows the results of the FE simulation results at different moments and the first 80 mm of the predicted load-displacement profile for both experiment and simulation. The results show a good correlation between the FE and experimental results. The two failure mechanisms including the transverse shearing and lamina bending (delamination) were observed which are illustrated in Figure  3. Table 3 indicates that the primary evaluation parameters calculated from the experimental and numerical analysis are close to each other.

Modeling Using GMDH Neural Network
The classical GMDH is a set of neurons in which different pairs of them in each layer are connected through a quadratic polynomial. Hence, this will produce new neurons for the next layer which can be used to map the inputs to output parameters. For example, consider a function f with multi-input X = (x1, x2, x3, , xn) and a single output yi and the ̂ is the function that approximates the predicted output value: The trained GMDH-type neural network that can predict the output values given M observations are: The goal is now to minimize the difference between the predicted value using GMDH-type and actual output: Based on this algorithm, the general form of the connection between the inputs and output variables are expressed with the following equation which is known as Ivankhnenko polynomial [32]: In most applications, the quadratic form was used only for two variables which can be written as: In Equation (8), all the coefficients ai are calculated using the regression method [35] to minimize the differences between the actual output y and approximated � for the input variables of xi and xj. This can be expressed with the following equation: All the combinations of the two independent variables from n input variables are taken into account for the basic form of the GMDH to construct the regression polynomial similar to Equation Consequently, using the quadratic sub-expression in the form of Equation (8) the following matrix can be obtained: where the variables are: To calculate the coefficients, the least square method from multiple regression analyses can be used which leads to the Equation (15). This allows us to determine the vector of the best coefficients for the quadratic equations for the whole set of M data. This procedure is repeated for each neuron in the next hidden layer based on the topology of the neural network. To prevent the singularity and round errors for the solutions, the singular value decomposition (SVD) method was implemented [32]: In this study, the generalized structure GMDH method (GS-GMDH) using a genetic algorithm was used to train the model which was proposed by Nariman-Zadeh, Darvizeh [31]. The genome or chromosome which presents the topology of the GMDH-type network consists of a symbolic string where the input variables are alphabetical names. This shows how the genetic algorithm was taken into account for the design of GMDH neural networks. Figure 4 shows an example of this network with four input variables and a single output with two hidden layers which corresponds to the length of 2 2+1 =8 genes. Further details of the mathematical model and the conditions of this method can be found in Atashkari, Nariman-Zadeh [36]. These two main concepts include a hybrid GA and SVD are involved to optimally design such a polynomial neural network [34]. The method that was used in that study was successfully used in this paper to obtain the polynomial models of the energy absorption, maximum force, and critical buckling force. The design variables in this study are the fiber angles (θ1, θ2) and the number of layers n. These variables range are from 0°-180° degrees for fiber angle and 5-9 for the number of layers. The GMDH-type models have shown a promising prediction of those outputs during the training process which will be presented in the following section.

GMDH Modeling of Crashworthiness Parameters
The validated numerical model was used to analyze the crashworthiness behavior of the FRP composite tube considering the fiber angles and the tube's thickness as important parameters. More than 100 FE simulations were conducted using different θ1 and θ2 and n (number of layers) ranged from θ1 = θ2 = [0-180] and n = [5][6][7][8][9] respectively. The outputs were the absorbed energy (E), maximum force (Pmax), and critical buckling force (Pcr). For each thickness, the polynomial function of each output was identified. It should be noted that, in some simulations where the direction of the fibers was θ1 and θ2 ≤ 40°, the simulation predicted a catastrophic failure mode which is not out interest for crashworthiness design. Also, by increasing the number of layers, the values of E, Pmax, and Pcr will go up. The architecture of the layers was not identified as the main factor that influences the output parameters. For example, the differences in the energy absorption between the tubes with [0/90]11 and [90/0]11 fiber direction were only 2%.
Here the results for the composite tube with the number of layers of n = 9 are presented. Figure  5 shows the feedforward GMDH-type networks that was used for three outputs. The GMDH type neural network that was used for the input-output data set generated from numerical simulations were presented in Table 4. The population of 13 with a crossover probability of 0.7 and a mutation probability of 0.07 was used for 1000 generations using the non-commercial code GEvoM.  Table 4 summarize the results inputs and outputs resulted from FE simulations and GMDHtype (for n = 9 layers). Figure 6 shows the comparison between the two model's outputs including critical buckling force, maximum force, and energy absorption. The results indicate that the GMDHtype neural network was adequate enough to predict the objective values for different inputs.  Figure 6. Example of comparison of the numerical FE results with predicted values using the GMDH neural network for three outputs for n = 9.
The error between the predicted values and numerical FE simulation for three outputs were quantified using three different methods and presented in Table 5. The root mean squared error (RMSE), the absolute fraction of variance (R 2 ), and the absolute percentage (MAPE) were measured based on the following equations where Yei and Ymi are the desired output from FE simulation and predicted GMDH respectively:

Multi-Objective Optimization
In this section, the multi-objective optimization has been conducted to find the vector of the decision variables that satisfies the constraints and provides the acceptable values for all the objectives function. The objectives are to maximize specific energy absorption (maximize energy E and minimize weight W) while minimizing the maximum and critical forces. The polynomial neural network models obtained from the previous section (shown in Appendix 1) were deployed to perform this optimization. Each output is defined as a set of polynomial functions of two variables. These objectives substantially influence the crashworthiness characteristics of the composite tube. The constraints of a minimum of 20 kN and a maximum of 60 kN were defined based on literature [39] for the safety of the vehicle and occupants respectively. It should be noted that the weight of the tube is only the function of the number of layers: Objectives: Constrains: The evolutionary process of Pareto multi-objective optimization has been conducted using the modified NSGA-II method [40]. In this method, the population size of 60, with the crossover probability (Pc) of 0.7, and the mutation probability (Pm) of 0.07 was defined. Figure 7 shows the results of all four objectives resulted from this optimization process. It also shows the best possible design point that has the minimum weight with a high capability of energy absorption. The plots are showing the existed relationship between the outputs. For example, energy absorption increased when the peak force went up. The highlighted point it the optimum decision vector which satisfies all the objectives and constraints defined in this analysis. The optimized values were obtained for the tube with five layers and a total mass of 19.6 g with the fiber direction of θ1 and θ2 equal to 0 and ±86 degrees, respectively. Furthermore, the optimum design was also simulated in LS-DYNA and Figure 8 shows the failure mechanism of this composite tube with the [0/±86]5 layers. The progressive crushing mechanism, which according to the literature [8,19] has the highest energy absorption capability, was observed from this simulation. Furthermore, the results from the FE simulation is very close to the values predicted by the GHMD-type neural network which confirms the reliability of this model for the optimization process (see Table 6). The FE simulation and GMDG prediction of crashworthy behavior of optimum design point showed that we can reach the same SEA (72.6 kJ/kg) as our experimental model, by using the aforementioned fiber direction and number of layers while decreasing the weight by 30%.
It should be noted that since many parameters contribute to the energy absorption of the tube such as a number of layers, fiber angle, their orientation, the geometry, and loading condition, it is impossible to generalize the results for CFRP tubes and expect to have the same crush behavior if those parameters change. However, several studies showed that the trigger mechanism, similar fiber direction will increase the energy absorption capability of the composite tube [8,15,19,21].

Conclusions
This study focuses on the characteristics of the failure mechanism for a CFRP composite tube under axial load using the appropriate material constitutive model and then optimizes its crashworthiness behavior. The FE model of the FRP composite tube and failure mechanism were successfully developed and validated in LS-DYNA. Promising results were provided by FE simulations which indicate that the Tsai-Wu failure criteria can effectively represent the crush characteristics of the carbon/epoxy composite tubes under a quasi-static axial load. Furthermore, the GS-GMDH-type neural network has been successfully used to derive the polynomial expression of the crashworthy output parameters include energy absorption, peak force, and critical buckling force as a function of fiber direction. These results were then used for the Pareto-based multi-objective optimization which provides some valuable insights into the optimum crashworthy design of the CFRP tube. The results of the optimum design point showed the same specific energy absorption (72.6 kJ/kg) while its weight decreased by 30% percent. Aside from the experimental and numerical limitations of this study, this approach can be applied in any structural design that requires crash safety such as either the automotive or aerospace industry.