Collaborative optimization of bushing installation angle and bushing stiffness of torsion beam suspension to improve the handing stability of the whole vehicle

Abstract

Both the bushing installation angle and bushing stiffness of torsion beam suspension have important and complex effects on suspension characteristics and vehicle handling stability. For improving the vehicle handling stability by optimizing the two, this paper systematically analyzes the influence mechanism of the bushing installation angle and bushing stiffness of the torsion beam suspension on the suspension C characteristic. The nonlinear relationship between the bushing installation angle and the steady-state characteristics index in the local design space is discussed. The influence mechanism of the bushing installation angle and bushing stiffness on the frequency response characteristics is deeply analyzed, the variation relationship of the frequency characteristic index at the low frequency of 0.5 Hz under different bushing installation angles and bushing stiffnesses is obtained. Finally, the bushing stiffness, which has a great influence on the frequency characteristics, and the bushing installation angle are taken as the optimization variables, and the frequency characteristic index is taken as the optimization goal, and the multi-objective collaborative optimization of the frequency characteristic is carried out with the help of genetic algorithm. From the comparative analysis of transient and steady-state characteristics before and after optimization, it can be seen the vehicle handling stability has been improved.

Introduction

In recent years, a kind of suspension between independent suspension and semi-independent suspension, namely torsion beam suspension, has the characteristics of simple structure, low manufacturing cost, and installation on the vehicle body. And then it is widely used in the rear suspension system of miniature and medium-sized cars. Although the torsion beam suspension structure is simple, many factors such as opening direction, bushing installation angle and bushing stiffness, cross-beam section shape, etc. will affect the suspension characteristics and vehicle handling stability.^1,2 In the design stage, if the above design parameters are different, the performance of the suspension and vehicle performance will also be different. At present, there are few systematic studies on the effects of torsion beam suspension bushing installation angle and bushing stiffness on suspension and vehicle performance. In addition, the bushing is an important bearing component of the suspension. The research on the bushing installation angle and bushing stiffness belongs to an important branch of the automotive field and is of great significance.^3,4

Until now, numerous scholars have executed lots of research on the rubber bushing of various suspension. Aydemir et al. ⁵ put forward a simplified frequency-dependent transfer function model to overcome the problem of lower accuracy due to limitations in representing the frequency-dependent dynamic characteristics of vibration isolators. The proposed method is validated on centain heavy commercial truck, and the proposed model was more accurate than that modeled by Voigt. Khaknejad et al. ⁶ optimized the installation angle and stiffness of the rubber bushing on the vehicle, improved the toe-in stiffness, and improved the handling stability of the vehicle through an example. In order to overcome the trade-off between handling stability and ride comfort in passenger cars equipped with traditional suspension, Qi et al. ⁷ created a new suspension configuration and compared it with the conventional suspension. The test results showed that vehicles equipped with hydraulic interconnected suspension and electronically controlled air springs can achieve higher performance in the light of handling stability and ride better comfort. Dzierzek ⁸ put forward an experimental model, which can simulate the damping and stiffness characteristics of the bushing well, and can be easily applied to the ADAMS model, but this model was limited to describing the axial and radial characteristics of the bushing. However, the torsion and swing characteristics of the bushing were invalid, and the coupling effect of the force in all directions of the rubber bushing was not considered. OK et al. ⁹ analyzed the bushing characteristics under various excitation conditions, and established a new nonlinear bushing characteristic coupling model using the Bouc-Wen hysteresis model. The results showed that the bushing force is well matched with the test, and the established model can be used in vehicle dynamics simulation. Fichera et al. ¹⁰ calculated and analyzed the compliance characteristics of the suspension including the mechanical properties of the bushing, and predicted and evaluated the effect of the rubber bushing on the NVH performance of the vehicle. The simulation results of the linear flexible body multi-body dynamic model and the nonlinear model are analyzed and compared, and the effects of different bushing stiffness characteristics on the performance of the suspension are explored. By view of on the quasi-static test results, Gracia et al. ¹¹ used an optimization algorithm to fit the model parameters required for the finite element calculation, and used the ABAQUS software to embed the Yeoh hyperelastic model and the multi-segment linear stress-strain elastic-plastic model into the cylindrical rubber lining set of finite element models, and the curve of reaction force and displacement of the action point of the rubber part was calculated. Wollscheid et al. ¹² took carbon black filled rubber vibration isolators as the research object, built a model in the pre-deformed elastic rubber material based on the storage modulus and loss modulus, and obtained the preload dependence and frequency dependence through static and dynamic of the test. The constitutive parameters required by the finite element calculation were identified by the optimization of the least squares method. The numerical simulation results show that it can describe the pre-deformation and frequency dependence of the vibration isolator well. Khajehsaeid et al. ¹³ established a hyperelastic-viscoelastic model of rubber bushing by view of on the VUMAT subroutine method of ABAQUS software, and calculated three working conditions of axial, torsional, and axial torsion. The calculation results showed the hyperelastic-viscoelastic model can not only accurately calculate the unidirectional condition, but also can better predict the coupled condition of the torsional axial direction. Liu et al. ¹⁴ established a dynamic mechanical model of rubber bushing by view of the principle of mesh superposition. The model includes a hyperelastic model, a viscoelastic model and a friction model. And through the user-defined program interface of SIMPACK software, the nonlinear dynamic characteristics such as superelasticity, viscoelasticity, and friction characteristics of the bushing were considered in detail.

The above scholars have executed research on the rubber bushing from various thinking angles, yet have not systematically explored the influence mechanism of the bushing installation angle and bushing stiffness of the torsion beam suspension on the C characteristics; The influence mechanism of the bushing installation angle and bushing stiffness of the torsion beam suspension on the frequency characteristics has not been deeply analyzed. Regarding the collaborative optimization of the bushing installation angle and bushing stiffness and the frequency response characteristics of the vehicle, there is no research in this area yet.

In view of the issues that have not yet been discussed, this paper explores the following four aspects: (1) The influence mechanism of the bushing installation angle and bushing stiffness of the torsion beam suspension on the C characteristics is analyzed in detail; (2) The influence mechanism of the bushing installation angle and bushing stiffness of the torsion beam suspension on the frequency response characteristics is revealed in depth; (3) The evaluation indexes of steady-state and transient handling stability are determined, the nonlinear relationship between the bushing installation angle within a certain range and the steady-state handling stability index is obtained. Meanwhile, The variation law of the frequency characteristic index under different bushing installation angles and bushing stiffness when the low frequency is around 0.5 Hz; (4) Taking the frequency characteristic index at 0.5 Hz as the optimization target, and the bushing stiffness, which has a great influence on the frequency characteristics, and the bushing installation angle are taken as the optimization variables, the multi-objective collaborative optimization is carried out, so as to improve the handing stability of vehicle.

Establishment of torsion beam suspension model

During the movement of the torsion beam suspension, the torsion beam will deform, so a flexible body is needed to establish the model. ¹⁰ First, the neutral file MNF of the flexible body of the torsion beam needs to be generated by the finite element software, and then the central file MNF is imported into the ADAMS multi-body dynamics software, so as to establish a rigid-flexible coupled multi-body dynamic model of the torsion beam suspension in Figure 1.

Figure 1.

Multibody dynamic model of torsion beam rear suspension.

The rubber bushing commonly used in suspension is composed of an intermediate rubber body and inner and outer metal sleeves. The rubber bushing is used as the connecting part of the two objects. By view of the relative displacement of the two objects at the connection point, the bushing model applies elastic force and damping force in six directions to the two objects. Therefore, the mechanical properties modeling of the bushing generally considers the stiffness properties in six directions, and the force during the working process is more complicated, which transmits three axial forces and three torsional moments between the body and the guide mechanism.

As shown in Figure 2, the modeling coordinate system takes the axial direction as the z-axis, and the two radial directions as the x-axis and the y-axis. Under the action of the six component forces, deformation in six directions will occur, and correspondingly there are stiffness in six directions. In Figure 2, $F_x,F_y$ represents the radial stiffness of the bushing, $F_z$ represents the axial stiffness of the bushing, $T_x,T_y$ represents the torsional stiffness of the bushing in the radial direction, $T_z$ represents the torsional stiffness in the axial direction.

Figure 2.

Schematic diagram of the force of the rubber bushing.

Figure 3 is schematic diagram of the installation angle of the rubber bushing. The bushing installation angle refers to the acute angle formed by the axial direction of the bushing and the axis of the torsion beam, and the bushing installation angle in the initial state of the suspension modeling is −24.5°. In Figure 3, $F_{x(y)(z)}$ represents the radial or axial stiffness of the bushing, and $T_{x(y)(z)}$ represents the torsional stiffness of the bushing around the radial or axial direction. $x(y)(z)$ represents the directions around the coordinate axes X, Y, and Z respectively. Usually, in order to effectively transmit the lateral force, the rubber bushing must be arranged diagonally. The purpose of changing the bushing installation angle can be achieved by adjusting the hard point coordinates of the bushing installation axis.

Figure 3.

Schematic diagram of bushing installation angle.

Calculation of roll stiffness of torsion beam suspension

The roll stiffness of the torsion beam suspension is mainly constitute of the roll stiffness contributed by the spring and the roll stiffness contributed by the torsion beam. So order to calculate the equivalent roll stiffness contributed by the torsion beam, it is necessary to first derive the torsion stiffness of the torsion beam.

Derivation of torsion stiffness of beam

Usually the stiffness of the trailing arm of the torsion beam suspension is much larger than that of the beam, so this paper mainly considers the torsion stiffness of the beam and the radial stiffness and axial stiffness of the bushing, ignoring the torsion stiffness of the trailing arm, so as to reflect the suspension deformation under various forces.

A pair of equal and opposite moments (torques) are applied to both ends of the beam member, and the torsion stiffness of the beam is defined as:

$$K_{\varphi \mathrm{rT}}=\frac{T}{\eta }$$

(1)

Where, $\eta$ is the torsion angle of the beam.

The V-section beam is a thin-walled member with a non-circular section opening, and the left and right ends are welded on the longitudinal arm. Therefore, the warping freedom of the beam is constrained when the beam is twisted, and the torsion stiffness of the beam is only calculated by the theory of restrained torsion. Before using the constrained torsion initial parameter solution method, the constrained torsion coefficient $k$ must be solved.

$$k=l_0\sqrt{\frac{G{J}_d}{E_1{J}_w}}$$

(2)

Among them, the converted elastic modulus: $E_1=\frac{E}{1-u^2}$ shear modulus: $G=\frac{E}{2(1+u)}$ , centerline length of section: $s=2({\vartheta }_1{r}_3+L_0)$ , free torsional moment of inertia of section: $J_d=\frac{1}{3}\sum _i{d}_i{L}_1^3=\frac{1}{3}s{L}_1^{3}s{L}_1^3°$

In equation (2): $l_0$ is the length of the beam, $r_3$ is the radius of the arc part of the midline of the section, $L_0$ is the length of the rectangular part, $L_1$ is the thickness of the beam, ${\vartheta }_1$ is the radian corresponding to the arc, $G$ is the shear modulus of the material, $E$ is the elastic modulus of the material, $u$ is the Poisson’s ratio, $E_1$ is the converted elastic modulus, $J_w$ is the main sector moment of inertia of the section; $J_d$ is the free torsional moment of inertia of the section.

Therefore, to obtain the constraint torsion coefficient $k$ , it is also necessary to solve the moment of inertia $J_w$ of the main sector of the section.

Since the beam is a symmetric section, the section intersects with the symmetry axis $z$ at point $M$ , and the circle centers $O$ and $M(z,s)$ are assumed to be the poles of the reference sector and the zero of the reference sector, respectively, as shown in Figure 4. Rotating point $M_0$ clockwise along the midline of the section by $\alpha$ angle to reach point $M$ , then the reference sector area $w_1$ of this point is:

$$w_1=\alpha r_3^2,0\le \alpha \le {\vartheta }_1$$

(3)

Figure 4.

Schematic diagram of the solution of the pole of the main sector.

The distance $y_1$ from $M$ point to the z-axis is:

$$y_1=r_3\sin \alpha$$

(4)

The main sectoral moment of inertia $J_{\mathrm{zw}1}$ of the arc part is:

$$J_{\mathrm{zw}1}=\underset{A}{\int }{y}_1{w}_1\mathrm{dA}=L_1{r}_3^4\underset{0}{\overset{{\vartheta }_1}{\int }}\alpha \sin \alpha d\alpha$$

(5)

The intersection point $c$ of the rectangle and the arc part moves along the centerline of the rectangle by the length $b$ to point $c_1$ , and the sector area $w_2$ of this point is:

$$w_2=b{r}_3+{\vartheta }_1{r}_3^2$$

(6)

The distance $y_2$ from point $c_1$ to the z-axis is:

$$y_2=r_3\sin ({\vartheta }_1)+b\sin (\eta )$$

(7)

The main sectoral moment $J_{\mathrm{zw}2}$ of inertia of the rectangular part is:

$$\begin{array}{c}\multicolumn{1}{c}{J_{\mathrm{zw}2}=\underset{A}{\int }{y}_2{w}_2\mathrm{dA}=L_1{\int }_0^l(b{r}_3+{\vartheta }_1{r}_3^2)}\\ \multicolumn{1}{c}{(\sin {\vartheta }_1{r}_3+\sin \eta b)L_1\mathrm{db}}\end{array}$$

(8)

The main sectoral moment $J_{\mathrm{zw}}$ of inertia of the entire section is:

$$J_{\mathrm{zw}}=2J_{\mathrm{zw}1}+2J_{\mathrm{zw}2}$$

(9)

The moment of inertia of the entire section with respect to the z-axis is:

$$\begin{array}{c}{J}_z=\underset{A}{\int }{y}^2\mathrm{dA}=\frac{2}{3}{L}_1{\cos }^2({\vartheta }_1){L_0}^3+r_3{L}_1{L_0}^2\sin (2{\vartheta }_1)\hfill \\ +\frac{2}{3}{L}_1{L}_0(3{r_3}^2+\frac{t^2}{4}){\sin }^2({\vartheta }_1)+L_1{r_3}^3\hfill \\ \multicolumn{1}{c}{({\vartheta }_1-\cos ({\vartheta }_1)\sin ({\vartheta }_1))}\end{array}$$

(10)

The value of $e$ can be obtained from equation $e=\frac{J_{\mathrm{zw}}}{J_z}$ , as shown in Figure 5. The position of the main sector pole is obtained by moving the reference sector pole $O$ to the left by a distance $e$ .

Figure 5.

Schematic diagram of the solution of the moment of inertia of the main sector.

After determining the poles and zeros of the main sector, the main sector area $w_M$ at point $M$ can be obtained as:

$$w_M=2(S_{\mathrm{\Delta }\mathrm{AOM}}-S_{\mathrm{\Delta }O{M}_{0}M})=e{r}_3\sin (\alpha )-\alpha r_3^2$$

(11)

Then the main sector moment $J_{w1}$ of inertia of the arc part is:

$$J_{w1}=\underset{A}{\int }{w_M^2\mathrm{dA}={\int }_0^{{\vartheta }_1}(e{r}_3\sin \alpha -\alpha r_3^2)}^2{L}_1{r}_{3}d\alpha$$

(12)

Analyzing the rectangular part, the main sector area $w_{c1}$ of point $c_1$ is:

$$\begin{array}{c}{w}_{c1}=S_{\mathrm{\Delta }\mathrm{AM}{N}_0}-S_{\mathrm{NC}{C}_1}=S_{\mathrm{\Delta }\mathrm{AO}{C}_1}-S_{\mathrm{OC}{M}_0}-S_{\mathrm{OC}{C}_1}\hfill \\ =e(r_3\sin {\vartheta }_1+b\sin \eta )-{\vartheta }_1{r}_3^2-b{r}_3\hfill \end{array}$$

(13)

Then the main sector moment $J_{w2}$ of inertia of the rectangular part is:

$$J_{w2}={\underset{A}{\int }{w}_{{\vartheta }_1}^{2}dA={\int }_0^l(e(r_3\sin {\vartheta }_1+b\sin \eta )-{\vartheta }_1{r}_3^2-b{r}_3)}^{2}db$$

(14)

Then the main sectoral moment $J_w$ of inertia of the entire section is:

$$J_w=2(J_{w1}+J_{w2})$$

(15)

After obtaining $J_w$ , the constraint torsion coefficient $k$ can be obtained from the above equation (2). Afterward, based on the initial parameter solution method of constrained torsion, ¹⁵ firstly, the coordinate system is selected at any section, and the initial torsion angle on the section with z = 0 is set as ${\eta }_0$ , the initial torsion rate is $\eta {\prime }_0$ , the double initial moment is $B_0$ , and the initial torque is $T_0$ . Then the initial state vector $\left\{Z_0\right\}$ of the constrained torsion of the thin-walled member can be obtained as:

$$\left\{z_0\right\}=\{0,0,-\frac{B_0}{G{J}_d},\frac{T_0}{G{J}_d}\}$$

(16)

When $z=l_0$ , the boundary conditions are:

$${{\eta }_{l_0}}^{\prime }=0;T_0=T$$

(17)

According to the constraint torsion equation, so there is:

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}{\eta }_{l_0}\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{-\frac{B_{l_0}}{G{J}_d}}\\ \multicolumn{1}{c}{\frac{T}{G{J}_d}}\end{array}\right\}=\left[\begin{array}{cccc}1\hfill & \frac{l_0}{k}\mathrm{shk}\hfill & \mathrm{chk}-1\hfill & l_0-\frac{l_0}{k}\mathrm{shk}\hfill \\ \multicolumn{1}{c}{0}& \mathrm{chk}\hfill & \frac{k}{l_0}\mathrm{shk}\hfill & 1-\mathrm{chk}\hfill \\ \multicolumn{1}{c}{0}& \frac{l_0}{k}\mathrm{shk}\hfill & \mathrm{chk}\hfill & -\frac{l_0}{k}\mathrm{shk}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]}\\ \multicolumn{1}{c}{\times \left\{\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{-\frac{B_0}{G{J}_d}}\\ \multicolumn{1}{c}{\frac{T_0}{G{J}_d}}\end{array}\right\}}\end{array}$$

(18)

In equation (18): ${\eta }_{l_0}$ is the torsion angle when the beam length is $l_0$ , $T_{l_0}$ is the torque on the beam when the beam length is $l_0$ , $B_{l_0}$ is the double initial moment on the beam when the beam length is $l_0$ , $\mathrm{sh},\mathrm{ch}$ respectively refer to the hyperbolic sine function and the hyperbolic cosine function.

According to the above equation, it is easy to solve:

$${\eta }_{l_0}=\frac{T{l}_0{(1-\mathrm{chk})}^2}{G{J}_d\mathrm{kshk}}+\frac{T_{l_0}}{G{J}_d}(1-\frac{\mathrm{shk}}{k})$$

(19)

Therefore, the torsion stiffness $K_{\varphi \mathrm{rT}}$ of the beam is obtained:

$$K_{\varphi \mathrm{rT}}=\frac{T}{{\eta }_{l_0}}=\frac{G{J}_d\mathrm{kshk}}{l_0(1-\mathrm{chk})+l_0\mathrm{shk}(k-\mathrm{shk})}$$

(20)

Calculation of the torsion angle of the beam

From the equation (1) in Section 3.1, it can be seen that in order to obtain the torsion stiffness of the beam, the torsion angle of the beam needs to be deduced. The motion of the wheel of the torsion beam suspension is equivalent to swing motion of the rotation axis with the line vector connecting the body connection point $O_1$ and the shear center of the beam symmetry plane (the main fan pole $O_2$ obtained above), as shown in Figure 6.

Figure 6.

Equivalent kinematics model.

Before calculating the torsion angle of the beam, an initial vector $P{\prime }_t$ should be defined in the local coordinate system, and solve the coordinate $\varsigma$ of the wheel center in the global coordinate system and the relationship between the wheel travel and the angle $\zeta$ of the rotation of one wheel relative to the rotation axis.

$$\begin{array}{c}{P\prime }_t=O-O_1=\left[\begin{array}{c}{P\prime }_{\mathrm{tx}}\hfill \\ \multicolumn{1}{c}{{P\prime }_{\mathrm{ty}}}\\ \multicolumn{1}{c}{{P\prime }_{\mathrm{tz}}}\end{array}\right]=\left[\begin{array}{c}{O}_x\hfill \\ \multicolumn{1}{c}{O_y}\\ \multicolumn{1}{c}{O_z}\end{array}\right]-\left[\begin{array}{c}{O}_{1x}\hfill \\ \multicolumn{1}{c}{O_{1y}}\\ \multicolumn{1}{c}{O_{1z}}\end{array}\right]\hfill \\ \varsigma =\left[\begin{array}{c}{\varsigma }_x\hfill \\ \multicolumn{1}{c}{{\varsigma }_y}\\ \multicolumn{1}{c}{{\varsigma }_z}\end{array}\right]=O_1+A{P\prime }_t\hfill \end{array}$$

(21)

$A$ in equation (21) is the change matrix that transforms the local coordinate system into the global coordinate system. From the equation (21), the component ${\varsigma }_z$ of the wheel center in the vertical direction in the global coordinate system is gained.

$$\begin{array}{c}{\varsigma }_z=O_{1z}+2({\lambda }_1{\lambda }_3-{\lambda }_2{\lambda }_0){P\prime }_{\mathrm{tx}}+2({\lambda }_2{\lambda }_3+{\lambda }_1{\lambda }_0){P\prime }_{\mathrm{ty}}\hfill \\ +({\lambda }_0^2-{\lambda }_1^2-{\lambda }_2^2+{\lambda }_3^2){P\prime }_{\mathrm{tz}}\hfill \end{array}$$

(22)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{d{\varsigma }_z}{d\zeta }=[w_x{w}_z\sin \zeta -w_y\cos \zeta ]{P\prime }_{\mathrm{tx}}}\\ \multicolumn{1}{c}{+[w_y{w}_z\sin \zeta +w_x\cos \zeta ]{P\prime }_{\mathrm{ty}}}\\ \multicolumn{1}{c}{+[w_z^2+\sin \zeta -\sin \zeta ]{P\prime }_{\mathrm{tz}}}\end{array}$$

(23)

From the same reasoning above, it can be seen that in the reverse travel, the wheels on both sides are equivalent to the line vector $t$ connecting the connection point of the body and the shear center of the cross-beam section, which is the swing motion of the rotating shaft. So a local coordinate system is established at point $O_1$ . Based on the method of Euler quaternion. ¹⁶ it is also assumed that the Euler angle of the swing axis of the suspension in the coordinate system is $w_x,w_y,w_z$ , and the Euler angle vector of the rotation axis is $t$ ,so there is:

$$t=[w_x,w_y,w_z{]}^T$$

(24)

The Euler four-element vector is:

$$\begin{array}{c}P″=[{\lambda }_0,{\lambda }_1,{\lambda }_2,{\lambda }_3{]}^T\hfill \\ =[\cos (\frac{\zeta }{2}),w_x\sin (\frac{\zeta }{2}),w_y\sin (\frac{\zeta }{2}),w_z\sin (\frac{\zeta }{2}){]}^T\hfill \end{array}$$

(25)

In equation (25), ${\lambda }_0,{\lambda }_1,{\lambda }_2,{\lambda }_3$ is Euler’s quaternion.

The change matrix $A$ is:

$$A=\left[\begin{array}{ccc}{\lambda }_0^2+{\lambda }_1^2-{\lambda }_2^2-{\lambda }_3^2\hfill & 2({\lambda }_1{\lambda }_2-{\lambda }_0{\lambda }_3)\hfill & 2({\lambda }_1{\lambda }_3+{\lambda }_0{\lambda }_2)\hfill \\ \multicolumn{1}{c}{2({\lambda }_1{\lambda }_2+{\lambda }_0{\lambda }_3)}& \lambda -{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2\hfill & 2({\lambda }_2{\lambda }_3-{\lambda }_0{\lambda }_1)\hfill \\ \multicolumn{1}{c}{2({\lambda }_1{\lambda }_3-{\lambda }_0{\lambda }_2)}& 2({\lambda }_2{\lambda }_3+{\lambda }_0{\lambda }_1)\hfill & {\lambda }_0^2-{\lambda }_1^2-{\lambda }_2^2+{\lambda }_3^2\hfill \end{array}\right]$$

(26)

Therefore, when the wheel on one side of the torsion beam suspension travels up by angle $\zeta$ around the rotation axis, the vector $P\prime$ in the local coordinate system can be converted into $P$ in the global coordinate system through the above change matrix, that is,

$$P=\mathrm{AP}\prime$$

(27)

A vertical unit vector is defined as $P{\prime }_z=[1,0,0{]}^T$ in the local coordinate system. And from the above change matrix, it is known that when the wheel rotates $\zeta$ angles around the rotation axis, the vector in the global coordinate system is $P_z=AP{\prime }_z$ , then the coordinate $P_z$ of the vertical unit vector is:

$$P_z=\left[\begin{array}{c}{p}_{\mathrm{zx}}\hfill \\ \multicolumn{1}{c}{P_{\mathrm{zy}}}\\ \multicolumn{1}{c}{P_{\mathrm{zz}}}\end{array}\right]=\left[\begin{array}{c}{\lambda }_0^2+{\lambda }_1^2-{\lambda }_2^2-{\lambda }_3^2\hfill \\ \multicolumn{1}{c}{2({\lambda }_1{\lambda }_2+{\lambda }_0{\lambda }_3)}\\ \multicolumn{1}{c}{2({\lambda }_1{\lambda }_3-{\lambda }_0{\lambda }_2)}\end{array}\right]$$

(28)

Combining the above equation and the geometric relationship in Figure 7, the torsion angle of the beam can be obtained:

$$\eta =2\mathrm{ac}\tan \left(\frac{P_{\mathrm{zz}}}{P_{\mathrm{zx}}}\right)=2\mathrm{ac}\tan \left[\frac{2({\lambda }_1{\lambda }_3-{\lambda }_2{\lambda }_0)}{{\lambda }_0^2+{\lambda }_1^2-{\lambda }_2^2-{\lambda }_3^2}\right]$$

(29)

Figure 7.

Schematic diagram of the geometric relationship of the torsion angle of the beam.

Calculation of the roll stiffness of the suspension

When the wheels on both sides bounce in the opposite direction, it is equivalent to have two forces with the same amplitude but opposite directions acting on the center point of the wheels. The energy generated is mainly absorbed by the deformation of the beam and the bushing, and since the stiffness of the trailing arm is much larger than the rest of the components, the deformation of the trailing arm can be ignored. Typically, rubber bushings have stiffness in six directions and are mounted on the body at an angle to increase the understeer tendency in Figures 2 and 3.

When considering the effect of the rubber bushing on the roll stiffness of the suspension, it is assumed that the torsion angle ${\theta }_F$ of the bushing about the axial direction and the torsion angle $\eta$ of the beam, and the torsion angle ${\theta }_T$ of the bushing about the radial direction and the camber angle $\gamma$ of the wheel have the following relationship. ¹⁷

$${\theta }_F=\frac{\eta }{2}$$

(30)

$${\theta }_T=\gamma$$

(31)

Assuming that the external force acting on the wheel center is $F_p$ , according to the law of conservation of energy, the energy $E_p$ required by the carriage should be the energy left after the beam and the rubber bushing absorb it, which can be expressed as follows:

$$E_p=\frac{1}{2}{K}_{\varphi \mathrm{rT}}{\eta }^2+\frac{5}{2}{K}_F{\theta }_F^2+\frac{5}{2}{K}_T{\theta }_T^2-2F_p{\varsigma }_z$$

(32)

Where, ${\varsigma }_z$ is the vertical displacement of the wheel center.

In the initial state, $E_p$ should satisfy the following equation:

$$\begin{array}{c}\multicolumn{1}{c}{\frac{d{E}_p}{d{\varsigma }_z}=K_{\varphi \mathrm{rT}}\eta \frac{d\eta }{d{\varsigma }_z}+2K_F{\theta }_F\frac{d{\theta }_F}{d{\varsigma }_z}+2K_T{\theta }_T\frac{d{\theta }_T}{d{\varsigma }_z}}\\ \multicolumn{1}{c}{-2F_{p}d{\varsigma }_z=0}\end{array}$$

(33)

When the wheel travel is not large, the torsion angle of the beam and the rubber bushing is still in the linear range, that is:

$$\eta =\frac{d\eta }{d{\varsigma }_z}{\varsigma }_z$$

(34)

The angle of rotation of the rubber bushing around the axial direction ${\theta }_F$ and around the radial direction ${\theta }_T$ is:

$$\{\begin{array}{c}{\theta }_F=\frac{1}{2}\frac{d\eta }{d{\varsigma }_z}{\varsigma }_z\hfill \\ \multicolumn{1}{c}{{\theta }_T=\frac{1}{2}\frac{d\gamma }{d{\varsigma }_z}{\varsigma }_z}\end{array}$$

(35)

Combining equations (30)–(35), the equivalent spring stiffness $K_{\varphi \mathrm{rs}}$ of the torsion beam suspension can be obtained as:

$$\begin{array}{c}\multicolumn{1}{c}{K_{\varphi \mathrm{rs}}=\frac{F_p}{{\varsigma }_z}=\frac{1}{2}}\\ \multicolumn{1}{c}{[K_{\varphi \mathrm{rT}}{\left(\frac{d\eta }{{\varsigma }_z}\right)}^2+2K_F{\left(\frac{d{\theta }_F}{{\varsigma }_z}\right)}^2+2K_T{\left(\frac{d{\theta }_T}{{\varsigma }_z}\right)}^2]}\end{array}$$

(36)

where $K_{\varphi \mathrm{rT}}$ is the torsion stiffness of the beam, which has been derived from Section 3.1.

Therefore, the roll stiffness $K_{\varphi r}$ of the torsion beam suspension can be expressed by the equivalent spring stiffness $K_{\varphi \mathrm{rs}}$ and wheel base $L$ , namely:

$$K_{\varphi r}=\frac{1}{2}{K}_{\varphi \mathrm{rs}}{L}^2$$

(37)

Bushing installation angle and bushing stiffness on suspension C characteristics

The C characteristic of the suspension, namely compliance kinematics characteristic, describes the change of the wheel alignment parameters caused by the deformation of the bushing and the connection point of the suspension in the process of transmitting the force and moment of the suspension. The suspension simulation analysis of the established torsion beam suspension is carried out, and the relationship between the toe angle and the roll stiffness under the action of lateral force under different bushing installation angles is obtained as shown in Figure 8.

Figure 8.

Influence of different bushing installation angles on suspension characteristics under lateral force conditions: (a) influence on the toe angle and (b) influence on the roll stiffness.

Under the lateral force condition of the torsion beam, the axle steering will occur, and the axle steering is mainly caused by the deformation of the rubber bushing at the connection between the trailing arm and the body, thus causing the C characteristics of the suspension to change, which in turn has a major impact on the handing stability of the rear axle. The change in toe angle due to bushing deformation under lateral force can be decomposed into the following two effects:

1. Due to the axial translation stiffness of the bushing, when a lateral force acts, the connection point between the torsion beam and the body moves in series along the bushing axis, resulting in a certain understeer of the rear axle. The lower the axial stiffness of the bushing, the more obvious the tendency of understeering.

2. Due to the radial translation stiffness of the bushing, when a lateral force acts, the connection point between the torsion beam and the body is compressed along the radial direction of the bushing, resulting in a certain understeer of the rear axle. The lower the axial stiffness of the bushing, the more obvious the tendency of oversteer.

In order to ensure the vehicle has better handling stability characteristics, when the rear suspension moves under lateral force conditions, the rear suspension should produce a change direction of the toe angle in the same direction as the lateral force, that is, appropriate understeer characteristics, which is called follow-up steering. ¹⁸ Therefore, the change gradient of the toe angle should be positive, that is, the slope of the curve is positive.

Usually for the rear suspension, if the change gradient of the toe angle is positive, it means that when the body rolls, the toe angle of the outer wheel increases, the wheel turns inward, the toe angle of the inner wheel decreases, and the wheel turns outward, then will beneficial to understeer. Conversely, if the change gradient of the toe angle is negative, it will be unfavorable for understeer.

It can be seen from Figure 8(a) that when the lateral force changes from −4000 N to 4000 N, when the bushing installation angle is within [−45°−45°], the toe angle changes regardless of whether the angle changes from small to large or from large to small, the change gradient of the toe angle is positive, which is conducive to understeer. It can be seen that it is consistent with the above analysis theory.

It can be known from equation (35) when the bushing installation angle changes (the torsion angle ${\theta }_F$ of the bushing about the axial direction and the torsion angle ${\theta }_T$ of the bushing about the radial direction), $2K_F{\left(\frac{d{\theta }_F}{{\varsigma }_z}\right)}^2+2K_T{\left(\frac{d{\theta }_T}{{\varsigma }_z}\right)}^2$ in equation (36) will change, and the stress on the bushing will change, which will cause the equivalent spring stiffness $K_{\varphi \mathrm{rs}}$ to change, which further and the stress on the bushing will change, which will cause the equivalent spring stiffness $K_{\varphi \mathrm{rs}}$ to change, which further causes the roll stiffness $K_{\varphi r}$ in equation (44) to change, and finally leads to the change in Figure 8(b).

Sensitivity analysis of influence of bushing stiffness on suspension characteristic index

Sensitivity analysis refers to the response quantity of the system dynamic change caused by the parameter change, that is, the influence of the parameter change on the time response, state vector, transfer function or other dynamic performance of the system. Since the suspension model is left and right symmetrical, the influence laws of the bushings on both sides are exactly the same. Therefore, the stiffness in six directions of the rubber bushing at the connection between the trailing arm and the vehicle body on one side of the model is used as the design variable to analyze the influence on the suspension characteristic index, and the sensitivity analysis results shown in Table 1 are obtained. Among them, the bushing stiffness in six directions is recorded as DV1-DV6 in turn, and the positive and negative signs of the sensitivity value indicate the changes of each variable and the performance index of the suspension are the same or opposite. The change gradient refers to the slope of the curve of the suspension performance index changing with the lateral force when the lateral force is [−4000 N, 4000 N].

Table 1.

Sensitivity analysis results of bushing stiffness.

Designvariable	Toe anglegradient	Camber anglegradient	Roll centerheight gradient	Roll stiffnessgradient
DV1	−2.53	2.25	2.86	1.49
DV2	−1.68	−2.88	2.17	0.93
DV3	2.07	−2.17	−1.64	2.74
DV4	1.46	−1.74	−2.23	2.20
DV5	0.79	1.43	−1.62	1.62
DV6	−2.05	1.80	1.52	0.79

From Table 1, the bushing stiffness in different directions that are most sensitive to the suspension characteristics can be screened out. According to the experience of engineering design, the bushing stiffness is reduced or increased to 10% or 20% of the initial value, respectively, and then the suspension simulation is carried out, and the relationship between the suspension characteristic index and the lateral force is obtained, as shown in Figure 9.

Figure 9.

Influence of different bushing stiffnesses on suspension characteristics under lateral force conditions: (a) influence on the toe angle, (b) influence on camber angle, (c) influence on the roll center height, and (d) influence on roll stiffness.

Under the force of the wheel, the bushing will be deformed, and the deformation of the bushing will lead to the change of the wheel alignment parameters, which will eventually affect the C characteristics of the suspension, resulting in the change relationship shown in Figure 9. It can be seen from Figure 9 that when the wheel is under the action of lateral force, the bush stiffness varying in the range of [−20%, 20%] has obvious effects on the wheel toe angle, camber angle, suspension roll center height, and roll stiffness.

Establishment of theoretical model of handling stability considering suspension roll

Figure 10 shows the established vehicle handling dynamics model, $X\prime -Y\prime -Z\prime$ is the absolute coordinate system. The origin of the centroid coordinate system $X-Y-Z$ is located at the intersection of the vertical line passing through the centroid and the roll axis.

Figure 10.

Three-degree-of-freedom manipulation model.

The magnitude of the absolute acceleration $a$ at the center of mass $O$ is:

$$a=|\stackrel{·}{\phi }\times V|=\stackrel{·}{\phi }V$$

(38)

In equation (38), the velocity at the center of mass $O$ is $V$ , and $\phi$ is the angle between $V$ and $X\prime$ .

The direction of the acceleration vector $a$ is perpendicular to the velocity vector $V$ , and the component on the y-axis $\stackrel{··}{Y}$ is:

$$\stackrel{··}{Y}=\stackrel{·}{\phi }V\cos \beta$$

(39)

where $\beta$ is the side slip angle. Since $\phi =\psi +\beta$ and $\cos \beta \approx 1$ , equation (39) can be transformed into:

$$\stackrel{··}{Y}=V\stackrel{·}{\beta }+\mathrm{Vr}$$

(40)

In equation (40), $r$ is the yaw rate.

The component ${\stackrel{··}{Y}}_s$ of the center of mass acceleration of the suspension sprung mass on the y-axis is:

$${\stackrel{··}{Y}}_s=V\stackrel{·}{\beta }+\mathrm{Vr}-h_s\stackrel{·}{p}$$

(41)

In equation (41), $\stackrel{·}{p}$ is the roll angle acceleration, $h_s$ is the distance from the sprung mass center to the roll axis.

The balance equation for the lateral force is:

$$\begin{array}{c}\multicolumn{1}{c}{(M-M_s)(r+\stackrel{·}{\beta })V+M_s[(r+\stackrel{·}{\beta })V-h_{s}p]}\\ \multicolumn{1}{c}{=\sum _{i=1}^n({F\prime }_{\mathrm{yi}}+{F″}_{\mathrm{yi}})}\end{array}$$

(42)

In equation (42), $M$ is the vehicle mass, $M_s$ is the sprung mass.

The balance equation for the moment about the z-axis is:

$$I_s\stackrel{·}{r}-I_{\mathrm{xz}}\stackrel{·}{p}=\sum _{i=1}^n{x}_i({F\prime }_{\mathrm{yi}}+{F″}_{\mathrm{yi}})$$

(43)

In equation (43), $I_s$ is the yaw moment of inertia, $I_{\mathrm{xz}}$ is the inertia product.

The balance equation for the moment about the x-axis is:

$$I_{\varphi }\stackrel{·}{p}-I_{\mathrm{xz}}r-M_s{h}_{s}V(\stackrel{·}{\beta }+r)=(-K_{\varphi }+M_{s}g{h}_s)\varphi -C_{\varphi }p$$

(44)

In equation (44), $I_{\varphi }$ is the roll moment of inertia, $C_{\varphi }$ is the roll angle damping.

The side slip angles $\beta {\prime }_i$ and ${\beta ″}_i$ of the inner and outer wheels are:

$$\beta {\prime }_i=\frac{V\beta +d_{i}r}{V-L_{i}r/2}-{\delta }_i=\frac{\beta +d_{i}r/V}{1-L_{i}r/2V}-{\delta }_i\approx \beta +d_{i}r/V-{\delta }_i$$

(45)

$${\beta ″}_i=\frac{V\beta +d_{i}r}{V+L_{i}r/2}-{\delta }_i=\frac{\beta +d_{i}r/V}{1+L_{i}r/2V}-{\delta }_i\approx \beta +d_{i}r/V-{\delta }_i$$

(46)

In equations (45) and (46), ${\delta }_i$ is the wheel steering angle, $d_i$ is the distance from the center of mass to the front and rear axles, and $L_i$ is the wheel track of front and rear wheels.

The relationship between the roll steering angle $\alpha {\prime }_i,{\alpha ″}_i$ of the front and rear wheels and the roll angle is. ¹⁹

$$\alpha {\prime }_i={\alpha ″}_i=\frac{\partial {\alpha }_i}{\partial \varphi }\varphi$$

(47)

In equation (47), $\partial {\alpha }_i/\partial \varphi$ is the roll steering angle corresponding to each unit body roll angle.

Supposing that the camber angle caused by the body roll is proportional to the body roll angle, the camber thrusts $Y{\prime }_{\mathrm{ci}}$ and ${Y″}_{\mathrm{ci}}$ generated by the front and rear wheels are:

$$Y{\prime }_{\mathrm{ci}}={Y″}_{\mathrm{ci}}=-K_{\mathrm{ci}}\frac{\partial {\sigma }_i}{\partial \varphi }\varphi$$

(48)

In equation (48), $K_{\mathrm{ci}}(i=f,r)$ is the camber lateral thrust coefficient of the front and rear tires, and $\partial {\sigma }_i/\partial \varphi$ is the camber angle generated per unit body roll angle.

Considering the roll steering and roll camber of the suspension, and combining equations (45)–(48), the tire lateral forces $F{\prime }_{\mathrm{yi}}$ and ${F″}_{\mathrm{yi}}$ acting on the wheels are:

$$F{\prime }_{\mathrm{yi}}={F″}_{\mathrm{yi}}=-K_i(\beta +d_{i}r/V-{\delta }_i-\frac{\partial {\alpha }_i}{\partial \varphi }\varphi )-K_{\mathrm{ci}}\frac{\partial {\sigma }_i}{\partial \varphi }\varphi$$

(49)

In equation (49), $K_i(i=f,r)$ is the side slip stiffness of the front and rear tires.

When the vehicle rolls, the front and rear wheels will produce corresponding slip angles due to the roll of the vehicle, which are denoted as $\frac{\partial {\alpha }_f}{\partial \varphi }\varphi ,\frac{\partial {\alpha }_r}{\partial \varphi }\varphi$ respectively. The front wheel itself has a rotation angle ${\delta }_f$ input by the steering wheel, thus the side slip angles of the front and rear wheels can be expressed as:

Side slip angle of front wheel:

$${\delta }_f+\frac{\partial {\alpha }_f}{\partial \varphi }\varphi -\beta -\frac{d_f\stackrel{·}{r}}{V}$$

(50)

Side slip angle of rear wheel:

$$\frac{\partial {\alpha }_r}{\partial \varphi }\varphi -\beta +\frac{d_r\stackrel{·}{r}}{V}$$

(51)

Among them, $\partial {\alpha }_f/\partial \varphi ,\partial {\alpha }_r/\partial \varphi$ means that the front wheel and the rear wheel produce the corresponding side slip angle when the vehicle produces a unit roll angle.

Equations (45) and (46) can calculate the roll steering angle and camber thrust, and after sorting out equations (42)–(44), the 3-DOF dynamic equations can be obtained as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{MV}\stackrel{·}{\beta }+2\sum _{i=1}^n{K}_i\beta +[\mathrm{MV}+\frac{2\sum _{i=1}^n{d}_i{K}_i}{V}]}\\ \multicolumn{1}{c}{r-M_s{h}_s\stackrel{·}{p}-2Y_{\varphi }\varphi =2K_f\delta }\end{array}$$

(52)

$$2\sum _{i=1}^n{d}_i{K}_i\beta +I_s\stackrel{·}{r}+\frac{2\sum _{i=1}^n{d}_i^2{K}_i}{V}r-I_{\mathrm{xz}}\stackrel{·}{p}-2N_{\varphi }\varphi =2d_f{K}_f\delta$$

(53)

$$\begin{array}{c}\multicolumn{1}{c}{-M_s{h}_{s}V\stackrel{·}{\beta }-I_{\mathrm{xz}}\stackrel{·}{r}-M_s{h}_s\mathrm{Vr}+I_{\varphi }\stackrel{·}{p}+C_{\varphi }p}\\ \multicolumn{1}{c}{+(K_{\varphi }-M_{s}g{h}_s)\varphi =0}\end{array}$$

(54)

In equations (52)–(54), there is:

$$\{\begin{array}{c}{Y}_{\varphi }=\sum _{i=1}^n{K}_i\frac{\partial {\alpha }_i}{\partial \varphi }-\sum _{i=1}^n{K}_{\mathrm{ci}}\frac{\partial {\sigma }_i}{\partial \varphi }\hfill \\ \multicolumn{1}{c}{N_{\varphi }=\sum _{i=1}^n{d}_i{K}_i\frac{\partial {\alpha }_i}{\partial \varphi }-\sum _{i=1}^n{d}_i{K}_{\mathrm{ci}}\frac{\partial {\sigma }_i}{\partial \varphi }}\end{array}$$

(55)

Then the Laplace change is performed on the equations (52)–(54), and the state equation system after the Laplace change can be obtained, and then the transfer function of the corresponding response indexes relative to the steering wheel angle can be obtained, and then then let $s=j\omega$ , where $j=\sqrt{-1}$ and finally the corresponding frequency characteristic index can be obtained.

Based on the above theories, it can be concluded that the influence mechanism of the bushing installation angle and bushing stiffness of the torsion beam suspension on the frequency response characteristics is as follows:

1. The difference in the bushing installation angle will affect the roll stiffness of the suspension, which in turn affects the roll stiffness of the suspension, thereby affecting the frequency response characteristics.

2. The camber angle of the front and rear wheels in the above equation (48) and the side slip angles of the front and rear wheels in the equations (50) and (51) only take into account the additional turning angle due to the suspension movement. Only the additional turning angle due to the suspension movement is considered, and the influence of the deformation of the suspension bushing, etc. on the side slip angle and camber angle of the front and rear wheels is not considered. It can be seen from the above analysis that when the bushing installation angle and the bushing stiffness are different, the C characteristics of the suspension will change, which will affect the side slip stiffness of the front and rear tires $(K_i)$ , the side slip angle of the front and rear wheels $({\delta }_f+\frac{\partial {\alpha }_f}{\partial \varphi }\varphi -\beta -\frac{d_f\stackrel{·}{r}}{V})$ , $(\frac{\partial {\alpha }_r}{\partial \varphi }\varphi -\beta +\frac{d_r\stackrel{·}{r}}{V})$ , and the camber angle of the front and rear wheels $(\frac{\partial {\sigma }_i}{\partial \varphi }\varphi )$ , which ultimately affects the frequency response characteristics.

Bushing installation angle on steady-state and transient handing stability of vehicle

The vehicle multi-body dynamics model of the front McPherson suspension and the rear torsion beam suspension is established from a real vehicle parameter. The parameters of the established vehicle model in Table 2.

Table 2.

Vehicle parameters.

Parameters	Value	Parameters	Value
Sprung mass (kg)	1097.00	Height from center of mass to ground (m)	0.37
Unsprung mass (kg)	551.98	Distance between front wheels (m)	1.59
Center of mass to front wheelbase (m)	0.71	Distance between rear wheels (m)	1.52
Center of mass to rear wheelbase (m)	1.85	Wheelbase (m)	2.56
Spring stiffness of front suspension (N/m)	2.10 × 104	Spring stiffness of rear suspension (N/m)	2.27 × 104

Analysis of transient handing stability index

The indicators of transient handing stability include yaw rate gain, roll angle gain, delay time of yaw rate (DT2), and delay time of lateral acceleration (DT1).

Among then, the yaw rate gain refers to the gain value of the yaw rate relative to the steering wheel angle, and the roll angle gain refers to the gain value of the roll angle relative to the lateral acceleration. The delay time of the yaw rate refers to the delay time of the yaw rate relative to the lateral acceleration, and the delay time of the lateral acceleration refers to the delay time of the lateral acceleration relative to the steering wheel angle. The resonance frequency refers to the frequency corresponding to maximum gain of yaw rate relative to steering wheel angle.

Figure 11(a) is the curve of the yaw rate gain. In Figure 11(a), when the low frequency is around 0.5 Hz, the change of the bushing installation angle has a significant impact on the yaw rate gain. The gain value does not change according to a certain law with the change of the bushing installation angle. Generally, the driver’s frequency does not exceed 1 Hz when manipulating the vehicle, therefore the transient stability index changes at the high frequencies are not considered for the time being.

Figure 11.

Influence of bushing installation angle on transient handing stability:(a) Influence on yaw rate gain, (b) Influence on roll angle gain, (c) Influence on the delay time of lateral acceleration, (d) Influence on delay time on yaw rate.

Figure 11(b) is the curve of the roll angle gain. In Figure 11(b), when the low frequency is around 0.5 Hz, the change of the bushing installation angle has little effect on the roll angle gain. The gain value does not show a certain regular change with the change of the bushing installation angle.

Figure 11(c) and (d) are curves of delay time of lateral acceleration and delay time of yaw rate. In Figure 11(c) and (d), when the low frequency is around 0.5 Hz, the change of the bushing installation angle has a significant impact on the delay time of lateral acceleration and the delay time of yaw rate. The delay time does not change according to a certain law with the change of the bushing installation angle.

Analysis of index of steady-state handing stability

Steady-state characteristic analysis adopts the constant radius cornering simulation, and the variation of steady-state handling index with bushing installation angle in the range of [−90° to 90°] is obtained through simulation calculation in Figure 12. The steady-state handling stability indicators include: maximum lateral acceleration, understeer, the side slip angle of lateral acceleration at 0.7g, and steering wheel angle of lateral acceleration at 0.7g.

Figure 12.

Influence of bushing installation angle on steady-state handling index: (a) influence on the maximum lateral acceleration, (b) influence on understeer, (c) influence on the side slip angle, and (d) influence on steering wheel angle.

The maximum lateral acceleration and the steering wheel angle of the lateral acceleration at 0.7g can be obtained from the slope of the curve of the steering wheel angle relative to the lateral acceleration; Understeer refers to the steering wheel angle of lateral acceleration at 0.4g, which can be obtained from the slope of the curve of the steering wheel angle relative to the lateral acceleration at 0.4g; The side slip angle of lateral acceleration at 0.7g can be obtained from the slope of the curve of the side slip angle of lateral acceleration at 0.7g.

In Figure 12(a), the maximum lateral acceleration generally shows a trend of first decreasing, then increasing, then decreasing and finally increasing with the increase of the bushing opening angle; In Figure 12(b) and (d), the understeer and the steering wheel angle of lateral acceleration shows a trend of first decreasing and then increasing with the increase of the bushing opening angle at 0.7g generally first decreases and then increases with the increase of the bushing opening angle. In Figure 12(c), the side slip angle of the lateral acceleration at 0.7g generally In addition, in Figure 10, when the bushing installation angle is around −22.5°, the understeer of the vehicle and the steering wheel angle of the lateral acceleration at 0.7g reach the maximum, while the maximum lateral acceleration and side slip angle of lateral acceleration at 0.7g reach the minimum, indicating that the understeering characteristics of the vehicle are the best, and the steady-state performance is optimal near this range.

It can be seen from the above analysis that when the wheel is subjected to external force, the bushing is deformed, which in turn affects the C characteristics of the suspension, which further affects the steady-state performance index of the vehicle, and finally affects the steady-state performance of the vehicle.

Collaborative optimization of bushing stiffness and bushing installation angle

Firstly, the sensitivity analysis of the handling stability of the vehicle is executed by taking the rubber bushing at the connection between the trailing arm and the body as the object of sensitivity analysis. The importance of each objective is comprehensively weighed and the bushing stiffness that has a greater impact on the frequency characteristics is screened out. The sensitivity analysis results are in Table 3. It can be known the radial stiffness of the bushing in the x- and z-direction is more sensitive $(F_x,F_z)$ . It can be seen from Table 3 that the radial stiffness in the Z direction of the bushing is the most sensitive to the yaw rate gain, and the radial stiffness in the X direction of the bushing is the most sensitive to the roll angle gain. The most sensitive effect on the delay time of the yaw rate and the delay time of the lateral acceleration is the torsion stiffness of the bushing around the Z direction.

Table 3.

Sensitivity analysis results of handing stability.

Variable	Resonance frequency at 0.5 Hz	Yaw rate gain at 0.5 Hz	Roll angle gain at 0.5 Hz	DT 1 at 0.5 Hz	DT 2 at 0.5 Hz
DV1	0.86	−2.51	−2.20	−1.93	−1.63
DV2	1.32	1.18	−0.83	1.39	1.39
DV3	0.78	−2.14	−1.43	1.54	1.36
DV4	−0.64	0.74	0.54	−0.67	−0.81
DV5	0.47	0.56	1.16	−0.42	−0.61
DV6	0.75	−0.76	−0.72	0.46	0.48

In addition, the swept-sine input simulation test is carried out on the whole vehicle, and the influence of the change of the bushing stiffness on the frequency response characteristics is obtained. Through the analysis, it is known the bushing stiffness that changes in the range of [−20%, 20%] has little effect on the frequency characteristics of the vehicle when the low frequency is around 0.5 Hz. This paper gives the effect of the change of bushing stiffness on the delay time of lateral acceleration in Figure 13. It can be seen from Figure 13 that the bushing stiffness varying in the range of [−20%, 20%] has little effect on the lateral acceleration delay time at low frequencies around 0.5 Hz, and the absolute value of the delay time of the lateral acceleration at 0.5 Hz gradually decreases as the bushing stiffness changes from large to small.

Figure 13.

Effect of variation in bushing stiffness on delay time of lateral acceleration.

According to the results of sensitivity analysis, the stiffness in the direction with higher sensitivity of the bushing at the connection between the suspension trailing arm and the body and the bushing installation angle are taken as optimization variables, and generally the frequency of the driver steering the steering wheel is about 0.5 Hz, so the frequency characteristic index of the whole vehicle at 0.5 Hz is taken as the optimization target (Gain and delay time at 0.5 Hz) and the multi-objective optimization of the frequency characteristics is carried out through the genetic algorithm.

From the following description of formula (55), it can be known that the transfer function of the frequency characteristic index of the vehicle can be obtained from the equations (52)–(54), and make $s=j\omega$ , where $j=\sqrt{-1}$ . $s$ is the Laplace operator. For example, the transfer function of the yaw rate relative to the steering wheel angle $\frac{\stackrel{·}{\psi }}{\delta }(j\omega )$ , the transfer function of the roll angle relative to the lateral acceleration $\frac{\varphi }{a_y}(j\omega )$ ,where ${\alpha }_y=V(\stackrel{·}{\beta }+r)-h_s\stackrel{·}{p}$ , the resonance frequency $\frac{\stackrel{·}{\psi }}{\delta }(j\omega {)}_{max}$ (frequency corresponding to maximum gain of yaw rate relative to steering wheel angle), the phase lag of the yaw rate relative to the lateral acceleration ${\zeta }_{\stackrel{·}{\psi }}$ , and the phase lag of the lateral acceleration relative to the steering wheel angle ${\zeta }_{a_y}$ . Where, the phase lag angle can be converted into the delay time by equation (56).

$$t=\frac{\zeta +180}{360·f}\mathrm{or}\frac{\zeta -180}{360·f}$$

(56)

In equation (56), $\zeta$ is the phase lag angle, $t$ is the delay time, $f$ is the frequency.

Therefore, the mathematical expression of the optimization objective can be expressed as follows:

$$\{\begin{array}{c}{f}_1=min\left({|\frac{\stackrel{·}{\psi }}{\delta }(j\omega )|}_{\omega =0.5\mathrm{Hz}}\right)\hfill \\ \multicolumn{1}{c}{f_2=max\left({|\frac{\stackrel{·}{\psi }}{\delta }(j\omega )|}_{max}\right)}\\ \multicolumn{1}{c}{f_3=min{\left(\frac{{\zeta }_{a_y}\pm 180}{360·f}\right)}_{\omega =0.5\mathrm{Hz}}}\\ \multicolumn{1}{c}{f_4=min{\left(\frac{{\zeta }_{\stackrel{·}{\psi }}\pm 180}{360·f}\right)}_{\omega =0.5\mathrm{Hz}}}\\ \multicolumn{1}{c}{f_5=min\left({|\frac{\varphi }{a_y}(j\omega )|}_{\omega =0.5\mathrm{Hz}}\right)}\end{array}$$

(57)

In equation (57), $f_1$ is the yaw rate gain, $f_2$ is the resonance frequency, and $f_3,f_4$ are the delay time of lateral acceleration and the delay time of yaw rate, respectively, $f_5$ is the roll angle gain.

The constraints are as follows:

$$\{\begin{array}{c}0.5\le x_i\le 2,i=1,2\hfill \\ \multicolumn{1}{c}{-{90}^{°}\le {\alpha }_0\le {90}^{°}}\end{array}$$

(58)

In equation (58), $x_i$ is the constraint condition of the bushing stiffness, ${\alpha }_0$ is the constraint condition of the bushing installation angle.

In this paper, the non-dominated sorting genetic algorithm NSGA-II is used to solve the constrained multi-objective optimization problem. ²⁰ The algorithm is an improved multi-objective algorithm based on the NSGA algorithm. The algorithm introduces fast non-dominated sorting, crowding, and elite strategies. Through fast non-dominated sorting, the complexity of computing sorting is effectively reduced; The crowding degree and the comparison factor are introduced as the comparison criteria of the individuals in the population to ensure the diversity of the population; Compared with other genetic algorithms, it has better convergence and robustness, and has a good optimization effect on multi-objective optimization problems; The elite strategy is introduced to expand the sampling space, improve the global search ability, and the accuracy of optimization results.

The crossover operation and sudden mutation operation between individuals in the NSGA-II algorithm use the SBX method, and then generate sub-individuals according to this method. Formula (59) is the crossover calculation formula:

$$\{\begin{array}{c}{x_i}^{(1,t+1)}=\frac{1+{\beta }_{\mathrm{qi}}}{2}{x_i}^{(1,t)}+\frac{1+{\beta }_{\mathrm{qi}}}{2}{x_i}^{(2,t)}\hfill \\ \multicolumn{1}{c}{{x_i}^{(2,t+1)}=\frac{1+{\beta }_{\mathrm{qi}}}{2}{x_i}^{(1,t)}+\frac{1+{\beta }_{\mathrm{qi}}}{2}{x_i}^{(2,t)}}\end{array}$$

(59)

where, $x_i^{(1,t+1)}$ is the $i$ th gene of the first offspring individual $x^{(1,t+1)}$ , $x_i^{(2,t+1)}$ is the $i$ th gene of the second offspring individual $x^{(2,t+1)}$ , $x_i^{(1,t)}$ is the $i$ th gene of the first parent individual $x^{(1,t)}$ , $x_i^{(2,t)}$ is the $i$ th gene of the first parent individual $x^{(2,t)}$ , ${\beta }_{\mathrm{qi}}$ is a random number.

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

$${x_i}^{(1,t=1)}={x_i}^{(1,t)}+{\delta }_q({x_i}^{\mathrm{UB}}-{x_i}^{\mathrm{LB}})$$

(60)

$${\delta }_q=\{\begin{array}{c}{[2u+(1-2u){(1-\delta )}^{{\eta }_m+1}]}_{{\eta }_m+1}^1-1,u\le 0.5\hfill \\ \multicolumn{1}{c}{1-{[2(1-u)+2(u-0.5){(1-\delta )}^{{\eta }_m+1}]}_{{\eta }_m+1}^1,u>0.5}\end{array}$$

(61)

$$\delta =min(x_i-{x_i}^{\mathrm{LB}},{x_i}^{\mathrm{UB}}-x_i)/({x_i}^{\mathrm{UB}}-{x_i}^{\mathrm{LB}}),u\in [0,1]$$

(62)

where, $x_i^{\mathrm{UB}}$ is the upper limit of the $i$ th gene of the parent individual, $x_i^{\mathrm{LB}}$ is the lower limit of the $i$ th gene of the parent individual, ${\delta }_q$ a is the offset calculated from the polynomial distribution, $u$ is a random number between (0, 1), ${\eta }_m$ is the variation distribution index, which is customized by the user.

After optimization, the iterative process of the optimization objective shown in Figure 14 and the Pareto solution set of the frequency characteristic index in Figure 15 are obtained. The abscissa in Figure 14 refers to the number of swept-sine input simulation. It is obvious that through the optimization process, the vehicle has undergone 240 times to obtain the optimal optimization scheme.

Figure 14.

Iterative diagram of the optimization objective: (a) iterative graph of yaw rate gain at 0.5 Hz, (b) iterative graph of roll gain at 0.5 Hz, (c) iterative graph of resonant frequency at 0.5 Hz, (d) iterative graph of resonant frequency at 0.5 Hz, and (e) iterative graph of delay time for lateral acceleration at 0.5 Hz.

Figure 15.

Pareto solution set of frequency characteristic index: (a) pareto solution set of yaw rate gain relative to roll angle gain and (b) pareto solution of the delay time of lateral acceleration relative to the delay time of yaw rate.

Generally, in order to make the vehicle achieve better frequency characteristics, the optimization objective should show the following changes during optimization: That is, the smaller the gain of the yaw rate within a certain range, the more favorable it is to improve the handling stability; The smaller the roll angle gain within a certain range, the better the roll stability; The smaller the delay time of lateral acceleration and the delay time of yaw rate within a certain range, the faster the vehicle reaction speed and the better the handling stability; The larger the resonance frequency within a certain range, the more stable the whole vehicle is. In addition, each optimization objective should gradually converge toward the optimization direction and finally reach a stable state.

Figure 14 is the optimization iteration diagram of the yaw rate gain at 0.5 Hz, the roll angle gain at 0.5 Hz, the resonance frequency at 0.5 Hz, the delay time of lateral acceleration at 0.5 Hz, and the delay time of the yaw rate at 0.5 Hz. In Figure 14, with the advancement of the optimization process, the value of yaw rate gain at 0.5 Hz, the value of roll angle gain at 0.5 Hz, the delay time of lateral acceleration at 0.5 Hz and the delay time of yaw rate at 0.5 Hz all show a decreasing trend, and gradually converge in the decreasing direction. However, the resonance frequency at 0.5 Hz generally shows a gradually increasing trend with the advancement of the optimization process, and gradually converges to the increasing direction. It can be seen that each optimization objective conforms to the requirements of the changing trend that the above optimization objectives should meet in the optimization process, which is favorable to improve the handling stability.

After optimization, the corresponding Pareto solution set of the optimization objective is also obtained. Figure 15(a) and (b) show the Pareto solution set of the yaw rate gain at 0.5 Hz and the roll angle gain at 0.5 Hz, the Pareto solution set of lateral acceleration delay time at 0.5 Hz relative to yaw rate delay time at 0.5 Hz respectively.

Comparative analysis before and after optimization

After optimizing the frequency response characteristics, so as to verify whether the optimization results are reliable and effective, the optimal Pareto solution is selected from the Pareto solution set and brought into the established model for simulation and compared with the original car. The following mainly analyzes the transient and steady-state characteristics. Among then, the values of the vehicle transient characteristics before and after optimization are in Table 4.

Table 4.

Comparison between before and after optimization of vehicle transient frequency characteristics.

	DV1 (N/m)	DV3 (N/m)	Bushing installation angle (°)	f 1 (s−1)	f 2 (Hz)	f 3 (s)	f 4 (s)	f 5 (s−1)
Original vehicle	2.5 $0\times$ 106	3.0 $0\times$ 105	−24.50	0.34	0.53	−0.23	0.13	6.29
Optimized vehicle	4.10 × 105	1.5 $0\times 10$ 4	−27.70	0.33	0.58	−0.20	0.12	6.27
Change percent	−83.60	−95.00	+13.06	−2.94	+9.43	−13.04	−7.69	−3.17

Analysis of transient characteristics

Figure 16 shows the yaw rate gain before and after optimization. In Figure 16, when the frequency is around 0.5 Hz at low frequency, the optimized yaw rate gain is smaller than that of the original vehicle, and a smaller yaw rate gain within a certain range is beneficial to the handling stability, which meets the changing trend of the optimization target requirements.

Figure 16.

Yaw rate gain.

Figure 17 is the roll angle gain before and after optimization. In Figure 17, when the frequency is near the low frequency of 0.5 Hz, the optimized roll angle gain is slightly reduced compared with the original car, and the amount of reduction is small, but it still meets the changing trend of the optimization target requirements.

Figure 17.

Roll angle gain.

Figures 18 and 19 show the delay time of lateral acceleration and the delay time of yaw rate before and after optimization. When the low frequency is around 0.5 Hz, the absolute value of the delay time of the two is reduced compared with the original car before and after optimization. The smaller the delay time within a certain range, the faster the reaction speed of the vehicle, the more beneficial to the handling stability.

Figure 18.

Delay time for lateral acceleration.

Figure 19.

Delay time of yaw rate.

Steady-state characteristic analysis

Steady-state characteristic analysis adopts the constant radius cornering simulation. In Figure 20, The optimized steering wheel angle increases with the increase of lateral acceleration compared with the original car, indicating that the car needs a larger steering wheel angle to maintain the turning radius, and the understeer of the car increases, which is conducive to improving the handling stability.

Figure 20.

Variation curve of steering wheel angle relative to lateral acceleration.

Figure 21 shows that the optimized roll angle remains basically unchanged as the lateral acceleration increases compared to the original car. Since it is a multi-objective optimization, it is normal for individual objectives to have no substantial improvement.

Figure 21.

Variation curve of roll angle relative to lateral acceleration.

Conclusion

In this paper, the influence mechanism of the bushing installation angle and bushing stiffness of the torsion beam suspension on the C characteristics and the performance of the whole vehicle is analyzed. The frequency characteristic index is derived by view of the equivalent handling stability model, and the influence mechanism of the bushing installation angle and bushing stiffness on the frequency characteristic is revealed. Finally, the multi-variable and multi-objective optimization is executed for the bushing stiffness and bushing installation angle, which are more sensitive to the frequency characteristics. The specific conclusions are as follows:

1. The relationship between the bushing installation angle and bushing stiffness and the C characteristics of the suspension is analyzed: Under the condition of lateral force, the bushing of the suspension is deformed due to the external force, which leads to the axle steering of the suspension, which further affects the performance of the suspension. In addition, the bushing installation angle and bushing stiffness will affect the roll stiffness and force of the suspension, which will affect the suspension C characteristics, which in turn will affect the steady-state and transient characteristics of the vehicle.

2. The relationship between the bushing installation angle and the steady-state characteristic index in the range of [−90° to 90°] is discussed, that is: when the bushing installation angle is around −22.5°, the understeer characteristic of vehicle is the best, and the steady-state performance is optimal. The relationship between the change of the bushing installation angle and bushing stiffness and the change of the frequency characteristic index when the low frequency is around 0.5 Hz is obtained, namely: when the low frequency is around 0.5 Hz, the change of the bushing installation angle will have a significant impact on the yaw rate gain, the delay time of lateral acceleration and the delay time of the yaw rate, while the effect on the roll angle gain is small. However, the change of the bushing stiffness has little effect on the various indicators of the frequency characteristics of the vehicle.

3. The multi-objective optimization is executed on the frequency response characteristic index. After optimization, the optimal Pareto solution is screened out for comparative analysis with the original car. The results show that: When the radial stiffness of the bushing at the connection between the trailing arm and the body is reduced by 83.6% in the x-direction, the axial stiffness in the z-direction is reduced by 95%, and the bushing installation angle is −27.7°, it is concluded that the transient characteristics of the vehicle have been substantially improved. For the steady-state characteristics of the vehicle, the understeer of the vehicle have increased, and the roll angle is basically unchanged compared with the original car.