The Nonlinear Dynamics Characteristics and Snap-Through of an SD Oscillator with Nonlinear Fractional Damping

Abstract

A smooth and discontinuous (SD) oscillator is a typical multi-stable state system with strong nonlinear properties and has been widely used in many fields. The nonlinear dynamic characteristics of the system have not been thoroughly investigated because the nonlinear restoring force cannot be integrated. In this paper, the nonlinear restoring force is represented by a piecewise nonlinear function. The equivalent coefficients of fractional damping are obtained with an orthogonal function. The influence of fractional damping on the transition set, the amplitude–frequency response and the snap-through of the SD oscillator are analyzed. The conclusions are as follows: The nonlinear piecewise function accurately mimics the nonlinear restoring force and maintains a nonlinearity property. Fractional damping can significantly affect the stiffness and damping property simultaneously. The equivalent coefficients of the fractional damping are variable with regard to the fractional-order power of the excitation frequency. A hysteresis point, a bifurcation point, a frequency island, pitchfork bifurcations and transcritical bifurcations were discovered in the small-amplitude resonant region. In the non-resonant region, the increase in the fractional parameters leads to the probability of snap-through declining by increasing the symmetry of the attraction domain or reducing the number of stable states.

1. Introduction

Geometric nonlinearity is one of the three major nonlinear problems in engineering applications. Many engineering applications, such as vehicle suspension and gear driving systems, typically exhibit nonlinear characteristics. Nonlinear systems are made equivalent to linear systems to facilitate analytical calculations. In this process, the accurate nonlinear dynamics behavior of the nonlinear system cannot be obtained. If the system has been working at high intensity for a long time, the model errors caused by this equivalent linearization will accumulate. Therefore, establishing an accurate geometric nonlinear system model and analyzing the nonlinear dynamic behavior of the model can ensure the safe and stable operation of the system.

The smooth and discontinuous (SD) oscillator is a typical geometric nonlinear system. The system, which consists of one concentrated mass and two slanted springs, was simplified by a simply supported beam model [1]. The prototypical model of an SD oscillator is shown in Figure 1. This model has been widely used to study snap-through phenomena [2], energy-harvesting mechanisms [3], quasi-zero vibration isolators [4] and applications in many other fields.

In Figure 1, $X(t)$ is the displacement of the oscillator, $m$ is the mass, $c$ is the damping coefficient, $k$ is the stiffness of the springs, $l$ is the constant distance and $L$ is the original length of the spring. The dynamics model of the prototypical SD oscillator is:

$$m\ddot{X}(t)+c\dot{X}(t)+2kX(t)\left(1-\frac{L}{\sqrt{X{(t)}^2+l^2}}\right)=0$$

(1)

This system (1) is made dimensionless by letting $\alpha =l/L$ ($\alpha$ is the smooth parameter) ${\omega }_0^2=2k/m$, $2\zeta =c/m$ and $x=X(t)/L$:

$$\ddot{x}+2\zeta \dot{x}+{\omega }_0^{2}x\left(1-\frac{1}{\sqrt{x^2+{\alpha }^2}}\right)=0$$

(2)

where the nonlinear restoring force is $F_{re}={\omega }_0^{2}x\left(1-1/\sqrt{x^2+{\alpha }^2}\right)$.

Notably, when the original length $L$ is greater than the distance $l$ ($0<\alpha <1$), the oscillator has two stable positions (position 1 and position 3, as shown in Figure 1b). When the displacement or the velocity of the oscillator is smaller than the threshold, the oscillator vibrates at a small amplitude in stable position 1 or 3. If the displacement or the velocity of the oscillator exceeds the threshold (the threshold is the boundary of the attract domain, which will be discussed in Section 4), the oscillator will quickly jump from stable position 1 to stable position 3, or from 3 to 1. These stable position changes are called “snap-through”. However, when the original length $L$ is smaller than the distance $l$ ($\alpha >1$), the oscillator has only one stable position (position 2). Thus, in this condition, the oscillator cannot change its stable position.

The dimensionless system parameters are selected as ${\omega }_0^2=2$, $\alpha$ = 0, 0.5, 1 and 1.5. Diagrams of the nonlinear restoring force $F_{re}$ and the stiffness of the SD oscillator are shown in Figure 2.

Figure 2a shows the nonlinear restoring force with different smooth parameters, and Figure 2b shows the stiffness of the SD oscillator. In Figure 2a, the nonlinear restoring force has a jump point at $x=0$ when $\alpha =0$. In Figure 2b, the stiffness is negative in the neighborhood of $x=0$ ($x_1<x<x_2$) when the smooth parameter is 0.5. The slope of $F_{re}$, i.e., the stiffness, approaches ${\omega }_0^2$ when $\left|x\right|$ approaches infinity. These phenomena indicate that the system is discontinuous when $\alpha =0$, and it has a negative stiffness property when the smooth parameter is in the interval $0<\alpha <1$. The system has a quasi-zero stiffness property when $\alpha =1$. Additionally, its stiffness is positive when $\alpha >1$.

SD oscillators have been widely used in many fields due to their special nonlinear properties. Researchers have conducted many in-depth studies on SD oscillators. Avramov et al. [5] connected an SD-like truss with a spring–mass system and investigated snap-through and nonlinear vibration characteristics. Brennan et al. [6] used an SD-like model with snap-through characteristics to represent the wings of dipteran insects. Their results show that the snap-through phenomenon resulted in a velocity jump in the system. When the operating frequency was below the natural frequency, the system had a great advantage in flight. Ichiro et al. [7] studied the nonlinear characteristics of an elastic–plastic system using a multi-folding microstructure system composed of SD structures. Waite et al. [8] investigated the attractor coexistence and the competing resonance of a truss with the SD structure. Hao et al. [9,10] designed a quasi-zero-stiffness vibration isolator based on an SD oscillator. Currier et al. [11] designed a mechanical fish that had an SD-like spine structure. The structure provided sufficient acceleration for the fish to swim. SD oscillators are also widely applied in nonlinear energy sinks [12,13] and energy harvesting [14,15]. Myongwon et al. [16] demonstrated a 1D lattice of bistable elements. The lattices exhibited energy-harvesting capabilities from transition waves, and energy was transmitted in the form of waves. Tan et al. [17] designed a triboelectric nanogenerator constituted by a bow-type Teng with a snap-through phenomenon. The equipment converted low-frequency vibrations into electricity with a high efficiency.

In terms of nonlinear dynamics, the trilinear piecewise function method [18] and the elliptic function method [19] have been proposed for investigating the nonlinear dynamics characteristics of SD oscillators. Subsequently, Tian et al. [20] studied the codimension-two bifurcation and the Hopf bifurcation of an SD oscillator. Cao et al. [21] analyzed the limit case response of an SD oscillator. In this case, the SD oscillator lost hyperbolicity due to a discontinuous characteristic that was different from a standard double-well system. Shen et al. [22] studied the bifurcation characteristics of symmetric subharmonic orbits and asymmetric subharmonic orbits of a discontinuous SD oscillator. Han et al. [23] used singularity theory to obtain the transition set of a coupled SD oscillator. The chaos threshold was obtained using the Melnikov method. [24] The multiple buckling and the codimension-three bifurcation of an SD oscillator were analyzed [25,26]. Tian et al. [27] studied the chaos threshold of an asymmetric SD oscillator subjected to constant excitation with the topological equivalence method. Yue et al. [28] investigated the random bifurcation of an SD oscillator based on the generalized cell-mapping method. In terms of the frequency domain, Santhosh et al. [29] found a chaotic solution for an SD oscillator caused by the symmetry breaking. Chen et al. [30,31] studied the global bifurcation characteristics of a discontinuous SD oscillator that was equivalent to a Filippov system. Wang et al. [32] represented fractional damping using a Markov chain. The random P bifurcation of an SD oscillator with fractional damping was studied. Their results show that both the smooth parameter and the excitation amplitude induced a random P-bifurcation of the SD oscillator. Through experiments, Chang et al. [33] proved that the periodic motion of an SD oscillator is highly sensitive to the initial displacement and the smooth parameter.

The nonlinear dynamic behaviors of SD oscillators with integral damping have been thoroughly studied. However, two limitations exist in the current research. One limitation is that the Taylor series expansion method is used to mimic the nonlinear restoring force in SD oscillators. The calculation accuracy of the equivalence transformation is acceptable when displacement responses are small. However, this can produce large errors if displacement responses are large. The other limitation is that the integer-damping model is used to represent the energy dissipation. However, in terms of long-term dynamic behaviors, the memory properties and frequency dependence characteristics are not well revealed. It is proven that these properties can be more accurately revealed by the fractional damping model [34,35]. Caputo et al. [36] proposed a fractional-order model to describe the dissipation characteristics, which correlate well with the experimental results of many materials. Padovan et al. [37] investigated the nonlinear response characteristics of a Duffing system with fractional damping. Chang et al. [38] proposed a nonlinear fractional damping model based on the Caputo fractional model [39]. The nonlinear fractional damping model described the energy-dissipation propertoes of metal rubber damping, and the theoretical solution fit well with the results of dynamic load experiments. Their results show that the fractional term simultaneously influenced the frequency and the amplitude of the response. However, research on the vibration characteristics of SD oscillators with fractional damping is not sufficient.

Combined with the current research status, it is necessary to study the nonlinear characteristics of SD oscillators with fractional damping. Therefore, in this article, an equivalent piecewise function is proposed for accurately representing the nonlinear restoring force of an SD oscillator. Moreover, a nonlinear fractional damping model is introduced into the SD oscillator to represent the energy dissipation and viscoelastic properties of a metal rubber damping. This article mainly focuses on the nonlinear amplitude–frequency response and transition characteristics in the resonant region and snap-through phenomena in the non-resonant region of the system. In Section 1, an introduction to this paper is presented. In Section 2, the nonlinear piecewise function is used to approximate the nonlinear restoring force of the SD oscillator. The equivalent damping coefficient and the equivalent stiffness coefficient of fractional damping are derived using the energy equivalent method. The piecewise differential equation of the SD oscillator with nonlinear fractional damping is obtained. In Section 3, the amplitude–frequency response function, the stable criteria and the transition set of the SD oscillator are obtained through calculations. In Section 4, the influence of fractional damping parameters on the nonlinear dynamic characteristics is analyzed in detail. The conclusions are provided in Section 5.

2. The Model of an SD Oscillator with Nonlinear Fractional Damping

The differential equation of the SD oscillator with nonlinear fractional damping is:

$$m\left(\ddot{X}(t)+\ddot{\overline{A}}\right)+c\dot{X}(t)+2k_{0}X(t)\left(1-\frac{L}{\sqrt{X{(t)}^2+l^2}}\right)+k_{1}X(t)+k_{3}X{(t)}^3+h{\mathrm{D}}^p[X(t)]=0$$

(3)

where $\overline{A}=A\cos \left(\omega t\right)$ is the external displacement excitation, $A$ is the excitation amplitude, $\omega$ is the excitation frequency and $X(t)$ is the displacement. $F_{frac}=k_{1}X(t)+k_{3}X{(t)}^3+h{\mathrm{D}}^p[X(t)]$ is the damping force of a metal rubber damping [38], $k_1$ is the linear stiffness coefficient, $k_3$ is the nonlinear stiffness coefficient, $h$ is the fractional coefficient, $p$ is the fractional order and $h{\mathrm{D}}^p[X(t)]$ is the Caputo fractional model [39,40]. The reason the Caputo model is selected is because its initial condition is the same form as the integer-order differential equation. It is convenient for addressing the problem of initial value in dynamics [41].

The new length scale $L_1$ and new time scale $t_1$ are selected as:

$$L_1=L,t_1=\sqrt{\frac{m}{k_1}}$$

(4)

The following dimensionless transformances are determined:

$$\begin{gathered} x(\tau )=\frac{X(t)}{L},\tau =\frac{t}{t_1},2\zeta =\frac{c{t}_1}{m},K_0=\frac{2k_0{t}_1^2}{m},K_1=\frac{k_1{t}_1^2}{m},K_3=\frac{k_3{t}_1^2{L}^2}{m}, \\ H=\frac{h{t}_1^2}{m},\alpha =\frac{l}{L},\mathrm{\Omega }=\omega t_1,F=\frac{A{\mathrm{\Omega }}^2}{L} \end{gathered}$$

By substituting the above dimensionless transformances into Equation (3), we can rewrite the equation as:

$$\ddot{x}(\tau )+2\zeta \dot{x}(\tau )-K_0{F}_n+K_{1}x(\tau )+K_{3}x{(\tau )}^3+H{\mathrm{D}}^p[x(\tau )]=F\cos (\mathrm{\Omega }\tau )$$

(5)

The nonlinear restoring force of Equation (5) is defined as:

$$F_n=-x(\tau )\left(1-\frac{1}{\sqrt{x{(\tau )}^2+{\alpha }^2}}\right)$$

(6)

Equation (6) is nonintegral; therefore, the amplitude–frequency response function of the system cannot be obtained. The nonlinear piecewise function $P_n$ is introduced to obtain the equivalent function of Equation (6):

$$P_n=\left\{\begin{array}{ll}-x(\tau )+B_1,\hfill & x(\tau )<-x_0\hfill \\ B_{2}x{(\tau )}^3+B_{3}x{(\tau )}^2+B_{4}x(\tau ),\hfill & -x_0\le x(\tau )<x_0\hfill \\ -x(\tau )+B_5,\hfill & x(\tau )\ge x_0\hfill \end{array}\right.$$

(7)

where $B_i\left(i=1,2,3,4,5\right)$ is undetermined constants and $x_0\left(x_0>0\right)$ is the piecewise point of Equation (7).

$F_n^{\prime }$ is the first derivative with respect to $\tau$. Let $F_n^{\prime }=0$, and two extreme points $({\tilde{x}}_{1,2},{\tilde{y}}_{1,2})$ of Equation (6) can be obtained:

$$({\tilde{x}}_{1,2},{\tilde{y}}_{1,2})=\left(\pm \sqrt{{\alpha }^{4/3}-{\alpha }^2},\pm \sqrt{{\alpha }^{4/3}-{\alpha }^2}\left(1-\frac{1}{{\alpha }^{2/3}}\right)\right)$$

(8)

Let the middle segment of Equation (7) satisfy Equation (8), and we can obtain:

$$\left\{\begin{array}{l}{B}_2=\frac{{\alpha }^{4/3}-{\alpha }^{2/3}}{2\left(-1+{\alpha }^{2/3}\right){\alpha }^{8/3}}\hfill \\ B_3=0\hfill \\ B_4=\frac{3\left(-1+{\alpha }^{2/3}\right)}{2{\alpha }^{2/3}}\hfill \end{array}\right.$$

(9)

Next, calculate the piecewise point $x_0$ and the corresponding function $y_0$ to obtain the constants $B_1$ and $B_5$. Obviously, $\underset{x\to \infty }{\lim }{F}_n^{\prime }=-1$, and the following can be set:

$$P_n^{\prime }(x_0)=-1$$

(10)

When $x=x_0$, the piecewise point $x_0$ and the function $y_0$ of Equation (7) can be obtained:

$$x_0=\frac{\sqrt{\frac{1}{3}-\frac{1}{{\alpha }^{2/3}}}}{\sqrt{\frac{-{\alpha }^{4/3}+{\alpha }^{2/3}}{\left(-1+{\alpha }^{2/3}\right){\alpha }^{8/3}}}}$$

(11)

$$y_0=\frac{\sqrt{1-\frac{3}{{\alpha }^{2/3}}}\left(-4{\alpha }^{2/3}+3{\alpha }^{2/3}\right)}{3\sqrt{3}{\alpha }^{10/3}{\left(\frac{-{\alpha }^{4/3}+{\alpha }^{2/3}}{\left(-1+{\alpha }^{2/3}\right){\alpha }^{8/3}}\right)}^{3/2}}$$

(12)

Let the first and the third parts of Equation (7) satisfy the following conditions:

$$P_n\left(x_0\right)=y_0$$

(13)

$$P_n\left(-x_0\right)=-y_0$$

(14)

We can obtain $B_1$ and $B_5$:

$$\left\{\begin{array}{c}{B}_1=\frac{\sqrt{1-\frac{3}{{\alpha }^{2/3}}}\left(-4{\alpha }^{4/3}+{\alpha }^{2/3}\left(3+{\alpha }^{4/3}\right)\right)}{3\sqrt{3}\left(-1+{\alpha }^{2/3}\right){\alpha }^{10/3}{\left(\frac{-{\alpha }^{4/3}+{\alpha }^{2/3}}{\left(-1+{\alpha }^{2/3}\right){\alpha }^{8/3}}\right)}^{3/2}}\\ B_5=\frac{\sqrt{1-\frac{3}{{\alpha }^{2/3}}}\left(4{\alpha }^{4/3}-{\alpha }^{2/3}\left(3+{\alpha }^{4/3}\right)\right)}{3\sqrt{3}\left(-1+{\alpha }^{2/3}\right){\alpha }^{10/3}{\left(\frac{-{\alpha }^{4/3}+{\alpha }^{2/3}}{\left(-1+{\alpha }^{2/3}\right){\alpha }^{8/3}}\right)}^{3/2}}\end{array}\right.$$

(15)

By substituting Equations (9) and (15) into Equation (7), the nonlinear piecewise function is derived. By selecting $\alpha =$ 0.1, 0.5 and 0.9, the results of $F_n$, $P_n$, a linear piecewise function and a three-order Taylor series-equivalent function, are depicted in Figure 3.

Figure 3 shows that using the three-order Taylor series-equivalent function to approach $F_n$ is not effective because the deviation between curves for three-order Taylor series and $F_n$ sharply increases. The linear piecewise function and $P_n$ are in good agreement with $F_n$. However, a larger error is generated at the extreme points of $F_n$, and the nonlinear characteristics of the system are eliminated in the interval $-x_0<x<x_0$, when the linear piecewise function is used to represent $F_n$. By analyzing the curve of $P_n$, it is found that $P_n$ correlates well with $F_n$, and the nonlinear characteristics are maintained. Thus, many novel nonlinear phenomena have been discovered, such as the hysteresis point, the bifurcation point and transcritical bifurcation, which will be discussed in Section 4.

By substituting Equation (7) into Equation (5) to replace Equation (6), the differential equation of the SD oscillator with nonlinear fractional damping is written as:

$$\ddot{x}(\tau )+2\zeta \dot{x}(\tau )-K_0{P}_n+K_{1}x(\tau )+K_{3}x{(\tau )}^3+H{\mathrm{D}}^p[x(\tau )]=F\cos (\mathrm{\Omega }\tau )$$

(16)

The Caputo model [39] is used to represent the viscoelastic characteristics of nonlinear fractional damping, i.e., the fractional term in Equation (16):

$${\mathrm{D}}^p[X]=\frac{1}{\mathrm{\Gamma }(1-p)}{\displaystyle {\int }_0^t\frac{X^{\prime }(u)}{{\left(t-u\right)}^p}\mathrm{d}u}$$

(17)

where $\mathrm{\Gamma }(\cdot )$ is a gamma function.

When setting $H{\mathrm{D}}^p[x(\tau )]=c_{eq}\dot{x}(\tau )+k_{eq}x(\tau )$, $\sin (\phi )$ and $\cos (\phi )$ are multiplied on both sides. Both sides of the equations are integrated in one vibration period because the sum of the energy dissipated by damping and the energy stored by springs is a constant.

$$\left\{\begin{array}{l}\frac{H}{\mathrm{\Gamma }(1-p)}\underset{\mathrm{T}\to \infty }{\lim }\frac{1}{T}{\int }_0^T{\int }_0^t\frac{-a\omega \sin (\omega u+\theta )}{{\left(\tau -u\right)}^p}\mathrm{d}u\sin (\omega \tau +\theta )\mathrm{d}t=\oint c_{eq}\left(a\cos \phi \right)\prime \times \sin \phi \mathrm{d}x\hfill \\ \frac{H}{\mathrm{\Gamma }(1-p)}\underset{\mathrm{T}\to \infty }{\lim }\frac{1}{T}{\int }_0^T{\int }_0^t\frac{-a\omega \sin (\omega u+\theta )}{{\left(\tau -u\right)}^p}\mathrm{d}u\cos (\omega \tau +\theta )\mathrm{d}t=\oint k_{eq}a\cos \phi \times \cos \phi \mathrm{d}x\hfill \end{array}\right.$$

(18)

where $x(\tau )=a\cos \left(\mathrm{\Omega }\tau +\theta \right)=a\cos \phi$.

Based on the orthogonality of trigonometric functions, the equivalent stiffness coefficient $k_{eq}$ and the equivalent damping coefficient $c_{eq}$ can be obtained as:

$$\left\{\begin{array}{l}{k}_{eq}=H{\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right)\hfill \\ c_{eq}=H{\mathrm{\Omega }}^{p-1}\sin \left(\frac{p\pi }{2}\right)\hfill \end{array}\right.$$

(19)

It can be found that the fractional order $p$ is the $p\mathrm{th}$ power or the $(p-1)\mathrm{th}$ power of the excitation frequency $\mathrm{\Omega }$ in Equation (19). Thus, $c_{eq}\dot{x}(\tau )+k_{eq}x(\tau )$ is a variable coefficient function with different values of power $\mathrm{\Omega }$ because the response frequency in $x(\tau )$ and $\dot{x}(\tau )$ is the same as the excitation frequency. Additionally, the fractional order $p$ has different influences on $k_{eq}$ and $c_{eq}$ because of different values of the power. In terms of the vibration theory, fractional order $p$ affects the stiffness property and the damping property of the system simultaneously.

Fractional order $p$ is selected as 0, 0.3, 0.7 and 1. The results of the equivalent coefficients are depicted in Figure 4.

Figure 4a shows the results of the equivalent stiffness coefficient $k_{eq}$. Figure 4b shows the results of the equivalent damping coefficient $c_{eq}$. Figure 4 illustrates that $k_{eq}$ increases with increases in $\mathrm{\Omega }$, while $c_{eq}$ decreases. When $\mathrm{\Omega }$ is above ${\mathrm{\Omega }}_1=e^{\pi \tan (p\pi /2)/2}$, $k_{eq}$ first increases, then decreases as $p$ increases (the variation tendency is marked by black arrows), while $k_{eq}$ monotonically decreases when $\mathrm{\Omega }$ is smaller than ${\mathrm{\Omega }}_1=e^{\pi \tan (p\pi /2)/2}$ (the variation tendency is marked by red arrows). When $\mathrm{\Omega }$ is above ${\mathrm{\Omega }}_2=e^{-\pi \cot (p\pi /2)/2}$, $c_{eq}$ monotonically increases as $p$ increases (the variation tendency is marked by red arrows), while $c_{eq}$ first increases, then decreases when $\mathrm{\Omega }$ is smaller than ${\mathrm{\Omega }}_2=e^{-\pi \cot (p\pi /2)/2}$ (the variation tendency is marked by black arrows). When $p$ equals 0, $k_{eq}$ and $c_{eq}$ equal $H$ and 0, respectively. When $p$ equals 1, $k_{eq}$ and $c_{eq}$ equal 0 and $H$, respectively. The above results illustrate that the fractional order and excitation frequency have a significant influence on the stiffness property and the damping property of the system simultaneously.

Substituting Equations (18) and (19) into Equation (16), the differential equation of the SD oscillator with a nonlinear fractional damping can be written as:

$$\ddot{x}(\tau )+\left(2\zeta +H{\mathrm{\Omega }}^{p-1}\sin \left(\frac{p\pi }{2}\right)\right)\dot{x}(\tau )-K_0{P}_n+\left(K_1+H{\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right)\right)x(\tau )+K_{3}x{(\tau )}^3=F\cos \left(\mathrm{\Omega }\tau \right)$$

(20)

3. The Bifurcation of the Amplitude–Frequency Response and the Stability Conditions of Steady-State Solutions

3.1. The Amplitude–Frequency Response Function of the Primary Resonance

The solution of Equation (20) is supposed as:

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

(21)

By substituting Equation (21) into Equation (20), and based on the average method, we can obtain:

$a<x_0$:

$$\left\{\begin{array}{l}\begin{array}{l}\dot{a}=-\frac{1}{2\mathrm{\Omega }}\left[2a\zeta \mathrm{\Omega }+aH{\mathrm{\Omega }}^p\sin \left(\frac{p\pi }{2}\right)+F\sin \left(\theta \right)\right]\hfill \\ a\dot{\theta }=\frac{1}{2\mathrm{\Omega }}\left[\frac{a}{4}\left(3a^2\left(B_2{K}_0+K_3\right)+4\left(B_4{K}_0+K_1-{\mathrm{\Omega }}^2\right)\right.\right.\hfill \end{array}\\ \left.\left.+4H{\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right)\right)-F\cos \left(\theta \right)\right]\end{array}\right.$$

(22)

$a>x_0$:

$$\left\{\begin{array}{l}\dot{a}=-\frac{1}{2\mathrm{\Omega }}\left[2a\zeta \mathrm{\Omega }+aH{\mathrm{\Omega }}^p\sin \left(\frac{p\pi }{2}\right)+F\sin \left(\theta \right)\right]\\ a\dot{\theta }=\frac{1}{2\mathrm{\Omega }}\left[\frac{K_0}{\pi }\left(-2B_1+x_0\right)\sqrt{1-{\left(\frac{x_0}{a}\right)}^2}+\frac{K_0}{\pi }\left(2B_5+x_0\right)\sqrt{1-{\left(\frac{x_0}{a}\right)}^2}\right.\\ -\frac{K_0{x}_0}{2\pi }\sqrt{1-{\left(\frac{x_0}{a}\right)}^2}\left(3a^2{B}_2+4B_4+2B_2{x}_0^2\right)+\frac{2}{\pi }a{K}_0\mathrm{acos}\left(\frac{x_0}{a}\right)\\ +\frac{a}{2\pi }{K}_0\mathrm{asin}(\frac{x_0}{a})\left(3a^2{B}_2+4B_4\right)+a{K}_1+\frac{3}{4}{a}^3{K}_3-a{\mathrm{\Omega }}^2\\ \left.+aH{\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right)-F\cos \left(\theta \right)\right]\end{array}\right.$$

(23)

By eliminating $\theta$, the amplitude–frequency response functions can be obtained as:

$a<x_0$:

$$\begin{array}{l}{\left(2a\zeta \mathrm{\Omega }+aH{\mathrm{\Omega }}^p\sin \left(\frac{p\pi }{2}\right)\right)}^2+\left(\frac{a}{4}\left(3a^2\left(B_2{K}_0+K_3\right)\right.\right.\\ {\left.\left.+4\left(B_4{K}_0+K_1-{\mathrm{\Omega }}^2\right)+4H{\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right)\right)\right)}^2=F^2\end{array}$$

(24)

$a>x_0$:

$$\begin{array}{l}{\left(2a\zeta \mathrm{\Omega }+aH{\mathrm{\Omega }}^p\sin \left(\frac{p\pi }{2}\right)\right)}^2+\left(\frac{K_0}{\pi }\left(-2B_1+x_0\right)\sqrt{1-{\left(\frac{x_0}{a}\right)}^2}\right.\\ +\frac{K_0}{\pi }\left(2B_5+x_0\right)\sqrt{1-{\left(\frac{x_0}{a}\right)}^2}-\frac{K_0{x}_0}{2\pi }\sqrt{1-{\left(\frac{x_0}{a}\right)}^2}\left(3a^2{B}_2+4B_4+2B_2{x}_0^2\right)\\ {\left.+\frac{2}{\pi }a{K}_0\mathrm{acos}\left(\frac{x_0}{a}\right)+\frac{a}{2\pi }{K}_0\mathrm{asin}\left(\frac{x_0}{a}\right)\left(3a^2{B}_2+4B_4\right)+a{K}_1+\frac{3}{4}{a}^3{K}_3\right)}^2\\ {\left.-a{\mathrm{\Omega }}^2+aH{\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right)\right)}^2=F^2\end{array}$$

(25)

3.2. The Stability Conditions of the Steady-State Solution

$a=\overline{a}+\mathrm{\Delta }a$ and $\theta =\overline{\theta }+\mathrm{\Delta }\theta$ are substituted into Equations (22) and (23). $\overline{a}$ and $\overline{\theta }$ are the singular points of Equations (22) and (23). $\mathrm{\Delta }a$ and $\mathrm{\Delta }\theta$ are small disturbances. By eliminating $\overline{\theta }$, the following can be obtained:

$a<x_0$:

$$\left\{\begin{array}{c}\frac{\mathrm{d}\mathrm{\Delta }a}{\mathrm{d}t}=-R_{11}\mathrm{\Delta }a-R_{12}\mathrm{\Delta }\theta \\ \frac{\mathrm{d}\mathrm{\Delta }\theta }{\mathrm{d}t}=R_{13}\mathrm{\Delta }a-R_{14}\mathrm{\Delta }\theta \end{array}\right.$$

(26)

$a>x_0$:

$$\left\{\begin{array}{c}\frac{\mathrm{d}\mathrm{\Delta }a}{\mathrm{d}t}=-R_{11}\mathrm{\Delta }a-R_{21}\mathrm{\Delta }\theta \\ \frac{\mathrm{d}\mathrm{\Delta }\theta }{\mathrm{d}t}=R_{22}\mathrm{\Delta }a-R_{14}\mathrm{\Delta }\theta \end{array}\right.$$

(27)

The detail of $R_{i,j}\left(i=1,2;j=1,2,3,4\right)$ is written in Appendix A.

Based on the Lyapunov theory, calculate the eigenvalues of the determinant, and the stability conditions of the steady-state solutions are:

$a<x_0$:

$$\left(R_{11}+R_{14}>0\right)\wedge \left(R_{12}{R}_{13}+R_{11}{R}_{14}>0\right)$$

(28)

$a>x_0$:

$$\left(R_{11}+R_{14}>0\right)\wedge \left(R_{21}{R}_{22}+R_{11}{R}_{14}>0\right)$$

(29)

3.3. The Transition Set of the Amplitude–Frequency Response Function

The unfolding functions are constructed by Equations (24) and (25).

$a<x_0$:

$$\begin{array}{l}G\left(a,\mathrm{\Omega },A,\alpha \right)={\left(2a\zeta \mathrm{\Omega }+aH{\mathrm{\Omega }}^p\sin \left(\frac{p\pi }{2}\right)\right)}^2+\left(\frac{a}{4}\left(3a^2\left(B_2{K}_0+K_3\right)\right.\right.\\ {\left.\left.+4\left(B_4{K}_0+K_1-{\mathrm{\Omega }}^2\right)+4H{\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right)\right)\right)}^2-F^2\end{array}$$

(30)

$a>x_0$:

$$\begin{array}{l}G\left(a,\mathrm{\Omega },A,\alpha \right)={\left(2a\zeta \mathrm{\Omega }+aH{\mathrm{\Omega }}^p\sin \left(\frac{p\pi }{2}\right)\right)}^2+\left(\frac{K_0}{\pi }\left(-2B_1+x_0\right)\sqrt{1-{\left(\frac{x_0}{a}\right)}^2}\right.\\ +\frac{K_0}{\pi }\left(2B_5+x_0\right)\sqrt{1-{\left(\frac{x_0}{a}\right)}^2}-\frac{K_0{x}_0}{2\pi }\sqrt{1-{\left(\frac{x_0}{a}\right)}^2}\left(3a^2{B}_2+4B_4+2B_2{x}_0^2\right)\\ +\frac{2}{\pi }a{K}_0\mathrm{acos}\left(\frac{x_0}{a}\right)+\frac{a}{2\pi }{K}_0\mathrm{asin}\left(\frac{x_0}{a}\right)\left(3a^2{B}_2+4B_4\right)+a{K}_1+\frac{3}{4}{a}^3{K}_3\\ {\left.-a{\mathrm{\Omega }}^2+aH{\mathrm{\Omega }}^p\cos \left(\frac{p\pi }{2}\right)-F\cos \left(\theta \right)\right)}^2-F^2\end{array}$$

(31)

The excitation frequency $\mathrm{\Omega }$ and the excitation amplitude $A$ are defined as bifurcation parameters. Based on the singularity theory, the bifurcation plane of the amplitude–frequency response is derived. It is too complex to obtain the explicit formulation of the transition set. Therefore, the transition set obtained is

$$\sum =Bif\cup Hys$$

(32)

where $Bif$ is the bifurcation set and $Hys$ is the hysteresis set.

$$Bif=\left\{\left(f,\alpha \right)\in {ℝ}^2\left|\exists \left(a,\mathrm{\Omega }\right),\exists G=G_a=G_{\mathrm{\Omega }}=0\right.\right\}$$

(33)

$$Hys=\left\{\left(f,\alpha \right)\in {ℝ}^2\left|\exists \left(a,\mathrm{\Omega }\right),\exists G=G_a=G_{aa}=0\right.\right\}$$

(34)

The specific formulations of Equations (33) and (34) are written in Appendix B.

4. The Nonlinear Characteristics Analysis of the SD Oscillator with Nonlinear Fractional Damping

4.1. The Nonlinear Characteristics in the Resonant Region

First, the correctness of Equations (24) and (25) should be verified. The ODE45 method is used to calculate the numerical solutions to obtain the amplitude–frequency response of the system. The calculation step $t_{step}=0.0025$. Only forced vibration phenomena are concerned, so that the initial displacement and the initial velocity are zero to exclude the influence of transient resonance. Thus, the results of the numerical solution, Equations (24) and (25) are shown in Figure 5.

Figure 5 displays the results of the amplitude–frequency response of Equations (24) and (25), and the numerical solution. The pink dashed line represents the vibration amplitude $a=x_0$. The blue and green solid lines represent the stable solutions calculated by Equations (24) and (25), respectively, while the unstable solution is drawn as a green dashed line. Circles are the numerical solution. The numerical solution is in good agreement with the theoretical solution. The correctness of Equations (24) and (25) is verified.

4.1.1. The Influence of the Smooth Parameter in the Transition Set of the System

Based on Equation (32), the influence of the smooth parameter and excitation parameters in the amplitude–frequency response is studied. The parameters selected are $\zeta =0.1443$, $K_0=20$, $K_1=1$, $K_3=0.0128$, $H=1$, and $\alpha \in [0.89,0.99]$. The transition set, the amplitude–frequency response curves and the distribution of the attractors in the frequency island are depicted in Figure 6.

Figure 6a,b show the results of the bifurcation set and the hysteresis set of Equations (33) and (34), respectively. The explicit function of the bifurcation set and the hysteresis set cannot be directly obtained because the expressions of Equations (33) and (34) are too complex. Thus, in Figure 6a, three surfaces, $G=0$ (the red surface), $G_a=0$ (the green surface) and $G_{\mathrm{\Omega }}=0$ (the blue surface), are drawn in the three-dimension space $(\mathrm{\Omega },A,a)$, and their intersection $Bi{f}_1=\left(1.6606,0.051901,0.38956\right)$ is obtained as the numerical solution of the bifurcation set. Using the same numerical method, in Figure 6b, the result of the hysteresis set $Hy{s}_1=\left(1.0149,0.06301,0.34944\right)$ can be obtained, and the red surface is $G=0$, the green surface is $G_a=0$, and the blue surface is $G_{aa}=0$. After obtaining the transition set, all of the bifurcation conditions of the amplitude–frequency response of the SD oscillator can be calculated.

As shown in Figure 6c, a hysteresis point is generated, as the smooth parameter $\alpha$ equals 0.89. The hysteresis phenomenon leads to the appearance of pitchfork bifurcation. Comparing the green curve and the red curve in Figure 6c, the number of the results of $a$ increases from 1 to 3 in the interval $0.5<\mathrm{\Omega }<1.5$, so that a new multi-solution coexistence region is generated. When $\alpha =$ 0.8594, a bifurcation point is generated. The bifurcation phenomenon leads to the appearance of transcritical bifurcation. Subsequently, a frequency island, which is depicted by closed black curves, occurs in the resonant region. There are three coexistent attractors in the entire region of the frequency island. The three attractors are depicted in Figure 6d by orange, purple and cyan circles, respectively. The attraction domains are drawn in the same colors as the attractors. The orange attractor is a stable attractor with the largest amplitude and the widest attraction domain, corresponding to the top black branch of the frequency island; the purple attractor is a stable attractor with the smallest amplitude and a narrow attraction domain corresponding to the bottom black solid line; and the other unstable attractor corresponds to the black dashed line in Figure 6c.

The above novel nonlinear phenomena in the small amplitude region are disclosed due to the fact that the nonlinear piecewise function $P_n$ maintains the nonlinear property of the system in the interval $-x_0<x<x_0$. Because of the frequency island, the vibration amplitude of the SD oscillator can be adjusted by imposing a disturbance to change the displacement or the velocity in order to make them exceed the boundary of the attraction domain.

Next, the frequency island is further analyzed in terms of the distributions of the attractors and the time history of vibration. Thus, the distribution of attractors in the frequency island, the time history of vibration of stable attractors and those without the frequency island are drawn in Figure 7 and Figure 8. To highlight the differences, the parameters in Figure 7 and Figure 8 are selected as $\mathrm{\Omega }=2$, $\alpha =0.85$ and $\mathrm{\Omega }=2$, $\alpha =0.9$, respectively.

Figure 7a shows that two stable attractors (the orange circle and the purple circle) with different amplitudes coexist with an unstable attractor (the cyan circle) in the frequency island. Figure 7b shows the time history of vibration. The stable motion in the black line corresponds to the orange attractor in Figure 7a, and the stable motion in the red line corresponds to the purple attractor in Figure 7a. In Figure 7a, the black arrows and black dots are schematic representations of the attractor-switching phenomenon, supposing a disturbance is imposed in the system, which is in the stable state of the orange attractor. The displacement and the velocity of the system are changed to the first black dot, which marks that the system has entered the purple attraction domain. After a period of time, the system goes through several unstable states (the other black dots) before finally returning to a stable state (the purple circle). The vibration characteristic also changes from the black line to the red line in Figure 7b, and the vibration amplitude significantly decreases. However, without the frequency island, this phenomenon cannot be generated because there is only one attractor, as shown in Figure 8.

In Figure 8a, there is only one stable attractor, and the attraction domain is full of the entire state space. Thus, no matter what the initial conditions are, for example $\left(1,0\right)$ or $\left(0,0\right)$, the stable state of the system can always be represented by the red line in Figure 8b. For most nonlinear systems, a narrow hysteresis region where the attractor-switching phenomenon can be generated, such as $2.5<\mathrm{\Omega }<4.1$ in Figure 6c, is found in the resonant region. It is rare that the frequency island results in the entire resonant region generating the attractor-switching phenomenon.

4.1.2. The Influence of the Fractional Damping Parameters in the Transition Set of the System

From the above investigations, when the smooth parameter $\alpha$ decreases and nonlinearity increases, a transition process is disclosed: a hysteresis point followed by a bifurcation point and a frequency island. Further, with the frequency island, the fractional SD oscillator has a novel attractor-switching phenomenon among the entire resonant region. Thus, it is worth investigating the influence of the fractional damping parameters on the transition set. The results of the transition sets with different fractional damping parameters are shown in Figure 9.

Figure 9a,b show the influence of the fractional damping parameters on the hysteresis set of the system. As the fractional parameters increase, the hysteresis point requires the generation of a larger excitation amplitude and higher excitation frequency. In other words, the hysteresis point moves to the high-frequency and large-amplitude region. Figure 9b,c exhibit the influence of the fractional damping parameters in the bifurcation set. As the fractional damping order $p$ increases, the bifurcation point moves to the high-frequency region, but the excitation amplitude has little change. As the fractional damping coefficient $H$ increases, the bifurcation point moves to the high-frequency and small-amplitude region. In other words, the system needs a small excitation amplitude to generate a bifurcation point.

After analyzing the influence of the fractional damping on the hysteresis set and the bifurcation set, the influence of these parameters in the frequency island is investigated. The smooth parameter is selected as $\alpha =0.7$. The amplitude–frequency response curves with different fractional damping parameters are drawn in Figure 10.

Figure 10 shows the amplitude–frequency response curves with differential fractional damping parameters. With the increase in the fractional parameters, both the frequency interval and amplitude of the stable solution in the frequency island decrease. When $p=0.9$ in Figure 10a and $H=3.5$ in Figure 10b, the frequency island disappears, while the stable solution with a small amplitude remains. The amplitude of the existing attractor decreases as the fractional order increases, while the amplitude increases as the fractional coefficient increases. The results correlate well with Figure 4, because with the fractional order equals 1 and $\mathrm{\Omega }>{\mathrm{\Omega }}_2=e^{-\pi \cot (p\pi /2)/2}$, the fractional damping has the equivalent of linear damping and is the strongest damping effect, causing the amplitude to decrease. As the fractional coefficient increases and $\mathrm{\Omega }>{\mathrm{\Omega }}_1=e^{\pi \tan (p\pi /2)/2}$, the fractional damping has a stiffness property and a damping property, but the stiffness property is dominated, causing the amplitude to slightly increase.

To further investigate the reason why the frequency island disappears, the distribution of attractors in the frequency island is depicted in Figure 11 and Figure 12.

In Figure 11 and Figure 12, the stable attractor with a large amplitude is drawn in green, the stable attractor with a small amplitude is drawn in red and the unstable attractor is drawn in dark blue. These figures show the influence of the fractional parameters on the attractors in the frequency island. Figure 11 shows that as the fractional order increases, the green attractor moves closer to the blue attractor; meanwhile, the red attraction domain expands. Then, the green attractor collides with the blue attractor and disappears, while the red attraction domain fills the entire state space since $p=0.7$. This is the reason why the frequency island disappears as the fractional order increases.

Figure 12 shows that as the fractional coefficient increases, the green attractor collides with the blue attractor and becomes a new chaotic attractor (the chaotic attraction domain and the chaotic attractor are multicolored in Figure 12c). Subsequently, the chaotic attractor disappears; meanwhile, the red attraction domain expands until it fills the entire state space as shown in Figure 12d. This is the reason why the frequency island disappears as the fractional coefficient increases.

4.2. Analysis of the Snap-Through Phenomenon in the Non-Resonant Region

In the non-resonant region, this article mainly focuses on the snap-through phenomenon that is commonly discovered in multi-well dynamics systems. A multi-well system which generates snap-through phenomena has been analyzed in terms of the nonlinear energy and the competing resonance. When the energy is sufficient to exceed the energy threshold of the potential well, two types of transition phenomena are generated: one is inner-well motion → inter-well motion → inner-well motion (it is a different potential well from the initial well), and the other is inner-well motion → inter-well motion [8,12].

In this section, the influence of the fractional parameters on the snap-through phenomenon in the non-resonant region is analyzed in terms of the global dynamics. The distribution of attractors is obtained when $\mathrm{\Omega }=$ 1 (the frequency is lower than the natural frequency) and 8.5 (the frequency is higher than the natural frequency). The initial displacement $x_s$ and the initial velocity ${\dot{x}}_s$ are selected as $\left(x_s,{\dot{x}}_s\right)=\left(-0.7,0\right)$ to keep the steady-state motion of the system corresponding to the same initial attractor. The random impulse disturbance, which has an amplitude of $x_{im}\in \left[-1,4\right]$, is applied 100 times, and the experiment is repeated 100 times. The probabilities that the system generates snap-through due to the random disturbance and the distribution of attractors in the non-resonant region are shown in Figure 13, Figure 14, Figure 15 and Figure 16.

Figure 13a,b show that the distribution of the attraction domain is helicoid. Two stable attractors (the cyan circle and the pink circle) and one unstable attractor (the purple circle) coexist. By increasing the fractional order, the continuity of the attraction domain increases, and the symmetry enhances. Thus, as shown in Figure 13c,d, the probability of the snap-through decreases from 62% to 40%, which means that the asymptotic stability of the steady-state solution increases.

Figure 14a,b show that by increasing the fractional coefficient, the continuity of the attraction domain of the stable solutions increases, and the symmetry enhances. The probability of the snap-through decreases from 54% to 28% in Figure 14c,d. The asymptotic stability of the steady-state solution increases.

When the excitation frequency is smaller than the natural frequency, the increase in the fractional parameters leads to the continuity and symmetry of the attraction domain of the stable solution increasing. Meanwhile, the asymptotic stability of the steady-state solution increases.

Figure 15a shows that when the excitation frequency is larger than the natural frequency, one stable period-three attractor (the green circles) and three period-one attractors coexist. The attraction domain of the period-one attractors is discontinuous and narrow, which means that the asymptotic stability of the period-one motion is weak. The probability of the snap-through is 77%, which is shown in Figure 15c. By increasing the fractional order, the continuity of the period-one attraction domain increases, and the period-three attractor disappears (the existence conditions of the period-three are broken), as shown in Figure 15b. The probability of the snap-through decreases from 77% to 47%. Note that the two period-one attraction domains are still asymmetric; therefore, the probability of the snap-through in Figure 15d is higher than that in Figure 13d.

Figure 16a,c,e show that by increasing the fractional coefficient, three period-one attractors merge and evolve to a period-one attractor, and the period-one attraction domain expands. The asymptotic stability of the period-one motion increases. Meanwhile, the period-three attractor disappears because the existence condition of the period-three motion is broken. As the fractional coefficient equals 4, the period-three attractor disappears, while the period-one attraction domain fills the entire state space, as shown in Figure 16e. The system evolves from a tristable state to a bistable state and finally to a monostable state. Figure 16b,d,f illustrate that as the continuity of the attraction domain increases, the probability of the snap-through decreases from 73% to 0%. No snap-through occurs when the system is monostable, as shown in Figure 16f.

When the frequency is greater than the natural frequency, caused by an increase in the fractional coefficient, the asymptotic stability of the steady-state solution increases, and the number of the coexisting attractors changes. Thus, the probability of the snap-through decreases. This influence is different from the increase in the fractional order.

5. Conclusions

A new nonlinear piecewise function is proposed to establish the differential equation of the SD oscillator with fractional damping. The Caputo fractional model was used to represent the fractional damping, and the equivalent stiffness coefficient and the equivalent damping coefficient were calculated. The nonlinear dynamics characteristics and snap-through phenomena were studied. The influence of fractional damping on nonlinear dynamic characteristics was analyzed, and some novel and interesting nonlinear phenomena were disclosed. The conclusions from the investigation are as follows:

• The nonlinear restoring force is accurately represented by the piecewise nonlinear function. The nonlinear characteristics of the restoring force in the interval $-x_0<x<x_0$ are retained, so that some novel nonlinear phenomena are found.
• The orthogonal function is used to calculate the equivalent fractional coefficients. The equivalent stiffness coefficient and the equivalent damping coefficient are variable with respect to the $p\mathrm{th}$ and the $(p-1)\mathrm{th}$ power of the excitation frequency. In the high-frequency region, the stiffness characteristic of the fractional model is dominant. The stiffness characteristic increases first and then decreases as the fractional order increases. In the low-frequency region, the damping characteristic of the fractional model is dominant. The damping characteristic increases first and then decreases as the fractional order increases. The fractional model affects the stiffness property and the damping property, simultaneously.
• Based on the amplitude–frequency response functions, a novel transition process is found. With the decrease in the smooth parameter, the nonlinearity of the system increases. A hysteresis point appears first, followed by a bifurcation point and a frequency island. There are three attractors (two stable attractors and one unstable attractor) in the frequency island. In addition, the variation in the number and the stable state of attractors means that pitchfork bifurcation and transcritical bifurcation are found. It is rare that the vibration amplitude of the system can be changed in the entire resonant region because of the frequency island.
• As the fractional parameters increase, the hysteresis point moves to the high-frequency and large-amplitude region, and the bifurcation point moves to the high-frequency and small-amplitude region. The frequency interval of the frequency island shortens. Finally, the frequency island disappears because the stable attractor and the unstable attractor collide.
• In the non-resonant region, the increase in the fractional parameters leads to the probability of the snap-through decreasing and the asymptotic stability of the steady-state solution increasing. When the excitation frequency is smaller than the natural frequency, the symmetry of the attraction domain enhances and the continuity increases. When the excitation frequency is larger than the natural frequency, the number of stable states of the system decreases. When the system is in a monostable state, no snap-through occurs.

Based on the above investigations, some novel nonlinear phenomena, such as a frequency island, of the SD oscillator with nonlinear fractional damping are revealed. This investigation may provide a theoretical basis for applying SD oscillators in energy harvesting, vibration isolation, and many other fields. In the future, some vibration isolation experiments and theoretical investigations for fractional SD oscillators utilizing the novel frequency island phenomenon may be useful for ultra-low-frequency vibration isolators.