Bifurcation phenomena and statistical regularities in dynamics of forced impacting oscillator

The paper is devoted to the study of harmonically forced impacting oscillator. The physical model for oscillator is a cart on a guide connected to the support with springs and excited by the stepper motor. The support also is provided with limiter of motion. The mathematical model for this system is defined with the second-order piecewise smooth differential equation. Model’s nonlinearity is connected with the incorporation of dry friction and generalized Hertz contact law. Analyzing the classical Poincare sections and inter-impact sequences obtained experimentally and numerically, the bifurcations and statistical properties of periodic, multi-periodic, and chaotic regimes were examined. The development of impact-adding regime as a new nonlinear phenomenon when the forcing frequency varies was observed.


Introduction
A large number of artificial devices and natural systems contain oscillating elements, dynamics of which essentially affect [1,2] the parent system. Therefore, one of the important problems is the classification of the dynamics types of oscillating systems, sources of their appearance, evolution scenario, etc. To construct vibration systems, springs and dash-pots as structural elements are used. To model their dynamics, the smooth mathematical models are applied, for studies of which the powerful analytical and numerical tools have been developed. But there are a wide range of observed phenomena [3][4][5][6][7] which cannot be described within the framework of smooth models.
Incorporation of dry friction and impact rheologic elements into vibration systems allows one to describe a wider range of physical processes. Corresponding mathematical models are governed by the nonsmooth equations of motions [5,[8][9][10]. The detailed analytical treatment is admitted mainly by piecewise linear models. If nonlinearity is taken into account, the nonsmooth models exhibit essentially nonlocal dynamics [11][12][13][14][15] for the study of which numerical methods are involved. The attractors inherent in the smooth models can be endowed by new properties. Moreover, such models can demonstrate special types of evolution including border-collision bifurcations [10,[15][16][17][18] and chattering phenomena [19].
The present studies are focused on revealing the peculiarities of nonlinear regimes occurred in nons-  [20] mooth models. In particular, it is investigated a onedegree-of-freedom oscillator incorporating dry and viscous resistance, the generalized Hertz contact law [14,20] during an impact and subjected to the harmonic loading. In [20], the multi-parametric mathematical model for the impact oscillator was developed. During the model validation, several sets of model parameters providing the best agreement of theoretical and experimental data were obtained. Validated mathematical submodels allowed one to estimate impact duration time and coefficient of restitution, and to classify roughly the vibration regimes based on the bifurcation diagrams. As shown in [21], additional information on models' regimes can be obtained from the inter-impact sequence analysis which was performed for the periodic and chaotic attractors.
The investigations mentioned in [20,21] do not cover the whole interval of forcing frequency variation, although a cursory glance at the bifurcation diagram indicates that no less interesting regimes are expected to be revealed in another part of the diagram.
Thus, the question arises: What else new regimes and their bifurcations can be caught in selected interval of forcing frequency? To answer the question, the classical Poincaré section technique, top Lyapunov exponents (TLE), inter-impact (spike) sequences, and information entropy are used.
The detailed description of the studies is provided below. Thus, the next section contains the experimental rig introduction and corresponding equation of motion. Section 3 is devoted to the presentation of experimental results when the frequency of forcing varies. Sections 4 and 5 are concerned with the numerical investigation of mathematical model and identification of observed regimes. The final section ends up with the concluding remarks.

Experimental stand and mathematical model
To begin with, let us give brief description of the experimental stand and proper mathematical model (the advance discussion can be found in [20]). In Fig. 1a, there is presented part of the larger stand used in the present work as an experimental model for investigations of dynamics of one-degree-of-freedom harmonically excited oscillator with one-sided obstacle [20]. It consists of a cart (1) on a guide (4) with the integrated linear ball bearing and Hall effect sensor (2). The cart is connected to the support (7) with the springs (3). A stepper motor (9) is mounted on the cart with a disk (10) and an unbalance (5) in order to realize harmonic forcing. Impacts are realized by limiters of motion in the form of stops (6) mounted on the support and the cart. A transoptor (8) is situated on the cart behind the disk in order to control absolute angular position of the stepper motor. Position of the cart is measured by the Hall effect sensor, while angular position of the disk is obtained from an encoder integrated with the stepper motor. Measurement data are then registered by National Instrument data acquisition system. Figure 1b exhibits physical concept of the experimental rig. As was shown in the work [20], it is described by the following differential equation of motion where In these equations, x denotes position of the cart; m-total mass of the cart; k-total stiffness of the springs; F R -resistance force acting on the cart; F I (x,ẋ)-impact force; m 0 -an unbalanced mass; eradius of the unbalanced mass; ω-angular speed of the disk; c-linear damping coefficient; T -magnitude of rolling resistance (assumed to not depend on the normal load); ε-special parameter of regularization of the sign function (ẋ/ √ẋ 2 + ε 2 ≈ sgnẋ); x I -position of the obstacle; k I -stiffness coefficient of the obstacle; b-damping coefficient of the obstacle; n 1 , n 2 , n 3parameters of the contact. Note that impact is modeled as compliant based on extended Hertzian contact stiffness and damping.
To achieve good agreement between experimental and simulated data, the estimation of the parameters based on the numerical optimization was carried out [20]. The procedure consists in minimization of objective function where x exp and x sim are the cart positions obtained experimentally and during numerical simulation on the time interval [t 1 , t 2 ], while μ = (k, c, T, k I , b, n 2,3 ) denotes the vector of the estimated parameters. Minimization of the function F 0 is provided by the Nelder-Mead multi-dimensional optimization algorithm realized with the help of the f minsearch MATLAB function. Using the parameters measured directly m = 8.735 kg, m 0 e = 0.01805 kg · m and fixed quantities ε = 10 −6 m/s and n 1 = 3/2, the application of F 0 minimization procedure results in several groups of parameters defining different submodels, concerning different mathematical descriptions of resistance in the linear bearing and various versions of compliant impact model, including different models of the related nonlinear damping. From these groups, the optimal (in the sense of simplicity and correspondence between experimental and simulation data) submodel (1) was chosen with rather small F 0 = 0.00997 mm 2 , for which the set of estimated parameters is as follows: k = 1418.9 N/m, c = 6.6511 N · s/m, T = 0.63133 N, k I = 2.3983 · 10 8 N/m 3/2 , b = 0.8485 m −n 3 s n 3 , n 2 = 3/2, n 3 = 0.18667. Note that small deviations from these derived values or certain small modifications of model of nonlinear damping during the impact (compare for example the models denoted as BC and BD, Table 2, [20]) do not affect essentially the results of simulations. In the reference [20], the correspondence between the mathematical model (1) and real behavior of the system for its different operation regimes is also investigated, including free and forced bifurcation vibrations for a wider range of forcing frequencies and different positions of the obstacle. In all cases, a good agreement between experiment and simulation is observed.
Further analysis presented in the paper corresponds to the obstacle position x I = −2.086 · 10 −3 m. Since its value is negative, at the starting moment the cart is located in the permanent contact with the obstacle. Therefore, considerations always deal with the two component equation (1) and do not encounter the purely linear oscillations. This essentially complicates the study of model when, for instance, the first impact one is interested in.

Experimental and preliminary simulation results
Let us briefly recall that the key point of paper [20] is the comprehensive analysis of admissible models accompanied by the comparison of bifurcation diagrams without studies of their fine structure. Therefore, it was constructed the bifurcation diagram for ω ∈ [0, 50] rad/s containing the interval ω ∈ [35, 50] rad/s where the classical period-doubling scenario leading to a chaotic attractor creation is observed. To find out the properties of these chaotic solutions, the inter-impact sequences were constructed numerically and their statistical properties were analyzed in [21]. But, besides the chaotic "window" at rather large ω, the bifurcation diagram contains another chaotic zone ω < 25 rad/s with oscillations whose amplitudes are five times smaller.
It turns out that at these values of ω the scenario of oscillations development is much more complicated and richer. In the present research, to consider the properties of the regimes, in addition to the classical Poincare section technique the new experimental data concerning the inter-impact events are extracted and compared with the results of simulation of system (1). Moreover, scaling features of bifurcational diagrams, the TLE and entropy index behavior and new impact-adding phenomenon are investigated. Estimation of the moment of first grazing bifurcation based on the approximate solution is discussed as well.
Let us consider each step of the studies in detail and start from the construction of classical Poincare sections, based on the points at time moments t = 2πi/ω, i = 1, 2, . . .. The time moments of impacts also are derived. These moments define the inter-impact (interspike) intervals or inter-impact sequence T j [22]. This sequence can be regarded as a sort of Poincare section [23,24].
Thus, Fig. 2 exhibits experimental bifurcation diagrams for increasing value of excitation frequency ω constructed on the base of classical Poincare sections (Fig. 2a) and inter-impact sequences (Fig. 2b). They are built for the continuous change of bifurcation parameter with rate 8.0905 × 10 −3 rad/s 2 . Since the impacts are very weak, a threshold below the obstacle is intro-duced and impacts are detected during crossing the level x I − 0.01 mm. Moreover, the experimental data are slightly smoothed by passing them through a filter of transfer function (θ s + 1) −1 , where θ = 10 −3 s. To clarify the observed regimes and prove physical and mathematical models consistency, the analogical bifurcation diagrams (Fig. 2c, d) are constructed numerically. Since experimental data are blurred because of certain additional oscillations not taken into account in the model and noise in the measurement system, it is not possible to observe some small details on experimental bifurcation diagrams. However, one can notice good general agreement between corresponding bifurcation diagrams.
Among others, the sequence of impact-adding events for frequencies ω ∈ [14; 17.5] rad/s is observed, both experimentally and numerically, while period of the attractor remains equal to the period of external forcing. It is the scenario related to transition between attractors corresponding to permanent contact between the cart and obstacle and solutions with impacts [25]. This case will be investigated in more detail in Sect. 5. It should be, however, noted that in the case of experimental bifurcation scenario, the impact-adding sequence is confirmed at least up to two-impact periodic trajec-tory. For higher frequencies, there are observed jumps between the extensions of the branches corresponding to the long inter-impact intervals (upper parts of bifurcation diagrams in Fig. 2c, d). It can be explained by appearance and loss of some impacts induced by some unpredictable factors present in the real system.
Analyzing Fig. 2, it is experimentally confirmed the scenario of impact-adding events detected later numerically, at least up to two-impact periodic trajectory. It should be, however, noted that in the case of experimental bifurcation diagram, there are observed jumps between the extensions of the branches corresponding to the long inter-impact intervals (upper part of bifurcation diagrams), for the periodic trajectories with more than two impacts. It can be explained by appearance and loss of some impacts induced by some unpredictable factors present in the real system.
When ω reaches 20 rad/s, the system goes out to the attractor with essentially different amplitude. At increasing ω, the system undergoes several grazing bifurcations which are identified clearly in the interimpact sequence diagram (Fig. 2c, d). There is a critical value of ω corresponding to the multi-periodic regime degeneration. Instead, a periodic trajectory describing the single-impact regime of oscillations develops.

Bifurcation analysis of the mathematical model for the oscillator with impacts
Let us now investigate the phenomena admitted by the mathematical model (1) in more detail. Figure 3 exhibits bifurcation diagrams presenting the system dynamics for frequencies of the external forcing lower than 25.0 rad/s. They are obtained based on numerical simulations for increasing (red points) and decreasing (green points) value of bifurcation parameter ω. Bifurcation diagrams consist of Poincare maps defined in standard way, i.e., they are based on sampling the position of the cart every period of external forcing (Fig. 3a) and inter-impact intervals (Fig. 3b).
Analysis of these diagrams (Fig. 3) for ω ∈ (20.5; 25) shows that several grazing bifurcations and attractor coexistence (ω ≈ 22.71) take place. Visual figure examination suggests the presence of some sort of self-similar diagram structure. In particular, the cascade of grazing bifurcations can manifest the scaling properties. To consider this, let us estimate the fre- No less interesting is the bifurcation diagram at ω ∈ [17.5; 20.5] (Fig. 4). The classical period-doubling cascade exists in the window ω ∈ (20; 20.5). It is also revealed the intervals near forcing frequencies which are equal to 17.9, 18.5, and 19.7 rad/s with the hysteretic behavior of bifurcation diagrams. This can imply the coexistence of at least two attractors.
Note that at the diagram in Fig. 4b several similar objects can be distinguished. They lay in the intervals Top Lyapunov exponents To consider the trajectory stability, the TLE [14,26] are derived. In Fig. 5, the dependence of the TLE λ 1 on the frequency ω is presented. One corresponding to the excitation phase, which is always equal to zero, is omitted. Period of normalization was chosen as equal to the period of external forcing. The presented bifurcation diagrams are performed for increasing (red line) and decreasing (green line) values of bifurcation parameter with a step equal to ±0.005 rad/s. 100 and 50 periods of external forcing were ignored as transient motion after starting simulation and after each further jump of bifurcation parameter, respectively. For each value of external forcing, the TLE is computed for 10 3 periods. In the diagram of Fig. 5 Fig. 5 The TLE dependence on the frequency ω ysis, it should be note that the inter-spike sequences can provide the possibility to extract the TLE estimation [27,28]. It shows promise in addressing the direct TLE measurement from impact sequences during experiments.
Statistical analysis Some additional information on the observed regimes can be obtained from the statistical analysis of interspike sequences [23,29]. The sequence T j at fixed ω is considered as a sample arranged in the normalized histogram. Such a sample allows one to estimate the probability distribution. This distribution is characterized within the framework of the Tsallis entropy approach [30,31]. This generalized entropy S T is defined as follows where the parameter q is the entropy index, N is the total number of microstates, k is Boltzmann's constant (k = 1 can be assumed). The value of q allows one to characterize the structure of distribution. Namely, when q < 1 the rare events are important, whereas at q > 1 the frequent events are more weight [31]. Therefore, the estimation of q is required. This can be done on the base of maximum entropy principle [32]. It can be proved that the entropy S T gains its maximum S T max at the equiprobable microstates when p i = 1/N .
The deviation of S T from the S T max can be characterized by the redundancy function R T = 1 − S T /S T max . This function depends on the entropy index q. The derivation of R T maximum with respect to q provides Fig. 6 The normalized histogram and the redundancy function R T . The R T maximum is achieved at q = 0.62 Fig. 7 The dependence of the entropy index q on the frequency ω the estimation of entropy index. Thus, consider the q estimation at fixed ω = 20.465 and arrange the sequence T j in the normalized histogram (Fig. 6a) at the interval (0.05; 0.3) with the step 0.001. Some zeroes in the frequencies sequence should be eliminated. Using resulting distribution, p i , N , and finally the redundancy function R T can be identified. The profile of R T is depicted in Fig. 6b. The numerical procedure of "Mathematica" packages gives the maximum estimation q = 0.62 < 1. This means that the distribution describes a large number of rare events which affect the system behavior. Doing in the same manner, the dependence of q on the frequency ω is studied. The ω interval (20.433; 20.515), when the chaotic regimes are developed, is chosen. Estimating the entropy indexes q with the step 0.005, the graph in Fig. 7 is plotted. From the analysis of this graph, it follows that chaotic attractors are characterized by q < 1 which are changed weakly. At the edges of frequency interval, when the chaotic attractors are reduced to the almost periodic trajectories, the entropy index grows abruptly to 1.
This means that the distribution contains fewer variants which occur more frequently.

The system behavior at low frequencies
The regimes created at low frequencies deserve to the more detailed analysis. For the best of our knowledge, this type of dynamics is not widely known [25]. It is worth noting that the obstacle grazing at frequencies from the interval ω ∈ [14; 18] does not cause the trajectory bifurcation as opposed to the phenomena at ω ∈ [20.5; 23]. Let us construct the phase portraits of typical orbits forming the so-called impact-adding scenario (do not confuse with the grazing-induced impact-adding bifurcation [10]) presented in Fig. 2d.
For ω = 14 rad/s, one can observe a periodic attractor (Fig. 8) corresponding to the cart being in permanent contact with the stop. Position of the obstacle is denoted by red line, while red dot corresponds to Poincare section t = 2πi/ω, i = 1, 2, . . . .
For ω = 14.1 rad/s, a phase portrait with a segment of trajectory detached from the obstacle and one impact is seen. The intersections of increasing position x of the cart with the level of obstacle x I (instances of impacts) are denoted by blue dots. Further increase in the forcing frequency leads to separation from the obstacle of further segments of the orbit and adding the consecutive impacts. For ω = 17.5 rad/s, one can see an orbit with five impacts and of period equal to the period of external excitation. Finally, an irregular orbit and a periodic stable solution of double period with five impacts for forcing frequency equal to 18 and 19 rad/s, respectively, are shown.
As it is shown in the first phase portrait of Fig. 8 at ω = 14 rad/s, the oscillations exist in the region, = m 0 eω 2 sin ωt.
Since k k I , then the nonlinear term prevails and equation (2) is closer to the strongly nonlinear equation. The profile (Fig. 9a) of periodic solution testifies, furthermore, about the strongly nonlinear character of oscillations.
In order to understand the nature of this regime, let us consider the Fourier spectrum of solution (Fig. 9b). It is seen that the signal contains only two dominant frequencies, i.e., ω and 2ω. Therefore, the main features of solution can be approximated by the following expression Using the numerical procedure from the "Mathematica" packages, the solution (3) can be derived in the form x = −2.0813 × 10 −3 − 1.3136 × 10 −6 cos ωt + 3.9328×10 −6 sin ωt+6.6538×10 −7 cos 2ωt+2.915× 10 −7 sin 2ωt.
But due to the presence of power and discontinuous functions, one cannot get the analytical expressions for the required coefficients of Fourier series. Instead, the first harmonics can be treated.
Consider the linearized counterpart [33] of equation (2). This equation can be considered in the vicinity of steady state x = X 0 = const which satisfies the algebraic equation k X 0 + k I (X 0 − x I ) n 1 = 0. From this, it follows X 0 = −2.0807 × 10 −3 .
Next let us suppose that the dissipation is small and consider the deviation from the steady state y = x +X 0 . Equation (2) reduces to the following one where is a small parameter prescribing formally the smallness of nonlinear friction. Then, the simple linear approximation gives the solution x = X 0 + A 1 sin ωt + A 2 cos ωt, where Thus, the solution at the parameters fixed above is This form of solution allows one to estimate the value of frequency at which grazing the x = x I is possible. Namely, let us derive ω when x = x I , i.e., From this, it follows that ω = 15.657. This value is overestimated in comparison with the exact value of first grazing.
To take into account the equation structure under manipulations more effectively, another method of approximation for the y amplitude is applied. Here, the approach used in [34,35] is referred. If the viscosity takes into account, then the approximate solution of Eq. (2) with the external force F = F 0 sin ωt reads as follows y = A sin(ωt − α). If instead the solution has the form y = A sin ωt, then the force can be written in the form F = F 0 sin(ωt + α). Since the equation holds for all ωt, and hence for ωt = π/2. But in this case y = A, dy/dt = 0, d 2 y/dt 2 = −Aω 2 , F = F 0 cos α.
The comparison of "exact" solution of system (2), numerically derived solution (3) and (6) at ω = 14, is presented in Fig. 9a. From the analysis of this figure, it can be concluded that the last method gives better results than the method of simple linearization. Using the solution (5), the estimation of grazing phenomenon is possible. Solving the system x I = X 0 + A(ω), one gets ω = 12.789. It is evident that this value is undervalued with respect to exact value of the first grazing. Analyzing Eq. (2), it is worth noting that the deviation from the fixed point X 0 is not small enough. Therefore, the omitting of high-order terms, appearing in the Taylor series for the power function (y + X 0 − x I ) n 1 expanded near X 0 − x I , contributes a notable inaccuracy in the approximated equation.
Between the first grazing (ω ≈ 14 rad/s) and chaotic mode (ω ≈ 18 rad/s) oscillations manifest some interesting features. Analyzing the profiles of solution (Fig. 10), one observes that trajectories contain small oscillations near stopper position. This regime looks like an incomplete chattering [11,19,24]. The finite number of small oscillations and the length of shelves (Fig. 10a) are caused by dry friction incorporation.
It can be checked the Fourier spectrum contains four dominant frequencies jω, j = 1 . . . 4. At ω ≈ Fig. 9 The numerical solution (black points) of system (2), numerically approximated solution (3) (green curve), the profile of function (4) (blue curve) and the graph of function (6) (brown curve). The Fourier spectrum for the solution of system (2) (Fig. 10b). The number of excited frequencies in Fourier spectrum reaches dozen of commensurate values. Exploring the micro-oscillations behavior in the chaotic mode at ω = 18 rad/s (Fig. 10c), it should be noted that three most prevailing peaks remain almost unchanged during time increasing, whereas the smallest peak undergoes the essential deviations, although four hits in the chattering window are still observed. The Fourier spectrum possesses clearly visible several dominant frequencies, and its smooth part is also essentially excited. Since these regimes possess the commensurate frequencies in Fourier spectra, one can consider them as nonlinear quasiperiodic oscillations forming in the phase space closed orbit on the torus surface.

Concluding remarks
Summarizing aforementioned results, the research presents new experimental and theoretical studies of the physical oscillating system with impacts. The observed dynamics of this system is managed to play with the help of nonsmooth mathematical model which takes into account the dry and viscous friction, generalized Hertz contact law during impacts. It should be emphasized that the system is considered under the conditions when pure linear dynamics is absent and the most of analytical methods are not applicable. Therefore, the information on the system properties is mainly obtained from the numerical integration, the results of which are used in the qualitative and statistical analysis. Dealing with impacts, it was shown that the inter-spike sequences analysis has significant capability to experimental and theoretical investigation of impacting systems.
The main results of the work can be enumerated as follows: -In the considered interval of frequency ω, the system possesses the periodic, multi-periodic, and chaotic attractors. Among bifurcation scenarios, the period-doubling cascades leading to chaotic attractor creation and stiff attractor appearance and destruction are revealed. The windows with hysteretic behavior and period-adding cascades are discovered as well. These dynamical phenomena are endowed by scaling properties. -It is shown that in the system the impact-adding scenario is realized which does not show up in the classic Poincare diagram, whereas the inter-impact sequence technique allows us to identify main its features. -The inter-impact sequences admit analysis via the q-entropy concept. This approach tells us that the distribution of impacts characterizes the set of rare events. The parameter of this distribution, i.e., entropy index q, is sensitive to the approaching to a bifurcation point. Therefore, it can be regarded as a natural identificator of chaotic attractor variability. -the estimation of the moment of first grazing bifurcation is carried out with the help of approximate solution derivation. In particular, despite the satisfactory estimation of bifurcation value of frequency, the improvement can be achieved by the incorporation of additional terms in the approximation. However, then the derivation of bifurcation value becomes more complex task. On the other hand, to improve the approximation procedure, instead of linear generating system the nonlinear one can be chosen. This way of studies is very promising but can be realized for that parameter values when the generating system possesses solutions written in the closed analytic form.
Theoretical studies thus clarified the experimentally observed phenomena, shed the light on the specific statistical properties of chaotic regimes, allowed one to identify novel bifurcation events and open up new challenges in the field of nonlinear dynamics.