Subharmonic resonance and fixed-range asymptotic stability of the fractional-order SD oscillator

Abstract

The smooth and discontinuous (SD) oscillator is a typical system with strong nonlinear characteristics, and it is widely used in low-frequency vibration isolation and energy harvesting. A fractional damping model denoted by the Caputo model is introduced into the SD oscillator to adjust the property of the secondary resonance and evaluate the stability of the system. The influence of the fractional damping model on the one-third subharmonic resonance and the fixed-range asymptotic stability is studied. Residue theory and the Laplace transform are used to solve the fractional damping model. The amplitude–frequency response function and the existence conditions are derived by means of the averaging method. Lyapunov theory is used to determine the stable criteria of steady-state solutions. The cell-mapping method is ameliorated and used to calculate the fixed-range asymptotic stability of the one-third subharmonic resonance. The main results are as follows: a gap in the excitation amplitude occurs in the region of the existence condition of the one-third subharmonic resonance when the smooth parameter is smaller than 1. The generation of one-third subharmonic resonance is totally avoided for all frequencies when the excitation amplitudes are within the gap. The width of the gap, as well as the amplitude of the one-third subharmonic resonance, is affected by the parameters of the fractional damping term. The fixed-range asymptotic stability of the one-third subharmonic resonance is weak when the fractional damping parameters are large, which indicates a low resistance of the one-third subharmonic resonance to the external disturbance. The tuning effects of the fractional damping model on the one-third subharmonic resonance and fixed-range asymptotic stability are beneficial for the applications of SD oscillators.

Introduction

Geometric nonlinearity, contact nonlinearity, and material nonlinearity are three major nonlinearity problems in the engineering field. The smooth and discontinuous (SD) oscillator is a typical geometric nonlinear system that was proposed by Cao et al. ¹ The SD oscillator has a special negative stiffness property and abundant nonlinear phenomena, such as chaotic motions, Hopf bifurcation, and attractor coexistence. ² This oscillator has been widely used in many fields, such as energy harvesting ³ and low-frequency vibration isolation.^4,5 Hao and Cao ⁶ designed a quasizero-stiffness vibration isolator based on the SD oscillator and demonstrated its local and global bifurcations using double-parameter bifurcation diagrams and the cell-mapping method. Li et al.^7,8 utilized the SD oscillator model to represent a kind of system with dry friction. Yang et al.^3,9,10 investigated the primary resonance phenomenon of a quasizero-stiffness SD oscillator with a time delay under the feedback control of velocity and displacement, as well as the stochastic resonance under harmonic excitation and Gaussian white noise. The researchers found that delayed feedback could enhance the phenomenon of stochastic resonance. ¹¹ The multistable phenomenon of the SD oscillator can be exploited to improve the performance of energy harvesting. ¹² By utilizing the snap-through characteristics of the SD-like system, Geng and Ding, ¹³ Ding and Chen, ¹⁴ and Geng et al. ¹⁵ comprehensively summarized the broadband vibration isolation characteristics and wide application prospects of a snap-through nonlinear energy sink. Based on the snap-through phenomenon generating higher acceleration, Currier et al. ¹⁶ designed a mechanical fish whose spine has an SD-like structure. The strong nonlinear characteristics of the SD oscillator give this kind of system enormous advantages in converting low-frequency vibrations into electricity with high efficiency.^17,18 The sketch model of the SD oscillator is shown in Figure 1.

Figure 1.

Model of the SD oscillator.

In Figure 1, m is the mass, c is the damping coefficient, k is the stiffness coefficient of two oblique springs, L is the original length of the spring, and $l$ is the vertical distance from the fixed points of the springs to the oscillator. The value of $l/L$ is defined as the smooth parameter α, which represents the nonlinear characteristics of the SD oscillator. Each side of the mass is supported by a spring that has two states: precompressed 0 < α = l/L < 1 and prestretched α_> 1. When $\alpha =0$ , the spring is equivalent to a rigid rod, and the system evolves into a discontinuous state. The sketch maps of the restoring force of the SD oscillator with different α are depicted in Figure 2.

Figure 2.

Geometric nonlinear restoring force of: (a) the restoring force and (b) the nonlinear part of the restoring force with $\alpha =0,0.6,1.2$ .

Figure 2 shows the variation in the geometric nonlinear restoring force with the displacement. The restoring force is discontinuous at x = 0 when $\alpha =0$ , and it can be used to study snap-through characteristics. ² In this case, the oscillator can be recognized as being supported by rigid rods. The slope at x = 0 is less than zero, which indicates that the system has a negative stiffness characteristic when the displacement is small and α = 0.6 (the springs are precompressed). No negative stiffness phenomenon occurs when α = 1.2 (the springs are prestretched). The limit of the slope as x approaches $\pm \infty$ is the square of the natural frequency of the SD oscillator, regardless of the value of α. Figure 2(b) shows that the slope is less than −1 at x = 0 when α < 1. When $\alpha =0$ , the nonlinear part of the restoring force is discontinuous, and its limit value is the square of the natural frequency as x approaches $\pm \infty$ , regardless of the value of $\alpha$ . When the smooth parameter equals 0 (the discontinuous state), the system experiences the attractor coexistence phenomenon, and the chaotic saddle point evolves into a chaotic attractor. ¹⁹

Many researchers intrigued by the special restoring force have extensively investigated the nonlinear phenomena of the SD oscillator. To tackle the problem that the nonlinear stiffness term of the SD oscillator cannot be integrated, new analytical methods were proposed: the piecewise approach method in which the nonlinear restoring force of the SD oscillator is substituted as a trilinear function, ²⁰ the four-dimensional averaging method, ²¹ and the generalized elliptic integral transformation method. ²² Moreover, Cao et al. ²³ proposed irrational elliptic functions and a hyperbolic function to study the analytical solution of the system. Many complex nonlinear dynamic behaviors of the SD oscillator have been proven by Tian et al., such as Hopf bifurcation, ²² codimension-2 bifurcation, ²⁴ complex KAM structures, and chaotic behaviors. ²⁵ Santhosh et al. ²⁶ analyzed the numerical solutions of an SD oscillator from the perspective of the frequency domain. Chen et al. ²⁷ found the sufficient and necessary existence conditions and the number of limit circles of grazing bifurcations, as well as the subharmonic solutions and harmonic solutions. ²⁸ The random bifurcation of the SD oscillator under random excitation was also investigated, which demonstrated the size and shape of the random attractor. ²⁹ Chen et al. ³⁰ comprehensively investigated the bifurcation behaviors of the SD oscillator, including pitchfork bifurcation, degenerate Hopf bifurcation, homoclinic bifurcation, double limit cycle bifurcation, Bautin bifurcation, and Bogganov–Takens bifurcation. Wang et al. ³¹ used the stochastic generalized cell mapping method to analyze the stochastic response and the bifurcation of an SD oscillator with fractional differential damping.

Although the abovementioned scholars conducted in-depth investigations on SD oscillators, many researchers used the viscous damping model to represent the energy dissipation characteristics. In fact, the dissipation characteristic cannot be accurately described using a simple viscous model. Fractional-order damping models have been proposed to describe the dissipation characteristics and replace the viscous damping model, which fit well with the experimental curves of many materials.^32,33 Rossikhin and Shitikova ³⁴ used a fractional model to describe the internal friction phenomenon of the free damped vibration of a suspension bridge. The results were in good agreement with the experimental data. Different from viscous damping, researchers found that the fractional order term affected the frequency and the amplitude of the response. ³⁵ Gao and Yu ³⁶ extended the fractional model to the complex domain and studied the chaotic behaviors of the fractional complex Duffing system under symmetric and asymmetric excitation. Sheu et al. ³⁷ studied the chaotic behavior of a Duffing system with fractional damping. Shen et al.^38–41 investigated the influence of fractional terms on different systems and uncovered a series of novel dynamic phenomena. Niu et al. ⁴² used the Melnikov method to study the chaotic threshold of a Duffing system with one fractional term and proposed chaotic detection conditions. Chang et al. ⁴³ used a fractional-order model to represent the hysteretic property of metal rubber damping, and the results fit well with the experimental data.

From the abovementioned investigations, the fractional damping model can represent the energy dissipation characteristics. Moreover, the one-third subharmonic resonance and the fixed-range asymptotic stability, which are important to the application of the SD oscillator in vibration isolation and energy harvesting, are rarely explored. The one-third subharmonic resonance, which typically occurs within the region of vibration isolation, significantly affects the effectiveness and safety of vibration isolation. The efficiency of energy harvesting can be improved by logically adjusting the existence region of the one-third subharmonic resonance. Thus, the influence of the fractional damping term on the one-third subharmonic resonance and fixed-range asymptotic stability of the SD oscillator is studied. In Section 2, the dynamic model of the fractional-order SD oscillator is established. The residue theorem and the Laplace transform are used to deal with the fractional-order term. The amplitude–frequency response function and the existence conditions of one-third subharmonic resonance are obtained. In Section 3, the influence of the fractional parameters on the nonlinear characteristics is analyzed. In Section 4, the stability condition of the steady-state solution is established by Lyapunov theory. The fixed-range asymptotic stability of the steady-state solutions is analyzed by the fractional cell-mapping method. The influence of the fractional parameters on the fixed-range asymptotic stability is explored in detail. The research results are discussed in Section 5.

Approximate analytical solution of the fractional-order SD oscillator

Motion differential equation of the fractional-order SD oscillator

Let ${\omega }_0^2=\sqrt{2k/m}$ and $c/m=2n$ the dynamic differential equation of the SD oscillator with viscous damping ¹ can be expressed as

$$\stackrel{··}{X}\left(t\right)+2n\stackrel{·}{X}\left(t\right)+{\omega }_0^{2}X\left(t\right)(1-\frac{1}{\sqrt{{\alpha }^2+X^2\left(t\right)}})=0$$ (1)

The restoring force of the SD oscillator is defined as

$$F_r={\omega }_0^{2}X\left(t\right)-\frac{{\omega }_0^{2}X\left(t\right)}{\sqrt{{\alpha }^2+X^2\left(t\right)}}$$ (2)

The nonlinear part of equation (2) is defined as

$$F_n=-\frac{{\omega }_0^{2}X\left(t\right)}{\sqrt{{\alpha }^2+X^2\left(t\right)}}$$ (3)

In equation (3), the limitation of $F_n$ with the displacement is ${\omega }_0^2$ , which is one of the features of SD oscillators. The fractional-order model is introduced to replace the viscous damping model in equation (1) to accurately describe the energy dissipation property for investigating the nonlinear response characteristics of the system. The motion differential equation of the fractional-order SD oscillator can be rewritten as

$$m\stackrel{··}{X}(t)+h{\mathrm{D}}^p[X(t)]+2\mathrm{kX}(t)(1-\frac{L}{\sqrt{l^2+X^2(t)}})=F\cos \left(\omega t\right)$$ (4)

where $F$ is the amplitude of external excitation, $\omega$ is the excitation frequency, $h$ is the fractional damping coefficient, $p$ is the fractional damping order which is within the interval $[0,1]$ . As $p$ approaches 0, the fractional term evolves into a linear stiffness term, while the fractional term performs viscous damping when $p$ increases to 1. ${\mathrm{D}}^p[·]$ is the Caputo fractional ³⁸ differential term, and $X$ is the displacement of the oscillator.

The time scale $t_1$ and the length scale $L_1$ are selected as

$$L_1=L,t_1=\sqrt{\frac{m}{k}}$$ (5)

Let $x=X/L_1$ , $\tau =t/t_1$ , $\alpha =l/L$ , $\epsilon \mu ={\mathrm{ht}}_1^{(2-p)}/m$ , ${\omega }_0^2=2{\mathrm{kt}}_1^2/m$ , $\mathrm{\Omega }=\omega t_1$ , and $f={\mathrm{Ft}}_1^2/\left(m{L}_1\right)$ . Equation (4) can be rewritten as

$$\stackrel{··}{x}(\tau )+\epsilon \mu {\mathrm{D}}^p[x(\tau )]+{\omega }_0^{2}x(\tau )-\frac{{\omega }_0^{2}x(\tau )}{\sqrt{{\alpha }^2+x^2(\tau )}}=f\cos \left(\mathrm{\Omega }\tau \right)$$ (6)

The restoring force of equation (6) is supposed to be

$$P={\omega }_0^{2}x(\tau )-\frac{{\omega }_0^{2}x(\tau )}{\sqrt{{\alpha }^2+x^2(\tau )}}$$ (7)

The nonlinear part of equation (7) is defined as

$$F_n=-\frac{{\omega }_0^{2}x(\tau )}{\sqrt{{\alpha }^2+x^2(\tau )}}$$ (8)

Notably, $F_n$ cannot be integrated and cannot be solved by the average method directly. Piecewise functions can be used to mimic $F_n$ but are not convenient for deriving the amplitude–frequency response function. Thus, $F_f=B_1(\alpha )\epsilon {\alpha }_{1}x+B_2(\alpha )\epsilon {\alpha }_2{x}^3$ is used to fit $F_n$ according to the least square method, where $\epsilon {\alpha }_1=-{\omega }_0^2/\alpha$ , $\epsilon {\alpha }_2={\omega }_0^2/(2{\alpha }^3)$ , $B_1(\alpha )$ , and $B_2(\alpha )$ are undetermined constants. The parameters are chosen as $m=5$ and $k=5$ . The fitting results of $F_f$ and three-order Taylor series are shown in Figure 3.

Figure 3.

Fitting results of $F_f$ and three-order Taylor series: (a) $\alpha =0.6$ and (b) $\alpha =1.2$ .

Figure 3 illustrates that $F_f$ fits well with $F_n$ . The third-order Taylor series will produce a large error when $x$ is large. In the follow-up investigation, only the amplitude of one-third subharmonic resonance within the fitting range is analyzed. Next, by substituting $F_f$ into equation (4), we can obtain

$$\begin{array}{c}\stackrel{··}{x}(\tau )+\epsilon \mu {\mathrm{D}}^p[x(\tau )]+{\omega }_0^{2}x(\tau )+B_1(\alpha )\epsilon {\alpha }_{1}x(\tau )\hfill \\ +B_2(\alpha )\epsilon {\alpha }_{2}x(\tau {)}^3=f\cos \left(\mathrm{\Omega }\tau \right)\hfill \end{array}$$ (9)

Subsequently, the amplitude–frequency response function of the one-third subharmonic resonance of the fractional-order SD oscillator is derived by the averaging method. To investigate the steady-state solutions, the following expression is introduced.

$$\frac{{\mathrm{\Omega }}^2}{9}={\omega }_0^2+\epsilon \sigma$$ (10)

where is the tuning parameter that measures the frequency proximity.

By substituting equation (10) into the left side of equation (9), we can obtain

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{··}{x}(\tau )+\epsilon \mu {\mathrm{D}}^p[x(\tau )]+\frac{{\mathrm{\Omega }}^2}{9}x(\tau )-\epsilon \sigma x(\tau )+B_1(\alpha )\epsilon {\alpha }_{1}x(\tau )}\\ \multicolumn{1}{c}{+B_2(\alpha )\epsilon {\alpha }_{2}x(\tau {)}^3=f\cos \left(\mathrm{\Omega }\tau \right)}\end{array}$$ (11)

Equation (11) can be rewritten as

$$\begin{array}{c}\stackrel{··}{x}(\tau )+\frac{1}{9}{\omega }^{2}x(\tau )\hfill \\ =\epsilon \{\sigma x(\tau )-\mu {\mathrm{D}}^p[x(\tau )]-B_1(\alpha ){\alpha }_{1}x(\tau )-B_2(\alpha ){\alpha }_{2}x{(\tau )}^3\}\hfill \\ +f\cos \left(\mathrm{\Omega }\tau \right)\hfill \end{array}$$ (12)

The solution of equation (11) is supposed to be

$$x(\tau )=a\cos \left(\phi \right)+b\cos \left(\mathrm{\Omega }\tau \right)$$ (13)

$$\stackrel{·}{x}(\tau )=-\frac{1}{3}a\mathrm{\Omega }\sin \left(\phi \right)-b\omega \sin \left(\mathrm{\Omega }\tau \right)$$ (14)

where $\phi =\frac{\mathrm{\Omega }}{3}\tau +\theta$ , $b=-\frac{9f}{8{\mathrm{\Omega }}^2}$ .

By substituting equations (13) and (14) into equation (12) and based on the averaging method, we can obtain

$$\stackrel{·}{a}=-\frac{3\sin \left(\phi \right)}{\mathrm{\Omega }}\epsilon \{P_1(a,\theta )+P_2(a,\theta )\}$$ (15)

$$a\stackrel{·}{\theta }=-\frac{3\cos \left(\phi \right)}{\mathrm{\Omega }}\epsilon \{P_1(a,\theta )+P_2(a,\theta )\}$$ (16)

where $P_1(a,\theta )=\sigma (a\cos \left(\phi \right)+b\cos \left(\mathrm{\Omega }\tau \right))-B_1(\alpha ){\alpha }_1(a\cos \left(\phi \right)+b\cos \left(\mathrm{\Omega }\tau \right))-B_2(\alpha ){\alpha }_2{(a\cos \left(\phi \right)+b\cos \left(\mathrm{\Omega }\tau \right))}^3$ $P_2(a,\theta )=-\mu {\mathrm{D}}^p[a\cos \left(\phi \right)+b\cos \left(\mathrm{\Omega }\tau \right)]$ .

Calculating the fractional-order term by using the residue theorem

The residue theorem is introduced into the average method to calculate $P_2(a,\theta )$ . The second part of $P_2(a,\theta )$ is taken as an example. Using the Laplace transform, we can obtain

$$\mathrm{L}[P_{22}(a,\theta )]=-\mu \mathrm{L}\left[{\mathrm{D}}^p[b\cos (\mathrm{\Omega }\tau )]\right]=-\mu b\times \frac{S^{p+1}}{S^2+{\mathrm{\Omega }}^2}$$ (17)

The inverse Laplace transform is performed on equation (17). The residue theorem is used to calculate the inverse Laplace transform, which can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{P_{22}(a,\theta )={\mathrm{L}}^{-1}[-\mu b\times \frac{S^{p+1}}{S^2+{\mathrm{\Omega }}^2}]}\\ =-\mu b\times (\underset{S=\mathrm{i}\omega }{\mathrm{Res}}\left[\frac{S^{p+1}}{S^2+{\mathrm{\Omega }}^2}\right]+\underset{S=-\mathrm{i}\omega }{\mathrm{Res}}\left[\frac{S^{p+1}}{S^2+{\mathrm{\Omega }}^2}\right])\hfill \\ =-\mu b\times ({\frac{S^{p+1}}{2}{\mathrm{e}}^{S\tau }|}_{S=\mathrm{i}\mathrm{\Omega }}+{\frac{S^{p+1}}{2}{\mathrm{e}}^{S\tau }|}_{S=-\mathrm{i}\mathrm{\Omega }})\hfill \\ \multicolumn{1}{c}{=-\mu b\times (\frac{{(\mathrm{i}\mathrm{\Omega })}^p}{2}{\mathrm{e}}^{\mathrm{i}\mathrm{\Omega }\tau }+\frac{{(-\mathrm{i}\mathrm{\Omega })}^p}{2}{\mathrm{e}}^{-\mathrm{i}\mathrm{\Omega }\tau })}\end{array}$$ (18)

where $S=\pm \mathrm{i}\mathrm{\Omega }$ is the singularity of $S^{p+1}/S^2+{\mathrm{\Omega }}^2$ .

To simplify equation (18), the following formula is introduced:

$${\mathrm{i}}^p={\left({\mathrm{e}}^{\frac{\mathrm{i}\pi }{2}}\right)}^p={\mathrm{e}}^{\frac{\mathrm{i}\pi p}{2}}$$ (19)

Equation (19) is performed on equation (18), and we can obtain

$$\begin{array}{c}\multicolumn{1}{c}{P_{22}(a,\theta )=-\mu \mathrm{b}\times (\frac{{\omega }^p}{2}{\mathrm{e}}^{\mathrm{i}(\mathrm{\Omega }\tau +\frac{p\pi }{2})}+\frac{{\omega }^p}{2}{\mathrm{e}}^{-\mathrm{i}(\mathrm{\Omega }\tau +\frac{p\pi }{2})})}\\ =-\mu b{\mathrm{\Omega }}^p\cos (\mathrm{\Omega }\tau +\frac{p\pi }{2})\hfill \end{array}$$ (20)

By the same process, $P_{21}(a,\theta )$ can be expressed as

$$P_{21}(a,\theta )=-\mu a{\left(\frac{\mathrm{\Omega }}{3}\right)}^p\cos (\frac{\mathrm{\Omega }}{3}\tau +\theta +\frac{p\pi }{2})$$ (21)

By summing equations (20) and (21), $P_{21}(a,\theta )$ can be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{P_{2L}(a,\theta )=P_{21}(a,\theta )+P_{22}(a,\theta )}\\ =-\mu a{\left(\frac{\mathrm{\Omega }}{3}\right)}^p\cos (\frac{\mathrm{\Omega }}{3}\tau +\theta +\frac{p\pi }{2})\hfill \\ -\mu b{\mathrm{\Omega }}^p\cos (\mathrm{\Omega }\tau +\frac{p\pi }{2})\hfill \end{array}$$ (22)

Substituting equation (22) into equations (15) and (16), they can be expressed as

$${\stackrel{·}{a}}_R=-\frac{3\sin \left(\phi \right)}{\mathrm{\Omega }}\epsilon \{P_1(a,\theta )+P_{2L}(a,\theta )\}$$ (23)

$$a{\stackrel{·}{\theta }}_R=-\frac{3\cos \left(\phi \right)}{\mathrm{\Omega }}\epsilon \{P_1(a,\theta )+P_{2L}(a,\theta )\}$$ (24)

where $P_{2L}(a,\theta )=-\mu a{\left(\frac{\mathrm{\Omega }}{3}\right)}^p\cos (\phi +\frac{p\pi }{2})-\mu b{\mathrm{\Omega }}^p\cos (\mathrm{\Omega }t+\frac{p\pi }{2})$ .

According to the averaging method, equations (23) and (24) are integrated over the time interval $[0,T]$ . By selecting $T=2\pi$ if $P_i(a,\theta )(i=1,2L)$ are periodic functions, they can be obtained as

$${\stackrel{·}{a}}_R=\frac{3a}{8\mathrm{\Omega }}(-4\times {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\epsilon \mu \sin \left(\frac{p\pi }{2}\right)+3\mathrm{ab}{B}_2(\alpha )\epsilon {\alpha }_2\sin (3\theta ))$$ (25)

$$\begin{array}{c}a{\stackrel{·}{\theta }}_R=\frac{3a}{8\mathrm{\Omega }}(4B_1(\alpha )\epsilon {\alpha }_1+(3a^2+6b^2)B_2(\alpha )\epsilon {\alpha }_2-4\epsilon \sigma \hfill \\ +4\times {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\epsilon \mu \cos \left(\frac{p\pi }{2}\right)+3\mathrm{ab}{B}_2(\alpha )\epsilon {\alpha }_2\cos (3\theta ))\hfill \end{array}$$ (26)

To obtain the steady-state solutions of the system, let $\stackrel{·}{a}=0$ and $a\stackrel{·}{\theta }=0$ . Equations (25) and (26) become

$$\frac{3a}{8\mathrm{\Omega }}(-4\times {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\epsilon \mu \sin \left(\frac{p\pi }{2}\right)+3\mathrm{ab}{B}_2(\alpha )\epsilon {\alpha }_2\sin (3\theta ))=0$$ (27)

$$\begin{array}{l}\frac{3a}{8\mathrm{\Omega }}(4B_1(\alpha )\epsilon {\alpha }_1+\left(3a^2+6b^2\right)B_2(\alpha )\epsilon {\alpha }_2-4\epsilon \sigma +4\times {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\\ \epsilon \mu \cos \left(\frac{p\pi }{2}\right)+3ab{B}_2(\alpha )\epsilon {\alpha }_2\cos (3\theta ))=0\end{array}$$ (28)

Amplitude–frequency response function and existence conditions

By eliminating $\theta$ from equations (27) and (28), the amplitude–frequency response function of the one-third subharmonic resonance is obtained as

$$\begin{array}{c}{\left[4\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\sin \left(\frac{p\pi }{2}\right)\right]}^2\hfill \\ +{[4\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\cos \left(\frac{p\pi }{2}\right)+4B_1(\alpha )\epsilon {\alpha }_1+3(a^2+2b^2)B_2(\alpha )\epsilon {\alpha }_2-4\epsilon \sigma ]}^2\hfill \\ ={(3\mathrm{ab}{B}_2(\alpha )\epsilon {\alpha }_2)}^2\hfill \end{array}$$ (29)

Let $A=a^2$ . Equation (29) can be rewritten in polynomial form as

$$W_1{A}^2+W_{2}A+W_3=0$$ (30)

where $W_1=9{(B_2(\alpha )\epsilon {\alpha }_2)}^2$ , $W_2=24B_1(\alpha )B_2(\alpha )\epsilon {\alpha }_1\epsilon {\alpha }_2+27b^2{(B_2(\alpha )\epsilon {\alpha }_2)}^2-24B_2(\alpha )\epsilon {\alpha }_2\epsilon \sigma +8\times 3^{1-p}{B}_2(\alpha )\epsilon {\alpha }_2\epsilon \mu {\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right)$ , $W_3=4\times 9^{-p}(9^p{(2B_1(\alpha )\epsilon {\alpha }_1+3b^2{B}_2(\alpha )\epsilon {\alpha }_2-2\epsilon \sigma )}^2+4{(\epsilon \mu )}^2{\mathrm{\Omega }}^{2p}$ $+4\times 3^p(2B_1(\alpha )\epsilon {\alpha }_1+3b^2{B}_2(\alpha )\epsilon {\alpha }_2-2\epsilon \sigma )\epsilon \mu {\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right))$

By solving equation (30), the existence conditions for real solutions to the quadratic equation are

$$-\frac{W_2}{2W_1}>0$$ (31)

$$W_2^2-4W_1{W}_3\ge 0$$ (32)

Equations (31) and (32) can be further simplified. Let $B=b^2$ ; then, we can obtain

$$B<\frac{8(-B_1(\alpha )\epsilon {\alpha }_1+\epsilon \sigma -{\left(\frac{\mathrm{\Omega }}{3}\right)}^p\epsilon \mu \cos \left(\frac{p\pi }{2}\right))}{9B_2(\alpha )\epsilon {\alpha }_2}$$ (33)

$$|2Q_{1}B+Q_2|<\sqrt{Q_2^2-4Q_1{Q}_3}$$ (34)

where $Q_1=-567{(B_2(\alpha )\epsilon {\alpha }_2)}^4<0$ , $Q_2=(-16\times 3^{3-p}{(B_2(\alpha )\epsilon {\alpha }_2)}^3(3^p\left(B_1(\alpha )\epsilon {\alpha }_1-\epsilon \sigma \right)+\epsilon \mu {\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right)))>0$ , $Q_3=32\times 9^{1-p}{(B_2(\alpha )\epsilon {\alpha }_2\epsilon \mu {\mathrm{\Omega }}^p)}^2(\cos (p\pi )-1)<0$

Equations (33) and (34) are the existence conditions for the one-third subharmonic resonance of the fractional-order SD oscillator.

Effects of the fractional parameters on the amplitude–frequency response function and the existence conditions

Effects of the fractional parameters when the smooth parameter is smaller than 1

The influence of the fractional damping coefficient $h$ on the existence condition (EC hereinafter) and the amplitude frequency response is studied when the springs are precompressed. The EC and the amplitude–frequency response curve (AFC) with different fractional damping coefficients $h$ are shown in Figures 4 and 5 when the smooth parameter $\alpha =0.6$ . The system parameters are selected as $m=5$ , $k=5$ , $l=0.096$ , and $L=0.16$ .

Figure 4.

Existence conditions when $p=0.5$ with $h$  = 1.5, 1.88, 2.5, 3.

Figure 5.

Amplitude–frequency response curves as $p=0.5$ : (a) $f=0.375$ and (b) $f=5$ .

Figure 4 depicts ECs with different fractional damping coefficients $h$ when the fractional damping order $p=0.5$ . ECs are depicted by solid lines and filled with the same colors as the lines. The one-third subharmonic resonance is generated as the parameters of excitation within the region surrounded by EC. The left line of EC indicates the start frequency of the one-third subharmonic resonance, while the right line indicates the end frequency. The range of EC shrinks as the fractional damping coefficient increases. A novel phenomenon is that a gap in the excitation amplitude $f$ occurs as the fractional damping coefficient increases (the light green area in Figure 4). The gap splits EC into two parts, which indicates that subharmonic resonance cannot be completely avoided by increasing or decreasing the excitation amplitude. The one-third subharmonic resonance is only generated in the high-frequency–high-amplitude region as the gap expands.

Figure 5 depicts AFCs with different fractional damping coefficients when the fractional damping order $p=0.5$ . The stable solutions are plotted as solid lines, while the unstable solutions are plotted as dashed lines. The amplitude of the stable solution decreases, while the start frequency of the one-third subharmonic resonance increases when the fractional damping coefficient rises. These phenomena illustrate that the damping effect of the system is reinforced, and the stiffness is hardened.

Next, the influence of the fractional damping order $p$ on the one-third subharmonic resonance is studied. The EC and AFC with different $p$ are shown in Figures 6 and 7 when the fractional damping coefficient $h=2.2$ .

Figure 6.

Existence conditions when $h=2.2$ : (a) $p=0.1$ , (b) $p=0.4$ , (c) $p=0.7$ , and (d) $p=1.0$ .

Figure 7.

Amplitude–frequency response curves when $h=2.2$ _: (a) $f=0.125$ and (b) $f=3.875$ .

Figure 6 depicts ECs with different fractional damping orders $p$ when the fractional damping coefficient $h=2.2$ . At the beginning, one EC area exists; then, the area shrinks, and the second part of it occurs as $p$ increases, as shown in Figure 6(b). Then, the upper EC area shrinks, while the lower EC area enlarges with the increase in the gap. As the fractional damping order approaches 1, the gap disappears, which results in that the two parts of EC merge, as shown in Figure 6(d). After the gap disappears, this system cannot avoid the one-third subharmonic resonance, but the start frequency can be adjusted by the excitation amplitude. The start frequency increases with the rise in the excitation amplitude.

Figure 7 depicts AFCs with different fractional damping orders when the fractional damping coefficient $h=2.2$ . As shown in Figure 7(a), the amplitude of the stable solution increases, while the start frequency of the one-third subharmonic resonance decreases as the fractional damping order $p$ rises. As shown in Figure 7(b), the start frequency of the one-third subharmonic resonance first increases and then decreases, while the amplitude rises as the fractional damping order $p$ improves. These phenomena illustrate that the damping effect is weakened, and the stiffness is softened when the excitation frequency $\mathrm{\Omega }$ is low. However, the stiffness first hardens and then softens when the excitation frequency $\mathrm{\Omega }$ is high.

Figure 8 depicts ECs with different fractional damping orders when the fractional damping coefficient $h=2.6$ . First, one EC area exists; then, the area shrinks, and the second part of it occurs as $p$ increases, as shown in Figure 8(c). As the fractional damping order continues to increase, the gap starts to decrease, but it does not disappear even if the fractional damping order equals 1. When the fractional damping order is smaller than or equal to 0.4, the one-third subharmonic resonance can be avoided by decreasing the excitation amplitude. However, the excitation amplitude should be constrained within a certain range to avoid the one-third subharmonic resonance after the gap appears.

Figure 8.

Existence conditions when $h=2.6$ _: (a) $p=0.1$ , (b) $p=0.4$ , (c) $p=0.7$ , and (d) $p=1.0$ .

Figure 9 depicts AFCs with different fractional damping orders when the fractional damping coefficient $h=2.6$ . As shown in Figure 9(a), the amplitude of the stable solution increases. The start frequency of the one-third subharmonic resonance decreases as the fractional damping order $p$ increases. As shown in Figure 9(b), the amplitude of the stable solution first decreases and then increases. The start frequency of the one-third subharmonic resonance first increases and then decreases as the fractional damping order $p$ rises. These phenomena illustrate that the damping effect is weakened, and the stiffness is softened when the excitation frequency $\omega$ is low. The damping effect first reinforces and then weakens, while the stiffness first hardens and then softens when the excitation frequency $\omega$ is high.

Figure 9.

Amplitude–frequency response curves when $h=2.6$ _: (a) $f=0.125$ and (b) $f=6.25$ .

Effects of the fractional parameters when the smoothing parameter is larger than 1

When $\alpha =1.2$ (the springs are prestretched), the EC and AFC with different fractional parameters are as shown in Figures 10 and 11.

Figure 10.

Existence conditions and amplitude–frequency response curves as $p=0.3$ , and the excitation amplitude $f$ in (b) is 7.5: (a) existence conditions and (b) amplitude–frequency response curves.

Figure 11.

Existence conditions and amplitude–frequency response curves as $h=1.5$ , and the excitation amplitude $f$ in (b) is 5: (a) existence conditions and (b) amplitude–frequency response curves.

Figure 10 depicts ECs and AFCs with different fractional damping coefficients when the fractional damping order $p=0.3$ . Figure 10(a) indicates that the area of EC shrinks as the fractional damping coefficient increases. Figure 10(b) shows that the amplitude of the stable solution decreases, while the start frequency of the one-third subharmonic resonance increases as the fractional damping coefficient $h$ rises. These phenomena illustrate that the damping effect is reinforced, and the stiffness is hardened. The system with a small excitation amplitude does not generate the one-third subharmonic resonance. Figure 11 depicts the ECs and AFCs with different fractional damping orders $p$ when the fractional damping coefficient $h=1.5$ . Figure 11(a) indicates that the area of EC shrinks as the fractional damping order increases. Notably, the left boundary of EC corresponding to different $p$ overlaps. Under this condition, the system cannot avoid the one-third subharmonic resonance, but the start frequency can be adjusted by the excitation amplitude. The start frequency increases with the rise in the excitation amplitude. Figure 11(b) shows that the amplitude of the stable solution increases, and the start frequency of the one-third subharmonic resonance slightly decreases as the fractional damping order $p$ rises. These phenomena illustrate that the damping effect is weakened, and the stiffness is softened. Under this condition, a gap in the excitation amplitude does not exist.

Stability analysis of the steady-state solution

Stability conditions of the steady-state solution

Section 4 focuses on the stability of the steady-state solution of the one-third subharmonic resonance. Substituting $a=\overline{a}+\mathrm{\Delta }a$ and $\theta =\overline{\theta }+\mathrm{\Delta }\theta$ into equations (25) and (26) yields the stability condition for the steady-state solution. The results can be obtained as

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\mathrm{d}\mathrm{\Delta }a}{\mathrm{d}t}=\mathrm{\Delta }a\times \frac{1}{2}\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^{p-1}\sin \left(\frac{p\pi }{2}\right)}\\ +\mathrm{\Delta }\theta \times \frac{9\overline{a}}{8\mathrm{\Omega }}(-4B_1(\alpha )\epsilon {\alpha }_1-3({\overline{a}}^2+2b^2)B_2(\alpha )\epsilon {\alpha }_2+4\epsilon \sigma -4\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\cos \left(\frac{p\pi }{2}\right))\hfill \end{array}$$ (35)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{d\mathrm{\Delta }\theta }{\mathrm{dt}}=\mathrm{\Delta }a\times \frac{3}{8\mathrm{\Omega }}(-4B_1\left(\alpha \right)\epsilon {\alpha }_1+3{\overline{a}}^2{B}_2\left(\alpha \right)\epsilon {\alpha }_2-6b^2{B}_2\left(\alpha \right)\epsilon {\alpha }_2+4\epsilon \sigma -4\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\cos \left(\frac{p\pi }{2}\right))}\\ -\mathrm{\Delta }\theta \times \frac{3\overline{a}}{2}\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^{p-1}\sin \left(\frac{p\pi }{2}\right)\hfill \end{array}$$ (36)

Equations (27) and (28) are substituted into equations (35) and (36). By removing $\overline{\theta }$ , the characteristic determinate of the system can be obtained as follows.

$$\left|\begin{array}{cc}\lambda -R_1\hfill & -R_2\hfill \\ \multicolumn{1}{c}{-R_3}& \lambda -R_4\hfill \end{array}\right|=0$$ (37)

where $R_1=\frac{1}{2}\epsilon \mu {\left(\frac{\omega }{3}\right)}^{p-1}\sin \left(\frac{p\pi }{2}\right)$ , $R_2=\frac{9\overline{a}}{8\mathrm{\Omega }}(-4B_1(\alpha )\epsilon {\alpha }_1-3({\overline{a}}^2+2b^2)B_2(\alpha )\epsilon {\alpha }_2+4\epsilon \sigma -4\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\cos \left(\frac{p\pi }{2}\right))$ , $R_3=\frac{3}{8\mathrm{\Omega }}(-4B_1\left(\alpha \right)\epsilon {\alpha }_1+3{\overline{a}}^2{B}_2\left(\alpha \right)\epsilon {\alpha }_2-6b^2{B}_2\left(\alpha \right)\epsilon {\alpha }_2+4\epsilon \sigma -4\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\cos \left(\frac{p\pi }{2}\right))$ , $R_4=-\frac{3\overline{a}}{2}\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^{p-1}\sin \left(\frac{p\pi }{2}\right)$

According to the Lyapunov stability theory, the stability condition of the steady-state solution can be obtained by calculating equation (37) as

$$\begin{array}{c}({\overline{a}}^2>\frac{2}{3B_2\left(\alpha \right)\epsilon {\alpha }_2}({(2B_1\left(\alpha \right)\epsilon {\alpha }_1+3b^2{B}_2\left(\alpha \right)\epsilon {\alpha }_2-2\epsilon \sigma )}^2+{\left(2\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^2\right)}^2\hfill \\ {+4(2B_1\left(\alpha \right)\epsilon {\alpha }_1+3b^2{B}_2\left(\alpha \right)\epsilon {\alpha }_2-2\epsilon \sigma )\epsilon \mu {\left(\frac{\mathrm{\Omega }}{3}\right)}^p\cos \left(\frac{p\pi }{2}\right))}^{-2})\wedge (a>\frac{1}{3})\hfill \end{array}$$ (38)

Equation (38) indicates that, in the solution of equation (29), the branch with the larger amplitude is asymptotically stable, while the other branch is unstable when two steady-state solutions appear simultaneously. Although equation (38) represents the asymptotic stability conditions of the steady-state solution, the resistance of the stable steady-state solution to external disturbances cannot be analyzed yet. The assumption is that, if the domain of attraction is narrow, then one small disturbance will lead to the mutation of the steady-state motion, and the system cannot be regarded as stable at this time. The anti-disturbance ability of steady-state solutions is weak.

The range of the stable region, that is, the domain of attraction corresponding to the stable attractor, needs to be studied. Only small disturbances need to be considered. Thus, we need to study the asymptotic stability of the steady-state solution under external perturbation in a certain range, that is, the fixed-range asymptotic stability. Then, the attraction domain of the steady-state solution is studied by the fractional cell-mapping method. Equation (9) is rewritten as an extended state equation set.

$$\{\begin{array}{c}{}_0^C{\mathrm{D}}_{t_n}^{p_1}[x_1(\tau )]=x_2(\tau )\hfill \\ {}_0^C{\mathrm{D}}_{t_n}^{p_2}[x_2(\tau )]=x_3(\tau )\hfill \\ {}_0^C{\mathrm{D}}_{t_n}^{p_3}[x_3(\tau )]=f\cos (\mathrm{\Omega }\tau )-\epsilon \mu {\mathrm{D}}^p[x_1(\tau )]\hfill \\ -{\omega }_0^{2}x(\tau )-B_1(\alpha )\epsilon {\alpha }_1{x}_1(t)-B_2(\alpha )\epsilon {\alpha }_2{x}_1^3(t)\hfill \end{array}$$ (40)

where ${}_0^C{\mathrm{D}}_{t_n}^{p_i}[·](i=1,2,3)$ is the fractional-order term represented by the Caputo form, $t_n$ is the iteration calculation time, and $p_i$ is the order of differentiation.

In the fractional cell-mapping method, the extended state equation method is used to calculate the equation. The iterative cell calculation number is 500, the period of iteration is 1, and the cell number is 90,000. Equation (40) regards the fractional-order term as a new dimension, and thus, it belongs to a three-dimensional system. In the iterative calculation, the results of the fractional-order term will be recorded and used as the initial value of the new iteration. The system parameters in verifying the accuracy of the method are selected as $m=5$ , $k=5$ , $l=0.096$ , $L=0.16$ , $f=5$ , $h=0.8$ , and $p=0.5$ . The iterative calculation step is 0.001, and the period of iteration $t_n$ is 500. The results of the fractional cell-mapping method, the Poincaré mapping, and the time history diagrams are shown in Figure 12.

Figure 12.

Numerical results of equation (40): (a) fractional cell-mapping diagram, (b) Poincaré mapping diagram, (c) time history diagram of the P1 point, and (d) time history diagram of the P3 point.

Figure 12(a) depicts the result of the fractional cell-mapping method. A group of period-3 (P3) points exist in the orange regions, and a period-1 (P1) point exists in the purple region. The orange and purple regions are the attraction domains corresponding to P3 and P1, respectively. Figure 12(b) is the result of the Poincaré mapping method. The time history diagrams corresponding to the P1 and P3 points are shown in Figure 12(c) and (d). The result of the fractional cell-mapping method is consistent with those of other methods, which indicates that the fractional cell-mapping method has good accuracy.

The one-third subharmonic resonance is mainly focused on the dimensionless region $S$

$$S=(0<\stackrel{-}{a}<1.5,1<\stackrel{-}{\theta }<2.5)$$ (41)

In the fractional cell-mapping method, $S$ is divided into at least cells. To measure the fixed-range asymptotic stability of the steady-state solutions, the proportion of the attraction domain $S$ is defined as

$$\mathrm{Por}=\frac{C_D}{C_S}$$ (42)

where $C_D$ is the number of cells in the attraction domain of the steady-state solution and $C_S$ is the number of cells in $S$ . As the proportion $\mathrm{Por}$ approaches 1, the fixed-range asymptotic stability of equation (9) is enhanced, which indicates that the resistance of the one-third subharmonic resonance to external disturbance is enhanced, and EC is difficult to break.

Effects of the fractional coefficient on the fixed-range asymptotic stability

The fixed-range asymptotic stability of the steady-state solution with fractional coefficients $h=$ 2, 3.5, and 4 is studied. The system parameters are selected as $m=5$ , $k=5$ , $l=0.096$ , $L=0.16$ , $p=0.5$ , $f\in [3,7.5]$ , and $\mathrm{\Omega }\in [2.8,6.2]$ . The results are shown in Figures 13–15. Graph (a) in each figure depicts the variation in $\mathrm{Por}$ with $(\mathrm{\Omega },f)$ . Graphs from (b) to (f) are typical results of $\mathrm{Por}$ . The attraction domain of the asymptotic stable steady-state solution is drawn in orange, while another domain in purple represents the steady-state solution for $a=0$ , that is, the steady-state response of the system has the same frequency as the external excitation, and EC is broken. The green region indicates the attraction domain of the unstable solution of the one-third subharmonic resonance. The attraction domain with a wider area implies that the corresponding approximate solution has a strong stability. Therefore, the attraction domain of the unstable solution of one-third subharmonic resonance contains the boundary between the orange and purple regions, which indicates that its stability is weak.

Figure 13.

Distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane and the results of the fractional cell-mapping method at typical points with the fractional damping coefficient $h=3$ _: (a) distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane, (b) $(3.9,7.5),\mathrm{Por}=0.3613$ , (c) $(4.4,7.5),\mathrm{Por}=0.4237$ , (d) $(6.2,7.5),\mathrm{Por}=0.2161$ , (e) $(4.9,7.5),\mathrm{Por}=0.3811$ , and (f) $(3.2,3.0),\mathrm{Por}=0.2899$ .

Figure 14.

Distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane and the results of the fractional cell-mapping method at typical points with the fractional damping coefficient $h=3.5$ _: (a) distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane, (b) $(3.9,7.5),\mathrm{Por}=0.2996$ , (c) $(4.4,7.5),\mathrm{Por}=0.3850$ , (d) $(6.2,7.5),\mathrm{Por}=0.1545$ , (e) $(4.9,7.5),\mathrm{Por}=0.3385$ , and (f) $(3.2,3.0),\mathrm{Por}=0.2382$ .

Figure 15.

Distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane and the results of the fractional cell-mapping method at typical points with the fractional damping coefficient $h=4$ _: (a) distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane, (b) $(3.9,7.5),\mathrm{Por}=0.2054$ , (c) $(4.4,7.5),\mathrm{Por}=0.3395$ , (d) $(6.2,7.5),\mathrm{Por}=0.0752$ , (e) $(4.9,7.5),\mathrm{Por}=0.2897$ , and (f) $(3.2,3.0),\mathrm{Por}=0.1708$ .

Figures 13(a), 14(a), and 15(a) depict the distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane, which indicates that $\mathrm{Por}$ decreases as the fractional damping coefficient $h$ increases. The fixed-range asymptotic stability of the steady-state solution reduces. In the three figures, $\mathrm{Por}$ first increases and then decreases with the rise in the excitation frequency $\omega$ , as marked by solid black arrows. Meanwhile, $\mathrm{Por}$ improves with the enhancement in the excitation amplitude $f$ , as shown by red arrows. Figures 13, 14, and 15(b)–(d) are the results of the fractional cell-mapping method with the increase in the excitation frequency $\omega$ , which indicate that the attraction domain of the steady-state solutions first rises and then decreases. The fixed-range asymptotic stability of the steady-state solutions first increases and then decreases. In Figures 13(e), (f), 14(e), (f), 15(e), and (f), the attraction domain of the steady-state solution increases with the rise in the excitation amplitude $f$ . The fixed-range asymptotic stability of the steady-state solutions increases.

Effects of the fractional differential order on the fixed-range asymptotic stability

The fixed-range asymptotic stability of the steady-state solution with $p=$ 0.7, 1.0 is studied. The system parameters are selected as $m=5$ , $k=5$ , $h=3$ , $l=0.096$ , $L=0.16$ , $f\in [3,7.5]$ , and $\mathrm{\Omega }\in [2.8,6.2]$ . The results are shown in Figures 16 and 17.

Figure 16.

Distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane and the results of the fractional cell-mapping method at typical points with $p=0.7$ : (a) distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane, (b) $(3.8,7.5),\mathrm{Por}=0.2640$ , (c) $(4.1,7.5),\mathrm{Por}=0.3799$ , (d) $(5.05,7.5),\mathrm{Por}=0.276$ , (e) $(4.5,7.5),\mathrm{Por}=0.3599$ , and (f) $(2.85,3),\mathrm{Por}=0.2848$ .

Figure 17.

Distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane and the results of the fractional cell-mapping method at typical points with $p=1$ : (a) distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane, (b) $(3.8,7.5),\mathrm{Por}=0.3164$ , (c) $(4.1,7.5),\mathrm{Por}=0.3649$ , (d) $(5.05,7.5),\mathrm{Por}=0.1701$ , (e) $(4.5,7.5),\mathrm{Por}=0.3090$ , and (f) $(2.85,3),\mathrm{Por}=0.3043$ .

Figures 16(a) and 17(a) depict the distribution of $\mathrm{Por}$ in the $(\mathrm{\Omega },f)$ plane, which indicates that $\mathrm{Por}$ decreases as the fractional damping order $p$ increases. The fixed-range asymptotic stability of the steady-state solution decreases. In the two figures, $\mathrm{Por}$ first increases and then decreases with the rise in the excitation frequency $\omega$ , as marked by solid black arrows. Meanwhile, $\mathrm{Por}$ improves with the enhancement in the excitation amplitude $f$ , as shown by red arrows. Figures 16(b)–(d), 17(a), (b), and (d) are the results of the fractional cell-mapping method with the rise in the excitation frequency $\omega$ , which illustrate that the attraction domain of the steady-state solutions first increases and then decreases. The fixed-range asymptotic stability of the steady-state solutions first increases and then decreases. In Figures 16(e), (f), 17(e), and (f), the attraction domain of the steady-state solution increases with the rise in the excitation amplitude $f$ . The fixed-range asymptotic stability of the steady-state solutions increases.

Conclusion

The one-third subharmonic resonance and the fixed-range asymptotic stability of a fractional-order SD oscillator under harmonic excitations are studied. The Caputo fractional model is introduced into the SD oscillator to represent the dissipation property of energy. Residue theory and the Laplace transform are used to calculate the fractional model. The average method is used to obtain the amplitude–frequency response function of the one-third subharmonic resonance. The cell-mapping method is improved and is used to calculate the fixed-range asymptotic stability. The influence of fractional damping on the nonlinear dynamic characteristics and the fixed-range asymptotic stability of the system is analyzed. The conclusions from the theoretical investigations presented in this paper are as follows:

In the existence region of the one-third subharmonic resonance, a novel amplitude gap exists, which is manipulated by the parameters of the fractional damping when the smooth parameter is smaller than 1. The existence region is divided into two parts because of a gap: one part is in the low-frequency region, and the other is in the high-frequency region. If the amplitude of the external excitation is within the gap, then the one-third subharmonic response is totally banned. The width of the gap can be enlarged with the increase in the fractional damping coefficient. However, the width shrinks when the fractional damping order is close to 0 or 1.

The amplitude of the one-third subharmonic resonance is restrained when the fractional damping coefficient increases. When the smooth parameter is larger than 1, the amplitude in the low-frequency part of the existence region is reinforced as the fractional damping order increases, while in the high-frequency region, the amplitude is impaired.

The fixed-range asymptotic stability of the one-third subharmonic resonance decreases, which indicates that the resistance of the one-third subharmonic resonance to the external disturbance degrades as the parameters of the fractional damping term increase. The fixed-range asymptotic stability is enhanced when the amplitude of the excitation is large and the excitation frequency is far from the boundary of the region.

Through the abovementioned theoretical investigations, the influence of fractional damping on the one-third subharmonic resonance and the fixed-range asymptotic stability of the SD oscillator can be understood more clearly. These investigations can be utilized to ameliorate the applications involving SD oscillators such as energy harvesters and low-frequency vibration isolators.