Dynamical response of fractional-order nonlinear system with combined parametric and forcing excitation

Hai Jun Xing1 , Ke Shi Xiao2 , Peng Sai Song3 , Lin Ru Li4

1, 2, 3School of Mechanical Engineering, Shijiazhuang Tiedao University, Shijiazhuang, 050043, China

4School of Information Technology, Hebei University of Economics and Business, Shijiazhuang, 050043, China

1Corresponding author

Journal of Vibroengineering, Vol. 20, Issue 1, 2018, p. 793-808. https://doi.org/10.21595/jve.2017.18863
Received 17 July 2017; received in revised form 12 December 2017; accepted 20 December 2017; published 15 February 2018

Copyright © 2018 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Table of Contents Download PDF Acknowledgements References
Cite this article
Views 138
Reads 52
Downloads 1142
WoS Core Citations 0
Abstract.

A dynamical analysis of a Mathieu-van der Pol-Duffing nonlinear system with fractional-order derivative under combined parametric and forcing excitation is studied in this paper. The approximate analytical solution is researched for 1/2 sub harmonic resonance coupled with primary parametric resonance based on the improved averaging approach. The steady-state periodic solution including its stability condition is established. The equivalent linear stiffness coefficient (ELDC) and equivalent linear damping coefficient (ELSC) for this nonlinear fractional-order oscillator are defined. Then, the numerical simulations are presented in three typical cases by iterative algorithms. The time history, phase portrait, FFT spectrum and Poincare maps are shown to explain the system response. Some different responses, such as quasi-periodic, multi-periodic and single periodic behaviors are observed and investigated. The results of comparisons between the numerical solutions and the approximate analytical solutions in three typical cases show the correctness of the analytical solutions. The influences of the fractional-order parameters on the system dynamical response are researched based on the ELDC and ELSC. Through analysis, it could be found that the increase of the fractional-order coefficient would result in the rightward and downward movements of the amplitude-frequency curves. The increase of the fractional-order coefficient will also move the bifurcation point rightwards and will make the existing range of steady-state solution larger. It could also be found that the ELSC will become larger and ELDC smaller when the fractional order is closer to zero, so that the decrease of the fractional order would make the response amplitude larger. At last, the detailed conclusions are summarized, which is beneficial to design and control this kind of fractional-order nonlinear system.

Keywords: fractional-order derivative, Mathieu-van der Pol-Duffing oscillator, amplitude-frequency curves, bifurcation, averaging approach.

1. Introduction

With the development of modern computation technology in recent years, the fractional calculus has caused wide concerns in many fields, especially in mechanical, civil, biological and other engineering fields. The fundamental characteristics of fractional calculus are studied, and great progress has been made in the basic theory[1-4]. For example, Petras [1-2] studied the modeling, analysis and simulation of the fractional-order nonlinear systems, and they also focused on the modeling and control of the fractional-order systems. Caputo et al. [3] proposed a new definition of the fractional derivative. Yang et al. [4] investigated the local fractional integral transforms and the applications problem, etc.

At present, the applications of the fractional calculus in engineering fields have been attracted much attention. Many engineering problems can be described by the fractional calculus equations to simulate some dynamical systems accurately. So, the fractional calculus on the dynamical system was important and interesting, which had been investigated by many researchers in some relevant fields [5, 6]. For example, Yang et al. [7] analyzed the response of a fractional-order Duffing system by the harmonic balance method, and concluded that the fractional-order damping could change the number of the steady stable states. Wang et al. [8] investigated the dynamical characteristics of a linear oscillator with fractional-order derivative under forcing excitation, and found that the system solution had two parts in the fractional-order linear oscillator based on the theory of the stability switches. Shen et al. [9, 10] investigated several fractional-order oscillators by the averaging method, and proposed the fractional-order derivative that had both the function of damping and the function of stiffness. Many mathematical theories of the fractional calculus were researched by Li et al. [11, 12]. They also compared the differences and connections between three typical definitions of fractional-order derivatives, and presented some characteristics of the Caputo fractional calculus. Chen et al. [13] studied the dynamical response of the fractional-order stochastic oscillator based on the stochastic averaging approach, and they mainly analyzed the influences of the fractional-order derivative on the dynamics and the control of the stochastic system. Leung et al. [14] raised the residue harmonic balance approach to solve the fractional-order van der Pol (VDP) system and the higher-order approximate solution could be deduced according to this method. You et al. [15] studied the optimal control of a fractional-order vehicle suspension system, and the viscoelastic characteristics of the vehicle suspension system were investigated. Shen et al. [16] studied the fractional-order PID controller on a non-smooth oscillator and found the fractional-order PID controller was better than the traditional PID controller. Wang et al. [17] investigated the chaos behavior of the fractional-order Lorenz systems by the active feedback approach. Then the fractional order of 0.98 was implemented and the result demonstrated the validity of this method. Yang et al. [18, 19] studied several fractional-order operators in the sense of Caputo definition, and discussed their characteristics based on the Laplace transform and Fourier transform. Those results were very efficient and helpful to investigate some complex phenomenon. Li et al. [20] investigated the stochastic bifurcations of Duffing-VDP system under colored noise by the stochastic averaging method, and found that the change of fractional order, noise intensity and correlation time would lead to bifurcation of this system. Niu et al. [21] studied the Duffing oscillator with delayed fractional-order PID control based on the KBM asymptotic method, and analyzed the influences of time delay. Giresse et al. [22] studied the synchronization of chaotic Mathieu-VDP and chaotic Duffing-VDP systems with the fractional-order derivative, and found fractional-order derivative would lead to faster synchronization than the integer-order derivative.

There are many nonlinear dynamics phenomena in engineering fields, such as non-smoothness, large deformation, backlash, and parametric excitation, etc., which will deteriorate the performance of engineering system [23, 24]. For example, Li et al. [25-26] found the bursting phenomenon in a mechanical system and analyzed the slow-fast effect. Belhaq et al. [27] analyzed the resonances of VDP-Mathieu-Duffing oscillator in the frequency-locking area of 2:1 and 1:1, and found that the fast harmonic excitation could change the spring behavior. Based on the authors’ knowledge, the published papers were mainly focused on the dynamical responses of the fractional-order system with forcing excitation and the integer-order nonlinear system. A few researchers studied the fractional-order complex nonlinear system. For example, Wen et al. [28, 29] investigated the dynamical characteristics of Mathieu-Duffing oscillator with the fractional-order delayed feedback, and found that the fractional-order delayed feedback had the characteristics of both delayed velocity feedback and displacement feedback. Xu et al. [30] analyzed the dynamical responses of fractional-order Duffing-Rayleigh oscillator with the Gauss white noise excitation. But the dynamical analysis of the complex fractional-order nonlinear system with combined parametric and forcing excitation hardly has been investigated.

In this paper, the dynamical responses of Mathieu-van der Pol-Duffing nonlinear system with the fractional-order derivative for 1/2 subharmonics resonance coupled with the primary parametric resonance are studied. The paper is arranged as follows: In Section 2, the approximate analytical solution of the fractional-order Mathieu-van der Pol-Duffing nonlinear system is obtained. Then the ELDC and ELSC are also defined in this section. In Section 3, the amplitude-frequency equation and the stability conditions of the steady-state nonzero periodic solution are presented. Then in Section 4 the numerical solutions in three typical cases are obtained. The comparisons between the analytical results and the numerical solutions in three cases are presented to prove the correctness of the approximate analytical solution. The effects of the fractional-order parameters are studied in Section 5. At last, the detailed conclusions are summarized and analyzed.

2. First-order approximate analytical solution

The considered fractional-order nonlinear system with 1/2 subharmonics resonance coupled with primary parametric resonance is:

(1)
m x ¨ ( t ) + [ K 1 + 2 B c o s ( 2 ω t ) ] x ( t ) + α 1 [ x ( t ) 2 - 1 ] x ˙ ( t )
            + α 2 x 3 t + K 2 D p x t = F c o s Ω t ,

where m, K1, α1, α2, F, Ω and 2Bcos2ωt are the system mass, linear stiffness coefficient, nonlinear damping coefficient, nonlinear stiffness coefficient, forcing excitation coefficient, forcing excitation frequency and the periodic time-varying stiffing coefficient respectively. Dp[x(t)] is the p-order derivative of x(t) to t with the fractional coefficient K2 (K2>0) and fractional order p (0p1). Here the Caputo definition is used below:

(2)
D p x t = 1 Γ 1 - p 0 t x ' u t - u p d u ,

where Γx is the Gamma function which satisfies Γx+1=xΓx.

The following transformations are used:

(3)
ω 0 = K 1 m ,         B m = ε β ,         α 1 m = ε α 1 ' ,         α 2 m = ε α 2 ' ,         K 2 m = ε k 2 ,         F m = ε f .

Eq. (1) becomes:

(4)
x ¨ ( t ) + ω 0 2 x ( t ) + 2 ε β c o s ( 2 ω t ) x ( t ) + ε α 1 ' [ x ( t ) 2 - 1 ] x ˙ ( t )
          + ε α 2 ' x 3 t + ε k 2 D p x t = ε f c o s Ω t .

Here, we focus on the 1/2 sub harmonic resonance coupled with the primary parametric resonance. It means the parametric frequency is close to the natural frequency and the forcing excitation frequency is close to the double of the parametric frequency, i.e. Ω=2ω, ωω0. The formula ω2=ω02+εσ is introduced to denote the approximate degree, where σ is the detuning factor, so that Eq. (4) becomes:

(5)
x ¨ t + ω 2 x t = ε σ x t + f c o s 2 ω t - 2 β cos 2 ω t x t + α 1 ' [ 1 - x ( t ) 2 ] x ˙ ( t ) - α 2 ' x 3 ( t ) - k 2 D p [ x ( t ) ] .

Supposing the solution of Eq. (5) is:

(6)
x t = a c o s φ ,         x ˙ t = - a ω s i n φ ,

where φ=ωt+θ. Based on the improved averaging approach [9, 10], we could integrate the right-hand part of Eq. (5) within the time interval [0, T] and then average it. The approximate solution of amplitude a and phase θ are:

(7)
a ˙ = - 1 T ω 0 T P 1 a , θ + P 2 a , θ , τ s i n φ d φ ,
a θ ˙ = - 1 T ω 0 T P 1 a , θ + P 2 a , θ , τ c o s φ d φ ,

where:

P 1 a , θ = ε { f c o s [ 2 φ - θ ] + σ a c o s φ - α 1 ' a ω s i n φ + α 1 ' a 3 ω s i n 3 φ 4 + α 1 ' a 3 ω s i n φ 4
            - 2 β c o s [ 2 φ - θ ] a c o s φ - α 2 ' a 3 c o s φ } ,
P 2 a , θ = - ε k 2 D p a c o s φ .

When Eq. (7) is expanded, it will yield:

(8a)
a ˙ = a ˙ 1 + a ˙ 2 ,
(8b)
a θ ˙ = a θ ˙ 1 + a θ ˙ 2 .

The first part of Eq. (8a) and Eq. (8b) are the periodic function with T=2π. Integrating the first part of Eq. (8), one could get:

(9)
a ˙ 1 = - 1 2 π ω 0 2 π P 1 a , θ s i n φ d φ = ε a β s i n 2 θ 2 ω + ε α 1 ' a 2 - ε α 1 ' a 3 8 ,
a θ ˙ 1 = - 1 2 π ω 0 2 π P 1 a , θ c o s φ d φ = ε a β c o s 2 θ 2 ω - ε σ a 2 ω + 3 ε α 2 ' a 3 8 ω .

The second part of Eq. (8a) and Eq. (8b) are the aperiodic functions, and the time terminal T is selected as :

(10)
a ˙ 2 = - l i m T 1 T ω 0 T P 2 a , θ s i n φ d φ = l i m T ε k 2 T ω 0 T D p a cos φ s i n φ d φ ,
a θ ˙ 2 = - l i m T 1 T ω 0 T P 2 a , θ c o s φ d φ = l i m T ε k 2 T ω 0 T D p a cos φ c o s φ d φ .

Considering Eq. (2), one can obtain:

(11)
a ˙ 2 = - ε a k 2 Γ ( 1 - p ) l i m T 1 T 0 T 0 t s i n ( ω u + θ ) ( t - u ) p d u s i n ( ω t + θ ) d t .

Supposing s=t-u, ds=-du, and when they are substituted into Eq. (11), it becomes:

(12)
a ˙ 2 = - ε a k 2 Γ ( 1 - p ) l i m T 1 T 0 T 0 t s i n ( ω t + θ - ω s ) s p d s s i n ( ω t + θ ) d t
            = - ε a k 2 Γ ( 1 - p ) l i m T 1 T 0 T 0 t c o s ω s s p d s s i n ( ω t + θ ) s i n ( ω t + θ ) d t
            + ε a k 2 Γ ( 1 - p ) l i m T 1 T 0 T 0 t s i n ω s s p d s c o s ( ω t + θ ) s i n ( ω t + θ ) d t .

Two important formulae in Refs [9-10] are introduced:

(13)
l i m T 0 T s i n ( ω t ) t p d t = ω p - 1 Γ 1 - p c o s p π 2 ,
l i m T 0 T c o s ( ω t ) t p d t = ω p - 1 Γ 1 - p s i n p π 2 .

In Eq. (12), defining the first part as a˙21 and the second part as a˙22, integrating it by parts respectively and considering Eq. (13), one can get:

(14)
a ˙ 21 = - ε a k 2 4 ω Γ ( 1 - p ) l i m T 1 T 2 ω t - s i n ( 2 ω t + 2 θ ) T 0 t c o s ω s s p d s 0 T
            + ε a k 2 4 ω Γ ( 1 - p ) l i m T 1 T 0 T [ 2 ω t - s i n ( 2 ω t + 2 θ ) ] c o s ω t t p d t = - ε a k 2 ω p - 1 2 s i n p π 2 .

Similarly, one can find the integral of the second part a˙22 in Eq. (12) approach 0 when T, so that:

(15a)
a ˙ 2 = - ε a k 2 ω p - 1 2 s i n p π 2 .

When the same method is used, it will yield:

(15b)
a θ ˙ 2 = ε a k 2 ω p - 1 2 c o s p π 2 .

Combining Eq. (9) and Eq. (15), one can obtain:

(16)
a ˙ = ε a β s i n 2 θ 2 ω + ε α 1 ' a 2 - ε α 1 ' a 3 8 - ε a k 1 ω p - 1 2 s i n p π 2 ,
a θ ˙ = ε a β c o s 2 θ 2 ω - ε σ a 2 ω + 3 ε α 2 ' a 3 8 ω + ε a k 1 ω p - 1 2 c o s p π 2 - ω τ .

When the original parameters are substituted into Eq. (16), it will yield:

(17)
a ˙ = B a s i n 2 θ 2 m ω - a 2 m C e p - α 1 a 3 8 m ,
a θ ˙ = B a c o s 2 θ 2 m ω + 3 α 2 a 3 8 m ω - ω a 2 + a 2 m ω K e p ,

where:

(18)
C e p = - α 1 + K 2 ω p - 1 s i n p π 2 ,
K e p = K 1 + K 2 ω p c o s p π 2 ,

are defined as the equivalent linear damping coefficient (ELDC) and the equivalent linear stiffness coefficient (ELSC) respectively.

Based on Eq. (17) and Eq. (18), it could be found that the fractional coefficient K2 and the fractional order p are all significant in system dynamics. The fractional coefficient K2 is linear with the ELDC and ELSC in this system. The fractional order p affects this system dynamical response with the ELDC and ELSC in the form of trigonometric function. When p0 and K20, the fractional-order derivative will possess both the function of stiffness and the function of damping simultaneously.

3. Steady-state periodic solutions and stability conditions

Letting a˙=0 and θ˙=0, Eq. (17) becomes:

(19)
B a - s i n 2 θ - 2 m ω - a - 2 m C e ( p ) - α 1 a - 3 8 m = 0 ,
B a - c o s 2 θ - 2 m ω + 3 α 2 a - 3 8 m ω - ω a - 2 + a - 2 m ω K e ( p ) = 0 .

When θ- is eliminated from Eq. (19), the amplitude-frequency equation is:

(20a)
a - 2 ω α 1 a - 2 4 + ω C e ( p ) 2 + 3 α 2 a - 2 4 - m ω 2 + K e ( p ) 2 = B 2 a - 2 .

That is:

(20b)
( α 1 2 ω 2 + 9 α 2 2 ) a - 4 - 8 { ω 2 α 1 C e   ( p ) + 3 α 2 [ K e ( p ) - m ω 2 ] } a - 2               + 16 { [ K e ( p ) - m ω 2 ] 2 + [ C e ( p ) ω ] 2 - B 2 } = 0 , a - = 0 .

Letting u=acosθ, v=asinθ, and when they are substituted into Eq. (17), it becomes:

(21)
u ˙ = - C e p 2 m + α 1 u 2 + v 2 8 m u + B 2 m w v - 3 α 2 u 2 + v 2 8 m ω - ω 2 + K e p 2 m ω v = d e f U u , v , v ˙ = - C e p 2 m + α 1 u 2 + v 2 8 m v + B 2 m w u + 3 α 2 u 2 + v 2 8 m ω - ω 2 + K e p 2 m ω u = d e f V u , v .

The Jacobian matrix for Eq. (21) is:

J = U u U v V u V v .

The Jacobian matrix of the nonzero periodic solution is:

J = - C e ( p ) 2 m - α 1 ( 3 u 2 + v 2 ) 8 m - λ - 3 α 2 u v 4 m ω B 2 m ω - α 1 u v 4 m - K e ( p ) 2 m ω - ω 2 + 3 α 2 ( u 2 + 3 v 2 ) 8 m ω B 2 m ω - α 1 u v 4 m + K e ( p ) 2 m ω - ω 2 + 3 α 2 ( 3 u 2 + v 2 ) 8 m ω - C e ( p ) 2 m - α 1 ( u 2 + 3 v 2 ) 8 m - λ + 3 α 2 u v 4 m ω .

Letting:

C 1 = C e ( p ) 2 m ,         C 2 = B 2 m ω ,         C 3 = K e ( p ) 2 m ω - ω 2 ,         C 4 = α 1 8 m ,         C 5 = 3 α 2 8 m ω ,

and using Eq. (19), the characteristic equation can be get as follows:

(22)
λ 2 + ( 2 C 1 + 4 C 4 a 2 ) λ - 4 ( C 1 C 4 + C 3 C 5 ) a 2 - 4 ( C 1 2 + C 3 2 - C 2 2 ) = 0 .

From the analysis of Eq. (22), it can be found that the stability conditions for the nonzero periodic solution are:

(23)
C 1 + 2 C 4 a 2 > 0 ,
4 ( C 1 2 + C 3 2 - C 2 2 ) + 4 ( C 1 C 4 + C 3 C 5 ) a 2 < 0 .

According to Eq. (20b), one could get the following results.

When [Ke(p)-mω2]2+[Ce(p)ω]2<B2, Eq. (20b) has the zero solution and one nonzero periodic solution as follows:

(24)
a 1 = 0 ,
a 2 = 4 ω 2 α 1 C e   p + 3 α 2 K e p - m ω 2 + [ α 1 2 ω 2 + 9 α 2 2 ] B 2 - { α 1 ω [ K e ( p ) - m ω 2 ] - 3 α 2 ω C e ( p ) } 2 1 / 2 α 1 2 ω 2 + 9 α 2 2 1 / 2 .

When:

[ K e ( p ) - m ω 2 ] 2 + [ C e ( p ) ω ] 2 > B 2 ,

and:

8 ω 2 α 1 C e ( p ) + 24 α 2 [ K e ( p ) - m ω 2 ] > 0 .

Eq. (20b) has the zero solution and two nonzero periodic solutions as follows:

(25)
a 1 = 0 ,
(25a)
a 2 = 4 ω 2 α 1 C e   p + 3 α 2 K e p - m ω 2 + [ α 1 2 ω 2 + 9 α 2 2 ] B 2 - { α 1 ω [ K e ( p ) - m ω 2 ] - 3 α 2 ω C e ( p ) } 2 1 / 2 α 1 2 ω 2 + 9 α 2 2 1 / 2 ,
(25b)
a 3 = 4 ω 2 α 1 C e   p + 3 α 2 K e p - m ω 2 - [ α 1 2 ω 2 + 9 α 2 2 ] B 2 - { α 1 ω [ K e ( p ) - m ω 2 ] - 3 α 2 ω C e ( p ) } 2 1 / 2 α 1 2 ω 2 + 9 α 2 2 1 / 2 .

According to the stability conditions of Eq. (23), it can be found that the nonzero solution given by Eq. (24) is stable and there is only one nonzero solution for Eq. (1) at this moment. When Eq. (1) has two nonzero solutions, one nonzero solution given in Eq. (25a) is stable and the other one given in Eq. (25b) is unstable.

4. Numerical simulation

Numerical simulation solutions of Eq. (1) are presented in the following part. The numerical iteration methods introduced in [1-3] are adopted and the numerical formula is:

(26)
D p x t l h - p j = 0 l C j p x t l - j ,
(27)
C 0 p = 1 ,           C j p = 1 - 1 + p j C j - 1 p ,

where Cjp is the fractional binomial coefficient which satisfies the iterative relationship of Eq. (27), tl=lh is the time sample points and h is the sample step.

Based on Eq. (26) and Eq. (27), the numerical iterative algorithm of Eq. (1) could be obtained:

(28)
Z 1 t l = Z 2 t l - 1 h - j = 1 l C j 1 Z 1 t l - j ,
Z 2 ( t l ) = { F c o s ( 2 ω t l ) - α 1 ( Z 1 ( t l - 1 ) 2 - 1 ) Z 2 ( t l - 1 ) - α 2 Z 1 ( t l - 1 ) 3
            - [ k + 2 b c o s ( 2 ω t l ) ] Z 1 ( t l - 1 ) - K 1 Z 3 ( t l - 1 ) } h - j = 1 l C j 1 Z 2 t l - j ,
Z 3 t l = Z 2 t l - 1 h 1 - p - j = 1 l C j 1 - p Z 3 t l - j ,

where Z1 is the displacement, Z2 is the velocity, and Z3 is the fractional-order derivative. Here the sample step h is selected as 0.001. The temporary response is omitted, and the maximum steady amplitude of the posterior steady-state response is adopted. Numerical results and subsequent analysis of Eq. (1) are presented in three typical topology cases.

The system parameters of the first illustrative example are selected as: m=1, K1=1, α1=0.1, α2=0.1, B=0.4, K2=0.1, p=0.5 and F=0.5. The numerical results are shown in Fig. 1. As the frequency is swept forward, the amplitude-frequency curve based on the numerical iterative algorithm is shown in Fig. 1 by small triangles and circles. From the observation of Fig. 1, it can be found that in the amplitude-frequency curve there are three sudden change points and four different response regions. The following parts will simulate the dynamical motions in different regions. In order to analyze the system motion, the time history figure only shows a partial steady state process in Figs. 2-7.

Fig. 1. Amplitude-frequency curves: (m=1, K1=1, α1=0.1, α2=0.1, B=0.4, K2=0.1, p=0.5, F=0.5)

 Amplitude-frequency curves:  (m=1, K1=1, α1=0.1, α2=0.1, B=0.4, K2=0.1, p=0.5, F=0.5)

(1) When ω= 0.55, the time history, phase portrait, FFT spectrum and Poincare maps are as shown in Fig. 2 respectively. It can be found that there is only one main frequency in the FFT spectrum of Fig. 2(c), and there is only one point in Poincare maps of Fig. 2(d), so that the system response is simple single-periodic motion in this case.

(2) When ω= 0.75, the time history, phase portrait, FFT spectrum and Poincare maps are shown in Fig. 3 respectively. It can be found that there are three irreducible frequency components in the FFT spectrum of Fig. 3(c), and there is a closed loop in Poincare maps of Fig. 3(d), so that the system response is quasi-periodic motion in this case.

(3) When ω= 0.8, the time history, phase portrait, FFT spectrum and Poincare maps are shown in Fig. 4 respectively. There are several frequency components in the FFT spectrum of Fig. 4(c) and several points in Poincare maps of Fig. 4(d). It can be found that the system response is multi-periodic motion in this case which is in the midst of quasi-periodic motion and single-periodic motion.

(4) When ω= 1.0, the time history, phase portrait, FFT spectrum and Poincare maps are shown in Fig. 5 respectively. The motion characteristics are similar to those shown in Fig. 2, so that the system response also is single-periodic motion in this case.

(5) When ω=1.5, the time history, phase portrait, FFT spectrum and Poincare maps are shown in Fig. 6 respectively. When ω= 1.55, the time history, phase portrait, FFT spectrum and Poincare maps are shown in Fig. 7 respectively. It can be found that the motion characteristics are similar to those shown in Fig. 4, so that the system responses are all multi-periodic motion in those cases.

Fig. 2. System response for ω= 0.55

 System response for ω= 0.55

a) Time history

 System response for ω= 0.55

b) Phase portrait

 System response for ω= 0.55

c) FFT spectrum

 System response for ω= 0.55

d) Poincare maps

Fig. 3. System response for ω= 0.75

 System response for ω= 0.75

a) Time history

 System response for ω= 0.75

b) Phase portrait

 System response for ω= 0.75

c) FFT spectrum

 System response for ω= 0.75

d) Poincare maps

Fig. 4. System response for ω= 0.8

 System response for ω= 0.8

a) Time history

 System response for ω= 0.8

b) Phase portrait

 System response for ω= 0.8

c) FFT spectrum

 System response for ω= 0.8

d) Poincare maps

Fig. 5. System response for ω= 1.0

 System response for ω= 1.0

a) Time history

 System response for ω= 1.0

b) Phase portrait

 System response for ω= 1.0

c) FFT spectrum

 System response for ω= 1.0

d) Poincare maps

We can find that there are quasi-periodic behaviors and multi-periodic behaviors (QMP) in Region II and IV denoted by small triangles, and single-periodic behaviors (SP) in Region I and III denoted by small circles in Fig. 1. With the frequency increase, the SP, QMP, SP and then QMP responses will take place in this system gradually. The SP response in Region I is the super harmonic resonance, and the SP response in Region III is the primary resonance.

The numerical results of the other two cases are shown in Fig. 8 and Fig. 9. The amplitude-frequency curves of the approximate analytical solution of Eq. (20) in three typical cases are also shown in Fig. 1, Fig. 8 and Fig. 9 respectively, where the solid line shows the stable solution and the dash line shows the unstable one. From the observation of Fig. 1, Fig. 8 and Fig. 9, one can find that the analytical solutions are all in a good agreement with the numerical ones at the regions ω1 in three topology cases. Here 1/2 sub harmonic resonance coupled with the primary parametric resonance is focused, so that the correctness and accuracy of the analytical solutions can be verified. Furthermore, we also find that the approximate analytical solution does not agree very well with the numerical one near the bifurcation points on the both sides in Region III of Fig. 1. In order to analyze the error reason, two points ω= 0.87 and ω= 0.875 near the bifurcation points are selected to simulate time histories in Fig. 10(b) and Fig. 10(a) respectively. From the comparison of Fig. 10, it can be found that the system response approaches the steady state after 100s when ω= 0.875, but the system response approaches the steady state after 500 s when ω= 0.87. Accordingly, it could be concluded that the closer ω is to the bifurcation point, the longer the system needs the convergence time. Because of the calculation and time troubles, the fittings on both sides of the amplitude-frequency curve are not realized.

Fig. 6. System response for ω= 1.5

System response for ω= 1.5

a) Time history

System response for ω= 1.5

b) Phase portrait

System response for ω= 1.5

c) FFT spectrum

System response for ω= 1.5

d) Poincare maps

Fig. 7. System response for ω= 1.55

System response for ω= 1.55

a) Time history

System response for ω= 1.55

b) Phase portrait

System response for ω= 1.55

c) FFT spectrum

System response for ω= 1.55

d) Poincare maps

5. Influences of fractional-order derivative parameters

In order to optimize the system design and to control the system response through choosing the fractional-order derivative parameters, the influences of the fractional-order derivative parameters on the system dynamical responses are investigated. When K2 is changed, the system amplitude-frequency curves of α1=0.4 and α1=0.8 are given in Fig. 11(a, b) respectively. The system basic parameters are selected as m=1, K1=1, α2=0.1, B=0.4, F=0.5 and p=0.5.

Fig. 8. Comparison between analytical results and numerical solutions: (m=1, K1=1, α1=0.42, α2=0.01, B=0.4, K2=0.1, p=0.5, F=0.5)

 Comparison between analytical results and numerical solutions:  (m=1, K1=1, α1=0.42, α2=0.01, B=0.4, K2=0.1, p=0.5, F=0.5)

Fig. 9. Comparison between analytical results and numerical solutions: (m=1, K1=1, α1=0.8, α2=0.01, B=0.4, K2=0.1, p=0.5, F=0.5)

 Comparison between analytical results and numerical solutions:  (m=1, K1=1, α1=0.8, α2=0.01, B=0.4, K2=0.1, p=0.5, F=0.5)

Fig. 10. System response for ω= 0.875 and 0.87: Time history

 System response for ω= 0.875 and 0.87: Time history

a)ω= 0.875

 System response for ω= 0.875 and 0.87: Time history

b)ω= 0.87

From the observation of Fig. 11, it could be found that the increase of K2 would result into the rightward and downward movements of the amplitude-frequency curves. The main reason is that the ELSC will become larger with the increase of K2. Besides the above characteristics of the dynamical system responses, the topological structures of the amplitude-frequency curves vary significantly too. The increase of fractional-order coefficient K2 would make the ovoid curve disappear, the open curve gradually turn up, and push the bifurcation point rightward. The proportion of steady-state solution in the amplitude-frequency curve is getting bigger and bigger with the increase of K2.

Fig. 11. Influences of K2

 Influences of K2

a)α1=0.4

 Influences of K2

b)α1=0.8

When p changes from 1 to 0, the amplitude-frequency curves of the fractional coefficient K2=0.1 and K2=0.5 are given in Fig. 12(a, b) respectively. The basic system parameters are selected as m=1, K1=1, α2=0.1, B=0.4, F=0.5 and α1=0.4. Through the analysis of Fig. 12, one can find that the decrease of p will result into the rightward and upward movements of the amplitude-frequency curves. The open curve disappears, the ovoid curve gradually turns up, and the bifurcation point will move rightwards with the decrease of p. Through the analysis of Eq. (17), one can find that the reason is that the ELSC will become larger, and the ELDC will become smaller when the fractional order p0. The ELSC will become smaller, and the ELDC will become larger when the fractional order p1. Accordingly, the decrease of p will make the maximum response amplitude larger.

Fig. 12. Influences of the fractional order p

 Influences of the fractional order p

a)K2= 0.1

 Influences of the fractional order p

b)K2= 0.5

6. Conclusions

The dynamical characteristics of a Mathieu-van der Pol-Duffing fractional-order nonlinear system under 1/2 subharmonic resonance coupled with primary parametric resonance are investigated. The approximate analytical system solution is presented by the improved averaging approach. The numerical solutions are also studied by an iterative approach in three typical cases. The time history, phase portrait, FFT spectrum and Poincare maps are presented to identify the system response. There are Quasi-periodic, multi-periodic or single periodic behaviors. As the frequency is swept forward, the QMP response, SP response and QMP response will takes place in this system gradually, which is observed and investigated in three cases. Then the numerical results are compared with the analytical results. From three comparisons between the numerical and analytical results, one can find that the stable periodic solutions are in a good agreement with the numerical solutions, which verify the correctness and high accuracy of the approximate analytical solution. Finally, the influences of the fractional-order derivative parameters on the system dynamical response are analyzed. Through the analysis of the fractional-order derivative parameters, it can be found that the ELSC and ELDC both will become larger with the increase of fractional-order coefficient K2. The ELSC will become larger, and the ELDC will become smaller when the fractional order p0. The ELSC will become smaller, and the ELDC will become larger if the fractional order p1. A change of the fractional-order derivative parameters would lead to the motion of the amplitude-frequency curve and to a change of the topological structure of the amplitude-frequency curve. Those conclusions are very useful for the system design and control of the system dynamical response by choosing suitable fractional-order derivative parameters.

Acknowledgements

The authors are grateful to the support by the National Natural Science Foundation of China (No. 11602152), and the department project in Hebei Province (14212202D).

References

  1. Petras I. Fractional-order Nonlinear Systems: Modeling, Analysis and Simulation. Higher Education Press, Beijing, 2011. [Publisher]
  2. Caponetto R., Dongola G., Fortuna L., Petras I. Fractional Order Systems: Modeling and Control Applications. World Scientific, New Jersey, 2010. [Publisher]
  3. Caputo M., Fabrizio M. A new definition of fractional derivative without singular kernel. Progress in Fractional Differentiation and Applications, Vol. 1, Issue 2, 2015, p. 73-85. [Search CrossRef]
  4. Yang X. J., Baleanu D., Srivastava H. M. Local Fractional Integral Transforms and Their Applications. Academic Press, London, 2015. [Search CrossRef]
  5. Agnieszka B., Malinowska, Delfim Torres F. M. Fractional calculus of variations for a combined Caputo derivative. Fractional Calculus and Applied Analysis, Vol. 14, Issue 4, 2011, p. 523-537. [Search CrossRef]
  6. Machado J. T., Kiryakova V., Mainardi F. Recent history of fractional calculus. Communications in Nonlinear Science and Numerical Simulation, Vol. 16, Issue 3, 2011, p. 1140-1153. [Publisher]
  7. Yang J. H., Zhu H. The response property of one kind of fractional-order linear system excited by different periodical signals. Acta Physica Sinica, Vol. 62, Issue 2, 2013, p. 24501, (in Chinese). [Search CrossRef]
  8. Wang Z. H., Hu H. Y. Stability of a linear oscillator with damping force of fractional-order derivative. Science in China G: Physics, Mechanics and Astronomy, Vol. 39, Issue 10, 2009, p. 1495-1502, (in Chinese). [Search CrossRef]
  9. Shen Y. J., Wei P., Yang S. P. Primary resonance of fractional-order van der Pol oscillator. Nonlinear Dynamics, Vol. 77, Issue 4, 2014, p. 1629-1642. [Publisher]
  10. Shen Y. J., Yang S. P., Sui C. Y. Analysis on limit cycle of fractional-order van der Pol oscillator. Chaos, Solitons and Fractals, Vol. 67, 2014, p. 94-102. [Publisher]
  11. Li C. P., Deng W. H. Remarks on fractional derivatives. Applied Mathematics and Computation, Vol. 187, Issue 2, 2007, p. 777-784. [Publisher]
  12. Deng W. H., Li C. P. The evolution of chaotic dynamics for fractional unified system. Physics Letters A, Vol. 372, Issue 4, 2008, p. 401-407. [Publisher]
  13. Chen L. C., Hu F., Zhu W. Q. Stochastic dynamics and fractional optimal control of quasi integrable Hamiltonian systems with fractional derivative damping. Fractional Calculus and Applied Analysis, Vol. 16, Issue 1, 2013, p. 189-225. [Publisher]
  14. Leung A. Y. T., Guo Z. J., Yang H. X. The residue harmonic balance for fractional order van der Pol like oscillators. Journal of Sound and Vibration, Vol. 331, Issue 5, 2012, p. 1115-1126. [Publisher]
  15. You H., Shen Y. J., Yang S. P. Optimal design for fractional-order active isolation system. Advances in Mechanical Engineering, Vol. 7, Issue 12, 2015, p. 1-11. [Publisher]
  16. Shen Y. J., Niu J. C., Yang S. P., Li S. J. Primary resonance of dry-friction oscillator with fractional-order PID Controller of velocity feedback. ASME Journal of Computational and Nonlinear Dynamics, Vol. 11, Issue 5, 2016, p. 51027. [Search CrossRef]
  17. Wang X. Y., Song J. M. Synchronization of the fractional order hyperchaos Lorenz systems with activation feedback control. Communications in Nonlinear Science and Numerical Simulation, Vol. 14, Issue 8, 2009, p. 3351-3357. [Publisher]
  18. Yang X. J. Fractional derivatives of constant and variable orders applied to anomalous relaxation models in heat-transfer problems. Thermal Science, Vol. 21, Issue 3, 2017, p. 1161-1171. [Publisher]
  19. Yang X. J., Machado J. A. T. A new fractional operator of variable order: application in the description of anomalous diffusion. Physica A: Statistical Mechanics and its Applications, Vol. 481, 2017, p. 276-283. [Publisher]
  20. Li W., Zhang M. T., Zhao J. F. Stochastic bifurcations of generalized Duffing-van der Pol system with fractional derivative under colored noise. Chinese Physics B, Vol. 26, 2017, p. 90501. [Search CrossRef]
  21. Niu J. C., Shen Y. J., Yang S. P. Analysis of duffing oscillator with time-delayed fractional-order PID controller. International Journal of Non-Linear Mechanics, Vol. 92, 2017, p. 66-75. [Publisher]
  22. Giresse T. A., Crépin K. T. Chaos generalized synchronization of coupled Mathieu-Van der Pol and coupled Duffing-Van der Pol systems using fractional order-derivative. Chaos, Solitons and Fractals, Vol. 98, 2017, p. 88-100. [Publisher]
  23. Nayfeh A. H. Nonlinear Oscillations. Wiley, New York, 1973. [Search CrossRef]
  24. Chen Y. S. Nonlinear Vibrations. Higher Education Press, Beijing, 2002. [Search CrossRef]
  25. Li X. H., Hou J. Y., Shen Y. J. Slow-fast effect and generation mechanism of Brusselator based on coordinate transformation. Open Physics, Vol. 14, Issue 1, 2016, p. 261-268. [Publisher]
  26. Li X. H., Hou J. Y. Bursting phenomenon in a piecewise mechanical system with parameter perturbation in stiffness. International Journal of Non-Linear Mechanics, Vol. 81, 2016, p. 65-176. [Publisher]
  27. Belhaq M., Fahsi A. 2:1 and 1:1 frequency-locking in fast excited van der Pol-Mathieu-Duffing oscillator. Nonlinear Dynamics, Vol. 53, Issues 1-2, 2008, p. 139-152. [Publisher]
  28. Wen S. F., Shen Y. J., Yang S. P., Wang J. Dynamical response of Mathieu-Duffing oscillator with fractional-order delayed feedback. Chaos, Solitons and Fractals, Vol. 94, 2017, p. 54-62. [Publisher]
  29. Wen S. F., Shen Y. J., Li X. H., Yang S. P. Dynamical analysis of Mathieu equation with two kinds of van der Pol fractional-order terms. International Journal of Non-Linear Mechanics, Vol. 84, 2016, p. 130-138. [Publisher]
  30. Xu Y., Li Y. G., Liu D., Jia W. T., Huang H. Responses of Duffing oscillator with fractional damping and random phase. Nonlinear Dynamics, Vol. 74, Issue 3, 2013, p. 745-753. [Publisher]

Cited By

International Journal of Nonlinear Sciences and Numerical Simulation
Jianhua Tang, Chuntao Yin
2021