Sage Journals: Discover world-class research

Abstract

The precision instruments and equipment are often utilized in low-frequency and micro-amplitude vibration systems, in which many vibration isolators of rubber materials are widely used. Ignoring the low-frequency amplitude will result in errors in the fatigue life design of the vibration isolators and predicting the dynamic response of each frequency band accurately becomes necessary. However, integer-order models cannot describe the frequency dependence of rubber materials, while the fractional-order models can describe it instead. On the other hand, the elastic restoring force is strongly nonlinear under large deformation, and the vibration of the nonlinear system contains multiple harmonic components. In order to solve those issues, the fractional nonlinear Nishimura model is used to characterize the constitutive relation of vibration isolators such as air springs, which are made of carbon black filled natural rubber. The high-order harmonic balance method is used to obtain the steady-state response of the vibration system, while the fourth-order Runge–Kutta method is applied to simulate the dynamic response of the system in the low-frequency region, and the Lyapunov exponent is used to determine the stability of the system. Furthermore, the influence of parameters on the amplitude–frequency characteristics of different frequency bands is also studied, and a method to solve the optimal damping coefficient is proposed based on the primary resonance amplitude–frequency curves. The results show that there is a diversity of periodic motions in the process of adjacent super-harmonic resonance transition. Numerical simulations also demonstrate that multi-periodic motions coexist in the system. The motion transition law between the polymorphic coexistence region and its adjacent regions is summarized.

Keywords

nonlinear Nishimura model fractional-order derivative super-harmonic resonance chaos transition law

Introduction

With the development of society, viscoelastic materials have been developed rapidly in the field of vibration control and have shown superior benefits in reducing noise and vibration. It is widely applied in aviation, machinery, civil engineering, as well as many other fields, such as metal rubber in spacecraft vibration reduction, and air springs in high-speed train bogies.

In order to accurately characterize the relaxation and creep of viscoelastic materials, the Kelvin model, Maxwell model, Zener model, and Nishimura model are applied to describe viscoelastic system.^1–5 Zheng et al.⁶ used the Kelvin model to describe the mechanical properties of air springs and analyzed the oscillation attenuation characteristics of the transmission force of air springs. Zhu et al.⁷ proposed a linear Zener model that can describe the dynamic stiffness of air springs with an auxiliary chamber based on engineering thermodynamics and the aerodynamic theory. Wu et al.⁸ investigated the accuracy of the dynamic stiffness of air springs described by the Zener model based on the thermodynamic point of view and found that air springs are dependent on nonlinear parameters and frequencies. Larsson et al.⁹ studied the micro-mechanically driven, quasi-static to dynamic failure response of fiber-reinforced unidirectional composite viscoelastic materials under finite deformation conditions based on the Zener model. Besides, the Nishimura model is also adopted to describe the mechanical properties of the air spring.¹⁰ Ryaboy¹¹ discussed a series of issues that need to be considered in the establishment of dynamic models from the perspective of air spring design and application. The rubber materials of the rail pads and air springs, mainly, are carbon black filled natural rubber.^12,13 The viscoelastic properties of rubber materials can be changed by changing the molecular structure or adjusting the filler proportion. The external force response of rubber is divided into the elastic and viscous parts, and the strain will fall behind the stress. The Nishimura model is more popular to be applied to describe the dynamic properties of air springs in practical situations.¹⁰ Compared with the Kelvin model and the Zener model, the Nishimura model is more popular to be applied to describe the mechanical properties of air springs in practical situations.

Fractional calculus has not been investigated and developed for a long time due to the lack of practical application background. In recent years, with the rapid development of computer technology, the definition, characteristics, and calculation of fractional calculus have been more and more applied in the field of engineering.^14,15 There are several definitions of fractional-order derivatives, such as Grünwald-Letnikov, Riemann-Liouville, Caputo definitions,¹⁶ He’s fractional definition,¹⁷ He’s fractal derivative,¹⁸ and so on.¹⁹ The local fractional derivative was first proposed by Yang, and it is an excellent tool to describe a lot of natural phenomena in fractal space.²⁰ Recently, the two-scale fractal theory appeared as a powerful mathematical tool to porous media or unsmooth boundary problems,²¹ and one is for the continuum model on a large scale, the other is figuring out the porous structure of the studied object on a small scale.²² Moreover, the nonlinear vibration in a porous medium can be modeled by the fractal vibration theory, which was first proposed by He et al.^23,24 Although traditional integer-order models can be applied to describe the mechanical properties of vibration isolation materials, the description of the frequency dependence phenomenon in vibration isolation systems is lack of investigations.⁸ From this point of view, the fractional-order model is more accurate in practical engineering situations.^12,25,26 Fractional-order models can be applied to not only describe the amplitude dependence and frequency dependence of rubber but also the stress relaxation and creep properties of viscoelastic materials.^27–29 However, the vibration isolation system can be approximated as a linear elastic system; only when the deformation is small, it behaves as a nonlinear system beyond this range.^30,31

Although the nonlinear characterization is closer to practical engineering, the dynamic response is more complex, even bifurcation and chaos occur. Among all kinds of approximate analytical methods, the harmonic balance method is the clearest in concept. For vibro-impact problems, the harmonic balance method is widely used in systems with nonlinear restoring forces, such as quasi-zero stiffness systems,^32,33 origami metamaterial systems,³⁴ and viscoelastic systems.³⁵ Besides, the harmonic balance method is also applied to the computation of periodic orbits for piecewise linear oscillators.³⁶ In addition to the harmonic balance method, many other effective methods for solving nonlinear systems exist. Wang et al.²⁰ used local fractional calculus and applied the variational iteration method to obtain the KdV-Burgers-Kuramoto equation approximate analytical solution. He et al.²¹ used the two-scale transform to convert the nonlinear Zhiber-Shabat oscillator with the fractal derivatives to the traditional model. The homotopy perturbation method has been demonstrated under a suitable transformation of the system containing several exponential nonlinear terms to the famous Helmholtz-Duffing oscillator. Chen et al.³⁷ established a nonlinear dynamic model of the suspension of the maglev system and solved the chaotic threshold by using the Melnikov method. Xu³⁸ proposed a nonlinear Kelvin model of a rolling lobe air spring and analyzed the relationship between geometric parameters and stiffness characteristics through experiments. Li et al.³⁹ developed a viscoelastic model of dielectric elastomer and studied the influence of the viscoelastic strain enhancement effect on the resonance frequency and oscillation amplitude. Sharma et al.⁴⁰ applied the Zener model to investigate the impact of visco-elasticity on the electromechanical response of complex dielectric elastomeric actuators involving inhomogeneous deformations. Qi et al.⁴¹ utilized the Lyapunov exponent to analyze the stability of fractional-order nonlinear time-varying systems. Shen et al.^16,42,43 studied the approximate analytical solutions of fractional Duffing oscillator and fractional van der Pol oscillator through the average method and proposed the concepts of the equivalent linear stiffness and the equivalent linear damping. Zhang et al.³² predicted the dynamic response of a geometrically nonlinear system with negative stiffness by using the higher harmonic balance method and analyzed its vibration isolation performance. Zhou et al.⁴⁴ studied the dynamic response and stability of a two-degree of freedom vehicle system under nonlinear external force. Jena et al.⁴⁵ analyzed the time-fractional fuzzy vibration equation of large membranes by using the double parametric-based Residual power series method. Matouk⁴⁶ studied the local stability of the system’s three equilibria when the fractional parameter belongs to (0,2]. According to Hopf bifurcation theory in fractional-order systems, approximations to the periodic solutions around the system’s three equilibria are explored.

In practical engineering, viscoelastic materials are neither ideal solid nor ideal fluid, but between them. Polyborosiloxanes were invented initially in the search for substitutes of natural rubbers, possess reversible physical cross-linking, and the viscoelasticity of rubber materials can be adjusted by changing the molecular structure of cross-linking.⁴⁷ However, the rubber materials exhibit frequency-dependent and nonlinear characteristics.⁷ Fortunately, fractional-order calculus can represent the frequency-dependent properties of viscoelastic materials, and the superelasticity of a material can be described by the Duffing equation.⁴⁸ Due to nonlinear factors, the dynamic response is more complex, especially in the low-frequency region, and even super-harmonic resonance occurs. The precision instruments and equipment are often utilized in low frequencies and micro-amplitude vibration environments. However, the vibration isolators of rubber materials are widely used in micro-vibration isolation. As we know, the fatigue failure of components is related to the amplitude of different frequencies. Ignoring low-frequency amplitude will result in errors in the fatigue life design of the vibration isolators, and predicting the dynamic response of each frequency band accurately becomes necessary.^49,50 The nonlinear system also has a diversity of periodic motions, the motion transition process is more complicated, and even polymorphism coexists, which affects the accurate prediction of the dynamic response in the low-frequency region.⁵¹ The investigations on the dynamic characteristics of integer-order nonlinear systems have been relatively mature. In contrast, the diversity of periodic motions and the transition laws of fractional-order nonlinear systems need to be revealed.^52,53

To reveal the dynamic response of the fractional nonlinear vibration isolation system, this paper focuses on the investigations of super-harmonic resonance and the diversity of periodic motions of the system. The high-order harmonic balance method is used to obtain the steady-state response of the vibration system, while the fourth-order Runge–Kutta method is applied to simulate the dynamic response of the system in the low-frequency region. Furthermore, the influence of parameters on the amplitude–frequency characteristics of different frequency bands is also studied. On the other hand, two Poincaré maps are established to describe the diversity of periodic motions and the transition law, and the Lyapunov exponent is used to determine the stability of the system. Finally, the relationship between super-harmonic resonance and the diversity of periodic motions of the system is analyzed. The system parameters that affect amplitude-frequency characteristics, bifurcation, and chaos are further revealed.

Analysis of fractional nonlinear Nishimura model

The air springs are made of carbon black filled natural rubber,¹³ and the mechanical properties of rubber materials are between ideal elastic bodies and ideal viscous bodies. However, integer-order models cannot describe the frequency dependence of rubber materials, while the memory recovery force is needed to be described by fractional differential operators.^15–17 The Nishimura model is more popular to be applied to describe the dynamic properties of air springs in practical situations, and model parameters are related to the upper and lower chamber stiffness of the air spring and the gas viscosity through the orifice plate.¹¹ The fractional nonlinear Nishimura model is applied to describe the viscoelasticity and nonlinear features of rubber and air springs as shown in Figure 1.

Figure 1.

The model of fractional nonlinear vibration isolation system.

Consider the following motion state of the vibration isolation system under periodic excitation. The equation of the dynamic model is

{\begin{array}{c} M {\ddot{X}}_{1} + K_{1} (X_{1} - X_{2}) + K_{3} X_{1} + K_{4} X_{1}^{3} + K D_{T}^{p} [X_{1} (T)] = F \cos (Ω T) \\ K_{1} (X_{1} - X_{2}) - K_{2} X_{2} - C {\dot{X}}_{2} = 0 \end{array}

(1)

where M is the mass of the system, F is the external excitation amplitude, X₁ is the displacement of the mass block, X₂ is the displacement of the node,

Ω

is the external excitation frequency, K₁ and K₂ are the stiffness coefficients of the linear elastic restoring force, C is the viscous damping coefficient, F_K is the nonlinear elastic restoring force, D is the fractional-order derivative, K is the fractional coefficient, and P is the fractional order.

The control parameter of external excitation amplitude F is introduced into Fs, and the following variables are replaced.

x_{i} = \frac{K_{1}}{F_{s}} X_{i}, t = T \sqrt{\frac{K_{1}}{M}}, ω = Ω \sqrt{\frac{M}{K_{1}}}, ξ = \frac{C}{2 \sqrt{K_{1} M}}, μ_{k 2} = \frac{K_{2}}{K_{1}}, μ_{k 3} = \frac{K_{3}}{K_{1}}, f = \frac{F}{F_{s}}, ε = \frac{K_{4}}{K_{1}^{3}} F_{s}^{2}

The dimensionless differential equation of the system can be rewritten as

{\begin{array}{c} {\ddot{x}}_{1} + (1 + μ_{k 3}) x_{1} - x_{2} + ε x_{1}^{3} + λ D_{t}^{p} [x_{1} (t)] = f \cos (ω t) \\ 2 ξ {\dot{x}}_{2} + x_{2} - x_{1} + μ_{k 2} x_{2} = 0 \end{array}

(2)

The fractional-order derivative defined by Caputo is conducive to the solution of the actual physical system and has more application value in practical engineering.^16,42,43 In order to represent the frequency-dependent properties of viscoelastic materials, here Caputo’s definition is adopted

D_{t}^{p} [x (t)] = \frac{1}{Γ (1 - p)} \int_{0}^{t} \frac{x^{'} (u)}{{(t - u)}^{p}} d u

(3)

Γ (y + 1) = y Γ (y)

(4)

where

Γ (n)

is the Gamma function satisfying equation (4). Supposing the steady-state response of the mass as

x_{1} = \sum_{n = 1}^{N} A_{2 n - 1} \cos [(2 n - 1) ω t] + \sum_{n = 1}^{N} B_{2 n - 1} \sin [(2 n - 1) ω t]

(5)

Substituting equation (5) into equation (2) and introducing the formula in^42,43

{\begin{array}{c} \int_{0}^{t} \frac{\sin (ω s)}{s^{p}} d s = ω^{p - 1} Γ (1 - p) \cos (\frac{p π}{2}) \\ \int_{0}^{t} \frac{\cos (ω s)}{s^{p}} d s = ω^{p - 1} Γ (1 - p) \sin (\frac{p π}{2}) \end{array}

(6)

Taking the first-order approximation simplification of the fractional term and equation (3) can be obtained as

D_{t}^{p} [x_{1} (t)] \approx ω^{p - 1} \sin (\frac{p π}{2}) {\dot{x}}_{1} (t) + ω^{p} \cos (\frac{p π}{2}) x_{1} (t)

(7)

The fractional term is obtained by equation (7). Viscoelastic materials are neither ideal solid nor ideal fluid, but between them. According to the Caputo definition, the fractional-order derivative can be reduced to the form of integer-order trigonometric functions, which can describe not only the viscoelasticity of materials, but also the frequency correlation of materials. Shen et al.⁴³ verified the rationality by using the first-order trigonometric function to represent the fractional-order derivative in practical engineering applications. This paper mainly studies the case where the amplitude of the higher harmonics of the vibration isolation system is much smaller than the amplitude of the fundamental wave, and the fractional-order derivative simplification ignores the high-order terms is reasonable.

Substituting equation (7) into equation (2), the dimensionless differential equation of the system can be rewritten as

{\ddot{x}}_{1} + (1 + μ_{k 3}) x_{1} - x_{2} + ε x_{1}^{3} + α {\dot{x}}_{1} + β x_{1} = f \cos (ω t)

(8)

where

α = λ ω^{p - 1} \sin (\frac{p π}{2})

and

β = λ ω^{p} \cos (\frac{p π}{2})

The equivalent linear stiffness coefficient Keq and equivalent damping coefficient Ceq of the system can be rewritten as

{\begin{array}{c} K e q = 1 + μ_{k 2} + μ_{k 3} + λ ω^{p} \cos (\frac{p π}{2}) \\ C e q = 2 ξ + λ ω^{p - 1} \sin (\frac{p π}{2}) \end{array}

(9)

It can be seen from equation (9) that the fractional coefficient $λ$ and the fractional order p influence the equivalent linear stiffness coefficient Keq and the equivalent damping coefficient Ceq. This subsection selects the following parameters

ξ = 0.1, μ_{k 2} = 1, μ_{k 3} = 1, ε = 0.2, f = 5

Figure 2 can be obtained by taking the fractional order $p = 0.5$ . As shown in Figure 2(a), the larger the fractional coefficient $λ$ is, the larger the equivalent linear stiffness coefficient Keq is, and it rises as the frequency increases. As shown in Figure 2(b), the less the fractional coefficient $λ$ is, the less the equivalent damping coefficient Ceq is. It also decreases rapidly in the low-frequency region as the frequency increases.

Figure 2.

The equivalent linear stiffness coefficient Keq and equivalent damping coefficient Ceq with different values of fractional coefficient $λ$ : (a) the equivalent linear stiffness coefficient Keq and (b) the equivalent damping coefficient Ceq.

Figure 3 can be obtained by taking the fractional coefficient $λ = 0.01$ . As shown in Figure 3(a), in the super low-frequency region, the less the fractional order p is, the larger the equivalent linear stiffness coefficient Keq is. Beyond the super low-frequency region, the less the fractional order p is, the less the effective linear stiffness coefficient Keq is, and the equivalent linear stiffness coefficient Keq rises as the frequency increases. As shown in Figure 3(b), in the super-low-frequency region, the less the fractional order p is, the more the equivalent damping coefficient Ceq is. Beyond the super-low-frequency region, the less the fractional order p is, the less the equivalent damping coefficient Ceq is. In general, the equivalent damping coefficient Ceq in the low-frequency region decreases rapidly as the frequency increases. There is no significant change in the equivalent damping coefficient Ceq beyond the low-frequency region.

Figure 3.

The equivalent linear stiffness coefficient Keq and equivalent damping coefficient Ceq with different values of fractional order p: (a) the equivalent linear stiffness coefficient Keq and (b) the equivalent damping coefficient Ceq.

It is observed that different values of fractional coefficient $λ$ and fractional order p will change the equivalent damping coefficient Ceq and the equivalent linear stiffness coefficient Keq. Even more, it will change the natural frequency of the system.

The steady-state response of the node could be supposed as

x_{2} = \sum_{n = 1}^{N} a_{2 n - 1} \cos [(2 n - 1) ω t] + \sum_{n = 1}^{N} b_{2 n - 1} \sin [(2 n - 1) ω t]

(10)

where (2n−1) indicates the order of harmonics. When n = 1, according to equations (5) and (10), the first-order harmonic responses of the mass block and the node can be obtained as

{\begin{array}{c} x_{1} = A_{1} \cos (ω t) + B_{1} \sin (ω t) = A_{11} \sin (ω t + θ_{1}) \\ x_{2} = a_{1} \cos (ω t) + b_{1} \sin (ω t) = a_{11} \sin (ω t + ϕ_{1}) \end{array}

(11)

Substituting equation (11) into equation (2) to perform first-order harmonic balance, the relationship between the mass block and the node of amplitude and phase can be obtained as

{\begin{array}{c} A_{11} = a_{11} \sqrt{1 + μ_{k 2} + {(2 ξ ω)}^{2}} \\ θ_{1} = φ_{1} + \arctan (\frac{2 ξ ω}{1 + μ_{k 2}}) \end{array}

(12)

Substituting equation (12) into equation (11), the first-order steady-state response of the mass block can be obtained as

x_{1} = a_{11} \sqrt{1 + μ_{k 2} + {(2 ξ ω)}^{2}} \sin [ω t + φ_{1} + \arctan (\frac{2 ξ ω}{1 + μ_{k 2}})]

(13)

Substituting equations (13) and (11) into equation (8) can be obtained as

{\begin{array}{c} ψ_{1} a_{11}^{3} + ψ_{2} a_{11} = f \sin (φ_{1}) \\ ψ_{3} a_{11}^{3} + ψ_{4} a_{11} = f \cos (φ_{1}) \end{array}

(14)

where equation (14) satisfies the following formula

{\begin{array}{c} ψ_{1} = (3 ξ^{2} ε μ_{k 2} + 3 ξ^{2} ε) ω^{2} + \frac{3 μ_{k 2}^{3} ε}{4} + \frac{9 μ_{k 2}^{2} ε}{4} + \frac{9 μ_{k 2} ε}{4} + \frac{3 ε}{4} \\ ψ_{2} = (- 2 α ξ - μ_{k 2} - 1) ω^{2} + β μ_{k 2} + μ_{k 2} μ_{k 3} + μ_{k 2} + μ_{k 3} + β \\ ψ_{3} = 6 ξ^{3} ε ω^{3} + (\frac{3}{2} μ_{k 2}^{2} ξ ε + 3 μ_{k 2} ξ ε + \frac{3}{2} ε ξ) ω \\ ψ_{4} = - 2 ξ ω^{3} + (α μ_{k 2} + 2 ξ β + 2 ξ μ_{k 3} + α + 2 ξ) ω \end{array}

(15)

Then the first-order amplitude and phase of the node can be solved as

f^{2} = {(ψ_{1} a_{11}^{3} + ψ_{2} a_{11})}^{2} + {(ψ_{3} a_{11}^{3} + ψ_{4} a_{11})}^{2}

(16)

φ_{1} = \arctan (\frac{ψ_{1} a_{11}^{3} + ψ_{2} a_{11}}{ψ_{3} a_{11}^{3} + ψ_{4} a_{11}})

(17)

The first-order harmonic balance method can only obtain the primary resonance of the system but cannot obtain the super-harmonic resonance. In order to further study the dynamic response in the low-frequency region, the high-order harmonic balance method is applied.⁵² The harmonic amplitude coefficients of each order of the node are solved at first, and the maximum displacement in the time history graph is selected as the amplitude. Finally, the high-order response of the mass block is obtained.⁵³ When n = 2, the third-order harmonic balance method is used to suppose the response of the node as

x_{2} = a_{1} \cos (ω t) + b_{1} \sin (ω t) + a_{3} \cos (3 ω t) + b_{3} \sin (3 ω t)

(18)

Performing the third-order harmonic balance on equation (18), the following equations can be obtained as

{\begin{array}{c} D_{13} ω^{3} + D_{12} ω^{2} + D_{11} ω + D_{10} = f \\ D_{23} ω^{3} + D_{22} ω^{2} + D_{21} ω + D_{20} = 0 \\ D_{33} ω^{3} + D_{32} ω^{2} + D_{31} ω + D_{30} = 0 \\ D_{43} ω^{3} + D_{42} ω^{2} + D_{41} ω + D_{40} = 0 \end{array}

(19)

where parameters D_1n,…, D_4n for the equation (19) have been described in the Appendix 1.

The values of the node amplitude coefficients a₁, b₁, a₃, and b₃ can be obtained with certain system parameters.⁵³ Then the amplitude coefficients of the node are substituted into equation (2), and the amplitude coefficients of the mass block can be obtained as

{\begin{array}{c} A_{1} = 2 ξ b_{1} ω + (1 + μ_{k 2}) a_{1} \\ B_{1} = - 2 ξ a_{1} ω + (1 + μ_{k 2}) b_{1} \\ A_{3} = 6 ξ b_{3} ω + (1 + μ_{k 2}) a_{3} \\ B_{3} = - 6 ξ a_{3} ω + (1 + μ_{k 2}) b_{3} \end{array}

(20)

The third-order steady response of the mass block can be obtained as

x_{1} = [2 ξ b_{1} ω + (1 + μ_{k 2}) a_{1}] \cos (ω t) + [- 2 ξ a_{1} ω + (1 + μ_{k 2}) b_{1}] \sin (ω t) + [6 ξ b_{3} ω + (1 + μ_{k 2}) a_{3}] \cos (3 ω t) + [- 6 ξ a_{3} ω + (1 + μ_{k 2}) b_{3}] \sin (3 ω t)

(21)

According to the solution process of the third-order harmonic balance method, the higher harmonic response of the system can be obtained by analogy.

In order to compare the solution results, a sample chosen system of particulars is given by the following parameters

λ = 0.1, ε = 0.1, f = 5, ξ = 0.1, p = 0.5, μ_{k 2} = 0.5, μ_{k 3} = 0.5, ω = 0.5

The fourth-order Runge–Kutta method is used to solve the system numerically, and the virtual experimental simulation of the vibration isolation system is carried out through the dynamic simulation software UM. Figure 4 illustrates that there is roughly consistency between numerical results obtained from the approximate and exact equations. The amplitude–frequency characteristics of the system can be obtained through sweeping the excitation frequency as shown in Figure 5. A small hump appears in the amplitude-frequency curve of the low-frequency range when the super-harmonic resonance is taken into consideration. It is evident that although the primary harmonic is the dominant frequency of the spectrum, it can also be observed that the super-harmonic response will occur with a small amplitude in the spectrum that affects the system behavior around the resonance peak.

Figure 4.

Time history diagram of the mass block: (a) transient response and (b) steady-state response.

Figure 5.

The amplitude–frequency curves of the mass block by different methods.

Amplitude–frequency characteristics

The fatigue failure of the parts of the vibration isolation system is related to the amplitudes of different frequencies. For a vibration isolation system, it is necessary to predict each frequency range dynamic response accurately. This subsection selects the following benchmark parameters

λ = 0.1, ε = 0.1, f = 5, ξ = 0.1, p = 0.5, μ_{k 2} = 0.5, μ_{k 3} = 0.5

The variation of parameters has great effects on the amplitude–frequency response, Figure 6 illustrates how the response curves behave.

Figure 6.

The amplitude-frequency curves of the mass block with different (a) parameter $μ_{k 3}$ , (b) parameter $μ_{k 2}$ , (c) parameter $ε$ , (d) parameter $p$ , (e) parameter $λ$ , and (f) parameter $f$ .

By comparing the amplitude–frequency curves of the system, it can be noticed from Figure 6(a) that the increase of parameter $μ_{k 3}$ will lead to lower values of the response amplitude. That is to say, the increase of parameter $μ_{k 3}$ will decrease the primary resonance peak and will flatten the super-harmonic resonance hump.

Figure 6(b) illustrates that as parameter $μ_{k 2}$ increases, the primary resonance peak decreases, and the frequency range reduces. Nevertheless, it can be seen that the peak value of the super-harmonic resonance decreases but will not flatten the super-harmonic resonance hump.

As shown in Figure 6(c), the primary resonance peak decreases from 13.01 to 10.94, while the resonance frequency increases from 3.644 to 5.196 as parameter $ε$ increases from 0.1 to 0.3. Additionally, the super-harmonic resonance peak decreases, while the resonance frequency increases.

According to Figure 6(d), the primary resonance peak decreases from 16.52 to 10.31, while the resonance frequency decreases from 4.556 to 1.791 as fractional order p increases from 0.1 to 0.9. Moreover, the frequency range of primary resonance polymorphic coexistence reduces from 2.765 to 1.165. It can be seen that the third-order super-harmonic resonance peak decreases from 4.911 to 4.667 as fractional order p increases from 0.1 to 0.5, and increases from 4.667 to 4.7 as fractional order p increases from 0.5 to 0.9.

Figure 6(e) shows that the primary resonance peak decreases from 13.01 to 9.042 and the resonance frequency decreases from 2.646 to 2.696 as fractional coefficient $λ$ increases from 0.1 to 0.3. It can be noticed that the super-harmonic resonance hump will be flattened.

As shown in Figure 6(f), the primary resonance peak increases from 7.207 to 17.89, and the resonance frequency increases from 2.187 to 4.922 as parameter f increases from 2 to 8. Moreover, the third-order super-harmonic resonance peak increases from 2.454 to 5.793, and the resonance frequency increases from 0.4107 to 0.5915.

As shown in Figure 7, the primary resonance peak decreases from 10.45 to 9.126 while resonance frequency decreases from 2.896 to 2.596 as parameter $ξ$ increases from 0.15 to 0.25. With the further increases of parameter $ξ$ from 0.25 to 0.35, the primary resonance peak increases from 9.126 to 10.4, and the resonance frequency increases from 2.596 to 2.975. To obtain the optimal parameter $ξ_{opt}$ , the amplitude–frequency curve of the primary resonance of the mass block is obtained from equations (12) and (16). The minimum point of the available curve $L$ can be obtained as

\frac{d A_{11} (ω, μ_{k 2}, μ_{k 3}, ξ, ε, α, β, f)}{d ω} = 0

(22)

Figure 7.

Effect of parameter $ξ$ on the amplitude–frequency curves of the mass block.

Then the optimal parameter $ξ_{opt}$ of primary resonance can be obtained as

\frac{\partial^{2} A_{11} (ω, μ_{k 2}, μ_{k 3}, ξ, ε, α, β, f)}{\partial ω \partial ξ} = 0

(23)

The damping coefficient $ξ$ is selected from 0.1 to 0.4, and calculate the optimal parameter $ξ_{opt}$ for the primary resonance of the system, as shown in Figure 8.

Figure 8.

Schematic diagram of optimal damping coefficient of principal resonance $ξ_{opt}$ . The diversity of periodic motions of super-harmonic response.

The diversity of periodic motions of super-harmonic response

The super-harmonic response is rigorous in the mathematical sense, and it has a great significance for engineering applications. In order to investigate the diversity and local bifurcation of the periodic motions of the system in the low-frequency region, the Poincaré mapping of the system is established. Select the mapping section $σ_{τ} = {(x, \dot{x}) \in R^{2} \times S, θ \mod 2 π = 0}$ to calculate the number of the periodic motions of the system, and $σ_{x_{1}} = {(x, \dot{x}) \in R^{2} \times S, {\dot{x}}_{1} = 0, {\ddot{x}}_{1} < 0}$ to calculate the number of harmonic components of the system. The symbol T-N is defined to describe the motion state, where T is the number of the periodic motions of the system and N is the number of simple harmonic vibrations in the periodic motions. In what follows, a sample system of particulars is given by the following parameters

λ = 0.1, ε = 0.1, f = 5, ξ = 0.1, p = 0.5, μ_{k 2} = 0.1, μ_{k 3} = 0.1

Taking n = 5 and using the high-order harmonic balance method to obtain the super-harmonic resonance of the mass block. Combined with Poincaré mapping, the motion state of the mass block can be obtained as shown in Figure 9.

Figure 9.

Dynamic characteristics of super-harmonic resonance.

It is observed that there is a diversity of periodic motions of super-harmonic resonance. SNB is used for the saddle junction bifurcation. In addition to the polymorphic coexistence and SNB in the primary resonance, the super-harmonic resonance in the low-frequency region also has polymorphic coexistence. The jumping direction of the amplitude-frequency curve is consistent with the SNB mutation direction in the bifurcation diagram. It is observed that decreasing the frequency, the periodic motions conform to the following rules: 1–1→1–3→SNB→1–3→1–5→1–7. The higher the order of the super-harmonic resonance is, the more the simple harmonic vibrations in the periodic motions are. Moreover, the order of the simple harmonic vibration in the periodic motions is approximately the same as the order of the super-harmonic resonance.

As shown in Figure 10, the dynamic response with different frequency bands is in the low-frequency region.

Figure 10.

The diversity of periodic motions of super-harmonic resonance: (a) 1–1 periodic motion at $ω = 0.6$ , (b) unstable 1–3 periodic motion at $ω = 0.5$ , (c) 1–3 periodic motion of upper curve at $ω = 0.5$ , (d) 1–3 periodic motion of under curve at $ω = 0.5$ , (e) 1–5 periodic motion at $ω = 0.3$ , and (f) 1–7 periodic motion at $ω = 0.2$ .

As frequency $ω$ decreases, the system always maintains period-1 motion (T = 1), but the number N of simple harmonic oscillations in the period increases and the kinetic responses are also different. In addition, the response of the system in period-1 (T = 1) contains not only primary harmonic but also multiple super-harmonics, which validates that the super-harmonic response does exist.

Analysis of bifurcation characteristics

Since the restoring force of the Duffing system is an odd function, the motion state of the system also presents an antisymmetric motion trajectory in the phase diagram.⁵¹ Moreover, the motion transition law is more complicated. The bifurcation characteristics and stability of periodic motions of the fractional nonlinear vibration isolation system are further analyzed. On the Poincaré mapping of the fixed phase plane, the response after a certain moment $t_{0}$ a cycle in the steady-state motion is taken as the initial value $X (t) = X (t_{0} + T)$ . Perturbations are applied to the points of the phase trajectories respectively at time t: ${\tilde{x}}_{1} = x_{1} + Δ x_{1}$ , ${\tilde{\dot{x}}}_{1} = {\dot{x}}_{1} + Δ x_{1}$ , ${\tilde{x}}_{2} = x_{2} + Δ x_{2}$ . The response $X (t_{0} + T)$ at the end of the period can be obtained by numerical calculation. The Poincaré mapping relationship from the beginning to the end is $X (t + T) = A \tilde{X} (t)$ . Where A is the Jocabi matrix at the fixed point and $X = {(x_{1} {\dot{x}}_{1} x_{2})}^{T}$ . The Lyapunov exponential spectrum is applied in this study as an effective means to determine whether the dynamic system is stable or chaotic.⁵⁴ The Poincaré mapping is used to transform the continuous dynamic system into a discrete dynamic system. The calculation formula is as follows

λ_{i} = \frac{1}{K} \sum_{k = 1}^{K} \ln ‖ z_{i}^{(k)} ‖

(24)

where K is the number of iterations,

z_{i}^{(k)}

is the vector processed by Gram-Schmidt orthogonalization and norm normalization after each iteration, and

‖ z_{i}^{(k)} ‖ (i = 1,2,3)

is the norm of

z_{i}^{(k)}

In order to investigate the distribution type of the periodic motion of the system and its transition law, the symbol $P_{AS (i)} - T$ is defined as the periodic motion state of the system in the i-th position. This subsection selects the system parameters

λ = 0.1, ε = 0.6, f = 5, ξ = 0.01, p = 0.5, μ_{k 2} = 0.01, μ_{k 3} = 0.01

The forward and backward sweeping bifurcation diagrams are depicted respectively in Figure 11, which indicates that there are coexisting motions under a given parameter condition. The Lyapunov exponent is adopted to determine the existence of the stable periodic motion and the chaotic motion. When the maximum Lyapunov exponent is close to 0, the system has bifurcation. When the maximum Lyapunov exponent is greater than 0, the system is judged to be in chaotic motion.

Figure 11.

Bifurcation and the Lyapunov exponent diagram: (a) forward sweeping and (b) backward sweeping.

The above bifurcation analysis also indicates that the system may respond in various motion states, either in stable periodic motion states or in unstable chaotic ones.

As shown in Figure 12, excitation frequency changes induce the appearance of various bifurcation types. PDB is used for period-doubling bifurcation, IPDB for inverse period-doubling bifurcation, PFB for the pitchfork bifurcation, IPFB for the inverse pitchfork bifurcation, CIC for the catastrophic bifurcation, and BC for the boundary crisis. As frequency $ω$ increases from 0.45, the period-3 motion (T = 3) rises first. When frequency $ω$ increases to 0.4605, the boundary crisis (BC) occurs and the system becomes chaotic.

Figure 12.

Polymorphic coexistence bifurcation diagram.

It is observed that for the backward sweeping of frequency, the response occurs first in P_AS2-1 motion when frequency $ω$ decreases from 0.6. As frequency $ω$ decreases to 0.561, the pitchfork bifurcation (PFB) occurs, makes the original P_AS2-1 solution unstable, and generates a pair of stable period-1 solution branches. (The forward sweeping solution follows the lower branch P_AS1 and the backward sweeping solution follows the upper branch P_AS2.) With the further decrease of frequency $ω$ , the saddle junction bifurcation occurs at $ω = 0.5192$ and the displacement $x_{1}$ occurs jump phenomenon. Then the period-doubling bifurcation (PDB) occurs at $ω = 0.4883$ , where the P_AS2-2 solution begins to exist. After a series of period doubling bifurcations occur, the system becomes chaotic. Figure 13(a)–(d) show that P_AS1-1 and P_AS2-1 motions coexist, P_AS1-2 and P_AS2-2 motions coexist, P_AS1-4 and P_AS2-4 motions coexist, and P_AS1-8 and P_AS2-8 motions coexist. Figure 13(e) shows that the chaotic attractors coexist. After the catastrophic bifurcation, the coexisting chaotic attractors are combined into a whole, as shown in Figure 13(f).

Figure 13.

Phase trajectories and attractors of different response types when (a) $ω = 0.55$ , (b) $ω = 0.48$ , (c) $ω = 0.4765$ , (d) $ω = 0.4755$ , (e) $ω = 0.4705$ , and (f) $ω = 0.4695$ .

When neighboring periodic motions transfer to each other, the transfer process is irreversible. Furthermore, the polymorphic coexistence region is composed of two neighboring periodic motions. In order to gain a better insight into the effect of frequency on the dynamics of the system, the cell mapping method is used with the selected initial state domain in $Ω_{a} = {(x, \dot{x}) | - 3 < x < 3, - 5 < \dot{x} < 5}$ with $500 \times 500$ state cells.⁵⁵ The basins of the attraction for different frequencies are shown in Figure 14. As the frequency decreases, the basin of attraction for P_AS1 (shown in green) also decreases and the basin of attraction for P_AS2 (shown in yellow) becomes the dominant solution. The ratio of the yellow area to the total area is depicted in Figure 14(a)–(c).

Figure 14.

The attraction domain on the phase-fixing surface when (a) $ω = 0.55$ , (b) $ω = 0.48$ , and (c) $ω = 0.4755$ .

Since there are various types of bifurcations in the system such as SNB, PDB, PFB, and BC, chaos can even occur under these bifurcations and the process is irreversible. The specific motion transition law is shown in Figure 15.

Figure 15.

Multi-state coexistence and motion transition law of adjacent regions.

Conclusion

In this study, the fractional nonlinear Nishimura model is used to characterize the constitutive relation of the fractional nonlinear vibration isolation system. The dynamic responses are obtained by using the high-order harmonic balance method, the super-harmonic resonance, and the diversity of periodic motions of the system are investigated. The Poincaré map is used to determine the characteristics of the periodic motions of the system, and the Lyapunov exponent is used to determine the stability of the system. Combined with the cell mapping method, the transition law between adjacent periodic motions of the system is summarized. The influence of system parameters on the amplitude–frequency characteristics is analyzed. The main conclusions are as follows

(1) There is a diversity of periodic motions in the low-frequency region of the system. The dynamic response contains not only primary harmonic component but also multiple super-harmonics components. As the excitation frequency decreases, the number of periodic motion simple harmonic vibrations increases, and the order of the simple harmonic vibration in the periodic motions is approximately the same as the order of the super-harmonic resonance. The frequency corresponding to SNB is consistent with the amplitude–frequency curve jumping frequency. Moreover, the jumping direction of the amplitude–frequency curve is consistent with the mutation direction of the saddle-junction bifurcation.

(2) There are coexisting motions under a given parameter condition. In the transition process, various bifurcation types such as PFB, SNB, PDF, and BC occur successively. Even more, chaos will occur under bifurcation inductions but the process is irreversible.

(3) The integer-order parameters of the system influence the super-harmonic resonance. As the damping coefficient $ξ$ increases, the super-harmonic resonance peak value decreases, but the primary resonance peak value first decreases and then increases. Parameter $μ_{k 3}$ can flatten the super-harmonic resonance as it increases, while parameter $μ_{k 2}$ can affect the peak of super-harmonic resonance as it increases as well.

(4) The fractional-order derivatives of the system influence the super-harmonic resonance. As fractional coefficient $λ$ increases, the super-harmonic resonance hump will be flattened. As fractional order p increases, the frequency range of primary resonance polymorphic coexistence reduces.

To sum up, this paper provides a fundamental insight into the design and optimization of the dynamic response of the viscoelastic vibration isolation system with a small amount of higher harmonics under micro-vibration. The research results will pay contributions to a better understanding of the effect of parameters on the super-harmonic resonance and periodic motion transition characteristics of fractional nonlinear vibration isolation system. When designing the vibration isolation systems in practical engineering, the appropriate system parameters should be selected to avoid the nonlinear jump phenomenon, bifurcation, and chaos in the dynamic response.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work has been supported by the Netural Science Foundation of China (11962013, 11732014, 12162020), Gansu Youth Science and Technology Fund (21JR7RA335).

ORCID iD

Minghe Qu

Appendix

The defined parameters for amplitude-frequency response in equation (19) are as follows

D_{13} = 12 ξ^{3} ε a_{1}^{2} b_{1} + 36 ξ^{3} ε a_{1} a_{3} b_{1} + 36 ξ^{3} ε b_{1}^{2} b_{3} + 108 ξ^{3} ε b_{1} b_{3}^{2} - 2 ξ b_{1}

\begin{array}{c} D_{12} = 6 ξ^{2} ε a_{1}^{3} μ_{k 2} + 15 ξ^{2} ε a_{1}^{2} a_{3} μ_{k 2} + 54 ξ^{2} ε a_{1} a_{3}^{2} μ_{k 2} + 60 ξ^{2} ε a_{1} b_{1} b_{3} μ_{k 2} + 6 ξ^{2} ε a_{1} b_{1}^{2} + 108 ξ^{2} ε a_{1} b_{3}^{2} \\ - 30 ξ^{2} ε a_{3} b_{1}^{2} - a_{1} μ_{k 2} + 54 ξ^{2} ε a_{1} b_{3}^{2} μ_{k 2} - 15 ξ^{2} ε a_{3} b_{1}^{2} μ_{k 2} + 6 ξ^{2} ε a_{1}^{3} + 30 ξ^{2} ε a_{1}^{2} a_{3} + 108 ξ^{2} ε a_{1} a_{3}^{2} \end{array}

D_{11} = - 6 ξ ε a_{1} a_{3} b_{1} μ_{k 2}^{2} + 3 ξ ε b_{1}^{3} μ_{k 2}^{2} - 12 ξ ε a_{1} a_{3} b_{1} μ_{k 2} + 6 ξ ε b_{1}^{3} μ_{k 2} + 3 ε a_{1}^{2} b_{1} ξ + 3 ξ ε a_{1}^{2} b_{3} + 3 ε b_{1}^{3} ξ + 6 ξ ε b_{1} b_{3}^{2} + 2 ξ b_{1} μ_{k 3} + α b_{1}

\begin{array}{c} D_{10} = \frac{3}{4} ε a_{1}^{3} μ_{k 2}^{3} + \frac{3}{2} ε a_{1} b_{1}^{2} μ_{k 2}^{3} + 3 ε a_{1} b_{1} b_{3} μ_{k 2}^{3} + \frac{9}{4} ε a_{1}^{3} μ_{k 2}^{2} + \frac{9}{2} ε a_{1} b_{1}^{2} μ_{k 2}^{2} + 2 a_{1} μ_{k 2} + 2 a_{1} μ_{k 3} + \frac{9}{2} ε a_{1}^{2} a_{3} μ_{k 2} \\ + 9 ε a_{1} b_{1} b_{3} μ_{k 2}^{2} + 9 ε a_{1} b_{3}^{2} μ_{k 2}^{2} - \frac{9}{2} ε a_{3} b_{1}^{2} μ_{k 2}^{2} + \frac{9}{4} ε a_{1}^{3} μ_{k 2} + 9 ε a_{1} b_{1} b_{3} μ_{k 2} + 9 ε a_{1} b_{3}^{2} μ_{k 2} - \frac{9}{2} ε a_{3} b_{1}^{2} μ_{k 2} \\ + 9 ε a_{1} a_{3}^{2} μ_{k 2} + \frac{9}{2} ε a_{1} b_{1}^{2} μ_{k 2} + \frac{3}{4} ε a_{1}^{3} + \frac{3}{4} ε a_{1}^{2} a_{3} + \frac{3}{2} ε a_{1} a_{3}^{2} + 3 ε a_{1} b_{1} b_{3} + β a_{1} μ_{k 2} + 2 a_{1} μ_{k 2} μ_{k 3} + β a_{1} \end{array}

D_{23} = - 6 ξ^{3} ε a_{1}^{3} + 18 ξ^{3} ε a_{1}^{2} a_{3} - 108 ξ^{3} ε a_{1} a_{3}^{2} - 6 ξ^{3} ε a_{1} b_{1}^{2} + 36 ξ^{3} ε a_{1} b_{1} b_{3} - 108 ξ^{3} ε a_{1} b_{3}^{2} - 18 ξ^{3} ε a_{3} b_{1}^{2} + 2 a_{1} ξ

\begin{array}{c} D_{22} = 3 ξ^{2} ε a_{1}^{2} b_{1} μ_{k 2} + 15 ξ^{2} ε a_{1}^{2} b_{3} μ_{k 2} - 30 ξ^{2} ε a_{1} a_{3} b_{1} μ_{k 2} + 54 ξ^{2} ε a_{3}^{2} b_{1} μ_{k 2} + 3 ξ^{2} ε b_{1}^{3} μ_{k 2} - b_{1} μ_{k 2} + 54 ξ^{2} ε b_{1} b_{3}^{2} \\ - 15 ξ^{2} ε b_{1}^{2} b_{3} μ_{k 2} + 54 ξ^{2} ε b_{1} b_{3}^{2} μ_{k 2} + 3 ξ^{2} ε a_{1}^{2} b_{1} + 15 ξ^{2} ε a_{1}^{2} b_{3} - 30 ξ^{2} ε a_{1} a_{3} b_{1} + 54 ξ^{2} ε a_{3}^{2} b_{1} - 15 ξ^{2} ε b_{1}^{2} b_{3} \\ + 3 ξ^{2} ε b_{1}^{3} - 2 α ξ b_{1} - b_{1} \end{array}

\begin{array}{c} D_{21} = - 3 ε a_{1} ξ a_{3}^{2} - 3 ξ ε a_{1} a_{3}^{2} μ_{k 2}^{2} - 3 ξ ε a_{1} b_{3}^{2} μ_{k 2}^{2} - \frac{3}{2} ε a_{1}^{3} ξ - \frac{3}{2} ξ a_{1}^{3} μ_{k 2}^{2} - 3 ξ a_{1}^{3} μ_{k 2} - \frac{3}{2} ε a_{1} b_{1}^{2} ξ - \frac{3}{2} ε a_{1}^{2} ε a_{3} \\ - 6 ξ ε a_{1} a_{3}^{2} μ_{k 2} - 6 ξ ε a_{1} b_{3}^{2} μ_{k 2} - 3 ξ ε a_{1} b_{1} b_{3} - a_{1} μ_{k 2} α - 2 a_{1} ξ β - \frac{3}{2} ξ ε a_{1}^{2} a_{3} μ_{k 2}^{2} - \frac{3}{2} ξ ε a_{1} b_{1}^{2} μ_{k 2}^{2} - α a_{1} \\ \begin{array}{c} + \frac{3}{2} ξ ε a_{3} b_{1}^{2} μ_{k 2}^{2} - 3 ξ ε a_{1}^{2} a_{3} μ_{k 2} - 3 ξ ε a_{1} b_{1}^{2} μ_{k 2} + 3 ξ ε a_{3} b_{1}^{2} μ_{k 2} + \frac{3}{2} ε b_{1}^{2} ξ a_{3} - 3 ξ ε a_{1} b_{3}^{2} - 6 ξ ε a_{1} b_{1} b_{3} μ_{k 2} \\ - 3 ξ ε a_{1} b_{1} b_{3} μ_{k 2}^{2} - 2 a_{1} ξ - 2 a_{1} μ_{k 3} ξ \end{array} \end{array}

\begin{array}{c} D_{20} = \frac{3}{4} ε b_{1}^{3} μ_{k 2}^{3} + \frac{9}{4} ε b_{1}^{3} μ_{k 2}^{2} + \frac{9}{4} ε b_{1}^{3} μ_{k 2} + \frac{3}{4} ε a_{1}^{2} b_{1} + \frac{3}{4} ε a_{1}^{2} b_{3} + \frac{3}{2} ε a_{3}^{2} b_{1} - \frac{3}{4} ε b_{1}^{2} b_{3} + \frac{3}{2} ε b_{1} b_{3}^{2} + β b_{1} μ_{k 2} \\ \begin{array}{c} + \frac{3}{4} ε a_{1}^{2} b_{1} μ_{k 2}^{3} + \frac{3}{4} ε a_{1}^{2} b_{3} μ_{k 2}^{3} + \frac{3}{2} ε a_{3}^{2} b_{1} μ_{k 2}^{3} - \frac{3}{4} ε b_{1}^{2} b_{3} μ_{k 2}^{3} + \frac{3}{2} ε b_{1} b_{3}^{2} μ_{k 2}^{3} + \frac{9}{4} ε a_{1}^{2} b_{1} μ_{k 2}^{2} + \frac{9}{4} ε a_{1}^{2} b_{3} μ_{k 2}^{2} \\ + \frac{9}{2} ε a_{3}^{2} b_{1} μ_{k 2}^{2} - \frac{9}{4} ε b_{1}^{2} b_{3} μ_{k 2}^{2} + \frac{9}{2} ε b_{1} b_{3}^{2} μ_{k 2}^{2} + \frac{9}{4} ε a_{1}^{2} b_{1} μ_{k 2} + \frac{9}{4} ε a_{1}^{2} b_{3} μ_{k 2} + \frac{9}{2} ε a_{3}^{2} b_{1} μ_{k 2} - \frac{9}{4} ε b_{1}^{2} b_{3} μ_{k 2} \\ - \frac{3}{2} ε a_{1} a_{3} b_{1} μ_{k 2}^{3} - \frac{9}{2} ε a_{1} a_{3} b_{1} μ_{k 2}^{2} - \frac{9}{2} ε a_{1} a_{3} b_{1} μ_{k 2} + b_{1} μ_{k 3} + \frac{3}{4} ε b_{1}^{3} + \frac{9}{2} ε b_{1} b_{3}^{2} μ_{k 2} - \frac{3}{2} ε a_{1} a_{3} b_{1} + β b_{1} \end{array} \\ + b_{1} μ_{k 2} + b_{1} μ_{k 2} μ_{k 3} \end{array}

D_{33} = - 6 ξ^{3} ε a_{1}^{2} b_{1} + 36 ξ^{3} ε a_{1}^{2} b_{3} + 162 b_{3} ξ^{3} ε a_{3}^{2} + 2 ξ^{3} b_{1}^{3} ε + 36 ξ^{3} ε b_{1}^{2} b_{3} + 162 ξ^{3} ε b_{3}^{3} - 54 ξ b_{3}

\begin{array}{c} D_{32} = - 3 ξ^{2} ε a_{1}^{3} μ_{k 2} + 6 ξ^{2} ε a_{1}^{2} a_{3} μ_{k 2} + 9 ε ξ^{2} a_{1} b_{1}^{2} μ_{k 2} + 27 ε ξ^{2} a_{3}^{3} μ_{k 2} + 6 ξ^{2} ε a_{3} b_{1}^{2} μ_{k 2} - 9 a_{3} \\ - 3 ξ^{2} ε a_{1}^{3} + 6 ξ^{2} ε a_{1}^{2} a_{3} + 9 ξ^{2} ε a_{1} b_{1}^{2} + 27 ε ξ^{2} a_{3}^{3} + 6 ξ^{2} ε a_{3} b_{1}^{2} + 27 ε ξ^{2} a_{3} b_{3}^{2} - 18 ξ α a_{3} \\ - 9 a_{3} μ_{k 2} + 27 ε ξ^{2} a_{3} b_{3}^{2} μ_{k 2} \end{array}

\begin{array}{c} D_{31} = \frac{9}{2} ξ ε a_{1}^{2} b_{1} μ_{k 2}^{2} + 9 ξ ε a_{1}^{2} b_{3} μ_{k 2}^{2} + \frac{9}{2} ξ ε a_{3}^{2} b_{3} μ_{k 2}^{2} + 9 ξ ε b_{1}^{2} b_{3} μ_{k 2}^{2} + 9 ξ ε a_{1}^{2} b_{1} μ_{k 2} + 18 ξ ε a_{1}^{2} b_{3} μ_{k 2} \\ + 9 ξ ε a_{3}^{2} b_{3} μ_{k 2} + 18 ξ ε b_{1}^{2} b_{3} μ_{k 2} + \frac{9}{2} ξ ε a_{3}^{2} b_{3} + 9 ξ ε b_{1}^{2} b_{3} + 9 ξ ε b_{3}^{3} μ_{k 2} + \frac{9}{2} ε a_{1}^{2} b_{1} ξ + 9 ε a_{1}^{2} b_{3} D_{30} = \frac{1}{4} ε \\ - \frac{3}{2} ξ ε b_{1}^{3} μ_{k 2}^{2} + \frac{9}{2} ξ ε b_{3}^{3} μ_{k 2}^{2} - 3 ξ ε b_{1}^{3} μ_{k 2} + 6 ξ b_{3} + 3 b_{3} α - \frac{3}{2} ε b_{1}^{3} ξ + 3 b_{3} μ_{k 2} α + \frac{9}{2} ξ ε b_{3}^{3} + 6 b_{3} ξ β + 6 b_{3} ξ μ_{k 3} \end{array}

\begin{array}{c} D_{30} = \frac{1}{4} ε a_{1}^{3} + \frac{3}{4} ε a_{3}^{3} + β a_{3} + a_{3} μ_{k 2} + a_{3} μ_{k 3} + \frac{1}{4} ε a_{1}^{3} μ_{k 2}^{3} + \frac{3}{4} ε a_{3}^{3} μ_{k 2}^{3} + \frac{3}{4} ε a_{1}^{3} μ_{k 2}^{2} + \frac{9}{4} ε a_{3}^{3} μ_{k 2}^{2} \\ \begin{array}{c} + \frac{3}{4} ε a_{1}^{3} μ_{k 2} + \frac{9}{4} ε a_{3}^{3} μ_{k 2} + \frac{3}{2} ε a_{1}^{2} a_{3} - \frac{3}{4} ε a_{1} b_{1}^{2} + \frac{3}{2} ε a_{3} b_{1}^{2} + \frac{3}{4} ε a_{3} b_{3}^{2} + β a_{3} μ_{k 2} + a_{3} μ_{k 2} μ_{k 3} \\ + \frac{3}{2} ε a_{1}^{2} a_{3} μ_{k 2}^{3} - \frac{3}{4} ε a_{1} b_{1}^{2} μ_{k 2}^{3} + \frac{3}{2} ε a_{3} b_{1}^{2} μ_{k 2}^{3} + \frac{3}{4} ε a_{3} b_{3}^{2} μ_{k 2}^{3} + \frac{9}{2} ε a_{1}^{2} a_{3} μ_{k 2}^{2} - \frac{9}{4} ε a_{1} b_{1}^{2} μ_{k 2}^{2} \end{array} \\ + \frac{9}{4} ε a_{3} b_{3}^{2} μ_{k 2}^{2} + \frac{9}{2} ε a_{1}^{2} a_{3} μ_{k 2} - \frac{9}{4} ε a_{1} b_{1}^{2} μ_{k 2} + \frac{9}{4} ε a_{3} b_{3}^{2} μ_{k 2} + \frac{9}{2} ε a_{3} b_{1}^{2} μ_{k 2} + \frac{9}{2} ε a_{3} b_{1}^{2} μ_{k 2}^{2} \end{array}

D_{43} = 2 ξ^{3} ε a_{1}^{3} - 36 ξ^{3} ε a_{1}^{2} a_{3} - 6 ξ^{3} ε a_{1} b_{1}^{2} - 162 ξ^{3} ε a_{3}^{3} - 36 ξ^{3} ε a_{3} b_{1}^{2} - 162 a_{3} ξ^{3} ε b_{3}^{2} + 54 ξ a_{3}

\begin{array}{c} D_{42} = - 9 ξ^{2} ε a_{1}^{2} b_{1} μ_{k 2} + 6 ξ^{2} ε a_{1}^{2} b_{3} μ_{k 2} + 27 ξ^{2} ε a_{3}^{2} b_{3} μ_{k 2} + 3 ξ^{2} ε b_{1}^{3} μ_{k 2} - 18 ξ α b_{3} - 9 b_{3} (μ_{k 2} + 1) \\ + 6 ξ^{2} ε a_{1}^{2} b_{3} + 3 ξ^{2} ε b_{1}^{3} + 6 ξ^{2} ε b_{1}^{2} b_{3} + 27 ξ^{2} ε b_{3}^{3} + 27 ξ^{2} ε b_{3}^{3} μ_{k 2} - 9 ξ^{2} ε a_{1}^{2} b_{1} + 6 ξ^{2} ε b_{1}^{2} b_{3} μ_{k 2} \\ + 27 ξ^{2} ε a_{3}^{2} b_{3} \end{array}

\begin{array}{c} D_{41} = - \frac{9}{2} ξ ε a_{3} b_{3}^{2} μ_{k 2}^{2} - 9 ξ ε a_{3} b_{3}^{2} μ_{k 2} - \frac{9}{2} ξ ε a_{3}^{3} μ_{k 2}^{2} - 9 ξ ε a_{3}^{3} μ_{k 2} - \frac{3}{2} ξ ε a_{1}^{3} μ_{k 2}^{2} - 18 ξ ε a_{3} b_{1}^{2} μ_{k 2} \\ - 3 ξ ε a_{1}^{3} μ_{k 2} - 9 ε a_{1}^{2} ξ a_{3} + \frac{9}{2} ε a_{1} b_{1}^{2} ξ - 9 ε b_{1}^{2} ξ a_{3} - 6 ξ a_{3} - 3 α a_{3} - 3 a_{3} μ_{k 2} α + 9 ξ ε a_{1} b_{1}^{2} μ_{k 2} \\ \begin{array}{c} - 6 a_{3} ξ β - 6 μ_{k 3} ξ a_{3} - \frac{9}{2} ε ξ a_{3}^{3} - 9 ξ ε a_{1}^{2} a_{3} μ_{k 2}^{2} - 9 ξ ε a_{3} b_{1}^{2} μ_{k 2}^{2} - 18 ξ ε a_{1}^{2} a_{3} μ_{k 2} \\ - \frac{3}{2} ε a_{1}^{3} ξ + \frac{9}{2} ξ ε a_{1} b_{1}^{2} μ_{k 2}^{2} \end{array} \end{array}

\begin{array}{c} D_{40} = - \frac{1}{4} ε b_{1}^{3} μ_{k 2}^{3} - \frac{3}{4} ε b_{1}^{3} μ_{k 2}^{2} - \frac{3}{4} ε b_{1}^{3} μ_{k 2} + \frac{3}{4} ε a_{1}^{2} b_{1} + \frac{3}{2} ε a_{1}^{2} b_{3} + \frac{3}{2} ε b_{1}^{2} b_{3} + \frac{3}{4} ε b_{3}^{3} + β b_{3} + b_{3} μ_{k 2} \\ + b_{3} μ_{k 3} - \frac{1}{4} ε b_{1}^{3} + \frac{3}{4} ε b_{3}^{3} μ_{k 2}^{3} + \frac{9}{4} ε b_{3}^{3} μ_{k 2}^{2} + \frac{9}{4} ε b_{3}^{3} μ_{k 2} + \frac{3}{4} ε a_{3}^{2} b_{3} + β b_{3} μ_{k 2} + b_{3} μ_{k 2} μ_{k 3} \\ \begin{array}{c} + \frac{3}{4} ε a_{1}^{2} b_{1} μ_{k 2}^{3} + \frac{3}{2} ε a_{1}^{2} b_{3} μ_{k 2}^{3} + \frac{3}{2} ε b_{1}^{2} b_{3} μ_{k 2}^{3} + \frac{9}{4} ε a_{1}^{2} b_{1} μ_{k 2}^{2} + \frac{9}{2} ε a_{1}^{2} b_{3} μ_{k 2}^{2} + \frac{9}{2} ε b_{1}^{2} b_{3} μ_{k 2}^{2} \\ + \frac{9}{4} ε a_{1}^{2} b_{1} μ_{k 2} + \frac{9}{2} ε a_{1}^{2} b_{3} μ_{k 2} + \frac{9}{2} ε b_{1}^{2} b_{3} μ_{k 2} + \frac{3}{4} ε a_{3}^{2} b_{3} μ_{k 2}^{3} + \frac{9}{4} ε a_{3}^{2} b_{3} μ_{k 2}^{2} + \frac{9}{4} ε a_{3}^{2} b_{3} μ_{k 2} \end{array} \end{array}

References

Okuka

Zorica

. Fractional burgers models in creep and stress relaxation tests. Appl Math Modell 2020; 77: 1894–1935. DOI: 10.1016/j.apm.2019.09.035.

Berg

. Nonlinear rubber spring model for vehicle dynamics analysis. Veh Syst Dyn 1998; 29: 723–728. DOI: 10.1080/00423119808969599.

Dzierzek

. Experiment-based modeling of cylindrical rubber bushings for the simulation of wheel suspension dynamic behavior. SAE Trans 2000; 109: 78–85.

Miyasato

Simionatto

VGS

Dias

. Simulation of a mass-on-belt dynamical model with the Zener viscoelastic support. J Sound Vib 2022; 534: 117025. DOI: 10.1016/j.jsv.2022.117025.

Zhou

Zhang

. Dynamic modeling and computer simulation analysis of the air spring suspension. J Mech Sci Technol 2022; 36: 1719–1727. DOI: 10.1007/s12206-022-0308-2.

Zheng

Shangguan

W-B

Rakheja

. Modeling and analysis of time-domain nonlinear characteristics of air spring with an auxiliary chamber. Mech Syst Sig Process 2022; 176: 109161. DOI: 10.1016/j.ymssp.2022.109161.

Zhu

Taiki

Heisha

, et al. Investigation on low-frequency micro-vibration model of air spring with auxiliary chamber for steel plate concrete composite structures. Front Mater 2022; 9: 896497. DOI: 10.3389/fmats.2022.896497.

Yin

, et al. A new dynamic stiffness model with hysteresis of air springs based on thermodynamics. J Sound Vib 2022; 521: 116693. DOI: 10.1016/j.jsv.2021.116693.

Larsson

Singh

Olsson

, et al. A micromechanically based model for dynamic damage evolution in unidirectional composites. Int J Sol Struct 2022; 238: 111368. DOI: 10.1016/j.ijsolstr.2021.111368.

10.

Oda

Nishimura

. Vibration of air suspension bogies and their design. Bull JSME 1970; 13: 43–50. DOI: 10.1299/jsme1958.13.43.

11.

Ryaboy

. Static and dynamic stability of pneumatic vibration isolators and systems of isolators. J Sound Vib 2014; 333: 31–51. DOI: 10.1016/j.jsv.2013.06.034.

12.

Yang

Wang

Wei

, et al. Investigation on nonlinear and fractional derivative Zener model of coupled vehicle-track system. Veh Syst Dyn 2019; 58: 864–889. DOI: 10.1080/00423114.2019.1602273.

13.

Mendia-Garcia

Gil-Negrete Laborda

Pradera-Mallabiabarrena

, et al. A survey on the modelling of air springs – secondary suspension in railway vehicles. Veh Syst Dyn 2020; 60: 835–864. DOI: 10.1080/00423114.2020.1838566.

14.

Sun

Zhang

Baleanu

, et al. A new collection of real world applications of fractional calculus in science and engineering. Commun Nonlinear Sci Numer Simul 2018; 64: 213–231. DOI: 10.1016/j.cnsns.2018.04.019.

15.

Machado

Kiryakova

Mainardi

. Recent history of fractional calculus. Commun Nonlinear Sci Numer Simul 2011; 16: 1140–1153. DOI: 10.1016/j.cnsns.2010.05.027.

16.

Shen

Y-J

Wen

S-F

X-H

, et al. Dynamical analysis of fractional-order nonlinear oscillator by incremental harmonic balance method. Nonlinear Dyn 2016; 85: 1457–1467. DOI: 10.1007/s11071-016-2771-8.

17.

Wang

K-L

Liu

S-Y

. He’s fractional derivative and its application for fractional Fornberg-Whitham equation. Therm Sci 2017; 21: 2049–2055. DOI: 10.2298/tsci151025054w.

18.

Wang

YAN

Deng

. Fractal derivative model for tsunami traveling. Fractals 2019; 27: 1950017. DOI: 10.1142/s0218348x19500178.

19.

C-H

Shen

YUE

F-Y

, et al. Taylor series solution for fractal Bratu-type equation arising in electrospinning process. Fractals 2020; 28: 2050011. DOI: 10.1142/s0218348x20500115.

20.

Wang

K-L

Wang

K-J

C-H

. Physical insight of local fractional calculus and its application to fractional kdv–burgers–kuramoto equation. Fractals 2019; 27: 1950122. DOI: 10.1142/s0218348x19501226.

21.

J-H

El-Dib

. A tutorial introduction to the two-scale fractal calculus and its application to the fractal Zhiber–Shabat oscillator. Fractals 2021; 29: 1–9. DOI: 10.1142/s0218348x21502686.

22.

C-H

Liu

J-H

, et al. Low frequency property of a fractal vibration model for a concrete beam. Fractals 2021; 29: 2150117. DOI: 10.1142/s0218348x21501176.

23.

C-H

Liu

. A modified frequency–amplitude formulation for fractal vibration systems. Fractals 2022; 30: 2250046. DOI: 10.1142/s0218348x22500463.

24.

J-H

Kou

S-J

C-H

, et al. Fractal oscillation and its frequency-amplitude property. Fractals 2021; 29: 2150105. DOI: 10.1142/s0218348x2150105x.

25.

Khajehsaeid

. A comparison between fractional-order and integer-order differential finite deformation viscoelastic models: effects of filler content and loading rate on material parameters. Int J Appl Mech 2019; 10: 1850099. DOI: 10.1142/s1758825118500990.

26.

Luo

Shi

Guo

, et al. A nonlinear rubber spring model for the dynamics simulation of a high-speed train. Veh Syst Dyn 2019; 58: 1367–1384. DOI: 10.1080/00423114.2019.1624788.

27.

Bhangale

Kachhia

Gómez‐Aguilar

. Fractional viscoelastic models with Caputo generalized fractional derivative. Math Methods Appl Sci 2021. DOI: 10.1002/mma.7229.

28.

Chen

J-J

Yin

Z-H

Yuan

X-J

, et al. A refined stiffness model of rolling lobe air spring with structural parameters and the stiffness characteristics of rubber bellows. Measurement 2021; 169: 108355. DOI: 10.1016/j.measurement.2020.108355.

29.

Chen

J-J

Yin

Z-H

Rakheja

, et al. Theoretical modelling and experimental analysis of the vertical stiffness of a convoluted air spring including the effect of the stiffness of the bellows. P Mech Eng D-J Aut 2017; 232: 547–561. DOI: 10.1177/0954407017704589.

30.

Zhu

Yang

Zhang

, et al. Nonlinear dynamic model of air spring with a damper for vehicle ride comfort. Nonlinear Dyn 2017; 89: 1545–1568. DOI: 10.1007/s11071-017-3535-9.

31.

. Research on nonlinear modeling and dynamic characteristics of lateral stiffness of vehicle air spring system. Adv Mech Eng 2020; 12: 168781402093045. DOI: 10.1177/1687814020930457.

32.

Zhang

Y-W

Ding

. Vibration control combining nonlinear isolation and nonlinear absorption. Nonlinear Dyn 2020; 100: 2121–2139. DOI: 10.1007/s11071-020-05606-6.

33.

Liu

Liao

, et al. Enhanced vibration isolation performance of quasi-zero-stiffness isolator by introducing tunable nonlinear inerter. Commun Nonlinear Sci Numer Simul 2021; 95: 105654. DOI: 10.1016/j.cnsns.2020.105654.

34.

Liu

Peng

Jin

. Design and characteristics of a novel QZS vibration isolation system with origami-inspired corrector. Nonlinear Dyn 2021; 106: 255–277. DOI: 10.1007/s11071-021-06821-5.

35.

Miguel

de Oliveira Teloli

da Silva

. Some practical regards on the application of the harmonic balance method for hysteresis models. Mech Syst Sig Process 2020; 143: 106842. DOI: 10.1016/j.ymssp.2020.106842.

36.

Pei

Chong

ASE

Pavlovskaia

, et al. Computation of periodic orbits for piecewise linear oscillator by harmonic balance methods. Commun Nonlinear Sci Numer Simul 2022; 108: 106220. DOI: 10.1016/j.cnsns.2021.106220.

37.

Chen

Xing

Wang

, et al. Study on chaos of nonlinear suspension system with fractional-order derivative under random excitation. Chaos Soliton Fract 2021; 152: 111300. DOI: 10.1016/j.chaos.2021.111300.

38.

. Theoretical modeling of the vertical stiffness of a rolling lobe air spring. Sci Prog 2020; 103: 36850420940898. DOI: 10.1177/0036850420940898.

39.

Chen

, et al. Nonlinear dynamic analysis and active control of visco-hyperelastic dielectric elastomer membrane. Int J Sol Struct 2018; 152–153: 28–38. DOI: 10.1016/j.ijsolstr.2018.05.006.

40.

Sharma

Joglekar

. A computationally efficient locking free numerical framework for modeling visco-hyperelastic dielectric elastomers. Comput Methods Appl Mech Eng 2019; 352: 625–653. DOI: 10.1016/j.cma.2019.04.035.

41.

Y-F

Peng

Y-H

. Stability analysis of fractional nonlinear dynamic systems with order lying in (1, 2). Chin Quart J Math 2019; 34: 188–195. DOI: 10.13371/j.cnki.chin.q.j.m.2019.02.006.

42.

Shen

Yang

Xing

, et al. Primary resonance of Duffing oscillator with fractional-order derivative. Commun Nonlinear Sci Numer Simul 2012; 17: 3092–3100. DOI: 10.1016/j.cnsns.2011.11.024.

43.

Shen

Y-J

Wei

Yang

S-P

. Primary resonance of fractional-order van der Pol oscillator. Nonlinear Dyn 2014; 77: 1629–1642. DOI: 10.1007/s11071-014-1405-2.

44.

Zhou

Song

, et al. Dynamic and steady analysis of a 2-DOF vehicle system by modified incremental harmonic balance method. Nonlinear Dyn 2019; 98: 75–94. DOI: 10.1007/s11071-019-05172-6.

45.

Jena

Chakraverty

Jena

, et al. Analysis of time‐fractional fuzzy vibration equation of large membranes using double parametric based residual power series method. Z Angew Math Mech 2020; 101: e202000165. DOI: 10.1002/zamm.202000165.

46.

Matouk

. A novel fractional-order system: chaos, hyperchaos and applications to linear control. J Appl Comput Mech 2021; 7: 701–714. DOI: 10.22055/JACM.2020.35092.2561.

47.

Zhao

Jiang

Zhang

, et al. Study on optimization of damping performance and damping temperature range of silicone rubber by polyborosiloxane gel. Polymers (Basel) 2020; 12: 1196. DOI: 10.3390/polym12051196.

48.

Liu

Niu

Shen

, et al. Stability and bifurcation analysis of two-degrees-of-freedom vibro-impact system with fractional-order derivative. Int J Non Linear Mech 2020; 126: 103570. DOI: 10.1016/j.ijnonlinmec.2020.103570.

49.

Hou

Shang

D-G

Zuo

L-X

, et al. Fatigue life prediction of needled ceramic matrix composite under variable amplitude loading. Int J Fatigue 2022; 156: 106690. DOI: 10.1016/j.ijfatigue.2021.106690.

50.

Zhang

Wei

. Prediction of variable amplitude fatigue crack growth life based on modified grey model. Eng Fail Anal 2022; 133: 105939. DOI: 10.1016/j.engfailanal.2021.105939.

51.

Liu

. Superharmonic resonance of the quasi-zero-stiffness vibration isolator and its effect on the isolation performance. Nonlinear Dyn 2020; 100: 95–117. DOI: 10.1007/s11071-020-05509-6.

52.

Yan

Liu

. Fractional-order harmonic resonance in a multi-frequency excited fractional Duffing oscillator with distributed time delay. Commun Nonlinear Sci Numer Simul 2021; 97: 105754. DOI: 10.1016/j.cnsns.2021.105754.

53.

Shahraeeni

Sorokin

Mace

, et al. Effect of damping nonlinearity on the dynamics and performance of a quasi-zero-stiffness vibration isolator. J Sound Vib 2022; 526: 116822. DOI: 10.1016/j.jsv.2022.116822.

54.

Ding

Wang

Ding

, et al. Periodic motion and transition of a vibro‐impact system with multilevel elastic constraints. Discrete Dyn Nat Soc 2021; 2021: 1–13. DOI: 10.1155/2021/6687887.

55.

G F

Sun

Ding

W C

. Dynamics of a vibro-impact system by the global analysis method in parameter-state space. Nonlinear Dyn 2019; 97: 541–557. DOI: 10.1007/s11071-019-04996-6.

Analysis of super-harmonic resonance and periodic motion transition of fractional nonlinear vibration isolation system

Abstract

Keywords

Introduction

Analysis of fractional nonlinear Nishimura model

Amplitude–frequency characteristics

The diversity of periodic motions of super-harmonic response

Analysis of bifurcation characteristics

Conclusion

Footnotes

Declaration of conflicting interests

Funding

ORCID iD

Appendix

References