Turbulent travel speeds in nonlinear vehicle dynamics

Quarter car models of vehicles rolling on wavy roads lead to limit cycles of travel speed and acceleration with period doublings and bifurcation effects for appropriate driving force parameters. In case of narrow-banded road excitations, speed jumps occur, additionally. This has the consequence that the driving speed becomes turbulent. Bifurcation and jump effects vanish with growing vehicle damping. The same happens for increasing bandwidth of road excitations when, e.g., on flat highways there are no big road waves but only small noisy slope processes generated by rough road surfaces. The paper derives a new stability condition in mean. Numerical time integrations are stabilized by means of polar coordinates. Equivalently, Fourier series expansions are introduced in the angle domain. Phase portraits of travel speed and acceleration show new period-doublings of limit cycles when speed gets stuck before resonance. The paper extends these investigations to the stochastic case that road surfaces are random generated by filtered white noise. By means of Gaussian closure, a nonlinear mean speed equation is derived which includes the extreme cases of wavy roads and road noise.


Introduction into nonlinear vehicle dynamics
shows the applied quarter car model rolling on wavy roads with level z and slope u that effects coupled vibrations of vertical displacement y and horizontal travel velocity v [1][2][3][4] described bÿ v = ω 2 1 (y − z) + 2Dω 1 (ẏ −ż) tanα + f /m. (2) Equations (1) and (2) are equations of motion where dots denote derivatives with respect to time t and ω 1 is the natural frequency of the vehicle given by ω 2 1 = c/m. The parameter D determines the vehicle damping by 2Dω 1 = b/m, and m is the vehicle mass driven by the force f assumed to be constant. When α is the slope angle in the frictionless contact point of road and vehicle, the vertical component Ncosα of the normal force N is determined by the linear forces of spring c and damper b so that the horizontal component Nsinα enters through tan α into the horizontal dynamic equilibrium, noted in Eq. (2). According to mass-point mechanics, the vehicle mass is assumed to be concentrated in the contact point. Consequently, there is no rotation, only two planar velocities: the horizontal travel velocity v and the vertical vibration velocityẏ. The latter can be substituted byẏ = ω 1 x what is later applied in order to rewrite the linear second-order equation (1) of motion into the state form of two first-order system equations.
Note that Eq. (2) is approximately investigated in [5,6] by means of averaging methods which are only To obtain exact solutions, the one-dimensional road surface is modeled, e.g., by means of the sinusoidal forms [7,8] where z 0 is the amplitude of the road wave, s (t) determines the current position of the contact point of the vehicle moving with the horizontal velocityṡ = v and denotes the road frequency measurable by the wave length L = 2π/ . From Eq. (3), the slope angle α(t) of the road is calculable by dz/ds that leads to tan α = dz/ds = u. (4) After the insertion of tan α = u into Eq. (2), both equations of motion have to be solved under the condition that the geometrical equations (3) are satisfied. This is typical for DA equation systems where the sinusoidal functions (3) in the given problem take the role of algebraic terms in DA equations [9]. To eliminate the growing motion coordinate s (t), the geometrical conditions (3) are replaced by the dynamical relationṡ which are derived from Eq. (3) by means of dz = − z 0 sin( s) ds and du = − z 0 cos( s) ds where −z 0 sin ( s) is replaced by u and z 0 cos ( s) by z and ds/dt is substituted by v. Altogether, Eq. (1), (2) and (5) describe a five-dimensional problem with five unknowns: the horizontal velocity v (t) of the vehicle, its vertical vibration by y (t) andẏ (t) = ω 1 x (t) and the road level z (t) and slope u (t) in the moving contact point. When v (t) is known, the moving contact coordinate s (t) can be calculated by the integration for any given initial position s 0 at the time t 0 . In the special case that the velocity v of the vehicle is constant withv = 0, the nonlinear dynamical equation (2) degenerates to a static one. Furthermore, Eq. (5) becomes linear for v = const. describing the motion of an oscillator with the constant time frequency v .

Limit cycles of travel speed and acceleration
For the control problem or inverse problem that the driving force is unknown and the vehicle is moved with given constant speed v, the road and vehicle equations (1) and (5) become linear and solvable by means of the sinusoidal solutions where V and ε are amplitude and phase, respectively, to be calculated from the linear equation (1) of motion. The insertion of the calculated solutions z (t) , u (t) and y (t) into Eq. (2) gives the explicit form of the sinusoidal driving force f (t) necessary to maintain a constant travel speed. The time average of f (t) leads to the mean force In Fig. 2, the dimensionless mean force (6) is plotted versus the dimensionless speed ν = v /ω 1 for the two dimensionless road factors z 0 marked in red and cyan together with two typical motor characteristics in yellow. The first one is valid for a constant driving force. The second one holds for a more realistic driving force which decreases with growing speed. The two cutting points of driving characteristics and mean forces represent possible stationary velocities marked by circles before and after the related resonance velocity ν = 1.
For the nonlinear case of fluctuating vehicle velocities, the simple result in Eq. (6) is confirmed in [5] by means of an averaging method. One finds the same result (6) in rotor dynamics [10][11][12][13] when the road factor z 0 is replaced by the mass ratio of an elastically supported unbalanced rotor and the translation velocity v is substituted by a rotation speed. In the following, numerical integration methods and nonlinear Fourier expansions are developed in order to extend the averaged solution in Eq. (6) to limit cycles of acceleration and travel velocity. First results are presented in Fig. 3 showing new limit cycles obtained by means of numerical integration of Eqs. (1), (2) and (5) applying Euler schemes for short dimensionless times τ = ω 1 t. The dimensionless acceleration ν = dν/dτ is plotted against the velocity ν = v /ω 1 . The obtained results are shown in an appropriately scaled form where the dimensionless acceleration is divided by five and lifted by adding the dimensionless force value, applied. Therewith, all drawings of the limit cycles are more separated in Fig. 3 and degenerate to small circles around the driving force when the road factor decreases. This is due to the fact that the time mean value of the calculated acceleration is zero in the stationary case. The results in Fig. 3 are calculated for eight different parameters f /c of the driving force marked by colored triangles on the red back bone curve of the driving force applied. The geometrical form of the road is chosen by the road factor z 0 = 1. This indicates that the applied wave amplitude of the road is z 0 = 10/π cm for the wave length given by L = 20 cm. The selected damping is D = 0.2. The simulations are started with initial values marked by squares. Subsequent motions are plotted by thin dashed lines. After a sufficiently long transient period, the simulations end in limit cycles marked by colored thick lines. Limit cycles of travel speeds with mean speeds on negative slopes are unstable and therefore not realizable, physically.

Mean travel speed and stability in mean
Provided that the vehicle is moving with velocity ν > 0, Eqs. (7) and (8) are averaged by taking the mean values z 2 = ū 2 = 1/2 and zu = 0 that leads to the mean drift equation (9) for the mean velocity m ν = v and to the averaged covariances matrix equation for the covariances of the vehicle vibration y (t) , x (t) multiplied by the road level z (t) and slope u (t). Since the system matrix in Eq. (10) is skew-symmetric, the determinant of Eq. (10) is positive definite determining the mean values and zy 0 = ux 0 /m ν and zx 0 = − uy 0 m ν . The insertion of Eq. (11) into Eq. (9) leads to by which the travel velocity mean m ν can be calculated for given driving force f /c, road factor z 0 and damping D. The result (12) formally coincides with Eq. (6) of the inverse or control problem Figure 4 shows simulation results of Eqs. (9) and (10) where the acceleration is lifted by adding the applied driving force (6) and plotted versus the vehicle velocity for different initials marked by squares. After sufficiently long times, all simulated curves end in red circles on the thick red line calculated by Eq. (12). Note that stationary mean accelerations and travel speeds end The stability of stationary mean velocities is investigated by means of the perturbations m ν = m 0 + m ν and αγ = αγ 0 + αγ for all state processes (α, γ ) ∈ {z, u, y, x}. Insertion of all perturbations into Eqs. (9) and (10) and linearization in -terms lead to the fifth-order stability equation All stationary covariances calculated in Eq. (11) are inserted into Eq. (13) what leads to the characteristic The determinant of Eq. (13) yields According to the Hurwitz Criterion, A 5 < 0 determines divergence and gives the boundaries of monotonous instability. The new instability condition A 5 < 0 determined by Eq. (14) coincides with the negative slope of the driving force calculable by differentiating Eq. (12) with respect to the mean travel speed m ν . The coincidence of both, negative slope and instability in mean, is plausible and physically explainable: A perturbation of stable speed into the negative direction on the left side of the dynamic equilibrium generates acceleration since the applied force is bigger than that one of the new driving force of the perturbed speed. The vehicle, however, is braked if the speed perturbation goes into positive speed direction on the right side of the dynamical equilibrium since then the driving force is smaller than that one necessary to maintain the perturbed dynamic equilibrium. Vice versa, driving forces with negative slope lead to monotonous instability with the effect that speed leaves the unstable branch of the driving force, applied.

Stabilized integration by means of polar coordinates
For constant speeds, the road equations (5) The insertion of the coordinates (15) into the dimensionless forms dz = νūdτ and dū = −νzdτ of the road equations (5) gives the transformed road equations where the related polar radius is integrated tor = 1 by which drift is eliminated and integration is stabilized.
Equations (15) and (16) are inserted into Eq. (7) that gives the new covariance matrix equation Note that in both Eqs. (17) and (18) the increment dτ of the dimensionless time is replaced by the polar angle increment dϕ given by Eq. (16). Consequently, the integration of both equations can be performed with respect to time or angle. In Fig. 5, the numerical integration of Eqs. (17) and (18) is performed by means of Euler schemes in the time domain applying the time step size τ = 10 −4 . The time integration is started in the blue square inside the limit cycle and ends in the yellowblack limit cycle around the yellow triangle on the driving force calculated by Eq. (6). The transient behavior from starting point up to the stationary limit cycle is marked by a blue dashed line. The applied road factor and driving force are z 0 = 0.7 and f /c = 0.3, respectively. In addition, Fig. 5 shows results of a second simulation marked by a dashed green line which represents the moving average of the simulation results. This line starts in the green square outside the limit cycle and ends in the yellow triangle on the driving force calculated by Eq. (12). Moving average results are already shown in Fig. 4 where the averaged equations (9) and (10) are solved, numerically. Both curves, the black one in Fig. 4 and the green in Fig. 5, are slightly different because the order of averaging and integration is reversed in both figures. Finally, Fig. 6 shows the period-doubling effect obtained for the parameters z 0 = 0.9, f /c = 0.6 and D = 0.2. Similar perioddoubling effects are obtained when the travel velocity is plotted against the vibration velocity. The doubleperiodic limit cycle in Fig. 6 degenerates to the oneperiodic limit cycle already shown in Fig. 5 when the road factor z 0 and the driving force f /c are decreasing. For analytical investigations, Eqs. (17) and (18) are solved by means of the Fourier expansions where k is the vector of the covariances (zy, uy, zx, uz). The insertion of the two Fourier series (19) into the drift equation (18) and covariances equation (17) gives the two first equations of expansion where the vectors p and q are determined by the righthand side of Eq. (17) and the matrices A and B by its left-hand side. Both equations lead to the same velocity mean already calculated by averaging in Eq. (12) when higher terms of the Fourier expansion are neglected.
Summarizing it is noted that the longitudinal road vehicle problem is solved according to the theory of DA equations by eliminating the infinitely growing vehicle motion coordinate. This is possible when the geometrical road equations are replaced by means of corresponding dynamical ones. The calculated stationary solutions are evaluated in phase portraits plotting the horizontal acceleration versus the travel velocity. The numerical simulations show stable limit cycles before and after the resonance velocity, respectively. The middle limit cycles on negative slopes of the applied driving force are unstable and only detectable by means of analytical methods, e.g., by means of nonlinear Fourier expansions. For increasing driving forces needed for growing road factor, the one periodic limit cycle bifurcates into a double periodic one. For further growing driving forces, the limit cycles become one periodic, again.

Random road profiles in space and time
In order to characterize wave length L = 2π/w and way frequency w of level and slope of stochastic road processes, measurements [14] of one-dimensional profiles are evaluated by means of power spectral densities [15][16][17] plotted versus the road frequency w in doublelogarithm scales. Random processes of road profiles are generated, e.g., by means of white way noise W s with intensity σ applied to the oscillator equation where δ > 0 denotes the bandwidth of the road process and the parameter denotes its center frequency. All processes in Eq. (20) are written in capital letters as set functions in dependence on the way coordinate s. Dashes denote derivatives with respect to s The input term in Eq. (20) is delta-correlated with zero mean. It is uniformly distributed in the frequency range [18] with the power spectrum S(w) = 1. The Fourier transform of Eq. (20) multiplied by its conjugate version gives the road power spectrum [19] It determines the frequency distribution of the vertical level of the road profile in way domain. In Figs. 7 and 8, evaluations of the above spectrum are marked by  smooth red lines plotted in double-logarithmic scales for two sets of the parameters and δ. The noise intensity in both figures is chosen to be σ = 1. Cyan lines belong to power spectra of the slope process U s . The ordinary differential equation of the oscillator (20) is rewritten into Itô differential equations [20,21] dZ s = U s ds, E dW 2 s = ds, where the Wiener increment dW s is approximated [22] by the square root of the way increment s multiplied by normally distributed numbers N n with zero mean E (N n ) = 0 and unit mean square E(N 2 n ) = 1.
The application of the Maruyama scheme [23] to Eqs. (22) and (23) leads to M samples of finite signals with length L = N s and N simulated values Z k to which the spectral analysis is numerically performed by means of the discrete Fourier transformū nm of all samples multiplied by its conjugate complex version u * nm . The product of both is averaged over all samples M in order to obtain the discrete power spectrum where c z is a scaling factor needed to satisfy Kolmogorov's consistency condition. The scaling is carried out applying the mean square σ 2 z to be additionally measured from all Z k and integrating over all frequencies w n = n w L with the step size w L = 2π/L on the right side of Eq. (25). In Figs. 7 and 8, the FFT averages are plotted versus the frequency w. They are marked by rough black lines overlaying the smooth analytical results. The rough results become smoother for growing number M of all samples.
When vehicles roll on road with nonnegative velocityv ≥ 0, the contact point between road and vehicle is moving with same velocity. Consequently, the road equations (22) and (23) must be transformed from way to time by means of the vehicle velocityv applying the classical way-time relation ds =vdt. According to Itô's calculus, the Wiener increments in way and time are scaled by dW s = √v dW t [19] where E dW 2 s = ds and E dW 2 t = dt. The introduction of both way-time relations into the road model (22) and (23) leads to transformed road processes Z t and U t described by the first order dynamical system where the parameters δ andv denote the bandwidth of the road and its middle time frequency, respectively. The application of the two-dimensional Fokker-Planck equation [21] to Eq. (26) yields This equation determines the joint density distribution of the two road processes in dependence on time. In the stationary case ∂ p/∂t = 0, the velocityv drops out with the consequence that the stationary density p(zu) is independent on velocity. Consequently, the road model (26) is applicable for all velocitiesv ≥ 0. This is exactly the same situation as, e.g., in the deterministic case where the amplitude of the harmonic wave road remains unchanged when the vehicle speeds up or slows down.

Quarter car models rolling on random road surfaces
Quarter car models on random road are excited to vertical vibrations described by the equation of motion where 2Dω 1 = b/m determines the car damping and ω 1 is the natural frequency of the vertical vibrations, given by ω 2 1 = c/m . In a first investigation, the velocityv of the vehicle is assumed to be constant with the consequence that the road vehicle system (26) and (27) becomes linear and can be solved, e.g., by means of classical spectral methods [24]. Subsequently, the obtained car spectra are integrated over all spectral frequencies in order to calculate all mean squares and covariances needed in the drift equatioṅ In the special caseV t = 0, the expectation E (F t ) of the driving force is subsequently calculated to the form where σ z = σ u are the standard deviations of the road level and slope process. Both standard deviations are calculable to σ u = σ/ √ 4δ by means of the linear filter equations, noted in Eq. (26).
In Fig. 9, the expected value (29) is plotted by a red line versus the related car velocity ν =v /ω 1 for the road intensity σ 2 = 1, the damping measure D = 0.15 and the bandwidth δ = 0.3. The  Fig. 9 shows probability distributions p ( f ) of the related driving process F t /c marked by blue lines and plotted for four speeds of the vehicle. The obtained densities have the exponential form of Gaussian product processes with extremely wide distributions indicating that the assumption of constant speed is physically not realizable.
Hence, it is more realistic to investigate the inverse problem that the velocity is fluctuating meanwhile the driving force is assumed to be constant or given by any characteristic depending on speed. This dynamic problem is numerically investigated by means of the nonlinear first-order system which is derived from Eq. (26), (27) and (28) introducing the dimensionless state processes (Z t , U t , Y t , X t ) = σ u Z t ,Ū t ,Ȳ t ,X t andF t = F t /c andV t = V t /c as well as the dimensionless time and noise by means of dτ = ω 1 dt and dW τ = √ ω 1 dW t , respec-  Fig. 10, the calculated speeds are narrowly distributed with one peak, only. For decreasing driving force, however, the velocity distribution becomes broader and bifurcates into bi-modal probability density with two peaks when the driving force approaches its critical value f /c = 1 where the vertical car vibrations become resonant. In this case, the global mean value is split into two separated mean values; i.e., the car velocity is fluctuating around two means. Simultaneously, the velocity is migrating from one mean to the next one and back again. For further decreasing driving forces into the under-critical range, the velocity regains its original distribution with one peak and one mean, only.

Probability bifurcation into turbulent travel speeds
The control problem, where the mean value of the driving force is calculated in Eq. (29), is now inverted in order to investigate the associated dynamic problem of fluctuating speeds when the driving force is constant. For these purposes, the velocity V t is centralized by The insertion of the mean-free velocity into the drift equation (28) gives Taking into account that E (Z t U t ) = 0 in the stationary case and that the velocityV t is uncorrelated to the road surface processes Z t and U t , the above differential equation leads to the stationary mean equation where E U 2 t = σ 2 u and the mean speed is introduced by m v = E(V t ) /ω 1 . The stationary drift equation (33) coincides with Eq. (8) when the dimensionless processes Ȳ t ,X t = (Y t , X t ) and Z t ,Ū t = (Z t U t )/z 0 are inserted. The covariances E (X t U t ) and E (Y t U t ) in Eq. (33) are determined by the vector equation is the system vector with four covariances and P t (m v ) is the driving vector given by (0, 0, P 3t , P 4t ) and the components t . The application of the expectation operator to Eq. (34) leads to the stationary covariances equation where the driving vector is calculated by means of E (Z t U t ) = 0 and E Z 2 t = E U 2 t = σ 2 u . The system matrices A and B in Eq. (35) are skew-symmetric with 4 × 4 elements. They are given, as follows: Equation (35) couples the covariances vector Q t with higher potencies of speedV t so that a sequence of moment's equations is obtained. It is approximately solved by means of the Gaussian closure [24] that E V t Q t = 0, when all zero mean processes are normally distributed Therewith, Eqs. (33) and (35) are solved by This result is nonlinear with respect to the unknown mean value m v = E(V t ) /ω 1 of the travel velocity. Note that Eq. (36) coincides with Eq. (12) of the deterministic wavy road when the bandwidth of the stochastic road process tends to zero (δ → 0) and its standard deviation is replaced by σ u = z 0 / √ 2. In Fig. 11, Eq. (36) is evaluated for the road bandwidth δ = 0.1 and the intensity σ u = 1. For driving forces f /c above or below the yellow range, Eq. (36) leads to one mean value m v and to velocity distributions with one maximum, only. They are approximated, e.g., by the density where C is the normalization constant and v is the density variable of the velocity processV t = V t /c with mean value m ν and standard deviation σ v For n = 2, the density p (v) in Eq. (37) is Gaussian. For n > 2, the density is more concentrated with shorter distribution tales and gives a better fitting to the measured densities, shown in the upper part of Fig. 10. For driving forces inside the yellow range in Fig. 11, Eq. 36 leads to three different mean values indicating that the velocity distribution possesses three extreme values with two probability concentrations approximated, e.g., by the non-normal density where a, b, c are three parameters of the cubic exponent to be determined by means of the three mean values m ν such that the calculated distribution density possesses two maxima and one minimum. Subsequently, the standard deviation σ v and the normalization C have to be determined. Since the stationary density is measurable applying the relative frequency and interval probability to all velocity trajectories realized in the whole speed range |V t | < ∞, it is physically not possible to introduce three different distributions, as, e.g., in the deterministic case. Hence, there is only one density of one stationary velocity process distributed over the whole speed space.
The obtained distribution density has the consequence that there is a permanent migration from the lower mean velocity to the upper one and back. The superposition of noisy velocity fluctuations and systematic migrations of vehicle velocities is typical for turbulences similar as in fluid dynamics. The velocity range in yellow in Fig. 11 is extended to the much bigger one marked in green when the damping decreases from D = 0.1 to D = 0.05. For almost vanishing damping marked by blue dashed lines, turbulences practically are no longer avoidable. In contrast to that, the yellow range of turbulences completely vanishes for stronger damping, e.g., for D = 0.4 where the cyanblack curve is obtained instead of the red-black curve. This new force velocity curve rises steadily without resonance since there is no longer a horizontal gradient. In the case of highly damped vehicles, high energy loss occurs in the dampers and one needs growing driving forces to maintain higher speeds of vehicles. Figure 12 shows resonance and travel velocities of vehicles for small car damping and growing road bandwidth. In this case, the resonance peak is lowered as well. However, the travel speed remains almost unchanged in the over-critical speed range. In the limiting case of white noise that the road bandwidth tends to infinity, the force-speed characteristic becomes linear marked by the thick black line in Fig. 12. This has the consequence that the turbulence speed region in yellow disappears, and Eq. 36 is simplified to the linear form ( f /c)/ (σ u ) 2 = 2Dm v .

Conclusions
The paper investigates strong nonlinear systems in longitudinal road vehicle dynamics. When the assumption of constant speed is removed, nonlinearities are introduced by multiplicative state variables. They lead to the effect that vehicle get stuck before resonance when driving on wavy road surfaces. To overcome this blockade, one has to speed up if the driving motor is strong enough.
Experiences on non-fixed wavy roads, e.g., dessert roads show that resonance speeds of vehicles are about 40 km/h due to plastic deformations by heavy vehicles of road maintenance. Turbulences are avoided by speeding up to over-critical speeds or by increas-ing the car damping. Both leads to stronger driving forces which are necessary to balance the energy loss in dampers.
Bifurcations in probability are explained by the fact that stationary density distributions of interest possess three extreme values: two maxima and one minimum indicating that solution trajectories are seldom in the minimum and often in the maxima with a permanent change between both. This is completely different to the linear case where the trajectories are Gaussian and distributed with one maximum, only.