Vibration of the Duffing Oscillator: Effect of Fractional Damping

We have applied the Melnikov criterion to examine a global homoclinic bifurcation and transition to chaos in a case of the Duffing system with nonlinear fractional damping and external excitation. Using perturbation methods we have found a critical forcing amplitude above which the system may behave chaotically. The results have been verified by numerical simulations using standard nonlinear tools as Poincare maps and a Lyapunov exponent. Above the critical Melnikov amplitude $\mu_c$, which is the sufficient condition of a global homoclinic bifurcation, we have observed the region with a transient chaotic motion.


Introduction
A nonlinear oscillator with the Duffing term but linear damping is one of the simplest systems leading to chaotic motion studied originally by Ueda [1,2,3]. The problem of its nonlinear vibrations has attracted researchers from various fields of research; from natural science [4,5] and mathematics [6] to mechanical [8,9,10] and electrical engineering [1,2,3]. This system, for a negative linear part of stiffness, shows homoclinic orbits, and the transition to chaotic vibration can be treated analytically by the Melnikov method. Note that such a treatment has been already performed successfully to selected problems with various potentials [6,7]. In spite of fact that a large bibliography is devoted to vibration of a single realization [1][2][3][4][5][6][7][8][9][10][11][12] or set of coupled Duffing oscillators [13,14,15,16,17] with numerous modifications to potential and forcing parts, the problem of nonlinear damping in chaotically vibrating system has been focused shortly [12]. The fractionally damped Duffing's equation, for a single well potential with positive linear part, has been studied in detail by Padovan and Sawicki [18]. They were especially focused on the determination of the influence of fractional damping on the frequency amplitude response. Some additional insight into this problem can be also found in the context of self excitation effects [11,13,16,17,19]. On the other hand, in the paper of Trueba et al. [12], the systematic discussion on square and cubic damping effects on global homoclinic bifurcations in the Duffing system was given given while Awrejcewicz and Holicke applied the Melnikov theory to Duffing system in a dry friction limit [20]. In the present paper we reexamine the problem of nonlinear damping including fractional damping cases by both; the analytical Melnikov method as well as numerical simulations. Namely we are going to study a single degree of freedom nonlinear oscillator with the Duffing potential and fractional damping: where x is displacement andẋ velocity, respectively, while the external potential force F x and corresponding potential V (x) (Fig. 1a) are defined; Here the non-linear damping term is defined by the exponent p: We have plotted the above function versus velocity (v =ẋ) in Fig. 1b for few values of p. Note that, the case p → 0 (see p = 0.1 in Fig. 1b for a relatively small velocity) mimics the dry friction phenomenon [21,22].

Melnikov Analysis and Numerical Simulations
We start our analysis from the unperturbed Hamiltonian H 0 : The potential function V (x) (Fig. 1a) has the local peak at the saddle point x = 0. Existence of this point with a horizontal tangent makes possible homoclinic bifurcations to take place. This includes transitions from regular to chaotic solutions. To study effects of damping and excitation, we apply small perturbations around the homoclinic orbits. This implies using a small parameter ǫ in the Eq. 1 with perturbation terms. Equation 1 can be simultaneously uncoupled into two differential equations of the first order: where ǫα = α and ǫμ = µ, respectively. At the saddle point x = 0, for an unperturbed system (Fig. 1a), the system velocity reaches zero in velocity v = 0 (for infinite time t = ±∞) so the total energy has only its potential part. Transforming Eqs. 2,4 for a chosen nodal energy (E = 0) and for δ < 0, γ > 0 we get the following expression for velocity: Now one can perform integration over x: where t 0 represents an integration constant. Finally, we get so called homoclinic orbits ( Fig. 2): where '+' and '−' signs are related to left-and right-side orbits, respectively (Fig. 2). Note, the central saddle point x 0 = 0 is reached in time t corresponding to +∞ and −∞, respectively. In case of perturbed orbits W S (a stable manifold) and W U (an unstable manifold) the distance between them is given by the Melnikov function M(t 0 ): where the corresponding differential form h means the gradient of unperturbed Hamiltonian (Eq. 4): while g is a perturbation form (Eq. 5) to the same Hamiltonian: All differential forms are defined on homoclinic orbits (x, v) = (x * , v * ) (Eq. 8).
Condition for a global homoclinic transition, corresponding to a horse-shoe type of stable and unstable manifolds cross-section (Fig. 2), can be written as: Evaluating the corresponding integral (Eq. 9) after some lengthly algebra the last condition (Eq. 12) yields to a critical value of excitation amplitude [4,12] µ c : where the B(r, s) is the Euler Beta function: and Γ(n) denotes the Euler Gamma function: In Fig. 3a we plotted the results of Melnikov analysis for a critical amplitude µ c /α for few values of p while δ = −1, γ = 1. For µ > µ c the system can transit to chaotic vibrations. In the next figure (Fig. 3b) we plotted µ c versus p. Here we have also compared the results of analytical investigations represented by a solid line with numerical simulations represented by circle points. The numerical results can give a local criterion of chaotic vibration appearance, namely a positive value of the largest Lyapunov exponent λ 1 [23]. Note that in all cases ( This indicates transition to chaotic vibrations. Note that for p = 0.1 we have got some small but positive values of λ 1 for weak excitation (µ close to 0). When we examined this case with more details we observed that amplitude of vibrations was very small. Appearance of small positive value of λ 1 seams to be a combined effect of physical phenomenon of dry friction and our assumed numerical integration procedure of constant time steps ∆t(in most cases we used ∆t = 2π/(500ω)) . It is known from the studies of a stick-slip phenomenon [24,25] that in such cases integration should be performed very carefully. Indeed, making time steps relatively smaller (50 times smaller) we checked that it is sufficient to reduce the positive value around µ = 0 leaving, in principle, unchanged that λ 1 which appears at µ ≈ 0.27. To reduce this artificial effect of small positive values of λ 1 and assuming an efficient algorithm (not so much time consuming) we decided to put a minimal threshold of λ 1 = 0.1 above which the system was identified in a chaotic state (see Figs. 3b, 4a and 4b).
For a given value of p (p = 0.1) we show (in Fig. 5a for µ = 0.27 and Fig. 5b for µ = 0.25) typical phase portraits (represented by lines) and simultaneously Poincare maps (represented by points). Note that in Fig. 5a we plotted trajectories for t ∈ [1200, 1500], in arbitrary units of 1/ω, while for Fig. 5b we found transient behaviour. To show differences we plotted corresponding trajectories for t ∈ [0, 1200] and t ∈ [1200, 1500] in various colours. The characteristic time t = 1200 has been chosen in such a way to reach by our system a steady state. It is easy to see that Fig. 5a represents a non-periodic vibration state (µ = 0.27) just after the system transition to a chaotic motion (see function λ 1 versus µ in Fig. 4a). The corresponding Lyapunow exponent λ 1 = 0.175 is positive indicating on exponential divergence. In spite of a small number of Poincare stroboscopic points we can see a general view of a strange chaotic attractor. On the other hand Fig. 5b shows system behaviour state for a smaller excitation amplitude (µ = 0.25). In this case the initial trajectories follow the strange attractor complex geometry but after long enough time (t = 1200) collapse to a steady state motion. Here the corresponding Lyapunov exponent λ 1 = −0.023 is negative which means a convergent motion.
For better clarity we also plotted the corresponding time histories (Figs. 6a and b) using the same colour as in Figs. 5a and b. One can easy note that for µ = 0.27 (Figs. 5a and 6a) we have chaotic motion of the system while for µ = 0.25 (Figs. 5b and 6b) after some time of chaotic transient motion with a large amplitude of oscillations and jumping between potential wells we got steady state regular vibrations of much smaller amplitude located in a single potential well.
We have examined criteria for transition to chaotic vibrations in the Duffing system with a damping term dpt(v) = v|v| p−1 described by a term fractional exponent p. We were interested especially in p < 1.0 (p = 0.1) representing an important case which can be associated with a dry friction phenomenon [22]. The critical value of excitation amplitude µ above which the system vibrate chaotically has been estimated, in the first step, by means of the Melnikov method and later confirmed by calculating corresponding Lyapunov exponent.
The Melnikov method, is sensitive to a global homoclinic bifurcation and gives a necessary condition for excitation amplitude µ = µ c1 system in its transition to chaos [6,7]. On the other hand the largest Lyapunov exponent [23], measuring the local exponential divergences of particular phase portrait trajectories gives a sufficient condition µ = µ c2 for this transition which has obviously a higher value of the excitation amplitude µ = µ c2 > µ c1 .
Above the Melnikov transition predictions (µ > µ c1 ) we have got transient chaotic vibrations [7,8,9] as we expected drifting to a regular steady state away the fractal attraction regions separation boundary. This is typical behaviour of the system which undergoes global homoclinic bifurcation. In the region of resonance the amplitude of vibration depends on the damping exponent in a nontrivial way [18]. In particalar amplitude response rised for small fractional exponents (p < 1). Namely, such increase can be identified as an analogues of increasing system response on stronger external forcing. Thus one may expect that a chaotic threshold would be reached easier in that case. In fact, comparing Figs. 4a and b, we observe larger Lyapunov exponents for a fractional p (p = 0.1 in Fig. 4a), however from our Melnikov theory analysis we draw a different conclusion. The cases for p < 1 (Figs. 3a-b) seem to be more stable against transition to chaos then the p = 1 case. This interesting contradiction should be clearyfied by further simultaneous studies of local bifurcations and resonances in the above system.
For small p and in the limit of a small velocity we have found a region of system parameters leading to an interesting dry friction effect. Note, in our approach this limit has been obtained naturally on the same footing as other realizations. This enabled us to compare results of various damping terms (Figs. 1b and 2b). However to get reliable results in that region one must improve a simple numerical algorithm. We have not studied this effect systematically in this paper leaving it for a future publication.