Bifurcation and Chaotic Behavior of Duffing System with Fractional-Order Derivative and Time Delay

Abstract

In this paper, the abundant nonlinear dynamical behaviors of a fractional-order time-delayed Duffing system under harmonic excitation are studied. By constructing Melnikov function, the necessary conditions of chaotic motion in horseshoe shape are detected, and the chaos threshold curve is obtained by comparing the results obtained through the Melnikov theory and numerical iterative algorithm. The results show that the trend of change is the same, which confirms the accuracy of the chaos threshold curve. It could be found that when the excitation frequency ω is larger than a certain value, the Melnikov theory is not valid for these values. Furthermore, by numerical simulation, some numerical results are obtained, including phase portraits, the largest Lyapunov exponents, and the bifurcation diagrams, Poincare maps, time histories, and frequency spectrograms at some typical points. These numerical simulation results show that the system exhibits some new complex dynamical behaviors, including entry into the state of chaotic motion from single period to period-doubling bifurcation and chaotic motion and periodic motion alternating under the necessary condition of chaotic occurrence. In addition, the effects of time delay, fractional-order coefficient, fractional order, linear viscous damping coefficient, and linear stiffness coefficient on the chaotic threshold curve are discussed, respectively. Those results reveal that there exist abundant nonlinear dynamic behaviors in this fractional-order system, and by adjusting these parameters reasonably, the system could be transformed from chaotic motion to non-chaotic motion.

1. Introduction

Many scholars have carried out in-depth research on the basic theory of the properties and characteristics of fractional calculus. Since the last century, fractional calculus theory has been gradually applied to the scientific research of engineering problems [1,2,3,4,5]. The main advantage of fractional calculus is memory and genetic properties. Therefore, fractional-order derivative has been applied in multiple fields, such as electrochemistry signal processing, anomalous diffusion process, viscoelastic damping, and so on [6,7,8]. Especially in some engineering fields, many researchers have applied the fractional-order model to study complicated dynamical problems [9,10].

Bifurcation and chaotic behaviors are complex vibration phenomena specific in nonlinear dynamical systems. Chaos is generated from deterministic systems, and it has complex dynamic characteristics; for example, its sensitivity to initial conditions and irregular and unpredictable behaviors for a long time could be employed in the application of secure communication and cryptography [11,12]. The Melnikov method is particularly suitable for the analysis of chaos threshold of multi-parameter nonlinear dynamical systems. Recently, many researchers have predicted chaotic motion by Melnikov’s theorem. Liu et al. [13] used the Melnikov method to study the chaotic motion of VEH system under random excitation based on fractional order physical properties. Battelli et al. [14] studied the chaotic behaviors in time-perturbed discontinuous systems based on the Melnikov method. Lian et al. [15] studied the chaotic motion and control in a tethered-sailcraft system orbiting an asteroid by comparing the Melnikov method with the numerical method. Tuwa et al. [16] used the Melnikov’s theorem to analyze the chaotic motion of nonlinear viscoelastic plate with fractional derivative model, and the numerical results demonstrate the validity of theoretical prediction results. Liang et al. [17] studied the Poincaré bifurcation of a planar piecewise near-Hamiltonian system. Farshidianfar et al. [18] used Melnikov’s method to analyze the global bifurcation and the transition to chaotic behavior of nonlinear gear systems. Tian et al. [19,20] used the extended Melnikov theory to determine the bifurcation criterion of inverted pendulum system under the condition of impulsive excitation and verified phase portraits, bifurcation diagrams, and Poincare maps by numerical methods.

In the field of nonlinear dynamics, Duffing system is a forced vibration system with nonlinear restoring force, which is widely representative in various fields of engineering. Duffing system is recognized as the prototype of systems with large deformation or similar properties in the field of physics and engineering and is one of the research hot spots. Based on Lei et al. [21]’s fractional deflection generalized Duffing system, the effects of color noise, green noise, and red noise on chaos initiation were investigated theoretically and verified numerically. Yang et al. [22] studied the Duffing equation under the action of single external force and parametric excitation and obtained abundant dynamic behaviors of bifurcation and chaos. Yang. et al. [23] studied the resonance of fractional order Duffing systems.

Time delay is a ubiquitous phenomenon in nature. It is more common and inevitable in control systems. It exists in two ways; one is in the internal state of the system, and the other is introduced as a control strategy to enhance system performance [24]. Many scholars have studied the fractional-order delayed dynamical systems. Leung. et al. [25] investigated a Duffing–van der Pol oscillator with fractional-order derivative and time delay. Based on the residue harmonic method, Wen et al. [26] studied the dynamical response of fractional-order time-delayed feedback for Mathieu–Duffing oscillator. Yang et al. [27] analyzed the resonance phenomena of fractional Duffing systems based on linear time-delay feedback with direct separation of slow and fast motion. Shen et al. [28,29] investigated the bifurcation and chaotic behaviors of a Duffing oscillator with delayed displacement and velocity feedbacks by Melnikov’s theorem. Mesbahia. et al. [30] used harmonic balance method to investigate the dynamical behaviors of Mathieu equation with fractional-order and damped time delay. Çelik et al. [31] used the largest Lyapunov exponents to study the existence of chaotic motion in nonlinear autonomous systems with fractional-order delay.

Through the analysis and summary of the above scholars’ research content, most of the existing research studies are focused on the numerical or qualitative analysis of the necessary condition for chaotic motion. At present, few researchers are working on the necessary condition for chaos and analyzing the complex bifurcation and chaotic dynamical behaviors of a fractional-order time-delayed Duffing system.

Based on the Melnikov theory, the chaos criterion of a typical fractional-order time-delay Duffing system under harmonic excitation is studied, and it exhibits some complex nonlinear dynamic phenomena, including bifurcation and chaotic behaviors. In this paper, Melnikov function is established by the Melnikov theorem, and the chaos threshold curve is obtained. Then, the necessary conditions are verified by numerical calculation. In addition, the bifurcation diagrams and the corresponding largest Lyapunov exponents, as well as the phase portraits, Poincare maps, time histories, and frequency spectrograms at typical system parameters by numerical solution, are investigated, which give full support to the theoretical analysis. The influences of the time-delay parameter, fractional-order coefficient, fractional order, linear stiffness coefficient, and linear viscous damping coefficient on chaotic threshold curve are analyzed. In the last section, the primary conclusions are analyzed and summarized.

2. Results

As we know, time delay is a common phenomenon when modeling some dynamic systems, such as mechanical dynamics, neural networks and biological systems, process control systems, etc. [32,33,34,35]. Fractional derivative is used to describe the dissipative behavior of the dynamic system. In this paper, the considered fractional-order time-delayed Duffing system is

$$m\ddot{x}-k_{1}x+k_3{x}^3+c_1\dot{x}+h{D}_t{}^p\left[x(t-\tau )\right]=F\cos (\omega t)$$

(1)

where m, c₁, k₁, and k₃ denote the system mass, linear viscous damping coefficient, linear stiffness coefficient, and nonlinear stiffness coefficient. Moreover, F and ω denote the excitation amplitude and excitation frequency of the system, respectively. τ is the time-delay parameter in fractional-order term, and the expression hD_t^p[x(t − τ)] is the p-order derivative of x(t − τ) with respect to t with fractional-order coefficient h, and the fractional order is restricted as 0 ≤ p ≤ 1. It should be mentioned that the fractional-order term in Equation (1) plays an important role that induces different bifurcation and chaotic motion [27]. The Duffing oscillator is considered as one of the most important nonlinear dynamics in physics, biology, and even economics. Many research subjects can be abstracted as Duffing systems. The Caputo’s definition is adopted in Equation (1) with the form as follows.

$$D_t^p\left[x\left(t\right)\right]=\frac{1}{\mathrm{\Gamma }\left(n-p\right)}{\displaystyle {\int }_0^t\frac{x^{\left(n\right)}\left(\tau \right)}{{\left(t-\tau \right)}^{p-n+1}}\mathrm{d}\tau }$$

(2)

where n − 1 < p < n and n ∈ N, and Γ(•) is the gamma function satisfying Γ(y + 1) = y Γ(y).

In practical engineering applications, there may be large deformation or harmonic excitation. The dynamic response of the foundation modes of simply supported Euler beams can be expressed in Formula (1). When the coefficient k₁ is positive, the compressive axial force exceeds the first Eulerian buckling load [36]. When k₁ = 1 and k₃ = 1, it is known as Holmes-type Duffing system, which was obtained by P. Holmes when he studied the vibration of magnetoelastic beams. Many mathematical models of mechanical problems could be summed up to Duffing equation, such as nonlinear oscillation of electronic circuit system, perturbation problem of conservative systems, and so on.

Using the following transformations

$$b_1=k_1/m,B=\epsilon b=c_1/m,b_2=k_3/m,{\mu }_1=\epsilon \mu =h/m,F_0=\epsilon f=F/m$$

where ε is a small real parameter, and it satisfies 0 < ε << 1. Equation (1) becomes

$$\ddot{x}-b_{1}x+b_2{x}^3=\epsilon \left\{f\cos (\omega t)-\mu D_t{}^p\left[x(t-\tau )\right]-b\dot{x}\right\}$$

(3)

Equation (3) could be transformed into the following form

$$\{\begin{array}{l}\dot{x}=y\\ \dot{y}=b_{1}x-b_2{x}^3+\epsilon \left\{f\cos (\omega t)-\mu D_t{}^p\left[x(t-\tau )\right]-by\right\}\end{array}$$

(4)

If ε is zero, system (1) is assumed to be an unperturbed system, then it can be written as

$$\ddot{x}-b_{1}x+b_2{x}^3=0$$

(5)

Here, the homoclinic orbits satisfy the formula

$$\frac{1}{2}{\dot{x}}^2-\frac{b_1}{2}{x}^2+\frac{b_2}{4}{x}^4=0$$

(6)

Supposing $\dot{x}$ = 0 at t = 0, the calculating result is

$$x_0=\pm \sqrt{\frac{2b_1}{b_2}}$$

(7)

Integrating Equation (6), it could yield

$${\displaystyle {\int }_{x_0}^x\frac{\mathrm{d}x}{\pm \sqrt{b_1{x}^2-\frac{b_2}{2}{x}^4}}}=t$$

(8)

Calculating Equation (8), one could obtain the homoclinic trajectory, and it is given as follows.

$$\left(x,y\right)=\left(\pm \sqrt{\frac{2b_1}{b_2}}\sec \mathrm{h}\left(\sqrt{b_1}t\right),\mp b_1\sqrt{\frac{2}{b_2}}\sec \mathrm{h}\left(\sqrt{b_1}t\right)\tanh (\sqrt{b_1}t)\right)$$

(9)

In the following part, the main purpose is to apply the Melnikov theory to system (1) to obtain the necessary condition under which the system may generate chaos. In Equation (1), the parameters F, h, and c are assumed to be small parameters. Hence, the system of Equation (2) could be rewritten as

$$\dot{\overrightarrow{x}}=f\left(\overrightarrow{x}\right)+\epsilon g\left(\overrightarrow{x},t\right)$$

(10)

where

$$\overrightarrow{x}=\left(\begin{array}{l}{x}_1\\ x_2\end{array}\right)=\left(\begin{array}{c}x\\ \dot{x}\end{array}\right),f\left(\overrightarrow{x}\right)=\left(\begin{array}{l}{x}_2\\ b_{1}x-b_2{x}^3\end{array}\right),g\left(\overrightarrow{x},t\right)=\left[\begin{array}{l}0\\ f\cos (\omega t)-\mu D_t{}^p\left[x(t-\tau )\right]-b\dot{x}\end{array}\right]$$

According to the Melnikov process [36], one could establish Melnikov function as follows.

$$\begin{array}{cc}M\left(t_0\right)\hfill & ={\displaystyle {\int }_{-\infty }^{\infty }f\left(u^0\left(s-t_0\right)\right)}\wedge g\left(u^0\left(s-t_0\right),s\right)\mathrm{d}s={\displaystyle {\int }_{-\infty }^{\infty }f\left(\overrightarrow{x}\right)\wedge g\left(\overrightarrow{x},s\right)\mathrm{d}s}\hfill \\ \hfill & =\epsilon {\displaystyle {\int }_{-\infty }^{\infty }{x}_2\left(s-t_0\right)\times \left[f\cos (\omega s)-\mu D_{s-t_0}^p\left(x\left(s-t_0-\tau \right)\right)-b{x}_2\left(s-t_0\right)\right]}\mathrm{d}s\hfill \end{array}$$

(11)

Letting t = s − t₀, Equation (11) becomes

$$\begin{array}{cc}M\left(t_0\right)\hfill & =\epsilon {\displaystyle {\int }_{-\infty }^{\infty }{x}_2\left(t\right)\times \left[f\cos (\omega t+\omega t_0)-\mu D_t^p\left(x\left(t-\tau \right)\right)-b{x}_2\left(t\right)\right]}\mathrm{d}t\hfill \\ \hfill & =M_1\left(t_0\right)+M_2\left(t_0\right)+M_3\left(t_0\right)\hfill \end{array}$$

(12)

where t₀ is the cross-section time of the Poincare map, which satisfies t₀ ∈ [0, T), and T is the period of the system motion.

Substituting the homoclinic orbits x and y into Equation (12) and evaluating the integral, one could obtain the Melnikov function. The detailed computational process would be presented separately in the following text.

2.1. Calculating M₁(t₀) and M₃(t₀)

Since M₁(t₀) and M₃(t₀) in Equation (12) do not contain fractional-order differential terms, they could be integrated directly. The calculation processes are as follows.

$$M_1\left(t_0\right)=\epsilon f{\displaystyle {\int }_{-\infty }^{\infty }{x}_2\left(t\right)\cos (\omega t+\omega t_0)\mathrm{d}t}$$

(13)

According to the odevity of function, given that x₂(t) is an odd function of t, Equation (13) becomes

$$\begin{array}{cc}{M}_1\left(t_0\right)\hfill & =\epsilon f{\displaystyle {\int }_{-\infty }^{\infty }{x}_2\left(t\right)\cos (\omega t+\omega t_0)\mathrm{d}t}\hfill \\ \hfill & =\epsilon f{\displaystyle {\int }_{-\infty }^{\infty }{x}_2\left(t\right)\left(\cos \omega t\cos \omega t_0-\sin \omega t\sin \omega t_0\right)\mathrm{d}t}\hfill \\ \hfill & =-\sin \left(\omega t_0\right)\epsilon f{\displaystyle {\int }_{-\infty }^{\infty }{x}_2\left(t\right)\sin (\omega t)\mathrm{d}t}\hfill \\ \hfill & =\pm b_1\sqrt{\frac{2}{b_2}}\sin \left(\omega t_0\right)\epsilon f{\displaystyle {\int }_{-\infty }^{\infty }\sec \mathrm{h}\left(\sqrt{b_1}t\right)\tanh (\sqrt{b_1}t)\sin (\omega t)\mathrm{d}t}\hfill \\ \hfill & =\pm \sqrt{\frac{2}{b_2}}\sin \left(\omega t_0\right)\epsilon f\pi \omega \sec \mathrm{h}\left(\frac{\pi \omega }{2\sqrt{b_1}}\right)\hfill \end{array}$$

(14)

$$\begin{array}{cc}{M}_3\left(t_0\right)& =-\epsilon b{\displaystyle {\int }_{-\infty }^{\infty }{x}_2^2\left(t\right)}\mathrm{d}t\hfill \\ \hfill & =-\frac{2\epsilon b{b}_1^2}{b_2}{\displaystyle {\int }_{-\infty }^{\infty }\sec {\mathrm{h}}^2\left(\sqrt{b_1}t\right){\tanh }^2(\sqrt{b_1}t)}\mathrm{d}t\hfill \\ \hfill & =-\frac{4\epsilon b{b}_1^2}{3\sqrt{b_1}{b}_2}\hfill \end{array}$$

(15)

2.2. Calculating M₂(t₀)

We can get

$$\begin{array}{cc}{M}_2\left(t_0\right)& =\epsilon {\displaystyle {\int }_{-\infty }^{+\infty }{x}_2\left(t\right)\mu D_t^p\left(x\left(t-\tau \right)\right)\mathrm{d}t}\hfill \\ \hfill & =\mp \epsilon b_1\mu \sqrt{\frac{2}{b_2}}{\displaystyle {\int }_{-\infty }^{+\infty }\sec \mathrm{h}\left(\sqrt{b_1}t\right)\tanh (\sqrt{b_1}t)D_t^p\left(\pm \sqrt{\frac{2b_1}{b_2}}\sec \mathrm{h}\left(\sqrt{b_1}\left(t-\tau \right)\right)\right)\mathrm{d}t}\hfill \\ \hfill & =-\frac{4\mu \epsilon b_1\sqrt{b_1}}{b_2}{\displaystyle {\int }_0^{+\infty }\sec \mathrm{h}\left(\sqrt{b_1}t\right)\tanh (\sqrt{b_1}t)D_t^p\left(\sec \mathrm{h}\left(\sqrt{b_1}\left(t-\tau \right)\right)\right)\mathrm{d}t}\hfill \\ \hfill & =-\frac{4\mu \epsilon b_1\sqrt{b_1}}{b_2}A\hfill \end{array}$$

(16)

where $A={\displaystyle {\int }_0^{+\infty }\sec \mathrm{h}\left(\sqrt{b_1}t\right)\tanh (\sqrt{b_1}t)D_t^p\left(\sec \mathrm{h}\left(\sqrt{b_1}\left(t-\tau \right)\right)\right)\mathrm{d}t}$.

Since there is a fractional-order differential term in Equation (16), it could not be integrated directly. Firstly, the fractional-order differential is calculated, and then the generalized integral is calculated. The detailed computational process of Equation (16) can be found in [37].

2.3. The Threshold Curve of the Chaotic Motion

According to the Melnikov theory [38], if there exists t₀ when M(t₀) = 0 and dM(t₀)/dt₀ ≠ 0, then lateral intersection for the stable manifold (W_s^±) and the unstable manifold (W_u^±) occurs, which means the system generates chaos in the sense of Smale horseshoes. Combining Equation (14) with Equation (16), one could obtain the Melnikov function

$$M\left(t_0\right)=\pm \sqrt{\frac{2}{b_2}}\sin \left(\omega t_0\right)\epsilon f\pi \omega \sec \mathrm{h}\left(\frac{\pi \omega }{2\sqrt{b_1}}\right)-\frac{4\mu \epsilon b_1\sqrt{b_1}}{b_2}A-\frac{4\epsilon b{b}_1^2}{3\sqrt{b_1}{b}_2}=0$$

(17)

$$\frac{\mathrm{d}M\left(t_0\right)}{\mathrm{d}{t}_0}=\pm \sqrt{\frac{2}{b_2}}\cos \left(\omega t_0\right)\epsilon f\pi {\omega }^2\sec \mathrm{h}\left(\frac{\pi \omega }{2\sqrt{b_1}}\right)\ne 0$$

(18)

Since |sin(ωt₀)| < 1, from Equations (17) and (18), one could obtain the necessary condition for the chaos

$$\left|\pm \sqrt{\frac{2}{b_2}}\epsilon f\pi \omega \sec \mathrm{h}\left(\frac{\pi \omega }{2\sqrt{b_1}}\right)\right|>\left|\frac{4\mu \epsilon b_1\sqrt{b_1}}{b_2}A+\frac{4\epsilon b{b}_1^2}{3\sqrt{b_1}{b}_2}\right|$$

(19)

Simplifying Equation (19), one could obtain

$$\sqrt{\frac{2}{b_2}}\epsilon f\pi \omega \sec \mathrm{h}\left(\frac{\pi \omega }{2\sqrt{b_1}}\right)>\left|\frac{4\mu \epsilon b_1\sqrt{b_1}}{b_2}A+\frac{4\epsilon b{b}_1^2}{3\sqrt{b_1}{b}_2}\right|$$

(20)

Substituting the original system parameters into Equation (20), one could yield the necessary condition for generating chaos in the sense of Smale horseshoes as follows.

$$\sqrt{\frac{2m}{k_3}}F\pi \omega \sec \mathrm{h}\left(\frac{\pi \omega }{2}\sqrt{\frac{m}{k_1}}\right)>\left|\frac{4h{k}_1}{k_3}\sqrt{\frac{k_1}{m}}A+\frac{4c_1{k}_1^2}{3m{k}_3}\sqrt{\frac{m}{k_1}}\right|$$

(21)

According to Equation (21), the necessary condition for transverse intersections of stable manifolds (W_s^±) and the unstable manifold (W_u^±) could be obtained. An illustrative example is investigated herein, as defined by system parameters m = 1, k₁ = 1, k₃ = 1, c₁ = 0.4, h = −0.5, p = 0.5, and τ = 0.5. After substituting those system parameters into Equation (19), the relationship curve between the critical excitation amplitude F and the excitation frequency ω is obtained and shown in Figure 1.

From the analysis of Figure 1, it could be found that the occurrence of chaotic motion is affected by the external excitation amplitude F and frequency ω. With the increase of excitation frequency ω, the chaotic threshold would be larger initially, then smaller, and finally towards zero. In region 1, the system may produce chaotic motion; in region 2, the system has periodic motion.

3. Dynamical Analysis of Fractional-Order Time-Delayed Duffing System

In this part, the numerical simulation is used to verify the above theoretical results and continue to study other complex nonlinear dynamic problems, including bifurcation and chaotic behaviors. In order to demonstrate the validity of the chaotic threshold curve obtained based on the Melnikov method in Figure 1, the numerical simulations of the system are carried out. The power series expansion method (PSE) [39] is applied to numerical simulations for Equation (1). The discretization process of fractional-order differential terms is as follows.

$$D^p\left[x\left(t_k\right)\right]\approx d^{-p}{\displaystyle \sum _{j=0}^l{C}_j^{p}x\left(t_{k-j}\right)}$$

(22)

where t_k = ld is the sampling points for time, d is the time step of calculation, and C_j^p is the binomial coefficients of fractional-order differential terms. It has the following iterative relationship

$$C_0^p=1,C_j^p=\left(1-\frac{1+p}{j}\right)C_{j-1}^p$$

(23)

Since the fractional-order derivative in Equation (1) involves time delay, one could not utilize Equation (23) to calculate directly. Therefore, in the process of numerical calculation, letting τ = i × d, (i is natural number), it could yield

$$D^p\left[x\left(t-\tau \right)\right]=D^p\left[x\left(t-id\right)\right]$$

(24)

Based on Equations (22)–(24), one could obtain the numerical iterative algorithm of Equation (1) as

$$x\left(t_k\right)=y\left(t_{k-1}\right)h-{\displaystyle \sum _{j=1}^l{c}_j^{1}x\left(t_{k-j}\right)}$$

(25)

$$y\left(t_k\right)=\left(F\cos \left(\omega t_k\right)+hz\left(t_{k-1-i}\right)+k_{1}x\left(t_{k-1}\right)-k_3{x}^3\left(t_{k-1}\right)-c_{1}y\left(t_{k-1}\right)\right)h-{\displaystyle \sum _{j=1}^l{c}_j^{1}y\left(t_{k-j}\right)}$$

(26)

$$z\left(t_k\right)=h^{q}y\left(t_{k-1}\right)-{{\displaystyle \sum _{j=1}^l{c}_j^{q}z\left(t_{k-j}\right)}}_{\left(q=1-p\right)}$$

(27)

where, x(t_k) is displacement, y(t_k) is velocity, and z(t_k) is the fractional-order derivative of displacement. In the process of numerical calculation, select time step d = 0.01, and the total computation time is generally 500 excitation periods. In total, 200 points are sampled in each period. The frontal transient response is omitted, and the posterior 100 periods (the steady response of the numerical results) are selected for analysis. The initial value is [x₀, y₀, z₀] = [1, 1, 1]. According to the above numerical iteration process, the following numerical simulation analysis is carried out for the system. In the following parts, the values for the system parameters are chosen as follows: m = 1, k₁ = 1, k₃ = 1, c₁ = 0.4, h = −0.5, p = 0.5, and τ = 0.5.

3.1. Numerical Solutions of the Chaotic Threshold Curve

The correctness of the approximate analytical solution is further verified by analyzing the numerical results. According to the numerical iterative algorithm of Equations (25)–(27), one could obtain the critical value of F at each external excitation frequency when the system generates chaos. The largest Lyapunov exponent σ is calculated by changing the amplitude F at each fixed frequency ω. When σ changes from negative to positive for the first time, the amplitude F is the chaotic threshold of the system. In the process of numerical calculation, the step size of the external excitation frequency ω is 0.1. The results are shown in Figure 2, where the small circles are the numerical results, and the solid linear line denotes the approximate analytical results.

From the analysis of Figure 2, it could be concluded that numerical results and approximate analytical results are in qualitative agreement. The chaotic threshold obtained by the two methods increases first and then decreases with the increase of external excitation frequency. Because Melnikov’s method is a first-order approximation method, there are numerical differences between the two results when solving chaotic motion. Even if this difference exists, it could be useful for reference in analyzing and predicting chaotic motion of similar dynamical systems.

In terms of the values of two methods, there are differences between them. When ω < A, the results obtained by the numerical method are smaller; when A < ω < B, the results obtained by the two methods are in good agreement; when B < ω, the results obtained by the numerical method are larger; when C < ω, it can be seen that the numerical results do not show chaotic motion, only single periodic motion and multiple periodic motion in the system. Moreover, when the value of 1/F tends to be zero, the external amplitude F tends to be infinite; this difference of those two methods means that the Melnikov method is not valid for these values.

Furthermore, it could be concluded that below the threshold curve of chaotic motion, there may be chaotic motion or periodic motion in the system. As we know, the chaos threshold curve calculated by the Melnikov method is a necessary condition but not a sufficient and necessary condition. That is to say that when the parameters of the system satisfy the conditions shown in Equation (21), chaos may occur in the system. When the excitation frequency ω stays constant, chaotic motion is easy to occur when the excitation amplitude is increased; when the excitation amplitude F stays constant, chaotic motion may occur when the excitation frequency ω is within a certain range. If the excitation frequency ω of the system is very large or the excitation amplitude F is very small, chaotic motion would not occur in the system. Chaotic motions mainly exist in the range of lower external excitation frequency ω and higher external excitation amplitude F. There results could provide a certain reference value and theoretical guidance for designing similar fractional-order systems.

3.2. Bifurcation Analysis and the Largest Lyapunov Exponents

In order to analyze the influence of external excitation amplitude F on the dynamic response of the system, as the amplitude F varies, the bifurcation diagrams and corresponding largest Lyapunov exponents of the system under external periodic excitation are calculated and obtained. In order to further study the bifurcation and chaos of the system, four cases, ω = 1.6, ω = 2.3, ω = 2.9, and ω = 3.5, are shown in Figure 3, Figure 4, Figure 5 and Figure 6, respectively.

Lyapunov index is generally used as a general quantitative index to evaluate the divergence of nearby orbits in nonlinear dynamics. When the maximum Lyapunov index is positive, the system is in a chaotic state [40]. The traditional method of determining the maximum Lyapunov index is based on the benetin-wolf algorithm. Detailed numerical simulation methods for the maximum Lyapunov exponent of fractional differential equations and some typical examples can be found in the relevant literature [41,42,43]. From the observation of Figure 3, Figure 4, Figure 5 and Figure 6, it could be found that the bifurcation diagrams are in agreement with the numerical results of the largest Lyapunov exponents graphs.

The bifurcation diagram of system (1) in (F, x) plane and the corresponding largest Lyapunov exponents for ω = 1.6 are given in Figure 3. From Figure 3, it could be found that when the excitation amplitude F < 0.7, the system dynamic response is period-1 motion. For the period-2 motion at 0.7 < F < 0.75, the chaos occurs at F = 0.75; that is to say that the chaotic threshold is 0.75; and the period-3 window appears at 1 < F < 1.3. When 1.3 < F, the system enters period-1 motion and period-doubling bifurcation to chaos for 5 < F. When the excitation amplitude is near F = 4, chaotic motion appears in little scope. Periodic or quasi-periodic windows and chaotic windows at 6 < F < 7.5 also exist. When F > 7.5, the system enters period-1 motion too. Therefore, the motion state of the system alternates between periodic motion and chaotic motion, with the excitation amplitude increasing.

The bifurcation diagram of system (1) in (F, x) plane and the corresponding largest Lyapunov exponents for ω = 1.6 are given in Figure 4. From Figure 4, it could be found that the motion state of the system alternates between periodic motion and chaotic motion with the excitation amplitude increasing, and period-doubling bifurcation behavior of the system does not exist.

The bifurcation diagram of system (1) in (F, x) plane and the corresponding largest Lyapunov exponents for ω = 2.9 are given in Figure 5. From Figure 5, it can be seen that when F < 1.5, the system dynamic response is period-1 motion. When 1.5 < F < 2.7, the system turns from period-doubling bifurcation period-2 and period-4 motion to chaotic motion. When 2.7 < F < 4.2, the system is in the chaotic motion state. Therefore, F = 2.7 is the chaotic threshold at ω = 2.9. This value is the same as the F value of point L in Figure 6b. With the increasing value of F, it could be found that the turn from an inverse period-doubling bifurcation to chaotic motion occurs when F = 3.9. With the increase of F, the system transforms from period-doubling motion to periodic-2 motion and finally to periodic-1 motion.

The bifurcation diagram of system (1) in (F, x) plane and the corresponding largest Lyapunov exponents for ω = 3.5 are given in Figure 6. From Figure 6, it can be seen that with the change of the amplitude of external excitation, period-1 motion and period-doubling motion state occur in the system response and the chaotic motion state does not exist. Moreover, as can be seen from the largest Lyapunov exponents graph, all the largest Lyapunov exponents are less than or equal to zero. This conclusion is the same as that in Figure 2.

3.3. Analysis of the System Dynamical Responses at Some Different Typical Points

From the observation of Figure 3, Figure 4, Figure 5 and Figure 6, it could be found that there exist abundant dynamical phenomena with the variations of different excitation amplitudes F. In order to further analyze these dynamical phenomena, the nonlinear dynamic analysis method of time history, phase portrait, frequency spectrogram, and Poincare maps at different typical points are presented.

Some typical points are shown in Figure 7. When the excitation frequency ω = 1.6, select points A, B, C, and D; when the excitation frequency ω = 2.3, select points E, F, and G; when the excitation frequency ω = 2.9, select points H, G, K, and L. In the following part, each point is presented and analyzed, respectively. From Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18, the four figures are (a) time history, (b) frequency spectrograms, (c) phase portraits, and (d) Poincare maps, respectively.

3.3.1. Numerical Simulation Results (ω = 1.6)

Point A: The system dynamical response is period-1 motion at F = 0.55 N in Figure 8.

Figure 8. The system response at F = 0.55 N. (a) Displacement response. (b) The amplitude–fre−quency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 8. The system response at F = 0.55 N. (a) Displacement response. (b) The amplitude–fre−quency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Point B: The system dynamical response is period-2 motion at F = 0.65 N in Figure 9.

Figure 9. The system response at F = 0.65 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 9. The system response at F = 0.65 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

From Figure 9, it could be found that the Poincare maps have two fixed points, and there is a dominant frequency, which is equal to the excitation frequency (ω = 1.6), half frequency division (ω = 0.8), and weak frequency multiplication in the frequency spectrogram. It is obvious that the system response is period-2 motion in this case.

Point C: The system dynamical response is chaotic motion at F = 0.75 N in Figure 10.

Figure 10. The system response at F = 0.75 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 10. The system response at F = 0.75 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

In Figure 10b, the frequency spectrogram exhibits a strong peak at the fundamental frequency together with a higher and lower frequency broad band. The Poincare maps show a collection of discrete points. It is obvious that the system response is chaotic motion in this case.

Point D: The system dynamical response is period-3 motion at F = 1.1 N in Figure 11.

Figure 11. The system response at F = 1.1 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 11. The system response at F = 1.1 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

3.3.2. Numerical Simulation Results (ω = 2.3)

Point E: The system dynamical response is period-1 motion at F = 1 N in Figure 12.

Figure 12. The system response at F = 1 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 12. The system response at F = 1 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Point F: The system dynamical response is intermittency chaotic motion at F = 1.38 N in Figure 13.

Figure 13. The system response at F = 1.38 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 13. The system response at F = 1.38 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

When F = 1.38, the time history, frequency spectrogram, phase portraits, and Poincare maps are shown in Figure 13a–d, respectively. It is obvious that the system response is intermittency chaotic motion in this case. Intermittency is a typical route to chaotic motion. Intermittent motion is called pseudo-random alternating motion between regular and irregular motion. The time of the irregular motion becomes longer and longer until the system enters the chaotic motion state completely.

Point G: The system dynamical response is chaotic motion at F = 5.5 N in Figure 14.

Figure 14. The system response at F = 5.5 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 14. The system response at F = 5.5 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

3.3.3. Numerical Simulation Results (ω = 2.9)

Point H: The system dynamical response is period-1 motion at F = 1.2 N in Figure 15.

Figure 15. The system response at F = 1.2 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 15. The system response at F = 1.2 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Point J: The system dynamical response is period-2 motion at F = 2 N in Figure 16.

Figure 16. The system response at F = 2 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 16. The system response at F = 2 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Point K: The system dynamical response is period-4 motion at F = 2.5 N in Figure 17.

Figure 17. The system response at F = 2.5 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 17. The system response at F = 2.5 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

The time history, frequency spectrogram, phase portrait, and Poincare maps are shown in Figure 17a–d, respectively. There are four fixed points in the Poincare maps, and frequency spectrogram exhibits a dominant frequency, quarter frequency division, half frequency division, and weak frequency multiplication. It could be found that the system response is period-4 motion in this case.

Point L: The system dynamical response is chaotic motion at F = 3 N in Figure 18.

Figure 18. The system response at F =3 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

Figure 18. The system response at F =3 N. (a) Displacement response. (b) The amplitude–frequency curves. (c) Phase trajectory. (d) Poncaret section diagram.

In this section, the system response for ω = 1.6 (A → B → C → D), ω = 2.3 (E → F → G), and ω = 2.9 (H → J → K → L) are analyzed. When ω = 1.6, the state of the system motion undergoes the transformation: period-1 → period-2 → chaotic motion, and a period-3 window appears in chaotic motion. When ω = 2.3, the state of the system motion undergoes the transformation: period-1 → intermittency chaotic motion → chaotic motion, and intermittency chaotic motion abruptly appears at F = 1.38. When ω = 2.9, the state of the system motion undergoes the transformation: period-1 → period-2 → period-4 → chaotic motion, and an inverse period-doubling bifurcation to period motion with the excitation amplitude F increases. From all the above analyses, it could be concluded that there are two routes to chaotic motion: period-doubling bifurcation and intermittency chaotic motion, and periodic or quasi-periodic windows exist in the region of chaotic motion. Moreover, the critical excitation amplitude of F when the chaotic motion appears is in agreement with the chaotic threshold based on the Melnikov method. With the change of external excitation parameters, the response of the system changes between chaotic motion and periodic motion. These results demonstrate the validity of the chaotic threshold curve obtained based on the Melnikov method in Figure 1.

Accordingly, the dynamic characteristics of the fractional-order time-delayed Duffing system can be achieved by adjusting the external excitation parameters of the system to change the motion state of the system.

4. Effect of System Parameters on the Chaotic Threshold Curve

In this section, the main purpose is to study the effects of fractional-order delayed feedback term parameters (including τ, h, and p), linear stiffness coefficient k₁, and linear viscous damping coefficient c₁ on the necessary condition for chaos in the sense of Smale horseshoes. According to Equation (23), the chaotic threshold curve could be obtained with variations of each parameter. The analysis on the effects of those parameters are presented, respectively, in the following text.

In Figure 19, the influence of the delay parameter τ on the chaotic threshold curve is given. From the observation of Figure 19, it could be found that the region of chaotic motion increases with the increase of time-delay parameter τ. That is to say, the threshold value of generating chaos decreases so that it is easy to find chaotic motion in this case, and the possibility of chaos of the system increases.

In order to further explain the influence of fractional-order time-delay parameter τ on chaotic threshold, a typical case is presented in the following part. The time-delay parameter τ is chosen as 0.1. When the system is at points C, G, and L, the time history and phase portrait of the displacement are shown in Figure 20, Figure 21 and Figure 22, respectively.

Through comparing Figure 10 with Figure 20, Figure 14 with Figure 21, and Figure 18 with Figure 22, it could be found that with the decrease of the fractional-order time-delay parameter τ, the dynamical response of the system would transform from chaotic motion to periodic or period-doubling motion.

In Figure 23, the influence of the fractional-order parameter p on the chaotic threshold curve is given. The graph shows that as the fractional-order p increases, the maximal value of the chaotic threshold curve decreases. That is to say, the possibility of generating chaos in the system decreases. The change of order p of fractional-order differential term mainly affects the peak value of the curve.

In Figure 24, the influence of the fractional-order coefficient h of fractional-order differential term on the chaotic threshold curve is given. The fractional-order coefficient h changes as 0, −0.25, −0.5, and −0.75, respectively. As the absolute value of fractional-order coefficient h increases, the region of system generating chaos decreases, and the possibility of chaos in the system decreases.

In Figure 25a, the influence of the linear stiffness coefficient k₁ on the chaotic threshold curve is given. From the observation of Figure 25, it could be found that as the linear stiffness coefficient k₁ decreases, the maximal value of the chaotic threshold curve increases, and the possibility of chaos in the system increases. Moreover, the excitation frequency corresponding to the maximal value of the chaotic threshold curve also decreases.

In Figure 25b, the influence of the linear viscous damping coefficient c₁ on the chaotic threshold curve is given. From the observation of Figure 25b, it could be found that the increase of linear viscous damping coefficient c₁ would make the threshold value of generating chaos decrease, especially the maximal value of the chaotic threshold curve, and the possibility of generating chaos in the system decreases.

Consequently, the system parameters of time delay τ, fractional-order coefficient h, and fractional order p, linear stiffness coefficient k₁ and linear viscous damping coefficient c₁ all affect the chaotic threshold. These parameters mainly affect the maximum value of curve in Figure 19 and Figure 23–25 but have no influence on the shape and trend of the curve. Since chaotic motion is unpredictable, in the application of practical engineering, by adjusting the parameters of the system, chaotic motion can be avoided.

5. Conclusions

The abundant nonlinear dynamic behaviors, which include bifurcation and chaotic motion of a fractional-order time-delayed Duffing system under harmonic excitation, are investigated in this paper.

Based on the Melnikov theory, the necessary condition for chaos in the sense of Smale horseshoes is obtained. The comparison between the results from the Melnikov method and numerical integration shows that the two methods have good qualitative and quantitative agreements, which demonstrate the validity of approximate analytical results directly. It could be found that the Melnikov method has certain limitations; when the excitation frequency of the system is larger than a certain value, the Melnikov theory is not valid. Furthermore, the bifurcation diagrams and the corresponding largest Lyapunov exponents are presented. Moreover, the phase portraits, Poincare maps, time histories, and frequency spectrograms at some typical points are analyzed. From the analysis of numerical simulations, it could be found that there exist two routes to chaotic motion: period-doubling bifurcation and intermittency chaotic motion, and alternation appears between periodic motion and chaotic motion with the variations of the system parameters. In addition, the effects of the system parameters on the chaotic threshold curve are analyzed in detail. The increase of fractional-order coefficient h, linear viscous damping coefficient c₁, linear stiffness coefficient k₁, and the decrease of time delay τ can suppress chaos. By adjusting these parameters, chaotic motion could be avoided.

In some fields of science or engineering, the analysis of chaotic dynamics is of great significance. The results shown in this paper might have a relevant value to better understanding the complexity of chaotic motion of the fractional-order time-delay Duffing system.