Speed Oscillations of a Vehicle Rolling on a Wavy Road

Abstract

Every driver knows that his car is slowing down or accelerating when driving up or down, respectively. The same happens on uneven roads with plastic wave deformations, e.g., in front of traffic lights or on nonpaved desert roads. This paper investigates the resulting travel speed oscillations of a quarter car model rolling in contact on a sinusoidal and stochastic road surface. The nonlinear equations of motion of the vehicle road system leads to ill-conditioned differential-algebraic equations. They are solved introducing polar coordinates into the sinusoidal road model. Numerical simulations show the Sommerfeld effect, in which the vehicle becomes stuck before the resonance speed, exhibiting limit cycles of oscillating acceleration and speed, which bifurcate from one-periodic limit cycle to one that is double periodic. Analytical approximations are derived by means of nonlinear Fourier expansions. Extensions to more realistic road models by means of noise perturbation show limit flows as bundles of nonperiodic trajectories with periodic side limits. Vehicles with higher degrees of freedom become stuck before the first speed resonance, as well as in between further resonance speeds with strong vertical vibrations and longitudinal speed oscillations. They need more power supply in order to overcome the resonance peak. For small damping, the speeds after resonance are unstable. They migrate to lower or supercritical speeds of operation. Stability in mean is investigated.

1. Introduction to the Problem

Vertical vibrations of a vehicle driven by a constant force and rolling on a sinusoidal road surface are coupled with its horizontal travel motion, affecting the vehicle speed which fluctuates around a mean value. The coupling between both planar motions is caused by the permanent direction change of the contact force to ground along the contour of the road profile. This paper explains the nonlinear model of this dynamic problem applying averaging methods to calculate stationary solutions before and after the resonance speed. Numerical integrations are applied to obtain limit cycles around the averaged solutions, plotting the fluctuating car acceleration against the true velocity. Stationary solutions are stable in mean when the slope of the driving force speed characteristic is positive. Vice versa, they are unstable for negative slopes. This leads to the so-called Sommerfeld effect [1] that for a given driving force the car becomes stuck before the resonance speed and can only pass over the resonance and the unstable velocity range after the resonance by considerably increasing the driving force [2]. First investigations of velocity jumps and turbulent speeds in nonlinear vehicle road dynamics are given by Wedig in [3,4,5,6,7] applying sinusoidal and random road models introduced by Robson et al. [8,9,10]. The first order road model in [11,12] is extended in [3,4] to a second order one which includes sinusoidal models. Blekhman and Kremer studied the same vehicle road system in [13,14] for the special case of small road excitations to calculate only the average response of driving cars (see also [15]). In [2,7], these investigations are extended to large sinusoidal road surfaces and double-periodic limit cycles in the phase plane of longitudinal acceleration and oscillating speed. In the present paper, new results for mean values and amplitudes of speed oscillations are calculated by means of Fourier expansions. Stability investigations by means of Floquet theory are proposed. Extensions to quarter car models with two degrees of freedom are made. The new speed amplitudes calculated in this paper show that the longitudinal speed oscillations of the vehicle are stable in the lower speed range before the resonance speed and in the upper higher speed range. In the middle range immediately after the resonance peak, stationary speeds are unstable and therefore physically not realizable. These stability properties correspond to the Duffing problem where vertical displacement vibrations of the vehicle possess three different amplitudes in the resonant speed range: the upper and lower displacement vibrations are stable and the middle range vibrations are unstable, as well.

2. Coupled Vertical and Longitudinal Vehicle Road Dynamics

Figure 1a shows the applied quarter car model rolling on a wavy road with vertical displacement $z$ and derivative $u$ taken along the travel way $s$. In the following, the quantities $z$ and $u$ are called road level and slope, respectively. They generate vertical car vibration displacement $y$ and velocity $\dot{y}$, which are coupled by the vehicle speed $v=\dot{s}$ and described by the two equations of motion

$$\dot{v}=\left[{\omega }_1^2\left(y-z\right)+2D{\omega }_1\left( \dot{y}-\dot{z}\right)\right]\tan \alpha +f/m ,$$

(1)

$$\ddot{y}+2D{\omega }_1(\dot{y}-\dot{z})+{\omega }_1^2(y-z)=0, \dot{y}={\omega }_{1}x$$

(2)

where s is the travel displacement of the vehicle and $x=\dot{y}/{\omega }_1$ denotes the coordinate of the vertical vibration velocity. In Equations (1) and (2), dots denote derivatives with respect to time t. The parameter ${\omega }_1^2=c/m$ determines the natural frequency ${\mathsf{\omega }}_1$ of the vertical vehicle vibrations, $2D{\omega }_1=b/m$ denotes the damping, and $f$ is the driving force which is constant or slightly decreasing with growing speed. In Figure 1b, both force characteristics are plotted in yellow-black. In the following, constant driving force is applied only. The nonlinear term in Equation (1) represents the damper and spring force multiplied by $\tan \alpha$ that takes the horizontal component of the contact force $N$ by means of $\tan \alpha =dz/ds$ One finds Equations (1) and (2) already in the literature in [13,14] and in [2,3,4,5,6,7].

In first investigations, the road level and slope are assumed to be sinusoidal as

$$z\left(s\right)=z_0\cos \left(\Omega s\right) \mathrm{and} u\left(s\right)=-z_0\sin \left(\Omega s\right),$$

(3)

where $z_0$ is the amplitude of the road surface and $\Omega$ is the angular spatial frequency determined by the wave length $L=2\pi /\Omega$ of the road. Equations (1)–(3) represent a DA equation system where the role of algebraic terms is taken by sinusoidal terms. To eliminate these terms, the road level $z$ and slope $u$ in Equation (3) are differentiated with respect to the way coordinate $s$ in order to obtain the increments $dz=-z_0\Omega \sin \left(\Omega s\right)ds$ and $du=-z_0\Omega \cos \left(\Omega s\right)ds$ that leads to the homogeneous nonlinear oscillator equations

$$\dot{z}=v\Omega u , \dot{u}=-v\Omega z ,$$

(4)

which are obtained when both increments above are divided by $dt$ and $ds/dt$ is replaced by the speed v. Furthermore, $dz/ds=-z_o\Omega \sin \left(\Omega s\right)=\Omega u$ holds so that Equation (1) reads as

$$\dot{v}=\left[{\omega }_1^2\left(y-z\right)+2D{\omega }_1\left( \dot{y}-\dot{z}\right)\right]\Omega u+f/m.$$

(5)

Equations (2), (4) and (5) describe a five-dimensional problem with five unknowns [4,7]: the horizontal travel speed $v\left(\mathrm{t}\right)$ of the vehicle, its vertical vibration by displacement $y\left(t\right)$ and velocity $\dot{y}\left(t\right)={\omega }_{1}x\left(t\right)$, and the road level $z\left(t\right)$ and slope $u\left(t\right)$. For analytical and numerical investigations, it is appropriate to introduce the dimensionless time $\tau ={\omega }_{1}t$ and the related speed $v=v\mathsf{\Omega }/{\omega }_1$, as well as the related coordinates $\left(\overline{z}, \overline{u}\right)=\left(z,u\right)z_0$ and $\left(\overline{y},\overline{x}\right)=\mathsf{\Omega }\left(y,x\right)$. Their insertion into Equations (2), (4) and (5) leads to

$$v^{\prime }+2D{(z_0\Omega )}^2{\overline{u}}^{2}v=z_0\Omega (\overline{uy}+2D\overline{ux})-{(z_0\Omega )}^2\overline{zu}+f\Omega /c,$$

(6)

$${\overline{z}}^{\prime }=v\overline{u}, {\overline{u}}^{\prime }=-v\overline{z}, \mathrm{IC}.: \overline{z}\left(0\right)=1, \overline{u}\left(0\right)=0,$$

(7)

$${\overline{y}}^{\prime }=x, {\overline{x}}^{\prime }=-\left(\overline{y}+2D\overline{x}\right)+z_o\Omega \left(\overline{z}+2Dv\overline{u}\right) . \tau ={\omega }_{1}t,$$

(8)

where prime denotes differentiation with respect to the dimensionless time $\tau ={\omega }_{1}t.$ To improve numerical integration [2,7] in Equations (6)–(8), the polar coordinates

$$\overline{z}=\overline{r} \cos \phi \mathrm{and} \overline{u}=\overline{r} \sin \phi$$

(9)

are introduced into the road Equation (7) that leads to the transformation equations

$${\overline{r}}^{\prime }\cos \phi -\overline{r}{\phi }^{\prime }\sin \phi =v\overline{r} \sin \phi ,$$

$${\overline{r}}^{\prime }\sin \phi +\overline{r}{\phi }^{\prime }\cos \phi =-v\overline{r}\cos \phi ,$$

which are solved by means of the determinant $∆={\cos }^2\phi +{\sin }^2\phi =1$ and Cramer’s rule to

$${\overline{r} }^{\prime }=0, {\phi }^{\prime }=-v.$$

Without loss of generality, the related polar radius is integrated to $\overline{r }=1.$ Note that the derivative of the polar angle is equal to the negative speed of the vehicle. According to the definition of polar coordinates, the polar angle turns counterclockwise into the mathematically positive direction. The applied oscillator, however, rotates clockwise. This is the reason why both quantities have opposite signs.

3. Speed Driving Force Characteristic of Traveling Vehicles

In order to derive the driving force speed characteristic, shown in Figure 2, the equations of motion are approximately investigated assuming that the oscillating speed of the vehicle can be averaged by $v\left(\tau \right)=v=\mathrm{const}.$ In this case, the sinusoidal solutions

$$\overline{y}\left(\tau \right)=y_c\cos \left(v\tau \right)+y_s\sin \left(v\tau \right), \overline{z}\left(\tau \right)=\cos \left(v\tau \right),$$

$$\overline{x}\left(\tau \right)=x_c\cos \left(v\tau \right)+x_s\sin \left(v\tau \right), \overline{u}\left(\tau \right)=-\sin \left(v\tau \right)$$

are applicable and inserted into Equation (8) in order to obtain the matrix equation

$$\left[\begin{array}{cccc}0& -v& -1& 0\\ v& 0& 0& -1\\ 1& 0& 2D& -v\\ 0& 1& v& 2D\end{array}\right]\left[\begin{array}{c}{y}_s\\ y_c\\ x_s\\ x_c\end{array}\right]=z_0\Omega \left[\begin{array}{c}0\\ 0\\ 2Dv\\ 1\end{array}\right]$$

applying the coefficient comparison. The first two rows of this matrix equation yield $x_c=v{y}_s$ and $y_c=-x_s/v$. Both relations are inserted into the last two rows, obtaining

$$\left(v^2-1\right) x_s+2D{v}^2{y}_s=z_o\Omega v,$$

$$-2D{x}_s+\left(v^2-1\right)y_s=z_o\Omega 2Dv.$$

This reduced equation system possesses the determinant $∆={\left(v^2-1\right)}^2+{\left(2Dv\right)}^2$. It is solved by Cramer’s rule in order to obtain first $x_s \mathrm{and} y_s$ and then $x_c \mathrm{and} y_c$, as follows:

$$x_s=z_o\Omega v\left[v^2-1-{\left(2Dv\right)}^2\right]/∆, x_c=z_o\Omega 2D{v}^4/∆,$$

(10)

$$y_s=z_o\Omega 2D{v}^3/∆, y_c=z_o\Omega [1+v^2\left(4D^2-1\right)]/∆.$$

(11)

These results coincide with the linear theory of a constant vehicle speed. They are inserted into the velocity Equation (6) to calculate the time dependent driving force

$$f\mathsf{\Omega }/c=2D{\left(z_0\Omega \right)}^{2}v{\sin }^{2}v\tau +z_0\Omega \left[\left(\overline{y}+2D\overline{x}\right)-z_0\Omega \cos \left(v\tau \right) \right]\sin \left(v\tau \right)$$

which is needed to keep the speed $v$ constant. Its mean value is calculated as

$$f\mathsf{\Omega }/c ={\left(z_0\Omega \right)}^{2}D{v}^5/[{\left(v^2-1\right)}^2+{\left(2Dv\right)}^2] , \mathrm{with} v=v\Omega /{\omega }_1.$$

(12)

The approximated force speed characteristic in Equation (12) is plotted in Figure 1b for the two road surface parameters $z_0\Omega =1.4 \mathrm{and} 0.8$, marked by red and cyan lines, respectively. The dashed line represents the asymptote $f\mathsf{\Omega }/c={\left(z_0\Omega \right)}^{2}Dv$, which is proportional to the speed, indicating that one needs a linearly growing driving force to reach higher speeds of operation. Obviously, the increasing driving force is needed to compensate the energy loss in the damper, which is growing with higher speeds, as well.

The force speed characteristic (12) represents a linear approximation obtained by averaging the oscillating driving speed. In order to check the validity and stability of the result in Equation (12), numerical simulations are performed by applying the Euler scheme to

$$v^{\prime }+2D{\left(z_0\Omega \right)}^2{\overline{u}}^{2}v=z_0\Omega \left(\overline{uy}+2D\overline{ux}\right)-{\left(z_0\Omega \right)}^2 \overline{zu}+f\mathsf{\Omega }/c,$$

(13)

$${\overline{y}}^{\prime }=x, {\overline{x}}^{\prime }=-\left(\overline{y}+2D\overline{x}\right)+z_o\Omega \left(\overline{z}+2Dv\overline{u}\right) . {\phi }^{\prime }=-v,$$

(14)

where the related level and slope of the road surface are given by $\overline{z}= \cos \phi$ and $\overline{u}= \sin \phi$ dependent on the polar angle $\phi$, the derivative of which is equal to the negative speed of the vehicle. The numerical results obtained are presented in Figure 2a by plotting the scaled and shifted acceleration of the vehicle against the true speed where squares denote initial values of acceleration and speed; triangles are mean values of acceleration and speed calculated after the initially transient period, providing a sufficiently long time for the averaging procedure. In Figure 2b, stationary limit cycles are shown for a stronger road excitation given, e.g., by $z_0\Omega =1.0$, which leads to the road amplitude $z_o\approx 3.2 \mathrm{mm}$, e.g., for the wave length $L=20 \mathrm{mm}.$ Note that phase portraits of velocity over displacement are not applicable for the travel kinematics since the horizontal displacement of the traveling vehicle is growing infinitely. Instead of phase portraits of displacement and velocity, Figure 2b shows limit cycles of velocity against acceleration. Obviously, there are two speed regions where the limit cycles are stable; the first is in the under-critical speed range before the resonance speed $v=1.$ The second is far beyond the resonance in the higher speed range of operation. In between both stable speed ranges, the slope of the speed driving force characteristic is negative. In this range, limit cycles are not stable and therefore not realizable. This instability is plausible and physically explained inside the yellow area shown in Figure 2a. Accordingly, a speed perturbation by means of $∆v<0$ into the negative speed direction on the left side of the dynamic equilibrium generates an acceleration back to the equilibrium since the applied force $∆f>0$ is larger than that one being necessary in the new perturbed situation. In this case the positive driving force difference is equal to the vertical distance between the green and red circle. However, the vehicle is braked if the speed perturbation goes into the positive speed direction on the right side of the dynamical equilibrium. In this case, the applied driving force is smaller than the one necessary to maintain the new perturbed dynamic equilibrium, marked by the right red circle. Vice versa, a speed driving force characteristic with negative slope leads to monotonous instability with the effect that speed leaves the unstable branch. In Section 5, it is shown that the negative slope condition coincides with the instability in mean, which is derived by applying the Hurwitz criterion to the variational equations of the averaged equations of motion.

Figure 2a shows three limit cycles in the stable under-critical speed range for the driving forces: $f\Omega /c=0.1, 0.2, \mathrm{and} 0.3$ marked by green, cyan, and yellow triangles, respectively. The limit cycles are obtained by plotting the scaled and shifted acceleration against the true travel speed. The selected initials are marked by colored squares. Thin black lines are transients which start in the squares and end in the thick colored lines of the one-periodic limit cycles. After this initial period, the simulation is continued a sufficiently long time in order to calculate the mean values of speed and acceleration, which are plotted marked by colored triangles. Obviously, the triangles showing speed and acceleration coincide with the initial values of the applied speed driving force characteristic. Since the accelerations are shifted upward by the applied value of the speed driving force characteristic, the acceleration mean value is vanishing when both triangles coincide. The same property is obtained in the stable region of super-critical speeds $2<v<4,$ where three one-periodic limit cycles are shown on the right side in Figure 2a. Again, the transients are marked by thin black lines. After a sufficiently long initial time, the transients go over to the stationary limit cycles marked by green, cyan, and yellow.

In Figure 2b, the applied wave parameter is doubled by $z_0\Omega =1$ with the consequence that the amplitudes of the speed oscillations become much broader. Moreover, the one-periodical limit cycles bifurcate into double-periodic ones when for growing driving force the vehicle becomes stuck in the under-critical speed range. This bifurcation scenario is demonstrated by the green limit cycle obtained (left side, Figure 2b) for the related driving force $f\Omega /c=0.3,$ which goes over to the yellow double-periodic one if the related force is increased to $f\Omega /c=0.6.$ For the further increased driving force $f\Omega /c=0.9$, the calculated limit cycle is still double-periodic marked by red and bifurcates back into the one-periodic limit cycle marked by cyan if the driving force is again increased to $f\Omega /c=1.2.$ For further growing driving force, the resonance speed $v=1$ is reached with strong oscillations in the vertical and horizontal direction and then passed up to high speed ranges where the driving force can be again decreased. Associated limit circles are single periodic as shown in the right side of Figure 2b. They are obtained for $f\Omega /c=2.7, 3.3, 3.9, \mathrm{and} 4.5,$ marked by green, yellow, red, and cyan, respectively. Correspondingly, colored triangles represent the averaged travel speeds together with the applied driving force. The applied vehicle damping is $D=0.2$ and the time step of integration is $∆\tau =2\times {10}^{-5}.$ Note that for $f\Omega /c=0.9$ two different limit cycles are obtained: the first in the under-critical speed range and the second one in the upper-critical speed range. Both are marked by red. In between both stable limit cycles, there is a middle one which is unstable und therefore not realizable. Note that one finds the same behavior in rotor dynamics, where the Sommerfeld rotor leads to the same driving characteristic [2] when in Equation (12) the road factor $z_0\Omega$ is replaced by the mass ratio of the unbalanced rotor, the translation velocity by the rotation speed ratio, and the driving force by the moment. More details on Sommerfeld effects in rotor dynamics are given in [16,17,18,19,20].

4. Stable and Unstable Oscillation Amplitudes of Vehicle Speeds

In Figure 3a, the vibration amplitudes of the vertical displacement (red) and velocity (pink) of the vehicle are plotted against the travel speed together with the speed driving force characteristic shown in blue. The displacement and velocity amplitude are

$$A_y=\sqrt{y_s^2+y_c^2}, \mathrm{and} A_x=\sqrt{x_s^2+x_c^2},$$

respectively, where the coefficients $y_s, y_c \mathrm{and} x_s,x_c$ are determined by Equations (10) and (11). The blue curve in Figure 3a represents the speed driving force characteristic dependent on speed. As already derived, stable travel speeds are only realizable in those speed regions where the slope of the characteristic is positive marked by thick blue lines in the under- and super-critical range. In the middle range, where the slope is negative, the blue curve is replaced by a thin black line, indicating that the calculated stationary speed is unstable and therefore not realizable. From this it follows that the calculated amplitudes $A_y \mathrm{and} A_x$ of the vertical vibrations of the vehicle are not realizable in the middle speed range, as well. In the super-critical speed range, however, all stationary vibration amplitudes are stabilized with growing speed and end into asymptotes marked by a dashed line, which are linearly growing, constant without increasing and decreasing for the force characteristic, the velocity, and the displacement amplitude, respectively. Obviously, the asymptote of the driving force corresponds to acceleration. Both quantities, force and acceleration, are linearly increasing with growing speed of the vehicle.

In order to introduce efficient solution methods, the time increment $d\tau$ in Equations (6) and (8) is eliminated by the angle increment $d\phi =-vd\tau$ that leads to the time-free equations

$$v\overline{y}°=-\overline{x}, v\overline{x}°=\left(\overline{y}+2D\overline{x}\right)-z_0\Omega \left(\overline{z}+2Dv\overline{u}\right),$$

(15)

$$vv°-2D{\left(z_0\Omega \right)}^2{\overline{u}}^{2}v={\left(z_0\Omega \right)}^2 \overline{zu}-z_0\Omega \left(\overline{uy}+2D\overline{ux}\right)-f\mathsf{\Omega }/c,$$

(16)

where prime in ${(.)}^{\prime }=d(.)/d\tau$ is replaced by a circle in $(.)°=d(.)/d\phi$ denoting differentiation with respect to the polar angle $\phi .$ In Equations (15) and (16), the related road level and slope are sinusoidal and dependent on the polar angle $\phi .$ Accordingly, the vertical vehicle vibrations and the longitudinal speed oscillations are calculable inside one period by means of shooting methods to obtain numerical solutions for $v_p\left(\phi \right),\overline{y}\left(\phi \right)$ $\mathrm{and} \overline{x}\left(\phi \right)$. The associated characteristic multipliers ${\rho }_i$ are investigated by means of the new equation system

$$v_p{\left[\begin{array}{c}∆\overline{y}\\ ∆\overline{x}\\ ∆v\end{array}\right]}^o+\left[\begin{array}{ccc}0& -1& {\overline{y}}_p^o\\ -1& -2D& {\overline{x}}_p^o+z_0\Omega 2D\overline{u}\\ z_0\Omega \overline{u}& 2D{z}_0\Omega \overline{u}& v_p^o-2D{\left(z_o\Omega \overline{u}\right)}^2\end{array}\right]\left[\begin{array}{c}∆\overline{y}\\ ∆\overline{x}\\ ∆v\end{array}\right]=0$$

which is obtained when the perturbations $\overline{y}={\overline{y}}_p+∆\overline{y}, \overline{x}={\overline{x}}_p+∆\overline{x},$ and $v=v_p+∆v$ are inserted into Equations (15) and (16) and the perturbation equations are linearized in the $∆$ quantities. According to the Floquet theory, the calculated periodic solutions are asymptotically stable with respect to $\phi$ if all multipliers ${\rho }_i$ satisfy $\left|{\rho }_i\right|<1.$ For known velocity $v_p\left(\phi \right)$, the applied polar angle can be integrated back to time by which the sequence of associated time distances is determined. Note that the four-dimensional time system (13) and (14) is not ergodic since the polar angle $\phi$ is a state variable that grows infinitely by permanent rotation. This disadvantage is avoided by eliminating the time variable by means of the polar angle, which leads to the three-dimensional angle system (15) and (16) where the polar angle now represents the independent integration variable restricted to one-periodic interval. The solutions of this new time-free equation system are ergodic, and multiplicative ergodic theorems are applicable in order to calculate characteristic numbers of the dynamic system of interest.

For the new angle Equations (15) and (16), new analytical solutions are derived by means of the introduction of the Fourier expansions

$$\overline{y}\left(\phi \right)=y_c\cos \phi +y_s\sin \phi +\cdots , \overline{z}\left(\phi \right)=\cos \phi$$

(17)

$$\overline{x}\left(\phi \right)=x_c\cos \phi +x_s\sin \phi +\cdots , \overline{u}\left(\phi \right)=-\sin \phi ,$$

(18)

$$v(\tau )=v_o+v_c\cos 2\phi +v_s\sin 2\phi +\cdots ,$$

(19)

where the zeroth Fourier coefficient $v_o$ of the speed expansion in Equation (19) takes the role of the averaged speed $v$ in Section 3. The insertion of the expansions (17) and (18) into Equation (15) and the coefficient comparison of the sinusoidal terms of $\cos \phi$ and $\sin \phi$ leads to the same result, already noted in Equations (10) and (11). All other terms are cutoff and can only be taken into account when higher expansions are introduced. The insertion of the expansions (17), (18), and (19) into the speed Equation (16) and the coefficient comparison leads to

$$v_s=\frac{D{\left(z_o\Omega \right)}^2{v}_o^2\left[{\left(z_o\Omega \right)}^2\left(1+4D^2{v}_o^2-v_o^4\right)-4v_o^4\right]}{2\left[{\left(v_o^2-1\right)}^2+{\left(2D{v}_o\right)}^2\right]\left[4v_o^2+{\left(z_o\Omega \right)}^4{D}^2\right]} ,$$

(20)

$$v_c=\frac{{\left(z_o\Omega \right)}^2{v}_o^3 \left\{1+v_o^2\left[4D^2-1+D^2{\left(z_o\Omega \right)}^2\right]\right\}}{\left[{\left(v_o^2-1\right)}^2+{\left(2D{v}_o\right)}^2\right]\left[4v_o^2+{\left(z_o\Omega \right)}^4{D}^2\right]} ,$$

(21)

which represent new results for the amplitudes of the travel speed oscillations. In Figure 3b, the resultant speed amplitude $A_v$ is plotted against the mean speed $v_o=v\Omega /{\omega }_1$ by

$$A_v=\sqrt{v_c^2+v_s^2} , \underset{v_o\to \infty }{\lim }{A}_v=D{\left(z_o\Omega \right)}^2[4+{\left(z_o\Omega \right)}^2]/8$$

(22)

and marked by a thick green line. Obviously, the speed amplitude vanishes if the vehicle slows. It is increasing up to the resonance speed and decreases again for further growing speed, up to the asymptote given in Equation (22). In addition to the above coefficient comparison of terms with $\cos 2\phi$ and $\sin 2\phi$, all terms with $\cos 0$ lead to the extended speed driving force characteristic

$$f\Omega /c=\frac{D{\left(z_o\Omega \right)}^2{v}_o^3}{{\left(v_o^2-1\right)}^2+{\left(2D{v}_o\right)}^2}\{v_o^2+{\left(z_o\Omega \right)}^2\frac{v_o^2\left[1-{\left(z_o\Omega \right)}^2{D}^2-4D^2\right]-1}{2\left[4v_o^2+{\left(z_o\Omega \right)}^4{D}^2\right]} \}$$

(23)

where the first part coincides with Equation (12) and the second part gives a correction of second order. The extended speed driving force characteristic in Equation (23) is plotted in Figure 3b. The red line marks the first approximation noted in Equation (12). Its second order approximation, noted in Equation (23), is marked by a thick cyan line. It is close to the red line of the first approximation. Thin black lines represent speed driving force characteristics with negative slope where the calculated solutions are unstable and therefore not realizable.

5. Stability in Mean and Robustness with Regard to Disturbances

Note that the velocity Equation (13) is determined by the products of road level z and slope u multiplied by the two vibration coordinates $y$ and $x=y\prime$. These products are called covariances [2]. They determine the dependency of the response of the system from the road excitation. The application of the product rule of differentiation to the road Equation (7) and the vehicle Equation (8) leads to the covariance equation system

$${\left[\begin{array}{c}\overline{zy}\\ \overline{uy}\\ \overline{zx}\\ \overline{ux}\end{array}\right]}^{\prime }+\left[\begin{array}{cccc}0& -v& -1& 0\\ v& 0& 0& -1\\ 1& 0& 2D& -v\\ 0& 1& v& 2D\end{array}\right]\left[\begin{array}{c}\overline{zy}\\ \overline{uy}\\ \overline{zx}\\ \overline{ux}\end{array}\right]=z_0\Omega \left[\begin{array}{c}0\\ 0\\ {\overline{z}}^2+2Dv\overline{zu}\\ \overline{zu}+2Dv{\overline{u}}^2\end{array}\right],$$

(24)

where the related speed $v$ of the vehicle is determined by Equation (13) as

$$v^{\prime }+2D{\left(z_0\Omega \right)}^2{\overline{u}}^{2}v=z_0\Omega \left(\overline{uy}+2D\overline{ux}\right)-{\left(z_0\Omega \right)}^2 \overline{zu}+f\mathsf{\Omega }/c .$$

(25)

Level and slope of the road are $\overline{z}\left(\tau \right)=\cos \left(v\tau \right)$ and $\overline{u}\left(\tau \right)=-\sin \left(v\tau \right),$ respectively. Again, this is a DA equation system where transcendentals take the role of algebraic terms. For numerical integrations, the sinusoidal terms in Equation (3) are generated by means of the road oscillator, noted in Equation (7). The oscillator is transformed by means of the polar coordinates

$$\overline{z}\left(\phi \right)=\cos \phi , \overline{u}\left(\phi \right)=-\sin \phi , {\phi }^{\prime }=-v,$$

(26)

which are inserted into Equations (24) and (25). The remaining equation ${\phi }^{\prime }=-v$ determines the derivative of the polar angle and couples the polar angle back to the speed of the vehicle.

In order to derive analytical approximations, the averaging method is applied by taking the time mean values $<{\overline{z}}^2>=<{\overline{u}}^2>=1/2$ and $<\overline{zu}>=0$ in Equations (24) and (25) that leads to

$$(1-v^2)\overline{uy}+2D\overline{ux}=z_o\Omega vD$$

$$-2{Dv}^2 \overline{\mathrm{uy}}+\left(1-v^2\right) \overline{\mathrm{ux}}={ z}_o\Omega v/2,$$

when in Equation (24) for ${(.)}^{\prime }=0$ the first two rows are applied to eliminate $\overline{zy}$ and $\overline{zx}.$ Subsequently, the covariances $\overline{uy}$ and $\overline{ux}$ are calculated and inserted into Equation (25)

$$D{\left(z_0\Omega \right)}^{2}v=z_0\Omega \left(\overline{uy}+2D\overline{ux}\right)+f\mathsf{\Omega }/c,$$

which finally leads to the same speed driving force characteristic

$$f\mathsf{\Omega }/c={\left(z_0\Omega \right)}^{2}D{v}^5/[{\left(v^2-1\right)}^2+{\left(2Dv\right)}^2] ,$$

already noted in Equation (12). Note that the above applied time mean values are exact when the speed of the vehicle is constant. For oscillating speeds, the time mean values are approximations.

Following [2], the averaged Equation (24) are applied to investigate the stability in mean of the averaged speed by means of the perturbation equation system

$$[\begin{array}{c}∆v\\ ∆\overline{zv}\\ ∆\overline{uy}\\ ∆\overline{zx}\\ ∆\overline{ux}\end{array}]+\left[\begin{array}{ccccc}D{\left(z_0\Omega \right)}^2& 0& -z_0\Omega & 0& -2D{z}_0\Omega \\ -\overline{u{y}_0}& 0& -v_0& -1& 0\\ \overline{z{y}_0}& v_0& 0& 0& -1\\ -\overline{u{x}_0}& 1& 0& 2D& -v_0\\ \overline{z{x}_0}-z_0\Omega D& 0& 1& v_0& 2D\end{array}\right][\begin{array}{c}∆v\\ ∆\overline{zv}\\ ∆\overline{uy}\\ ∆\overline{zx}\\ ∆\overline{ux}\end{array}]=0$$

(27)

which is obtained by means of the perturbations $v=v_0+∆v$ and $\overline{\alpha \gamma }=\overline{\alpha {\gamma }_0}+∆\overline{\alpha \gamma }$ for $\left(\alpha ,\gamma \right)\in \left\{z,u,y,x\right\}$. The insertion of the perturbations into the averaged Equations (24) and (25) and linearization in the $∆$-terms yields the fifth order stability Equation (27) with the characteristic equation $∆\left(\lambda \right)=$$A_0{\lambda }^5+A_1{\lambda }^4+\cdots A_5=0.$ The determinant of Equation (27) yields

$$A_5=\frac{{\left(z_0\Omega \right)}^{2}D{v}_0^4\left[v_0^4-6\left(1-2D^2\right)v_0^2+5\right]}{{\left(1-v_0^2\right)}^2+{\left(2D{v}_0\right)}^2} .$$

(28)

According to the Hurwitz criterion, $A_5<0$ determines divergence that gives the boundary of monotonous instability. Obviously, the stability condition $A_5>0$ coincides with the positive slope condition of the speed driving force characteristic calculable by differentiating Equation (12) with respect to the travel speed $v$ that gives

$$\frac{d}{dv}\left(\frac{f\mathsf{\Omega }}{c}\right)={\left(z_0\Omega \right)}^2{v}^{4}D\frac{v^4-6\left(1-2D^2\right)v^2+5}{{\left[{\left(1-v_0^2\right)}^2+{\left(2Dv\right)}^2\right]}^2}>0 .$$

This result coincides with Equation (28) except that the positive definite denominator in Equation (28) is squared. Hence, the negative slope condition and the instability in mean lead to the same instability boundary. This is plausible and physically explainable, as already done in the yellow area in Figure 2a. The instable speed range vanishes with increasing damping.

In addition to the stability behavior, the robustness of the limit cycle calculation with regard to disturbances is of interest. When disturbances are generated by brief shocks, the eigenvalues of the stability matrix in Equation (27) must be calculated. They determine the growth behavior with which the disturbances increase or decrease, respectively. In the case of stationary disturbances, noise models are introduced. Figure 4a shows a stochastic limit cycle obtained for the quarter car model in the case that the angle motion on the sinusoidal road form is perturbed by additive noise given by

$$d{\varphi }_{\tau }=-v_{\tau }d\tau +\sigma d{W}_{\tau }, \mathsf{\Delta }{W}_n=\sqrt{∆\tau } N_n.$$

(29)

In Equation (29), capital letters with index $\tau$ denote set functions [21,22] dependent on time. Noise is generated by normally distributed numbers $N_n$ with zero mean [23]. The stochastic angle perturbation takes into account that the road surface is no longer sinusoidal but more realistically irregular and noisy with bounded realizations. For small noise intensities $\sigma ,$ this leads to response realizations which are bounded, as well.

Initial results are shown in Figure 4a where trajectories of a stochastic limit cycle are plotted in the phase plane of travel velocity and acceleration scaled by 0.3 and shifted by the applied driving force. The realizations are calculated by means of Equations (24)–(26) where the polar angle in Equation (26) is replaced by Equation (29). The applied damping is given by $D=0.2,$ the driving force by $f\Omega /c=0.6,$ the road level by $z_o\Omega =0.9,$ and the noise intensity by $\sigma =0.03.$ The Euler scheme is applied with the time step $∆\tau ={10}^{-4}.$ The mean value of the shifted acceleration and true velocity is marked by a yellow triangle on the red curve of the speed driving force characteristic, indicating that the mean acceleration is vanishing and the mean travel speed coincides with Equation (12). The comparison with the double-periodic yellow limit cycle in Figure 2b shows that its sharp line is widened to a bundle of nonperiodic realizations, the boundaries of which are double periodic with two loops and one node of two crossing limit flows. Figure 4b shows the associated double-crater-like probability distribution density on the phase plane of velocity and acceleration calculated for the stronger noise intensity $\sigma =0.05.$ Clockwise rotation in the phase plane of acceleration and true speed is slow when the probability density is high and vice versa fast for low densities. Note that Figure 4a,b present new results in nonlinear stochastic vehicle dynamics, obtained by applying noise perturbations that are bounded by means of the applied sinusoidal terms. For small intensity σ, the limit flow in Figure 4a possesses periodic side limits in spite of the fact that all limit cycle realizations are random. This is a new effect presented in this paper. For growing intensity σ, the inner side limits of the limit cycle shown in Figure 4a become broader and finally disappear, such that only the outer side limits remain, and the whole phase plane in Figure 4a is covered by realizations of velocity and acceleration.

6. Stable Travel Speeds of Road Vehicle Systems

To extend the above investigations to higher order models, consider the quarter car model with two degrees of freedom shown in Figure 5a. The motions of the car body of mass $M$ and the wheel mass $m$ are described by the equations of motion

$$\ddot{y}+2D{\omega }_1\left(\dot{y}-\dot{x}\right)+{\omega }_1^2\left(y-x\right)=0, z\left(s\right)=z_0\cos \left(\Omega s\right),$$

(30)

$$\ddot{x}+2B{\omega }_2\left[\left(1+\beta \right)\dot{x}-\dot{z}-\beta \dot{y}\right]+{\omega }_2^2\left[\left(1+\gamma \right)x-z-\gamma y\right]=0,$$

(31)

where the system frequencies ${\omega }_1$, ${\omega }_2$ and the damping coefficients $D, B$ are introduced together with the stiffness and damping ratios $\gamma \mathrm{and} \beta ,$ and defined as follows:

$${\omega }_1^2=c/M, 2D{\omega }_1=d/M, \gamma =c/k,$$

$${\omega }_2^2=k/m, 2B{\omega }_2=b/m, \beta =d/b.$$

The parameter $\gamma$ denotes the stiffness ratio of the car and wheel spring $c$ and $k$, respectively. Correspondingly, $\beta$ is the damping ratio of the car and wheel damper $d$ and $b$, respectively. The reference frequencies ${\omega }_2$ and ${\omega }_1$ describe the decoupled vibrations of the wheel and car body. In addition to Equations (30) and (31), the dynamic balance in horizontal direction gives a third equation of motion that determines the travel speed, as follows:

$$\left(M+m\right) \dot{v}=f+\left[k\left(x-z\right)+b\left(\dot{x}-\dot{z}\right)\right] \tan \mathsf{\alpha } , \tan \alpha =dz/ds.$$

(32)

Note that Equation (32) is of first order with respect to the car speed $v,$ and s denotes the longitudinal coordinate of the travel path. Note that both masses are assumed to be concentrated in the contact point of road and vehicle so that only planar translations are considered. Rotations are excluded.

It is appropriate to introduce the dimensionless vibration and road coordinates by means of $\left(y,x\right)=(\overline{y},\overline{x})/\Omega$ and $\left(z,u\right)=z_o\left(\overline{z},\overline{u}\right),$ respectively. The insertion of these coordinates into Equations (30) and (31) leads to the dimensionless equations of motion

$$\ddot{\overline{y}}+2D{\omega }_1\left(\dot{\overline{y}}-\dot{\overline{x}}\right)+{\omega }_1^2\left(\overline{y}-\overline{x}\right)=0, v=v\Omega /{\omega }_1,$$

(33)

$$\ddot{\overline{x}}+2\beta {\omega }_2\left[\left(1+\beta \right)\dot{\overline{x}}-\beta \dot{\overline{y}}\right]+{\omega }_2^2\left[\left(1+\gamma \right)\overline{x}-\gamma \overline{y}\right]=z_o\Omega {\omega }_2^2\left(\overline{z}+2Bv\overline{u}\right),$$

(34)

where $v$ is the related speed of the vehicle rolling on road with level $\overline{z}=\cos \Omega s$ and slope $\overline{u}=-\sin \Omega s.$ In order to derive a first approximation, it is assumed that the oscillating speed of the vehicle can be averaged by $v\left(\tau \right)=v=\mathrm{const}.$ In this case, the travel path is s $=\mathrm{vt}$ and the equations of motion become linear. They are solved by the set-up

$$\overline{y}\left(t\right)=y_c\cos \left(v{\omega }_{1}t\right)+y_s\sin \left(v{\omega }_{1}t\right), \overline{z}\left(t\right)=\cos \left(v{\omega }_{1}t\right),$$

$$\overline{x}\left(t\right)=x_c\cos \left(v{\omega }_{1}t\right)+x_s\sin \left(v{\omega }_{1}t\right), \overline{u}\left(t\right)=-\sin \left(v{\omega }_{1}t\right)$$

In the stationary case, the insertion of these set-ups into Equations (33) and (34) and the coefficient comparison leads to the linear matrix equation

$$\left[\begin{array}{cccc}1& 2Dv& v^2-1& -2Dv\\ -2Dv& 1& 2Dv& v^2-1\\ 1+\gamma -v^2{\kappa }^2& 2B\left(1+\beta \right)v\kappa & -\gamma & -2B\beta v\kappa \\ -2B\left(1+\beta \right)v\kappa & 1+\gamma -v^2{\kappa }^2& 2B\beta v\kappa & -\gamma \end{array}\right]\left[\begin{array}{c}{x}_c\\ x_s\\ y_c\\ y_s\end{array}\right]=z_o\Omega \left[\begin{array}{c}0\\ 0\\ 1\\ -2Bv\end{array}\right]$$

(35)

where $\kappa$ is the reference frequency ratio given by $\kappa ={\omega }_2/{\omega }_1.$ The coefficients $x_c, x_s, y_c, y_s$ of the sinusoidal set-up are calculated by means of Equation (35) and inserted into

$$A_y=\sqrt{y_s^2+y_c^2}, \mathrm{and} A_x=\sqrt{x_s^2+x_c^2},$$

which are plotted in Figure 5b against the related speed of the vehicle. The amplitude $A_y$ of the car body is marked by red and the wheel amplitude $A_x$ by pink. In unstable speed ranges, both amplitudes are marked by thin black lines.

The vertical vibration amplitudes of both masses are calculated for the reference frequencies ${\omega }_1=1/s,$ ${\omega }_2=2/s, \left(\kappa =2\right),$ for the stiffness-ratio $\gamma =0.5$ and for the damping values $D=0.3, B=0.2,$ and $\beta =0.5 .$ The applied road wave is chosen by $z_o\Omega =1.$ This is equivalent to the height–length ratio $z_o/L=0.159$ when the wave frequency $\Omega =2\pi /L$ is inserted into $z_o\Omega =1$. The frequency responses in Figure 5b possess two resonance speeds: the first at $v_1=0.27$ and the second at $v_1=1.29.$ Both eigenvalues are valid for vanishing damping. In order to decide the realizability of the stationary vehicle speeds applied in Equation (35), the new speed Equation (32) is investigated in the dimensionless form

$$\frac{m+M}{m}\dot{v}/{\omega }_2+2B{\left(z_o\Omega \right)}^2{\overline{u}}^{2}v=f\Omega /k+z_o\Omega (\overline{xu}+2B\overline{u}\dot{\overline{x}}/{\omega }_2)-{\left(z_o\Omega \right)}^2\overline{uz} .$$

(36)

As already shown in Section 3, Equation (36) can approximately be evaluated by averaging the oscillating speed with $\dot{v}=0$ and introducing the time mean values $\overline{uz}=0,$ $\overline{xu}=-x_s/2,$ and $\overline{u}\dot{\overline{x}}/{\omega }_2=\kappa v{x}_c/2$ into Equation (36). This gives the new approximated force characteristic

$$f\Omega /k=B{\left(z_o\Omega \right)}^{2}v+z_o\Omega \left(x_s-2B\kappa v{x}_c\right)/2 ,$$

(37)

where $x_s$ and $x_c$ are calculated by means of Equation (35). New numerical evaluations of Equations (35) and (37) are plotted in Figure 5b and marked by green for stable speeds when the slope of the speed driving force characteristic is positive and by black-thin lines when the slope is negative. In the latter case, the speed applied in Equation (35) is unstable and physically not realizable. Finally, it is noted that both instability speed ranges are vanishing for growing damping coefficients.

The instability behavior is numerically investigated by means of simulations applied to a slightly modified quarter car model with two equal masses m and a third damper $B$, as shown in Figure 6a. The car is rolling with velocity $v$ on a wavy road described by $z\left(s\right)=z_{o }\cos \Omega s$ in dependence on the travel way coordinate $s$. The speed frequency $v\Omega$ is related to the reference frequency ${\omega }_1$ given by ${\omega }_1^2=c/m$. The vertical vibrations of this car model are described by the two nonlinear equations of motion

$${\ddot{y}}_1+2d_1{\omega }_1\left({\dot{y}}_1-{\dot{y}}_2\right)+2D{\omega }_1\left({\dot{y}}_1-\dot{z}\right)+{\omega }_1^2\left(y_1-y_2\right)=0,$$

(38)

$${\ddot{y}}_2+2d_2{\omega }_1\left({\dot{y}}_2-\dot{z}\right)-2d_1{\omega }_1\left({\dot{y}}_1-{\dot{y}}_2\right)+{\omega }_1^2\left(2y_2-y_1-z\right)=0,$$

(39)

where $\dot{z}=v\Omega u$ couples both vertical vibrations to the horizontal velocity described by

$$2\dot{v}+{\omega }_1[{\omega }_1\left(z-y_2\right)+2d_2(\dot{z}-{\dot{y}}_2)+2D\left(\dot{z}-{\dot{y}}_1\right)]\Omega u=f/m.$$

(40)

The damping coefficients are introduced by $B/m=2D{\omega }_1 \mathrm{and} b_i/m=2d_i{\omega }_i$ for $i=1,2.$ For numerical integration, the dimensionless time $\tau ={\omega }_{1}t$ and the velocity $v=v\Omega /{\omega }_1,$ together with the related coordinates $\left(\overline{z}, \overline{u}\right)=\left(z,u\right)/z_o$ and $\left({\overline{y}}_i,\overline{x_i}\right)=\left(y_i,x_i\right)\mathsf{\Omega },$ are inserted into the velocity Equation (40) and vehicle Equations (38) and (39).

Figure 6b shows the average vehicle speed for a given driving force $f\Omega /c$, marked by a blue thick line. It is calculated for the road level $z_o\Omega =1$ and the damping values $D=d_1=d_2=0.1$. In the middle range of the related driving force $\left(0.60<f\Omega /c<0.82\right)$, there are five branches of stationary solutions: two solution branches are unstable and three are stable. The stable speeds possess positive slopes in the driving force speed characteristic; meanwhile, speeds started in negative slope ranges do not remain in this range. They migrate to higher or lower speed ranges. Figure 6b shows three limit cycles (red, yellow, green) before the first resonance, two further limit cycles (red, cyan) between the first and second resonance and the last three limit cycles for overcritical speeds near the asymptote $f\Omega /c=\left(d_2+D\right)v,$ which is marked by a dashed red line. Mean values of shifted acceleration and true speed are plotted by colored triangles on positive slopes.

7. Summary and Remarks on Main New Results

When vehicles are rolling on uneven roads and driven by a constant force, vertical vibrations of the vehicle are induced that are coupled with the horizontal motions of the vehicle, affecting the travel velocity, which fluctuates around a stationary speed with zero mean acceleration. The coupling between both planar motions is caused by the permanent direction change of the contact force to ground along the contour of the road profile. In the case of sinusoidal road profiles, the vehicle becomes stuck before the resonance speed and needs more power supply to overcome the resonance in order to reach higher speeds. The speed range immediately after the resonance where the slope of the speed driving force characteristic is negative, is not realizable. In this range, the calculated vertical vibration amplitudes are unstable. This is completely different from the results in linear road vehicle dynamics. To describe the dynamic effects obtained, the following main new results are presented:

• Multiple Sommerfeld effects are shown before and between all resonances. Regular speeds of operation are far beyond. The periodic limit cycles in Figure 6b are new results showing how the car speed becomes stuck on positive slopes of the driving force speed characteristic.
• In Figure 4a,b, new stochastic limit cycles are shown. For small intensities of white noise, a bundle of realizations is obtained described by a periodic probability distribution over the phase plane of the oscillating true speed and zero mean travel acceleration.
• Fourier expansions are introduced into the nonlinear equations of motion in order to calculate new analytical approximations for the mean travel speed and the amplitude of the speed oscillations. These expansions depend on the polar angle which takes the role of time.

In future, the one wheel vehicle is extended to half car models by which additional effects are introduced by road excitations with time delays when there is a finite distance between the front and rear wheel. Analytical investigations by means of Fourier expansions are applied to higher order vehicle models. For numerical solutions, shooting methods restrict the integration range to one period of the polar angle by which unstable period limit cycles can be calculated, as well. The stability of the periodic solutions is investigated by means of the characteristic multipliers calculated by means of the Floquet theory.