Optimal design of leaf springs for vehicle suspensions under cyclic conditions

The suspension systems are designed to provide good performance in terms of comfort and maneuverability and to satisfy other requirements such as fatigue strength. This study focuses on the leaf springs, a classic mechanism; leaf springs are still being extensively used in several types of vehicles because of their high load capacity and low manufacturing and maintenance costs. Its dynamic behavior assesses stationary nonlinear preload state components to provide considerable added value to this suspension type. This assessment considers the contact condition of the suspension’s components and the large deflections and tightening torques observed in the whole assembly. Furthermore, the components of the non-suspended mass subsystems, such as tires, shock absorbers, and stabilizer bars, are characterized according to the simplified models for reducing their computational cost. In addition, a commercial test vehicle is used for simulating the complete system using three-dimensional modeling for describing its most relevant components in terms of their mass and rigid connection. The vehicle is additionally analyzed using multibody system simulations (MBS) coupled with the finite element method (FEM) in an implicit nonlinear transient environment using the ANSYS APDL solver. This dynamic simulation is parameter-driven for obtaining the experimental design and determining the optimal suspension stiffness and damping features required for transporting suitable load sizes.


RESUMEN
Los sistemas de suspensión para un vehículo se diseñan para brindar un mejor rendimiento en cuanto a confort y maniobrabilidad, además de otros requerimientos como la resistencia a la fatiga. Este trabajo centra su interés en los resortes tipo ballesta, los cuales son ampliamente utilizados en una gran variedad de vehículos, deben su vigencia a su gran capacidad de carga y los bajos costos de manufactura y mantenimiento. Para dar un mayor valor agregado a este tipo suspensión se ha analizado su comportamiento dinámico desde el punto de vista de componente desde un estado estacionario de precarga no-lineal. Este incluye la condición de contacto de sus propios componentes, las grandes deflexiones y los torques de apriete de ensamble de todo el conjunto. Los componentes del subsistema de masa no suspendida tales como neumáticos, amortiguadores y barra estabilizadora se caracterizan conforme a modelos simplificados para disminuir el costo computacional. Para llevar a cabo la simulación del sistema completo se define

INTRODUCCIÓN
The operation of freight vehicles on highways imposes several demanding requirements due to the topography, orography, climate, and road infrastructure conditions, which vary significantly in different countries [1,2]. Several studies have been carried out to improve the efficiency and functionality of the vehicles, the researcher Rocha-Hoyos and associates [3] presents an alternative to reduce pollutant emissions from go-kart 7.5 kW gasoline vehicles in the city of Baños in Ecuador, considered before making the EV conversion, the battery and electric motors, to achieve maximum efficiency.
Other research [4] indicates that the functional and morphological dimensions require a better approach, more precision, and mutual interaction in the design space; diminishing the uncertainty in the handling and characterization of conceptual ideas is an important task in the stage of conceptual design. Thus, models must be described in the most straightforward manner to optimize and facilitate the making of the design process both in engineering and product design.
Some researchers indicate suspension is a subsystem that considerably impacts the response of land vehicles to terrain conditions, which are directly related to various aspects such as dynamic comfort, maneuverability, and performance. These aspects are because the suspension is responsible for maintaining contact between tires and road, for providing insulation against the vibrations produced by irregularities in the terrain, for responding to angles of inclination, for steering to avoid road irregularities, and for resisting the chassis rotation in its longitudinal and lateral axes [5,6].
On the research [7] has been evaluated the dynamic behavior of suspension-type monorail (STM) system, a coupled dynamic model of the STM vehicle-bridge system has been developed based on multibody dynamics and finite element (FE) theory, finding that the simulation with the proposed dynamic model is in good agreement with the listed test data. Some discrepancies were resolved when the bridge is treated as a rigid body and flexible body, respectively, showing the importance of considering the flexible bridge and demonstrating several modeling advantages of the proposed coupled dynamic model.
Recent studies have focused on achieving maximum dynamic and structural performances for suspension of these types of springs based on the characterization of various factors that influence the dynamic stiffness and a set of simulation and optimization techniques. Further, the influence of friction, frequency, and amplitude is investigated to improve vehicle comfort [8,9].
This study intends to model and optimize the suspension system of cargo vehicles by considering the application of the vehicle, operating conditions to which it will be subjected, different geometry constraints, expected dynamics, structural behavior, and robustness.
This study comprises an initial stage where preliminary calculations have been performed for a given commercially available vehicle. A workflow in V cascade type must be proposed, where the importance of each step is specified; however, it should allow for a certain level of automation. Figure 1  transmission. The information required for component modeling and characterization was extracted from the technical literature [10], which was originally in the form of a two-dimensional AutoCAD file. In addition, this study considers the information that has been extracted from the technical data sheet [11], which defines some aspects, including the centers of gravity of the permissible loads, cabin class and its turning radius, structural chassis profiles, disposition of some transversal reinforcements, and placement of staples and anti-vibration springs.
Modeling using the finite element method and multibody systems (MBS) for the test vehicle is simulated considering the complete structure (imported from the CAD system) [12,13]. The latter includes the suspension, chassis, and cabin subsystems for obtaining relevant results for a time interval.
Further, transient dynamic analysis is used to determine the dynamic response of a structure under the action of loads against time [14,15]. The basic expression of movement for solving transient analysis is defined by equation (1).
Where [M] denotes the mass matrix, [C] denotes the damping matrix, [K] denotes the stiffness matrix, denotes the nodal acceleration vector, denotes the nodal velocity vector, denotes the nodal displacement vector, and {F(t)} denotes the force vector. After every t time, this equation is divided into a series of equations in a dynamic equilibrium considering both inertial and viscous damping forces.
In the case of nonlinear structural dynamics, the internal forces are not proportionally linear to the nodal displacement; further, the stiffness matrix is dependent on the current displacement [10]. Instead of equation (1), any integration scheme must be applied to the nonlinear semi-discrete equation ( Where is a vector of internal forces and is a vector of the applied forces. The Hilbert-Huang transform (HHT) method is selected because it is considered stable and precise and provides better numerical dissipation when compared to that provided by the Newmark method. Further, the conditionally stable integration algorithms are affected by the size of the time interval Δt. When an algorithm is unconditionally stable, the size of Δt may be selected regardless of the stability.
According to a previous study [11], the iterative Newton-Raphson process can solve the nonlinear equation, which is defined by equation (3).  [10].
Where {u i } denotes the displacement in iteration i and {u i+1 } denotes the displacement vector in iteration i + 1.

MATERIALS AND METHODS
There are several software applications for simulating vehicle behavior. Programs, such as ADAMS/CAR, include different vehicle libraries, terrains, and driving conditions that facilitate the dynamic simulation and automation of the design process for any suspension and chassis component. In this study, the ANSYS software performs a dynamic vehicle simulation, which does not contain any special modules for simulating specific events. This simulation includes the parameters of the previously characterized components and the properties of the dampers and suspensions of the driver's cabin. The simulation also includes a test dummy, which contains an accelerometer for measuring the increase in speed to which it will be subjected.
The generated mesh comprises 438 elements, among which 190 elements are BEAM189, and 248 are MASS21 elements. The stiffness matrix for this system will be symmetric and in leaf, which is also applicable to the mass and damping matrices.

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.
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.
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.

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 [12]. 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 [1]. 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.

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 [10]. 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.
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 [12].

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 "inputoutput" 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)(7)(8).
min A rms = A rms k 1 = k 2 ,C 1 = C 2 , M PL ( ) subject to : A rms ≤ 0,12g, min f n = f n k 1 = k 2 , where A z is the average of peak accelerations during the simulation, k 1 is the stiffness constant for the right front spring, k 2 is the stiffness constant for the left front spring, C 1 is the damping constant for the right side, C 2 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)(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][19][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][22][23] to reduce the mass of the spring. This reduction is achieved through a leaf spring thickness function using equation (12).
where t 0 is the thickness of the sheet at a distance x, P is the applied force, b is the width of the sheet, and σ x is the effort at distance x.
However, when the stress remains constant throughout the beam [24], the spring constant can be obtained using equation (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 postprocessed, as depicted in Figure 7 and Figure 8. 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).
subject to : k s = k opt. dyn. sist. (15) subject to: (15) 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.

A. Transient simulation results for the test vehicle
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 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.

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 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 11. Vertical forces on the leaf springs for the left side and the three axles. 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 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.
"Handling" rotations Figure 16 shows the roll rotation, which is the primary "handling" simulation objective. Figure 17 depicts the roll rotation speed when the vehicle changes lanes.

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 2 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.
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 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.   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. Figure 22 depicts the displacement results for a two-sheet parabolic spring at maximum load and its kinematic interaction with the suspension components. 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. Table 1 presents three possible candidates for the optimal design according to the objectives that have been set using the final geometric adjustments.

H. Candidate selection
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.

CONCLUSIONS
The applied methodology allows for variations in the CAD model of a vehicle that may be quickly   components because these can be easily replaced by the MASS21 elements and MPC184 rigid connections by adjusting their mass centers to the field measurements that have been previously conducted.
The proposed methodology allows for the direct integration of kinematic and commitment analysis in the vehicle model as well as in the geometric parameter definition and position of all elements, thereby enabling evaluations for the entire system.
The proposed methodology also allows for the integration of "ride" and "handling" simulations for performing a joint assessment under an experimental design and for the optimization of not only the springs but also any other suspension component.
Further, the virtual FEM characterization of the suspension components supports proper modeling of its static-dynamic characteristics in the vehicle,  adjusted in the CAE model to assess different thickness segments quickly.
The CAD models of vehicles do not necessarily need to contain the same level of detail in all its thereby allowing for a high computational saving in terms of future experiments.
The optimization process for dynamic events should be integrated with the leaf spring design to generate a response surface related to fatigue life and to include it in the Pareto boundary.
Furthermore, the design of full factorial experiments is useful in automated simulation processes because it avoids any manual user intervention to change the parameters. If an experimental design needs to be manually created, a fragmented or Taguchi design must be applied.
The vertical dynamics responses in the design of the experiment conflict with the lateral dynamics responses, thereby limiting the stiffness constant when lateral maneuvers are performed.
The durability objective is to find a design boundary as the stiffness constant decreases when the durability is considerably low for the proposed restriction.