Prediction of frequency response of sub-frame bushing and study of high-order fractional derivative viscoelastic model

This paper presents experimental and dynamic modeling research on the rubber bushings of the rear sub-frame. The Particle Swarm Optimization algorithm was utilized to optimize a Backpropagation (BP) neural network, which was separately trained and tested across two frequency ranges: 1–40 Hz and 41–50 Hz, using wideband frequency sweep dynamic stiffness test data. The testing errors at amplitudes of 0.2 mm, 0.3 mm, and 0.5 mm were found to be 1.03%, 3.05%, and 1.96%, respectively. Subsequently, the trained neural network was employed to predict data within the frequency range of 51–70 Hz. To incorporate the predicted data into simulation software, a dynamic model of the rubber bushing was established, encompassing elastic, friction, and viscoelastic elements. Additionally, a novel model, integrating high-order fractional derivatives, was proposed based on the frequency-dependent model for the viscoelastic element. An enhanced Particle Swarm Optimization algorithm was introduced to identify the model's parameters using the predicted data. In comparison to the frequency-dependent model, the new model exhibited lower fitting errors at various amplitudes, with reductions of 3.84%, 3.61%, and 5.49%, respectively. This research establishes a solid foundation for subsequent vehicle dynamic modeling and simulation.

In simulation software such as MSC Adams, rubber bushing is expressed in the form of a dynamic model, so a more accurate rubber bushing dynamic model can help improve the accuracy of the simulation model in the simulation software.Further improve the reliability of pre-product development and product optimization.
This paper focuses on the rubber bushing of the rear subframe of a vehicle.The rubber bushing is shown in Fig. 1, where the X-direction represents the radial solid direction of the rubber bushing, the Y-direction represents the radial hollow direction of the rubber bushing, and the Z-direction represents the axial direction of the rubber bushing.The stiffness of different directions is different, so this paper focuses on Y-direction.
The rubber bushing was subjected to experimental analysis, with a focus on the Y-direction.The PSO-BP (Particle Swarm Optimization-Backpropagation) neural network was trained and tested across two frequency ranges: 1-40 Hz and 41-50 Hz.The trained PSO-BP neural network was then used to predict data in the frequency range of 51-70 Hz.In addition, to improve the accuracy of the rubber bushing's dynamic model, a high-order fractional derivative new model was proposed based on the Frequency-dependent model.The new model aimed to enhance the overall model accuracy.Then parameter identification was performed on the dynamic model, and a modified particle swarm optimization algorithm was proposed for parameter identification.

Rubber bushing experiment
Experiments are the most effective and intuitive method for studying the mechanical properties of rubber bushings.In this study, the rear-point rubber bushing of a vehicle rear suspension sub-frame was selected as the experimental object.Dynamic and static loading tests were conducted on the rubber bushing to obtain experimental data.The LETRY dynamic stiffness testing platform, as shown in Fig. 2, was used for the experiments.

Rubber bushing static loading test
The elastic and friction units of the bushing model established in this study are used to simulate the static behavior of the bushing.
The loading range in the X/Y direction is ± 12000N, and in the Z direction is ± 6000N.The experimental results are shown in Fig. 3. Due to the anisotropic nature of the rubber bushing studied in this paper, the static characteristics in the X/Y/Z directions exhibit significant differences.

Rubber bushing dynamic loading test
The dynamic loading test of the rubber bushing involves conducting wideband frequency sweep tests with sinusoidal excitations of different amplitudes.To fully investigate the dynamic characteristics of the rubber bushing, dynamic loading tests were performed in the frequency range of 1-50 Hz with amplitudes of 0.2 mm, 0.3 mm, and 0.5 mm.The relationship curve between the dynamic stiffness of the rubber bushing and the sweep frequency was obtained.
The experimental results for the X, Y, and Z directions are shown in Figs. 4, 5, 6.In the case of constant amplitude, the dynamic stiffness of the rubber bushing in the X/Y/Z directions increases with increasing frequency.Conversely, under constant frequency, the dynamic stiffness decreases with increasing amplitude.There is significant variation in the dynamic stiffness of the rubber bushing in different directions and amplitudes, particularly noticeable between the X and Y directions.In this study, the focus was on modeling the rubber bushing in the Y-direction, which exhibits higher dynamic stiffness.

Prediction of Dynamic Stiffness in the Y-Direction of Rubber Bushing using PSO-BP Neural Network
In the testing of rubber bushings, wideband frequency sweep tests are required at different amplitudes.This is especially crucial in the research process of NVH characteristics, where more frequency test samples are needed.However, the increased demand for test samples leads to higher costs and longer testing cycles.Additionally, during the testing process, the resonance of the testing machine itself can cause abrupt changes in the test data of the rubber bushings 14 .
To reduce the testing cost and cycle, BP neural networks are used to predict the test data of rubber bushings.In order to improve the prediction accuracy of the BP neural network, this study combines it with the PSO algorithm 15 .The BP neural network adjusts its network weights and thresholds based on the prediction error 16 .The PSO algorithm, proposed by Dr. Eberhart and Dr. Kennedy in 1995 17 , is an intelligent optimization algorithm inspired by birds searching for food.It searches for the optimal solution based on fitness by updating the position and velocity of particles.
The iteration formulas for updating velocity and position in the PSO algorithm is calculated as follows: (1)  www.nature.com/scientificreports/In Eq. 1, w represents the inertia weight; v t and x t denote the current particle's velocity and position; p b and g b respectively indicate the positions associated with the individual best fitness value and the global best fitness value; r 1 and r 2 represent random numbers within the range (0,1); and c 1 and c 2 are the learning factors; the veloc- ity and position of the particle have ranges of [v min , v max ] and [x min , x max ] , respectively.
The fitness function for the PSO-BP neural network can be expressed as follows: In the formula, n represents the number of output nodes in the neural network, y i denotes the desired output of the i-th node in the BP neural network, and o i represents the predicted output of the i-th node.The coefficient k is set to 1 in this study.
The parameter settings are provided in Table 1, and the PSO-BP neural network process is illustrated in Fig. 7.
The prediction results based on PSO-BP neural network test data can be shown by Fig. 8.According to Fig. 9, the comparison between the prediction results and experimental data results of the PSO-BP neural network can be observed. (2) Table 1.Parameter Settings for PSO-BP Neural Network.The particle dimension refers to the sum of the threshold and weight count of the entire neural network.The weight count is calculated as follows: 1 × 5 + 5 × 1 = 10, and the threshold count is 5 + 1 = 6.Therefore, the particle dimension is 10 + 6 = 16.

Establishment of the parameterized model for rubber bushing.
The dynamic model of the bushing, as shown in Fig. 11, can be constructed by parallel connections of elastic, frictional, and viscoelastic elements.In the figure, F e represents the elastic force in units of N ; F f represents the force of the frictional hysteresis element in units of N;F v and represents the viscoelastic force in units of N ; F represents the response force of the entire parameterized model in units of N.
Since the elastic, frictional, and viscoelastic elements are connected in parallel, the combined force of the three elements represents the response force of the entire bushing, as expressed in Eq. 3: (3)

Elastic element
The static characteristics of the bushing are caused by its elastic deformation.Constitutive models commonly used to describe the static mechanical behavior include the Mooney-Rivlin model 18 , Neo-Hookean model 19 , Yeoh model 20 , Ogden model 21 , etc.However, considering the elastic deformation characteristics of the bushing, in order to more flexibly accommodate the nonlinearity of the elastic element, a polynomial spring model is used to represent the static characteristics of the bushing.The polynomial spring can adjust the highest degree or coefficients to adapt to the nonlinearity of the elastic element.Its mechanical expression is as follows 22 : F e represents the force of the elastic element, measured in units of N .Under the influence of a sinusoidal excitation with an amplitude of x 0 , the amplitude of the elastic module is given by: The elastic element does not consider friction, so there is no energy loss.( 4)  www.nature.com/scientificreports/

Frictional element
The hysteresis effect of a rubber bushing becomes more pronounced with increasing deformation, and the nonlinearity becomes more evident.The expression for the smooth friction force model is as follows 23 : Among them, F f represents the frictional force, x represents the displacement of the loading, measured in units of mm; F fmax is the maximum frictional force, measured in units of N; x 2 is the displacement at which the frictional force increases from 0 to F fmax /2 , measured in units of mm;(x s , F fs ) represents a reference point on the force-displacement curve obtained from static loading tests.Under the influence of a sinusoidal excitation with an amplitude of x 0 , the amplitude of the frictional hysteresis module is given by: In the formula,u = F f 0 /F fmax , and E f represents the energy dissipation per cycle, measured in units of N•mm.

Viscoelastic element
For the dynamic model of the viscoelastic element in a rubber bushing, following the approach proposed by Liu Guo jia et al., a high-order fractional derivative model is derived based on the Frequency-dependent model.
The Frequency-dependent model is then modified to develop a new model called the High-Order Fractional Derivative Frequency-dependent model.The structure of the Frequency-dependent model and the new model is shown in Table 3.

Frequency-dependent model
The mechanical expression of the Frequency-dependent model is as follows: In the equations:x represents the loading displacement of the rubber bushing, measured in units of mm; z represents the displacement of the elastic element k 2 and the damping element c 2 , measured in units of mm; k 1 , k 2 , c 1 , c 2 are the elastic coefficients and damping coefficients of the model.( 6) www.nature.com/scientificreports/Substituting Eqs. 9 and 10 into Eq.11: In the formula, ż represents: Equations 14 and 15 are obtained from 12 and 13 through Laplace transformation: From Eq. 15, it can be concluded that: Equation 16 is substituted into Eq.14: From Eq. 17, the formula for calculating the complex stiffness can be derived: The complex stiffness converted to the frequency domain yields the following equation: By further deriving from the above equation, the amplitude of the real part and imaginary part of the response force, denoted as F v0Re and F v0lm respectively, under a sinusoidal excitation of amplitude x 0 , can be obtained:

New model
Because the Frequency-dependent model cannot accurately describe the viscoelastic properties of rubber, this paper proposes a Frequency-dependent model based on the High-Order Derivative Frequency-dependent model.The relationship between force and displacement in this new model is given by: In the equation, α and β represent the order of the fractional derivatives, which range (0,1) ; k 1 , k 2 , c 1 , c 2 are the elastic modulus and viscosity coefficients of the model, respectively; and 0 D α t is the Riemann-Liouville fractional derivative operator.

D α t can be definition by Eq. 23:
In other words, f (t) is first to do (n − α) fractional integration, and then take the n derivative, n is 1.
Based on the Frequency-dependent model, changing the damping to a sticky pot with fractional derivative of displacement can better describe the viscoelastic properties of rubber.When α and β are both 1, the new model will be equal to the Frequency-dependent model, so the mechanical properties of the new model already include the mechanical properties that the Frequency-dependent model can represent.From Eq. 23, the formula for calculating the complex stiffness can be derived: The complex stiffness of the new model in the frequency domain can be derived as: Setting n = 0 as the principal root, the result is as follows: From Euler's formula, we can obtain: Substituting Eq. 26 into Eq.23, we have: By further deriving from the above equations, the magnitudes of the real part and imaginary part of the response force, denoted as F v0Re and F v0lm respectively, under a sinusoidal excitation of amplitude x 0 can be obtained: )) (30) Vol.:(0123456789) www.nature.com/scientificreports/ The above equations will be used for the subsequent parameter identification of the viscoelastic element.

Parameter identification of the parameterized model for rubber bushing
In the parameter identification process of the rubber bushing model, the elastic unit and friction unit are first identified using quasi-static loading data.Subsequently, the viscoelastic unit is identified by combining the dynamic stiffness data.

Parameter identification of the elastic and friction units
Parameter identification is performed using static loading test data.The static elastic stiffness K e of the elastic unit, as shown in the Fig. 12, can be approximate by the slope of the curve near the limit position of displacement.The maximum friction force F fmax in the friction model expressed by half the vertical distance between the upper and lower limits of the hysteresis loop.The maximum slope of the curve is K max .The parameter x 2 in the friction unit can be determined using Eq. 32.
By aligning the upper and lower boundary curves of the hysteresis loop in Fig. 12 through translation, the force-displacement test curve for the elastic component is obtained.Using the data from this curve, a 3rd-degree polynomial spring model is fitted as shown in Fig. 13.
Results of parameter identification for the elastic unit and friction unit are shown in Table 4.

Parameter identification of the viscoelastic unit
The parameter identification of the viscoelastic unit involves a large number of parameters, resulting in a significant computational burden.Parameter identification for the viscoelastic unit is typically performed using algorithms such as least squares 24 , genetic algorithms 25 , and particle swarm optimization.
Because the model established in this paper has many parameters and strong nonlinear.The PSO algorithm is often used to solve the optimization problem with many parameters, wide range and strong nonlinear.Therefore, this paper selects PSO algorithm for optimization.However, the PSO algorithm is prone to premature convergence, meaning it may get trapped in local optima and fail to explore the entire search space.Due to the premature convergence problem of the particle swarm optimization algorithm, where it gets trapped in a local (32) www.nature.com/scientificreports/optimal solution, genetic algorithm has the ability of mutation.The proposed improved particle algorithm is proposed based on the genetic and mutation ideas.When the particle swarm optimization algorithm gets trapped in a local optimal solution, a new particle swarm is generated by mutating it to seek a better solution and thus avoid premature convergence.
To improve the speed of optimization, a random particle is selected from the particle swarm during the velocity update process.By controlling the particle's velocity update in three directions, the speed of particle optimization is enhanced, and it helps prevent getting trapped in local optima.The velocity update equation is as follows: In the Eq.33: c 3 is the learning factor; r 3 is a random number in the range (0,1); p s is the randomly selected particle from the current particle swarm.
To prevent getting trapped in local optima, the results of each optimization iteration are evaluated.If the historical best fitness of the particle swarm remains unchanged after the current iteration is completed, the entire particle swarm undergoes crossover and mutation operations similar to those in genetic algorithms.This generates new particles and changes the search direction, thereby avoiding local optima.
The crossover operation is performed in a real-valued encoding format.For particles that meet the crossover condition, one random particle is selected for the crossover operation.The specific method is as follows: In the Eq.34: x k represents the particle that satisfies the crossover condition, x l is the randomly selected particle, and σ is a random number in the range (0,1).
The mutation operation applies different mutation probabilities to different particles.Therefore, the particle swarm is sorted in ascending order based on their fitness values, where particles with higher fitness values have higher mutation probabilities.The specific method is as follows:  In the Eq.35: P m represents the mutation probability, and i represents the index of the particle in the popula- tion, ranging from 1 to n.
The mutation operation selects the j-th dimension of the i-th particle for mutation.The mutation method is as follows: In the Eq.36: max(j) represents the upper bound of the j-th dimension of the particle, min(j) represents the lower bound of the j-th dimension of the particle, and r is a random number in the range (0, 1).
Fitness function of the algorithm: In the Eq.37: n represents the number of operating conditions being considered; k i dyn_t represents the experimentally measured dynamic stiffness data; k i dyn represents the dynamically calculated dynamic stiffness for the i-th operating condition.
During the identification process, it is necessary to ensure that the data fitted during the model calculation does not have significant errors.Therefore, a constraint is established: Calculate the dynamically calculated dynamic stiffness of the bushing using Eq.39: The two types of model parameters will be identified using the MPSO (Modified Particle Swarm Optimization) and PSO (Particle Swarm Optimization) algorithms separately.The MPSO algorithm will follow the process outlined in Fig. 14.
The specific parameter settings are provided in Table 5 and Table 6.Both models are selected for parameter identification using dynamic stiffness data at frequencies of 1, 10, 20, 30, 40, 50, 60, and 70 Hz, with an amplitude of 0.2 mm.
To verify the reliability of the MSPO algorithm proposed in this paper, the optimization effects of the PSO algorithm, the Adaptive chaotic particle swarm optimization (ACPSO) algorithm, and the MPSO were compared.From Fig. 15 to Fig. 16, it can be observed that both for the new model and the frequency-dependent model, the MPSO algorithm demonstrates stronger optimization capabilities compared to the PSO algorithm and the ACPSO optimization algorithm.
For the new model, the particle is represented as x = (k 1 , k 2 , c 1 , c 2 , α, β) .Both MPSO and PSO algorithms will have a maximum of 300 iterations.
For the frequency-dependent model, the particle is represented as x = (k 1 , k 2 , c 1 , c 2 ) .Both MPSO and PSO algorithms will have a maximum of 300 iterations.
From Fig. 17, it can be observed that at a 0.2 mm amplitude, the new model exhibits better fitting performance compared to the frequency-dependent model.The error results are shown in Table 7.
The error results also indicate that the new model exhibits better fitting performance, thereby improving the model accuracy.
The identification results of unknown parameters of the two dynamic models are shown in Table 8.Both models were selected for parameter identification using dynamic stiffness data at frequencies of 1, 10, 20, 30, 40, 50, 60, and 70 Hz, at amplitudes of 0.3 mm and 0.5 mm.The results are shown in Fig. 18.
From Fig. 18, it can be observed that for both 0.3 mm and 0.5 mm amplitudes, the new model exhibits better fitting performance compared to the frequency-dependent model.The error results are shown in Table 9.
The error results also indicate that the new model exhibits better fitting performance, thereby improving the model accuracy.

Conclusion
This study conducted experiments on the rear suspension sub-frame rubber bushing of a certain electric vehicle model.The Y-direction of the rubber bushing was selected as the research object, and the PSO-BP neural network was used to predict the dynamic stiffness test data of the rubber bushing.To improve the accuracy of the bushing model, a new model was proposed based on the Frequency-dependent model when establishing the dynamic model.The parameter identification of the models and a comparison of the fitting accuracy between the two models were performed, leading to the following conclusions: (1) The proposed MPSO algorithm for parameter identification demonstrated stronger optimization capability compared to the PSO algorithm, highlighting its practical value.

Figure 1 .
Figure 1.(a) is an actual vehicle; (b) is an actual vehicle chassis; (c) is the rubber bushing form the actual vehicle.

Figure 2 .
Figure 2. (a) is the LETRY dynamic stiffness testing platform; (b) is the rubber bushing X loading.

Figure 3 .
Figure 3. Rubber bushing static loading test data in X/Y/Z-directions.

Figure 4 .
Figure 4. Relationship between dynamic stiffness and frequency in X-direction for different amplitudes of rubber bushing.

Figure 5 .
Figure 5. Relationship between dynamic stiffness and frequency in Y-direction for different amplitudes of rubber bushing.

Figure 6 .
Figure 6.Relationship between dynamic stiffness and frequency in Z-direction for different amplitudes of rubber bushing.

Figure 9 .
Figure 9.Comparison of prediction results from PSO-BP neural network with experimental data.

Figure 10 .
Figure 10.Difference between predicted values and experimental values for 1-50 Hz.

Figure 13 .
Figure 13. the curve of the elastic unit.

Figure 16 .( 2 )
Figure 16.Comparison of optimization algorithms for the frequency-dependent model.

Table 2 .
Prediction errors of the PSO-BP neural network.

Table 4 .
Parameter identification results for the elastic and friction units.

Table 5 .
MPSO and PSO parameter settings for the new model.

Table 8 .
Identification Results of Rubber Bushing Parameters at 0.2 mm Amplitude.Fitting results of models mm and 0.5 mm data.