Mechanism of time-delay feedback control of suspension damping with an annular vibration-absorbing structure

2.1. The mechanical model

In order to research the effect of the new suspension with the annular vibration-absorbing structure on the body vertical vibration better intuitively, the structure schematic shown in Fig. 1 is constructed, in which the annular vibration-absorbing mass block is connected to the wheel hub and the swing arm support bracket through the spring, damping, and mass block support bracket. For the problem to be researched in this paper, the structure is simplified to a 3-DOF system vibration model including the body, the annular vibration-absorbing structure, and the wheel, as shown in Fig. 2(a). The tire displacement time-delay feedback control is introduced to this model to form the main research object of this paper, as shown in Fig. 2(b), and is analyzed in comparison with the conventional passive suspension, as shown in Fig. 2(c). The physical significance and values of the parameters in the model are shown in Table 1, and due to the small damping of the tires, all the researches in this paper do not consider the tire damping.

Fig. 1A new suspension model with an annular vibration-absorbing structure

Separate the system of degrees of freedom and establish the differential equations of motion of the system according to Newton’s second law:

where, $u=g{x}_t(t-\tau )$ denotes the tire displacement time-delay feedback control force, $g$ is the time-delay feedback gain coefficient, and $\tau$ is the time-delay amount.

According to the vehicle suspension performance index and passenger ride comfort performance index, for this purpose, selected body displacement $x_s$, body velocity ${\dot{x}}_s$, annular vibration-absorbing structure displacement $x_h$, annular vibration-absorbing structure velocity ${\dot{x}}_h$, tire displacement $x_t$, tire velocity ${\dot{x}}_t$ as state space variables, the differential equation of motion of the system Eq. (1) can be written in the form of the following state space equation:

where, the state variable $\mathbf{X}\left(t\right)=[x_s,{\dot{x}}_s,x_h,{\dot{x}}_h,x_t,{\dot{x}}_t{]}^T$, $\dot{\mathbf{X}}\left(t\right)=[{\dot{x}}_s,{\ddot{x}}_s,{\dot{x}}_h,{\ddot{x}}_h,{\dot{x}}_t,{\ddot{x}}_t{]}^T$; the input variables $\mathbf{U}\left(t\right)=\left[x_t\right(t-\tau \left)\right]$, $\mathbf{h}\left(t\right)=\left[x_g\right]$; the output variable $\mathbf{Y}\left(t\right)=[{\ddot{x}}_s,x_s-x_t,x_t-x_g{]}^T$; and the coefficient matrix of the state-space equation is respectively:

Fig. 21/4 mathematical model

a) PVASS

b) TDVASS

c) CS

2.2. Suspension inherent characteristics analysis

In order to analyze the vibration characteristics of the new suspension system controlled by tire displacement time-delay feedback in the frequency domain, the Fourier transform of Eq. (1) is performed separately to obtain:

which, $A_{11}=-m_s{\omega }^2+k_s+i\omega c_s$, $A_{12}=-k_s-i\omega c_s$,$A_{13}=0$,$A_{21}=-k_s-i\omega c_s$, $A_{22}=-m_h{\omega }^2+k_s+i\omega c_s+k_h+i\omega c_h$, $A_{23}=-k_h-i\omega c_h+g{e}^{-i\omega \tau }$, $A_{31}=0$, $A_{32}=-k_h-i\omega c_h$, $A_{33}=-m_t{\omega }^2+k_t+k_h+i\omega c_h-g{e}^{-i\omega \tau }$.

Thus, the frequency response function of the body droop displacement for the pavement excitation displacement input and the amplitude function of the body droop acceleration for the pavement excitation displacement input can be obtained as:

2.3. Annular vibration-absorbing structure parameters determination

The design of new suspension parameters with annular vibration-absorbing structure includes $m_h$, $k_h$ and $c_h$. The three parameters have a direct impact on the damping performance of the body. As for the design requirements of the suspension system, one is to achieve compact structure and occupy small space size; the other is to ensure the reliability of all force and moment transmission between the body and the wheels, and to ensure sufficient strength while guaranteeing the mass of components to be small. Therefore, the mass $m_h$ of the annular vibration-absorbing structure is chosen to be 0.1 times the wheel mass $m_t$, i.e., $m_h=0.1m_t$ [15]. The stiffness coefficient $k_h$ and damping coefficient $c_h$ of the annular vibration-absorbing structure were obtained using a particle swarm optimization algorithm [16]. In order to prevent the optimized parameters from not conforming to the actual characteristics of the annular vibration-absorbing structure, the optimization ranges of $k_h$ and $c_h$ were set to $(0,k_s)$ and $(0,c_s)$, respectively.

2.3.2. Particle swarm optimization algorithm (PSO)

PSO is a classical swarm intelligence algorithm, which mainly uses the sharing of information of individuals in the group, so that the motion of the whole group can produce the evolution process from disorder to order in the problem solution space, and thus obtain the optimal solution of the problem, so the PSO is used to optimize the solution of the stiffness coefficient $k_s$ and damping coefficient $c_s$ of the annular vibration-absorbing structure. The specific steps of the particle swarm optimization algorithm are shown below,

Step 1: Initialize the parameters of the optimization algorithm and generate $m$ particles randomly within the parameter optimization range, giving each particle a random initial position $X_{id}$ and an initial velocity $V_{id}$.

Step 2: Design of the fitness function, and in order to optimize the damping performance of the system, the fitness function is therefore designed as shown in Eq. (6).

Step 3: Calculate the fitness value of the generated particles and select the parameter that minimizes the fitness value as the individual best value and the objective function that minimizes the fitness value as the global best value.

Step 4: Update the position and velocity of each particle as shown in Eq. (7):

7

$\left\{\begin{array}{l}{V_{id}}^{k+1}=w^k{V_{id}}^k+c_1{r}_1\left({P_{id}}^k-{X_{id}}^k\right)+c_2{r}_2\left({P_{gd}}^k-{X_{id}}^k\right),\\ {X_{id}}^{k+1}={X_{id}}^k+{V_{id}}^{k+1},\end{array}\right.$

where $k$ is the number of current iterations, $w$ is the inertia weight, $c_1$ and $c_2$ are learning factors, $r_1$ and $r_2$ are random values of $\left[\mathrm{0,1}\right]$, $P_{id}$ is the individual optimal particle position, $P_{gd}$ is the global optimal particle position, $V_{id}$ is the particle flight speed, and $X_{id}$ is the current position of the particle.

Step 5: Optimize the second set of parameters, get the individual best and global best according to steps 2 and 3, and compare them with the first set, and update the individual best and the global best.

Step 6: According to the above steps, the individual optimum and the global optimum are updated continuously until the set number of iterations $N$ is satisfied, and the system will output the optimized optimum parameters.

The flow chart of the algorithm steps is shown in Fig. 3.

Fig. 3Flow chart of algorithm steps

By the above particle swarm optimization algorithm, the stiffness and damping coefficients of the annular vibration-absorbing structure are optimized, and the optimal parameters of the annular vibration absorbing structure are obtained as $k_h=$ 6444 N/m, $c_h=$348 N⋅s/m. Thus, the vehicle model parameters are shown in Table 1.

Table 1Vehicle model parameters meaning and values

Vehicle parameters	Physical meaning	Values	Units
$m_s$	Body mass	345	kg
$m_h$	Annular vibration-absorbing structure mass	4.05	kg
$m_t$	Wheel mass	40.5	kg
$k_s/k_t$	Vehicle suspension, tire stiffness coefficient	16000/190000	N/m
$c_s$	Vehicle suspension damping coefficient	1500	N⋅s/m
$k_h$	Annular vibration-absorbing structure stiffness coefficient	6444	N/m
$c_h$	Annular vibration-absorbing structure damping coefficient	348	N⋅s/m
$x_s/x_h/x_t$	Body, annular vibration-absorbing structure, wheel vertical displacement	–	m
$x_g$	Road surface vertical input displacement	–	m

3.1. Time-delay independent stability

Whether the time-delay feedback control can ensure that the system is stable depends on the value of the feedback gain $g$ and the time-delay amount $\tau$. When any value of the time-delay amount is taken, the stability of the system is not changed with the time-delay amount. The stability of the system does not change with the time-delay amount when the time-delay amount takes any value, but only changes with the feedback gain, which is called the full time-delay stable state, i.e., time-delay independent stability. For the time-delay independent stability of the system, it is necessary to satisfy the following two conditions,

(1) When $\tau =0$, the characteristic Eq. (8) of the system satisfies Routh-Hurwitz stability.

(2) When $\tau >0$, the characteristic Eq. (8) of the system has no positive real roots ${\omega }_0$.

Condition (1): when $\tau =0$, the characteristic equation of the system polynomial Eq. (8) is expressed as:

It is concluded from the Routh-Hurwitz stability criterion that for the system to be in a stable state, firstly, all coefficients of the characteristic Eq. (10) of the system are required to be greater than zero; secondly, all terms in the first column of the Routh-Hurwitz table are required to be positive, and the Routh-Hurwitz discriminant table constructed from the polynomial coefficients can be obtained according to Eq. (10), as shown in Fig. 4, which, $d_1=a_4+b_4$, $d_2=a_3+b_3$, $d_3=a_2+b_2$, $e_1=\frac{a_5{d}_1-a_6{d}_2}{a_5}$, $e_2=\frac{a_5{d}_3-a_6{a}_1}{a_5}$， $f_1=\frac{e_1{d}_2-a_5{e}_2}{e_1}$, $f_2=\frac{e_1{a}_1-a_5{a}_0}{e_1}$, $l_1=\frac{f_1{e}_2-e_1{f}_2}{f_1}$, $p_1=\frac{l_1{f}_2-f_1{a}_0}{l_1}$.

Fig. 4Routh-Hurwitz discriminant table

Substituting the parameters of Table 1 into Fig. 4, the stability conditions $a_6$, $a_5$, $d_1$, $d_2$, $d_3$, $a_1$, and $a_0$ can be found to be constant greater than zero. Other stability regions drawn according to the Routh-Hurwitz stability conditions are shown in Fig. 5.

Fig. 5Routh-Hurwitz stability conditions

From Fig. 5, it is clear that when $g\in (-\infty ,171870.8]$, the polynomial coefficients of the characteristic equation of the system and the first column coefficients in the Routh-Hurwitz discriminant table are greater than zero, which satisfies the Routh-Hurwitz stability criterion and indicates that the no time-delay system is stable in this interval.

Condition (2): When $\tau >0$, it is known from stability theory that the asymptotic stability of the time-delay system needs to satisfy that the characteristic roots of the polynomial of the characteristic equation don't contain positive real roots. When the characteristic roots contain pure imaginary numbers, the system is at the boundary of instability and there exists stability switching as well as critical time-delay stability interval. Assuming the existence of imaginary roots $s=i{\omega }_0$ in the characteristic equation, the characteristic polynomial of the system Eq. (8) can be expressed as $D(s,\tau )=D(i{\omega }_0,\tau )=0$ and simplified by using Euler’s formula $e^{-i{\omega }_0\tau }=\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_0\tau \right)-\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_0\tau \right)$ for $D(i{\omega }_0,\tau )=0$, which can be obtained:

where:

Separating the real and imaginary parts shows that if the polynomial Eq. (11) is zero, it must satisfy that both the real and imaginary parts are zero, i.e.:

Further derivation leads to:

With $\mathrm{s}\mathrm{i}{\mathrm{n}}^2\left({\omega }_0\tau \right)+\mathrm{c}\mathrm{o}{\mathrm{s}}^2\left({\omega }_0\tau \right)=1\text{,}$ the trigonometric term is eliminated to obtain a polynomial equation for ${\omega }_0$, i.e.:

Such that $F\left({\omega }_0\right)=P_R({\omega }_0{)}^2+P_I({\omega }_0{)}^2-\left[Q_R\right({\omega }_0{)}^2+Q_I\left({\omega }_0{)}^2\right]$, the characteristic root discriminant polynomial of the system can be obtained, i.e.:

where:

After a series of derivations, $D(i{\omega }_0,\tau )=0$ can finally be reduced to $F\left({\omega }_0\right)=0$. By solving the value and number of positive real roots of the polynomial equation $F\left({\omega }_0\right)=0$, the range of its feedback gain $g$ can be determined when the time-delay independent stability is achieved. Fig. 6 shows the time-delay independent stability region of the new suspension under the tire displacement time-delay feedback control.

According to Fig. 6, when $g\in [-\mathrm{24427.7,24427.7}]$, the characteristic equation of the system has no positive real roots, and the range of values of this gain coefficient $g$ is in the stability interval of condition Eq. (1), which means that the system is time-delay independent stable in this interval range. When the value range of this feedback gain is exceeded, the polynomial equation has positive real roots, which means that the stability of the system will switch a certain number of times and needs further discussion and analysis.

Fig. 6Time-delay independent stability region

3.2. Stability switching

When the time-delay feedback control parameter $\tau$ is outside the region of time-delay independent stability, the stability of the system is not fixed, and the sign of the real part of the characteristic root changes with time-delay, which means that the stability of the system changes with time-delay, and this situation is called stability switching. Therefore, finding the critical time-delay is crucial to analyze whether the system is stable or not.

Assuming that the polynomial Eq. (15) has positive real roots ${\omega }_0$, all the critical time-delays corresponding to the positive real roots ${\omega }_0$ can be solved by the following equation:

where, $\theta ={\omega }_0\tau$, $\theta \in \left[\mathrm{0,2}\pi \right]$.

After all the critical time-delays are determined by Eq. (16), the change of the real part of the characteristic root corresponding to each critical time-delay is then judged to further determine whether the system is stable or not. The characteristic equation polynomial Eq. (8) can be regarded as a function of $s$ with respect to $\tau$. Its derivative can be obtained from the stability switching discriminant formula:

Then, the change of the real part of the characteristic root can be determined according to the positive or negative sign of $S$. In the case of $S=+1$, as the time-delay $\tau$ increases, it means that the real part of the characteristic root crosses the imaginary axis from left to right in the complex plane, thus switching the system from stable to unstable state; in the case of $S=-1$, it means that the real part of the characteristic root crosses the imaginary axis from right to left in the complex plane, and the system switches from unstable to stable state. The following is a concrete example to further illustrate the switching principle of stability.

When $g=$ 50000 Ν/m, the number of positive real roots of the characteristic equation is two, ${\omega }_1^{*}=$ 62.3487, ${\omega }_2^{*}=$ 75.1147, according to Eq. (15), as shown in Fig. 6. Substituting these two positive real roots into Eq. (17), we can get $\dot{F}\left({\omega }_1^{*}\right)<0$ and $\dot{F}\left({\omega }_2^{*}\right)>0$, and thus the critical time-delay corresponding to these two positive real roots can be found according to Eq. (16) as:

The obtained critical time-delays are sorted to obtain:

Since $\dot{F}\left({\omega }_1\right)<0$, i.e., $S=-1$, when the time-delay $\tau$ crosses the time-delay ${\tau }_{1,j}^{*}$$(j=\mathrm{1,2},3,\cdots )$, the real part of the characteristic root changes from positive to negative, and the system switches from unstable to stable state; similarly, when $\dot{F}\left({\omega }_2\right)>0$, i.e., $S=+1$, the time-delay $\tau$ crosses ${\tau }_{2,j}^{*}(j=\mathrm{1,2},3,\cdots )$, and the real part of the characteristic root changes from negative to positive, then the system switches from stable to unstable state. From the above analysis, it can be seen that when $\tau \in [0,{\tau }_{14}^{*})$, the system is asymptotically stable; when $\tau \in [{\tau }_{14}^{*},{\tau }_{15}^{*}]$, the system is unstable. From Eq. (19), it can be seen that the critical time-delay ${\tau }_{14}^{*}$ is followed by ${\tau }_{15}^{*}$, no switching occurs, for two consecutive unstable critical time-delays, then it shows that after ${\tau }_{15}^{*}$, the system is always in an unstable state. The final time-delay interval that puts the system in a stable state when $g=$ 50000 Ν/m is obtained is $[0,{\tau }_{11}^{*}]\bigcap [{\tau }_{21}^{*},{\tau }_{12}^{*}]\bigcap [{\tau }_{22}^{*},{\tau }_{13}^{*}]\bigcap [{\tau }_{23}^{*},{\tau }_{14}^{*}]$.

When $g$ = –50000 Ν/m, the positive real roots of the characteristic equation are obtained as ${\omega }_3^{*}=$ 62.3487, ${\omega }_4^{*}=$ 75.1147, respectively, as shown in Fig. 6. Further solving for $\dot{F}\left({\omega }_3\right)>0$ and $\dot{F}\left({\omega }_4\right)<0$, then the critical time-delays corresponding to these two positive real roots are obtained as:

The obtained critical time-delays are sorted to obtain:

It can be concluded that no switching occurs after the critical time-delay ${\tau }_{34}^{*}$, which directly becomes ${\tau }_{35}^{*}$, for two consecutive unstable critical time-delays, then it is derived that after ${\tau }_{35}^{*}$ the system are in an unstable state. Therefore, when $g$ = –50000 N/m, the critical time-delay interval that keeps the system stable is $[0,{\tau }_{31}^{*}]\bigcap [{\tau }_{41}^{*},{\tau }_{32}^{*}]\bigcap [{\tau }_{42}^{*},{\tau }_{33}^{*}]\bigcap [{\tau }_{43}^{*},{\tau }_{34}^{*}]$.

Fig. 7 shows the critical time-delay curves of $(g,\tau )$ under tire displacement time-delay feedback control. The intersection of the purple and light blue curves on the right side of Fig. 7 is before the critical time-delay point, and the system is stable, while the intersection of the green and purple curves is after the critical time-delay point, and the system is unstable, and the system state has changed, i.e., the stability switching occurs between $[{\tau }_{14}^{*},{\tau }_{15}^{*}]$ shown in the figure. Similarly, the intersection of the purple and light blue curves on the left side of Fig. 7 is on the outside of the critical time-delay point, and the system is stable, while the intersection of the green and purple curves is on the inside of the critical time-delay point, and the system is unstable, which indicates that the system has switched stability between $[{\tau }_{34}^{*},{\tau }_{35}^{*}]$.

Fig. 7(g,τ) critical time-delay curves graph

In the case that the gain coefficient is determined, the polynomial solution of the characteristic equation is used to get the positive real roots. If there are no positive real roots in the equation, the system is stable regardless of the value of $\tau$; if there are positive real roots in the equation, the system will undergo a certain stability transformation and finally become unstable. According to the change of the number of positive real roots, we can get the full time-delay stability region and stability switching region of the system, and then judge whether the system corresponding to the critical time-delay point is in the stable state according to the stability switching discriminant formula, and finally we can get the critical stability region in the $(g,\tau )$ plane under the tire displacement time-delay feedback control, as shown in Fig. 8.

Fig. 8Critical time-delay stability region

3.3. Stability verification

In order to verify the correction of the region in Fig. 8, three points are randomly selected for verification in each stability interval and instability interval of the system, where points $A_1^{*}$ ($g=\text{40000}$, $\tau =0.\text{03}$), $A_4^{*}$($g=$ –40000, $\tau =$ 0.25), and $A_5^{*}$ ($g=$ 40000, $\tau =$ 0.3) correspond to the stability interval of the system; points $A_2^{*}$ ($g=$–40000, $\tau =$ 0.11), $A_3^{*}$ ($g=$ 40000, $\tau =$ 0.15), and $A_6^{*}$ ($g=$ –40000, $\tau =$ 0.4) correspond to the instability interval of the system. The coordinates of each point are substituted into Eq. (1), and the initial condition is taken as $X\left(t\right)={\left[0.0100000\right]}^T$. The numerical verification results of the BA at each point are shown in Fig. 9.

Fig. 9Numerical verification results

a)$A_1^{*}$

b)$A_2^{*}$

c)$A_3^{*}$

d)$A_4^{*}$

e)$A_5^{*}$

f)$A_6^{*}$

From the numerical verification results in Fig. 9, it can be seen that when the initial excitation of the system is given, the BA vibration amplitude gradually converges to zero, which indicates that the system is asymptotically stable; the BA vibration amplitude gradually diverges to infinity, which indicates that the system is unstable. The verification results in Fig. 9 are consistent with the corresponding system stability cases in Fig. 8, and the numerical simulation results verify the correctness of the critical stability time-delay stability region map.

4.1. Method of solving the time-delay equation

The solution is difficult because the characteristic equation contains a transcendental function term. Therefore, in this paper, the fine integration method of time-delay differential equations is used to solve the numerical solution of each response quantity of the 1/4 vehicle model vibration.

At first, treating the time-delay feedback control term $U\left(t\right)$ in Eq. (2) as the excitation term, the solution of the equation can be written as:

Next, the numerical discretization of Eq. (22) is performed by assuming that the simulation time is $T$. Taking the time step $\mathrm{\Delta }t=t_{k+1}-t_k(k\in [0,n\left]\right)$, then $n=T/\mathrm{\Delta }t$, the solution of the equation is discretized into the following stepwise integral equation in recursive form:

The above equation can be further derived as:

where, $m=\tau /\mathrm{\Delta }t$, is expressed as the number of nodal steps that the time-delay quantity $\tau$ is divided into during the solution process. Further solving the solution of Eq. (24), assuming the initial condition $\mathbf{X}\left(t\right)={\left[000000\right]}^T$ of the equation, the system vibration response of each mass block (body, annular vibration-absorbing structure, wheels) at each time point $k$ is found as:

where, $x_{sk}$, ${\dot{x}}_{sk}$ represents the vibration displacement response and vibration velocity response of the body mass $m_s$ at time $t_k$; $x_{hk}$, ${\dot{x}}_{hk}$ represents the vibration displacement response and vibration velocity response of the annular vibration-absorbing structure mass $m_h$ at time $t_k$; $x_{tk}$, ${\dot{x}}_{tk}$ represents the vibration displacement response and vibration velocity response of the tire mass $m_t$ at time $t_k$, respectively.

Finally, by substituting Eq. (25) into Eq. (1), the vibration acceleration of the body, the annular vibration-absorbing structure and the tire in the time domain are obtained respectively:

4.2. Optimization of feedback control parameters

For the optimization of the time-delay feedback control parameters the PSO described above is still used, and its objective function is shown in Eq. (6). With 100 randomly generated particles and 200 iterations, the fine integration method for solving the time-delay and the pavement displacement are introduced into the PSO to optimize the objective function in the time domain and solve the appropriate time-delay feedback control parameters for different pavements. The simple harmonic excitation and random excitation are selected for the pavement displacement, and the variation curves of the adaptation degree values with the number of iterations for different pavement excitations are shown in Fig. 10. From Fig. 10, it can be seen that the optimal time-delay feedback control parameters are basically obtained at the later stage of the iteration number, and the optimal time-delay feedback control parameters under simple harmonic excitation and random excitation are finally obtained as $g_1=$ 5386 N/m, ${\tau }_1=$ 2×10^-4 s and $g_2=$ 2528 N/m, ${\tau }_2=$1.28×10^-4 s, respectively. The optimized two sets of parameters are in the stability region in Fig. 8, which can make the system in a stable state, indicating that the optimized time-delay feedback control parameters are reasonable.

Fig. 10Adaptation function change curves

a) Simple harmonic excitation

b) Random excitation

5.1. Simulation analysis under simple harmonic excitation

The above optimized tire displacement time-delay feedback control parameters $g_1$ and ${\tau }_1$ are substituted into Eq. (1), and the simple harmonic excitation $x_g=0.05\mathrm{s}\mathrm{i}\mathrm{n}\left(5t\right)$ is simulated as the tire road input. The effects of the CS, the PVASS, the TDVASS on the system damping performance in the time domain are analyzed. The simulation curves of each response of the system were obtained as shown in Fig. 11, and the RMS values of each response quantity of the system were shown in Table 2.

As can be seen from Fig. 11, under the optimal feedback control parameters, the BA, SDD and TDD of the TDVASS are reduced compared with the PVASS, and the amplitude is lower in the full frequency domain than that of the PVASS; while compared with the CS, the BA and TDD are significantly reduced, and the SDD compared with the CS, the BA and TDD are significantly reduced, while the SDD is slightly deteriorated, but the increase is still within the design requirement (the suspension dynamic travel of ±100 mm is selected in this paper [17]), which doesn’t exceed the limit travel of the suspension dynamic deflection, i.e., it doesn’t increase the probability of hitting the limit block.

Fig. 11Simulation comparison under simple harmonic excitation

a) Body droop acceleration

b) Suspension dynamic deflection

c) Tire dynamic displacement

In Table 2, it can be obtained that compared with the CS, the RMS values of BA and TDD of the TDVASS are reduced from 1.6081, 3.1109 to 0.3068, 0.6251, respectively, and the optimization ratios are 80.92 % and 79.91 %; compared with the PVASS, the RMS values of BA, SDD, TDD of the TDVASS are reduced from 0.9715, 0.0693, 1.6606 to 0.3068, 0.0388, 0.6251, respectively, and the optimized ratios are 68.42 %, 44.01 %, 62.36 %, respectively. The above analysis shows that the tire displacement time-delay damping control method can effectively alleviate the impact of road excitation on the body vibration, so that its comfort can be significantly improved. At the same time, the grounding of the tires is also improved, which guarantees the safety of the vehicle.

5.2. Simulation analysis under random excitation

The simple harmonic excitation is only for a certain fixed frequency pavement, which is more desirable for a road, while the random excitation is more in line with the actual road conditions. According to the filtered white noise method described in ISO 8608:1995 as shown in Eq. (27), a standard C-level random road surface [18] is constructed as shown in Fig. 12(a), which is used to further verify the damping effect of the TDVASS on the body. To illustrate the effectiveness of the time domain random pavement excitation model generated based on the filtered white noise method, a comparative power spectral density plot of the C-level pavement is established, as shown in Fig. 12(b), which shows that the simulated power spectral density of the pavement unevenness is in good agreement with the standard power spectral density, indicating the correctness of the generated pavement. This pavement can be used for simulation analysis of the model:

where, $f_0$ is the under cut-off frequency, generally take $f_0=0.1$; $G_q\left(n_0\right)$ is the A-H level pavement unevenness coefficient; $n_0$ is the reference spatial frequency, $n_0=$ 0.1 m^-1; $u$ is the car driving speed, this paper takes $u=$ 20 m/s; $\omega \left(t\right)$ is the Gaussian white noise signal with zero mean value.

Fig. 12C-level pavement

a) C-level pavement unevenness

b) Power spectrum density comparison chart

Using Fig. 12 (a) as the random excitation of the pavement and substituting the optimal feedback control parameters $g_2$ and ${\tau }_2$ into Eq. (1), the simulation curves of each response quantity of the system are shown in Fig. 13, and the root mean square values of each response quantity of the system are shown in Table 3.

Fig. 13Simulation comparison under random excitation

a) Body droop acceleration

b) Suspension dynamic deflection

c) Tire dynamic displacement

According to Fig. 13 and Table 3, it can be concluded that the TDVASS has greater advantages compared to the CS and the PVASS. Compared with the CS, the RMS values of BA and TDD are reduced from 0.3184, 0.6399 to 0.1014, 0.5289, respectively, with the optimized ratio of 81.72 % and 19.41 %, and the SDD is slightly deteriorated, but its vibration amplitude is still within the dynamic travel of the suspension and does not exceed its limit travel. Compared with the PVASS, the RMS values of BA, SDD, and TDD were reduced from 0.1512, 0.0105, and 0.5607 to 0.1014, 0.0073, and 0.5289, respectively, with optimization ratios of 61.51 %, 28.57 %, and 8.03 %. As a result, it can be derived that the time-delay feedback control can also effectively reduce the body vibration and improve the comfort of the vehicle and the grounding of the tires under the random excitation with the selection of suitable feedback parameters.

5.3. Amplitude and frequency characteristics analysis

Through frequency domain analysis of CS, PVASS and TDVASS, the response quantities of each vibration characteristic of the suspension system at different frequencies in the studied frequency band can be obtained, and the effect of different time-delay feedback control on vehicle damping can be further investigated. With Eq. (4) and (5), the amplitude frequency characteristic functions of the system can be obtained. The optimal time-delay feedback control parameters $g_2$ and ${\tau }_2$ obtained under random road excitation are substituted into the system to obtain the amplitude frequency characteristic curves of BD and BA under tire displacement time-delay feedback, as shown in Fig. 14.

Fig. 14Amplitude frequency characteristic curves

a) Body droop displacement

b) Body droop acceleration

It can be seen that due to the addition of the under-spring mass annular vibration-absorbing structure to form the PVASS on the basis of the CS , and the introduction of the tire displacement time-delay feedback control to form the PVASS, which results in the low frequency resonance frequencies of the PVASS and the TDVASS being smaller than those of the CS, and the high frequency resonance frequencies are almost unchanged. It can also be concluded that compared with the CS, the TDVASS only has a certain degree of deterioration in the low frequency resonance region, and the optimization effects in the rest of the frequency bands are substantially improved. With PVASS, the amplitude frequency characteristic curves of BD and BA in TDVASS show substantial improvement in both the low frequency region and the medium frequency region, and the curves almost overlap in the high frequency region with no obvious optimization effect. The above analysis shows that TDVASS has obvious effect on reducing body vibration in the low and medium frequency bands, and the optimization effect is not obvious in the high frequency band.