Nonlinear Dynamics Study of Giant Magnetostrictive Actuators with Fractional Damping

Abstract

Featured Application

This study contributes to a better understanding of the dynamics of GMA systems and provides a new perspective for controlling the stable operation of GMA systems in engineering practice.

Abstract

Since the structural mechanics of the super magnetostrictive actuator (GMA) system involves problems related to viscoelastic damping materials, the fractional order is more accurate than the integer order calculus to characterize the viscoelastic features in the structure. In order to further investigate the intrinsic mechanism and dynamical characteristics of the GMA dynamical system, the dynamical equations of the nonlinear GMA system containing fractional damping terms are established and the main resonance of the system is analyzed using the averaging method. The mechanism of the influence of some parameters on the GMA system is analyzed by MATLAB numerical simulation to study the bifurcation and chaotic motion phenomena of the system from the qualitative and quantitative perspectives. The results show that the fractional damping coefficient, external excitation amplitude and fractional order have significant effects on the amplitude-frequency characteristics of the system; the fractional order has a greater influence on the bifurcation and chaotic behavior of the system; the dynamic behavior of the system caused by the change of external excitation amplitude and fractional damping coefficient at different damping orders is similar but the chaotic region is different.

1. Introduction

Giant magnetostrictive material (GMM) is a strategic rare earth functional material [1]. The giant magnetostrictive actuator (GMA) made from GMM rods as the core component has the advantages of a high magneto-mechanical coupling coefficient, good low frequency response, high magnetostrictive coefficient, and a high output force and wide vibration response band. It has been widely used in the fields of micro-displacement drive, sensor, precision positioning control and vibration control [2,3,4,5,6,7]. However, due to the influence of multiple coupled nonlinear factors during the operation of the GMA system [8,9], the nonlinearity of GMA structural materials, i.e., the stress-strain relationship does not satisfy Hooke’s law, and the nonlinearity of the system’s construction, is reflected in the nonlinearity of deformation and external forces, and the recovery force and damping are nonlinear functions of structural vibration displacement and velocity, respectively. In addition, the GMM rod also has strong hysteresis characteristics, frequency doubling characteristics, and eddy current effect nonlinear characteristics; it is easy to fall into nonlinear instability, resulting in large output error and low control accuracy of the system. At present, studies on enhancing and improving the stability of hysteresis nonlinearities in GMA systems have mainly focused on the hysteresis compensation control of the system [10,11], and the mechanism of system dynamics instability and dynamics characteristics have not been studied in depth.

The existence of nonlinear instability and even chaotic motion in GMA systems makes it difficult to accurately predict and control their output response, which seriously hinders the application of GMA in precision positioning control and active vibration isolation, etc. Zeng et al. [12] established a high-power GMA nonlinear mathematical model and showed that increasing the system damping coefficient helps to improve output stability, and chaotic phenomena occur when the system’s stiffness is low. Gao et al. [13] established a mathematical model of GMA electromagnetic coupling dynamics and showed that when the system stiffness coefficient and damping coefficient are low, it leads to system instability. Yan et al. [14] established a mathematical model of GMA hysteresis nonlinearity and showed that the system has a complex motion pattern under different structural parameters. Sylvain et al. [15] established a mathematical model containing a five times nonlinear thin magnetostrictive actuator mathematical model, and it was shown that the system exhibits rich and complex dynamical properties, such as multistability and anti-monotonicity. In order to study the structural vibrations [16,17,18,19] in a system, it is essential to establish an accurate theoretical model. Nam et al. [20] established a new model of beam, this new beam model is free of shear-locking for both thick and thin beams, is easy to apply in computation, and has efficiency in simulating the variable thickness beams. Gebrel et al. [21] developed a suitable theoretical model to generate nonlinear electrostatic forces acting on the MEMS ring structure and examined the nonlinear dynamic response in the drive and sensing directions by time response, phase diagram and Poincare mapping when both input angular motion and nonlinear electrostatic forces were considered. Although some progress has been made in the study of chaos in GMA systems, there is a lack of more accurate dynamical model descriptions for GMA systems generating chaotic phenomena.

By studying the relationship between fractional order calculus and integer order calculus, and the properties of fractional order calculus, it can be found that the fractional order calculus integrates the effects of historical as well as nonlocal distributions and does not reflect only local and some point properties. It is this long memory property and global nature that the fractional order calculus more aptly describes; the damping properties of viscoelastic materials [22,23,24,25,26]. Coccolo et al. [22] considered nonlinear Duffing oscillators in the presence of fractional damping and showed that oscillations exist in the fractional order derivatives and that their amplitude and asymptotic time change abruptly with small changes in the fractional parameter. Li et al. [23] proposed a model of oblique crack dynamics based on fractional order damping and experimentally verified that the proposed fractional order model compensates for the deficiencies of the integer order model. Zarraga et al. [24] investigated the distribution of poles and the dynamic response to several excitations for different values of the parameters of the fractional damping model, and showed that the fractional model differs from the conventional viscous damping system by observing specific behaviors that cannot be replicated by the classical behavior. Sylvain et al. [25] studied the dynamics of fractional-order thin magnetostrictive actuators and showed that, compared to integer order systems, fractional order systems have more complex chaotic properties and extreme multistability. Yang et al. [26] studied a frictionally damped system with fractional order derivative damping under Gaussian white noise excitation, analyzed the stochastic bifurcation phenomenon caused by system parameters such as fractional order and fractional order coefficients, by setting critical conditions, and finally gave the differences between the bifurcation regions of the fractional order model and the integer order model. So far, most of the studies on GMA dynamical systems are based on the mechanical behavior described by integer-order calculus, which has some limitations for studying the stability of the system. The above studies show that it is important to consider fractional-order calculus in the model of dynamical systems and that GMM is a viscoelastic material with memory [25,27], and for systems with viscoelastic materials, the damping term depends not only on the current time and position but also on the previous state, and the memory properties of the viscoelastic material are characterized by the power-law kernel function associated with the fractional-order derivatives. The use of the fractional order model to describe the damping characteristics of GMA systems is essential to provide theoretical guidance for the stability design of GMA structures, not only to further reveal the intrinsic mechanism during the nonlinear motion of the system, but also to study the stability of GMA systems from a new perspective.

In this paper, based on the working mechanism of the GMA system based on previous research, the GMA dynamical system model is extended to fractional order based on fractional order calculus theory, the dynamical equations of nonlinear GMA system containing fractional damping terms are established, and the nonlinear vibration response of the system is analyzed based on the averaging method; the numerical solution of the system is solved using the algorithm defined by Grunwald–Letnikov fractional order derivatives, and the effects of fractional damping coefficient $\mathsf{\delta },$ and different external excitation amplitudes f, and fractional order $\mathsf{\alpha },$ on the main resonant output response of the GMA system are studied through MATLAB numerical simulations. By varying the fractional order $\mathsf{\alpha }$ and the magnitude of fractional damping coefficient $\mathsf{\delta }$ and excitation amplitude f at different orders, the bifurcation and chaos characteristics of the GMA system are investigated from qualitative and quantitative perspectives using bifurcation diagrams, Lyapunov exponent diagrams, phase trajectory diagrams, amplitude spectrum diagrams, Poincaré diagrams and time domain waveform diagrams.

2. Traditional Mechanical Model of the GMA System

Figure 1 shows the structure of the GMA system. The GMA working principle: the permanent magnet generates bias magnetic field, so that the driver can avoid the problem of frequency doubling; the disc spring and output rod sleeve form the preload structure, which generates the axial preload to make the GMM bar work under the compressive stress condition; the magnetic guide ring, magnetic guide sleeve and magnetic guide sheet can form the closed magnetic circuit, which can improve the magnetic field utilization and reduce the leakage, and can also reduce the external environment on the internal magnetic field. The influence of the external environment on the internal magnetic field can also be reduced. The drive coil generates a magnetic field through the current, which causes the GMM bar to change its length and push the output rigid rod to move in the form of output shaft force to realize the displacement and force output.

Considering the connection stiffness of the output part and the stiffness characteristics of the pre-pressed disc spring, the output of the GMA system is equated to a single degree of freedom k-m-c system, and its simplified model is shown in Figure 2. The upper half of the simplified model is a cross-shaped output rod, and the lower half is a GMM rod with a lower end face displacement of 0. Where N is the output axis force generated by the GMM rod under the excitation magnetic field, F and $\mathsf{\omega }$ are the amplitude and frequency of the external excitation signal on the output rod, c is the system damping factor, k is the coefficient of elasticity of the disc spring, u and i are the voltage and current of the drive coil, respectively, and $\ddot{\mathrm{x}}$ is the acceleration of the output rod. The spring module consists of the pre-pressed disc spring and the GMM bar in the GMA system, and the system damping includes the internal damping of the GMM bar and the damping of the output bar.

The forces on the output rod in the GMA system are analyzed as shown in Figure 3. Where K(x) is the prestress applied by the disc spring, and C($\dot{\mathrm{x}}$) is the damping force on the output rod. Under the action of static uniform magnetic field, the upper face of the GMM rod and the output rod have the same displacement x, velocity $\dot{\mathrm{x}}$, and acceleration $\ddot{\mathrm{x}}$ between them.

Under the action of a static stable magnetic field, according to the first piezomagnetic equation, the magnetostrictive effect, the linear piezomagnetic equation along the length direction of the GMM bar can be obtained as follows:

(1)$\left\{\begin{array}{c}{\mathrm{B}=\mathrm{d}}_{33}{\mathrm{T}}_3{+\mathsf{\mu }}_{33}^{\mathrm{T}}\mathrm{H}\\ {\mathsf{\zeta }=\mathrm{s}}_{33}^{\mathrm{T}}{\mathrm{T}}_3{+\mathrm{d}}_{33}\mathrm{H}\end{array}\right.$

where B is the magnetic induction intensity, d₃₃ is the piezomagnetic constant, T₃ is the stress on the GMM bar, ${\mathsf{\mu }}_{33}^{\mathrm{T}}$ is the magnetic permeability, H is the magnetic field strength inside the GMM bar, and ${\mathrm{s}}_{33}^{\mathrm{T}}$ is the flexibility factor.

According to Newton’s second law, the force analysis of the output rod, the dynamics equation of the GMA system, is obtained as follows:

(2)$\mathrm{m}\ddot{\mathrm{x}}=\mathrm{N}+\mathrm{Fcos}\mathsf{\omega }\mathrm{t}-\mathrm{K}-\mathrm{C}$

where m is the mass of the output rod.

The displacement of the output lever of the GMA system is as follows:

(3)$\mathrm{x}=\mathsf{\zeta }\mathrm{l}$

where $\mathsf{\zeta }$ is the strain of the GMM bar, and l is the length of the GMM bar.

The output axial force of GMA in a static uniform magnetic field without considering the effect of magnetic field variation on each parameter and the total magnetoresistance of the magnetic circuit is as follows:

(4)$\mathrm{N}=\frac{{\mathrm{E}}_{\mathrm{y}}^{\mathrm{H}}\mathrm{S}}{\mathrm{l}}(1-\frac{{\mathrm{E}}_{\mathrm{y}}^{\mathrm{H}}{\mathrm{d}}_{33}^2}{{\mathsf{\mu }}_{33}^{\mathrm{T}}})\mathrm{x}$

where ${\mathrm{E}}_{\mathrm{y}}^{\mathrm{H}}$ is the modulus of elasticity of the GMM bar, and S is the cross-sectional area of the GMM bar.

The GMA output lever is subjected to a system damping force as follows:

(5)$\mathrm{C}=\mathrm{c}\dot{\mathrm{x}}$

According to the disc spring design theory and the magnetostriction theory, so that the disc spring deformation ${\mathrm{f}}_{\mathrm{z}}=\mathrm{x}$, the preload stress K(x) of the preloaded disc spring in GMA can be expressed as follows:

(6)$\mathrm{K}\left(\mathrm{x}\right)=\frac{{\mathrm{h}}_0^2{\mathrm{h}+\mathrm{h}}^3}{{\mathsf{\chi }\mathrm{D}}^2}\mathrm{x}+\frac{{3\mathrm{h}}_0}{{2\mathsf{\chi }\mathrm{D}}^2}{\mathrm{x}}^2+\frac{\mathrm{h}}{{2\mathsf{\chi }\mathrm{D}}^2}{\mathrm{x}}^2$

where $\mathsf{\chi }$ is the correction value obtained from the table of D/d value, h is the thickness of the disc spring, D is the outer diameter of the disc spring, and ${\mathrm{h}}_0$ is the maximum deformation of the disc spring.

Simplify Equation (6) as follows:

(7)$\mathrm{K}\left(\mathrm{x}\right){=\mathrm{k}}_0{\mathrm{x}+\mathsf{\beta }}_1{\mathrm{x}}^2{+\mathsf{\gamma }}_1{\mathrm{x}}^2$

where ${\mathrm{k}}_0=\frac{{\mathrm{h}}_0^2{\mathrm{h}+\mathrm{h}}^3}{{\mathsf{\chi }\mathrm{D}}^2}$, ${\mathsf{\beta }}_1=\frac{{3\mathrm{h}}_0}{{2\mathsf{\chi }\mathrm{D}}^2}$, ${\mathsf{\gamma }}_1=\frac{\mathrm{h}}{{2\mathsf{\chi }\mathrm{D}}^2}$.

Substituting Equations (4), (5) and (7) into Equation (2), the magneto-mechanical coupling dynamics equation of the GMA system is obtained as:

(8)$\mathrm{m}\ddot{\mathrm{x}}{+\mathrm{c}\dot{\mathrm{x}}+\mathrm{k}}_1{\mathrm{x}+\mathsf{\beta }}_1{\mathrm{x}}^2{+\mathsf{\gamma }}_1{\mathrm{x}}^2 -\mathrm{Fcos}\mathsf{\omega }\mathrm{t}=0$

where ${\mathrm{k}}_1{=\mathrm{k}}_0+\frac{{\mathrm{E}}_{\mathrm{y}}^{\mathrm{H}}\mathrm{S}}{\mathrm{l}}(\frac{{\mathrm{E}}_{\mathrm{y}}^{\mathrm{H}}{\mathrm{d}}_{33}^2}{{\mathsf{\mu }}_{33}^{\mathrm{T}}} -1)$.

3. Mechanical Model of the GMA System with Fractional Damping Term

It is well known that the step response depends on the structural parameters of the system itself, which can well express the influence of the structural parameters of the system on its dynamic performance. Wang [28] derived the system transfer function of the GMA system with the current as the input and the displacement as the output, and the specific transfer function expression was obtained by substituting the system model parameters in this literature as follows:

(9)$\mathrm{G}\left(\mathrm{s}\right)=\frac{\mathrm{x}\left(\mathrm{s}\right)}{\mathrm{I}\left(\mathrm{s}\right)}=\frac{137,170}{{19\mathrm{s}}^2+2{.48\times 10}^5\mathrm{s}+2{.22\times 10}^9}$

Figure 4 shows the step response curves of the GMA system measured by the experiment and the theoretical curve of the simulation. From the Figure, it can be seen that the step response curves of the system obtained from the experiment and the theory not only have some errors in the overshoot and peak time, but also have large differences in the curve overlap.

Since the fractional order derivatives of the Grunwald–Letnikov form are convenient for discretization and easy to perform numerical operations, it is appropriate to include the fractional damping term defined in its form for a more accurate study of the stability of the GMA system. The defining equation of the Grunwald–Letnikov derivative of the function f(t) is as follows:

(10)$\left\{\begin{array}{c}{\mathrm{P}=\mathsf{\delta }}_1\frac{{\mathrm{d}}^{\mathsf{\alpha }}}{{\mathrm{dt}}^{\mathsf{\alpha }}}\\ {\left.\frac{{\mathrm{d}}^{\mathsf{\alpha }}\mathrm{f}\left(\mathrm{t}\right)}{{\mathrm{dt}}^{\mathsf{\alpha }}}\right|}_{{\mathrm{t}=\mathrm{k}}_{\mathrm{n}}\Delta \mathrm{t}}=\underset{\Delta \mathrm{t}\to 0}{\lim }\frac{1}{∆{\mathrm{t}}^{\mathsf{\alpha }}}{\displaystyle {\displaystyle \sum }_{\mathrm{j}=0}^{{\mathrm{k}}_{\mathrm{n}}}}{(-1)}^{\mathrm{j}}\left(\begin{array}{c}\mathsf{\alpha }\\ \mathrm{j}\end{array}\right){\mathrm{f}(\mathrm{k}}_{\mathrm{n}}\Delta \mathrm{t}-\mathrm{j}\Delta \mathrm{t})\end{array}\right.$

where $\mathsf{\alpha }$ is the order of derivation, $∆\mathrm{t}$ is the calculated time step, ${\mathsf{\delta }}_1$ is the fractional damping coefficient, considering the actual engineering context of the GMA system, and $\mathsf{\alpha }$ is taken as $0<\mathsf{\alpha }<2$. The binomial coefficients are as follows:

(11)$\left(\begin{array}{c}\mathsf{\alpha }\\ \mathrm{j}\end{array}\right)=\frac{\mathsf{\Gamma }(\mathsf{\alpha }+1)}{\mathsf{\Gamma }(\mathrm{j}+1)\mathsf{\Gamma }(\mathsf{\alpha }+1-\mathrm{j})}$

where $\mathsf{\Gamma }(·)$ is the gamma function, defined as follows:

(12)$\mathsf{\Gamma }(\mathsf{\alpha }+1)={{\displaystyle \int }}_0^{+\infty }{\mathrm{e}}^{-\mathrm{t}}{\mathrm{t}}^{\mathrm{z}-1}\mathrm{dt},\mathrm{Re}\left(\mathrm{z}\right)>0$

and satisfies the equation: $\mathsf{\Gamma }(\mathrm{z}+1)=\mathrm{z}\mathsf{\Gamma }\left(\mathrm{z}\right)$.

When fractional damping is added to the GMA system, the linear damping term describes the damping of the output rod and the fractional damping term describes the internal damping of the GMM bar, and the above two terms, together, constitute the damping term of the system. Equation (9) is transformed to obtain the transfer function of the system containing the fractional order calculus as follows:

(13)$\mathrm{G}\left(\mathrm{s}\right)=\frac{\mathrm{x}\left(\mathrm{s}\right)}{\mathrm{I}\left(\mathrm{s}\right)}=\frac{137,170}{{19\mathrm{s}}^2+2.48{\times 10}^5{(\mathrm{s}-\mathrm{s}}^{\mathsf{\alpha }})+2.22{\times 10}^9}$

Figure 5 shows the step response curves and experimental curves of the GMA system at different fractional orders $\mathsf{\alpha }$. From the Figure, it can be seen that when the order $\mathsf{\alpha }$ gradually changes from 0.84 to 0.89, the overshoot and peak time of the system gradually approach to the experimental value. Compared with the theoretical curve of the integer-order GMA system, the step response curve of the fractional-order GMA system gradually tends to the experimental curve with a higher degree of overlap. Therefore, adding the fractional order calculus to describe the damping characteristics of the GMA system can more accurately characterize the physical change process of the GMA system, which is more consistent with the actual situation.

To further investigate the effect of fractional damping on the GMA system, fractional damping is added to the simplified model, as shown in Figure 6.

After substituting Equation (10) into Equation (8) and normalizing it, the differential equation for the dynamics of the GMA system with fractional damping terms is obtained as follows:

(14)$\ddot{\mathrm{x}}{+\mathsf{\omega }}_0^2\mathrm{x}=\mathrm{fcos}\mathsf{\omega }\mathrm{t}-2\mathsf{\xi }\dot{\mathrm{x}} -\mathsf{\delta }\frac{{\mathrm{d}}^{\mathsf{\alpha }}}{{\mathrm{dt}}^{\mathsf{\alpha }}}{\mathrm{x}-\mathsf{\beta }\mathrm{x}}^2{-\mathsf{\gamma }\mathrm{x}}^3$

where ${\mathsf{\omega }}_0=\sqrt{\frac{\mathrm{k}}{\mathrm{m}}}$, $\mathrm{f}=\frac{\mathrm{F}}{\mathrm{m}}$, $2\mathsf{\xi }=\frac{\mathrm{c}}{\mathrm{m}}$, $\mathsf{\delta }=\frac{{\mathsf{\delta }}_1}{\mathrm{m}}$, $\mathsf{\beta }=\frac{{\mathsf{\beta }}_1}{\mathrm{m}}$, $\mathsf{\gamma }=\frac{{\mathsf{\gamma }}_1}{\mathrm{m}}$.

When the external excitation frequency of the GMA system is close to the intrinsic frequency of the system, the system will undergo resonant response. To facilitate the micro-amplitude vibration analysis of the nonlinear GMA system, adding the weak nonlinear factor $\mathsf{\epsilon }$ in the system before the nonlinear and damping terms and introducing the tuning parameter $\mathsf{\sigma }$, and ${\mathsf{\omega }}^2{=\mathsf{\omega }}_0^2+\mathsf{\epsilon }\mathsf{\sigma }$ is substituted into Equation (14) and organized as follows:

(15)$\ddot{\mathrm{x}}{+\mathsf{\omega }}^2\mathrm{x}=\mathsf{\epsilon }(\mathrm{fcos}\mathsf{\omega }\mathrm{t}-\mathsf{\sigma }\mathrm{x}-2\mathsf{\xi }\dot{\mathrm{x}} -\mathsf{\delta }\frac{{\mathrm{d}}^{\mathsf{\alpha }}}{{\mathrm{dt}}^{\mathsf{\alpha }}}{\mathrm{x}-\mathsf{\beta }\mathrm{x}}^2{-\mathsf{\gamma }\mathrm{x}}^3)$

The solution of Equation (14) can be set as:

(16)$\mathrm{x}=a\cos \mathsf{\varphi },\dot{\mathrm{x}}=-a\mathsf{\omega }\sin \mathsf{\varphi }$

where a is the amplitude. $\mathsf{\varphi }$ ($\mathsf{\varphi }=\mathsf{\omega }\mathrm{t}+\mathsf{\theta }$) is the phase, which is a slow-varying function of time.

Differentiating x and $\dot{\mathrm{x}}$ in Equation (16) with respect to t yields the following:

(17)$\left\{\begin{array}{c}\dot{\mathrm{x}}=-\dot{a}\cos \mathsf{\varphi }-(\dot{\mathsf{\theta }}+\mathsf{\omega })a\sin \mathsf{\varphi }\\ \ddot{\mathrm{x}}=-\dot{a}\mathsf{\omega }\sin \mathsf{\varphi }-(\dot{\mathsf{\theta }}+\mathsf{\omega })a\mathsf{\omega }\cos \mathsf{\varphi }\end{array}\right.$

Combining Equations (14) and (17) yields the following:

(18)$\dot{a}\sin \mathsf{\varphi }-a\dot{\mathsf{\theta }}\cos \mathsf{\varphi }=-\frac{1}{\mathsf{\omega }}{(\mathrm{E}}_1{+\mathrm{E}}_2)$

where ${\mathrm{E}}_1=\mathsf{\epsilon }\left[\mathrm{fcos}\right(\mathsf{\varphi }-\mathsf{\theta })+\mathsf{\sigma }\mathrm{acos}\mathsf{\varphi }+2\mathsf{\xi }\mathrm{a}\mathsf{\omega }\sin \mathsf{\varphi }{-\mathsf{\beta }\mathrm{a}}^2{\cos }^2\mathsf{\varphi }{-\mathsf{\gamma }\mathrm{a}}^3{\cos }^3\mathsf{\varphi }]$, ${\mathrm{E}}_2=-\mathsf{\epsilon }\mathsf{\delta }\frac{{\mathrm{d}}^{\mathsf{\alpha }}}{{\mathrm{dt}}^{\mathsf{\alpha }}}(\mathrm{acos}\mathsf{\varphi })$.

The averaging method [29] was used over the time interval [0, T] as follows:

(19)$\left\{\begin{array}{c}\dot{a}=-\frac{\mathrm{f}}{2\mathsf{\omega }}\sin \mathsf{\theta }-a\mathrm{A}\left(\mathsf{\alpha }\right)\\ a\dot{\mathsf{\theta }}=-\frac{\mathrm{f}}{2\mathsf{\omega }}\cos \mathsf{\theta }+\frac{{3\mathsf{\gamma }\mathrm{a}}^3}{8\mathsf{\omega }} -\frac{\mathsf{\omega }a}{2}+\frac{a}{2\mathsf{\omega }}\mathrm{B}\left(\mathsf{\alpha }\right)\end{array}\right.$

where $\mathrm{A}\left(\mathsf{\alpha }\right)=\mathsf{\xi }\frac{{\mathsf{\delta }\mathsf{\omega }}^{\mathsf{\alpha }-1}}{2}\sin (\frac{\mathsf{\alpha }\mathsf{\pi }}{2})$, ${\mathrm{B}\left(\mathsf{\alpha }\right)=\mathsf{\omega }}_0^2{+\mathsf{\xi }\mathsf{\delta }\mathsf{\omega }}^{\mathsf{\alpha }}\cos (\frac{\mathsf{\alpha }\mathsf{\pi }}{2})$.

To obtain the steady-state solution of the GMA system, let $\dot{a}=a\dot{\mathsf{\theta }}=0$ when $\mathsf{\theta }$ is eliminated, the amplitude-frequency response equation of the nonlinear GMA system is obtained as follows:

(20)${[\mathrm{A}}^2\left(\mathsf{\alpha }\right)+(\frac{{3\mathsf{\gamma }\mathrm{a}}^2}{8\mathsf{\omega }}-\frac{\mathsf{\omega }}{2}+\frac{{\mathrm{B}}^2\left(\mathsf{\alpha }\right)}{2\mathsf{\omega }}{)}^2{]a}^2=(\frac{\mathrm{f}}{2\mathsf{\omega }}{)}^2$

4. Amplitude and Frequency Characteristics of the GMA System

To solve the numerical solution of the GMA system, Equation (14) is deformed in reduced order using the superposition property of fractional order calculus as follows:

(21)$\left\{\begin{array}{c}\frac{{\mathrm{d}}^{\mathsf{\alpha }}\mathrm{x}\left(\mathrm{t}\right)}{{\mathrm{dt}}^{\mathsf{\alpha }}}=\mathrm{y}\left(\mathrm{t}\right)\\ \frac{{\mathrm{d}}^{1-\mathsf{\alpha }}\mathrm{y}\left(\mathrm{t}\right)}{{\mathrm{dt}}^{1-\mathsf{\alpha }}}=\mathrm{z}\left(\mathrm{t}\right)\\ \frac{\mathrm{dz}\left(\mathrm{t}\right)}{\mathrm{dt}}{=\mathrm{fcos}\mathsf{\omega }\mathrm{t}-\mathsf{\sigma }\mathrm{y}\left(\mathrm{t}\right)-2\mathsf{\xi }\mathrm{z}\left(\mathrm{t}\right)-\mathsf{\omega }}_0^2{\mathrm{x}\left(\mathrm{t}\right)-\mathsf{\beta }\mathrm{x}}^2{\left(\mathrm{t}\right)-\mathsf{\gamma }\mathrm{x}}^3\left(\mathrm{t}\right)\end{array}\right.$

The discretization of Equation (21) using the algorithm defined by the Grunwald–Letnikov fractional order calculus results in the following [30].

(22)$\left\{\begin{array}{c}{\mathrm{x}(\mathrm{t}}_{\mathrm{n}})=-{\displaystyle {\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}-1}}{\mathrm{w}}_{\mathrm{j}}^{\left(\mathsf{\alpha }\right)}{\mathrm{x}(\mathrm{t}}_{\mathrm{n}-\mathrm{j}})+∆{\mathrm{t}}^{\mathsf{\alpha }}{\mathrm{y}(\mathrm{t}}_{\mathrm{n}-1})\\ {\mathrm{y}(\mathrm{t}}_{\mathrm{n}})=-{\displaystyle {\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}-1}}{\mathrm{w}}_{\mathrm{j}}^{(1-\mathsf{\alpha })}{\mathrm{y}(\mathrm{t}}_{\mathrm{n}-\mathrm{j}})+∆{\mathrm{t}}^{1-\mathsf{\alpha }}z{(\mathrm{t}}_{\mathrm{n}-1})\\ {\mathrm{z}(\mathrm{t}}_{\mathrm{n}})=-{\displaystyle {\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}-1}}\begin{array}{c}{\mathrm{w}}_{\mathrm{j}}^{\left(1\right)}{\mathrm{x}(\mathrm{t}}_{\mathrm{n}-\mathrm{j}}{)+[\mathrm{fcos}\mathsf{\omega }\mathrm{t}}_{\mathrm{n}-1}{-\mathsf{\omega }}_0^2{\mathrm{x}(\mathrm{t}}_{\mathrm{n}-1}{)-\mathsf{\beta }\mathrm{x}}^2{(\mathrm{t}}_{\mathrm{n}-1}{)-\mathsf{\gamma }\mathrm{x}}^3{(\mathrm{t}}_{\mathrm{n}-1})\\ {-2\mathsf{\xi }\mathrm{z}(\mathrm{t}}_{\mathrm{n}-1}{)-\mathsf{\delta }\mathrm{y}(\mathrm{t}}_{\mathrm{n}-1}\left)\right]\end{array}∆\mathrm{t}\end{array}\right.$

where t_n = n∆t, and n is a large enough positive integer, which generally means that the system runs for n cycles after the transient response disappears. ${\mathrm{w}}_{\mathrm{j}}^{\left(\mathsf{\alpha }\right)}$ is the quadratic term coefficient and satisfies the following iterative relation:

(23)${\mathrm{w}}_0^{\left(\mathsf{\alpha }\right)}{=1,\mathrm{w}}_{\mathrm{j}}^{\left(\mathsf{\alpha }\right)}=(1-\frac{1+\mathsf{\alpha }}{\mathrm{j}}{)\mathrm{w}}_{\mathrm{j}-1}^{\left(\mathsf{\alpha }\right)}$

The basic parameter of the GMA system is taken as ${\mathsf{\omega }}_0^2=0.4$, $\mathsf{\gamma }=0.1$, $\mathsf{\beta }=-0.4$, $2\mathsf{\xi }=0.2$, $\mathsf{\delta }=0.5$, and $\mathrm{f}=0.1$. During the calculation, the initial conditions x, y and z take the values [0, 0, 0.5], the calculation time is 500 s, the time step $∆\mathrm{t}$ = 0.01, and the first 400 s are rounded off to ensure stability of the numerical values. Figure 7 shows the relationship between the fractional order $\mathsf{\alpha }$ and the external excitation frequency $\mathsf{\omega }$ and the response amplitude a of the system. From the Figure, it can be seen that the relationship between the order $\mathsf{\alpha }$ and the maximum value of the amplitude a is not a monotonic one in the interval $\mathsf{\alpha }\in (0,2)$, but the amplitude a tends to decrease and then increase as the order $\mathsf{\alpha }$ increases. This Figure and Equation (19) both show that the system response amplitude a is not only related to the external excitation amplitude f and frequency $\mathsf{\omega }$, but also related to the fractional order $\mathsf{\alpha }$.

Figure 8 shows the main resonance amplitude-frequency characteristic curves of the GMA system at different fractional orders $\mathsf{\alpha }$ for the above values of the basic parameters. From the Figure, it can be seen that the nonlinear curve bends to the right and is hard nonlinear. When $0<\mathsf{\alpha }<1$ the order $\mathsf{\alpha }$ is 0.4, 0.5 and 0.6, respectively, with the increase of order $\mathsf{\alpha }$, the resonant frequency and the main resonant peak of the system tend to decrease, and the hysteresis region also decreases, because the fractional damping coefficient $\mathsf{\delta }$ in Equation (11) is less than 0, which also explains the phenomenon that the overshoot of the system increases with the increase of the fractional order in Figure 5. When $1<\mathsf{\alpha }<2$, the order $\mathsf{\alpha }$ is 1, 1.2 and 1.4, respectively, with the increase of order $\mathsf{\alpha }$, the main resonance peak of the system tends to increase; the left side of the resonance region of the curve is a single-value solution, and the right side of the resonance bifurcation appears, the left side of the resonance curve is a stable solution, and the right side of the resonance curve is an unstable solution, and the resonance curve “jumps” phenomenon occurs. This also shows the trend of decreasing and then increasing with the order $\mathsf{\alpha }$ when the amplitude a is at the resonance frequency in Figure 7.

Figure 9 shows the resonance amplitude and frequency characteristics of the GMA system with different external excitation amplitude f when the basic parameter of the GMA system is ${\mathsf{\omega }}_0^2=0.4$, $\mathsf{\gamma }=0.1$, $\mathsf{\beta }=-0.4$, $2\mathsf{\xi }=0.2$, $\mathsf{\delta }=0.5$, $\mathsf{\alpha }=0.88$. As can be seen from the Figure, the resonance amplitude of the GMA system shows a trend of increasing and then decreasing with the increase of external excitation frequency $\mathsf{\omega }$. Along with the increase of external excitation amplitude, the resonance peak of the system also becomes larger accordingly. The resonance curve shows the phenomenon of outward expansion, the skeleton curve is deflected to the right, and it is hard nonlinear, and the resonance region and instability region are gradually widened.

Figure 10 shows the main resonance amplitude-frequency characteristic curves of the GMA system with different fractional damping coefficients $\mathsf{\delta }$ when the basic parameters of the GMA system are taken as ${\mathsf{\omega }}_0^2=0.4$, $\mathsf{\gamma }=0.1$, $\mathsf{\beta }=-0.4$, $2\mathsf{\xi }=0.2$, $\mathrm{f}=0.1$, and $\mathsf{\alpha }=0.88$. As can be seen from the Figure, along with the absolute value of the damping coefficient $\mathsf{\delta }$ becomes larger, the resonance peak of the GMA system becomes smaller accordingly, the resonance curve inwardly shrinks, and its resonance region and instability region gradually become narrower.

5. Bifurcation and Chaotic Behavior of GMA Fractional Damping Systems

In order to further study the instability of GMA system dynamics, this paper adds the fractional order calculus to describe the damping characteristics of the system, and studies the bifurcation and chaos characteristics of the GMA system containing fractional damping terms.

5.1. Chaotic Properties of the System When the Fractional Order $\alpha$ Changes

In order to further reveal the mechanism and dynamics of the effect of fractional damping on the motion of the nonlinear GMA system, the system parameters in Equation (14) are taken as ${\mathsf{\omega }}_0^2=-0.3$, $\mathsf{\omega }=1.1$, $\mathsf{\gamma }=0.5$, $\mathsf{\beta }=-0.4$, $2\mathsf{\xi }=0.3$, $\mathsf{\delta }=0.2$, and $\mathrm{f}=0.4$. The numerical solution of Equation (14) is solved using the algorithm defined by the Grunwald–Letnikov fractional order derivative, and its numerical solution is simulated by simulation, and the steady-state response part is selected for the analysis of the GMA system. The fractional order $\mathsf{\alpha }$ is varied from 0 to 2, the step size is taken as 0.01, and the initial condition is [0, 0, 0].

Figure 11 shows the bifurcation diagram and the Lyapunov exponent of the GMA system with the fractional order as the control parameter. The bifurcation and chaos occur in the GMA system as the order changes, and the path to chaos in the GMA system when the order $\mathsf{\alpha }$ changes is as follows: period 3 → chaotic motion → period 4 → chaotic motion → period 4 → period 2 → period 2 → chaotic motion → period 2 → period 4 → chaotic motion.

From Figure 11, the GMA system enters into chaotic motion from period 3 when in the $\mathsf{\alpha }\in (0,0.04)$ interval, at which time the Lyapunov exponent in Figure 11b changes from negative to positive. When in the $\mathsf{\alpha }\in (0.04,0.22)$ interval, the GMA system degenerates from chaotic motion to period 4. When in the $\mathsf{\alpha }\in (0.22,0.61)$ interval, the GMA system enters from period 4 to chaotic motion. As shown in Figure 12 ($\mathsf{\alpha }$ = 0.41), the time domain waveform of the system shows that the waveform is a non-periodic irregular reciprocating motion; the phase trajectory curve fills the phase diagram space area, and the trajectory curve is an irregular multi-turn coil and cannot be closed for a long time; there are multiple frequency peaks in the amplitude spectral response that cannot be metric; the Poincaré diagram shows a “horseshoe” attractor, and the GMA system can be judged to be in a state of chaotic motion.

When in the $\mathsf{\alpha }\in (0.61,0.65)$ interval, the GMA system degenerates into a period 4 motion via chaotic motion, as shown in Figure 13 ($\mathsf{\alpha }$ = 0.63), and the system waveform diagram behaves as a superposition combination of four waveform signals; the system phase trajectory curve is a closed coil of four multi-turn phase sets; the spectral components in the amplitude spectrogram are mainly concentrated in the four frequency regions of 0.28, 0.57, 1.1 and 1.69; and four discrete attractors appear in the Poincaré diagram, when the Lyapunov exponent corresponding to Figure 11b also changes from positive to negative in this local range, so that the system is in a period of 4 motion.

When in the $\mathsf{\alpha }\in (0.65,0.94)$ interval, the GMA system bifurcates into the period 2 motion via the period 4 inverse times period, as shown in Figure 14 ($\mathsf{\alpha }$ = 0.83), the time domain waveform of the GMA system is an iterative combination of two period waveform signals; the closed-loop smooth curve in the phase trajectory diagram forms two mutually nested ellipses; the spectral components in the amplitude spectrum diagram are relatively single, mainly concentrated in the frequency region of 0.57 and 1.1; the Poincaré diagram is presented as two mutually isolated attractors, at which time the Lyapunov exponent corresponding to Figure 11b is less than 0, and thus, it can be determined that the GMA system is in period 2 motion at this point.

When in the $\mathsf{\alpha }\in (0.94,1.7)$ interval, the GMA system degenerates into stable periodic motion by a period 2 motion inverse times period bifurcation, as shown in Figure 15 ($\mathsf{\alpha }$ = 0.97), the time domain waveform of the GMA system shows a single and regular reciprocal motion; the closed-loop smooth curve in the phase trajectory forms an ellipse; the frequency component in the amplitude spectrum is relatively single, mainly concentrated in the frequency region of 1.1; the Poincaré diagram shows a single isolated attractor, which corresponds to the Lyapunov index greater than 0 in Figure 11b, and thus, it can be concluded that the GMA system has a typical stable periodic motion at this point.

When in the $\mathsf{\alpha }\in (1.7,1.87)$ interval, the GMA system enters period 2 motion through the period motion times period bifurcation. When in the $\mathsf{\alpha }\in (1.87,2)$ interval, the GMA system enters the chaotic motion after bifurcating into the period 4 motion through the period 2 motion times period, which corresponds to the change of Lyapunov exponent from negative to positive in Figure 11b.

5.2. Chaotic Characteristics of the GMA System When the Damping Coefficient $\delta$ Varies at Different Fractional Order $\alpha$

To further study the influence law of fractional order $\mathsf{\alpha }$ on the GMA system, the fractional damping coefficient $\mathsf{\delta }$ is taken as the control parameter, and different orders $\mathsf{\alpha }$, respectively, are selected to study the dynamic behavior of the GMA system. When $\mathsf{\alpha }$ = 1 is the traditional integer order GMA system, in order to facilitate the observation of the relationship between the global motion form of the system and the damping coefficient $\mathsf{\delta }\in (-0.2,0.2)$, the damping coefficient $\mathsf{\delta }$ is selected, the step size 0.005 is chosen, and the initial condition is [0, 0, 0], and the bifurcation diagram and Lyapunov exponent diagram of the GMA system are obtained, as shown in Figure 16.

From Figure 16, it can be seen that the GMA system appears as bifurcation and chaos with the change of fractional damping coefficient $\mathsf{\delta }$. The path of the GMA system to chaos when the damping coefficient $\mathsf{\delta }$ changes is: paroxysmal chaotic motion → period 2 → paroxysmal chaotic motion → period motion → jump motion → period 2 → period motion.

When in the $\mathsf{\delta }\in (-0.2,-0.12)$ interval, the system degenerates from paroxysmal chaotic motion into period 2 motion, at which time the Lyapunov exponent corresponding to Figure 16b changes repeatedly from positive to negative up and down to negative values, and the side-response system transforms from chaotic motion to periodic motion. As shown in Figure 17 ($\mathsf{\delta }$ = −0.18), the waveform of the system shows a reciprocal motion with neither regularity nor periodicity; its phase trajectory curve is a long-term unclosable, irregular multi-turned coil; a continuous wider amplitude spectrum appears in the amplitude spectrum; the Poincaré mapping is neither a finite point nor a closed curve, and the Lyapunov index corresponding to Figure 16b is positive at this time. All the above indicators show that the GMA system is in a typical chaotic motion at this time. When in the interval $\mathsf{\delta }\in (-0.12,-0.01)$, the system performs period 2 motion, and the Lyapunov exponent in Figure 16b is less than 0. When in the interval $\mathsf{\delta }\in (-0.01,0.055)$, the system enters from the period 2 motion into the paroxysmal chaotic motion. When in the interval $\mathsf{\delta }\in (0.055,0.085)$, the system degenerates from paroxysmal chaotic motion to periodic motion, and the Lyapunov exponent in Figure 16b changes from positive to negative at this time. At $\mathsf{\delta }$ = 0.09, after jumping motion times, the period bifurcation changes into the period 2 motion. When in the $\mathsf{\delta }\in (0.09,0.2)$ interval, the system degenerates from the period 2 motion inverted times period bifurcation to the periodic motion, and then the GMA system enters into the periodic motion.

To facilitate the distinction between the details of the bifurcation diagrams of the GMA system caused by the fractional damping coefficient $\mathsf{\delta }$ at different fractional orders $\mathsf{\alpha }$. Figure 18 shows the local enlarged plots in the interval of damping coefficient $\mathsf{\delta }\in (-0.3,-0.15)$ when the order $\mathsf{\alpha }$ is taken as 0.55, 0.75, 1, 1.25 and 1.45, respectively. It is easy to see from the plots that the GMA system has similar dynamical behaviors at different orders $\mathsf{\alpha }$ but the regions where chaotic motions are generated are different. At $0<\mathsf{\alpha }<1$, the area of the chaotic region in the system decreases significantly with the increase of order $\mathsf{\alpha }$ and enters the chaotic region later; at $1<\mathsf{\alpha }<2$, the area of the chaotic region in the system increases significantly with the increase of order $\mathsf{\alpha }$ and enters the chaotic region earlier.

5.3. Chaotic Characteristics of the GMA System When the Excitation Amplitude f Varies at Different Fractional Orders $\alpha$

Similarly, the external excitation amplitude f is taken as the control parameter, and different orders $\mathsf{\alpha }$ are chosen to study the dynamic behavior of the GMA system, respectively. To facilitate the observation of the relationship between the global motion form of the system and the amplitude f, the conventional GMA system at $\mathsf{\alpha }$ = 1 is still chosen first, and the external excitation amplitude $\mathrm{f}\in (0,1.5)$ is chosen with a step size of 0.01 and initial conditions of [0, 0, 0], and the bifurcation diagram and Lyapunov exponent diagram of the GMA system are obtained, as shown in Figure 19.

From Figure 19, it can be seen that the GMA system bifurcates and chaos occurs with the change of external excitation amplitude f. The path to chaos when the amplitude f changes is as follows: periodic motion → period 2 → paroxysmal chaotic motion → period 2 → period 4 → chaotic motion → proposed periodic motion → chaotic motion → period 3 → chaotic motion → period 4 → period 2 → periodic motion → jump motion → periodic motion.

When in the $\mathrm{f}\in (0,0.4)$ interval, the system changes from a period motion times period bifurcation to period 2 motion, when the Lyapunov exponent in Figure 19b corresponds to less than zero. When in the $\mathrm{f}\in (0.4,0.76)$ interval, the system transforms from period 2 motion to paroxysmal chaotic motion, when the Lyapunov exponent in Figure 19b changes from negative to repeatedly changing above and below the positive and negative values, and the side reaction system transforms from periodic motion to paroxysmal chaotic motion. In the interval $\mathrm{f}\in (0.76,0.84)$, the system degenerates from chaotic motion to period 2 motion. In the $\mathrm{f}\in (0.84,1.07)$ interval, the system changes from period 2 times bifurcation to period 4 motion, and then to chaotic motion, and when f = 0.97, the system is in the proposed periodic motion, at which time the Lyapunov exponent in Figure 19b appears as a spike. At the interval of $\mathrm{f}\in (1.07,1.11)$, the system changes from chaotic motion to period 3 motion. As shown in Figure 20 (f = 1.1), the time domain waveform of the GMA system is an iterative combination of three periodic waveform signals; the closed-loop smooth curve in the phase trajectory diagram forms three mutually nested ellipses; the spectral components in the amplitude spectrum diagram are relatively single, and mainly concentrated in three frequency regions; the Poincaré mapping diagram is presented as three mutually isolated attractors, at which time the corresponding Figure 19b Lyapunov exponent is negative, thus, it can be determined that the GMA system is in period 3 motion at this point. In the interval $\mathrm{f}\in (1.11,1.27)$, the system changes from period 3 to chaotic motion. In the $\mathrm{f}\in (1.27,1.43)$ interval, the system degenerates from chaotic motion to period 4, then to period 2 motion through the inverse period bifurcation, and then to period motion. When in the $\mathrm{f}\in (1.43,1.5)$ interval, the system goes to the periodic motion after the jumping motion.

To facilitate the distinction between the details of the bifurcation diagrams of the GMA system caused by the external excitation amplitude f at different fractional orders $\mathsf{\alpha }$. Figure 21 shows the local enlarged plots in the interval of amplitude $\mathrm{f}\in (0,0.7)$ when $\mathsf{\alpha }$ is taken as 0.5, 0.8, 1 and 1.3, respectively. As shown in the Figure, the GMA system produces different regions of chaotic motion at different orders $\mathsf{\alpha }$. As the order $\mathsf{\alpha }$ increases, the area of the chaotic region appears significantly less in the system and enters the chaotic region later.

From the above analysis, it can be seen that with the change of external excitation amplitude f, fractional damping coefficient $\mathsf{\delta }$ and fractional order $\mathsf{\alpha }$, the output response of the GMA system will be accompanied by complex forms of motion such as periodic, apropos-periodic and chaotic motion. However, the proposed periodic and chaotic motions are unstable and difficult to predict and control accurately. As a result, uncontrollable and unpredictable chaotic phenomena can occur during the operation of the GMA system due to the changes of the excitation force and the fractional order damping coefficient.

6. Conclusions

In this paper, according to the working principle of the GMA system, the fractional order calculus is invoked to describe the damping characteristics of the system, the dynamics equations of the GMA system containing fractional damping are established, and the nonlinear vibration response of the GMA system is studied; the numerical solution of the system is solved by using the algorithm defined by Grunwald–Letnikov fractional order derivatives, and the effects of fractional order $\mathsf{\alpha }$ as well as the external excitation amplitude f and fractional damping coefficient $\mathsf{\delta }$ on the bifurcation and chaotic characteristics of the GMA system at different orders are studied from qualitative and quantitative perspectives through MATLAB simulations, and the results are shown as follows:

• (1). Compared with the integer-order GMA system, the fractional-order GMA system can more accurately express the physical change process of the system and is closer to the actual situation.

• (2). When $0<\mathsf{\alpha }<1$, with the increase of fractional order $\mathsf{\alpha }$, the resonance frequency and the main resonance peak of the system tend to decrease, and the lag region also decreases; when $1<\mathsf{\alpha }<2$, the main resonance peak of the system tends to increase. With the increase of excitation frequency, the system amplitude tends to increase and then decrease, along with the increase of excitation amplitude, the resonance curve shows the phenomenon of outward expansion, the skeleton curve deflects to the right, the hard nonlinearity, the resonance region and the instability region gradually widen, and the hysteresis phenomenon appears. With the increase of absolute value of damping coefficient $\mathsf{\delta }$, the resonance peak of the system decreases accordingly, the resonance curve inwardly shrinks, and the resonance region and the instability region gradually become narrower.

• (3). The order $\mathsf{\alpha }$ has a more significant effect on the dynamical behavior of the GMA system. When the order $\mathsf{\alpha }$ is less than 0.61 and greater than 1.88, complex dynamical behaviors such as paroxysmal chaos and chaotic motion appear in the GMA system.

• (4). At different orders $\mathsf{\alpha }$, chaos exists in the system, as the external excitation amplitude f and damping coefficient $\mathsf{\delta }$ vary; at $0<\mathsf{\alpha }<1$, the fractional order GMA system enters into chaotic motion earlier than the integer order GMA system, showing richer and different instability conditions of the dynamical system.

Author Contributions

Conceptualization, Q.M.; methodology, Q.M.; software, X.F.; formal analysis, Q.M.; investigation, Q.M.; resources, J.W. and H.Y.; data curation, J.Y.; writing—original draft preparation, Q.M.; writing—review and editing, Q.M.; supervision, H.Y.; project administration, H.Y.; funding acquisition, H.Y. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Some or all of the data, models, or code generated or used during the study are available from the corresponding author on reasonable request.

Acknowledgments

The authors gratefully acknowledge the support of the Inner Mongolia Natural Science Foundation Project (2020LH05023).

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. GMA structure: 1. output rod; 2. disc spring; 3. drive coil; 4. GMM bar; 5. permanent magnet; 6. lower guide ring; 7. base; 8. upper end cap; 9. upper permeable ring; 10. permeable sleeve; 11. outer wall sleeve; 12. permeable sheet; 13. coil.

Figure 2. Simplified model of the GMA system.

Figure 3. Equivalent mechanical model of the GMA output rod.

Figure 4. Experimental and theoretical curves of step response of the GMA.

Figure 5. Experimental step response curves of GMA and step response curves of different fractional orders.

Figure 6. Simplified model of GMA fractional damping system.

Figure 7. The three-dimensional curve of the amplitude and frequency characteristics of the GMA system.

Figure 8. Amplitude and frequency characteristic curves of the system under different [Forumla omitted. See PDF.]: (a) [Forumla omitted. See PDF.]; (b) [Forumla omitted. See PDF.].

Figure 9. Amplitude and frequency characteristic curves of the system at different f.

Figure 10. Amplitude and frequency characteristic curves of the system under different [Forumla omitted. See PDF.]: (a) [Forumla omitted. See PDF.]; (b) [Forumla omitted. See PDF.].

Figure 11. Bifurcation diagram and Lyapunov exponent diagram of the system when [Forumla omitted. See PDF.] varies: (a) Bifurcation diagram; (b) Lyapunov exponent diagram.

Figure 12. Nonlinear response of the GMA system at [Forumla omitted. See PDF.] = 0.41: (a) Time domain waveform diagram; (b) Phase trajectory curve; (c) Amplitude spectrum; (d) Poincaré diagram.

Figure 13. Nonlinear response of the GMA system at [Forumla omitted. See PDF.] = 0.63: (a) Time domain waveform diagram; (b) Phase trajectory curve; (c) Amplitude spectrum; (d) Poincaré diagram.

Figure 14. Nonlinear response of the GMA system at [Forumla omitted. See PDF.] = 0.83: (a) Time domain waveform diagram; (b) Phase trajectory curve; (c) Amplitude spectrum; (d) Poincaré diagram.

Figure 14. Nonlinear response of the GMA system at [Forumla omitted. See PDF.] = 0.83: (a) Time domain waveform diagram; (b) Phase trajectory curve; (c) Amplitude spectrum; (d) Poincaré diagram.

Figure 15. Nonlinear response of the GMA system at [Forumla omitted. See PDF.] = 0.97: (a) Time domain waveform diagram; (b) Phase trajectory curve; (c) Amplitude spectrum; (d) Poincaré diagram.

Figure 16. Bifurcation diagram and Lyapunov exponent diagram of the system with [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] = 1: (a) Bifurcation diagram; (b) Lyapunov exponent diagram.

Figure 17. Nonlinear response of the GMA system at [Forumla omitted. See PDF.] = −0.18: (a) Time domain waveform diagram; (b) Phase trajectory curve; (c) Amplitude spectrum; (d) Poincaré diagram.

Figure 18. Local bifurcation diagram of the system with [Forumla omitted. See PDF.] at different orders [Forumla omitted. See PDF.]: (a) [Forumla omitted. See PDF.] = 0.55; (b) [Forumla omitted. See PDF.] = 0.75; (c) [Forumla omitted. See PDF.] = 1; (d) [Forumla omitted. See PDF.] = 1.25; (e) [Forumla omitted. See PDF.] = 1.45.

Figure 19. Bifurcation diagram and Lyapunov exponent diagram of the system with f at [Forumla omitted. See PDF.] = 1: (a) Bifurcation diagram; (b) Lyapunov exponent diagram.

Figure 20. Nonlinear response of the GMA system at f = 1.1: (a) Time domain waveform diagram; (b) Phase trajectory curve; (c) Amplitude spectrum; (d) Poincaré diagram.

Figure 21. Local bifurcation diagram of the system with f at different orders [Forumla omitted. See PDF.]: (a) [Forumla omitted. See PDF.] = 0.5; (b) [Forumla omitted. See PDF.] = 0.8; (c) [Forumla omitted. See PDF.] = 1; (d) [Forumla omitted. See PDF.] = 1.3.

Figure 21. Local bifurcation diagram of the system with f at different orders [Forumla omitted. See PDF.]: (a) [Forumla omitted. See PDF.] = 0.5; (b) [Forumla omitted. See PDF.] = 0.8; (c) [Forumla omitted. See PDF.] = 1; (d) [Forumla omitted. See PDF.] = 1.3.

References

1. Clark, A.E.; Belson, H.S. Giant Room-Temperature Magnetostriction in TbFe2 and DyFe2. Phys. Rev. B; 1972; 5, pp. 3642-3644. [DOI: https://dx.doi.org/10.1103/PhysRevB.5.3642]

2. Liu, X.H.; Gao, L.; Wu, Y.; Song, H.R.; He, Y. Micro-displacement amplifier of giant magnetostrictive actuator using flexure hinges. J. Magn. Magn. Mater.; 2022; 556, 169415. [DOI: https://dx.doi.org/10.2139/ssrn.4036635]

3. Stachowiak, D.; Demenko, A. Finite Element and Experimental Analysis of an Axisymmetric Electromechanical Converter with a Magnetostrictive Rod. Energies; 2020; 13, 1230. [DOI: https://dx.doi.org/10.3390/en13051230]

4. Yang, Y.; Yang, B.; Niu, M. Adaptive trajectory tracking of magnetostrictive actuator based on preliminary hysteresis compensation and further adaptive filter controller. Nonlinear Dyn.; 2018; 92, pp. 1109-1118. [DOI: https://dx.doi.org/10.1007/s11071-018-4112-6]

5. Zhu, Y.C.; Li, Y.S. Development of a deflector-jet electrohydraulic servovalve using a giant magnetostrictive material. Smart Mater. Struct.; 2014; 23, 115001. [DOI: https://dx.doi.org/10.1088/0964-1726/23/11/115001]

6. Duo, T.; Li, Y.T. Finite Element Solutions for Magnetic Field Problems in Terfenol-D Transducers. Sensors; 2020; 20, 2808.

7. Fang, Z.W.; Zhang, Y.W.; Li, X.; Ding, H.; Chen, L.Q. Integration of a nonlinear energy sink and a giant magnetostrictive energy harvester. J. Sound Vib.; 2017; 391, pp. 35-49. [DOI: https://dx.doi.org/10.1016/j.jsv.2016.12.019]

8. Braghin, F.; Cinquemani, S.; Resta, F. A low frequency magnetostrictive inertial actuator for vibration control. Sens. Actuators A Phys.; 2012; 180, pp. 67-74. [DOI: https://dx.doi.org/10.1016/j.sna.2012.03.015]

9. Feng, X.U.; Zhang, H.; Jiang, C.B.; Hu, H.B. Designing and performance research of giant magnetostrictive actuator. Acta Aeronaut. Astronaut. Sin.; 2002; 23, pp. 552-555.

10. Yi, S.C.; Yang, B.T.; Meng, G. Ill-conditioned dynamic hysteresis compensation for a low-frequency magnetostrictive vibration shaker. Nonlinear Dyn.; 2019; 96, pp. 535-551. [DOI: https://dx.doi.org/10.1007/s11071-019-04804-1]

11. Zhu, W.; Bian, L.X.; Cheng, L.; Rui, X.T. Non-linear compensation and displacement control of the bias-rate-dependent hysteresis of a magnetostrictive actuator. Precis. Eng.; 2017; 50, pp. 107-113. [DOI: https://dx.doi.org/10.1016/j.precisioneng.2017.04.018]

12. Zeng, H.Q.; Zeng, G.X. Nonlinear Behaviors of Giant Magnetostrictive High Power Ultrasonic Transducer. Appl. Mech. Mater.; 2012; 128, pp. 918-922. [DOI: https://dx.doi.org/10.4028/www.scientific.net/AMM.128-129.918]

13. Gao, X.; Liu, Y. Research of giant magnetostrictive actuator’s nonlinear dynamic behaviours. Nonlinear Dyn.; 2018; 92, pp. 793-802. [DOI: https://dx.doi.org/10.1007/s11071-018-4061-0]

14. Yan, H.B.; Gao, H.; Yang, G.; He, H.; Niu, Y.; Liu, P. Bifurcation and chaos characteristics of hysteresis vibration system of giant magnetostrictive actuator. Chin. Phys. B; 2020; 29, pp. 194-205. [DOI: https://dx.doi.org/10.1088/1674-1056/ab65b4]

15. Sylvain, Z.N.; Victor, K.T.; Kisito, T.P. Hysteretic dynamics inducing coexistence of attractors in a thin magnetostrictive actuator system with quintic nonlinearity. J. Magn. Magn. Mater.; 2020; 507, pp. 858-878.

16. Dung, N.V.; Tho, N.C.; Ha, N.M.; Hieu, V.T. On the Finite Element Model of Rotating Functionally Graded Graphene Beams Resting on Elastic Foundation. Math. Probl. Eng.; 2021; 2021, 1586388.

17. Javani, M.; Kiani, Y.; Eslami, M.R. Free vibration analysis of FG-GPLRC L-shaped plates implementing GDQE approach. Thin Walled Struct.; 2021; 162, 107600. [DOI: https://dx.doi.org/10.1016/j.tws.2021.107600]

18. Nematollahi, M.S.; Mohammadi, H.; Dimitri, R.; Tornabene, F. Nonlinear vibration of functionally graded graphene nanoplatelets polymer nanocomposite sandwich beams. Appl. Sci.; 2020; 10, 5669. [DOI: https://dx.doi.org/10.3390/app10165669]

19. Bedon, C. Issues on the vibration analysis of in-service laminated glass structures: Analytical, experimental and numerical investigations on delaminated beams. Appl. Sci.; 2019; 9, 3928. [DOI: https://dx.doi.org/10.3390/app9183928]

20. Nam, V.H.; Vinh, P.V.; Chinh, N.V.; Thom, D.V.; Hong, T.T. A New Beam Model for Simulation of the Mechanical Behaviour of Variable Thickness Functionally Graded Material Beams Based on Modified First Order Shear Deformation Theory. Materials; 2019; 12, 404. [DOI: https://dx.doi.org/10.3390/ma12030404]

21. Gebrel, I.F.; Asokanthan, S.F. Influence of System and Actuator Nonlinearities on the Dynamics of Ring-Type MEMS Gyroscopes. Vibration; 2021; 4, pp. 805-821. [DOI: https://dx.doi.org/10.3390/vibration4040045]

22. Coccolo, M.; Seoane, J.M.; Lenci, S.; Sanjuán, A.F. Fractional damping effects on the transient dynamics of the Duffing oscillator. Commun. Nonlinear Sci. Numer. Simul.; 2022; 117, 106959. [DOI: https://dx.doi.org/10.1016/j.cnsns.2022.106959]

23. Li, Z.N.; Li, Y.L.; Wang, D.; Peng, Z.K.; Wang, H.F. Dynamic Characteristics of Rotor System with a Slant Crack Based on Fractional Damping. Chin. J. Mech. Eng.; 2021; 34, pp. 89-106. [DOI: https://dx.doi.org/10.1186/s10033-021-00543-w]

24. Zarraga, O.; Sarría, I.; García-Barruetabeña, J.; Cortés, F. An Analysis of the Dynamical Behaviour of Systems with Fractional Damping for Mechanical Engineering Applications. Symmetry; 2019; 11, 1499. [DOI: https://dx.doi.org/10.3390/sym11121499]

25. Sylvain, Z.N.; Victor, K.T.; Bruno, N.G.; Kisito, T.P. Extreme multistability in a fractional-order thin magnetostrictive actuator (TMA). SeMA J.; 2021; 78, pp. 347-365. [DOI: https://dx.doi.org/10.1007/s40324-020-00238-7]

26. Yang, Y.G.; Sun, Y.H.; Xu, W. Stochastic bifurcation analysis of a friction-damped system with impact and fractional derivative damping. Nonlinear Dyn.; 2021; 105, pp. 3131-3138. [DOI: https://dx.doi.org/10.1007/s11071-021-06806-4]

27. Rodríguez, C.; Barrio, A.; Orue, I.; Vilas, J.L.; León, L.M.; Barandiarán, J.M.; Ruiz, M.F. High magnetostriction polymer-bonded Terfenol-D composites. Sens. Actuators A Phys.; 2008; 142, pp. 538-541. [DOI: https://dx.doi.org/10.1016/j.sna.2007.05.021]

28. Wang, A.M. Research on Structural Design, System Modeling and Parameter Identification of Giant Magnetostrictive Exciter. Ph.D. Thesis; Lanzhou Jiaotong University: Lanzhou, China, 2020.

29. Shen, Y.; Yang, S.; Xing, H.; Ma, H. Primary resonance of Duffing oscillator with two kinds of fractional-order derivatives. Int. J. Non-Linear Mech.; 2012; 47, pp. 975-983. [DOI: https://dx.doi.org/10.1016/j.ijnonlinmec.2012.06.012]

30. Petrá, I. Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2011.