Transition to chaos and escape phenomenon in two-degrees-of-freedom oscillator with a kinematic excitation

We study the dynamics of a two-degreesof-freedom (two-DOF) nonlinear oscillator representing a quarter-car model excited by a road roughness profile. Modeling the road profile by means of a harmonic function, we derive the Melnikov criterion for a system transition to chaos or escape. The analytically obtained estimations are confirmed by numerical simulations. To analyze the transient vibrations, we used recurrences.

strategy, the damping coefficient switched between two different values. Gobbi and Mastinu [2] considered a two-DOF model to derive a number of analytical formulae describing the dynamic behavior of passively suspended vehicles running on randomly profiled roads. Their linear model approach was generalized by Von Wagner [3], who determined a corresponding high-dimensional probability density by solving the Fokker-Planck equations. Finally, in papers [4][5][6] dynamics, bifurcations, and appearance of chaotic solutions were discussed.
However, the main issues of vehicle dynamics studies were unwanted and harmful vibration responses generated by vehicle as an effect of a rough surface road profile kinematic excitation [1,4,[7][8][9][10][11]. Thus, the efficient reduction of them is still a subject of research among automotive manufacturers and research groups [14,15]. Turkay and Akcay [7] considered constraints on the transfer functions from the road disturbance to the vertical acceleration, the suspension travel, and the tire deflection are derived for a quarter-car active suspension system using the vertical acceleration and/or the suspension travel measurements for feedback. The recent experimental and theoretical studies involved many new applications of active and semi-active control procedures [12][13][14][15][16]. Consequently, previous mechanical quarter-car models [3,10,11,15] have been reexamined in the context of active damper applications. Dampers based on a magnetorheological fluid with typical hysteretic char- acteristics have a significant expectation for effective vibration damping in many applications [4,[17][18][19][20][21][22][23].
Recent efforts have been also focused on studies of the excitation of an automobile by a road surface profile with harmful noise components [3,24,25]. Noise-like chaotic vibrations, appearing due to the system nonlinearities, have been investigated in simplified single-DOF models [4,6,19]. These papers follow the rich literature on an escape phenomenon in the symmetric and nonsymmetric Duffing or Helmholtz potentials [26][27][28][29][30], where critical system parameters were determined. Note that the condition for escape could be applied as a criterion of fractality in basins of attraction and also as a transition to chaos [29]. On the other hand, more sophisticated models of vehicle dynamics, namely half-car and full-car models, in the context of nonlinear response including chaotic solutions have been studied by Zhu and Ishitobi and also by Wang et al. [31][32][33].
To study a transition to the chaotic region and the corresponding critical parameters in a low-dimensional dynamical system analytically, the Melnikov theory [34][35][36][37][38] is often advocated. The application of this approach to a simple quarter-car model has been recently proposed by Li et al. [4] and Litak et al. [19]. In the above papers, a single-DOF model was used because of its simplicity. Consequently, the analytic consideration included also multiple-scale analysis and harmonic balance [20,39]. The present paper is a continuation of the previous studies with an extension to a more realistic two-DOF model, which includes the sprung and unsprung masses (Fig. 1).
Our article is organized in four sections. After introduction in the present section (Sect. 1) we present the model and discuss possibility of a global homoclinic bifurcation (Sect. 2). By the reduction of dimension we introduce the foundations of the Melnikov approach to the model. In this section (Sect. 2) we obtain the principal result as a critical curve which defines the system parameter regions of regular and nonperiodic (chaos, transient chaos, or escape) behavior. The simulation results illustrating those transition results are also shown. Further results including recurrence analysis confirming the theoretical predictions are provided in Sect. 3. Finally, in Sect. 4 we end up with conclusions and final remarks.

The model and global bifurcations
We start the analysis from the two-DOF model presented in Fig. 1. The dynamics of vehicle excited by the road profile quarter-car model is governed by nonlinear suspension characteristics. To examine a transition to the chaotic regime of vibrations, the Melnikov theory has been recently proposed [4,19]. In this perturbation approach the authors of previous works used a single-DOF model.
In the present note we go beyond this assumption by considering an extension of a vehicle model with the defined unsprung and sprung masses (Fig. 1). The differential equations of motion for both masses have the following form: where m 1 and m 2 denote the corresponding sprung and unsprung masses (Fig. 1), and x 0 = a cos ωt describes the harmonic corrugation of a road profile. c i , α i (for i = 1, 2) and β 1 are damping and stiffness coefficients, respectively. Now we define the new variable y = x 1 − x 2 of the relative motion: To adjust the above equations to application of the higher-dimensional Melnikov approach [40], the above equations can be approximated by introducing a small parameter : Consequently, the second equation (Eq. (2)) gets the new form Using the Melnikov approach [34] in y coordinate equation (Eq. (3)), we study the Hamiltonian system of the nodal kinetic energy defined at the saddle points. In this way, time variability of y andẏ is negligible. The above equations can be rewritten in the dimensionless form: where For further consideration, we assumed that the system parameters are: The excitation frequency and the corresponding amplitudes used in the analysis have been fixed to ω = 1.5, a = 0.08, and a = 0.12. In the numerical simulations we used the sampling time δt = 0.00418.
To analyze a homoclinic bifurcation in the sprung mass vibration, we make further approximation. The role of the small parameter is to determine the heteroclinic trajectory and decouple the equations of motion into separate equations for sprung and unsprung masses. Thus, after above normalizations the equations can be expressed aṡ Interestingly, in the limit of small limit, x 2 can be approximated as Note that in the above expression, x st , generated by slowly changing terms with y andẏ, is playing a role of the static displacement, while A is an amplitude. Thus, the unperturbed equations (for = 0) can be obtained from the gradient of the unperturbed Hamiltonian H 0 (y, v): where H 0 is defined as follows: The corresponding effective potential (Fig. 2a) is given by the expression Following the standard Melnikov theory [34][35][36][37][38], we get the heteroclinic orbits connecting the two hyperbolic saddle points (coinciding with the maxima of the potential V (y) (Eq. (12)): y = ± A −1 3 .
They can be expressed analytically as After perturbations of the heteroclinic orbits (Fig. 2b), the stable and unstable manifolds are calculated. Because of perturbations, they are detracted (Fig. 3). As the characteristic distance between them d → 0, the system possibilities of mixed solutions (regular and escape) appear, which is equivalent to irregular chaotic and chaotic transient solutions. Thus, d = 0 is the ideal criterion for chaos appearance.
Formally, the distance d between perturbed stable and unstable manifolds are proportional to the Melnikov integral (Fig. 3 [34][35][36], which can be written as where ∧ denotes a wedge product, the differential form h is the gradient of the unperturbed Hamiltonian, g is related to the perturbation part, and t 0 is an integration constant. Both forms are defined on unperturbed heteroclinic orbit stable and unstable manifolds W st(unst) = (y * st(unst) , v * st(unst) ): Thus, after substitution the above forms into the Melnikov function M(t 0 ) (Eq. (14)), we get: The condition for a global homoclinic transition, corresponding to a possible horse-shoe-type crosssection of stable and unstable manifolds, can be written as Consequently, the critical parameter η = a/C 1 is now The results of the above analysis are presented in Fig. 4a. The curve separates the region of regular solutions from the chaotic and escape ones. The black where the two indicated points symbolize the parameters used for numerical calculation points above and below the critical curve for a = 0.08 at the point "1", and a = 0.12115 at the point "2" (ω = 1.5 for both cases)  Note that Fig. 5a shows the mono-frequency, while Fig. 5b shows more complex responses (multifrequency or irregular). This is clearly visible in Figs. 6a and 6b, where we show the corresponding phase diagrams. In Fig. 6a one can see a single line, while in Fig. 6b lines are split into characteristic threeline patterns. It is worth mentioning that the Melnikov criterion (Eq. (18) and Fig. 4) specifies the global transition associated with the destruction of borders between basins of attractions belonging to different solutions. In our case, one of basins is related to an escape from the potential well. The irregular solution denoted as no. 2 in Fig. 4 (see also Figs. 5b and 6b) must be related to such an escape. Indeed, continuing the calculations for long enough time interval, one observes the escape in the plot on time series (Fig. 7a) and on the phase portrait (Fig. 7b), respectively.
Focusing on the above solutions, we discuss the recurrence properties of numerical results with more details in the next section.

Recurrence plot analysis
The numerical solutions of the regular and transient nature can be analyzed more carefully by the recurrence plots [41,42]. This method is based on the statistics of recurrences and can be described by the matrix form R m,ε with corresponding 0 and 1 elements: where x i and x j are usually defined in the embedding space of dimension m, and ε is the threshold value.
Here indices i and j denote the sampling instants. In our case we decided to use the m = 2 and the embedding space including the displacements of sprung and unsprung masses: x j = [y, x 2 ] (Eqs. (7-8) with = 1).
Having 0 and 1 values, they are translated into the recurrence diagram as an empty place and colored dots, respectively (see Fig. 9). Here w denotes the Theiler window used to exclude identical and neighboring points from the analysis [44]. Webber and Zbilut [45] and later Marwan and collaborators [43,44] developed the recurrence quantification analysis (RQA) for recurrence plots.
Shortly after its invention, RQA was addressed to the biological and physiologic systems [45,46]. Recently this method has been applied for several technical systems [47][48][49]. The first parameter of the RQA Fig. 8 The idea of states summation idea in the embedded space and the sphere of radius ε defining the correlation function is the recurrence rate RR, which calculates the number of recurrences. In our case, the Theiler window w = 0 in order to get the consistency with correlation sum [44]. In this language, ε expresses the correlation length of the characteristic sphere radius in the embedded space (Fig. 8). Note that the correlation sum is an important tool which could be used to derive correlation dimension D 2 . For small enough ε, where ε 0 is the arbitrary length. Furthermore, the RQA can be used to identify topological structures of diagonal and vertical lines. In its frame, RQA provides us with the probability p(l) or p(v) of line distribution according to their lengths l or v (for diagonal and vertical lines). Practically they are calculated as where x = l or v, depending on diagonal or vertical structures in the specific recurrence diagram. P ε (x) denotes the unnormalized probability for a given threshold value ε. In this way, the Shannon information entropies (L ENTR ) can be defined for diagonal line collections: Other properties, the determinism DET and laminarity LAM, are defined as Fig. 9 Recurrence plots for regular (a) and irregular (b) solutions for the same system parameters as point "1" and "2" in Figs. 4a and 4b, respectively (ε = 0.01) where l min and v min denote the minimal values which should be chosen for a specific dynamical system. In our calculations we have assumed that l min = v min = 2.
The determinism DET is a measure of the predictability of the examined time series and gives the ratio of recurrent points formed in diagonals to all recurrent points. Note that in a periodic system all points should be included in the lines. On the other hand, the laminarity LAM is a similar measure which corresponds to points formed in vertical lines. This measure indicates the dynamics behind sampling point changes.
The results of our analysis calculated for time series "1" and "2" (see Figs. 4, 5 and 6) are presented in Fig. 9, where we present the results of RP for ε = 0.01. Note that both plots show regular features. However, Fig. 9a is composed of full diagonals lines only, while in Fig. 9b each fifth line is full, and between them one observes short line pieces or even insulated points. This would indicate a multi-frequency modulated solution. There can be also a transient with basic regular and superimposed chaotic solutions. To shine this difference with more lights, we show some estimated RQA parameters in Table 1. Obviously, RR is smaller for "2" (Fig. 9b) as we have broken lines instead of full ones (Fig. 9a). Moreover, the determinism and laminarity (DET and LAM) are smaller, telling us that the system is less regular. Consequently the more peculiar distribution of line lengths is confirmed by L ENTR , which is larger for the solution "2". Additionally, in Fig. 10 we present the results of RR versus ε (for relatively small ε). One can see a significant difference between both solutions. Fig. 10 Recurrence rate RR versus threshold value ε for regular "1" and irregular "2" solutions ("1" and "2" as in Figs. 4a and 4b, respectively)

Summary and conclusions
We have analyzed the two-DOF quarter-car model, assuming that damping and suspension through the unsprung mass excited by the road profile corrugation can act as a perturbation on the main sprung mass. The obtained Melnikov criterion was latter confirmed by numerical simulations. The main conclusion coming from that point would be loss of stability of the system appearing as the chaotic or transient chaotic or escape solution. The present investigation is going beyond research dealing with a single-DOF quarter-car model [4,18,19,24]. In particular, the single-DOF model assumes that the unsprung mass is significantly smaller than the sprung mass.
One should note that the model used in this paper, although more realistic than the previous singledegree-of-freedom ones, is relatively simple and would not be sufficient to simulate a detailed response of a vehicle or compare to experimental results from real vehicles. Unfortunately, more sophisticated half-car and full-car models [31][32][33] cannot be used in the frame of the presented approach as the heteroclinc trajectories could not be defined reliably. Furthermore, in higher-DOF systems the analytic perturbation calculations are not possible.
However, the present quarter-car model is able to capture the major nonlinear effects that occur in vehicle dynamics and has demonstrated the transition to chaotic vibrations and synchronization phenomena [25,32,39]. Interestingly, the resulting critical amplitude curve (Fig. 4a, Eq. (19)) has the maximum for ω ≈ 3.1 and the minimum for ω ≈ 2.2. The minimum is obviously related to the resonance region of the decoupled unsprung solution x 2 (Eqs. (8) and (9)).
It should be also noted that the recurrence plot technique appeared to be very useful to study the transient signals. Thus, the conclusions that came from the analytic approach have been confirmed. Indeed, this method is designed for the short time series [44,49] Interestingly, it also works for nonstationary signals [44,48]. The recurrences for RP and RQ analyses have been obtain using the available command line code written by Marwan [50]. The recurrence approach can be also used to higher-DOF models of vehicle dynamics. The corresponding results on half-and full-vehicle models will reported in a separate paper.