A. Boundary conditions

The input excitation can be assessed by considering a three-dimensional terrain condition, including the contact model with the wheels; however, this analysis is computationally expensive and unproductive concerning the optimization process. Instead, the test vehicle is used for achieving a vertical dynamic condition. This vehicle is supported by a platform of vertical actuators that simulate a vehicle's contact with the ground. Figure 2 displays a representation of the simplification that has been used.

Figure 2 Vehicle placed on a virtual mechanical system of dynamic tests. 

The mechanism presented in Figure 3 is modeled to emulate the interaction between the wheel and the surface because the COMBIN14 elements do not control the freedom degrees between the interacting components.

Figure 3 Simplification of the tire contact problem. 

Similarly, for the front and rear leaf spring suspensions, models with sufficient degrees of freedom are included to allow their movement, as depicted in Figure 4 and Figure 5. The bars control the vertical movement of the rear axles through spherical joints. These joints allow the angular displacement of the axles against the cardan and control spring contact in any terrain condition.

Figure 4 Multibody system simplification of the front leaf springs. 

Figure 5 Multibody system simplification of the rear leaf springs. 

B. Force conditions

The total time that has been defined for performing the simulation is 7 s. For this time-lapse, static balance and test track route are the two states that are considered.

For the static balance, an interval of 2 s is allocated where the vehicle mechanically adapts to the standard gravity acceleration. Furthermore, a 5 s interval is assigned to the second state wherein the excitation functions begin to vibrate on the test vehicle according to the gap, depending on the simulation's speed. In this case, an average speed of 45 km/h is defined following previous research ¹². The gap includes the same excitation function with a time delay for the rear axles (axes 2 and 3).

The excitation function for the hydraulic actuators is randomly adopted by considering adverse conditions ¹. Figure 6 depicts the excitation function for the right and left terrain profiles.

The conditions for performing the analysis are controlled using Δt, which starts with ITS = 1e -2 for all steps and MTS = 1 -e5 in case of divergences that may occur when the input parameters change in an experimental design. Additionally, the APDL solver uses the HHT method to resolve the transient equation system.

Figure 6 Off road terrain profile in the random function. 

C. Design of experiments (DOE)

The design of the experiments used in this study adopts full factorial by considering that the computational cost of a single run is relatively low compared with that of a dynamic MBS + FEM simulation for the same vehicle, including the discrete model of leaf springs. The latter could increase the computational cost by up to 28,800 s (8 h) per run.

The experiment involves executing a series of transient ride/handling of dynamic simulations by continuously adjusting the leaf springs, dampers, and payload features according to the vehicle specification ¹⁰. For performing the lateral dynamics simulation, only the stiffness and damping parameters of the front leaf springs have been defined. Further, the payload remains constant throughout the experiment.

In addition, the entire factorial arrangement provides a high resolution of the results that may be expressed practically on the response surface.

D. Experimental factor levels

The virtual experimental arrangement is defined in equation (5) for a total of 160 treatments, where the factor is the average constant of the front spring assessed when 160 N/mm < k < 528 N/mm. The factor is the average damping factor of the front damper that is assessed when 1 Ns/mm < C < 15 Ns/mm. Further, the factor is the vehicle payload assessed when 1475 kg < M < 17000 kg. An experimental arrangement based on multiple levels of the experimental factors is preferred to estimate better the non-linear effects expected in a suspension system.

(5)

The k and C factors are used similarly for performing the "handling" simulation. Further, the payload value remains constant at 17000 kg for performing this last simulation.

The design of the experiment is established considering the most influential independent factors in the dynamic response of a vehicle suspension. Furthermore, with the different ranges covered by the independent factors, an attempt is made to cover a wide set of operational possibilities for trucks with the selected characteristics ¹².

E. Objective function

A metamodel is used to define the type of response surface. This standard regression model provides a statistical analysis of the relation between two or more quantitative variables, allowing us to estimate the dependent variable from others ^12,14.

The regression analysis assumes that there are n samples with a corresponding known output value for each sample. Further, the regression analysis determines the relation between the input and output parameters based on the sample points. Typically, a second-order polynomial is used for the regression model. This model approximates a real "input-output" relation. An exact relation is reached only in special cases ^(12,16,17.

The optimization function for the "ride" simulation depends on the independent variables with restrictions, and it is expressed in equations (6-8).

(6)

(7)

(8)

where is the average of peak accelerations during the simulation, k ₁ is the stiffness constant for the right front spring, k ₂ is the stiffness constant for the left front spring, C₁ is the damping constant for the right side, C₂ is the damping constant for the left side, M _PL is the load transported by the vehicle, g is the acceleration due to gravity, A _rms is the RMS acceleration, and f _n is the natural frequency of the vehicle.

The optimization function for the "handling" simulation depends on the independent variables without any restrictions, and it is expressed in equations (9-11).

(9)

(10)

(11)

Where θ is the angle of lateral rotation of the suspended mass, is the rotation speed of the said mass, and A _y is the maximum lateral acceleration measured at the head of the test dummy.

The objective functions are defined following the ranges established in studies performed to investigate the passenger and cargo transport vehicles by defining the amplitudes and oscillation frequencies of the vehicle in accordance with ISO 2631-1.

The MOGA optimization algorithm is a variant of the NSGA-II algorithm based on elitist control concepts, supporting multiple objectives and restrictions for finding the optimum global value. The number of samples for the generation of candidates is 500, with the maximum number of candidates being 6 because many responses have to be evaluated. In comparison with the direct optimization methods, the optimization process is quick, mainly because the response surfaces have already been generated as per the design of the experiment.

F. Detailed parabolic leaf spring design

Currently, there is a tendency to use a parabolic spring instead of a traditional laminated spring of constant width to improve the comfort and efficiency of transportation ^18-20. For the parabolic spring, the flat distribution of the stress allows the mass to decrease. Further, if the parabolic spring constant is the same as the constant of the traditional spring, the former will have an advantage in decreasing the dynamic constant, thereby improving the vehicle's comfort. The main stress must remain the same throughout the whole beam ^21-23 to reduce the mass of the spring. This reduction is achieved through a leaf spring thickness function using equation (12).

(12)

where t ₀ is the thickness of the sheet at a distance x, P is the applied force, b is the width of the sheet, and o _x is the effort at distance x.

However, when the stress remains constant throughout the beam ²⁴, the spring constant can be obtained using equation (13).

(13)

where k is the spring stiffness constant, E is the material's elastic modulus, Ix is the moment of inertia of the cross-section, b is the width of the sheet, t is the thickness of the sheet, and l is the length of the sheet.

According to equation (12) and equation (13), a preliminary calculation is performed for sizing the CAD model using 100-mm intervals. The mechanical properties of SAE 5160H steel are included by considering its fatigue behavior through the S-N curve. The S-N curve is constructed according to equation (14) ^(25,26, where the parameters consider the material and the manufacturing process. For example, the hardness obtained through the heat treatment and residual compressive stresses are derived from the blasting process.

(14)

where S _n is the alternating stress, N is the number of cycles, a is the value [MPa] composed by the average stress, and b is an exponent composed by the reduction factor of the fatigue strength.

The load history is extracted from the transient nonlinear simulations of the test vehicle in which the compression forces of the front springs are post-processed, as depicted in Figure 7 and Figure 8.

Figure 7 Dynamic load history for the "ride" test of the front springs. 

Figure 8 Dynamic load history for the "handling" test of the front springs. 

For reliability of 3σ or 99.73%, the maximum loads for both the events are 50,000 N, whereas the minimum loads are 10,000 N. These loads along with the average load (25,851 N) and empty vehicle load (23,975 N) are incrementally introduced in ANSYS for performing analysis together with an unloading step.

The pre-stress loads are also introduced onto the U-Bolt staples according to their SAE grades and diameters (74464 N) for assembling without misalignments during the application of test loads. The rainflow-counting method establishes the frequency of load values and their amplitude using the data from Figure 7 and Figure 8 as input information.

Further, the durability equivalence of the spring is defined according to the distance traveled by the vehicle during the ride/handling events (7 s) at a speed of 45 km/h, indicating that a full load cycle will be repeated for every 87.5 m traveled.

The stress component that has to be assessed is the main maximum stress, which occurs on the stress side of the material because it is the leading cause of fatigue failures in a leaf spring. The material's S-N curve is used as the failure criterion, and the fatigue of any other component within the assembly is not considered because they are not studied. In the actual work, it is decided not to use an objective function for fatigue resistance based on a safety factor, but rather it is preferred to use a defined objective function to maximize the fatigue life of the leaf spring.

The target function defined for the parameter optimization of the spring is expressed through equation (15).

(15)

subject to: (15)

A. Transient simulation results for the test vehicle

Figure 9 depicts that the acceleration peaks of the test dummy are higher than the peaks recorded in the engine support. These peaks may be attributed to the absence of damping in the cabin suspension system and the travel restriction presented by the suspension on the vertical axis, which produces a rigid stop.

Figure 9 Vertical accelerations measured in the motor support and in the head of the test dummy. 

The sinusoidal signal observed in a time interval of 0-2 s corresponds to the static equilibrium obtained by the vehicle because of acceleration due to gravity, which exerts influence from the first sub-step.

Figure 10 depicts the lateral acceleration measured in the dummy's head placed in the driver's seat, which undergoes a lateral movement at a moment when the vehicle changes lanes at a constant speed.

Figure 10 Lateral acceleration in the test dummy's head. 

B. Resulting forces

Figure 11 shows the operation ranges of the fluctuating loads exerted by each suspension spring. There is a considerable rear suspension load due to the location of the gravity center of the payload.

Figure 11 Vertical forces on the leaf springs for the left side and the three axles. 

C. Rotation results "Ride" rotations

Rotations are measured about the mass center of the vehicle because they represent the global movement of the structure. Figure 12 shows the rotation in the longitudinal direction against the transverse axis.

Figure 12 Pitch rotation about the suspended mass center of the vehicle. 

Figure 13 depicts the rotation of the vehicle on the longitudinal axis. This rotation is a consequence of the difference between amplitude and disparity owing to terrain irregularities.

Figure 13 Roll rotation about the suspended mass center of the vehicle. 

Figure 14 depicts the rotation angle for the vertical axis. In this case, the magnitude is considerably smaller than the angles described in Figure 16 and Figure 17 because the vehicle does not perform maneuvers using the steering wheel.

Figure 14 Yaw rotation about the suspended mass center of the vehicle. 

Figure 15 depicts a plot of the rotation speed vs. time on the lateral axis. This chart shows a significant magnitude when compared to those of other rotations owing to terrain irregularities.

Figure 15 Roll rotation speed for the "ride" simulation. 

"Handling" rotations

Figure 16 shows the roll rotation, which is the primary "handling" simulation objective.

Figure 16 Roll rotation for the "handling" simulation. 

Figure 17 depicts the roll rotation speed when the vehicle changes lanes.

Figure 17 Roll rotation speed for the "handling" simulation. 

D. DOE results

The results of the design parameters of the suspension elements and the loading conditions are presented. The adjustment is used to measure the quality of the response surface, and it is linked to the regression model that has been used, which is a second-order standard in this case. The correlation coefficient R² is the percentage of variation of an output parameter that can be explained by the regression equation of the response surface ^27,28. Figure 18 shows the normalized adjustment quality for the mean vertical accelerations and RMS. All points are included for the plot, considering both experimental factors under study. Although good results can be obtained for the mean accelerations, the same is not applicable for RMS due to the results' high dispersion. The high dispersion of the results may be due to the non-proportionality of some acceleration peaks derived from the unexpected behavior of the vehicle's movement obtained through the adjustments of different parameters.

Figure 18 Normalized quality adjustment for the mean vertical acceleration and the square root of the mean results. 

Moreover, the response surfaces are presented, where the output variable is on the vertical Z-axis. Figure 19 depicts the influence of input parameters, elastic constant, and damping coefficient on the deflection of the left front leaf springs. The damping variation is considered to be important for the low stiffness constants. For the front spring on the right side, similar conditions are presented. However, Figure 20 depicts the response surface for the vertical acceleration that the test dummy has perceived. Both response surfaces are obtained for a constant payload of 17,000 kg. As can be seen in these last two figures, both experimental factors studied have a significant influence on the vehicle's dynamic response.

Figure 19 Response surface for the deflection in front left leaf springs. 

Figure 20 Response surface for the mean vertical acceleration on the test dummy in a loaded condition. 

Figure 21 depicts the lateral acceleration of the test dummy that has been placed in the driver's seat. The different lines denote the influence of damping, which does not seem considerably strong, especially for low stiffness.

E. Leaf spring modeling results

Also, the FEM and optimization results (structural dynamic fatigue) are presented by considering the results of the vehicle optimization model, which is selected according to the response surfaces and the influence factor as the input parameters. The spring constant parameter is a restriction parameter that the optimization process must comply. It can be concluded that the stiffness factor of the spring extensively influences the vertical and lateral behaviors of the vehicle.

Figure 21 Response surface for lateral acceleration measured in the head of the test dummy. 

The FEM results are part of the preliminary nonlinear model. This study defines the smallest number of possible leaf sheets for the spring based on the functional and manufacturing conditions. Functional conditions adopt a better spring behavior with less friction between sheets, less weight, and less space within the suspension system; however, there is a hazard that a single sheet bearing full responsibility will fail to owe to overload or fatigue. The manufacturing conditions refer to low manufacturing times and costs because of the low number of sheets; therefore, a smaller number of assembly accessories are required.

F. Displacement results

Figure 22 depicts the displacement results for a two-sheet parabolic spring at maximum load and its kinematic interaction with the suspension components.

Figure 22 FEM displacements for parabolic springs under maximum load. 

G. Main stress results

Figure 23 depicts the linear stresses for both sheets at different load levels. The dashed lines represent the control points for the experiment design in which the mean stress and its standard deviation are controlled to maintain the stress balance. In total, seven control points are generated for each spring sheet.

Figure 23 Linear stress for each sheet under different load conditions. 

H. Candidate selection

Table 1 presents three possible candidates for the optimal design according to the objectives that have been set using the final geometric adjustments.

Table 1 Optimal spring design candidates by thickness (dimensions in mm). 

The durability, natural frequency, k-stiffness constant, mass results, and their variation against candidate 1 are reported for each candidate. The latter turns out to be the optimal design alternative based on the proposed objectives.

Therefore, candidate 1 is selected to close the optimum leaf spring design cycle. Figure 24 depicts the parabolic spring model with the final dimensions obtained through the optimization process for candidate 1 and specifies the dimensions subjected to modification.

Figure 24 Final dimensions for candidate 1 (in mm).