Crashworthiness design of novel CFRP functionally-graded conical square structures with variable circumscribed diameters

Abstract

Carbon fiber reinforced plastic (CFRP) thin-walled structures have been widely acknowledged as superior energy absorbers due to comprehensive advantages on crashworthiness and light weight. However, few studies have been reported concerning CFRP thin-walled structures with nonlinearly variable sectional dimensions along the longitudinal direction. Hence, in the present study, a novel CFRP functionally-graded conical square thin-walled structure with variable circumscribed diameter (functionally-graded conical square tube, FGCST) is proposed and comprehensively investigated. Specifically, finite element (FE) models of FGCST under multi-angle crushing loading are first built and validated. Second, a comprehensive crashworthiness comparison among the FGCST, the conventional CFRP conical square tube (CST) and the CFRP straight square tube (SST) under multi-angle crushing loading is executed to explore the relative merits of FGCST over others. Third, a parameter study is implemented on FGCST to explore the action laws of its structural parameters on its overall crashworthiness under multi-angle crushing loading. Finally, the FGCST is multi-objective crashworthiness optimized for maximizing overall specific energy absorption (SEA_θ) and minimizing overall peak crushing force (PCF⁰) simultaneously through Taguchi method coupled with grey relational analysis (GRA). The optimized FGCST obtains 6.06% higher of SEA_θ and meanwhile 40.58% lower of PCF⁰ compared with the baseline design, thus demonstrating good potential as a superior candidate for vehicle collision energy absorbers.

Introduction

Research, design, and application of thin-walled structures with large energy absorption capacity yet light weight, are of great significance for synergistically improving the crashworthiness and lightweight level of the vehicle body composed of a large number of thin-walled structures (Figure 1), thereby promoting safety, energy conservation and environmental protection of vehicles.^1,2 Traditional thin-walled structures of the vehicle body are mostly made of metal materials, such as steel, aluminum, etc., which are not so high as carbon fiber reinforced plastic (CFRP) in terms of specific strength and specific energy absorption (SEA) capacity to the weight. Therefore, CFRP thin-walled structures demonstrate good potential to replace traditional metal thin-walled structures in vehicle bodies, whose research has become a current research hotspot particularly in the field of automotive safety and lightweight design.^3–5

Figure 1.

Thin-wall structures of vehicle body under multi-angle crushing loading.

To better understand the crushing mechanism of CFRP thin-walled structures and fully exploit their potentials as energy absorbers, extensive experimental, theoretical and numerical studies have been conducted during past decades on crashworthiness of CFRP thin-walled structures with various sections, such as circular,^6–10 square,^11–13 polygonal,^14,15 multi-cell,^16,17 top-hat, ¹⁸ etc. For instance, Sun et al. ¹⁹ conducted experimental study on crashworthiness of CFRP circular tubes under axial quasi-static crushing. Mamalis et al. ²⁰ conducted experimental study on crashworthiness of CFRP thin-walled square tubes subjected to static and dynamic axial compressive loading. Jiang et al., ²¹ Liu et al. ²² and Song et al. ²³ conducted research on energy absorption characteristics of perforated CFRP square tube structures. Xiao et al. ²⁴ carried out experimental and numerical study of CFRP hat-shaped structures under quasi-static axial crushing. Liu et al. ²⁵ investigated the energy absorption of CFRP bio-inspired multi-cell square tubes under quasi-static axial crushing. Zhu et al. ¹¹ executed numerical study on CFRP thin-walled tubes with different sectional profiles subjected to quasi-static crushing, including triangle, square, pentagon, hexagon and circle, etc.

Nevertheless, to date, the crashworthiness research of CFRP thin-walled structures mainly focused on structures with constant or linearly variable sectional dimensions along the longitudinal direction, such as conventional straight tubes (STs) or conical tubes (CTs). To be the authors’ best knowledge, basically there have not been any related reports concerning CFRP thin-walled structures with variable sectional dimensions along the longitudinal direction in a nonlinear gradient, which however possess better designability and thus greater potential as superior energy absorbers due to introduction of functional gradient.^26–29 In addition, most existed studies concerning crashworthiness of CFRP thin-walled structures focused merely on a single load direction, ignoring the uncertainty of load angles which affects crashworthiness largely, as exampled in Figure 1 with thin-wall structures of the vehicle body under multi-angle crushing loading in actual vehicle collisions.^29–31

Hence to make up for the above deficiency, the present study first proposes a novel CFRP functionally-graded conical square tube (FGCST) with variable circumscribed diameters along the longitudinal direction, for which a comprehensive crashworthiness comparison, a parameter study and a multi-objective crashworthiness optimization are carried out in succession. Specifically, the finite element (FE) models of FGCST under multi-angle crushing loading are first built and validated. Second, a comprehensive crashworthiness comparison among the FGCST, the conventional CFRP conical square tube (CST) and the CFRP straight square tube (SST) under multi-angle crushing loading is executed to explore the relative merits of FGCST over others. Third, a parameter study is implemented on FGCST to explore the action laws of its structural parameters on its overall crashworthiness under multi-angle crushing loading. Finally, the FGCST is multi-objective crashworthiness optimized for maximizing overall specific energy absorption (SEA_θ) and minimizing overall peak crushing force (PCF⁰) simultaneously through Taguchi method coupled with grey relational analysis (GRA).

Geometric configuration

Unlike the conventional straight square tube (SST) or conical square tube (CST) with simply straight walls, the proposed novel FGCST is characterized by functionally-graded circumscribed diameter along its longitudinal direction, thereby presenting either straight or curved tube wall depending on the gradient exponent. Figure 2(a) illustrates the geometrical configuration and boundary conditions of FGCST under multi-angle (axial and oblique) crushing loading. In the axial (or longitudinal) direction, the circumscribed diameter of FGCST varies continuously according to a predefined gradient function, while the wall thickness keeps uniform. Specifically, the circumscribed diameter variation of FGCST along the axial direction can be defined through the following power-law function^32,33:

$$R\left(x\right)=R_t+(R_b-R_t){\left(\frac{x}{L}\right)}^n$$

(1)

where $R\left(x\right)$ denotes the circumscribed radius of the tube at current section; $x$ is the distance from current section to the top end face of the tube; $R_t$ and $R_b$ denote the top and bottom circumscribed radius of the tube, respectively; L is the total length of the tube; $n$ is the gradient exponent controlling the variation of the circumscribed radius, which varies between 0.1 and 10.

Figure 2.

FGCST under multi-angle crushing loading: (a) geometric diagram of FGCST and (b) variation modes of circumscribed radius.

To improve load uniformity and resist oblique loads, the circumscribed radius at the top end face of FGCST is set less than or equal to that at the bottom end face ( $R_t\le R_b$ ). Obviously, with the gradient exponent n increases from less than 1 to larger than 1, the circumscribed radius increases continuously along the length with a gradient function changing form convexity to concavity, as shown in Figure 2(b). To better demonstrate the influence of gradient exponent n on the geometric configuration of FGCST, the geometric configuration diagrams of FGCSTs with different gradient exponent n are shown in Figure 3, herein all the displayed FGCSTs have the same structural parameters as the following: $R_t$  = 40 mm, $R_b$  = 60 mm, L = 150  mm, and t = 1 mm. Note that the FGCST transforms to conventional CST when n equals 1, while it transforms to conventional SST when n equals 0 ( $R_t=R_b$ ). In other words, the conventional CST and SST could be treated as special cases of FGCST.

Figure 3.

FGCSTs with different gradient exponent n.

Finite element modeling and verification

Constitutive models

The finite element (FE) models in the present study are developed using FE code ABAQUS. To capture intra-laminar and inter-laminar failure details of CFRP thin-walled structures simultaneously with a relatively simple constitutive relation yet high computational efficiency, the macro-mechanical modeling approach of stacked shell is adopted to monitor the crushing behaviors of FGCST under multi-angle crushing loading. ¹¹ Specifically, the intra-laminar failure mechanism is relevant to the fiber damage and matrix cracking, while the inter-laminar failure mechanism is connected with the delamination failure between the layers. The constitutive models applied for modeling the intra-laminar and inter-laminar failure mechanisms are outlined as follows. ¹¹

Intra-laminar failure

The intra-laminar failure model of CFRP material is actualized through a user defined material subroutine VUMAT, in which the fiber damage under tensile or compressive loading and plastic deformation under shear loading are involved. The model is composed of several material properties and damage parameters. For plain woven fabric reinforcements with orthogonal fiber directions considered in the present study, the intra-laminar stress-strain relation can be expressed as^11,35:

$$\left\{\begin{array}{c}{\epsilon }_{11}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{22}}\\ \multicolumn{1}{c}{{\epsilon }_{12}^{\mathrm{el}}}\end{array}\right\}=\left[\begin{array}{ccc}\frac{1}{(1-d_1)E_1}\hfill & -\frac{v_{12}}{E_1}\hfill & 0\hfill \\ \multicolumn{1}{c}{-\frac{v_{21}}{E_2}}& \frac{1}{(1-d_2)E_2}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{1}{(1-d_{12})2G_{12}}\hfill \end{array}\right]\left\{\begin{array}{c}{\sigma }_{11}\hfill \\ \multicolumn{1}{c}{{\sigma }_{22}}\\ \multicolumn{1}{c}{{\sigma }_{12}}\end{array}\right\}$$

(2)

where subscript 1 and subscript 2 represent the fiber longitudinal direction and transversal direction, respectively; $\mathbf{\sigma }={\left\{\begin{array}{ccc}{\sigma }_{11}\hfill & {\sigma }_{22}\hfill & {\sigma }_{12}\hfill \end{array}\right\}}^T$ and $\mathit{\epsilon }={\left\{\begin{array}{ccc}{\epsilon }_{11}\hfill & {\epsilon }_{22}\hfill & {\epsilon }_{12}^{\mathrm{el}}\hfill \end{array}\right\}}^T$ denote the stress vector and elastic strain vector, respectively; $E_1$ and $E_2$ denote the elastic modulus in the longitudinal and transversal direction, respectively, while $G_{12}$ denote the in-plane shear modulus; $v_{12}$ and $v_{21}$ denote the in-plane Poisson’s ratios; $d_1$ and $d_2$ denote the damage parameter associated with fiber fracture in the longitudinal and transversal direction, respectively, while $d_{12}$ stands for the damage parameter related to matrix micro-cracking under shear loading, value of which varying between 0 and 1 stands for the stiffness degradation caused by micro-damage in the material.

The initiation and evolution of the individual damage mechanisms are described by means of the corresponding effective stress $\hat{\sigma }$ , as ( $i=1,2$ )^11,35:

$${\hat{\sigma }}_{i+}=\frac{〈{\sigma }_{\mathrm{ii}}〉}{1-d_{i+}}{\hat{\sigma }}_{i-}=\frac{〈-{\sigma }_{\mathrm{ii}}〉}{1-d_{i-}}{\hat{\sigma }}_{12}=\frac{{\sigma }_{12}}{1-d_{12}}$$

(3)

where $d_{i+}$ and $d_{i-}$ denote the tensile and compressive damage variable of the fiber in the ith direction, respectively; angle bracket < > in the equation represents the Macquarie operator. ¹¹

In any given time, the elastic domain can be defined in terms of the damage activation function F _i+ , F _i− and F₁₂ as follows ( $i=1,2$ )^11,35 :

$$F_{i+}=\frac{{\hat{\sigma }}_{\mathrm{ii}}}{X_{i+}}{F}_{i-}=-\frac{{\hat{\sigma }}_{\mathrm{ii}}}{X_{i-}}{F}_{12}=\left|\frac{{\hat{\sigma }}_{12}}{S}\right|$$

(4)

where $X_{i+}$ and $X_{i-}$ denote the tensile and compressive strength in the ith direction, $S$ denotes in-plane shear strength, F _i+ , F _i− and F₁₂ denote the failure coefficients corresponding to each failure mode, any value of which equal to or greater than 1.0 indicates the corresponding failure criterion has been met.

Once the fiber appears initial damage, then the degradation of fiber properties would take place according to the damage evolution law, which could be expressed as follows^11,35:

$$d_{\alpha }=1-\frac{1}{r_{\alpha }}\exp [-\frac{2g_0^{\alpha }{L}_c}{G_{\mathrm{fc}}^{\alpha }-g_0^{\alpha }{L}_c}(r_{\alpha }-1)](a=1\pm ,2\pm )$$

(5)

where $L_c$ is the characteristic length of the element, $G_{\mathrm{fc}}^{\alpha }$ denotes the fracture energy per unit area under uniaxial tensile or compressive loading, $g_0^{\alpha }$ denotes the elastic energy per unit volume at the point of damage initiation, $r_{\alpha }$ denotes the damage threshold for fiber damage and can be defined as:

$$r_{i+}=max(1,F_{i+})$$

(6)

$$r_{i-}=max(1,F_{i-})$$

(7)

As for the evolution of shear damage, it can be defined as:

$$d_{12}=min[{\beta }_{12}\ln \left(r_{12}\right),d_{12}^{\max }]$$

(8)

where ${\beta }_{12}$ and $d_{12}^{\max }$ denote the material properties, which can be obtained through in-plane shear test; $r_{12}$ denotes the damage threshold for shear:

$$r_{12}=max(1,F_{12})$$

(9)

The damage thresholds, $r_{i+}$ , $r_{i-}$ , and $r_{12}$ , indicate the historical damage of material. They were initially install to 1 before the initial damage and equal to each corresponding failure coefficient when damage evolution starts.

The in-plane shear response of the material is dominated by the plastic behavior due to the matrix micro-cracking. To account for this, a plasticity model with a yield function and a hardening law is employed as follows^11,35:

$$F=\left|{\hat{\sigma }}_{12}\right|-{\hat{\sigma }}_0({\overline{\epsilon }}^{\mathrm{pl}})\le 0$$

(10)

$${\hat{\sigma }}_0({\overline{\epsilon }}^{\mathrm{pl}})={\hat{\sigma }}_{y0}+C{\left({\overline{\epsilon }}^{\mathrm{pl}}\right)}^p$$

(11)

where ${\overline{\epsilon }}^{\mathrm{pl}}$ is the equivalent plastic strain; ${\hat{\sigma }}_{y0}$ is the initial effective yield shear stress; $C$ and $p$ denote the hardening coefficients of the material, which can be obtained through calibration.

Table 1 presents the material parameters of CFRP, including the elasticity constants, the tensile and compression strengths that generally obtained by uniaxial tensile and compressive tests, as well as the in-plane shear response parameters that normally obtained through in-plane shear test. Table 2 gives the critical fracture energy for intra-laminar damage evolution. In the present study, all these above parameters to establish the intra-laminar failure model are referred to Zhu et al., ³⁴ where the similar material is used.

Table 1.

Material parameters of CFRP. ³⁴

Var.	Description	Value	Var.	Description	Value
E1	Longitudinal modulus (GPa)	57	X2−	Transverse compressive strength (MPa)	512
E2	Transverse modulus (GPa)	57	S	In-plane shear strength (MPa)	71
v12	Poisson’s ratio	0.067	${\hat{\sigma }}_{y0}$	Initial effective shear yield stress (MPa)	115
G12	Shear modulus (GPa)	8.4	${\beta }_{12}$	Coefficient in the shear damage	0.154
X1+	Longitudinal tensile strength (MPa)	679	$C$	Coefficient in the hardening equation	3080
X1−	Longitudinal compressive strength (MPa)	679	p	Power term in hardening equation	0.8106
X2+	Transverse tensile strength (MPa)	512

Table 2.

The critical fracture energy for intra-laminar damage evolution. ³⁴

Variable	Description	Value
$G_{\mathrm{fc}}^{1+}$	Longitudinal tensile fracture energy (kJ/m2)	155
$G_{\mathrm{fc}}^{1-}$	Longitudinal compression fracture energy (kJ/m2)	255
$G_{\mathrm{fc}}^{2+}$	Transverse tensile fracture energy (kJ/m2)	155
$G_{\mathrm{fc}}^{2-}$	Transverse compression fracture energy (kJ/m2)	255

Inter-laminar failure

In order to monitor the adhesion and debonding of adjacent plies of CFRP, a bonding model based on the generalized traction-separation law including corresponding damage criterion and evolution is adopted to define the cohesive zone. The behavior of its elastic phase can be described by the following formula^11,35:

$$\mathit{t}=\left\{\begin{array}{c}{t}_n\hfill \\ \multicolumn{1}{c}{t_s}\\ \multicolumn{1}{c}{t_t}\end{array}\right\}=\left[\begin{array}{ccc}{K}_n\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& K_s\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & K_t\hfill \end{array}\right]\left\{\begin{array}{c}{\epsilon }_n\hfill \\ \multicolumn{1}{c}{{\epsilon }_s}\\ \multicolumn{1}{c}{{\epsilon }_t}\end{array}\right\}=\mathit{K}\mathit{\epsilon }$$

(12)

where $\mathit{t}$ and $\mathit{\epsilon }$ denote traction stress vector and strain vector, respectively; ${\epsilon }_i$ and $t_i(i=n,s,t)$ represent the traction displacement and traction stress in normal (n), first shear ( $s$ ) and second shear (t) directions, and $\mathit{K}$ represents the stiffness matrix of interaction.

To describe the initiation of delamination damage, a quadratic interaction function is employed as^11,35:

$${\left(\frac{〈t_n〉}{t_n^0}\right)}^2+{\left(\frac{t_s}{t_s^0}\right)}^2+{\left(\frac{t_t}{t_t^0}\right)}^2=1$$

(13)

where $t_i^0(i=n,s,t)$ represents the peak value of traction stress in each direction. When the function reaches a value of 1, the damage occurs and evolves according to the following rule^11,35:

$$G_n^C+(G_s^C-G_n^C){\left(\frac{G_s^C+G_t^C}{G_n^C+G_s^C+G_t^C}\right)}^{\mathrm{\eta }}=G^C$$

(14)

where $G_n^C,G_s^C,G_n^C$ demote the fracture energy in the normal, first shear and second shear direction, respectively, η is the cohesive layer attribute parameter equal to 2.284. Table 3 demonstrates the inter-laminar damage material parameters. ⁴

Table 3.

Inter-laminar damage parameters. ⁴

Variable	Description	Value
$G_n^C$	Mode I fracture energy (kJ/m2)	0.202
$G_s^C,G_n^C$	Mode II fracture energy (kJ/m2)	0.891
$t_n^0$	Normal strength (MPa)	4.8
$t_s^0,t_t^0$	Shear direction strength (MPa)	14.8

Finite element (FE) modeling

The FGCST is placed on a fixed lower plate with the bottom end fully constrained, while the top end is compressed by a moving upper plate at an initial velocity of 2 m/s from variable multiple loading angles (0°, 10°, 20°, 30°). Note that in the numerical simulation a relatively higher loading rate is used in order to balance computational efficiency and accuracy. Both the lower and upper plate are modeled as rigid surfaces. To accurately simulate the crushing behaviors of FGCST, the multi-layer shell element (S4R) is employed to model the CFRP layers, and the cohesive zone element (COH3D8) is utilized to model the bonding between the layers. To improve computational efficiency, the FGCST is meshed with an element size of 1.0 × 1.5 mm, which is sufficient to predict the crushing behaviors of CFRP tube. The material properties and parameters for defining the intra-laminar and inter-laminar failure models are listed in Tables 1 to 3, which are employed by a built-in subroutine of ABAQUS to simulate the crushing behaviors of the CFRP tube. The contact between the tube and the pressing plate is simulated using the general contact with friction coefficient of 0.3. ³⁴ Figure 4 demonstrates the finite element models of FGCSTs under multi-angle crushing loading.

Figure 4.

Finite element models of FGCST under multi-angle crushing loading.

Verification of FE models

As mentioned earlier, the FGCST will degenerate to a conventional straight square tube (SST) when the grading exponent n equals 0, in turn, the SST can be regarded as a special case of FGCST. Therefore, the SST could be selected to perform experimental verification of FE models for FGCST. Toward this end, Figure 5 demonstrates the geometric configuration of experimental specimen of SST by Zhu et al. ³⁴ and corresponding FE model built in the present study. The SST is composed of six layer of plain woven CFRP laminates with each layer a uniform thickness of 0.235 mm, and possesses a height of 100 mm and a side length of 55.89 m. Figure 6(a) compares the crushing force versus displacement curves of the experimental test ³⁴ and FE simulation, in which the simulation curve can basically captures the experimental counterpart concerning amplitude and trend throughout the crushing stroke. Besides, the comparison of final deformation mode between experiment and simulation is also presented in Figure 6(b), from which it is easily seen that the final deformation mode of simulation basically agrees with that of experiment. To sum up, the FE models developed in the present study are of acceptable accuracy and provide considerable confidence for us to further investigate and optimize the crashworthiness of FGCST under multi-angle crushing loading.

Figure 5.

Experimental specimen and FE model of SST. ³⁴

Figure 6.

Crashworthiness comparison between experiment ³⁴ and simulation: (a) force-displacement curve and (b) final deformation mode.

Crashworthiness analysis and comparison

Crashworthiness indicators

To evaluate the crashworthiness of FGCST under multi-angle crushing loading, the following conventional crashworthiness indicators are used, including energy absorption (EA), specific energy absorption (SEA), peak crushing force (PCF) and mean crushing force (MCF). Specifically, EA denotes the absorbed energy of structures via deformation, which demonstrates the total capacity for energy absorption. Specific energy absorption (SEA) denotes the energy absorption per unit mass, which measures the energy absorption capacity of material. PCF denotes the maximum instantaneous crushing force along the crushing process, which generally measures the load severity of crushing. MCF denotes the average instantaneous crushing force along the crushing process, which could be employed to evaluate the average energy absorption capability on the past crushing stroke. The above crashworthiness indicators, as shown in Figure 7, could normally be calculated as follows:

$$\mathrm{EA}\left(d\right)={\int }_0^{d}F\left(x\right)\mathrm{dx}$$

(15)

$$\mathrm{SEA}\left(d\right)=\frac{\mathrm{EA}\left(d\right)}{M}$$

(16)

$$\mathrm{PCF}\left(d\right)=\underset{d}{max}F\left(x\right)$$

(17)

$$\mathrm{MCF}\left(d\right)=\frac{\mathrm{EA}\left(d\right)}{d}$$

(18)

where $d$ is the crushing distance, which is taken as 0.75 L in the present study, $F\left(x\right)$ is the instantaneous crushing force, M is the total mass of tube.

Figure 7.

Schematic diagram of crashworthiness indicators.

Nevertheless, due to the uncertainty of vehicle collision direction, vehicle body thin-walled structures may withstand crushing loads from different angles, which could largely affect the value of energy absorption. Hence, taking into account the effect of different loading angles, two overall crashworthiness indicators, namely $E{A}_{\theta }$ and $\mathrm{SE}{A}_{\theta }$ , are defined and adopted herein as follows:

$$E{A}_{\theta }=\sum _{j=1}^N{\omega }_{{\theta }_j}E{A}_{{\theta }_j}$$

(19)

$$\mathrm{SE}{A}_{\theta }=\sum _{j=1}^N{\omega }_{{\theta }_j}\mathrm{SE}{A}_{{\theta }_j}$$

(20)

where $E{A}_{{\theta }_j}$ and $\mathrm{SE}{A}_{{\theta }_j}$ denote the EA and SEA under specific loading angle ${\theta }_j$ ( j = 1,2,…N ), N stands for the total number of crushing loading angles considered; ${\omega }_{{\theta }_j}$ denotes the relative weight of ${\theta }_j$ indicating importance or occurrence probability, and the sum of ${\omega }_{{\theta }_j}$ under all crushing loading angles is 1.

Generally, under multi-angle crushing loading scenarios, the overall or weighted EA ( $E{A}_{\theta }$ ) and SEA ( $\mathrm{SE}{A}_{\theta }$ ) is “larger-the-better (LTB)” to dissipate crushing kinetic energy, while the PCF under loading angle of 0° (PCF⁰), which is generally the maximum among PCFs under all loading angles, is “smaller-the-better (STB)” to control crushing load severity.

Crashworthiness analysis of FGCST

Numerical crashworthiness analyses are carried out based on the validated FE models to investigate the crushing and energy absorption behaviors of FGCSTs under multi-angle crushing loading. Herein, four typical loading angles, namely 0°, 10°, 20°, and 30° are considered. For FGCSTs with gradient exponent n greater than or equal to 1, a chamfer with an unilateral angle of 45° is added at the edges of their top end face in order to trigger a relatively more stable form of fragmentation. In contrast, for FGCSTs with gradient exponent n less than 1, no chamfer is added due to they present an arc structure, which is already conducive to triggering stable form of fragmentation. The crushing force versus displacement curves of FGCSTs with different n under multi-angle (0°, 10°, 20°, and 30°) crushing loading are depicted in Figure 8, and the corresponding final deformation modes of FGCSTs are demonstrated in Table 4.

Figure 8.

Force-displacement curves of FGCSTs under multi-angle crushing loading: (a) 0°, (b) 10°, (c) 20°, and (d) 30°.

Table 4.

Final deformation modes of FGCSTs with different n.

Gradient exponent	Initial states	Final deformation modes
		0°	10°	20°	30°
n = 0.1
n = 0.6
n = 1
n = 6
n = 10

It can be seen that: (1) For most FGCSTs with different gradient exponent n, the crushing forces present not so large initial peak values, followed by alternating oscillations with relatively large global mean values, namely MCFs, which are different from conventional meta-tubes yet beneficial for energy absorption. (2) The FGCSTs with n less than 1 present relatively larger global MCFs compared with that of counterparts with n greater than or equal to 1, which could also be reflected by the final deformation modes in Table 4. It can be seen that progressive fragmentation accompanied by collapse folding and local bending relatively beneficial for energy absorption are dominant for the former, while progressive fragmentation accompanied by overall bending and four-corner cracking relatively adverse to energy absorption are dominant for the latter. (3) For all FGCSTs with different gradient exponent n, the MCF decreases along with increasing crushing loading angles, indicating the deterioration of energy absorption. This could also be reflected by the final deformation modes, which change from progressive fragmentation dominated deformation to overall bending and four-corner cracking dominated deformation with increasing crushing loading angles.

Crashworthiness comparison

In this section, a comprehensive crashworthiness comparison among the novel FGCST, the conventional CST and SST under multi-angle crushing loading is carried out to explore the relative merits and demerits of FGCST over others.

Configuration relationship of different tubes with equal mass

To objectively compare the crashworthiness of FGCST, CST and SST, it is necessary to ensure equal mass of them. To this end, it is critical to derive the geometrical relationship among them to obtain equal mass (or volume). Herein, for all above CFRP tubes, including FGCST, CST, and SST, the tube walls are composed of four layer of plain woven carbon fiber laminates with each layer a thickness of 0.25 mm, and have an identical length (L) of 150 mm and thickness (t) of 2 mm. Thus, based on the geometric parameters of FGCST shown in Figure 9, the volume of the FGCST (shadowed) $V_F$ is calculated by integration and can be given as:

$$\begin{array}{c}\multicolumn{1}{c}{V_F={\int }_0^L[{\left(\sqrt{2}(R\left(x\right)+\frac{\sqrt{2}}{2}t)\right)}^2-{\left(\sqrt{2}(R\left(x\right)-\frac{\sqrt{2}}{2}t)\right)}^2]\mathrm{dx}}\\ =4\sqrt{2}\mathrm{Lt}[R_t+\frac{(R_b-R_t)}{n+1}]\hfill \end{array}$$

(21)

where $R\left(x\right)$ is the circumscribed radius of the FGCST defined in equation (1). Herein, keep $R_t$  = 30 mm and $R_b$  = 40 mm, then the volume changes with variable gradient exponent n.

Figure 9.

Illustration of cross section area for FGCST.

Likewise, based on geometric parameters of CST and SST shown in Figure 10, the volume of CST ( $V_C$ ) and SST ( $V_S$ ) could be respectively derived as:

$$\begin{array}{c}\multicolumn{1}{c}{V_C={\int }_0^L[{\left(\sqrt{2}(R\left(x\right)+\frac{\sqrt{2}}{2}t)\right)}^2-{\left(\sqrt{2}(R\left(x\right)-\frac{\sqrt{2}}{2}t)\right)}^2]\mathrm{dx}}\\ =4\sqrt{2}\mathrm{Lt}[R_c+\frac{(40-R_c)}{2}]\hfill \end{array}$$

(22)

$$\begin{array}{c}\multicolumn{1}{c}{V_S={\int }_0^L[{\left(\sqrt{2}(R_S+\frac{\sqrt{2}}{2}t)\right)}^2-{\left(\sqrt{2}(R_S-\frac{\sqrt{2}}{2}{t}_S)\right)}^2]\mathrm{dx}}\\ =4\sqrt{2}\mathrm{Lt}{R}_S\hfill \end{array}$$

(23)

where $R\left(x\right)$ is the circumscribed radius of CST, $R_C$ is the equivalent top circumscribed radius of CST (assume $R_b$  = 40 mm unchanged), $R_S$ is the equivalent circumscribed radius of SST.

Figure 10.

Illustration of equivalent circumscribed radius for CST and SST: (a) CST and (b) SST.

Let $V_F=V_C=V_S$ , the following relationships can be obtained:

$$R_C=\frac{(2n+2)R_t-(n+1)R_b}{n+1}$$

(24)

$$R_S=R_t+\frac{(R_b-R_t)}{n+1}$$

(25)

According to equations (24) and (25), given specific geometrical parameters of FGCST and gradient exponent n, the equivalent top circumscribed radius ( $R_C$ ) of CST and equivalent circum-scribed radius ( $R_S$ ) of SST that guarantee equal mass of them can be determined uniquely, which varies along with the variation of gradient exponent n, as shown in Figure 11. Herein, gradient exponent $n$ are taken as 0.1, 0.2, 0.4, 0.6, 0.8, 1, 2, 4, 6, 8, and 10 to exectute the comparison.

Figure 11.

Equivalent circumscribed radius of CST and SST with different n.

Comparison between FGCST and CST

Figure 12 compares the crashworthiness indicators of FGCST and CST under multi-angle (0°, 10°, 20°, 30°) crushing loading, including specific energy absorption (SEA) in Figure 12(a) and peck crushing force (PCF) in Figure 12(b). It can be seen that when n is less than 1, with all values of gradient exponent n, the SEA under all loading angles of FGCST are relatively larger than that of CST, meanwhile, the PCF under axial (0°) crushing loading (PCF⁰) of FGCST are relatively smaller than that of CST. When n is greater than 1, with all values of gradient exponent n, the SEA under most loading angles of FGCST are relatively larger than that of CST, while with most values of gradient exponent n (n = 6,8,10), the PCF under axial (0°) crushing loading (PCF⁰) of FGCST are even relatively larger than that of CST. Of course, the FGCST is CST when n is equal to 1, which needs no discussion. To conclude, on the whole, the novel FGCST outperforms the CST in terms of multi-angle crashworthiness when n is less than 1, while there is no absolute advantage or disadvantage between FGCST and CST in terms of multi-angle crashworthiness when n is larger than 1, which depends on the value of n and the emphasis of SEA or PCF.

Figure 12.

Comparison of SEA and PCF between FGCST and CST with equal mass: (a) SEA and (b) PCF.

Comparison between FGCST and SST

Figure 13 compares the crashworthiness indicators of FGCST and SST under multi-angle (0°, 10°, 20°, 30°) crushing loading, including specific energy absorption (SEA) in Figure 13(a) and peck crushing force (PCF) in Figure 13(b). It can be seen that when n is less than 1, with all values of gradient exponent n, the SEA under most loading angles of FGCST are relatively larger than that of SST, meanwhile, the PCF under axial (0°) crushing loading (PCF⁰) of FGCST are relatively smaller than that of SST. In contrast, when n is greater than or equal to 1, with all values of gradient exponent n, the SEA under all loading angles of FGCST are relatively smaller than that of SST, meanwhile, the PCF under axial (0°) crushing loading (PCF⁰) of FGCST with partial values of gradient exponent n (n = 8, 10) are even relatively larger than that of SST. To conclude, on the whole, the novel FGCST outperforms the SST in terms of multi-angle crashworthiness when n is less than 1, while it underperforms the SST in terms of multi-angle crashworthiness when n is greater than or equal to 1.

Figure 13.

Comparison of SEA and PCF between FGCST and SST with equal mass: (a) SEA and (b) PCF.

Parameter study

In order to better understand the energy absorption characteristics of FGCST, the effects of structural parameters of FGCST on its crashworthiness under multi-angle crushing loading are further studied. Considering the effect of gradient exponent n has been analyzed, another two structural parameters, namely the difference between top and bottom circumscribed radius (Δ $R$ ) and the number of (45°/−45°) layers (N), are selected to carry out parameter study in this section. Note that in this section and subsequent sections, assume the baseline (or initial) design of FGCST has such structural parameters: R_t = 30 mm, R_b = 40 mm (namely Δ $R$  = 10 mm), n = 2 and N = 0.

Effect of difference between top and bottom circumscribed radius (ΔR)

For FGCST, set the gradient exponent n to 2 and keep the circumscribed radius of the bottom end face ( $R_b$ ) to 40 mm, then vary the circumscribed radius of the top end face ( $R_t$ ) within the range of 25–35 mm with a sample point every 2.5 mm, leading to a group of Δ $R$ (5, 7.5, 10, 12.5, and 15 mm). Figure 14 compares the crashworthiness indicators (EA and PCF) of FGCSTs with different Δ $R$ under multi-angle crushing loading, while Table 5 compares the corresponding final deformation modes. It can be seen that for FGCST with any Δ $R$ , both EA and PCF decrease along with increasing crushing loading angles on the whole, which once again confirms the change of crushing loading angle could significantly affect the energy absorption. In addition, for FGCST under any specific crushing loading angle, except for some local deviations, both the EA and PCF approximately decrease with increasing Δ $R$ on the whole, indicating the decrease of stiffness and energy absorption capacity.

Figure 14.

Crashworthiness indicators of FGCST with different Δ $R$ : (a) EA and (b) PCF.

Table 5.

Final deformation modes of FGCSTs with different ΔR.

Different ΔR	Initial mode	Final deformation modes
		0°	10°	20°	30°
ΔR = 15 mm
ΔR = 12.5 mm
ΔR = 10 mm
ΔR = 7.5 mm
ΔR = 5 mm

These might also be reflected by the final deformation modes in Table 5, where it can be seen that for all FGCSTs with different Δ $R$ , the proportion of local yielding in deformation increases with increasing loading angles, leading to the final deformation modes transform from approximately progressive fragmentation under axial (0°) crushing loading to approximately overall bending with four-corner splitting under 30° crushing loading, which inevitably results in a decrease in energy absorption due to local yielding is not so conducive to energy absorption as fragmentation. As for FGCSTs with different Δ $R$ under any specific crushing loading angle, the effects of Δ $R$ on the final deformation modes are relatively not so obvious.

Effect of number of (45°/−45°) layers (N)

The layer structure of CFRP laminates, particularly the layer angle, could has a significant impact on the crashworthiness of CFRP thin-walled structures like FGCST under multi-angle crushing loading. To explore this, assume the FGCST has such layer structure as (0°/90°)_(4-N) (45°/−45°) _N , where N denotes the number of (45°/−45°) layers and is taken as parameter variable varying between 0 and 4. Except for this, the remaining structural parameters are consistent with the baseline design. Then the crashworthiness indicators of FGCST with different layer structures, namely N, under multi-angle crushing loading are compared as shown in Figure 15, and corresponding final deformation modes are also compared in Table 6. It can be seen that for FGCST with any number of (45°/−45°) layers (N), both EA and PCF decrease along with increasing crushing loading angles on the whole, which once again confirms the change of crushing loading angle could significantly affect the energy absorption. In contrast, the effects of layer structure on the EA and PCF under multi-angle crushing loading are more sophisticated. On the whole, the FGCST with larger number of (45°/−45°) layers (N) tends to present relatively higher EA and PCF under crushing loading angles of 0° and 30°, while it tends to present relatively lower EA and PCF under crushing loading angles of 10° and 20°, indicating the effect of number of (45°/−45°) layers (N) on crashworthiness varies with crushing loading angle.

Figure 15.

Crashworthiness indicators of FGCSTs with different N: (a) EA and (b) PCF.

Table 6.

Final deformation modes of FGCSTs with different N.

Different N	Final deformation modes
	0°	10°	20°	30°
N = 0
N = 1
N = 2
N = 3
N = 4

These might also be reflected by the final deformation modes in Table 6, where it can be seen that for all FGCSTs with different N, the final deformation modes transform from approximately progressive fragmentation under axial (0°) loading to approximately overall bending with four-corner splitting under 30° crushing loading, which inevitably results in a decrease in energy absorption similarly. As for FGCSTs with different N under any specific crushing loading angle, the effects of N on the final deformation modes vary with crushing loading angles and are relatively not so obvious. To sum up, the selection of the number of (45°/−45°) layers (N) in CFRP thin-walled structure designs depend on actual needs.

Multi-objective crashworthiness optimization of FGCST

To fully tap into the potential of FGCST for better overall crashworthiness under multi-angle crushing loading, it is necessary to explore optimization methods for it. In general, for energy absorbing thin-walled structures like FGCST, SEA is “larger-the-better (LTB)” to improve energy absorbing efficiency and capacity, while PCF is “smaller-the-better (STB)” to control load severity. However, as this article and many studies have revealed, SEA and PCF are generally competing or conflicting indicators, hence multi-objective optimization is suitable and necessary.^36,37 Compared with many other proposed optimization methods, the Taguchi method, which adopts orthogonal experimental design and processes experimental data based on signal-to-noise ratio (SNR) that to be maximized, is relatively more efficient and simple to assess robustness of experiments and effects of factors on responses, as well as to explore optimal combination of factor levels. ³⁸ Nevertheless, the traditional Taguchi method is normally limited to handling single-objective optimization issues, which is not suitable to address multi-objective crashworthiness optimization of FGCST in the present study. However, this limitation could be eliminated by coupling Taguchi method with grey relational analysis (GRA), due to GRA could convert multiple responses to a single response called grey relational grade (GRG), which represents the relative closeness of a candidate to the ideal referential alternative and generally is “larger-the-better (LTB).” Consequently, this paper carries out multi-objective crashworthiness optimization of FGCST under multi-angle crushing loading using Taguchi method coupled with grey relational analysis (GRA).

Taguchi-GRA method

Step 1: Taguchi experimental design

For multi-objective crashworthiness optimization of FGCST in the present study, the aforementioned three structural parameters, including the gradient exponent (n), the difference between top and bottom circumscribed radius (Δ $R$ ) and the number of (45°/−45°) layers (N) could be selected as three design variables. Nevertheless, for convenience in this section, two of the above parameters are slightly adjusted. Specifically: (1) keep R_b constant thus the top circumscribed radius R_t is selected as one design variable instead Δ $R$ ; (2) considering the value dimensions of the gradient exponent n between n < 1 (0.1,1) and n > 1 (1,10) are inconsistent, a new parameter m called normalized gradient exponent, is selected as another design variable calculated as:

$$m=\mathrm{lo}{\mathrm{g}}_{10}n$$

(26)

Obviously, m varies between (−1, 1) along with n varying between (0.1, 10). Therefore, herein the top circumscribed radius (R_t or TCR), the normalized gradient exponent (m or NGE) and the number of (45°/−45°) layers (N or NL) are selected as three design factors and for each factor four levels are assigned, thus a L₁₆(4³) orthogonal array (OA) could be built. Table 7 demonstrates design variables and corresponding levels, which are set according to parameter values of baseline design of FGCST, as boldly marked in the table and abbreviated as R2-m3-N1, namely R_t= 30 mm, m = 0.301 (or n = 2), N = 0. Considering multi-angle crushing scenarios, the overall SEA (SEA_θ) is taken as one response (or objective), while the PCF under 0° crushing loading (PCF⁰) is taken as an another based on the fact that PCF of FGCST generally obtains maximum under pure axial (0°) crushing loading, as disclosed before. Table 8 displays the constituted table for Taguchi experimental design and corresponding experimental results.

Table 7.

Design variables and corresponding levels.

Variables	Symbols	Levels
		1	2	3	4
Top circumscribed radius (TCR)	Rt	27.5	30	32.5	35
Normalized gradient exponent (NGE)	m	−0.301	0	0.301	0.602
Number of (45°/−45°) layers (NL)	N	0	1	2	3

Table 8.

Experimental results and corresponding SNRs.

No.	Experiment layout			Experimental results
	Rt (mm)	m	N	SEAθ (kJ/kg)(LTB)	SNR1 (LTB)	PCF0 (kN) (STB)	SNR2 (LTB)
1	27.5	−0.301	0	14.60	23.28	14.55	−23.25
2	27.5	0	1	15.43	23.76	22.62	−27.08
3	27.5	0.301	2	14.31	23.11	25.2	−28.02
4	27.5	0.602	3	14.75	23.37	22.24	−26.94
5	30	−0.301	1	14.94	23.48	12.17	−21.70
6	30	0	0	14.53	23.24	24.59	−27.81
7	30	0.301	3	15.06	23.55	25.1	−27.99
8	30	0.602	2	14.40	23.16	19.08	−25.61
9	32.5	−0.301	2	15.19	23.63	15.82	−23.98
10	32.5	0	3	14.63	23.30	24.29	−27.70
11	32.5	0.301	0	13.59	22.66	24.43	−27.75
12	32.5	0.602	1	13.70	22.73	29.85	−29.49
13	35	−0.301	3	12.83	22.16	13.13	−22.36
14	35	0	2	14.60	23.28	29.09	−29.27
15	35	0.301	1	14.13	23.00	26.51	−28.46
16	35	0.602	0	13.19	22.40	25.11	−27.99

Step 2: Signal-to-noise ratio (SNR)

After the Taguchi experimental design has been accomplished and the experimental response for each experiment has been counted, the signal-to-noise ratio (SNR), which measures the magnitude of the mean of a process compared to its variation, is used to assess the response performance according to the response characteristics. In general, for a response with “larger-the-better (LTB)” characteristic, the corresponding SNR could be calculated using equation (27), while for a response with “smaller-the-better (STB)” characteristic, the corresponding SNR could be calculated using equation (28) ³⁸ :

$$\mathrm{SN}{R}_{\mathrm{LTB}}=-10{\log }_{10}\left(\frac{1}{\mathrm{num}}\sum _{i=1}^{\mathrm{num}}\frac{1}{{y_i}^2}\right)$$

(27)

$$\mathrm{SN}{R}_{\mathrm{STB}}=-10{\log }_{10}\left(\frac{1}{\mathrm{num}}\sum _{i=1}^{\mathrm{num}}{y_i}^2\right)$$

(28)

where $y_i$ is the experimental result of response y at the ith measurement, $\mathrm{num}$ is the total number of measurements per experiment.

For multi-objective crashworthiness optimization of FGCST in the present study, two competing or conflicting objectives, namely SEA_θ and PCF⁰, are selected to measure and optimize the overall crashworthiness of FGCST under multi-angle crushing loading. Evidently, the SEA_θ is “larger-the-better (LTB)” to improve overall energy absorption, while PCF⁰ is “smaller-the-better (STB)” to control overall load severity. Hence, the SNR of SEA_θ (called SNR1) could be calculated using equation (27), while SNR of PCF⁰ (called SNR2) could be calculate using equation (28). The SNR1 for SEA_θ and SNR2 for PCF⁰ are calculated and demonstrated in Table 8.

Step 3: Grey relational analysis (GRA)

To overcome the hinder of traditional Taguchi method on handling multi-objective optimization, the GRA is integrated to transform the multiple responses to a single grey relational grade (GRG) that to be maximized. The specific steps are as follows:

• (1) Grey relational generation.

The experimental data sequence are first transformed to a uniform range of [0,1] through normalization, which is called grey relational generation. In general, the normalization methods vary according to response characteristics. For a response with “larger-the-better (LTB)” characteristic ³⁹ :

$$y_i\left(k\right)=\frac{x_i\left(k\right)-\min x_i\left(k\right)}{\max x_i\left(k\right)-\min x_i\left(k\right)}$$

(29)

Else, for a response with “smaller-the-better (STB)” characteristic:

$$y_i\left(k\right)=\frac{\max x_i\left(k\right)-x_i\left(k\right)}{\max x_i\left(k\right)-\min x_i\left(k\right)}$$

(30)

where $x_i\left(k\right)$ and $y_i\left(k\right)$ respectively denote the original and normalized value of kth response at ith experiment; $\max x_i\left(k\right)$ and $\min x_i\left(k\right)$ respectively denote the maximum and minimum value of $x_i\left(k\right)$ for the kth response among all experiments.

• (2) Grey relational coefficient (GRC)

Then, the GRC $r_i\left(k\right)$ is calculated from the normalized experimental data to measure the correlation between the current experiment and the desired solution as:

$$r_i\left(k\right)=\frac{\underset{i}{min}\underset{k}{min}|y_0\left(k\right)-y_i\left(k\right)|+\delta \underset{i}{max}\underset{k}{max}|y_0\left(k\right)-y_i\left(k\right)|}{|y_0\left(k\right)-y_i\left(k\right)|+\delta \underset{i}{max}\underset{k}{max}|y_0\left(k\right)-y_i\left(k\right)|}$$

(31)

where $y_0\left(k\right)$ is the ideal normalized value, that is, the maximum normalized SNR for the kth response among all experiments. $\delta$ denotes the identification coefficient varying in interval (0,1), δ = 0.5 is adopted in this paper as generally recommended.

• (3) Grey relational grade (GRG)

At last, the grey relational grade (GRG) measuring the correlation of current experiment and the ideal solution can be obtained by averaging the GRC of multiple responses as:

$${\epsilon }_i=\left(\frac{1}{Q}\right){\sum }_{k=1}^Q{r}_i\left(k\right)$$

(32)

Nevertheless, taking the difference of response preference or significance into account, different weights are assigned to different responses to calculate the weighted GRG:

$${\epsilon }_{i\omega }={\sum }_{k=1}^Q{\omega }_k{r}_i\left(k\right)$$

(33)

where ${\epsilon }_i$ and ${\epsilon }_{i\omega }$ represent the GRG and weighted GRG for ith experiment, ${\omega }_k$ represents the weight of the kth response and $Q$ represents the total number of responses.

In the present study, $Q$ =2 as SEA_θ and PCF⁰ are selected as two optimization responses, thus the GRG could be treated as a single comprehensive response representing the overall crashworthiness of FGCST. Larger value of ${\epsilon }_i$ or ${\epsilon }_{i\omega }$ signifies that the corresponding combination of factor levels is closer to the optimal condition and therefore is preferred.

Step 4: Main effect and interaction effect analyses

With the help of GRA, the multi-objective crashworthiness optimization of FGCST has now been transformed to a single-objective optimization of GRG, which then could be addressed by Taguchi analysis to explore the optimal combination of factor levels. To be more specific, the average response value of GRG for each factor level is first calculated, then the main effect value of each factor on GRG is obtained by computing the difference between the maximum and minimum average GRG of all factor levels, which could be described as follows:

$${\overline{f}}_p=\left(\frac{1}{U}\right){\sum }_{q=1}^U{\epsilon }_{\mathrm{pq}}$$

(34)

$$\Delta f=\max ({\overline{f}}_1,{\overline{f}}_2\dots ,{\overline{f}}_V)-\min ({\overline{f}}_1,{\overline{f}}_2\dots ,{\overline{f}}_V)$$

(35)

where ${\overline{f}}_p$ represents the average GRG value for factor $f$ at pth level, ${\epsilon }_{\mathrm{pq}}$ represents the GRG of $q$ th experiment with factor $f$ at $p$ th level, U denotes the total number of experiments with factor $f$ at pth level, V represents the total number of levels of factor $f$ and $\Delta f$ represents the main effect value of factor $f$ .

Commonly, the average response value of each factor level ( ${\overline{f}}_p$ ), the main effect value of each factor on response ( $\Delta f$ ) 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 then 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 a response. Parallel lines in an interaction plot indicates no interaction, and the greater the difference in slope between the lines, the higher the degree of interaction.

Results and discussion

Table 8 demonstrates the Taguchi experimental results, as well as calculated SNR of SEA_θ (SNR1) using equation (27) and SNR of PCF⁰ (SNR2) using equation (28) for each experiment. Table 9 demonstrates the response table for SNR1 of SEA_θ, while Table 10 demonstrates that for SNR2 of PCF⁰. For a certain factor, the mean of SNR at a certain level stands for the overall average value of SNR for all experiments with this specific level, while the Delta which is the maximum difference among different levels, represents the main effect value of this factor on SNR. In general, a factor with a higher Delta value among the factors signifies a relatively larger influence on the response over others. Consequently, according to Table 9, the most influential factor on SNR1 (or SEA_θ) is the top circumscribed radius (R_t), followed by the normalized gradient exponent (m) and finally the number of (45°/−45°) layers (N). In contrast, according to Table 10, the most influential factor on SNR2 (or PCF⁰) is normalized gradient exponent (m), followed by the top circumscribed radius (R_t) and at last the number of (45°/−45°) layers (N). These could also be concluded from the main effect plots of factors on SNR1 of SEA_θ and SNR2 of PCF⁰, as shown in Figure 16, from which the rank of influence of factors on SNR of responses could be determined by comparing the slopes of the lines. In addition, the main effect plots could be employed to determine the optimal combination of factor levels for each response. Specifically, for maximizing SNR1 (namely maximizing SEA_θ), the optimal combination of factor levels is R1-m2-N3 (namely R_t = 27.5 mm, m = 0 (or n = 1) and N = 2), as boldly marked in Table 9. In contrast, for maximizing SNR2 (namely minimizing PCF⁰), the optimal combination of factor levels is R2-m1-N4 (namely R_t = 30 mm, m = −0.301 (or n = 0.5) and N = 3), as boldly marked in Table 10.

Table 9.

Response table for SNR1 of SEA_θ.

Means	Level	Rt	m	N
SNR1 (LTB)	1	23.39	23.15	22.90
	2	23.37	23.40	23.25
	3	23.09	23.09	23.30
	4	22.72	22.92	23.10
	Delta	0.67	0.48	0.40
	Rank	1	2	3
	Total mean	23.14

Table 10.

Response table for SNR2 of PCF⁰.

Means	Level	Rt	m	N
SNR2 (LTB)	1	−26.33	−22.83	−26.71
	2	−25.78	−27.97	−26.69
	3	−27.24	−28.06	−26.72
	4	−27.03	−27.51	−26.25
	Delta	1.46	5.23	0.47
	Rank	2	1	3
	Total mean	−26.59

Figure 16.

Main effect plots of factors on SNR1 of SEA_θ and SNR2 of PCF⁰: (a) mean of SNR1 and (b) mean of SNR2.

Figures 17 and 18 demonstrate the interaction effect plots of design factors on SNR1 of SEA_θ and SNR2 of PCF⁰, respectively, based on which whether the interactions among design factors exert any influence on each response could be explored. 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 used to assess the two-way interactions among factors on the responses. Specifically, if the lines are parallel then there is no interaction, otherwise there is, and the greater the lines deviate from being parallel, the larger the interaction. Obviously, no matter on SNR1 or SNR2, interaction effects are observed between either two of the three factors (R_t, m, and N), and they differ from each other.

Figure 17.

Interaction effect plots of factors on SNR1.

Figure 18.

Interaction effect plots of factors on SNR2.

As disclosed above, the optimal combination of factor levels for maximizing SEA_θ (R1-m2-N3) entirely differ with that for minimizing PCF⁰ (R2-m1-N4), while the FGCST is desired to be manufactured with a specific level of design parameters to obtain optimal SEA_θ and PCF⁰ simultaneously, which can not be addressed by traditional Taguchi method limited to optimize a single response at a time. Thus, to promote the multi-objective crashworthiness optimization of FGCST, the grey relational analysis (GRA) is employed herein to convert the two competing responses (SEA_θ and PCF⁰) to a single one, that is, grey relational grade (GRG) that to be maximized (LTB). For this purpose, the gray relational generation for SEA_θ using equation (29) and PCF⁰ using equation (30), the grey relational coefficient (GRC) using equation (31) and grey relational grade (GRG) using equations (32) or (33) are obtained and demonstrated in Table 11. Likewise, the response table for GRG is calculated and demonstrated in Table 12, from which it can been seen the most influential factor on GRG is normalized gradient exponent (m), followed by top circumscribed radius (R_t) and finally the number of (45°/−45°) layers. In addition, the main effect and interaction effect plots for GRG are shown in Figures 19 and 20, respectively. It can be seen from Figure 19 that the optimal combination of factor levels for GRG is R2-m1-N2 (namely R_t= 30 mm, m = −0.301(n = 0.5), N = 1), as boldly marked in Tables 11 and 12. Besides, strong interaction effect on GRG is similarly observed among either two of the three factors.

Table 11.

Grey relational coefficient (GRC) and grey relational grade (GRG).

No.	Experiment levels			Grey relational generation		Grey relational coefficient (GRC)		Grey relational grade (GRG) (LTB)	Rank
	Rt	m	N	SNR1 (LTB)	SNR2 (LTB)	SNR1 (LTB)	SNR2 (LTB)
Weight				0.500	0.500	0.500	0.500
Ideal	/	/	/	1.000	1.000	1.000	1.000
1	1	1	1	0.800	0.700	0.715	0.625	0.670	4
2	1	2	2	0.309	1	0.418	1	0.709	3
3	1	3	3	0.188	0.591	0.381	0.550	0.465	12
4	1	4	4	0.328	0.755	0.426	0.671	0.549	7
5	2	1	2	1	0.825	1	0.740	0.870	1
6	2	2	1	0.216	0.674	0.389	0.605	0.497	10
7	2	3	4	0.193	0.868	0.382	0.791	0.587	6
8	2	4	3	0.498	0.625	0.499	0.571	0.535	8
9	3	1	3	0.707	0.915	0.631	0.854	0.742	2
10	3	2	4	0.229	0.711	0.393	0.634	0.513	9
11	3	3	1	0.223	0.311	0.391	0.420	0.406	14
12	3	4	2	0	0.355	0.333	0.436	0.385	15
13	4	1	4	0.915	0	0.855	0.333	0.594	5
14	4	2	3	0.028	0.700	0.339	0.625	0.482	11
15	4	3	2	0.132	0.523	0.365	0.511	0.438	13
16	4	4	1	0.192	0.149	0.382	0.370	0.376	16

Table 12.

Response table for GRG.

Means	Level	Rt	m	N
GRG (LTB)	1	0.5988	0.7195	0.4876
	2	0.6227	0.5510	0.6010
	3	0.5120	0.4745	0.5567
	4	0.4730	0.4616	0.5611
	Delta	0.1497	0.2579	0.1134
	Rank	2	1	3
	Total mean	0.5516

Figure 19.

Main effect plots of factors on GRG.

Figure 20.

Interaction effect plots of factors on GRG.

Verification of optimization results

Since the optimal combination of factor levels, as well as corresponding structural parameters of FGCST for better overall crashworthiness (SEA_θ and IPCF⁰) is obtained, the next step is to verify the effectiveness of the optimization results. Specifically, a detailed crashworthiness comparison among the optimized design (R2-m1-N2, namely R_t= 30 mm, m = −0.301 (or n = 0.5), N = 1) and the baseline design (R2-m3-N1, namely R_t= 30 mm, m = 0.301 (or n = 2), N = 0) of FGCST under multi-angle crushing loading is conducted.

Table 13 demonstrates the comparison of crashworthiness indicators before and after optimization, including SEA and PCF under individual angle crushing loading (0°, 10°, 20°, 30°), and the overall SEA (SEA_θ) and PCF(PCF⁰) under multi-angle crushing loading. Evidently, the SEAs under individual angle crushing loading have been improved by different percentages, among which the largest is 32.01% (30°), followed successively by 2.78% (0°), 2.17% (20°) and 1.62% (10°). Meanwhile, the PCFs under individual angle crushing loading present various change, with decrease by 40.58% (0°) and 5.05% (10°), as well as increase by 1.72% (20°) and 24.25% (30°). These changes under individual angle crushing loading could not only be more clearly observed according to Figure 21, but also be reflected by the comparison of final deformation modes in Table 14. For instance, it can be seen that under 30° angle crushing loading, the optimized design presents local bending accompanied by progress collapse on the upper part of the structure, which is more conducive to energy absorption compared with the baseline design with overall bending accompanied by four-corner cracking. Except for this, there is minor difference in final deformation modes between optimized design and baseline design under other angle crushing loading (0°, 10°, 20°), just as reflected by comparison of crashworthiness indicators in Table 13 and Figure 21. Nevertheless, taking multi-angle crushing loading (0°∼30°) into account, the overall SEA (SEA_θ) is increased by 6.06%, meanwhile the overall PCF (PCF⁰) is decreased by 40.58%, indicating that the multi-objective crashworthiness optimization for FGCST is correct and effective.

Table 13.

Comparison of crashworthiness indicators before and after optimization.

Crushing scenario	Crushing angle (°)	Original	Optimized	Change (%)	Original	Optimized	Change (%)
Individual angle		SEA			PCF
	0	17.748	18.242	2.78	20.48	12.17	−40.58
	10	17.022	17.297	1.62	12.48	11.85	−5.05
	20	14.309	14.620	2.17	10.49	10.67	1.72
	30	7.298	9.634	32.01	5.03	8.25	24.25
		SEAθ			PCF0
Multi-angle	0–30	14.094	14.948	6.06	20.48	12.17	−40.58

Figure 21.

Comparison of crashworthiness indicators before and after optimization: (a) SEA and (b) PCF.

Table 14.

Comparison of final deformation modes before and after optimization.

Designs	Initial states	Final deformation modes
		0°	10°	20°	30°
Baseline design
Optimized design

To sum up, the optimized FGCST presents more superior crashworthiness characteristics under multi-angle crushing loading compared with the baseline (or original) counterpart, thus demonstrating good potential as a superior candidate for energy absorbers especially under multi-angle crushing loading.

Conclusions

The present study proposes a novel CFRP functionally-graded conical square tube (FGCST), for which a comprehensive crashworthiness comparison, a detailed parameter study and a multi-objective crashworthiness optimization under multi-angle crushing loading are carried out in succession. Main findings are summarized as follows:

• (1) The constructed finite element (FE) models of FGCST are of good accuracy with simulation results, namely crushing force versus displacement and final deformation modes, basically consistent with that of experiments. Therefore, the subsequent crashworthiness comparison and multi-objective optimization for FGCST based on FE models are reasonable and credible.

• (2) Compared with conventional CST and SST, the novel FGCST demonstrates relative merits on crashworthiness. In particular, the FGCST with arc-shape (n < 1) exhibits better crashworthiness with lower PCF and similar SEA, demonstrating good potential as superior candidates for energy absorbers.

• (3) According to parameter study, both the EA and PCF under any specific crushing loading angle approximately decrease with increasing Δ $R$ . Meanwhile, on the whole, the FGCST with larger value of number of (45°/−45°) layers (N) tends to present relatively higher EA and PCF under crushing loading angles of 0° and 30°, while it tends to present relatively lower EA and PCF under crushing loading angles of 10° and 20°, indicating the effect of number of (45°/−45°) layers (N) on crashworthiness varies with crushing loading angles.

• (4) Compared with the baseline (or original) design, the optimized FGCST obtains 6.06% higher of SEA_θ and 40.58% lower of PCF⁰ simultaneously, thus demonstrating good potential as a superior candidate for energy absorbers especially under multi-angle crushing loading. In addition, Taguchi method coupled with GRA could be an effective approach for handling similar multi-objective optimization issues.