Sage Journals: Discover world-class research

Abstract

The dynamic model of the front double wishbone suspension and the rear multi-link suspension of the vehicle are established. On the basis of detailed analysis of suspension kinematics, calculation method of wheel alignment angle and force calculation of suspension bushing, the influence mechanism of suspension bushing on the vehicle transient state is clarified, and the vehicle transient characteristic index is derived from the vehicle three-free dynamic model. The sensitivity analysis of the suspension bushing is carried out, and the bushing stiffness which has a great influence on the transient state of the vehicle is obtained. The bushing stiffness scale factor is used as the optimization variable, the vehicle transient characteristic index is used as the optimization target, and the NSGA-II optimization algorithm is used for multi-objective optimization. After optimization, one Pareto solution is selected to compare with the original vehicle, the comparison results show that the yaw rate gain, resonance frequency and delay time of yaw rate in the vehicle transient characteristic index are all improved, other optimization targets change less. In the steady-state comparison, the understeer tendency of the vehicle increases, and the roll angle of the vehicle increases but is within an acceptable range.

Keywords

Suspension kinematics suspension compliance bushing deformation vehicle transient characteristics multi-objective optimization

Introduction

The K (Kinematics) characteristic and C (Compliance) characteristic of the suspension are important performances of the suspension, which largely determines the handling stability of the vehicle. During the movement of the vehicle, the suspension and the suspension bushing play an important role. The deformation of the bushing will affect the wheel alignment angle of the vehicle’s suspension, and will also affect the force of the vehicle’s tires, and ultimately affect the vehicle’s dynamic characteristics. It is of great significance to carry out research on suspension and bushing, especially the optimization of the stiffness of the suspension bushing. At present, many experts and scholars have carried out a lot of research on suspension bushing, suspension movement and suspension force.

García et al. studied the frequency and amplitude dynamics of a commercial bushing with carbon black filled rubber bushings in the axial and radial directions. Based on the measured data, the axial and radial dynamic stiffness models of the rubber bushing are established. The accuracy of the new engineering model describing the axial and radial dynamic characteristics of carbon black filled rubber bushing is verified by experiments. This work allows the bushing characteristics in a complex system to be calculated quickly and easily.¹ Gibanica proposed a new method of calculating the stiffness of the bushing. In the experiment, four bushings were used for comparative analysis to quantify the uncertainty of the bushing parameters, and an uncertainty quantization method based on bootstrapping linear alternative model within parameters was proposed. For three identical rear subframes, finally a model update program using damping equalization has been successfully used to update all bushes of the rear subframe and find the optimal parameters of the bushes.² Plagge et al. proposed a model that is based on physical ideas and reasonable assumptions about the microstructure of materials and is designed for high efficiency and robustness in finite element applications.³ Sedlaczek et al., proposed an advanced phenomenological Daimler model of rubber bushings. The modeling complexity of this model can be easily adapted to the specific problem requirements of various simulation tasks. Examples of two typical chassis bushings under durable load conditions demonstrate the practical applicability and superior performance of the modular and adaptable Daimler model. The modular approach allows the basic building blocks of a specific combination of problems, representing stiffness, static friction and frequency-dependent dynamics. They can be easily parameterized using Excel’s graphical preprocessing, and the achievable model accuracy can be directly checked and adjusted to current requirements. In order to overcome the problem of low accuracy due to the limitation of expressing the frequency-dependent dynamic characteristics of the vibration isolator.⁴ Aydemir and Sendur⁵ developed and proposed a simplified transfer function model that varies with frequency, by generating complex stiffness and damping varying with frequency from the horizontal measurement of the vehicle. The proposed method has been verified on a heavy commercial truck, and the accuracy of the proposed model is better than that of Voigt modeling method. Zenglian et al. first derives the mathematical equations of the Berg model and Dzierzek model, and discusses the theoretical methods of parameter identification. Secondly, the static and dynamic characteristics of the rubber sleeve were tested, and the parameters of the two models were determined based on the experimental data, and the prediction accuracy of the two models was analyzed. On the basis of detailed analysis of dynamic stiffness and damping, a new empirical model is established using Berg model and Dzierzek model. The model can achieve a reasonable compromise between prediction accuracy, identification difficulty and calculation amount, and is an ideal tool for vehicle chassis dynamics simulation and analysis.⁶ Sancibrian et al. use a spatial synthesis method to design the kinematics of the vehicle’s double wishbone suspension system, using a multi-angle synthesis technology to establish the RSSR-SS spatial link, and express it through a constraint equation. Numerical examples show that this method can effectively improve the design of double wishbone suspension system and has good robustness. This study can achieve high accuracy in the definition of functional parameters during suspension movement, and provides a better reference for suspension design.⁷ Balike et al. proposed a comprehensive power quarter vehicle model to study the kinematics and dynamics of the link suspension and the influence of link geometry on selected performance measures. Combined with the kinematics of the double wishbone suspension, a plane two-degree-of-freedom model was established.⁸

The above experts combined different bushing structures and materials to put forward an equivalent model of bushing stiffness, and analyzed the characteristics of the bushing. From different perspective, use the different methods to study vehicle suspension. However, the above experts still have some problems that have not been solved. First of all, experts did not analysis the relationship between bushing deformation and force in detail, and did not calculate the amount of bushing deformation and force during suspension movement. Secondly, the experts did not analysis the impact of suspension kinematics characteristics on the transient and steady-state characteristics of the vehicle, and the impact of the suspension compliance characteristics caused by the deformation of the suspension bush on the transient and steady-state characteristics of the vehicle. The experts did not conduct sensitivity analysis on the stiffness of the suspension bushings, and did not find the bushes that have a greater impact on the transient characteristics of the vehicle. The third is that the experts did not match and optimize the bushing stiffness to improve the transient performance of the vehicle.

Aiming at the problems that experts have not solved yet, this article will start our research from three aspects. The first will analysis the impact of suspension movement on the wheel alignment angle and calculate the value of the wheel alignment angle, and analysis the displacement and force of the specific suspension, especially the displacement and force analysis of the multi-link suspension. In addition, the deformation and force of the bushing during suspension motion will be calculated, and the influence of bushing deformation on the compliance of the suspension will be analyzed in detail. The second aspect will analysis the impact of the bushing on the transient characteristics of the vehicle, focusing on explaining the influence mechanism of the bushing on the transient and steady-state characteristics of the vehicle. Suspension bushings will be analyzed for sensitivity, and suspension bushings that have a greater impact on the vehicle’s transient conditions will be selected. In the third aspect, the selected bushing stiffness is used as the optimization variable, and the frequency characteristic index of the vehicle under the sine swept steering input is used as the optimization target to match and optimize the suspension bushing, thereby improving the dynamic characteristics of the vehicle.

Establish vehicle suspension model and calculation of rubber bushing’s stiffness

In order to analysis the kinematics characteristics and compliance characteristics of the vehicle suspension and the influence of the suspension bushings on the suspension performance, the multi-body dynamic model of the vehicle suspension was established according to a real vehicle,⁹ as shown in Figure 1.

Figure 1.

Multi-body dynamic model of vehicle’s suspension: (a) front double wishbone suspension and (b) rear multi-link suspension.

The rubber bushing is a rubber metal piece composed of metal and rubber, the outer cylinder and the inner cylinder are composed of metal, and the rubber material is in the middle, the rubber bushing is integrated after vulcanization. In order to accurately analysis the dynamic characteristics of the vehicle suspension system, the mechanical properties of the rubber components need to be analyzed accurately. The structure diagram of the cylindrical rubber bushing is shown in Figure 2.

Figure 2.

Schematic diagram of rubber bushing’s structure.

The rubber bushing is fixed on the metal inner sleeve, and a radial ( $y$ or $z$ axis direction) load $W$ is applied to the metal outer sleeve to cause the outer sleeve producing a displacement $S$ . In theory, the dimensionless radial stiffness coefficient $n$ is often used to characterize the ability of the bush to resist radial deformation. The mechanical model of the rubber bushing is shown in Figure 3 and the stiffness curve of the bush is shown in Figure 4.

Figure 3.

Schematic diagram of bushing mechanical model.

Figure 4.

Diagram of rubber bush’s stiffness.

The radial stiffness coefficient of the bushing is defined as:

n = \frac{W}{GSL}

(1)

Among them, $G$ is the shear modulus of rubber. The radial stiffness is:

K_{ty} = \frac{W}{S} = nGL

(2)

The approximate calculation formula of radial stiffness coefficient is:

n = \frac{π}{\frac{1}{4} [\ln (\frac{r_{2}}{r_{1}}) - \frac{{r_{2}}^{2} - r_{1}^{2}}{{r_{2}}^{2} + r_{1}^{2}}] + \frac{1}{10 [c_{1} {(\frac{L}{r_{1}})}^{2} + c_{2}]}}

(3)

among them,

{\begin{matrix} c_{1} = - \frac{1}{60 s_{1}} \\ c_{2} = \frac{s_{2}}{s_{1}^{2}} \\ s_{1} = \frac{{r_{2}}^{2} \ln^{2} (\frac{r_{2}}{r_{1}})}{2 ({r_{2}}^{2} - {r_{1}}^{2})} + \frac{{r_{2}}^{2} \ln (\frac{r_{2}}{r_{1}})}{4 ({r_{2}}^{2} + {r_{1}}^{2})} - \frac{3 ({r_{2}}^{2} - {r_{1}}^{2})}{16 {r_{1}}^{2}} \\ s_{2} = \frac{{r_{2}}^{4} \ln^{3} (\frac{r_{2}}{r_{1}})}{4 {({r_{2}}^{2} - {r_{1}}^{2})}^{2}} + \frac{{r_{2}}^{4} \ln^{2} (\frac{r_{2}}{r_{1}})}{8 ({r_{2}}^{4} - {r_{1}}^{4})} - \frac{11 {r_{2}}^{2} \ln (\frac{r_{2}}{r_{1}})}{96 {r_{1}}^{2}} + \frac{({r_{2}}^{2} - {r_{1}}^{2}) (7 {r_{2}}^{2} - 2 {r_{1}}^{2} {r_{2}}^{2} + 7 {r_{1}}^{4})}{576 {r_{1}}^{4} ({r_{2}}^{2} + {r_{1}}^{2})} \end{matrix}

(4)

When the rubber bush is subjected to an axial ( $x$ axial direction) load, the elastic theory can be used to derive the axial stiffness calculation formula:

K_{tx} = \frac{- 2 π LG {(1 - \ln \frac{r_{2}}{r_{1}})}^{2}}{[\frac{3}{r_{1}} - 2 (\frac{1}{r_{1}} - \frac{1}{r_{2}}) \ln \frac{r_{2}}{r_{1}} - (1 + \ln \frac{r_{2}}{r_{1}} + \ln^{2} \frac{r_{2}}{r_{1}}) \frac{2}{r_{2}} - \frac{1}{r_{2}} (1 + \ln^{2} \frac{r_{2}}{r_{1}})]}

(5)

The calculation of torsional stiffness (around the x-axis) is:

K_{rx} = 4 π LG {(\frac{1}{{r_{1}}^{2}} - \frac{1}{{r_{2}}^{2}})}^{- 1}

(6)

Kinematics analysis of vehicle’s suspension

The influence from suspension’s movement on the vehicle’s handling stability is specifically reflected in the influence on the wheel alignment angle. Controlling the change of the wheel alignment angle within a reasonable range can maintain good dynamic characteristics of the vehicle during the movement. After reasonable simplification of the suspension, the wheel alignment angle of the suspension can be obtained through geometric relation and vector calculation of suspension.

Analysis and calculation of vector ${\vec{e}}_{k}$

In order to calculate the wheel alignment angle of the suspension, the vectors ${\vec{e}}_{k}$ and ${\vec{e}}_{hub}$ in Figure 5 need to be calculated.¹⁰ In order to calculate the vector ${\vec{e}}_{k}$ , the upper control arm, lower control arm and wheel hub of the suspension need to be taken into consideration. Because the mechanical structure of the double wishbone suspension can be seen as consisting of RSSR-SSP (R: revolute, S: spherical, P: prismatic) links,¹¹ the suspension in Figure 5 can perform two independent complete loop analysis of OBDAO and OBNHJO.

Figure 5.

Schematic diagram of double wishbone suspension.

First consider the RSSR closed-loop connection rod composed of OABCO. The ball hinge position at $J$ is considered to be fixed in a specific position. In order to analysis the RSSR connecting rod, the required structural parameters are: $l_{l}$ , $l_{u}$ , $l_{k}$ , $ξ$ , and $Φ$ . When these parameters are determined and given a set of lower control arm angle $θ$ , the RSSR connection rod can be used for kinematic analysis. The vector equation of this loop is: $| \overset{⇀}{BD} | = | \overset{⇀}{OD} - \overset{⇀}{OB} |$ . The vector on the right side of the equation can be written as:

\begin{matrix} | \overset{⇀}{OB} | = l_{l} \cos θ \overset{⇀}{i} + l_{l} \sin θ \overset{⇀}{k} \\ | \vec{OD} | = (d + l_{u} \cos χ \cos ξ + l_{u} \sin χ \sin Φ \sin ξ) \overset{⇀}{i} \\ + (h + l_{u} \cos χ \sin ξ - l_{u} \sin χ \sin Φ \cos ξ) \overset{⇀}{j} \\ + (q + l_{u} \sin χ \cos Φ) \overset{⇀}{k} \end{matrix}

(7)

Usually the displacement equation of the RSSR connecting rod is:

\begin{matrix} k_{1} + (k_{2} - k_{3} \cos θ) \cos χ - (k_{4} + k_{5} \sin θ - k_{6} \cos θ) \sin χ \\ - k_{7} \cos θ - k_{8} \sin θ = 0 \end{matrix}

(8)

among them,

{\begin{matrix} k_{1} = h^{2} + d^{2} + q^{2} + l_{l}^{2} + l_{u}^{2} - l_{k}^{2} \\ k_{2} = 2 h l_{u} \sin ξ + 2 d l_{u} \cos ξ \\ k_{3} = 2 l_{l} l_{u} \cos ξ \\ k_{4} = 2 h l_{u} \cos ξ \sin Φ - 2 q l_{u} \sin ξ \sin Φ \\ k_{5} = 2 l_{l} l_{u} \cos Φ \\ k_{6} = 2 l_{l} l_{u} \sin ξ \sin Φ \\ k_{7} = 2 d l_{l} \\ k_{8} = 2 q l_{l} \end{matrix}

(9)

For a given set of lower control arm angle $θ$ and combine formula (8), use the half-tangent formula can get the value of $χ$ :

χ = 2 \tan^{- 1} (\frac{- b \pm \sqrt{b^{2} - 4 ac}}{2 a})

(10)

among them,

{\begin{matrix} a = - k_{1} + k_{2} - k_{3} \cos θ + k_{7} \cos θ + k_{8} \sin θ \\ b = 2 k_{4} + 2 k_{5} \sin θ + 2 k_{6} \cos θ \\ c = - k_{1} - k_{2} - k_{3} \cos θ + k_{7} \cos θ + k_{8} \sin θ \end{matrix}

(11)

The position vector of the wheel knuckle can be defined as: $\overset{⇀}{BD} = l_{k} \overset{⇀}{e_{k}}$ , the unit vector of the wheel knuckle is defined as:

\overset{⇀}{e_{k}} = \frac{1}{l_{k}} [\begin{matrix} d + l_{u} \cos χ \cos ξ + l_{u} \sin χ \sin Φ \sin ξ - l_{l} \cos θ \\ h + l_{u} \cos χ \sin ξ - l_{u} \sin χ \sin Φ \cos ξ \\ q + l_{u} \sin χ \cos Φ - l_{l} \sin θ \end{matrix}]

(12)

The unit vector of the lower control arm can be defined as:

\overset{⇀}{e_{1}} = {[\cos θ 0 \sin θ]}^{T}

(13)

By formulas (10)–(13), a point position can be obtained when the $θ$ value of a set of lower control arm is given.

Analysis and calculation of vector ${\vec{e}}_{hub}$

In order to calculate vector ${\vec{e}}_{hub}$ , now discuss the SSP connection rod system. Due to the needs of analyzing, the OBNHJO connecting rod system was considered. For this analysis, the required parameters are: the initial coordinates of the spherical hinge points H and J. The length of the steering rod is $t_{r} = | \overset{⇀}{HJ} |$ , the effective steering arm length is $s_{a} = | \overset{⇀}{HB} \times \overset{⇀}{HD} | / | \overset{⇀}{DB} |$ . The change in toe angle of the double wishbone suspension due to the input angle $θ$ can be determined by constraining OBNHJO. The first constraint is that for each position of the mechanism, the length of the tie rod remains unchanged, so there is:

{t_{r}}^{2} = {(H_{x} - J_{x 0})}^{2} + {(H_{v} - J_{v 0})}^{2} + {(H_{z} - J_{z 0})}^{2}

(14)

The second constraint is that the length of the steering arm remains constant for each position of the mechanism, so there is:

{s_{a}}^{2} = {(N_{x} - H_{x})}^{2} + {(N_{y} - H_{y})}^{2} + {(N_{z} - H_{z})}^{2}

(15)

among them,

\overset{⇀}{ON} = l_{l} \overset{⇀}{e_{1}} + l_{a} \overset{⇀}{e_{k}}, l_{a} = \sqrt{{| \overset{⇀}{BH} |}^{2} - s_{a}^{2}}

(16)

The third constraint is that for each position of the mechanism, the steering arm vector and the steering knuckle vector are perpendicular to each other $(\overset{⇀}{HN} ⊥ \overset{⇀}{BN})$ . Because the scalar product of two vertical vectors is equal to zero, so: $\overset{⇀}{HN} \cdot \overset{⇀}{BN} = 0$ .

(l_{l} \overset{⇀}{e_{1}} + l_{a} \overset{⇀}{e_{k}} - (H_{x} \overset{⇀}{i} + H_{y} \overset{⇀}{j} + H_{z} \overset{⇀}{k})) \cdot (l_{a} \overset{⇀}{e_{k}}) = 0

(17)

In equation (17), it can get:

c_{3} H_{x} + c_{4} H_{y} + c_{5} H_{z} + c_{6} = 0

(18)

among them,

\begin{matrix} c_{3} = - l_{a} \overset{⇀}{e_{k_{x}}}, c_{4} = - l_{a} \overset{⇀}{e_{k_{y}}}, c_{5} = - l_{a} \overset{⇀}{e_{k_{z}}}, \\ c_{6} = {l_{a}}^{2} + l_{l} l_{a} (\overset{⇀}{e_{k_{x}}} \overset{⇀}{e_{1_{x}}} + \overset{⇀}{e_{k_{y}}} \overset{⇀}{e_{1_{y}}} + \overset{⇀}{e_{k_{z}}} \overset{⇀}{e_{1_{z}}}) \end{matrix}

(19)

In equations (17)–(19), the unknown parameters $H_{x}$ , $H_{y}$ , and $H_{z}$ are the new coordinates of the steering rod ball hinge point $H$ . These three equations can solve three unknowns, at this stage, the positions of all joints are determined. In the static position of the suspension system, the unit vector GW of the hub can be determined as:

\overset{⇀}{e_{hu b_{0}}} = {[\cos ε_{0} 0 - \sin ε_{0}]}^{T}

(20)

$ε_{0}$ is the position of the positive static surface of the wheel, as shown in Figure 6(a). The static position of the wheel center can be determined:

\overset{⇀}{l_{hu b_{0}}} = c_{0} {[0 - 1 0]}^{T} + l_{hub} \overset{⇀}{e_{hu b_{0}}}

(21)

Figure 6.

Vector transformation: (a) schematic diagram of wheel static surface position vector, (b) schematic diagram of Angle change of vector e_k, and (c) schematic diagram of the rotation of wheel knuckle relative to e_k.

In order to get the any position of the vector on the wheel knuckle, the rotation matrix can be used and take the following steps: starting from the initial position of the steering arm, determine the rotation matrix of the unit vector of the steering arm from $\overset{⇀}{e_{a 0}}$ to any position $\overset{⇀}{e_{a}}$ . Among them, the rotation of the wheel knuckle is considered as the superposition of two different rotations. Initially, the rotation of the joint around its axis is not considered. Since the unit matrix $\overset{⇀}{e_{k}}$ of any position of the mechanism can be calculated by formula (13), the rotation amount of $\overset{⇀}{e_{k}}$ can be determined by the outer product as:

λ = \sin^{- 1} [(\overset{⇀}{e_{k 0}} \times \overset{⇀}{e_{k}}) \cdot \overset{⇀}{e_{k λ}}]

(22)

$\overset{⇀}{e_{k λ}}$ is the vector on the rotation axis perpendicular to the plane formed by vectors $\overset{⇀}{e_{k}}$ and $\overset{⇀}{e_{k 0}}$ ,as shown in Figure 6(b):

\overset{⇀}{e_{k λ}} = (\overset{⇀}{e_{k 0}} \times \overset{⇀}{e_{k}}) / | \overset{⇀}{e_{k 0}} \times \overset{⇀}{e_{k}} |

(23)

The second rotation analysis requires a complete analysis of the rotation of wheel knuckle with respect to the unit vector $\overset{⇀}{e_{k}}$ , as shown in Figure 6(c):

λ_{p} = \sin^{- 1} [(\overset{⇀}{e_{a 0 λ}} \times \overset{⇀}{e_{a}}) \cdot \overset{⇀}{e_{k}}]

(24)

The matrix with a rotation angle of $θ$ around a certain axis in the direction of the unit vector $\overset{⇀}{u}$ is:

\begin{matrix} R (\overset{⇀}{u}, θ) = [\begin{matrix} \cos θ + u_{x}^{2} (1 - \cos θ) & u_{x} u_{y} (1 - \cos θ) - u_{z} \sin θ & u_{x} u_{z} (1 - \cos θ) + u_{y} \sin θ \\ u_{y} u_{x} (1 - \cos θ) + u_{z} \sin θ & \cos θ + u_{y}^{2} (1 - \cos θ) & u_{y} u_{z} (1 - \cos θ) - u_{x} \sin θ \\ u_{z} u_{x} (1 - \cos θ) - u_{y} \sin θ & u_{z} u_{y} (1 - \cos θ) + u_{x} \sin θ & \cos θ + u_{z}^{2} (1 - \cos θ) \end{matrix}] \end{matrix}

(25)

Because of the related rotation matrix has obtained, the unit vector $\overset{⇀}{e_{hub}}$ from the initial position to the new position can be determined as:

\overset{⇀}{e_{hub}} = R (\overset{⇀}{e_{k}}, λ_{p}) R (\overset{⇀}{e_{k λ}}, λ) \overset{⇀}{e_{hu b_{0}}}

(26)

Calculation of wheel alignment angle

In order to analysis and calculate the wheel alignment angle of the double wishbone suspension, the double wishbone suspension is simplified to obtain a schematic diagram as shown in Figure 7. The camber angle, caster angle, kingpin angle, and toe angle in Figure 7 can all be calculated by vectors.¹²

Figure 7.

Schematic diagram of wheel alignment angle of double wishbone suspension.

In Figure 7, $\vec{k}$ represents the unit vector in the z-axis direction, $\vec{j}$ represents the unit vector in the y-axis direction, ${\vec{e}}_{k}$ represents the unit vector of the wheel hub, ${\vec{e}}_{hub}$ represents the unit vector of the wheel axis, ${\vec{e}}_{hubp}$ represents the unit vector obtained by projecting the vector ${\vec{e}}_{hub}$ in the XY plane, ${\vec{e}}_{caster}$ is the unit vector perpendicular to vector ${\vec{e}}_{hubp}$ in the XY plane, $α$ is the toe angle of the wheel, $φ$ is the camber angle, $τ$ is caster angle, $σ$ is kingpin angle. The toe angle of the suspension system can be calculated from the vectors ${\vec{e}}_{hubp}$ and $\vec{j}$ :

α = π / 2 - \cos^{- 1} (\vec{e_{hubp}} \cdot \vec{j})

(27)

The camber angle shown in Figure 7 can be calculated by vectors ${\vec{e}}_{hub}$ and $\vec{k}$ :

φ = \cos^{- 1} (\overset{⇀}{e_{hub}} \cdot \overset{⇀}{k}) - π / 2

(28)

The caster angle shown in Figure 7 can be calculated by vectors ${\vec{e}}_{k}$ and ${\vec{e}}_{caster}$ :

τ = \cos^{- 1} (\overset{⇀}{e_{k}} \cdot \overset{⇀}{e_{caster}}) - π / 2

(29)

among them,

\begin{matrix} {e_{caster}}^{⇀} = \\ {[\mp \frac{e_{huby}}{e_{hubx}} \sqrt{\frac{1}{{(\frac{e_{huby}}{e_{hubx}})}^{2} + 1}} \pm \sqrt{\frac{1}{{(\frac{e_{huby}}{e_{hubx}})}^{2} + 1}} 0]}^{T} \end{matrix}

(30)

The kingpin angle of the suspension system can be calculated by vector ${\vec{e}}_{k}$ and vector ${\vec{e}}_{hubp}$ :

σ = \cos^{- 1} (\overset{⇀}{e_{k}} \cdot \overset{⇀}{e_{hubp}}) - π / 2

(31)

among them,

\overset{⇀}{e_{hubp}} = {[e_{hub p_{x}} e_{hub p_{y}} 0]}^{T} / \sqrt{e_{hub p_{x}}^{2} + e_{hub p_{y}}^{2}}

(32)

After analyzing the influence of the suspension motion on the wheel alignment angle, the opposite travel simulation on the front suspension and the rear suspension was carried out, the changes of toe angle and camber angle with wheel travel as shown in Figures 8 and 9 were obtained.

Figure 8.

Variation diagram of toe angle with wheel travel.

Figure 9.

Variation diagram of camber angle with wheel travel.

In Figure 8, during the wheel travel downward, the toe angle of the front suspension wheel (the inner wheel when turning) increases in the positive direction, and the toe angle of the rear suspension wheel (the inner wheel when turning) increases in the negative direction. The toe angle of the corresponding outer wheel of the front suspension increases negatively, and the toe angle of the outer wheel of the rear suspension increases positively. For the front suspension, the toe angle of the inner wheel increases positively, and the toe angle of the outer wheel increases negatively, which will increase the understeer characteristics of the vehicle. For the rear suspension, the toe angle of the inner wheel increases negatively, and the toe angle of the outer wheel increases positively, which will also increase the understeer characteristics of the vehicle.

In Figure 9, when the wheel travels downward, the camber angle of the front and rear suspension wheel (the inner wheel when turning) increases in the positive direction. When the wheel travels upward, the camber angle of the wheel (the outer wheel when turning) of the front suspension and the rear suspension both increases in the negative camber direction. In terms of adhesion, the front suspension’s and the rear suspension’s outer wheel camber angle increase in the negative camber direction, the tire’s adhesion will be improved. For the perspective of understeer characteristics of the vehicle, when the vehicle is turning, the camber angle of the front suspension’s outer wheel increases to the negative camber direction, the understeer characteristic of the vehicle will decrease, and the camber angle of the rear suspension’s outer wheel increases to the negative camber direction, the understeer characteristics of the vehicle will increase.

Statics analysis of double wishbone suspension

When the wheels of the double wishbone suspension are excited by the road and jump, the hinge point of the suspension link and the link rod will move with the travel of the wheel, so the position of each hinge point will change accordingly. Through the theory of coordinate transformation, after selecting a reference point, the coordinates of other hinge points can be calculated accordingly. The force balance and torque balance equations can be used to calculate the force of the vehicle suspension hinge points and bushing.

Displacement analysis of double wishbone suspension’s connecting point

The schematic diagram of the double wishbone suspension is shown in Figure 10. In order to analysis the force and moment of suspension’s ball joint and bushing, firstly, the position transformation of each connection point of the double wishbone suspension is analyzed.¹³ The double wishbone suspension is composed of upper control arm, lower control arm, wheel knuckle, springs, shock absorbers and tie rods. The four bushes A₁ and A₂, A₃ and A₄ in Figure 10 connect the upper and lower control arms to the vehicle body respectively. The ball joint D and the ball joint G respectively connect the wheel knuckle with the upper and lower control arms. The two ends of the coil spring and the shock absorber are installed between the lower control arm and the vehicle’s body, which are installed between the ball joints C and E, respectively.

Figure 10.

Schematic diagram of double wishbone suspension system.

If an object rotates from one position to another without displacement, the coordinates of the object in the global coordinate system $o - xyz$ can be determined by the coordinate system fixed on the rigid body. Take the upper control arms A₁D and A₂D of the double wishbone suspension as an example, assuming that the rotation angles around the $x$ , $y$ , and $z$ axis are $η_{1}$ , $ρ_{1}$ , and $ζ_{1}$ respectively, the new coordinates ${A_{1}^{'}}$ , ${A_{2}^{'}}$ and ${D^{'}}$ of points A₁, A₂ and D on the rigid body A₁D and A₂D in the coordinate system fixed to itself can be obtained.

{A_{1}^{'}} = [X (η_{1})] [Y (ρ_{1})] [Z (ζ_{1})] {A_{1}^{0}}

(33)

{A_{2}^{'}} = [X (η_{1})] [Y (ρ_{1})] [Z (ζ_{1})] {A_{2}^{0}}

(34)

{D^{'}} = [X (η_{1})] [Y (ρ_{1})] [Z (ζ_{1})] {D^{0}}

(35)

$[X (η_{1})]$ , $[Y (ρ_{1})]$ , and $[Z (ζ_{1})]$ are the $(3 \times 3)$ rotation matrices of the upper control arm around the $x$ , $y$ , and $z$ axis on their own coordinate system and the rotation angles are $η_{1}$ , $ρ_{1}$ , and $ζ_{1}$ . ${A_{1}^{0}}$ , ${A_{2}^{0}}$ , and ${D^{0}}$ are the initial coordinates of points $A_{1}$ , $A_{2}$ , and $D$ .

When the rigid body rotates and moves at the same time, select point $D$ as the reference point, then the new coordinates ${A_{1}^{'}}$ and ${A_{2}^{'}}$ of points $A_{1}$ and $A_{2}$ in the global coordinate system can be obtained by formula (33) subtracting formula (35) and formula (34) subtracting formula (35):

\begin{matrix} {\begin{matrix} {A_{1}^{'}} \\ {A_{2}^{'}} \end{matrix}} = [\begin{matrix} [R (η_{1}, ρ_{1}, ζ_{1})], 0 \\ 0, [R (η_{1}, ρ_{1}, ζ_{1})] \end{matrix}] \\ {\begin{matrix} {{A_{1}^{0}} - {D^{0}}} \\ {{A_{2}^{0}} - {D^{0}}} \end{matrix}} + {\begin{matrix} {D^{'}} \\ {D^{'}} \end{matrix}} \end{matrix}

(36)

The 0 in formula (36) represents a zero matrix of $3 \times 3$ , ${A_{1}^{0}}$ , ${A_{2}^{0}}$ and ${D^{0}}$ are the initial coordinates of the connection points $A_{1}$ , $A_{2}$ and $D$ in the global coordinate system; ${D^{'}}$ is the new coordinate value of hinge point $D$ , matrix $[R (η_{1}, ρ_{1}, ζ_{1})]$ is a rotation matrix:

[R (η_{1}, ρ_{1}, ζ_{1})] = [X (η_{1})] [Y (ρ_{1})] [Z (ζ_{1})]

(37)

In the same way, taking point $G$ as the reference point, the new coordinates ${A_{3}^{'}}$ , ${A_{4}^{'}}$ , and ${C^{'}}$ of the connection points $A_{3}$ , $A_{4}$ , and $C$ on the lower control arm can be obtained:

\begin{matrix} {\begin{matrix} {A_{3}^{'}} \\ {A_{4}^{'}} \\ {C^{'}} \end{matrix}} = [\begin{matrix} [R (η_{2}, ρ_{2}, ζ_{2})], 0 0 \\ 0, [R (η_{2}, ρ_{2}, ζ_{2})], 0 \\ 0 0 [R (η_{2}, ρ_{2}, ζ_{2})] \end{matrix}] \\ {\begin{matrix} {A_{3}^{0}} - {G^{0}} \\ {A_{4}^{0}} - {G^{0}} \\ {C^{0}} - {G^{0}} \end{matrix}} + {\begin{matrix} {G^{'}} \\ {G^{'}} \\ {G^{'}} \end{matrix}} \end{matrix}

(38)

${A_{3}^{0}}$ , ${A_{4}^{0}}$ , ${C_{0}}$ , and ${G_{0}}$ are the initial positions of the hinge points $A_{3}$ , $A_{4}$ , $C$ , and $G$ . ${G^{'}}$ is the new coordinate of the hinge point $G$ ; Matrix $[R (η_{2}, ρ_{2}, ζ_{2})]$ is a rotation matrix:

[R (η_{1}, ρ_{1}, ζ_{1})] = [X (η_{2})] [Y (ρ_{2})] [Z (ζ_{2})]

(39)

$[X (η_{2})]$ , $[Y (ρ_{2})]$ , and $[Z (ζ_{2})]$ are the $3 \times 3$ rotation matrices of the lower control arm around the $x$ , $y$ , and $z$ axis on its own coordinate system and the rotation angles are $η_{2}$ , $ρ_{2}$ , $ζ_{2}$ . In the same way, with point $D$ as the reference point, the new coordinates ${G^{'}}$ , ${H^{'}}$ , and ${P_{1}^{'}}$ of the hinge points $G$ , $H$ , and $P_{1}$ can be obtained:

\begin{matrix} {\begin{matrix} {G^{'}} \\ {H^{'}} \\ {P_{1}^{'}} \end{matrix}} = [\begin{matrix} [R (η_{3}, ρ_{3}, ζ_{3})], 0, 0 \\ 0, [R (η_{3}, ρ_{3}, ζ_{3})], 0 \\ 0, 0, [R (η_{3}, ρ_{3}, ζ_{3})] \end{matrix}] \\ {\begin{matrix} {{G^{0}} - {D^{0}}} \\ {{H^{0}} - {D^{0}}} \\ {{P_{1}^{'}} - {D^{0}}} \end{matrix}} + {\begin{matrix} {D^{'}} \\ {D^{'}} \\ {D^{'}} \end{matrix}} \end{matrix}

(40)

$H^{0}$ and $P_{1}^{0}$ are the initial coordinates of points $H$ and $P_{1}$ , and the matrix $[R (η_{3}, ρ_{3}, ζ_{3})]$ can be calculated as follows:

[R (η_{3}, ρ_{3}, ζ_{3})] = [X (η_{3})] [Y (ρ_{3})] [Z (ζ_{3})]

(41)

$[X (η_{3})]$ , $[Y (ρ_{3})]$ , and $[Z (ζ_{3})]$ are the $3 \times 3$ rotation matrices of the wheel knuckle around the $x$ , $y$ , and $z$ axis on its own coordinate system and the rotation angles are $η_{3}$ , $ρ_{3}$ , and $ζ_{3}$ . According to equations (36), (38), (40), when external force is applied to the suspension, the new position of each connection point and bushing can be obtained. The displacement of the bushing and spring can also be easily estimated.

Calculate the forces and moments acting on the bushing

If a static external force ${F_{P 1}^{g}}$ is applied to the contact point between the tire and the ground, the bushing will rotate and translate. The displacement of the bushing $S$ ( $S$ means $A_{1}$ , $A_{2}$ , $A_{3}$ , and $A_{4}$ ) can be calculated:

{Q_{S}^{g}} = {\begin{matrix} {S^{'}} - {S_{0}} \\ {ϕ_{S}} \end{matrix}} + {\begin{matrix} {Δ S_{0}^{g}} \\ {Δ ϕ_{S 0}^{g}} \end{matrix}}

(42)

Superscript $g$ means in the global coordinate system. ${S^{'}}$ represents the translational position of the bush after deformation. ${ϕ_{S}}$ is the rotation angle of the bush $S$ in the global coordinate system when the external force is input. ${Δ S_{0}^{g}}$ and ${Δ ϕ_{S 0}^{g}}$ in equation (42) are the initial displacement and initial angle of the bushing in the global coordinate system under the action of the vehicle’s own weight. In the global coordinate system, the force and moment applied to the bush due to the displacement of bush $S$ are:

{\begin{matrix} {F_{S}^{g}} \\ {M_{S}^{g}} \end{matrix}} = [\begin{matrix} [K_{St}^{g}], 0 \\ 0, [K_{Sr}^{g}] \end{matrix}] {Q_{S}^{g}}

(43)

${F_{S}^{g}}$ and ${M_{S}^{g}}$ are the force matrix and moment matrix of $3 \times 1$ , respectively.

{\begin{matrix} [K_{St}^{g}] = {[R_{S}]}^{T} [K_{St}^{l}] [R_{S}] \\ [K_{Sr}^{g}] = {[R_{S}]}^{T} [K_{Sr}^{l}] [R_{S}] \end{matrix}

(44)

$T$ represents the transpose of the matrix. $[K_{St}^{g}]$ and $[K_{Sr}^{g}]$ are the translational movement stiffness matrix and the rotation stiffness matrix respectively of the rubber bushing $S$ in the global coordinate system. $[K_{St}] = diag (k_{Stx}, k_{Sty}, k_{Stz})$ and $[K_{Sr}] = diag (k_{Srx}, k_{Sry}, k_{Srz})$ are the translational stiffness matrix and rotational stiffness matrix of the rubber bushing $S$ in the local coordinate system. $[R_{S}]$ is the $3 \times 3$ transformation matrix of the bushing $S$ from the local coordinate system to the global coordinate system. In the global coordinate system, the force and moment applied to the upper control arm by the bush $S_{u}$ ( $S_{u}$ refers to the bush A₁ and the bush A₂) are:

{\begin{cases} {\begin{cases} {R F_{S_{u}}^{g}} \\ {R M_{S_{u}}^{g}} \end{cases}} = [\begin{array}{l} I, & 0 \\ [{\tilde{R}}_{S_{u} D}], & I \end{array}] {\begin{cases} - {F_{S_{u}}^{g}} \\ - {M_{S_{u}}^{g}} \end{cases}} \\ {\begin{cases} {R F_{S_{l}}^{g}} \\ {R M_{S_{l}}^{g}} \end{cases}} = [\begin{array}{l} I & 0 \\ [{\tilde{R}}_{S_{l} G}], & I \end{array}] {\begin{cases} - {F_{S_{l}}^{g}} \\ - {M_{S_{l}}^{g}} \end{cases}} \end{cases}

(45)

$I$ is the unit matrix of $(3 \times 3)$ , $[{\tilde{R}}_{S_{u} D}]$ is an oblique symmetric matrix of the position vector of the bush $S_{u}$ in the global coordinate system, and point $D$ is the reference point. $[{\tilde{R}}_{S_{l} G}]$ is an oblique symmetric matrix of the position vector of the bush $S_{l}$ in the global coordinate system. The oblique symmetric matrix $[{\tilde{R}}_{S_{u} D}]$ and $[{\tilde{R}}_{S_{l} G}]$ are:

{\begin{matrix} [{\tilde{R}}_{S_{u} D}] = [\begin{matrix} 0 & - z_{S_{u} D}^{g} & y_{S_{u} D}^{g} \\ z_{S_{u} D}^{g} & 0 & - x_{S_{u} D}^{g} \\ - y_{S_{u} D}^{g} & x_{S_{u} D}^{g} & 0 \end{matrix}] \\ [{\tilde{R}}_{S_{l} D}] = [\begin{matrix} 0 & - z_{S_{l} D}^{g} & y_{S_{l} D}^{g} \\ z_{S_{l} D}^{g} & 0 & - x_{S_{l} D}^{g} \\ - y_{S_{l} D}^{g} & x_{S_{l} D}^{g} & 0 \end{matrix}] \end{matrix}

(46)

$x_{S_{u} D}^{g}$ , $y_{S_{u} D}^{g}$ , $z_{S_{u} D}^{g}$ are the displacement of the bush $S_{u}$ along the $x$ , $y$ , and $z$ directions in the global coordinate system. $x_{S_{l} D}^{g}$ , $y_{S_{l} D}^{g}$ , $z_{S_{l} D}^{g}$ are the displacement of the bush $S_{l}$ along the $x$ , $y$ , and $z$ directions in the global coordinate system.

The upper control arm is in a balanced state under the action of the reaction force and the reaction moment of the bushing, the spring, and the ball joint of the wheel knuckle. Therefore, the sum of each force component of all forces, as well as the sum of each force and each moment component of the connection point D in the global coordinate system is zero. According to this principle, the static balance equation of the upper control arm can be obtained and expressed as:

{\begin{matrix} {{RF}_{A 1}^{g}} \\ {{RM}_{A 1}^{g}} \end{matrix}} + {\begin{matrix} {{RF}_{A 2}^{g}} \\ {{RM}_{A 2}^{g}} \end{matrix}} + {\begin{matrix} {F_{D}^{g}} \\ 0 \end{matrix}} = {0}

(47)

${F_{D}^{g}}$ is the combined force applied to the upper control arm from the wheel knuckle and the connection point D. Similar to obtaining the static balance equation of the upper control arm, the static balance equation of the lower control arm can be obtained, expressed as:

\begin{matrix} {\begin{matrix} {{RF}_{A 3}^{g}} \\ {{RM}_{A 3}^{g}} \end{matrix}} + {\begin{matrix} {{RF}_{A 4}^{g}} \\ {{RM}_{A 4}^{g}} \end{matrix}} \\ + {\begin{matrix} {F_{C}^{g}} \\ {M_{C}^{g}} \end{matrix}} + {\begin{matrix} {F_{G}^{g}} \\ 0 \end{matrix}} = {0} \end{matrix}

(48)

${F_{G}^{g}}$ is the reaction force applied to the lower control arm from the wheel knuckle G. ${F_{C}^{g}}$ is the reaction force of the spring acting on the lower control arm. Under the action of the external load of the upper and lower control arms, tie rods and the contact point between the tire and ground, the wheel knuckle is in a balanced state,¹⁴ and the static balance equation is obtained:

{\begin{matrix} {F_{P 1}^{g}} \\ {M_{P 1}^{g}} \end{matrix}} - {\begin{matrix} {F_{D}^{g}} \\ {M_{D}^{g}} \end{matrix}} - {\begin{matrix} {F_{H}^{g}} \\ {M_{H}^{g}} \end{matrix}} - {\begin{matrix} {F_{G}^{g}} \\ 0 \end{matrix}} = {0}

(49)

${F_{P 1}^{g}}$ is the external force input to the tire from the road, ${M_{P 1}^{g}}$ is the moment of ${F_{P 1}^{g}}$ about the connection point $D$ . ${F_{H}^{g}}$ is the reaction force acting on the wheel knuckle from the tie rod, expressed as:

{F_{H}^{g}} = ν_{s 1} {R_{JH}^{*}}

(50)

$ν_{s 1}$ is the scale factor, ${R_{JH}^{*}}$ is the vector from point $H$ to point $J$ . ${M_{H}^{g}}$ is the moment of force ${F_{H}^{g}}$ relative to point $G$ . The length, displacement and initial length of the spring after deformation satisfy the following formula:

L_{EC}^{'} + Δ L_{EC} - L_{EC} = 0

(51)

$L_{EC}^{'}$ is the new length of the spring after deformation, $L_{EC}$ is the initial length of the spring, $Δ L_{EC}$ is the deformation of the spring:

Δ L_{EC} = (ν_{s 2} L_{EC}^{'} - ν_{s 20} L_{EC}) / K_{d}

(52)

$ν_{s 20}$ is the initial scale factor of the spring force, $K_{d}$ is the stiffness of the spring. Due to the length of the tie rod does not change during the movement of the suspension, there is:

L_{JH}^{'} - L_{JH} = 0

(53)

$L_{JH}^{'}$ and $L_{JH}$ are the length of the tie rod after and before the suspension change, respectively. Through the above analysis, the overall balance equation can be obtained as:

{\begin{matrix} {\begin{matrix} {{RF}_{A 1}^{g}} \\ {{RM}_{A 1}^{g}} \end{matrix}} + {\begin{matrix} {{RF}_{A 2}^{g}} \\ {{RM}_{A 2}^{g}} \end{matrix}} + {\begin{matrix} {F_{D}^{g}} \\ 0 \end{matrix}} \\ {\begin{matrix} {{RF}_{A 3}^{g}} \\ {{RM}_{A 3}^{g}} \end{matrix}} + {\begin{matrix} {{RF}_{A 4}^{g}} \\ {{RM}_{A 4}^{g}} \end{matrix}} + {\begin{matrix} {F_{C}^{g}} \\ {M_{C}^{g}} \end{matrix}} + {\begin{matrix} {F_{G}^{g}} \\ 0 \end{matrix}} \\ {\begin{matrix} {F_{P 1}^{g}} \\ {M_{P 1}^{g}} \end{matrix}} - {\begin{matrix} {F_{D}^{g}} \\ {M_{D}^{g}} \end{matrix}} - {\begin{matrix} {F_{H}^{g}} \\ {M_{H}^{g}} \end{matrix}} - {\begin{matrix} {F_{G}^{g}} \\ 0 \end{matrix}} \\ {\begin{matrix} L_{EC}^{'} + Δ L_{EC} - L_{EC} \\ L_{JH}^{'} - L_{JH} \end{matrix}} \end{matrix}} = {0}

(54)

Take the lateral force in the interval $[- 3000, 3000]$ , take a value into the formula (54), the equation (54) can be solved by using the Newton iteration method, the force and moment situations shown in Figures 11 and 12 can be obtained.

Figure 11.

The force in the x direction of the rubber bushing.

Figure 12.

The moment of the rubber bushing around the x-axis.

Statics analysis of multi-link suspension

When an external load acts on the contact point between the tire and the ground, the wheels, wheel knuckle and connection rods will move to new position due to the elastic deformation of the bushes and springs in the suspension. In this way, a new static equilibrium is reached. The new position of each bush in this load can then be calculated.

Displacement analysis of multi-link suspension’s connection point

The schematic diagram of the multi-link suspension is shown in Figure 13.

Figure 13.

Schematic diagram of multi-link suspension.

For the wheel knuckle, when point $B_{5}$ is selected as the reference point, the new position ${B_{i}^{'}}$ for other points $B_{i}$ ( $i = 5 ~ 8$ ) can be calculated:

{B_{i}^{'}} = [X (η_{4})] [Y (ρ_{4})] [Z (ζ_{4})] {{B_{i}^{0}} - {P_{2}^{0}}} + {P_{2}^{'}}

(55)

${P_{2}^{0}}$ and ${B_{i}^{0}}$ are the initial coordinates of points $P_{2}$ and ${B_{i}}$ respectively; ${P_{2}^{'}}$ is the new coordinate of point $P_{2}$ ; $[X (η_{4})]$ , $[Y (ρ_{4})]$ , and $[Z (ζ_{4})]$ are the $(3 \times 3)$ rotation matrices when the wheel knuckle rotates around the $x$ , $y$ , and $z$ axis respectively, and the rotation angles are $η_{4}$ , $ρ_{4}$ , and $ζ_{4}$ .

For connection rod $A_{i} B_{i}$ ( $i = 5 ~ 8$ ), when point $A_{i}$ is selected as the reference point, the new coordinates of point $B_{i}$ in the global coordinates system is:

{B_{i}^{'}} = [X (η_{i})] [Y (ρ_{i})] [Z (ζ_{i})] {{B_{i}^{0}} - {A_{i}^{0}}} + {A_{i}^{'}}

(56)

Where ${A_{i}^{0}}$ and ${B_{i}^{0}}$ are the initial coordinates of points $A_{i}$ and $B_{i}$ respectively; ${A_{i}^{'}}$ is the new position of point $A_{i}$ ; $[X (η_{i})]$ , $[Y (ρ_{i})]$ , and $[Z (ζ_{i})]$ are the $(3 \times 3)$ rotation matrix when connecting rod $A_{i} B_{i}$ rotates around $x$ , $y$ , and $z$ axis respectively, and the rotation angles are $η_{i}$ , $ρ_{i}$ , and $ζ_{i}$ .

According to equations (55) and (56), when external force acts on the suspension, the new position and displacement of each bush can be obtained.

Calculate the forces and moments acting on the rubber bushing

When an external force is applied to the wheel knuckle, the bushing will rotate and translate.¹⁵ For the bushing $w$ ( $w$ means $A_{5}$ , $A_{6}$ , $A_{7}$ , and $A_{8}$ ) connecting the connection rod and the body, the displacement ${Q_{w}^{g}}$ in the global coordinate system is:

{Q_{w}^{g}} = {\begin{matrix} {w^{'}} - {w_{0}} \\ {ϕ_{w}} \end{matrix}} + {\begin{matrix} {Δ w_{0}^{g}} \\ {Δ ϕ_{w 0}^{g}} \end{matrix}}

(57)

The superscript $g$ in equation (57) means in the global coordinate system; ${ϕ_{w}}$ is the rotation angle matrix of the bushing $w$ in the global coordinate system; ${Δ w_{0}^{g}}$ and ${Δ ϕ_{w 0}^{g}}$ are the initial displacement and rotation angle of bushing $w$ in the global coordinate system under the weight of the vehicle. The force ${F_{w}^{g}}$ and moment ${M_{w}^{g}}$ due to the displacement of the bushing are:

{\begin{matrix} {F_{w}^{g}} \\ {M_{w}^{g}} \end{matrix}} = [\begin{matrix} [K_{wt}^{g}], & 0 \\ 0, & [K_{wr}^{g}] \end{matrix}] {Q_{w}^{g}}

(58)

0 is the zero matrix of $(3 \times 3)$ .

{\begin{matrix} [K_{wt}^{g}] = {[R_{w}]}^{T} [K_{wt}^{l}] [R_{w}] \\ [K_{wr}^{g}] = {[R_{w}]}^{T} [K_{wr}^{l}] [R_{w}] \end{matrix}

(59)

$T$ represents the transpose of the matrix. $[K_{wt}^{g}]$ and $[K_{wr}^{g}]$ are the translational movement stiffness matrix and the rotation stiffness matrix respectively of the bush $w$ in the global coordinate system. $[K_{wt}^{l}] = diag (k_{wtx}, k_{wty}, k_{wtz})$ and $[K_{wr}^{l}] = diag (k_{wrx}, k_{wry}, k_{wrz})$ are the translational stiffness matrix and rotational stiffness matrix of the bush $w$ in the local coordinate system. $[R_{w}]$ is the $3 \times 3$ transformation matrix of the bushing $w$ from the local coordinate system to the global coordinate system. According to Newton’s third law and equation (58), the reaction force and moment applied by the rubber bushing to the control arm in the global coordinate system are:

{\begin{matrix} {{RF}_{w}^{g}} \\ {{RM}_{w}^{g}} \end{matrix}} = - {\begin{matrix} {F_{w}^{g}} \\ {M_{w}^{g}} \end{matrix}} = [\begin{matrix} - [K_{wt}^{g}], & 0 \\ 0, & - [K_{wr}^{g}] \end{matrix}] {Q_{w}^{g}}

(60)

Mark the forces and moments on the hinge points $B_{5}$ , $B_{6}$ , $B_{7}$ and $B_{8}$ as ${F_{Bi}}$ and ${M_{Bi}}$ , then for the rigid body $A_{5} B_{5}$ , $A_{6} B_{6}$ , $A_{7} B_{7}$ , and $A_{8} B_{8}$ , the reaction force and moment at the $B_{i}$ $(i = 5 ~ 8)$ end are denoted as ${R F_{Bi}}$ and ${R M_{Bi}}$ .

Force balance and moment balance analysis

Under the action of the reaction force and moment of the hinge point $B_{i}$ and the external load, the wheel knuckle is in a balanced state. Therefore, the sum of each reaction force (RF), reaction moment (RM) and external force applied to the wheel knuckle is zero in the global coordinate system. The static balance equation of the wheel knuckle is expressed as:

{F_{P 2}^{g}} + {\begin{matrix} F_{B 0}^{g} \\ M_{B 0}^{g} \end{matrix}} + {\begin{matrix} \sum_{j = 5}^{8} {{RF}_{Bj}^{g}} \\ \sum_{j = 5}^{8} {{RM}_{Bj}^{g}} \end{matrix}} = {0}

(61)

${F_{P 2}^{g}}$ is the external force applied to the wheel at point $P_{2}$ .The rod is in equilibrium under the reaction force and moment of bushing $A_{i}$ and hinge point $B_{i}$ , therefore, in the global coordinate system, the sum of each force generated by the elastic deformation of the bush and each moment generated by the force relative to point $A_{i}$ is zero. Therefore, the static balance equation of rod $A_{i} B_{i}$ (i = 5–8) can be obtained and expressed as:

{\begin{matrix} {{RF}_{w}^{g}} \\ {{RM}_{w}^{g}} \end{matrix}} + {\begin{matrix} {{RF}_{Bi}^{g}} \\ {{RM}_{Bi}^{g}} \end{matrix}} = {0}

(62)

The initial length, deformed length and displacement of the spring satisfy the following formula:

L_{A_{0} B_{0}}^{'} + Δ L_{A_{0} B_{0}} - L_{A_{0} B_{0}} = 0

(63)

$L_{A_{0} B_{0}}^{'}$ is the length of the spring after deformation; $L_{A_{0} B_{0}}$ is the initial length of the spring; $Δ L_{A_{0} B_{0}}$ is the deformation of the spring:

Δ L_{A_{0} B_{0}} = (ν_{s 3} L_{A_{0} B_{0}}^{'} - ν_{s 30} L_{A_{0} B_{0}}) / K_{s}

(64)

$ν_{s 30}$ is the initial scale factor of the spring force, $K_{s}$ is the stiffness of the spring.

On the basis of the above analysis, the balance equations of the suspension system can be obtained. By combining the balance equation of the wheel knuckle (equation (61)), the balance equation of the rod (equation (62) and the constraint equation of spring (equation (63)), the equation of the multi-link suspension is expressed as:

{\begin{matrix} {F_{P 2}^{g}} + {\begin{matrix} F_{B 0}^{g} \\ M_{B 0}^{g} \end{matrix}} + {\begin{matrix} \sum_{j = 5}^{8} {{RF}_{Bj}^{g}} \\ \sum_{j = 5}^{8} {{RM}_{Bj}^{g}} \end{matrix}} \\ {\begin{matrix} {{RF}_{w}^{g}} \\ {{RM}_{w}^{g}} \end{matrix}} + {\begin{matrix} {{RF}_{Bi}^{g}} \\ {{RM}_{Bi}^{g}} \end{matrix}} \\ {L_{A_{0} B_{0}}^{'} + Δ L_{A_{0} B_{0}} - L_{A_{0} B_{0}}} \end{matrix}} = {0}

(65)

Figure 14 shows the change of the force on the rear multi-link suspension bushing with the lateral force on the wheels. Figure 15 shows the change of the torque on the rear multi-link suspension bushing with the lateral force on the wheels.

Figure 14.

The change of the bush’s force in the x-direction with the lateral force of the wheel.

Figure 15.

The torque of the bushing around the x-axis changes with the lateral force of the wheel.

Analysis of the influence of suspension kinematics and compliance characteristic on vehicle dynamics

The actual deformation of the rubber bush during the movement of the suspension is more complicated, including the axial movement of the inner ring (metal inner surface) relative to the outer ring (metal outer surface), the movement of the inner and outer rings in different radial directions, and the rotation of the inner and outer rings around the axis etc. The suspension kinematics characteristic is the change characteristic of the wheel alignment angles caused by the geometric movement of the suspension. The kinematics characteristic of the suspension include wheel bouncing in the same direction and opposite direction. The kinematics characteristics of the suspension are analyzed in Figures 8 and 9. The characteristic of the tire force and the change of the wheel alignment angle caused by the deformation of the elastic elements such as the bushing in the suspension is called the compliance characteristic of the suspension. The compliance characteristic of the suspension includes the characteristic of the tires subjected the same direction lateral force, reverse lateral force, longitudinal force, same direction align torque, and reverse align torque. The impact of the bushing on the dynamic characteristic of the vehicle is mainly in two aspects, one is that due to the deformation of the bushing, when the tire is subjected to lateral force, longitudinal force and align torque, it will cause the wheel to produce additional steering angle, which will affect the steering characteristic of the vehicle. At the same time, the bushing will affect the force on the vehicle wheels, that is, it affects the dynamic characteristic of the vehicle by influencing the tire force.

The influence of suspension motion and bushing on the change of toe Angle and camber Angle

From the analysis of the front two parts, it can be seen that the force of the vehicle wheel will cause the deformation of the suspension bushing, and the deformation of the bushing will cause the change of the wheel alignment angle. Based on the force analysis about front double wishbone suspension and the rear multi-link suspension, the compliance characteristic simulation analysis of the double wishbone suspension and the rear multi-link suspension of the vehicle is carried out, the changes of the toe angle and the camber angle of the wheel with the lateral force and the longitudinal force as shown in Figures 16 to 19 are obtained.

Figure 16.

The toe angle and camber angle of the front double wishbone suspension vary with the lateral force of the wheel.

Figure 17.

The toe angle and camber angle of the front double wishbone suspension vary with the longitudinal force of the wheel.

Figure 18.

Changes in toe angle and camber angle of rear multi-link suspension with wheel’s lateral force.

Figure 19.

Changes in toe angle and camber angle of rear multi-link suspension with wheel’s longitudinal force.

Figure 16 shows the change of toe angle and camber angle of the front double wishbone suspension with the lateral force of the wheel. In Figure 16, both the toe angle and the camber angle increase in the positive direction as the reverse lateral force increases. The toe angle and camber angle both increase in the negative direction as the positive lateral force increases. Figure 17 shows the change of toe angle and camber angle of the front double wishbone suspension with the longitudinal force of the wheel. In Figure 17, when the longitudinal force received by the wheels is the braking force, the toe angle increases in the negative direction with the increase of the braking force, and the camber angle increases in the positive direction with the increase of the braking force. When the longitudinal force received by the wheels is the driving force, the toe angle increases in the positive direction as the driving force increases, and the camber angle decreases in the negative direction as the driving force increases.

Figure 18 shows the change of the toe angle and camber angle of the rear multi-link suspension with the lateral force received by the wheels. In Figure 18, when the lateral force received is the reverse lateral force, the toe angle of the wheels changes less, and the camber angle increases in the positive direction as the reverse lateral force increases. When the lateral force received by the vehicle is a positive lateral force, the toe angle of the wheels changes less, and the camber angle increases in the negative direction. Figure 19 shows the change of the toe angle and camber angle of the rear multi-link suspension wheels when subjected to longitudinal force. In Figure 19, whether the wheel receives the driving force or the braking force, the toe angle changes little. For the camber angle, when the wheel receives the braking force, the camber angle increases in the negative direction, when the wheel receives the driving force, the camber angle of the wheel increases in the positive direction.

Effect of toe angle change caused by suspension movement and bushing deformation on tire side slip stiffness

When there is toe angle of the wheel, the side slip angle of the outer wheel is $α_{a} = α + δ_{0}$ , and the side slip angle of the inner wheel is $α_{a} = α - δ_{0}$ , there will be a deviation as shown in Figure 20. At this time, the sum of the side slip stiffness at point 1 is too large. The reason is: due to the higher wheel load at point 2, the side slip stiffness increases more than at point 3. The result is the toe angle reduces the axle’s side slip angle when driving in a curve. This means that the toe angle of the front wheels and the toe angle change caused by the driving force of the front wheels reduce the side slip angle of the front axle, the toe of the rear wheels and the change of the toe angle caused by the driving force of the rear wheels will increase the side slip angle of the rear axle. As a result, the vehicle tend to oversteer.

Figure 20.

The effect of toe angle on wheel’s side slip angle.

The average side slip stiffness of the inner and outer wheels is denoted as $K_{f}$ , without considering the toe angle of the wheel, the side slip stiffness of the inner wheel of the front wheel is denoted as $K_{yi}$ , and the side slip stiffness of the outer wheel of the front wheel is denoted as $K_{ya}$ . When considering the toe angle of the wheel, the side slip stiffness of the inner wheel of the front wheel is denoted as $K_{yi}^{'}$ , and the side slip stiffness of the outer wheel of the front wheel is denoted as $K_{ya}^{'}$ . Regardless of toe angle, the total side slip stiffness $K_{yi} + K_{ya}$ of the inner and outer wheels is greater than the average side slip stiffness $K_{f}$ , when considering the toe angle of the wheel, the total side slip stiffness $K_{yi}^{'} + K_{ya}^{'}$ of the inner and outer wheels is greater than the side slip stiffness $K_{yi} + K_{ya}$ of the inner and outer wheels without considering the toe angle, that is

\frac{K_{yi} + K_{ya}}{2} < \frac{K_{yi}^{'} + K_{ya}^{'}}{2} < K_{f}

(66)

The toe angle of the wheel itself and the change of the toe angle will affect the tire side slip angle. The change of the wheel’s toe angle caused by the wheel lateral force, longitudinal force and align torque is expressed as $Δ α (F_{y})$ , $Δ α (F_{x})$ , and $Δ α (T)$ , denote the change of toe angle caused by the suspension roll as $Δ α (κ)$ , therefore, the toe angle of the vehicle’s wheel can be expressed as:

{\begin{matrix} α_{v} = \pm α_{V 0} + Δ α_{V} (κ) + Δ α_{V} (F_{yV}) + Δ α_{V} (F_{xV}) + Δ α_{V} (T) \\ α_{H} = \pm α_{H 0} + Δ α_{H} (κ) + Δ α_{H} (F_{yH}) + Δ α_{H} (F_{xH}) + Δ α_{H} (T) \end{matrix}

(67)

$α_{V}$ is the toe angle of the front wheels, $α_{H}$ is the toe angle of the rear wheels, $α_{V 0}$ is the initial toe angle of the front wheels of the vehicle, and $α_{H 0}$ is the initial toe angle of the rear wheels of the vehicle.

For the toe angle, when designing the suspension, it requires to increase the side slip angle of the wheel during the change of the toe angle of the front axle caused by the deformation of the bushing, the toe angle of the rear axle wheels requires to reduce the wheel’s side slip angle during the change of the toe angle of the rear axle, so that the vehicle tends to understeer.

The effect of suspension motion and the change of camber caused by bushing on vehicle’s steering characteristic

The variation characteristics of wheel’s side slip Angle with relative centripetal acceleration are shown in Figure 21. When the camber is positive, the side slip angle is greater than when there is no camber, and when the camber is negative, it is smaller than when there is no camber. The understeer/oversteer trend of the vehicle will be affected by changes in camber as the function of centripetal acceleration (or the function of roll or the function of lateral force). When the centripetal acceleration is low, the front wheel’s camber changes in positive direction/rear wheel no change, or front wheel’s camber no change/rear wheel’s camber change in negative direction, the understeer tendency of the vehicle will increase; otherwise, the understeer will decrease. When the lateral acceleration is high, the front wheels have negative camber/the rear wheels have no camber, there is a slight oversteer of the vehicle; while the front wheels without camber/the rear wheels have a positive camber, there is a strong oversteer characteristic of the vehicle.

Figure 21.

Effect of camber on understeer/oversteer characteristic.

The camber of the wheel will have an effect on the side slip angle of the wheel. The change of camber angle caused by the lateral force, longitudinal force and the align torque of the wheel is expressed as $Δ φ (F_{y})$ , $Δ φ (F_{x})$ , and $Δ φ (T)$ , and the change of camber angle caused by the suspension roll is expressed as $Δ φ (κ)$ , so the camber angle of the vehicle’s wheel can be expressed as:

{\begin{matrix} φ_{v} = \pm φ_{V 0} + Δ φ_{V} (δ_{VL}) + Δ φ_{V} (κ) + Δ φ_{V} (F_{yV}) + Δ φ_{V} (F_{xV}) \\ + Δ φ_{V} (T) \\ φ_{H} = \pm φ_{H 0} + Δ φ_{H} (κ) + Δ φ_{H} (F_{yH}) + Δ φ_{H} (F_{xH}) + Δ φ_{H} (T) \end{matrix}

(68)

$φ_{V}$ is the camber angle of the front wheels, $φ_{H}$ is the camber angle of the rear wheels, $φ_{V 0}$ is the initial camber angle of the front wheels of the vehicle, and $φ_{H 0}$ is the initial camber angle of the rear wheels of the vehicle.

When designing the suspension, for the wheel’s camber angle, when the centripetal acceleration is low, the deformation of the bushing is required to make the front wheel camber change in positive direction/rear wheel no change, or front wheel camber no change/rear wheel’s camber changes in negative direction, which will increase the understeer tendency of the vehicle. When the lateral acceleration is high, let the front wheels have negative camber/the rear wheels have no camber, so that the vehicle has a slight oversteer to avoid strong oversteer characteristic.

Lateral force on vehicle’s tires

During the vehicle’s movement, the lateral force received by a vehicle’s tire can be calculated by two parts, one is the product of the tire’s side slip stiffness and the cornering angle, and the other is the product of the lateral camber thrust coefficient and the camber angle. The side slip angle and camber angle of the tire are obtained by the suspension movement and suspension bushing deformation. The force obtained by the product of the tire side slip stiffness and the tire side slip angle is:

{\begin{matrix} Y_{f} = Y_{f 1} + Y_{f 2} = k_{fL} α_{fL} + k_{fR} α_{fR} \\ Y_{r} = Y_{r 1} + Y_{r 2} = k_{rL} α_{rL} + k_{rR} α_{rR} \end{matrix}

(69)

The lateral force obtained by product of the tire’s lateral thrust coefficient and the camber angle is:

{\begin{matrix} 2 Y_{cf} = - 2 K_{cf} ϕ_{f} \\ 2 Y_{cr} = - 2 K_{cr} ϕ_{r} \end{matrix}

(70)

The total lateral force produced by the suspension movement and the bushing’s deformation is:

\sum F_{y} = 2 Y_{f} + 2 Y_{r} + 2 Y_{cf} + 2 Y_{cr}

(71)

Multi-objective optimization of vehicle’s handling stability

There are multiple evaluation parameters for the transient characteristics of vehicle. In order to explain these evaluation parameters in the form of formulas, the three degree of freedom vehicle dynamics model is used. Since there are multiple optimization targets and multiple optimization variables in the optimization, therefore, optimizing the bushing of the suspension based on the sine swept input to improve the transient dynamics of the vehicle requires multi-objective optimization. Therefore, on the basis of the three degree of freedom vehicle’s dynamics, the optimization targets of vehicle’s transient characteristic are derived, the multi-objective optimization program are established.

Optimize objective function

For parameter-oriented vehicle’s dynamics models, the three degree of freedom vehicle dynamics model is mostly used at present, because the three degrees of freedom vehicle’s dynamic model is simpler than other dynamic models, but it can describe the main motion characteristic and force characteristic of the vehicle. The vehicle’s three degree of freedom model considers the vehicle’s lateral, yaw and roll these three degrees of freedom. For the transient characteristic of the vehicle, the three degree of freedom vehicle dynamics model is also easy to describe the frequency characteristic. The main parameters such as gain and delay time can also be derived with the help of the three degree of freedom vehicle dynamics model. Because the sine swept steering input simulation is used in this article, the vehicle’s roll, lateral and yaw these three degrees of freedom are also mainly considered. The three degree of freedom vehicle dynamics model is:

{\begin{matrix} mV (\frac{d β}{dt} + \overset{\cdot}{ψ}) - m_{s} h_{s} \frac{d^{2} κ}{d t^{2}} = 2 K_{f} (δ + \frac{\partial α_{f}}{\partial κ} κ - β - \frac{a}{V} \overset{\cdot}{ψ}) + 2 K_{r} (\frac{\partial α_{r}}{\partial κ} κ - β + \frac{b}{V} \overset{\cdot}{ψ}) \\ - 2 (K_{cf} \frac{\partial κ_{f}}{\partial κ} ϕ + K_{cr} \frac{\partial κ_{r}}{\partial κ} κ) \\ I \frac{dr}{dt} - I_{xz} \frac{d^{2} κ}{d t^{2}} = 2 K_{f} (δ + \frac{\partial α_{f}}{\partial κ} κ - β - \frac{a}{V} \overset{\cdot}{ψ}) a - 2 K_{r} (\frac{\partial α_{r}}{\partial κ} κ - β + \frac{b}{V} \overset{\cdot}{ψ}) b \\ - 2 (a K_{cf} \frac{\partial κ_{f}}{\partial κ} κ - b K_{cr} \frac{\partial κ_{r}}{\partial κ} κ) \\ I_{κ} \frac{d^{2} κ}{d t^{2}} - I_{xz} \frac{dr}{dt} - m_{s} h_{s} V (\frac{d β}{dt} + \overset{\cdot}{ψ}) = (- K_{κ} + m_{s} g h_{s}) κ - C_{κ} \frac{d κ}{dt} \end{matrix}

(72)

After Laplace transform of equation (72), the characteristic equation is:

[\begin{matrix} mVs + 2 (K_{f} + K_{r}) & mV + \frac{2 (a K_{f} - b K_{r})}{V} & - m_{s} h_{s} s^{2} - 2 Y_{κ} \\ 2 (a K_{f} - b K_{r}) & Is + \frac{2 a^{2} K_{f} + b^{2} K_{r}}{V} & - I_{xz} s^{2} - 2 N_{κ} \\ - m_{s} h_{s} Vs & - I_{xz} s - m_{s} h_{s} V & I_{κ} s^{2} + C_{κ} s + (K_{κ} - m_{s} g h_{s}) \end{matrix}] = 0

(73)

Based on the formula (73), the objective function based on the sine swept steering input simulation is:

\frac{β (s)}{δ (s)} = \frac{| \begin{matrix} 2 K_{f} & mV + \frac{2 (a K_{f} - b K_{r})}{V} & - m_{s} h_{s} s^{2} - 2 Y_{κ} \\ 2 a K_{f} & Is + \frac{2 a^{2} K_{f} + b^{2} K_{r}}{V} & - I_{xz} s^{2} - 2 N_{κ} \\ 0 & - I_{xz} s - m_{s} h_{s} V & I_{κ} s^{2} + C_{κ} s + (K_{κ} - m_{s} g h_{s}) \end{matrix} |}{| \begin{matrix} mVs + 2 (K_{f} + K_{r}) & mV + \frac{2 (a K_{f} - b K_{r})}{V} & - m_{s} h_{s} s^{2} - 2 Y_{κ} \\ 2 (a K_{f} - b K_{r}) & Is + \frac{2 a^{2} K_{f} + b^{2} K_{r}}{V} & - I_{xz} s^{2} - 2 N_{κ} \\ - m_{s} h_{s} Vs & - I_{xz} s - m_{s} h_{s} V & I_{κ} s^{2} + C_{κ} s + (K_{κ} - m_{s} g h_{s}) \end{matrix} |}

(74)

\frac{\overset{\cdot}{ψ} (s)}{δ (s)} = \frac{| \begin{matrix} mVs + 2 (K_{f} + K_{r}) & 2 K_{f} & - m_{s} h_{s} s^{2} - 2 Y_{κ} \\ 2 (a K_{f} - b K_{r}) & 2 a K_{f} & - I_{xz} s^{2} - 2 N_{κ} \\ - m_{s} h_{s} Vs & 0 & I_{κ} s^{2} + C_{κ} s + (K_{κ} - m_{s} g h_{s}) \end{matrix} |}{| \begin{matrix} mVs + 2 (K_{f} + K_{r}) & mV + \frac{2 (a K_{f} - b K_{r})}{V} & - m_{s} h_{s} s^{2} - 2 Y_{κ} \\ 2 (a K_{f} - b K_{r}) & Is + \frac{2 a^{2} K_{f} + b^{2} K_{r}}{V} & - I_{xz} s^{2} - 2 N_{κ} \\ - m_{s} h_{s} Vs & - I_{xz} s - m_{s} h_{s} V & I_{κ} s^{2} + C_{κ} s + (K_{κ} - m_{s} g h_{s}) \end{matrix} |}

(75)

\frac{κ (s)}{δ (s)} = \frac{| \begin{matrix} mVs + 2 (K_{f} + K_{r}) & mV + \frac{2 (l_{f} K_{f} - l_{r} K_{r})}{V} & 2 K_{f} \\ 2 (l_{f} K_{f} - l_{r} K_{r}) & Is + \frac{2 l_{f}^{2} K_{f} + l_{r}^{2} K_{r}}{V} & 2 l_{f} K_{f} \\ - m_{s} h_{s} Vs & - I_{xz} s - m_{s} h_{s} V & 0 \end{matrix} |}{| \begin{matrix} mVs + 2 (K_{f} + K_{r}) & mV + \frac{2 (l_{f} K_{f} - l_{r} K_{r})}{V} & - m_{s} h_{s} s^{2} - 2 Y_{κ} \\ 2 (l_{f} K_{f} - l_{r} K_{r}) & Is + \frac{2 l_{f}^{2} K_{f} + l_{r}^{2} K_{r}}{V} & - I_{xz} s^{2} - 2 N_{κ} \\ - m_{s} h_{s} Vs & - I_{xz} s - m_{s} h_{s} V & I_{κ} s^{2} + C_{κ} s + (K_{κ} - m_{s} g h_{s}) \end{matrix} |}

(76)

In equations (74)–(76), $Y_{κ}$ and $N_{κ}$ are:

{\begin{matrix} Y_{κ} = (\frac{\partial α_{f}}{\partial κ} K_{f} + \frac{\partial α_{r}}{\partial κ} K_{r}) - (\frac{\partial κ_{f}}{\partial κ} K_{cf} + \frac{\partial κ_{r}}{\partial κ} K_{cr}) \\ N_{κ} = (\frac{\partial α_{f}}{\partial κ} a K_{f} - \frac{\partial α_{r}}{\partial κ} b K_{r}) - (\frac{\partial κ_{f}}{\partial κ} a K_{cf} - \frac{\partial κ_{r}}{\partial κ} b K_{cr}) \end{matrix}

(77)

In equation (72), the side slip stiffness $K_{f}$ and $K_{r}$ of the vehicle’s tires do not consider the influence of suspension movement and bushing deformation. The front wheels’ side slip angle $(δ + \frac{\partial α_{f}}{\partial κ} κ - β - \frac{a}{V} \overset{\cdot}{ψ})$ and camber angle $\frac{\partial κ_{f}}{\partial κ} κ$ of the vehicle, as well as the rear wheels’ side slip angle $(\frac{\partial α_{r}}{\partial κ} κ - β + \frac{b}{V} \overset{\cdot}{ψ})$ and camber angle $\frac{\partial κ_{r}}{\partial κ} κ$ , only consider the additional corners generated by the movement of the suspension, and do not take into account the impact of the deformation of the suspension bushing on the wheel side slip angle and camber angle. When considering the deformation of the suspension bushing, the deformation of the bushing will cause the compliance characteristic of the suspension to change. It can be seen from equation (69) that the sum of the side slip stiffness of the left and right tires will change accordingly, and the side slip stiffness of the front axle is no longer $K_{f}$ , the side slip stiffness of the rear axle is no longer $K_{r}$ . From equations (70) and (71), it can be seen that the side slip angle and camber angle will also change further, which will cause the transient characteristics of the vehicle to change.

In equations (74)–(76), make $s = j ω$ , where $j = \sqrt{- 1}$ , the frequency response characteristic function of the vehicle can be obtained, including the amplitude-frequency characteristic and the phase-frequency characteristic. Among them, $| \frac{\overset{\cdot}{ψ}}{δ} (ω) |$ , $| \frac{κ}{δ} (ω) |$ , and $| \frac{β}{δ} (ω) |$ are the gains of the yaw rate, roll angle and the side slip angle relative to the steering wheel angle respectively. The frequency corresponding to the maximum yaw rate gain is the resonance frequency of the vehicle. The resonance frequency of the vehicle is set as one of the objective functions. $θ_{ψ} (ω)$ , $θ_{κ} (ω)$ , and $θ_{β} (ω)$ are the corresponding phase lag angles. Since the general driver turns the steering wheel at a frequency of about 0.5 Hz, the gain of the yaw rate, roll angle, and side slip angle at 0.5 Hz relative to the steering wheel angle and the corresponding delay time at 0.5 Hz are taken as optimize targets. The phase lag angle can be converted to delay time by the following formula:

t = \frac{θ + 180}{360 • f} or \frac{θ - 180}{360 • f}

(78)

The overall optimization objective functions are:

{\begin{matrix} f_{1} = {| \frac{\overset{\cdot}{ψ}}{δ} (ω) |}_{ω = 0.5 Hz} \\ f_{2} = {| \frac{κ}{δ} (ω) |}_{ω = 0.5 Hz} \\ f_{3} = {| \frac{β}{δ} (ω) |}_{ω = 0.5 Hz} \\ f_{4} = ω {{| \frac{\overset{\cdot}{ψ}}{δ} (ω) |}_{max}} \\ f_{5} = \frac{θ_{ψ} + 180}{360 • f} or \frac{θ_{ψ} - 180}{360 • f} \\ f_{6} = \frac{θ_{κ} + 180}{360 • f} or \frac{θ_{κ} - 180}{360 • f} \end{matrix}

(79)

Sensitivity analysis of the influence of suspension bushing’s stiffness on vehicle handling stability

The vehicle’s parameters used in the optimization are shown in Table 1:

Table 1.

Vehicle parameters.

Parameter	Value	Parameter	Value
Sprung mass (kg)	1067	Height from center of mass to ground (mm)	330
Unsprung mass (kg)	595	Distance between rear wheels (mm)	1487.8
Center of mass to front wheelbase (mm)	1452.6	Distance between front wheels (mm)	1508.3
Center of mass to rear wheelbase (mm)	1249.5	Wheelbase (mm)	2702.1

There are eight bushings connecting the connection rods and the vehicle body in Figures 10 and 13, four for the front suspension and four for the rear suspension, and each bushing has stiffness in six directions. Taking the upper rear control arm on the front suspension as an example, the bushing is named $K_{2}$ , then the displacement stiffness of the bushing in the three directions $x$ , $y$ , $z$ , and the rotational stiffness around the $x$ , $y$ , and $z$ axis are denoted as $K_{2}^{x}$ , $K_{2}^{y}$ , $K_{2}^{z}$ , $K_{2}^{tx}$ , $K_{2}^{ty}$ , and $K_{2}^{tz}$ .

There are eight bushings in Figures 10 and 12, so there are 48 stiffness values, but some stiffness values have little effect on the handling stability of the vehicle, and some stiffness values have a greater influence on the handling stability of the vehicle. In order to find out the bushing stiffness values that have a greater impact on the handling stability of the vehicle, the sensitivity analysis of these 48 stiffness values is carried out first, and then several stiffness values with higher sensitivity are selected for optimization analysis. After the sensitivity analysis, five stiffness values with greater sensitivity are selected as shown in Table 2.

Table 2.

Bushing stiffness with larger sensitivity value.

Stiffness	Resonance frequency ( $f_{4}$ )	Yaw rate gain ( $f_{1}$ )	Max yaw rate gain ( $f_{3}$ )	Roll angle gain ( $f_{2}$ )	Delay time of lateral acceleration ( $f_{6}$ )	Delay time of yaw rate ( $f_{5}$ )
$K_{A 1}^{z}$	0.97	1.69	−0.19	1.62	2.08	2.16
$K_{A 3}^{tz}$	2.22	1.99	−1.13	1.77	1.24	1.60
$K_{A 4}^{y}$	1.86	2.38	−0.44	2.24	0.83	1.45
$K_{A 7}^{z}$	−2.39	−1.96	1.43	−2.53	−1.50	2.00
$K_{A 8}^{tz}$	1.63	1.69	−1.18	1.96	2.05	2.28

Optimization algorithm

By optimizing the stiffness of the suspension bushing to improve the handling stability of the vehicle, the multi-objective optimization method and algorithm are adopted. The NSGA-II algorithm is a relatively successful algorithm used in multi-objective optimization algorithms.¹⁶ In non-dominated sorting, individuals close to the Pareto frontier are selected, which enhances Pareto’s forward ability. Combined with the analysis of section 7.2, the mathematical expression of multi-objective optimization is as follows:

\begin{matrix} F (X) = {min F_{1} {f_{1} (X), f_{2} (X), f_{3} (X), f_{5} (X), f_{6} (X)} \\ + max F_{2} {f_{4}}} \end{matrix}

(80)

{\begin{matrix} X = {x_{1}, x_{2}, \dots, x_{5}} \\ 0 \leq x_{i} \leq 2.5, i = 1, 2, \dots, 5 \end{matrix}

(81)

$F (X)$ is the objective function, $f_{k} (X)$ is the sub-objective function, $X$ is the optimized variable. In the NSGA-II algorithm, the calculation mechanism of crossover and sudden mutation is called SBX method. The crossover operation in generating sub-individuals according to the SBX method is:

{\begin{matrix} x_{i}^{(1, t + 1)} = \frac{1 + β_{qi}}{2} x_{i}^{(1, t)} + \frac{1 + β_{qi}}{2} x_{i}^{(2, t)} \\ x_{i}^{(2, t + 1)} = \frac{1 + β_{qi}}{2} x_{i}^{(1, t)} + \frac{1 + β_{qi}}{2} x_{i}^{(2, t)} \end{matrix}

(82)

The mutation operation in generating sub-individuals according to the SBX method is:

x_{i}^{(1, t = 1)} = x_{i}^{(1, t)} + δ_{q} (x_{i}^{UB} - x_{i}^{LB})

(83)

δ_{q} = {\begin{matrix} {[2 u + (1 - 2 u) {(1 - δ)}^{η_{m} + 1}]}_{η_{m} + 1}^{1} - 1, u \leq 0.5 \\ 1 - {[2 (1 - u) + 2 (u - 0.5) {(1 - δ)}^{η_{m} + 1}]}_{η_{m} + 1}^{1}, u > 0.5 \end{matrix}

(84)

δ = min (x_{i} - x_{i}^{LB}, x_{i}^{UB} - x_{i}) / (x_{i}^{UB} - x_{i}^{LB}), u \in [1, 0]

(85)

At present, the optimization objective function used in the optimization of vehicle handling stability is mainly realized by vehicle dynamic parameter modeling. In the process of establishing parameter model, some secondary factors and degrees of freedom of the vehicle are often ignored, which leads to the optimized model compared with the actual vehicle, the error is larger. In order to make the vehicle’s model as close as possible to the real vehicle model, this paper adopts the multi-body dynamics model, that is, the multi-objective optimization of vehicle’s handling stability is carried out for the multi-body dynamics model.

Since the sine swept steering input test is used in this paper to evaluate the handling stability of a vehicle, the parameters of the sine swept steering input simulation to evaluate the vehicle’s handling stability are taken as the optimization objectives, and the suspension bushing stiffness with higher sensitivity are taken as the optimization variables to perform multi objective optimization. After the multi objective optimization of the whole vehicle, the optimization iterative diagram of the optimization targets is shown in Figures 22 to 27.

Figure 22.

Iterative diagram of the maximum gain value of the yaw rate relative to the steering wheel angle.

Figure 23.

Iterative graph of yaw rate gain at 0.5 Hz relative to the steering wheel angle.

Figure 24.

Iterative diagram of resonance frequency.

Figure 25.

Iteration diagram of roll angle gain relative to the lateral acceleration.

Figure 26.

Iterative diagram of lateral acceleration delay time.

Figure 27.

Iterative diagram of yaw rate delay time.

Figures 22 and 23 are iterative graphs of yaw rate gain with the steering wheel angle. Figure 22 is the maximum value of the yaw rate gain relative to the steering wheel angle, and Figure 23 is the yaw rate gain value at 0.5 Hz relative to the steering wheel angle. In the performance of the frequency characteristics of the vehicle, the smaller the maximum gain value of the yaw rate relative to the steering wheel angle, the more stable the vehicle. Therefore, it is hoped that the maximum gain value of the yaw rate relative to the steering wheel angle will be minimized during optimization. The change trend of the iterative graphs in Figures 22 and 23 meets the requirements of the optimization goal.

Figure 24 is an iterative diagram of the resonance frequency. The larger the resonance frequency, the more stable the vehicle. The resonance frequency in Figure 24 is continuously increasing as the optimization progresses, which meets the requirements of the optimization goal. Figure 25 shows the roll angle gain relative to the lateral acceleration. In Figure 25, the roll angle gain relative to the lateral acceleration is increased compared to the original vehicle, but the increased value is small, so it has little effect on the handling stability of the vehicle.

Figures 26 and 27 show the delay time of lateral acceleration relative to the steering wheel angle and the delay time of yaw rate relative to lateral acceleration. The smaller the delay time, the faster the response of the vehicle. Figures 26 and 27 reflect that the delay time is decreasing, which is beneficial to improve the response speed of the vehicle, indicating that the optimization is carried out in accordance with the required target.

Optimization results

After the multi-objective optimization for the transient characteristic of the vehicle, the corresponding Pareto solution sets are obtained, as shown in Figures 28 to 33. Among them, the black point represents the general solution, the blue point represents the Pareto solution, and the green point represents the relatively optimal solution. Delay time 1 refers to the delay time of the lateral acceleration relatives to the steering wheel angle, and delay time 2 refers to the delay time of yaw rate relatives to the lateral acceleration.

Figure 28.

Pareto solution sets of yaw rate gain at 0.5 Hz with resonance frequency.

Figure 29.

Pareto solution sets of roll angle gain at 0.5 Hz relative to yaw rate gain at 0.5 Hz.

Figure 30.

Pareto solution sets of delay time 1 with yaw rate gain at 0.5 Hz.

Figure 31.

Pareto solution sets of delay time 1 relative to resonance frequency.

Figure 32.

Pareto solution sets of delay time 1 relative to delay time 2.

Figure 33.

Pareto solution sets of delay time 2 relative to yaw rate gain.

Through optimization analysis, Pareto solution sets and relative optimal solution are finally obtained. The comparison of parameters before and after optimization of the vehicle are shown in Table 3.

Table 3.

Comparison of parameters before and after optimization.

	$f_{1}$ (s⁻¹)	$f_{2}$ (s⁻¹)	$f_{3}$ (s⁻¹)	$f_{4}$ (Hz)	$f_{5}$ (s)	$f_{6}$ (s)	$x_{1}$	$x_{2}$	$x_{3}$	$x_{4}$	$x_{5}$
Original vehicle	0.47	0.46	4.40	0.63	0.20	0.14	1	1	1	1	1
Optimized vehicle	0.44	0.42	4.52	0.83	0.17	0.14	1.31	0.001	0.03	0.18	0.21
Percentage change (%)	−6.4	−8.7	2.7	31.7	−15	0	+31	−99.9	−97	−82	−79

Comparative analysis of optimization results of vehicle handling and stability

In order to analysis the optimization effect, a Pareto solution is selected from the optimization results and compared with original vehicle. In order to ensure that the handling stability of the vehicle is optimized and improved, comparisons are made from two aspects of the vehicle’s transient characteristics and steady-state characteristics.

Comparative analysis of vehicle’s transient characteristics

Bring the selected Pareto solution into the original vehicle for simulation analysis, and compare the results with the original vehicle, the result comparison chart in Figures 34 to 39 are obtained. Figure 34 shows the yaw rate gain relative to the steering wheel angle. In Figure 34, in the range of 0 to 1 Hz, especially at 0.5 Hz, the curve of the optimized yaw rate gain relative to the steering wheel angle is smaller than that of the original vehicle, which meets the requirements of the optimization goal. When the frequency exceeds 1 Hz, the difference of the yaw rate gain before and after optimization is not obvious. Since the frequency of the driver in the process of driving the vehicle rarely exceeds the range of 1 Hz, the high frequency range is not considered.

Figure 34.

Yaw rate gain relative to the steering wheel angle.

Figure 35.

Delay time of yaw rate relatives to the lateral acceleration.

Figure 36.

Roll angle gain relative to the lateral acceleration.

Figure 37.

An enlarged view of roll angle gain relative to the lateral acceleration.

Figure 38.

Delay time of lateral acceleration relative to the steering wheel angle.

Figure 39.

The enlarged view of the lateral acceleration’s delay time relative to the steering wheel angle.

Figure 35 is the delay time of yaw rate relatives to the lateral acceleration. In Figure 35, at 0.5 Hz, the difference to the delay time of yaw rate between optimized and original vehicle is obvious, and it is reduced after optimization compared with the original vehicle, which meets the requirements of optimization goals. Between 0 and 1 Hz, the optimized delay time of yaw rate is basically smaller than the original vehicle. When the frequency exceeds 1 Hz, the difference of yaw rate delay time between the optimized vehicle and the original vehicle is small. Since the frequency rarely reaches above 1 Hz during the operation of the vehicle, the situation above 1 Hz will not be discussed too much.

Figure 36 shows the roll angle gain relative to the lateral acceleration. In Figure 36, from low frequency to high frequency, as the frequency increases, the difference of roll angle gain between the optimized vehicle and original vehicle is relatively small. After enlarging the curve at 0.5 Hz, the situation shown in Figure 37 is obtained. In Figure 37, the roll angle gain relative to the lateral acceleration after optimization is increased than original vehicle, but the increased value is relatively small and within an acceptable range. Because it is a multi-objective optimization, it is normal for individual optimization objectives to fail to reach the optimum.

Figure 38 shows the delay time of the lateral acceleration relative to the steering wheel angle, and after enlarging the range around 0.5 Hz in Figure 38, the delay time of the lateral acceleration around 0.5 Hz relative to the steering wheel angle as shown in Figure 39 is obtained.

In Figures 38 and 39, the absolute value of the lateral acceleration’s delay time of the optimized vehicle is reduced, but in general, the difference of the lateral acceleration’s delay time between the optimized vehicle and original vehicle in the entire frequency range is smaller. Because it is a multi-objective optimization, it is difficult to achieve all the optimization objectives to reach the optimal state, the lateral acceleration’s delay time after optimizing didn’t significantly reduced but it is acceptable.

Comparison and analysis of steady state with constant radius cornering simulation

In order to further verify the improvement of the vehicle’s dynamics characteristics after the vehicle optimization, the constant radius cornering simulation of the whole vehicle before and after the optimization is compared and analyzed. The curve of the steering wheel angle with the lateral acceleration shown in Figure 40 and the curve of the roll angle with the lateral acceleration shown in Figure 41 are obtained.

Figure 40.

The curve of steering wheel angle with the lateral acceleration.

Figure 41.

The curve of roll angle with the lateral acceleration.

In Figure 40, at the same lateral acceleration, the steering wheel angle of the optimized vehicle is larger, indicating that the understeer characteristic of the vehicle is increased, and the vehicle tends to be more stable. In Figure 41, under the same lateral acceleration, the optimized vehicle’s roll angle has increased, but the increase value is not much and within in the acceptable range.

Conclusion

The stiffness calculation of vehicle suspension bushing is analyzed and the multi-body dynamic model of vehicle’s front and rear suspension is established according to the actual vehicle. The influence mechanism of the bushing on the suspension wheel alignment angle and the suspension force is analyzed in detail and systematically. Through sensitivity analysis, five bushes that have a greater impact on the transient characteristics of the vehicle are selected, and multi-objective optimization is performed for the selected bushings. The specific conclusions obtained are:

Through the calculation of the vector, the link position and link angle during the movement of the double wishbone suspension are analyzed in detail, and the wheel alignment angle of the double wishbone suspension is finally calculated. Through the simulation calculation of the multi-body dynamics model, the toe angle and camber angle change with the wheel travel are obtained.

The bushing deformation and bushing force of double wishbone suspension and multi-link suspension are analyzed, and the relationship between bushing deformation and suspension compliance characteristics is explained.

Sensitivity analysis shows that for the double wishbone suspension, the torsional stiffness of the lower front control arm bushing around the z axis, the translational stiffness of the lower rear control arm bushing along the y direction, and the translational stiffness of upper front control arm bushing along the z direction have a greater impact on the vehicle’s transient characteristic. For the rear multi-link suspension, the stiffness of the lower front control arm bushing along the z direction and the torsional stiffness of the lower rear control arm bushing around the z-axis have a greater impact on the vehicle’s transient characteristics.

For the front suspension, when the torsional stiffness of the lower front control arm bushing around the z axis increases by 31%, the translational stiffness of the lower rear control arm bushing along the y direction is reduced by 99.9%, and the translational stiffness of the upper front control arm bushing in the z direction is reduced by 97%. For the rear suspension, when the stiffness of the lower front control arm bushing along the z direction is reduced by 82% and the torsional stiffness of the lower rear control arm bushing around the z axis is reduced by 79%, the yaw rate gain at 0.5 Hz is reduced by 6.4%, and the roll angle gain is increased by 2.7%, the delay time of yaw rate is reduced by 15%, the delay time of roll angle is unchanged, and the resonance frequency is increased by 31.7%. From the perspective of steady-state conditions, the understeer trend of the vehicle increases, and the roll angle increases slightly, but it is within an acceptable range.

Footnotes

Appendix

Handling Editor: James Baldwin

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This project was supported by National Natural Science Foundation of China (NSFC) (No.51965026), Yunnan Province Applied Basic Research Foundation (No. 2018FB097) and Scientific Research Fund Project of Yunnan Provincial Education Department (No. 2018JS022). The authors are greatly appreciated for the financial support.

ORCID iD

Jin Gao

References

García Tárrago

Kari

Vinolas

, et al. Frequency and amplitude dependence of the axial and radial stiffness of vehiclebon-black filled rubber bushings. Polymer Test 2007; 26: 629–638.

Gibanica

Abrahamsson

TJS

Olsson

. Model updating of multiple nominally identical vehicle components. SAE Technical Paper 2018-01-0554, 2018.

Plagge

Ricker

Kröger

, et al. Efficient modeling of filled rubber assuming stress-induced microscopic restructurization. Int J Eng Sci 2020; 151: 1021–1028.

Sedlaczek

Dronka

Rauh

. Advanced modular modelling of rubber bushings for vehicle simulations. Veh Syst Dyn 2011; 49: 741–759.

Aydemir

Sendur

. Simplified transfer function approach for modeling frequency dependency of damping characteristics of rubber bushings. Proc Inst Mech Eng D J Automobile Eng 2019; 233: 2518–2531.

Zenglian

Lljun

Zhuoping

. A novel empirical model of rubber bushing in automotive suspension system. Int Conf Motion Vib Control 2010; 2010: _3A12 to 1_.

Sancibrian

Garcia

Viadero

, et al. Kinematic design of double-wishbone suspension systems using a multiobjective optimisation approach. Veh Syst Dyn 2010; 48: 793–813.

Balike

Rakheja

Stiharu

. Development of kineto-dynamic quarter-vehicle model for synthesis of a double wishbone suspension. Veh Syst Dyn 2011; 49: 107–128.

Attia

. Dynamic modelling of the double wishbone motor-vehicle suspension system. Eur J Mech 2002; 21: 167–174.

10.

Tanik

Parlaktaş

. On the analysis of double wishbone suspension. J Adv Mech Des 2016; 10: 1–11.

11.

Shen

Lee

W-T

Russell

. A small-scale optimization model for RSSR-SS motion generation with branch and order defect elimination in Matlab. J Adv Mech Des Syst Manuf 2014; 8: 1–8.

12.

Reddy

Kodati

Chatra

, et al. A comprehensive kinematic analysis of the double wishbone and MacPherson strut suspension systems. Mech Mach Theory 2016; 105: 441–470.

13.

Tang

Shangguan

W-B

Dai

. A calculation method of joint forces for a suspension considering nonlinear elasticity of bushings. Proc Inst Mech Eng K J Multi-body Dyn 2012; 226: 281–297.

14.

Nabawy

Abdelrahman

Abdalla

, et al. Analysis of the dynamic behavior of the double wishbone suspension system. Int J Appl Mech 2019; 11: 1–19.

15.

Wang

X-L

Dai

Shangguan

W-B

. Calculation of joint forces of a multi-link suspension for strength and fatigue analysis of bushings and control arms. Int J Veh Des 2014; 66: 217–234.

16.

Zhao

Liu

Zhang

, et al. An improved adaptive NSGA-II with multi-population algorithm. Appl Intell 2019; 49: 569–5801.

Analysis and optimization of the vehicle handling stability with considering suspension kinematics and compliance characteristics

Abstract

Keywords

Introduction

Establish vehicle suspension model and calculation of rubber bushing’s stiffness

Kinematics analysis of vehicle’s suspension

Analysis and calculation of vector e → k

Analysis and calculation of vector e → hub

Calculation of wheel alignment angle

Statics analysis of double wishbone suspension

Displacement analysis of double wishbone suspension’s connecting point

Calculate the forces and moments acting on the bushing

Statics analysis of multi-link suspension

Displacement analysis of multi-link suspension’s connection point

Calculate the forces and moments acting on the rubber bushing

Force balance and moment balance analysis

Analysis of the influence of suspension kinematics and compliance characteristic on vehicle dynamics

The influence of suspension motion and bushing on the change of toe Angle and camber Angle

Effect of toe angle change caused by suspension movement and bushing deformation on tire side slip stiffness

The effect of suspension motion and the change of camber caused by bushing on vehicle’s steering characteristic

Lateral force on vehicle’s tires

Multi-objective optimization of vehicle’s handling stability

Optimize objective function

Sensitivity analysis of the influence of suspension bushing’s stiffness on vehicle handling stability

Optimization algorithm

Optimization results

Comparative analysis of optimization results of vehicle handling and stability

Comparative analysis of vehicle’s transient characteristics

Comparison and analysis of steady state with constant radius cornering simulation

Conclusion

Footnotes

Appendix

Declaration of conflicting interests

Funding

ORCID iD

References

Analysis and calculation of vector ${\vec{e}}_{k}$

Analysis and calculation of vector ${\vec{e}}_{hub}$