Multi-objective crashworthiness design optimization of a rollover protective structure by an improved constraint-handling technique

1.   Introduction

Rollover protective structure (ROPS) is a frame or a cab installed on tractors or construction machines. ROPS can absorb a certain proportion of impact energy (generated by the weight of the machine) to protect the operator wearing a seat belt from being injured in a rollover accident. For a ROPS certified by ISO 12117-2:2008 ^[1], on the one hand, it should be stiff enough to protect the deflection-limiting volume (DLV, a survival or safety space to mimic the male operator wearing a protective helmet ^[2]) from the deformed components. On the other hand, it should absorb/dissipate enough impact energy so as to avoid further rolls which will cause additional damage to the operators. Commonly, these two goals are conflicting ^[3]. Therefore, determining an appropriate design for the ROPS to achieve these two goals simultaneously is extremely important.

The crashworthiness of the protective structure is easily influenced by the design methods. One way to improve the crashworthiness of ROPS is to equip the ROPS with additional energy absorbing devices ^[4,5]. The devices contact the ground ahead of the ROPS when the machine rolls over and absorb a certain amount of impact energy by plastic deformation. Thus, the total amount of energy absorption is increased. However, these additional devices either affect the aesthetics of the structure or increase the complexity and price of the structure.

Another way to improve the crashworthiness of ROPS is to use structural optimization design methods to find an appropriate structural configuration ^[6,7,8,9]. However, to our best knowledge, few studies focus on the crashworthiness optimization of ROPS. Therefore, the well-studied crashworthiness optimization of the automotive bodies is used as an example here since most of them are space frame structures too ^{[10,11,12,13]}. The crashworthiness of a structure is greatly affected by its structural configuration ^[14,15,16]. Thus, before the optimization process, the relationship between the crashworthiness and the structural configuration should be well studied by the structural response analysis ^[17,18,19]. Ghadianlou and Abdullah ^[20] use the material type and geometry of the door components as design variables to study the effects of different structural configurations on the door deformation and crushing force during a side impact. Zhu et al. ^[21] use extensive numerical simulations to study the relationship between structural responses (occupant safety space and peak acceleration) and design variables during a frontal crash. Then, the surrogate model-based crashworthiness optimization is performed to design a lighter body-in-white. Yu et al. ^[22] analyze the influence of the front end components made out of uniform thickness and variable thickness sheet-metals on the structural response (acceleration and energy absorption) of pure electric vehicles in the front collision scenario. A multi-objective optimization (MOO) method is utilized to determine the thickness distribution of the structure. In general, for ROPS, the second way is more practical.

As mentioned above, ROPS should find a trade-off between stiffness and energy absorption and fulfill certain performance constraints set by regulations. Thus, the crashworthiness optimization of a ROPS belongs to constrained MOO problem. The commonly used constraint handling method is the penalty function method (PFM). PFM can be categorized into static PFM ^[23], dynamic PFM ^[24] and the adaptive PFM ^[25] based on the determination of penalty weights and the main drawback of PFM is the requirement of fine tuning of the penalty weights. The design space will be narrowed if inappropriate penalty weights are selected, resulting poor solutions.

This study aims to conduct crashworthiness optimization of a ROPS by combining the structural response analysis and the MOO with an improved adaptive PFM. First, an experimental test under the ISO 12117 criteria is conducted and is used to verify the developed finite element (FE) model of the ROPS. Then, the energy absorption capacity and deformation behaviors of the components in the ROPS are analyzed based on the amount of energy absorption and the cross-sectional forces. Subsequently, components that are sensitive to energy absorption or weight reduction are identified and the thicknesses of these components are considered as the design variables. A constrained MOO problem, which aims to maximize the safety distance between DLV and the ROPS and minimize the total mass without sacrificing the energy absorption capacity, is formulated. To improve the convergence rate and obtain optimal solutions with good diversity of distribution and fully explore the design space, a constrained MOO algorithm is proposed based on an improved adaptive PFM on the basis of the non-dominated sorting genetic algorithm (NSGA-Ⅱ).

The remainder of the paper is organized as follows: the experimental ROPS test, the FE modeling, the validation of the FE model by the experimental test and the structural response analysis are introduced in Section 2. The proposed constrained MOO algorithm is introduced and is validated by benchmark problems in Section 3. The crashworthiness optimization of the ROPS is presented in Section 4. Section 5 gives a conclusion.

2.   Experimental test and finite element modeling

The ROPS, installed on a hydraulic excavator (Figure 1(a)), is studied based on ISO 12117. ISO 12117 defines a quasi-static laboratory test for ROPS which includes a sequence of three consecutive loads: 1) a lateral load applied from the side which places the severest requirements on the ROPS (left side in this study), 2) a horizontal rear (longitudinal) load, 3) a vertical load (Figure 1(b)). The ROPS passes the test if it tolerates a predefined force in the side loading test and the vertical loading test and absorbs particular energy in the side loading test and the longitudinal loading test. Moreover, during the test, the DLV should not be violated by the ROPS. To accurately predict the deformation behaviors of the ROPS during the laboratory test by numerical simulation, since the above three loading tests are carried out consecutively, one at a time, the accumulated plastic deformation by the former loading test should not be ignored ^[26,27,28]. However, researches ^[5,26] show that if the ROPS can pass the side loading test, it can also pass the longitudinal and the vertical loading tests in most cases. That is, the side loading test is the most critical loading test. Thus, in this study, to simplify the numerical simulation and reduce the computational cost of the crashworthiness optimization, the numerical model only simulates the side loading test.

2.1.   Experimental test

The laboratory test was carried out at an authorized station for testing ROPS. The ROPS is fastened to the fixture which is totally constrained on the testbed in the experimental test. The side load is applied by a hydraulic actuator at a displacement rate that less than 5 mm/s. To avoid the local penetration at the load application point, a load distribution device (LDD, made up of thick steel plate) is welded on the horizontal part of the A pillar as is shown in Figure 2. As introduced earlier, besides no infringement to DLV is compulsory, the ROPS has to satisfy the resistance force and energy criteria in the side loading test. The lower bounds of these two criteria are defined by Eq (1) ^[1].

where M is the operating mass of the excavator, 9000 kg in this study. F_s and E_s denote the minimum resistance force and absorbed energy that have to be satisfied, respectively. That is, the ROPS in question has to tolerate a resistance force of 30.8 kN and absorb an impact energy of 11,395.8 J at least.

2.2.   Finite element modeling

In this study, the numerical simulation/virtual ROPS test is conducted using the explicit FE software LS-DYNA. The ROPS in question is mainly composed of panels and thin-walled beams. Thus, components are all discretized by Belytschko-Tsay (B-T) shell finite elements, 67,796 four-node elements and 234 three-node elements. The fixture is also discretized by the B-T shell finite elements, 12,485 four-node shell finite elements. The element size is chosen as 10 mm. The B-T elements use the reduced integration algorithm which can significantly reduce the computational cost ^[29,30]. The disadvantage of this type of element is that it would result in zero-energy mode. To avoid this issue, the hourglass control is necessary. It can be considered that the numerical results are not affected by the hourglass control when the ratio between the hourglass energy and total energy is less than 5% ^[26,31]. The DLV is molded in null beam elements (*ELEMENT_PLOTEL) to monitor the deformations of the ROPS.

Two grade steel materials are adopted in the ROPS, DC04 and Q345. The roof and right panel are made out of DC04 and the rest components are made out of Q345. The piecewise linear plasticity model is implemented. The densities, Young's moduli and the Poisson's ratios of DC04 and Q345 are 7900 kg/m³, 210 GPa and 0.3, respectively. The yield stresses of DC04 and Q345 are of 180 MPa and 345 MPa, respectively ^[26]. Due to the inevitable plastic deformation during the ROPS test, to obtain a more accurate numerical simulation result, the detailed stress-strain relationship is necessary. In this study, the specimens taken from components of the ROPS are tested based on ASTM E8:2011 ^[32] to obtain the engineering stress-strain curves (Figure 3).

The degrees of freedom of nodes on the bottom face of the fixture are totally constrained which is consistent with the experimental test. A virtual hydraulic actuator modeled by *RIGIDWALL_GEOMETRIC_FLAT_MOTION is used to push the LDD. ISO 12117 defines that in the quasi-static test the displacement rate should less than 5 mm/s ^[1]. However, the computational cost would be unaffordable if the numerical model adopts such a speed. Researches ^[26,30,33] show that if the amount of kinetic energy only accounts for a small portion of the total energy throughout the virtual ROPS test, it can be considered that improving the displacement rate will not affect the energy numerical result. Thus, the virtual hydraulic actuator moved at a relatively high displacement rate, 1 mm/ms and continued 220 ms in this study. After 220 ms, the load was held for 5 ms. With this improvement and a large outputs' time interval, the simulation can be done in 1.5 h (Intel Core i7 2.3GHz, 16GB RAM). Due to the significantly reduced computational cost, the surrogate model-based crashworthiness optimization which needs extensive FE simulations will be more practical. The automatic surface to surface contact algorithm is utilized for modeling the contact between each pair of connected components. Both the static and dynamic coefficients of friction are set 0.15 ^[26,34].

2.3.   Results and discussion

First, the reliability of the numerical model is verified by the energy conservation. The energy history of the numerical simulation is shown in Figure 4. It is observed that the summation of the kinetic energy, internal energy, hourglass energy, damping energy and sliding energy equals to the total energy. Thus, the law of conservation of energy is satisfied. In addition, during the virtual ROPS test, the maximal proportions of kinetic energy and hourglass energy to the total energy are 0.143 and 0.397%, respectively. Thus, the great increase in displacement rate and the utilization of B-T element have little impact on energy. At 220 ms, the value of absorbed energy is 11,800.5 J, i.e., the energy criterion is satisfied.

The comparison of resistance force histories between the numerical simulation and the experimental test is shown in Figure 5(a). It is noted that the numerical values are in good agreement with the experimental values. In addition, at 20 ms, both the virtual ROPS and the real ROPS meet the force criterion (>30.8 kN). That is, the force criterion is satisfied much earlier than the energy criterion. Therefore, compared to the energy criterion, the force criterion is easier to be satisfied, which is consistent with engineering experience. The external work history by the applied load can be calculated based on Eq (2). The maximal difference between the numerical result and the experimental result during the test is of 0.53%, which is negligible (Figure 5(b)).

where F denotes the external force applied on the ROPS and u denotes the displacement load.

The distance (d_s shown in Figure 6(a)) in x-axis between DLV and a specified node on the door post which has same z-axis coordinate with the shoulder of DLV is monitored. The minimal value is of 20.86 mm when the energy criterion is satisfied. Thus, during the virtual ROPS test, DLV is not violated by the ROPS which is consistent with the experimental test. Limited by the test conditions, d_s in the experimental test didn't be monitored. However, the deformation behaviors during the virtual ROPS test agree well with those in the experimental ROPS test (Figure 6(b, c)). For more details about the experimental test, the reader is referred to ^[26].

The above comprehensive comparisons of structural responses for numerical simulation and the experimental test demonstrate that the ROPS meets the requirements set by ISO 12117 and pass the ROPS test. Moreover, it can be considered that the parameter settings for the FE model are reasonable. The accuracy of the numerical simulation is high enough to predict the deformation behavior and performance and the FE model can be used to replace the expensive and time-consuming experimental test.

It can be noticed from the aforementioned discussion the ROPS meets the energy and force criteria. Hereinafter, a further structural response analysis on the energy absorption and cross-sectional force of each component is carried out to assess its effect on the crashworthiness of the ROPS.

Figure 7(a) shows the amounts of energy absorption of the main energy absorbing components and where they locate in the ROPS (Figure 7(b)). It is noted that the main energy absorbing components are floor, posts and roof which suffer large local plastic deformation. The floor, for example, its energy absorption accounts for 18.38% of the total which should be attributed to the plastic deformation around the fixing holes (enlarged view in Figure 8) though the thickness of the floor is of 8 mm. The energy absorptions of the left C pillar, front panel, left A pillar, right A pillar and the roof are 14.62, 8.23, 8.03, 7.27 and 6.09% of the total, respectively. Moreover, it can be observed from Figure 8 that the areas suffer plastic deformation are only a small part of the corresponding components.

Figure 9(a) shows the values of cross-sectional force (values at 220 ms) of the main load-carrying components along x-axis and where they locate in the ROPS (Figure 9(b)). From the cross-sectional force analysis, it can be seen that the main load-carrying components are rear plates and the rear lower beam. The deformation-resistant substructure composed of the rear plates, the rear lower beam and the posts (including the doorpost and the left and right C pillar) can effectively resist to the side impact force. The ratio of the cross-sectional force (a compressive load) of the rear plates to the corresponding external load is 0.574. The components that tolerate the compressive load also include the rear lower beam, roof reinforced plate, upper front reinforced plate and the front panel. The ratios of the cross-sectional forces of these components to the external load are 0.408, 0.261, 0.152, 0.062, respectively. Meanwhile, the rear middle and rear upper beams tolerate tensile loads. The upper part of the right C pillar and the doorpost also tolerate tensile force (z-axis). The load path distribution can be seen in Figure 9(b).

Combined with the analysis on energy absorption, it can be seen that the box-shaped components located along the x-axis mainly subject to axial loads while suffer to few plastic deformations. Thus, the amounts of the energy absorption of these box-shaped components are also very small. The roof and front panel have large energy absorption due to the obvious plastic deformation. The roof, for example, cannot withstand the external load after 30 ms while the energy absorption is significantly increased due to large plastic deformation. However, studies ^[35,36] have shown that premature plastic deformation does not achieve the desirable energy absorption capacity, that is, the energy absorption capacity is not fully exploited. The premature plastic deformation can be avoided by increasing the stiffness of the corresponding components (e.g., increasing the thicknesses of the components) and is beneficial to increase the absorbed energy. Numerical experiment shows that when the thickness of the front panel is set to 5 mm (initial value 3 mm), the energy absorption of the front panel is increased by 33.94%. However, the mass is increased by 40%, which is not conducive to reduce the total mass of the ROPS. Because the energy absorption is mainly contributed by the plastic deformation of the component and the plastic deformation area only occupies a small part of it. Therefore, to improve the energy absorption capacity and avoid premature plastic deformation, this study uses the tailor welded blanks (TWB) technique to make several components so that components with variable thickness can be obtained which improves the material utilization in the workshop and reduces the total mass of the ROPS.

3.   Methodology

3.1.   NSGA-Ⅱ algorithm

The crashworthiness optimization problem, in this study, is a constrained MOO problem. The response values are results of numerical simulations. The heuristic algorithms are usually used for such a black-box MOO problem ^{[10,37,38,39]}. NSGA-Ⅱ is one of the most popular multi-objective heuristic algorithms. The main characteristics of NSGA-Ⅱ are the non-dominated sorting scheme and the crowding distance estimation. These two characteristics are used to improve the convergence of genetic algorithm and maintain the diversity of candidate solutions, respectively. Thus, NSGA-Ⅱ is used in this study and the optimization parameters are presented in Table 1.

The steps of NSGA-Ⅱ can be summarized as follows:

Step 1. Initialize the parent population randomly.

Step 2. Generate offspring population through genetic operations, i.e., crossover and mutation.

Step 3. Evaluate all objectives.

Step 4. Rank all individuals based on the non-dominated sorting and the crowding distances.

Step 5. Choose next parent population by environmental selection.

Step 6. If the maximum iteration is achieved, the optimization should be terminated, or return to Step 2.

3.2.   Modified optimization algorithm

Compared to the unconstrained MOO problem, constrained MOO problem has to balance the satisfaction of constraints and the optimization of objective functions ^{[40,41,42,43]}. Thus, it is more challenging to solve the constrained MOO problems than the unconstrained MOO problems.

Without loss of generality, the constrained MOO problem can be formulated by:

where x is the decision vector represents a candidate solution. x^L and x^U are the corresponding upper bound and lower bound of x. f_i(x) is the ith objective function. m, p, q denote the number of objective functions, inequality constraint and all constraint functions, respectively. g_j(x) and h_j(x) denote the jth inequality constraint and jth equality constraint functions, respectively. To handle the constrained MOO problem, objectives are usually normalized ^[41], i.e.,

where $f_i^{\max }$ and $f_i^{\min }$ denote the maximal and minimal values of the ith objective. $f_i^{\text{N}}({\bf{x}})$ is the normalized ${f_i}({\bf{x}})$. Then, the constraints have to be normalized too:

where δ is a very small positive value to relax the equality constraints into inequality one since the equality constraint is too strict to be satisfied, δ = 10^-6 in this study. c_j(x) is the constraint violation of jth constraint of x. cmax j is the maximal constraint violation of jth constraints. cN j(x) is the normalized c_j(x). v(x) is the average of the constraint violations of x.

Finally, v(x) is added to the objective as a penalty term (Eq (9)). That is, the constrained MOO problem is transformed into an unconstrained MOO problem and the common NSGA-Ⅱ can be used as the optimizer ^[41].

where $f_i^{{\text{mod}}}({\bf{x}})$ denotes the ith modified objective. γ_f, the so-called penalty weight, is a coefficient equals the ratio of number of feasible solutions to that of the population size ^[41]. However, by this method ^[41], all the solutions in current population have the same penalty weight. If the penalty weight is high, due to the preference of feasible solution, the feasible solutions close to the constraint boundary will not be fully explored ^[40]. While, if the penalty weight is small, the difference between the feasible solution and infeasible solution close to the constraint boundary may be negligible. Thus, more sophisticated non-dominated sorting is needed, leading to a low convergence rate. To remedy this issue, a feasible-guiding strategy was proposed ^[41] to fully explore solutions close to the constraint boundary and improve the convergence rate on constrained MOO problems.

In this study, an improved adaptive PFM is proposed. Solutions will be assigned different penalty weights based on the amounts of constraint violation. The modified objective function is presented in Eq (10):

where v_j(x) is the jth constraint violation. d_j(x) denotes the value of jth constraint function of x. [d_j(x)] is the corresponding upper bound. γ_j is the penalty weight for the jth constraint of x and can be defined as:

where τ_m is a small positive value, τ_m = 0.01 in this study. The relationship between the penalty weight and the constraint violation can be illustrated in Figure 10. It is noted that the greater the constraint violation, the greater the penalty weight assigned to the constraint and vice versa. When the constraint is inactive, the penalty weight is negative, i.e., the objective function can be further penalized. Thus, the feasible solutions have a higher probability to be retained to next generation and infeasible solutions with large constraint violations are hard to be retained. In addition, for the solutions close to the constraint boundary, the penalty weight is small which is beneficial to save the potential optimal solutions so as to improve the convergence rate and diversity of distribution of optimal solutions.

3.3.   Feasibility analysis of optimization algorithm

To verify the effectiveness of the proposed MOO algorithm, two benchmark problems, the Srinivas (SRN) test problem (Srinivas and Deb ^[44], Eq (12)) and the Tanaka (TNK) test problem (Tanaka et al. ^[45], Eq (13)) are solved by the proposed method and method proposed by ^[41], respectively.

The Pareto optimal fronts are presented in Figure 11. It is noted that the Pareto fronts obtained by the proposed method (indicated by nonlinear weight in Figure 11) is more superior on diversity of distribution than those obtained by the method of ^[41] (indicated by same weight in Figure 11). Solutions obtained by the proposed MOO algorithm converge to the true Pareto optimal front and several solutions obtained by the method of ^[41] violate the constraints. It is clear that the optimal solutions obtained by the proposed method is also in good agreement with the true Pareto optimal front though the TNK that has a disconnected true Pareto optimal front. Again, infeasible solutions are obtained by the method of ^[41].

To quantitatively compare the effectiveness of the proposed method, the inverted generational distance (IGD) ^[46] is utilized as the performance metric to assess the convergence and diversity of population of the proposed and the method of ^[32]. IGD is a dimensionless value represents the average distance from solutions in the true Pareto optimal front to the nearest solution in the Pareto front obtained by a certain algorithm. That is, the smaller the IGD, the better the convergence and diversity in population of the algorithm. IGD can be defined as:

where Pa_R and Pa denote the true Pareto optimal front and the Pareto optimal front obtained by a certain algorithm, respectively. p_j is the jth individual in Pa_R. d(p_j, Pa) is the minimal Euclidean distance between p_j and Pa. |Pa_R| is the number of solutions in the true Pareto optimal front.

To reduce the computational errors as much as possible, the above two benchmark problems are solved in 50 separate runs by the proposed method and the method of ^[41], respectively, and the average of the IGDs in 50 runs are considered as the final IGD. Table 2 presents the IGDs. Note that the IGDs obtained by the proposed method are much smaller than those obtained by the method of ^[41]. Thus, in summary, the optimal solutions obtained by the proposed method show superiority on convergence rate and diversity of distribution.

4.   Multi-objective crashworthiness optimization for the ROPS

4.1.   Multi-objective optimization process

In this section, the crashworthiness optimization problem for the ROPS is solved by the finite element analysis (FEA)-based surrogate model. First, the design variables and output responses are determined. The output responses are considered as objective or constraint functions based on the optimization model. Then, sampling points are generated by using the design of experiment method. FEA is utilized to calculate the output responses of the sampling points which can be used to construct surrogate models to approximate the responses. An error analysis is conducted to determine the most appropriate surrogate model. Finally, the optimal solutions are obtained by solving the optimization problem. The flowchart of the multi-objective crashworthiness optimization is described in Figure 12.

As mentioned above, energy absorption, the resistance force and d_s are set as criteria to determine whether the ROPS meets the performance requirements. However, the force criterion is easily to be satisfied. Thus, in this study, energy absorption and d_s are selected as the output responses and the resistance force is not selected. In addition, the lower the total mass, the lower the manufacturing cost. Thus, the total mass of the ROPS will also be selected as the output response.

As revealed in Section 2, the floor panel, C pillar, A pillar, front panel and roof are the main contributors to the energy absorption of the ROPS. However, the thicknesses of the floor panel, C pillar and A pillar are custom made and are applied to many other ROPSs. New molds and adjustments to assembly process have to be made if their thicknesses are changed, leading to a significantly increase in manufacturing cost. Thus, thicknesses of these components will not be selected as the design variables. The thicknesses of the front panel and roof will be selected as the design variables and intended to be made by the TWB technique. These two components will be partitioned into two separate parts, respectively, i.e., the front panel_1, front panel_2, roof_1 and roof_2 which correspond to the design variables x₁~x₄, respectively. In addition, the thicknesses of the upper front reinforced plate (x₅), roof reinforced plate (x₆) and the rear plates (x₇~x₉) are insensitive to energy absorption but have large values which have great potentials to be reduced. Thus, they also be selected as the design variables (Figure 13). From the viewpoint of manufacturability, the range of each design variable is determined by commercially available steel thicknesses which is presented in Table 3.

4.2.   Definition of optimization problem

To improve the safety distance of the operator and reduce the total mass of the ROPS, d_s and the total mass of the ROPS (mass) are considered as the objective functions. Ideally, the value of d_s should be the minimum distance, along the x-axis, between the shoulder of DLV and the specified node on the doorpost when the ROPS meets the energy criterion, denoted as d_{s, R} in this study. However, different solutions meet energy criterion at different moments. Some solutions even failed to meet the energy criterion until the end of numerical simulations. Thus, approximating d_{s, R} accurately is difficult. Therefore, to construct a more accurate surrogate model of d_s and facilitate the comparisons between solutions, d_{s, 220} (d_s at 220 ms) is selected to construct the surrogate model. Meanwhile, the amount of energy absorption at 220 ms, E_s is considered as the constraint function, [E_s] (E_s of the baseline design) is set as the lower bound. The optimization problem is described as:

4.3.   Error analysis of multiple surrogate models

The optimal Latin hypercube sample technique is used to generate 60 sampling points in this study. FE models are established based on these sampling points and are analyzed. The true output responses are extracted and are used to construct the corresponding surrogate models.

The fitness to response depends on the surrogate model utilized. To choose an appropriate surrogate model, four types of surrogate model are utilized, i.e., the first and second order polynomial response surface models (denoted as RSM₁ and RSM₂, respectively), the radial basis function (RBF) neural network and the Kriging model. Two metrics, the R-squared (R²) and the root mean square error (RMSE), are used to evaluate the fitness of the surrogate models. In general, the closer the R² values are to 1, the more accurate the surrogate model. While, the closer the RMSE are to 0, the better the surrogate model. The definitions of R² and RMSE are presented in Eqs (16) and (17):

where n is the number of samples. y_i is the true output response. ${\hat y_i}$ denotes the corresponding predicted values and ${\bar y_i}$ denotes the average of all y_i. The comparison on R² and RMSE for the above four surrogate models are presented in Table 4. Note that, Kriging model provides the worst approximation on all three responses. RSM₂ is superior to RSM₁. With regard to the comparison between RSM₂ and RBF, RSM₂ has lower R² on d_s while the other two responses are fitted better than those by RBF. Thus, RSM₂ is chosen to approximate the responses.

4.4.   Optimization results and discussion

For a black-box constrained MOO problem, no true Pareto optimal fronts exists and can be used to evaluate the effectiveness of a certain algorithm. Thus, to demonstrate the effectiveness of the proposed constrained MOO method, the archive-based micro genetic algorithm (AMGA) and the multi-objective particle swarm optimization (MOPSO) are also used as the optimizer. The parameter settings are described in Table 5. Each algorithm is performed by 50 separate runs and the optimal solutions by these 50 runs are considered as the final Pareto front.

The Pareto fronts by the three algorithms are illustrated in Figure 14. It is noted that the Pareto fronts obtained by the NSGA-Ⅱ and AMGA are in the lower right position, which indicates that solutions obtained by these two algorithms have higher safety distance in the same mass or have lighter mass in the same safety distance. Moreover, the solutions obtained by NSGA-Ⅱ are evenly distributed than those obtained by the AMGA. To quantitatively compare these three algorithms, the minimum distance selection method ^[47] is used to select the compromise solutions (indicated by the red solid marks in Figure 14).

The comparison on performances and configuration among the compromise solutions are presented in Table 6. Note that, all these solutions obtained variable thickness front panel and roof. The determined values of x₁ and x₃ are all higher than those of x₂ and x₄. This is because that the plastic deformation mainly occurs on the front panel_1 (x₁) and roof_1 (x₃). Therefore, the front panel and roof improved their material utilization by optimizing the distribution of thicknesses of plates using the TWB technique. Moreover, it is clear that the compromise solutions obtained by NSGA-Ⅱ and AMGA, solution N and solution A for short, achieve a significant weight reduction. However, the amount of energy absorption of the solution A almost activates the constraint. From the comparison between the optimal values of solutions by NSGA-Ⅱ and MOPSO, the compromise solution by MOPSO, solution M for short, is not only 1.62 kg higher in mass, but it is also 102.78 J lower in energy absorption.

The true responses are calculated by FEA and solution N has the highest value of energy absorption, which is 621.75 and 290.04 J higher than those of solution A and solution M, respectively. d_{s, 220} of solution N is the lowest among three compromise solutions. However, the moments that solution N, solution A and solution M meet the energy criterion are 211,220 and 215 ms, respectively and the corresponding d_{s, R}s are 26.58, 20.93 and 24.04 mm, respectively. Therefore, solution N has the largest occupant safety distance. The resistance force and energy absorption histories of the three solutions are illustrated in Figure 15(a, b). It can be seen from Figure 15(a) that all these solutions meet the force criterion. During the virtual ROPS test, most of the time, the resistance force of solution N is larger than those of solution A and solution M which indicates a better ability to resist to the external load. It can be seen from Figure 15(b) that although all three solutions satisfy the energy criterion, the amount of energy absorption of solution N is always higher than those of solution A and solution M during the simulation process.

After comparing the performance indicators, i.e., the resistance force, energy absorption, occupant safety distance and the total mass, solution N is chosen as the optimal solution of the multi-objective crashworthiness optimization. Compared with the baseline design, the total mass of the optimal solution is reduced by 12.39 kg (7.06% of the baseline design). The optimal solution absorbs 2.45% more energy than that of the baseline design at 220 ms. The d_{s, R} is 26.58 mm, which is 27.42% higher than that of the baseline design (20.86 mm). Therefore, the optimal solution improves occupant safety distance, absorbs more impact energy and reduces the total mass simultaneously. That is, the crashworthiness and the total mass of the ROPS are optimized simultaneously.

The variations of mass and amount of energy absorption of the optimized components before and after the crashworthiness optimization are presented in Table 7. It can be seen from Table 7, except the front panel, the masses of the rest components are all decreased, especially the rear plates which are insensitive to energy absorption. The mass of these three rear plates is reduced by 8.87 kg, contributing 71.59% of the total weight reduction. After the crashworthiness optimization, the amounts of energy absorption of most components have been increased (except for the upper front reinforced plate and the rear plate_3). The amount of energy absorption of rear plate_1 has been increased 209.65% after reducing its thickness and the amounts of energy absorption of the front panel and roof have been increased 18.24 and 37.68%, respectively, which indicates that the baseline design has unreasonable configuration.

The comparison of plastic strain is plotted in Figure 16. Noticeably, new plastic deformation appeared in the rear plate as is shown in Figure 16(b) because its thickness is reduced. The maximal plastic strain is 0.293 (appeared in the right joint which connects the right C pillar and rear upper beam), which is less than the rapture strain of Q345, 0.33 ^[48]. Thus, the ROPS can resist to the external load and is stiff enough to protect the operator.

Due to reduce the manufacturing cost, nine design variables are selected. In the future, more design variables are expected to be considered so that a better design of ROPS can be obtained.

5.   Conclusions

In this study, the deformation behavior and energy absorption capacity of a ROPS are studied by experiment and numerical methods, and a multi-objective crashworthiness optimization strategy for ROPS is proposed. The main contributions and conclusions are described as follows:

1) The main components that have a significant influence on crashworthiness are identified. The main energy absorbing components that suffer large plastic deformation include posts, panel-like components between posts. The thicknesses of the main energy absorbing components have a significant influence on energy absorption capacity. The deformation-resistant substructure composed of the doorpost, C pillar, rear plates and the corresponding box-shaped beams among them contributes few to the energy absorption and resists the side load to maintain a survival space around the operator.

2) The proposed MOO algorithm is developed based on the modified objective function method and the improved adaptive PFM. The mathematical benchmark problems and the real-world engineering problem demonstrate the proposed MOO algorithm is superior to existing methods in terms of convergence rate and diversity of distribution.

3) The multi-objective crashworthiness optimization results show that the safety distance can be increased and the total mass can be reduced without sacrificing the energy absorption capacity of the ROPS. Compared with the baseline design, a more reasonable configuration is determined. It is noted that rear plate_1 and roof absorb much more impact energy when their thicknesses are reduced. The safety distance of the optimal solution by the proposed method is increased by 27.42% and the total mass is reduced by 7.06%. The method proposed in this study can be used as an alternative design method for the design of other protective structures in construction machinery and crash-resistant structures in the automotive industry.

Acknowledgements

This work was supported by the Basic General Scientific Research Program of Higher Education of Jiangsu Province [22KJD460003], the Science Foundation of Jiangsu Vocational Institute of Architectural Technology [JYJBZX22-05] and the National Natural Science Foundation of China [52805123, 52165029].

Conflict of interest

The authors declare that they have no conflict of interest.