Nonlinear dynamics and static bifurcations control of the 12-pole magnetic bearings system utilizing the integral resonant control strategy

Abstract

In this study, the Integral Resonant Controller (IRC) is presented along with the Proportional-Derivative (PD) controller as a novel control technique to control the dynamical behaviors of the 12-pole rotor active magnetic bearing system. According to the proposed control strategy, the system model has been derived as a two-degree-of-freedom nonlinear dynamical system coupled with two first-order filters. The obtained mathematical model has been analyzed utilizing the asymptotic analysis. The nonlinear algebraic system that governs the steady-state vibration amplitudes and corresponding phase angles of the considered system has been extracted. The influence of the IRC control parameters on the rotor dynamics has been explored by plotting the different bifurcation diagrams. The analytical investigations demonstrated that the vibration suppression efficiency of the combined controller (i.e., PD and IRC) is proportional to the product of control and feedback gains of the IRC as well as the derivative gain of the PD-controller. In addition, it is found that the controller efficiency is inversely proportional to both the square of the internal loop feedback gain of the IRC and the position gain of the PD-controller. Accordingly, an objective function has been derived to design the best control gains of the proposed controller strategy. The analytical and numerical simulations confirmed that the suggested control method can suppress the system vibration and eliminate the catastrophic bifurcation behaviors if the control gains are selected according to the proposed objective function. Finally, the effect of failing one of the coupled integral resonant controllers on the rotor dynamics has been explored as a precautionary procedure. It is found that the failure of one of the coupled integral resonant controllers may breakdown the bifurcation symmetry of the rotor system, but the system does not lose its stability.

Introduction

The magnetic bearings system is one of the best suspension systems for rotating shafts. Its working principle is based on generating a balanced electromagnetic force to suspend the rotor in its hovering attitude without any contact with the stationary parts. The suspension of the shafts without physical contact provides the magnetic bearings system with many advantages over the other conventional suspension systems such as operation with high speed, frictionless operation, lubrication is not required, low failure rate, less maintenance, high reliability, and high durability. As a result, numerous research papers have been devoted to investigate the dynamical characteristics of various magnetic bearing configurations, as well as the many control methodologies that have been introduced. Ji et al. ¹ studied the four-pole configuration of the magnetic bearings system. They investigated the considered model as an autonomous and non-autonomous system. They found that the investigated model shows a variety of nonlinear behaviors, including Hopf, saddle-connection, and saddle-node bifurcations. The six-pole configuration of the magnetic bearing system has been investigated with two different control methodologies by Saeed et al.^2,3 They presented two control techniques based on the Cartesian displacement and velocity of rotor lateral oscillations. ² In addition, they introduced a proportional-derivative control algorithm based on the polar displacement and polar velocity of the rotating shaft. According to the studied control methods, they reported that the Cartesian-based controller has high efficiency than the polar-based one in reducing the system oscillations. However, the polar-based control method is more stable than the Cartesian-based one. The eight-pole configuration has received many analyses and control techniques.^4–14 Ji and Hansen ⁴ used a traditional PD-controller to reduce the undesired vibrations of an eight-pole system at the primary resonance condition. The suggested configuration model has been derived as a nonlinear two-degree-of-freedom dynamical system. The asymptotic analysis has been applied to explore the static bifurcation behaviors of the rotor system. They concluded that the system may exhibit simultaneous bi-stable, or tri-stable solutions depending on the rotor angular velocity. Ji and Leung ⁵ studied the static bifurcation for the same model studied in Ref.[4] but at the superharmonic resonance condition. Yang et al. ⁶ discussed the oscillation modes of the eight-pole system using the energy and phase angle methods. They found that the system can perform one of two vibration modes, which are the quasiperiodic and elliptic motions. El-Shourbagy et al. ⁷ investigated the vibration control of the eight-pole system utilizing a nonlinear form of the PD-controller. They illustrated analytically and numerically that the nonlinear PD-controller is more efficient than the linear one. Saeed et al. ⁸ investigated a discontinuous dynamical model for the eight-pole system, where the rub-impact force among the rotor and the electromagnetic poles are included in the proposed model. They reported that the system can vibrate with periodic, quasiperiodic, or chaotic oscillations depending on the impact stiffness coefficient as well as the dynamic friction coefficient. Zhang et al.^9–14 introduced a PD-controller with time-varying proportional gain to control the eight-pole magnetic bearings system. The system mathematical model has been derived as a time-varying dynamical system with two degrees of freedom based on the proposed control method. Then, they investigated the proposed model with different analytical and numerical methods. The introduced analysis showed that the studied time-varying model may perform multipulse chaotic behavior. Recently, El-Shourbagy et al. ¹⁵ studied the twelve-pole magnetic bearings configuration. They applied the traditional linear position-velocity controller as an active control strategy to suppress the system’s lateral oscillations. According to the introduced analysis, they found that the position control gain has a significant impact on the dynamical characteristics of the system. In addition, it is found that the system can exhibit chaotic and quasiperiodic oscillations besides the periodic ones. Saeed and Kandil ¹⁶ explored the steady-state bifurcation behaviors of the sixteen-pole configuration. They proposed two control algorithms based on the PD-controller to improve the rotor dynamics. The first control strategy has been established based on the rotor Cartesian vibrations in the plane, while the second control method has been designed according to the rotor radial oscillations. According to the introduced analysis, they reported that the system may have simultaneous bi-stable, tri-stable, or quadric-stable periodic oscillations depending on the rotor angular velocity in the case of the first control strategy, while the system can perform simultaneous bi-stable solutions only in the second control method. Zhang et al.^17–21 studied a time-varying dynamical model for the sixteen-pole configuration, where the position-velocity controller with time-varying coefficients has been utilized to control the system’s undesired vibrations. The controllability, durability, and reliability of the active magnetic bearings system have encouraged many scientists and engineers to utilize such systems with some adaptive control techniques to serve as smart actuators.^22–27 Ishida and Inoue ²² utilized a four-pole magnetic bearings system with the push-pull mechanism to control the undesired oscillations of a Jeffcott rotor system supported by the traditional bearings. Saeed et al.^23–27 applied various control algorithms to control the lateral oscillation of conventionally supported rotating shafts employing the magnetic bearings system as an actuator. Finally, Srinivas et al. ²⁸ summarized the different usages of the magnetic bearings system regarding the flexible rotor dynamics.

One of the efficient control approaches that have been widely used to reduce the nonlinear oscillation of various dynamical systems is the Integral Resonant Controller (IRC). IRC is a first-order filter that is linearly coupled to the targeted oscillatory system to dampen these unwanted oscillations.^29–36 Diaz et al., ²⁹ and Al-Mamun et al. ³⁰ employed the IRC to reduce the nonlinear vibrations of a light-weight smart structure, where the authors reported that the IRC has vibration suppression performance higher than the conventional control methods. MacLean and Sumeet ³¹ and Omidi and Mahmoodi ³² applied a nonlinear version of the IRC to reduce the transversal oscillation of a cantilever beam. Omidi and Mahmoodi^33,34 utilized a combination of both the IRC and PPF controller to enhance the control efficiency of the PPF controller. Saeed et al. ³⁵ utilized a time-delayed nonlinear IRC to control the parametric oscillation of a cantilever beam system. The authors concluded that the optimal design for the control parameter can not only mitigate the system oscillations to an acceptable limit but also it can eliminate the system vibrations. Despite the IRC’s high efficiency in suppressing unwanted oscillatory motion in many dynamical systems, it has never been used to suppress the nonlinear oscillations of a 12-pole active magnetic bearings system.

Based on the above discussions, one can notice that the nonlinear vibration control of the twelve-pole magnetic bearing system has not been studied before (except in Ref.[15]). However, this type of magnetic bearing has several advantages over other configurations, including high dynamic stiffness, low power consumption, a negligible cross-coupling effect, and better suspension behaviors. ³⁷ As a result, a new control algorithm consists of both the IRC and PD-controllers is proposed within this article to control the unwanted lateral oscillations of the twelve-pole system. According to the proposed controller, the system mathematical model has been established as two nonlinear second-order differential equations coupled to two linear first-order differential equations. Then, the obtained dynamical system has been investigated by applying perturbation methods. Many Perturbation techniques have been applied in the literature to derive the amplitude-frequency relationship of different nonlinear dynamical systems. For example, the multiple time-scales perturbation technique^38,39 has been applied extensively to investigate nonlinear oscillatory systems with ordinary derivatives. In addition, Chinese mathematician Prof. He et al.^40–42 introduced modified versions of the homotopy perturbation method to analyze different classes of nonlinear oscillators. Recently, Prof. He et al. ⁴³ introduced a simple analytical technique relying on the homotopy perturbation method to analyze the fractal oscillators. In this study, the multiple time-scale perturbation method is applied to investigate the dynamical behaviors of the considered system. The corresponding nonlinear algebraic equations that govern the steady-state motions are extracted. Then, the influence of all control parameters on the static bifurcation behaviors of the rotor system is explored through the different response curves. In addition, the system transient response is also investigated by solving the equations of motion numerically. The obtained results demonstrated that the introduced control method is more efficient than the traditional PD-controller, where a new objective function has been introduced to design the optimum control gains.

Nonlinear dynamical equations of the 12-pole system

This section is dedicated to derive the mathematical model that governs the dynamics of the 12-pole system controlled by the IRC and PD controllers. It is supposed that the considered rotor performs a planner motion in the $X-Y$ plane as illustrated in Figure 1(a). Accordingly, the equations of motion in $X$ and $Y$ directions can be expressed as follows^44,45

$$m\ddot{x}\left(t\right)=me{y}^2\cos \left(yt\right)+F_X$$

(1)

$$m\ddot{y}\left(t\right)=me{\psi }^{2}sin\left(\psi t\right)+F_Y$$

(2)

where $m$ denote the rotor mass, $e$ represents the rotor eccentricity, $\psi$ is the angular velocity, and $F_X$ and $F_Y$ are the net attractive control forces that are applied in the $X$ and $Y$ directions using the 12 electromagnetic poles as shown in Figure 1(b).

Figure 1.

(a) 12-pole system schematic diagram and (b) closed-loop control system.

It is considered that each opposite two poles are designed to produce push-pull electromagnetic force to suspend the rotating shaft within the origin $O$ as shown in Figure 1(a) regardless of the external loads (i.e., the centrifugal force $me{\psi }^2$ ). Accordingly, the push–pull control force of each pair of the opposite poles can be written such that ⁴⁵

$$f_i=H\left[\frac{{\left(I_0-I_j\right)}^2}{{\left(c_0-d_j\right)}^2}-\frac{{\left(I_0+I_j\right)}^2}{{\left(c_0+d_j\right)}^2}\right]\hspace{0.17em};\hspace{0.17em}j=1,\dots ,6$$

(3)

where $H$ is constant ( $H=\frac{1}{4}{\mu }_0{N}^{2}Acos\left(\phi \right)$ ), $I_0$ is a constant permanentized electrical current, $I_j$ is the control current of the j^th pole (It will be defined later according to the proposed control law), and $c_0$ denotes the air-gap size when the geometric center of the rotor coincides with the origin $O$ . $c_0\mp {\delta }_j$ is the instantaneous air-gap of the j^th pole and its opposite one when the rotor system performs Cartesian displacements $x\left(t\right)$ and $y\left(t\right)$ . According to Figure 1(a), ${\delta }_j, (j=1,\dots ,6)$ can be defined as follows

$\begin{array}{c}{\delta }_1=x\cos \left(\alpha \right)-y\sin \left(\alpha \right),{\delta }_2=x,{\delta }_3=x\cos \left(\alpha \right)+y\sin \left(\alpha \right)\\ {\delta }_4=x\sin \left(\alpha \right)+y\cos \left(\alpha \right),{\delta }_5=y,{\delta }_6=-x\sin \left(\alpha \right)+y\cos \left(\alpha \right)\end{array}\}$

(4)

It is supposed that the control forces $f_1, f_2,$ and $f_3$ are functions of the rotor horizontal displacement $x\left(t\right)$ and the corresponding velocity $\dot{x}\left(t\right)$ , while the forces $f_4, f_5,$ and $f_6$ are a function of the rotor vertical displacement $y\left(t\right)$ and the corresponding velocity $\dot{y}\left(t\right)$ . Therefore, we can write

$$I_1=I_2=I_3=I_x,\hspace{0.17em}{I}_4=I_5=I_6=I_y$$

(5)

where $I_x$ denotes the control current in $X$ direction and $I_y$ is the control current in $Y$ direction. These control currents are designed in this article to be generated via combined IRC and PD controllers. So, we can express the control currents ( $I_x$ and $I_y$ ) as follows

$$I_x=k_{1}x+k_2\dot{x}+k_3{q}_x,\hspace{0.17em}{I}_y=k_{1}y+k_2\dot{y}+k_4{q}_y$$

(6)

where $k_1$ and $k_2$ represent the position and velocity control gains, respectively, $k_3$ denotes the control signal gain of the IRC in $X-$ direction, and $k_4$ is the control signal gain of the IRC $Y-$ direction. The IRC is a first-order linear dynamical system that is coupled to the targeted system linearly.^29–36 Accordingly, the IRC dynamical model can be expressed as

$${\dot{q}}_x+{\omega }_x\hspace{0.17em}{q}_x=k_{5}x$$

(7)

$${\dot{q}}_y+{\omega }_y\hspace{0.17em}{q}_y=k_{6}y$$

(8)

where $k_5$ and $k_6$ are the IRC feedback signal gains, and ${\omega }_x$ and ${\omega }_y$ are the internal loop feedback gain. Substituting equations (4)–(6) into equation (3), yields

$$f_1=H\left(\frac{{\left(I_0-k_{1}x-k_2\dot{x}-k_3{q}_x\right)}^2}{{\left(c_0-x\cos \left(\alpha \right)+y\sin \left(\alpha \right)\right)}^2}-\frac{{\left(I_0+k_{1}x+k_2\dot{x}+k_3{q}_x\right)}^2}{{\left(c_0+x\cos \left(\alpha \right)-y\sin \left(\alpha \right)\right)}^2}\right)$$

(9)

$$f_2=H\left(\frac{{\left(I_0-k_{1}x-k_2\dot{x}-k_3{q}_x\right)}^2}{{\left(c_0-x\right)}^2}-\frac{{\left(I_0+k_{1}x+k_2\dot{x}+k_3{q}_x\right)}^2}{{\left(c_0+x\right)}^2}\right)$$

(10)

$$f_3=H\left(\frac{{\left(I_0-k_{1}x-k_2\dot{x}-k_3{q}_x\right)}^2}{{\left(c_0-x\cos \left(\alpha \right)-y\sin \left(\alpha \right)\right)}^2}-\frac{{\left(I_0+k_{1}x+k_2\dot{x}+k_3{q}_x\right)}^2}{{\left(c_0+x\cos \left(\alpha \right)+y\sin \left(\alpha \right)\right)}^2}\right)$$

(11)

$$f_4=H\left(\frac{{\left(I_0-k_{1}y-k_2\dot{y}-k_4{q}_y\right)}^2}{{\left(c_0-x\sin \left(\alpha \right)-y\cos \left(\alpha \right)\right)}^2}-\frac{{\left(I_0+k_{1}y+k_2\dot{y}+k_4{q}_y\right)}^2}{{\left(c_0+x\sin \left(\alpha \right)+y\cos \left(\alpha \right)\right)}^2}\right)$$

(12)

$$f_5=H\left(\frac{{\left(I_0-k_{1}y-k_2\dot{y}-k_4{q}_y\right)}^2}{{\left(c_0-y\right)}^2}-\frac{{\left(I_0+k_{1}y+k_2\dot{y}+k_4{q}_y\right)}^2}{{\left(c_0+y\right)}^2}\right)$$

(13)

$$f_6=H\left(\frac{{\left(I_0-k_{1}y-k_2\dot{y}-k_4{q}_y\right)}^2}{{\left(c_0+x\sin \left(\alpha \right)-y\cos \left(\alpha \right)\right)}^2}-\frac{{\left(I_0+k_{1}y+k_2\dot{y}+k_4{q}_y\right)}^2}{{\left(c_0-x\sin \left(\alpha \right)+y\cos \left(\alpha \right)\right)}^2}\right)$$

(14)

According to the geometry of Figure 1(a), The net attractive control forces $F_X$ and $F_Y$ can be expressed as

$$F_X=f_2+(f_1+f_3)\cos \left(\alpha \right)+(f_4-f_6)\sin \left(\alpha \right)$$

(15)

$$F_Y=f_5+(f_4+f_6)\cos \left(\alpha \right)+(f_3-f_1)\sin \left(\alpha \right)$$

(16)

Substituting equations (9)–(16) in equations (1) and (2), with introducing the dimensionless variables: $t^{\ast }=\zeta t, x^{\ast }=\frac{x}{c_0}, y^{\ast }=\frac{y}{c_0}, q_x^{\ast }=\frac{q_x}{c_0}, q_y^{\ast }=\frac{q_y}{c_0}, p=\frac{c_0}{I_0}{k}_1, d=\frac{c_0\zeta }{I_0}{k}_2,{\eta }_1=\frac{c_0}{I_0}\left(1+2\cos \left(\alpha \right)\right)k_3,{\eta }_2=\frac{c_0}{I_0}\left(1+2\cos \left(\alpha \right)\right)k_4 , {\eta }_3=\frac{k_5}{\zeta }, {\eta }_4=\frac{k_6}{\zeta } ,{\omega }_1=\frac{{\omega }_x}{\zeta },{\omega }_2=\frac{{\omega }_y}{\zeta }, \mathrm{\Omega }=\frac{\psi }{\zeta }, f=\frac{e}{c_0}, \zeta =\sqrt{H/m{c}_0^3}$ with and omitting the asterisks for simplicity, we get the whole system model from equations (1), (2), (7), and (8) as follows

$$\begin{array}{l}\ddot{x}+2\mu \dot{x}+{\omega }^{2}x-\left({\alpha }_1{x}^3+{\alpha }_{2}x{y}^2+{\alpha }_3{x}^2\dot{x}+{\alpha }_4\dot{x}{y}^2+{\alpha }_{5}x{\dot{y}}^2+{\alpha }_{6}x{\dot{x}}^2+{\alpha }_{7}xy\dot{y}\right)-({\beta }_1{x}^2{q}_x+{\beta }_{2}x\dot{x}{q}_x\\ +{\beta }_{3}x{q}_x^2+{\beta }_{4}xy{q}_y+{\beta }_{5}x{q}_y^2+{\beta }_{6}x\dot{y}{q}_y+{\beta }_7{q}_x{y}^2)={\mathrm{\Omega }}^{2}f\cos \left(\mathrm{\Omega }t\right)+{\eta }_1{q}_x\end{array}$$

(17)

$$\begin{array}{l}\ddot{y}+2\mu \dot{y}+{\omega }^{2}y-\left({\alpha }_1{y}^3+{\alpha }_{2}y{x}^2+{\alpha }_3{y}^2\dot{y}+{\alpha }_4\dot{y}{x}^2+{\alpha }_{5}y{\dot{x}}^2+{\alpha }_{6}y{\dot{y}}^2+{\alpha }_{7}yx\dot{x}\right)-({\gamma }_1{y}^2{q}_y+{\gamma }_{2}y\dot{y}{q}_y\\ +{\gamma }_{3}y{q}_y^2+{\gamma }_{4}yx{q}_x+{\gamma }_{5}y{q}_x^2+{\gamma }_{6}y\dot{x}{q}_x+{\gamma }_7{q}_y{x}^2)={\mathrm{\Omega }}^{2}f\sin \left(\mathrm{\Omega }t\right)+{\eta }_2{q}_y\end{array}$$

(18)

$${\dot{q}}_x+{\omega }_1\hspace{0.17em}{q}_x={\eta }_{3}x$$

(19)

$${\dot{q}}_y+{\omega }_2\hspace{0.17em}{q}_y={\eta }_{4}y$$

(20)

where, $\mu =d\left[\cos \left(\alpha \right)+\frac{1}{2}\right],\omega =\sqrt{2p\cos \left(\alpha \right)+p-3}$ , ${\alpha }_1=6-6p{\cos }^3\left(\alpha \right)-3p+2p^2{\cos }^2\left(\alpha \right)+p^2+8{\cos }^4\left(\alpha \right)-8{\cos }^2\left(\alpha \right)$ , ${\alpha }_2=2p^2-2p^2{\cos }^2\left(\alpha \right)+24{\cos }^2\left(\alpha \right)-24{\cos }^4\left(\alpha \right)-18p\cos \left(\alpha \right)+18p{\cos }^3\left(\alpha \right)$ , ${\alpha }_3=-3d+2pd+4pd{\cos }^2\left(\alpha \right)-6d{\cos }^3\left(\alpha \right),\hspace{0.17em}{\alpha }_4=6d{\cos }^3\left(\alpha \right)-6d\cos \left(\alpha \right)$ , ${\alpha }_5=2d^2-2d^2{\cos }^2\left(\alpha \right),{\alpha }_6=2d^2{\cos }^2\left(\alpha \right)+d^2,{\alpha }_7=12d{\cos }^3\left(\alpha \right)-12d\cos \left(\alpha \right)+4pd-4pd{\cos }^2\left(\alpha \right)$ , ${\beta }_1=\frac{{\eta }_1}{1+2\cos \left(\alpha \right)}\left[-3+2p-6{\cos }^3\left(\alpha \right)+4p{\cos }^2\left(\alpha \right)\right],{\beta }_2=\frac{{\eta }_1}{1+2\cos \left(\alpha \right)}\left[2d+4d{\cos }^2\left(\alpha \right)\right]$ , ${\beta }_3=\frac{{\eta }_1^2}{{(1+2\cos \left(\alpha \right))}^2}\left[1+2{\cos }^2\left(\alpha \right)\right],{\beta }_4=\frac{{\eta }_2}{1+2\cos \left(\alpha \right)}\left[4psi{n}^2\left(\alpha \right)-12\cos \left(\alpha \right)si{n}^2\left(\alpha \right)\right]$ , ${\beta }_5=\frac{2{\eta }_2^{2}si{n}^2\left(\alpha \right)}{{(1+2\cos \left(\alpha \right))}^2},{\beta }_6=\frac{4d{\eta }_{2}si{n}^2\left(\alpha \right)}{1+2\cos \left(\alpha \right)},{\beta }_7=\frac{-6{\eta }_1\cos \left(\alpha \right)si{n}^2\left(\alpha \right)}{1+2\cos \left(\alpha \right)}$ , ${\gamma }_1=\frac{{\eta }_2}{1+2\cos \left(\alpha \right)}\left[-3+2p-6{\cos }^3\left(\alpha \right)+4p{\cos }^2\left(\alpha \right)\right],{\gamma }_2=\frac{{\eta }_2}{1+2\cos \left(\alpha \right)}\left[2d+4d{\cos }^2\left(\alpha \right)\right]$ , ${\gamma }_3=\frac{{\eta }_2^2}{{(1+2\cos \left(\alpha \right))}^2}\left[1+2{\cos }^2\left(\alpha \right)\right],{\gamma }_4=\frac{{\eta }_1}{1+2\cos \left(\alpha \right)}\left[4psi{n}^2\left(\alpha \right)-12\cos \left(\alpha \right)si{n}^2\left(\alpha \right)\right]$ , ${\gamma }_5=\frac{2{\eta }_1^{2}si{n}^2\left(\alpha \right)}{{(1+2\cos \left(\alpha \right))}^2},{\gamma }_6=\frac{4d{\eta }_{1}si{n}^2\left(\alpha \right)}{1+2\cos \left(\alpha \right)},{\gamma }_7=\frac{-6{\eta }_2\cos \left(\alpha \right)si{n}^2\left(\alpha \right)}{1+2\cos \left(\alpha \right)}$

Based on the introduced dimensionless parameters before equations (17)–(20), we can deduce that ${\eta }_1=\frac{c_0}{I_0}\left(1+2\cos \left(\alpha \right)\right)k_3$ and ${\eta }_2=\frac{c_0}{I_0}\left(1+2\cos \left(\alpha \right)\right)k_4$ denote the dimensionless forms of the IRCs control gains ( $k_3$ and $k_4$ ). In addition, ${\eta }_3=\frac{k_5}{\zeta }$ and ${\eta }_4=\frac{k_6}{\zeta }$ represent the dimensionless forms of the IRCs feedback gains ( $k_5$ and $k_6$ ). Moreover, ${\omega }_1=\zeta {\omega }_x$ and ${\omega }_2=\zeta {\omega }_y$ are the dimensionless forms of the IRCs internal loop feedback gains ( ${\omega }_x$ and ${\omega }_y$ ). Also, $p=\frac{c_0}{I_0}{k}_1$ and $d=\frac{c_0\zeta }{I_0}{k}_2$ are the dimensionless forms of the proportional ( $k_1$ ) and derivative ( $k_2$ ) gains of the PD-controller. To explore the effects of the various controller gains (i.e., ${\eta }_1, {\eta }_2, {\eta }_3, {\eta }_4,{\omega }_1,{\omega }_2, p,$ and $d$ ) on the oscillatory motion of the considered system, an analytical approximate solution for the derived dynamical system (i.e., equations (17)–(20)) is obtained using as given below.

Analytical investigation and autonomous amplitude-phase equations

Relying on the perturbation theory, a first-order approximate solution for the derived equations of motions (i.e., equations (17)–(20)) can be expressed as follows ^{38, 39}

$$x\left(t,\epsilon \right)=x_0\left(T_0,T_1\right)+\epsilon x_1\left(T_0,T_1\right)$$

(21)

$$y\left(t,\epsilon \right)=y_0\left(T_0,T_1\right)+\epsilon y_1\left(T_0,T_1\right)$$

(22)

$$q_x\left(t,\epsilon \right)=\epsilon q_{x0}\left(T_0,T_1\right)+{\epsilon }^2{q}_{x1}\left(T_0,T_1\right)$$

(23)

$$q_y\left(t,\epsilon \right)=\epsilon q_{y0}\left(T_0,T_1\right)+{\epsilon }^2{q}_{y1}\left(T_0,T_1\right)$$

(24)

where $\epsilon$ is the perturbation parameter, $T_0=t$ , and $T_1=\epsilon t$ . According to the time scales $T_0$ and $T_1$ , the derivatives $\frac{d}{dt}$ and $\frac{d^2}{d{t}^2}$ can be expressed as

$$\frac{d}{dt}=D_0+\epsilon D_1,\frac{d^2}{d{t}^2}=D_0^2+2\epsilon D_0{D}_1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{D}_j=\frac{\partial }{\partial T_j},j=\mathrm{0,1}$$

(25)

To perform the solution procedure, the system parameters may be scaled such that

$$\mu =\epsilon \tilde{\mu },\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\alpha }_j=\epsilon {\tilde{\alpha }}_j,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}f=\epsilon \tilde{f}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\eta }_3=\epsilon {\tilde{\eta }}_3,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\eta }_4=\epsilon {\tilde{\eta }}_4,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,\mathrm{...},7$$

(26)

Inserting equations (21)–(26) into equations (17)–(20), and then equating the coefficients that have the same power of $\epsilon$ , yields O ( ${\epsilon }^0$ )

$$(D_0^2+{\omega }^2)x_0=0$$

(27)

$$(D_0^2+{\omega }^2)y_0=0$$

(28)

O ( ${\epsilon }^1$ )

$$\left(D_0+{\omega }_1\right)q_{x0}={\tilde{\eta }}_3{x}_0$$

(29)

$$\left(D_0+{\omega }_2\right)q_{y0}={\tilde{\eta }}_4{y}_0$$

(30)

$$\begin{array}{c}\left(D_0^2+{\omega }^2\right)x_1=-2D_0{D}_1{x}_0-2\tilde{\mu }{D}_0{x}_0+{\tilde{\alpha }}_1{x}_0^3+{\tilde{\alpha }}_2{x}_0{y}_0^2+{\tilde{\alpha }}_3{x}_0^2{D}_0{x}_0+{\tilde{\alpha }}_4{y}_0^2{D}_0{x}_0+{\tilde{\alpha }}_5{x}_0{\left(D_0{y}_0\right)}^2\\ +{\tilde{\alpha }}_6{x}_0{\left(D_0{x}_0\right)}^2+{\tilde{\alpha }}_7{x}_0{y}_0{D}_0{y}_0+{\beta }_1{x}_0^2{q}_{x0}+{\beta }_2{x}_0{D}_0{x}_0{q}_{x0}+{\beta }_3{x}_0{q}_{x0}^2+{\beta }_4{x}_0{y}_0{q}_{y0}+{\beta }_5{x}_0{q}_{y0}^2\\ +{\beta }_6{x}_0{D}_0{y}_0{q}_{y0}+{\beta }_7{y}_0^2{q}_{x0}+{\mathrm{\Omega }}^2\tilde{f}\cos \left(\mathrm{\Omega }t\right)+{\eta }_1{q}_{x0}\end{array}$$

(31)

$$\begin{array}{c}\left(D_0^2+{\omega }^2\right)y_1=-2D_0{D}_1{y}_0-2\tilde{\mu }{D}_0{y}_0+{\tilde{\alpha }}_1{y}_0^3+{\tilde{\alpha }}_2{y}_0{x}_0^2+{\tilde{\alpha }}_3{y}_0^2{D}_0{y}_0+{\tilde{\alpha }}_4{x}_0^2{D}_0{y}_0+{\tilde{\alpha }}_5{y}_0{\left(D_0{x}_0\right)}^2\\ +{\tilde{\alpha }}_6{y}_0{\left(D_0{y}_0\right)}^2+{\tilde{\alpha }}_7{y}_0{x}_0{D}_0{x}_0+{\gamma }_1{y}_0^2{q}_{y0}+{\gamma }_2{y}_0{D}_0{y}_0{q}_{y0}+{\gamma }_3{y}_0{q}_{y0}^2+{\gamma }_4{y}_0{x}_0{q}_{x0}+{\gamma }_5{y}_0{q}_{x0}^2\\ +{\gamma }_6{y}_0{D}_0{x}_0{q}_{x0}+{\gamma }_7{x}_0^2{q}_{y0}+{\mathrm{\Omega }}^2\tilde{f}\sin \left(\mathrm{\Omega }t\right)+{\eta }_2{q}_{y0}\end{array}$$

(32)

The solution of equations (27) and (28) is

$$x_0\left(T_0,T_1\right)=A\left(T_1\right)e^{i\omega T_0}+\overline{A}\left(T_1\right)e^{-i\omega T_0}$$

(33)

$$y_0\left(T_0,T_1\right)=B\left(T_1\right)e^{i\omega T_0}+\overline{B}\left(T_1\right)e^{-i\omega T_0}$$

(34)

Substituting equations (33) and (34) into the right-hand side of equations (29) and (30), and then solving the resulting equation, we have

$$q_{x0}\left(T_0,T_1\right)={\lambda }_{1}A\left(T_1\right)e^{i\omega T_0}+{\overline{\lambda }}_1\overline{A}\left(T_1\right)e^{-i\omega T_0}$$

(35)

$$q_{y0}\left(T_0,T_1\right)={\lambda }_{2}B\left(T_1\right)e^{i\omega T_0}+{\overline{\lambda }}_2\overline{B}\left(T_1\right)e^{-i\omega T_0}$$

(36)

where $i=\sqrt{-1}, {\lambda }_1=\frac{{\omega }_1-i\omega }{{\omega }_1^2+{\omega }^2}{\tilde{\eta }}_3, {\lambda }_2=\frac{{\omega }_2-i\omega }{{\omega }_2^2+{\omega }^2}{\tilde{\eta }}_4$ , and the unknown functions $A\left(T_1\right)$ , $B\left(T_1\right)$ will be determined later. Now, inserting equations (33)–(36) into the right-hand side of equations (31) and (32) into equations (31) and (32), yields

$\begin{array}{c}(D_0^2+{\omega }^2)x_1=[-2i\omega \left(D_{1}A\right)-2i\tilde{\mu }\omega A+3{\tilde{\alpha }}_1{A}^2\overline{A}+2{\tilde{\alpha }}_{2}AB\overline{B}+{\tilde{\alpha }}_2\overline{A}{B}^2+i{\tilde{\alpha }}_3\omega A^2\overline{A}+2i{\tilde{\alpha }}_4\omega AB\overline{B}\\ -i{\tilde{\alpha }}_4\omega \overline{A}{B}^2+2{\tilde{\alpha }}_5{\omega }^{2}AB\overline{B}-{\tilde{\alpha }}_5{\omega }^2\overline{A}{B}^2+{\tilde{\alpha }}_6{\omega }^2{A}^2\overline{A}+i{\tilde{\alpha }}_7\omega \overline{A}{B}^2+2{\beta }_1{\lambda }_1{A}^2\overline{A}+{\beta }_3{\lambda }_1^2{A}^2\overline{A}\\ +{\beta }_4({\lambda }_{2}AB\overline{B}+{\lambda }_2\overline{A}{B}^2)+{\beta }_5{\lambda }_2^2\overline{A}{B}^2+{\beta }_6(-i\omega {\lambda }_{2}AB\overline{B}+i\omega {\lambda }_2\overline{A}{B}^2)+2{\beta }_7{\lambda }_{1}AB\overline{B}\\ +{\eta }_1{\lambda }_{1}A]e^{i\omega T_0}+[{\widehat{\alpha }}_1{A}^3+{\tilde{\alpha }}_{2}A{B}^2+i{\tilde{\alpha }}_3\omega A^3+i{\tilde{\alpha }}_4\omega A{B}^2-{\tilde{\alpha }}_5{\omega }^{2}A{B}^2-{\tilde{\alpha }}_6{\omega }^2{A}^3\\ +i{\tilde{\alpha }}_7\omega A{B}^2+{\beta }_1{\lambda }_1{A}^3+i{\beta }_2\omega {\lambda }_1{A}^3{e}^{3i\omega T_0}+{\beta }_3{\lambda }_1^2{A}^3{e}^{3i\omega T_0}+{\beta }_4{\lambda }_{2}A{B}^2+{\beta }_5{\lambda }_2^{2}A{B}^2\\ +i{\beta }_6\omega {\lambda }_{2}A{B}^2+{\beta }_7{\lambda }_{1}A{B}^2]e^{3i\omega T_0}+\frac{1}{2}{\mathrm{\Omega }}^2\tilde{f}{e}^{i\Omega T_0}+cc\end{array}$

(37)

$\begin{array}{c}(D_0^2+{\omega }^2)y_1=[-2i\omega \left(D_{1}B\right)-2i\tilde{\mu }\omega B+3{\tilde{\alpha }}_1{B}^2\overline{B}+2{\tilde{\alpha }}_{2}BA\overline{A}+{\tilde{\alpha }}_2\overline{B}{A}^2+i{\tilde{\alpha }}_3\omega B^2\overline{B}+2i{\tilde{\alpha }}_4\omega BA\overline{A}\\ -i{\tilde{\alpha }}_4\omega \overline{B}{A}^2+2{\tilde{\alpha }}_5{\omega }^{2}BA\overline{A}-{\tilde{\alpha }}_5{\omega }^2\overline{B}{A}^2+{\tilde{\alpha }}_6{\omega }^2{B}^2\overline{B}+i{\tilde{\alpha }}_7\omega \overline{B}{A}^2+2{\gamma }_1{\lambda }_2{B}^2\overline{B}+{\gamma }_3{\lambda }_2^2{B}^2\overline{B}\\ +{\gamma }_4({\lambda }_{1}BA\overline{A}+{\lambda }_1\overline{B}{A}^2)+{\gamma }_5{\lambda }_1^2\overline{B}{A}^2+{\gamma }_6(-i\omega {\lambda }_{1}BA\overline{A}+i\omega {\lambda }_1\overline{B}{A}^2)+2{\gamma }_7{\lambda }_{2}BA\overline{A}\\ +{\eta }_2{\lambda }_{2}B]e^{i\omega T_0}+[{\tilde{\alpha }}_1{B}^3+{\tilde{\alpha }}_{2}B{A}^2+i{\tilde{\alpha }}_3\omega B^3+i{\widehat{\alpha }}_4\omega B{A}^2-{\tilde{\alpha }}_5{\omega }^{2}B{A}^2-{\tilde{\alpha }}_6{\omega }^2{B}^3\\ +i{\tilde{\alpha }}_7\omega B{A}^2+{\gamma }_1{\lambda }_2{B}^3+{\gamma }_{2}i\omega {\lambda }_2{B}^3{e}^{3i\omega T_0}+{\gamma }_3{\lambda }_2^2{B}^3+{\gamma }_4{\lambda }_{1}B{A}^2{e}^{3i\omega T_0}+{\gamma }_5{\lambda }_1^{2}B{A}^2\\ +{\gamma }_{6}i\omega {\lambda }_{1}B{A}^2+{\gamma }_7{\lambda }_{2}B{A}^2]e^{3i\omega T_0}-\frac{1}{2}i{\mathrm{\Omega }}^2\tilde{f}{e}^{i\Omega T_0}+cc\end{array}$

(38)

where $cc$ denotes the complex conjugate term. To investigate system oscillatory behavior at the primary resonance conditions (i.e., when $\mathrm{\Omega }=\omega$ ), it is convenient to present the detuning parameter $\sigma$ to describe the closeness $\mathrm{\Omega }$ to $\omega$ as follows

$$\mathrm{\Omega }=\omega +\sigma =\omega +\epsilon \tilde{\sigma }$$

(39)

Inserting equation (39) into equations (37) and (38), we can extract the following secular terms^38,39

$$\begin{array}{l}-2i\omega \left(D_{1}A\right)-2i\tilde{\mu }\omega A+3{\tilde{\alpha }}_1{A}^2\overline{A}+2{\tilde{\alpha }}_{2}AB\overline{B}+{\tilde{\alpha }}_2\overline{A}{B}^2+i{\tilde{\alpha }}_3\omega A^2\overline{A}+2i{\tilde{\alpha }}_4\omega AB\overline{B}-i{\tilde{\alpha }}_4\omega \overline{A}{B}^2\\ +2{\tilde{\alpha }}_5{\omega }^{2}AB\overline{B}-{\tilde{\alpha }}_5{\omega }^2\overline{A}{B}^2+{\tilde{\alpha }}_6{\omega }^2{A}^2\overline{A}+i{\tilde{\alpha }}_7\omega \overline{A}{B}^2+2{\beta }_1{\lambda }_1{A}^2\overline{A}+{\beta }_3{\lambda }_1^2{A}^2\overline{A}+{\beta }_4{\lambda }_{2}AB\overline{B}\\ +{\beta }_4{\lambda }_2\overline{A}{B}^2+{\beta }_5{\lambda }_2^2\overline{A}{B}^2-i{\beta }_6\omega {\lambda }_{2}AB\overline{B}+i{\beta }_6\omega {\lambda }_2\overline{A}{B}^2+2{\beta }_7{\lambda }_{1}AB\overline{B}+{\eta }_1{\lambda }_{1}A+\frac{1}{2}{\left(\omega +\sigma \right)}^2\tilde{f}{e}^{i\epsilon \tilde{\sigma }{T}_0}=0\end{array}$$

(40)

$$\begin{array}{l}-2i\omega \left(D_{1}B\right)-2i\tilde{\mu }\omega B+3{\tilde{\alpha }}_1{B}^2\overline{B}+2{\tilde{\alpha }}_{2}BA\overline{A}+{\tilde{\alpha }}_2\overline{B}{A}^2+i{\tilde{\alpha }}_3\omega B^2\overline{B}+2i{\tilde{\alpha }}_4\omega BA\overline{A}-i{\tilde{\alpha }}_4\omega \overline{B}{A}^2\\ +2{\tilde{\alpha }}_5{\omega }^{2}BA\overline{A}-{\tilde{\alpha }}_5{\omega }^2\overline{B}{A}^2\hspace{0.17em}+{\tilde{\alpha }}_6{\omega }^2{B}^2\overline{B}+i{\tilde{\alpha }}_7\omega \overline{B}{A}^2+2{\gamma }_1{\lambda }_2{B}^2\overline{B}+{\gamma }_3{\lambda }_2^2{B}^2\overline{B}+{\gamma }_4{\lambda }_{1}BA\overline{A}\\ +{\gamma }_4{\lambda }_1\overline{B}{A}^2+{\gamma }_5{\lambda }_1^2\overline{B}{A}^2-i{\gamma }_6\omega {\lambda }_{1}BA\overline{A}+i{\gamma }_6\omega {\lambda }_1\overline{B}{A}^2+2{\gamma }_7{\lambda }_{2}BA\overline{A}+{\eta }_2{\lambda }_{2}B-\frac{1}{2}i{(\omega +\sigma )}^2\tilde{f}{e}^{i\epsilon \tilde{\sigma }{T}_0}=0\end{array}$$

(41)

To obtain the slow-flow modulating equations of the considered 12-pole system, we can express the unknown functions $A\left(T_1\right)$ and $B\left(T_1\right)$ in their polar form as^38,39

$$A=\frac{1}{2}a{e}^{i{\theta }_1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}B=\frac{1}{2}b{e}^{i{\theta }_2}\Rightarrow D_{1}A=\frac{d}{\epsilon dt}A=\frac{1}{2\epsilon }\left(\dot{a}{e}^{i{\theta }_1}+ia{\dot{\theta }}_1{e}^{i{\theta }_1}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{D}_{1}B=\hspace{0.17em}\hspace{0.17em}\frac{d}{\epsilon dt}B=\frac{1}{2\epsilon }\left(\dot{b}{e}^{i{\theta }_2}+ib{\dot{\theta }}_2{e}^{i{\theta }_2}\right)$$

(42)

Inserting equation (42) into equations (40) and (41), with separating the real and imaginary parts, one can obtain the following nonlinear autonomous dynamical system

$$\begin{array}{c}\dot{a}=-\left[\mu +\frac{{\eta }_1{\eta }_3}{2\left({\omega }_1^2+{\omega }^2\right)}\right]a+\frac{1}{8}\left[{\alpha }_3-\frac{2{\beta }_1{\eta }_3}{{\omega }_1^2+{\omega }^2}-\frac{2{\beta }_3{\eta }_3^2{\omega }_1}{{\left({\omega }_1^2+{\omega }^2\right)}^2}\right]a^3+\frac{1}{8}[2{\alpha }_4-\frac{{\beta }_4{\eta }_4}{{\omega }_2^2+{\omega }^2}-\frac{{\beta }_6{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}\\ -\frac{2{\beta }_7{\eta }_3}{{\omega }_1^2+{\omega }^2}]a{b}^2+\frac{1}{8}\left[-{\alpha }_4+{\alpha }_7-\frac{{\beta }_4{\eta }_4}{{\omega }_2^2+{\omega }^2}+\frac{2{\beta }_5{\eta }_4^2{\omega }_2}{{\left({\omega }_1^2+{\omega }^2\right)}^2}+\frac{{\beta }_6{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}\right]a{b}^2\cos (2{\varphi }_2-2{\varphi }_1)\\ -\frac{1}{8}\left[\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega +\frac{{\beta }_4{\eta }_4{\omega }_2}{\omega \left({\omega }_1^2+{\omega }^2\right)}+\frac{{\beta }_5{\eta }_4^2\left({\omega }_1^2+{\omega }^2\right)}{\omega {\left({\omega }_1^2+{\omega }^2\right)}^2}-\frac{{\beta }_6{\eta }_4\omega }{{\omega }_2^2+{\omega }^2}\right]a{b}^2\sin \left(2{\varphi }_2-2{\varphi }_1\right)+\frac{f}{2\omega }{\left(\omega +\sigma \right)}^2\sin \left({\varphi }_1\right)\end{array}$$

(43)

$$\begin{array}{c}\dot{b}=-\left[\mu +\frac{{\eta }_2{\eta }_4}{2\left({\omega }_2^2+{\omega }^2\right)}\right]b+\frac{1}{8}\left[{\alpha }_3-\frac{2{\gamma }_1{\eta }_4}{{\omega }_2^2+{\omega }^2}-\frac{2{\gamma }_3{\eta }_4^2{\omega }_2}{{\left({\omega }_2^2+{\omega }^2\right)}^2}\right]b^3+\frac{1}{8}[2{\alpha }_4-\frac{{\gamma }_4{\eta }_3}{{\omega }_1^2+{\omega }^2}-\frac{{\gamma }_6{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}\\ -\frac{2{\gamma }_7{\eta }_4}{{\omega }_2^2+{\omega }^2}]a^{2}b+\frac{1}{8}\left[-{\alpha }_4+{\alpha }_7-\frac{{\gamma }_4{\eta }_3}{{\omega }_1^2+{\omega }^2}-\frac{2{\gamma }_5{\eta }_3^2{\omega }_1}{{\left({\omega }_2^2+{\omega }^2\right)}^2}+\frac{{\gamma }_6{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}\right]a^{2}b\cos \left(2{\varphi }_2-2{\varphi }_1\right)\\ +\frac{1}{8}\left[\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega +\frac{{\gamma }_4{\eta }_3{\omega }_1}{\omega \left({\omega }_2^2+{\omega }^2\right)}+\frac{{\gamma }_5{\eta }_3^2({\omega }_1^2-{\omega }^2)}{\omega {\left({\omega }_2^2+{\omega }^2\right)}^2}+\frac{{\gamma }_6{\eta }_3\omega }{{\omega }_1^2+{\omega }^2}\right]a^{2}b\sin \left(2{\varphi }_2-2{\varphi }_1\right)-\frac{f}{2\omega }{\left(\omega +\sigma \right)}^2\cos \left({\varphi }_2\right)\end{array}$$

(44)

$$\begin{array}{c}{\dot{\varphi }}_1=\left[\sigma +\frac{{\eta }_1{\eta }_3{\omega }_1}{2\omega \left({\omega }_1^2+{\omega }^2\right)}\right]+\frac{1}{8\omega }\left[3{\alpha }_1+{\alpha }_6{\omega }^2+\frac{2{\beta }_1{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}+\frac{{\beta }_3{\eta }_3^2({\omega }_1^2-{\omega }^2)}{{({\omega }_1^2+{\omega }^2)}^2}\right]a^2\\ +\frac{1}{8\omega }\left[2{\alpha }_2+2{\alpha }_5{\omega }^2+\frac{{\beta }_4{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}-\frac{{\beta }_6{\eta }_4{\omega }^2}{{\omega }_2^2+{\omega }^2}+\frac{2{\beta }_7{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}\right]b^2+\frac{1}{8\omega }[{\alpha }_2-{\alpha }_5{\omega }^2+\frac{{\beta }_4{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}\\ +\frac{{\beta }_5{\eta }_4^2\left({\omega }_2^2+{\omega }^2\right)}{{\left({\omega }_2^2+{\omega }^2\right)}^2}+\frac{{\beta }_6{\eta }_4{\omega }^2}{{\omega }_2^2+{\omega }^2}]b^2\cos \left(2{\varphi }_2-2{\varphi }_1\right)+\frac{1}{8}[-{\alpha }_4+{\alpha }_7-\frac{{\beta }_4{\eta }_4}{{\omega }_2^2+{\omega }^2}-\frac{2{\beta }_5{\eta }_4^2{\omega }_2}{{\left({\omega }_2^2+{\omega }^2\right)}^2}\\ +\frac{{\beta }_6{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}]b^2\sin \left(2{\varphi }_2-2{\varphi }_1\right)+\frac{f}{2a\omega }{\left(\omega +\sigma \right)}^2\cos \left({\varphi }_1\right)\end{array}$$

(45)

$$\begin{array}{c}{\dot{\varphi }}_2=\left[\sigma +\frac{{\eta }_2{\eta }_4{\omega }_2}{2\omega \left({\omega }_2^2+{\omega }^2\right)}\right]+\frac{1}{8\omega }\left[3{\alpha }_1+{\alpha }_6{\omega }^2+\frac{2{\gamma }_1{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}+\frac{{\gamma }_3{\eta }_4^2\left({\omega }_2^2+{\omega }^2\right)}{{\left({\omega }_2^2+{\omega }^2\right)}^2}\right]b^2\\ +\frac{1}{8\omega }\left[2{\alpha }_2+2{\alpha }_5{\omega }^2+\frac{{\gamma }_4{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}-\frac{{\gamma }_6{\eta }_3{\omega }^2}{{\omega }_1^2+{\omega }^2}+\frac{2{\gamma }_7{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}\right]a^2+\frac{1}{8\omega }[{\alpha }_2-{\alpha }_5{\omega }^2+\frac{{\gamma }_4{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}\\ +\frac{{\gamma }_5{\eta }_3^2\left({\omega }_1^2+{\omega }^2\right)}{{\left({\omega }_1^2+{\omega }^2\right)}^2}+\frac{{\gamma }_6{\eta }_3{\omega }^2}{{\omega }_1^2+{\omega }^2}]a^2\cos (2{\varphi }_2-2{\varphi }_1)-\frac{1}{8}[-{\alpha }_4+{\alpha }_7+\frac{{\gamma }_4{\eta }_3}{{\omega }_1^2+{\omega }^2}-\frac{2{\gamma }_5{\eta }_3^2{\omega }_1}{{\left({\omega }_1^2+{\omega }^2\right)}^2}\\ -\frac{{\gamma }_6{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}]a^2\sin \left(2{\varphi }_2-2{\varphi }_1\right)+\frac{f}{2b\omega }{\left(\omega +\sigma \right)}^2\sin \left({\varphi }_2\right)\end{array}$$

(46)

ϕ₁ = σt − θ₁, ϕ₂ = σt − θ₂ where Now, by inserting equations (42) and (33)–(36) into equations (21)–(24) taking into account equation (39), we can obtain the periodic solution of equations (17)–(20) as follows

$$x\left(t\right)=a\left(t\right)\cos (\mathrm{\Omega }t-{\varphi }_1\left(t\right))$$

(47)

$$y\left(t\right)=b\left(t\right)\cos (\mathrm{\Omega }t-{\varphi }_2\left(t\right))$$

(48)

$$q_x\left(t\right)=\frac{{\eta }_3}{\left({\omega }_1^2+{\omega }^2\right)}a\left(t\right)\left[{\omega }_1\cos \left(\mathrm{\Omega }t-{\varphi }_1\right)+\omega \sin \left(\mathrm{\Omega }t-{\varphi }_1\right)\right]$$

(49)

$$q_y\left(t\right)=\frac{{\eta }_4}{\left({\omega }_2^2+{\omega }^2\right)}b\left(t\right)\left[{\omega }_2\cos \left(\mathrm{\Omega }t-{\varphi }_2\right)+\omega \sin \left(\mathrm{\Omega }t-{\varphi }_2\right)\right]$$

(50)

It is clear from equations (47)–(50) that $a\left(t\right)$ and $b\left(t\right)$ are the vibration amplitudes of the 12-pole rotor in $X$ and $Y$ directions, respectively, and ${\varphi }_1\left(t\right)$ and ${\varphi }_2\left(t\right)$ are the motion phase angles. The nonlinear autonomous system (43)–(46) governs the vibration amplitudes ( $a\left(t\right),$ $b\left(t\right)$ ) and phase angles ( ${\varphi }_1\left(t\right), {\varphi }_2\left(t\right)$ ). Accordingly, the steady-state amplitudes and phase angles of equations (17)–(20) can be obtained via setting $\dot{a}=\dot{b}={\dot{\varphi }}_1={\dot{\varphi }}_2=0$ into equations (43)–(46), to get the following nonlinear algebraic system of equations

$$\begin{array}{l}{F}_1(a,b,{\varphi }_1,{\varphi }_2)=-\left[\mu +\frac{{\eta }_1{\eta }_3}{2({\omega }_1^2+{\omega }^2)}\right]a+\frac{1}{8}\left[{\alpha }_3-\frac{2{\beta }_1{\eta }_3}{{\omega }_1^2+{\omega }^2}-\frac{2{\beta }_3{\eta }_3^2{\omega }_1}{{({\omega }_1^2+{\omega }^2)}^2}\right]a^3+\frac{1}{8}[2{\alpha }_4-\frac{{\beta }_4{\eta }_4}{{\omega }_2^2+{\omega }^2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\frac{{\beta }_6{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}-\frac{2{\beta }_7{\eta }_3}{{\omega }_1^2+{\omega }^2}]a{b}^2+\frac{1}{8}[-{\alpha }_4+{\alpha }_7-\frac{{\beta }_4{\eta }_4}{{\omega }_2^2+{\omega }^2}+\frac{2{\beta }_5{\eta }_4^2{\omega }_2}{{({\omega }_2^2+{\omega }^2)}^2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{{\beta }_6{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}]a{b}^2\cos (2{\varphi }_2-2{\varphi }_1)-\frac{1}{8}[\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega +\frac{{\beta }_4{\eta }_4{\omega }_2}{\omega ({\omega }_2^2+{\omega }^2)}+\frac{{\beta }_5{\eta }_4^2({\omega }_2^2-{\omega }^2)}{\omega {({\omega }_2^2+{\omega }^2)}^2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\frac{{\beta }_6{\eta }_4\omega }{{\omega }_2^2+{\omega }^2}]a{b}^2\sin (2{\varphi }_2-2{\varphi }_1)+\frac{f}{2\omega }{(\omega +\sigma )}^2\sin \left({\phi }_1\right)=0\end{array}$$

(51)

$$\begin{array}{l}{F}_2(a,b,{\varphi }_1,{\varphi }_2)=-\left[\mu +\frac{{\eta }_2{\eta }_4}{2({\omega }_2^2+{\omega }^2)}\right]b+\frac{1}{8}\left[{\alpha }_3-\frac{2{\gamma }_1{\eta }_4}{{\omega }_2^2+{\omega }^2}-\frac{2{\gamma }_3{\eta }_4^2{\omega }_2}{{({\omega }_2^2+{\omega }^2)}^2}\right]b^3+\frac{1}{8}[2{\alpha }_4-\frac{{\gamma }_4{\eta }_3}{{\omega }_1^2+{\omega }^2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\frac{{\gamma }_6{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}-\frac{2{\gamma }_7{\eta }_4}{{\omega }_2^2+{\omega }^2}]a^{2}b+\frac{1}{8}[-{\alpha }_4+{\alpha }_7-\frac{{\gamma }_4{\eta }_3}{{\omega }_1^2+{\omega }^2}-\frac{2{\gamma }_5{\eta }_3^2{\omega }_1}{{({\omega }_1^2+{\omega }^2)}^2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{{\gamma }_6{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}]a^{2}b\cos (2{\varphi }_2-2{\varphi }_1)+\frac{1}{8}[\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega +\frac{{\gamma }_4{\eta }_3{\omega }_1}{\omega ({\omega }_1^2+{\omega }^2)}+\frac{{\gamma }_5{\eta }_3^2({\omega }_1^2-{\omega }^2)}{\omega {({\omega }_1^2+{\omega }^2)}^2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{{\gamma }_6{\eta }_3\omega }{{\omega }_1^2+{\omega }^2}]a^{2}b\sin (2{\varphi }_2-2{\varphi }_1)-\frac{f}{2\omega }{(\omega +\sigma )}^2\cos \left({\varphi }_2\right)=0\end{array}$$

(52)

$$\begin{array}{l}{F}_3(a,b,{\varphi }_1,{\varphi }_2)=\left[\sigma +\frac{{\eta }_1{\eta }_3{\omega }_1}{2\omega ({\omega }_1^2+{\omega }^2)}\right]+\frac{1}{8\omega }\left[3{\alpha }_1+{\alpha }_6{\omega }^2+\frac{2{\beta }_1{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}+\frac{{\beta }_3{\eta }_3^2({\omega }_1^2-{\omega }^2)}{{({\omega }_1^2+{\omega }^2)}^2}\right]a^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{1}{8\omega }\left[2{\alpha }_2+2{\alpha }_5{\omega }^2+\frac{{\beta }_4{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}-\frac{{\beta }_6{\eta }_4{\omega }^2}{{\omega }_2^2+{\omega }^2}+\frac{2{\beta }_7{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}\right]b^2+\frac{1}{8\omega }[{\alpha }_2-{\alpha }_5{\omega }^2+\frac{{\beta }_4{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{{\beta }_5{\eta }_4^2({\omega }_2^2-{\omega }^2)}{{({\omega }_2^2+{\omega }^2)}^2}+\frac{{\beta }_6{\eta }_4{\omega }^2}{{\omega }_2^2+{\omega }^2}]b^2\cos (2{\varphi }_2-2{\varphi }_1)+\frac{1}{8}[-{\alpha }_4+{\alpha }_7-\frac{{\beta }_4{\eta }_4}{{\omega }_2^2+{\omega }^2}-\frac{2{\beta }_5{\eta }_4^2{\omega }_2}{{({\omega }_2^2+{\omega }^2)}^2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{{\beta }_6{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}]b^2\sin (2{\varphi }_2-2{\varphi }_1)+\frac{f}{2a\omega }{(\omega +\sigma )}^2\cos \left({\varphi }_1\right)=0\end{array}$$

(53)

$$\begin{array}{l}{F}_4(a,b,{\varphi }_1,{\varphi }_2)=\left[\sigma +\frac{{\eta }_2{\eta }_4{\omega }_2}{2\omega ({\omega }_2^2+{\omega }^2)}\right]+\frac{1}{8\omega }\left[3{\alpha }_1+{\alpha }_6{\omega }^2+\frac{2{\gamma }_1{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}+\frac{{\gamma }_3{\eta }_4^2({\omega }_2^2-{\omega }^2)}{{({\omega }_2^2+{\omega }^2)}^2}\right]b^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{1}{8\omega }\left[2{\alpha }_2+2{\alpha }_5{\omega }^2+\frac{{\gamma }_4{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}-\frac{{\gamma }_6{\eta }_3{\omega }^2}{{\omega }_1^2+{\omega }^2}+\frac{2{\gamma }_7{\eta }_4{\omega }_2}{{\omega }_2^2+{\omega }^2}\right]a^2+\frac{1}{8\omega }[{\alpha }_2-{\alpha }_5{\omega }^2+\frac{{\gamma }_4{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{{\gamma }_5{\eta }_3^2({\omega }_1^2-{\omega }^2)}{{({\omega }_1^2+{\omega }^2)}^2}+\frac{{\gamma }_6{\eta }_3{\omega }^2}{{\omega }_1^2+{\omega }^2}]a^2\cos (2{\varphi }_2-2{\varphi }_1)-\frac{1}{8}[-{\alpha }_4+{\alpha }_7+\frac{{\gamma }_4{\eta }_3}{{\omega }_1^2+{\omega }^2}-\frac{2{\gamma }_5{\eta }_3^2{\omega }_1}{{({\omega }_1^2+{\omega }^2)}^2}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\frac{{\gamma }_6{\eta }_3{\omega }_1}{{\omega }_1^2+{\omega }^2}]a^2\sin (2{\varphi }_2-2{\varphi }_1)+\frac{f}{2b\omega }{(\omega +\sigma )}^2\sin \left({\varphi }_2\right)=0\end{array}$$

(54)

By solving the above nonlinear algebraic system (i.e., $F_j\left(a, b, {\varphi }_1, {\varphi }_2\right)=0, j=\mathrm{1,2,3,4}$ ) in terms of one of the parameters $\sigma , f, {\eta }_1, {\eta }_2, {\eta }_3, {\eta }_4, {\omega }_1, {\omega }_2, p,$ and $d$ , we can establish the various bifurcation diagrams as given in the following section. In addition, the stability of the steady-state solution can be explored by applying the Hartman–Grobman theorem. ⁴⁶ This can be done via letting ( $a_0, b_0, {\varphi }_{10}, {\varphi }_{20}$ ) be the solution of equations (51)–(54) and assuming that ( $a_1, b_1, {\varphi }_{11}, {\varphi }_{21}$ ) is a variation about that solution. According to these assumptions, we can write

$$a=a_0+a_1,\hspace{0.17em}\hspace{0.17em}b=b_0+b_1,\hspace{0.17em}\hspace{0.17em}{\varphi }_1={\varphi }_{10}+{\varphi }_{11},\hspace{0.17em}\hspace{0.17em}{\varphi }_2={\varphi }_{20}+{\varphi }_{21}$$

(55)

By differentiating equation (55), we have

$$\dot{a}={\dot{a}}_1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\dot{b}={\dot{b}}_1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\dot{\varphi }}_1={\dot{\varphi }}_{11},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\dot{\varphi }}_2={\dot{\varphi }}_{21}$$

(56)

Substituting equations (55) and (56) into equations (43)–(46), expanding the resulting equations for the small deviations $a_1, b_1, {\phi }_{11},$ and ${\phi }_{21}$ keeping the linear terms only, yield

$$\left(\begin{array}{c}{\dot{a}}_1\\ {\dot{b}}_1\\ {\dot{\varphi }}_{11}\\ {\dot{\varphi }}_{21}\end{array}\right)=\left(\begin{array}{cccc}\frac{\partial F_1}{\partial a_1}& \frac{\partial F_1}{\partial b_1}& \frac{\partial F_1}{\partial {\varphi }_{11}}& \frac{\partial F_1}{\partial {\varphi }_{21}}\\ \frac{\partial F_2}{\partial a_1}& \frac{\partial F_2}{\partial b_1}& \frac{\partial F_2}{\partial {\varphi }_{11}}& \frac{\partial F_2}{\partial {\varphi }_{21}}\\ \frac{\partial F_3}{\partial a_1}& \frac{\partial F_3}{\partial b_1}& \frac{\partial F_3}{\partial {\varphi }_{11}}& \frac{\partial F_3}{\partial {\varphi }_{21}}\\ \frac{\partial F_4}{\partial a_1}& \frac{\partial F_4}{\partial b_1}& \frac{\partial F_4}{\partial {\varphi }_{11}}& \frac{\partial F_4}{\partial {\varphi }_{21}}\end{array}\right)\left(\begin{array}{c}{a}_1\\ b_1\\ {\varphi }_{11}\\ {\varphi }_{21}\end{array}\right)$$

(57)

According to the Hartman–Grobman theorem ⁴⁶ , the linear dynamical system (57) is equivalent to the nonlinear system (43)–(46). Therefore, the stability of the solution of equations (51)–(54) can be investigated by examining the linear system’s eigenvalues given by equation (57).

Bifurcation diagrams and the performance of the integral resonant controller

This section investigates the effects of various controller parameters on rotor steady-state vibration amplitudes ( $a,b$ ) and phase angles ( ${\varphi }_1, {\varphi }_2$ ). Solving the derived nonlinear algebraic system equations (51)–(54) (numerically using the Newton–Raphson algorithm ⁴⁷ ), one can obtain the response curves as shown in Figure 2, Figure 3, Figure 4, etc. The solution’s stability of equations (51)–(54) can be explored by investigating the eigenvalues of the linear system given by equation (57), where the stable solution is plotted as a solid line and the unstable one is illustrated as a dotted line. In addition, to confirm that the derived nonlinear algebraic system (i.e. equations (51)–(54)) describes with high accuracy the steady-state vibration amplitudes of the nonlinear system governed by equations (17)–(20), the steady-state oscillation amplitudes of equations (17)–(20) are obtained numerically using Rung–Kutta algorithm. The numerical solution of equations (17)–(20) is plotted as a small circle when the bifurcation parameter is being swept forward and as a big dot when it is being swept backward. The following parameter values are used to simulate the considered system. Unless otherwise stated, the analytical and numerical findings were obtained by using the following values for the system parameters: $P=1.5, d=0.005, \alpha ={30}^o, {\eta }_1={\eta }_2= {\eta }_3={\eta }_4={\omega }_1={\omega }_2=0.0,$ $f=0.012$ , and $\mathrm{\Omega }=\omega +\sigma , \sigma =0$ [15].

Figure 2.

12-pole system $\sigma -$ response curve when $f=0.004, 0.008, 0.012,$ $P=1.5, d=0.005, {\eta }_1={\eta }_2= {\eta }_3={\eta }_4={\omega }_1={\omega }_2=0.0$ : (a) vibration amplitude in $X-$ direction, (b) vibration amplitude in $Y-$ direction, (c) the steady-state phase angles at $f=0.004$ , and (d) the steady-state phase angles at $f=0.012$ .

Figure 3.

12-pole system $\sigma -$ response curve at $P=1.35, f=0.012, d=0.005,$ and ${\eta }_1={\eta }_2= {\eta }_3={\eta }_4={\omega }_1={\omega }_2=0.0$ : (a) vibration amplitude in $X-$ direction, (b) vibration amplitude in $Y-$ direction, and (c) steady-state phase angles.

Figure 4.

12-pole system $\sigma -$ response curve at $P=1.2, f=0.012, d=0.005,$ and ${\eta }_1={\eta }_2= {\eta }_3={\eta }_4={\omega }_1={\omega }_2=0.0$ : (a) vibration amplitude in $X-$ direction, (b) vibration amplitude in $Y-$ direction, and (c) steady-state phase angles.

12-pole rotor system dynamics without integral resonant controller

The dynamical behavior of the 12-pole system without IRC control is discussed within this subsection via solving equations (51)–(54) when ${\eta }_1={\eta }_2= {\eta }_3={\eta }_4={\omega }_1={\omega }_2=0.0$ . Before proceeding further it should be remembered that the dimensionless parameters of the rotor system and PD-controller are defined before equations (17) so that $f=\frac{e}{c_0}, \mathrm{\Omega }=\frac{\psi }{\zeta }, P=\frac{c_0}{I_0}{k}_1,$ and $d=\frac{c_0\zeta }{I_0}{k}_2$ . Accordingly, one can deduce that $f$ represents the dimensionless form of the actual eccentricity $e$ , $\mathrm{\Omega }$ is the normalized form of the actual angular velocity $\psi$ , $P$ denotes the dimensionless form of the actual proportional gain $k_1$ , and $d$ is the dimensionless form of the actual derivative gain $k_2$ . In addition, the parameter $\sigma$ represents the closeness of $\mathrm{\Omega }$ to the system natural frequency $\omega$ as given in equation (39) (i.e., $\sigma =\mathrm{\Omega }-\omega$ ). Therefore, the parameter $\sigma$ is employed as a bifurcation parameter to describe the system oscillation amplitudes ( $a, b$ ) and phase angles ( ${\varphi }_1$ , ${\varphi }_2$ ) of the rotor system at the primary resonance case (i.e., $\mathrm{\Omega }\to \omega$ ) as shown in Figures 2–4 and so on.

Figure (2) illustrates the 12-pole system oscillation amplitudes and phase angles at different values of the eccentricity $f$ when $P=1.5$ and $d=0.005$ . Figures 2(a) and (b) show that the 12-pole system vibrations when $f=0.004, 0.008,$ and $0.012$ . The figures depict that system lateral vibration amplitudes are a monotonic growing function of $f$ , where the nonlinearities dominate the system motion when $f=0.008$ and $0.012$ . Figures 2(a) and (b) depict that the rotor system behaves as a hardening Duffing oscillator at large disc electricity with symmetric motion in $X$ and $Y$ directions. Based on equations (47) and (48), the rotor system can execute forward whirling motion (i.e., whirling counterclockwise) as long as ${\varphi }_2-{\varphi }_1>0$ , while the system exhibits backward whirling (i.e., whirling clockwise) if ${\varphi }_2-{\varphi }_1<0$ . In addition, the system oscillates along a straight line with a slope $\varphi$ if ${\varphi }_2={\varphi }_1=\varphi$ . Therefore, one can determine the whirling direction of the rotor system by plotting the phase angles ( ${\varphi }_1, {\varphi }_2$ ) against the targeted bifurcation parameter. Figures 2(c) and (d) illustrate the evolution of the system phase angles against the detuning parameter $\sigma$ when $f=0.004,$ and $0.012$ , respectively. Figures 2(c) and (d) show that the phase difference ${\varphi }_2-{\varphi }_1=\frac{\pi }{2}$ , which demonstrated that the rotor system performs circular forward whirling motion when $P=1.5$ and $d=0.005$ regardless of the magnitude of the eccentricity.

To confirm the precision of the findings shown in Figure 2, the rotor system’s steady-state temporal oscillations and whirling motion have been numerically simulated as shown in Figure 5 via solving the system dynamical model given by equations (17)–(20) using ODE45 MATLAB function according to Figure 2 when $f=0.012, \sigma =0.05$ (i.e., when $f=0.012, P=1.5, d=0.005, \mathrm{\Omega }=\omega +0.05,$ and ${\eta }_1={\eta }_2= {\eta }_3={\eta }_4={\omega }_1={\omega }_2=0.0$ ). It is clear from Figures 2(a) and (b) that the system has two stable solutions at $f=0.012$ and $\sigma =0.05$ . Therefore, equations (17)–(20) have been solved numerically using two trial initial conditions (i.e. $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=q_x\left(0\right)=q_y\left(0\right)=0$ and $x\left(0\right)=0.8, \dot{x}\left(0\right)=0.3, y\left(0\right)=-0.8, \dot{y}\left(0\right)=0.3, q_x\left(0\right)=q_y\left(0\right)=0$ ). By examining Figure 5, one can confirm that the rotor may vibrate by one of two stable periodic solutions according to the initial position. In addition, the system performs forward circular whirling motion because $a$ and $b$ have the same magnitude as well as the phase difference ${\varphi }_2-{\varphi }_1=\frac{\pi }{2}$ as is clear in Figure 2(d).

Figure 5.

The steady-state temporal vibrations of the rotor system according to Figure 2 when $\sigma =0.05, f=0.012$ at the initials $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=q_x\left(0\right)=q_y\left(0\right)=0$ and $x\left(0\right)=0.8, \dot{x}\left(0\right)=0.3, y\left(0\right)=-0.8, \dot{y}\left(0\right)=0.3, q_x\left(0\right)=q_y\left(0\right)=0$ : (a, b) the vibration in $X$ and $Y$ directions, (c) the whirling motions.

The effect of decreasing the proportional gain from $P=1.5$ as in Figure 2 to $P=1.35$ or $P=1.2$ is illustrated in Figures 3 and 4. The figures show the complex static bifurcation behaviors when decreasing the proportional gain $P$ beyond $1.5$ ., where the system may respond with mono-stale, bi-stable, and tri-stable solutions depending on the angular velocity $\mathrm{\Omega }=\omega +\sigma$ . In addition, Figures 3(c) and 4(c) show that the system can whirl backward or forward according to the value of the initial conditions at the same spinning speed. The rotor temporal oscillations and the whirling orbit according to Figure 3 at $\sigma =0.05$ have been illustrated in Figure 6 via solving equations (17)–(20) numerically utilizing the two initials $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=q_x\left(0\right)=q_y\left(0\right)=0$ and $x\left(0\right)=0.5, \dot{x}\left(0\right)=0.3, y\left(0\right)=\dot{y}\left(0\right)= q_x\left(0\right)=q_y\left(0\right)=0$ , where the figure confirms the elliptic forward or backward whirling motion of the system.

Figure 6.

The steady-state temporal vibrations of the rotor system according to Figure 4 when $\sigma =0.0$ five at the initial conditions $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=q_x\left(0\right)=q_y\left(0\right)=0$ and $x\left(0\right)=0.5, \dot{x}\left(0\right)=0.3, y\left(0\right)= \dot{y}\left(0\right)=q_x\left(0\right)=q_y\left(0\right)=0$ : (a, b) the vibration in $X$ and $Y$ directions, (c) the whirling motions.

Based on Figure 3, Figure 4, and Figure 6, one can demonstrate that the dynamical characteristics of the 12-pole system are particularly sensitive to the proportional gain, where for the large $P$ values (i.e., $P=1.5$ ), the system responds with simple bifurcation behaviors like the hardening spring Duffing oscillator but with strong oscillation amplitudes. On the other hand, at the small values of the proportional gain (i.e. $P=1.2$ ), the system responds with small oscillation amplitudes, but it may exhibit complex bifurcation behaviors. Accordingly, the IRC has been introduced along with the PD-controller to improve the dynamical characteristics of the considered system as discussed in The 12-pole rotor dynamics with asymmetric IRC controller and The 12-pole rotor dynamics with symmetric IRC controller.

Before going any further, it’s crucial to keep in mind that the IRC control parameters have been defined in their dimensionless form before equation (17) such that ${\eta }_1=\frac{c_0}{I_0}\left(1+2\cos \left(\alpha \right)\right)k_3, {\eta }_2=\frac{c_0}{I_0}\left(1+2\cos \left(\alpha \right)\right)k_4 , {\eta }_3=\frac{k_5}{\zeta }, {\eta }_4=\frac{k_6}{\zeta } , {\omega }_1=\frac{{\omega }_x}{\zeta },$ and ${\omega }_2=\frac{{\omega }_y}{\zeta }$ . According to these definitions and based on equations (6)–(8), one can simply deduce that ${\eta }_1$ and ${\eta }_2$ represent the dimensionless forms of the control signal gains ( $k_3$ and $k_4$ ) of the IRCs that are connected to the horizontal and vertical oscillation modes, respectively. In addition, ${\eta }_3$ and ${\eta }_4$ are the dimensionless forms of the feedback signal gains ( $k_5$ and $k_6$ ) of the IRCs. Moreover, ${\omega }_1$ and ${\omega }_2$ denote the dimensionless forms of the internal feedback signal gains ( ${\omega }_x$ and ${\omega }_y$ ) of the integral resonant controllers. Secondly, it is clear from the derived autonomous dynamical system given by equations (43)–(46) that the integration of the integral resonant controllers to the 12-pole system has modified the linear damping (i.e., $\mu$ ) to the equivalent damping ${\mu }_x=\mu +\frac{{\eta }_1{\eta }_3}{2({\omega }_1^2+{\omega }^2)}$ (in the horizontal direction) and ${\mu }_y=\mu +\frac{{\eta }_2{\eta }_4}{2({\omega }_2^2+{\omega }^2)}$ (in the vertical direction). In addition, the detuning parameter $\sigma$ has been modified to the equivalent detuning parameters ${\sigma }_x=\mu +\frac{{\eta }_1{\eta }_3{\omega }_1}{2\omega ({\omega }_1^2+{\omega }^2)}$ (in the horizontal direction) and ${\sigma }_y=\mu +\frac{{\eta }_2{\eta }_4{\omega }_2}{2\omega ({\omega }_2^2+{\omega }^2)}$ (in the vertical direction). Accordingly, ${\mu }_x$ is proportional to the multiplication of both control and feedback gains (i.e., ${\eta }_1\times {\eta }_3$ ) of the IRC that is connected to the horizontal oscillation mode only, and inversely proportional to the square of ${\omega }_1$ as well as the square of the system’s natural frequency ( $\omega$ ). In addition, ${\mu }_y$ is proportional to the multiplication of both control and feedback gains ( ${\eta }_2\times {\eta }_4$ ) of the IRC that is connected to the vertical oscillation mode, and inversely proportional to the square of both ${\omega }_2$ and $\omega$ . Therefore, the product of the control and feedback gains ${\eta }_1$ and ${\eta }_3$ as well as the internal loop feedback gain ${\omega }_1$ can be utilized to reshape the rotor dynamics in the horizontal direction via adjusting the values of ${\mu }_x$ and ${\sigma }_x$ , while the product of the control and feedback gains ${\eta }_2$ and ${\eta }_4$ as well as the internal loop feedback gain ${\omega }_2$ can be utilized to control the rotor dynamics in the vertical direction via controlling the values of ${\mu }_y$ and ${\sigma }_y$ .

Based on the above explanations, the influence of the different control parameters (i.e., ${\eta }_1, {\eta }_2, {\eta }_3, {\eta }_4 , {\omega }_1,$ and ${\omega }_2$ ) on the rotor dynamics has been investigated in two-stage of analyses. In the first stage, the rotor dynamics have been explored when activating the IRC controller that is connected to the $X$ direction only when turning off the other controller (i.e., ${\text{ η}}_2={\text{η}}_4={\text{ω}}_2=0$ ) in The 12-pole rotor dynamics with asymmetric IRC controller. In the second stage, the two controllers are treated as one controller by letting ${\eta }_1={\eta }_2={\eta }_c, {\eta }_3={\eta }_4={\eta }_f$ and ${\omega }_1={\omega }_2={\omega }_c$ as discussed in The 12-pole rotor dynamics with symmetric IRC controller.

12-pole rotor dynamics with asymmetric IRC controller

The 12-pole rotor dynamics have been investigated within this section when activating the IRC that is coupled to the horizontal direction with turning off the controller in the vertical direction (i.e., ${\eta }_2={\eta }_4={\omega }_2=0$ ). Figure 7 shows the $\sigma -$ response curve of the 12-pole rotor system at three different values of ${\eta }_1\times {\eta }_3$ (i.e., at ${\eta }_1{\eta }_3=0.04, 0.08, 0.16$ ) when ${\omega }_1=1.0$ . It is clear from Figures 7(a) and (c), and (e) that the maximum oscillation amplitude of the rotor system in the horizontal direction is $a\approx 0.55$ when ${\eta }_1{\eta }_3=0.04$ , $a\approx 0.3$ when ${\eta }_1{\eta }_3=0.08$ , and $a\approx 0.2$ when ${\eta }_1{\eta }_3=0.16$ . Accordingly, one can deduce that the system oscillation amplitude in the horizontal direction is a monotonic decreasing function of the product of both the control and feedback signal gains ${\eta }_1{\eta }_3$ , which agrees with the mathematical explanation given above, where increasing of ${\eta }_1{\eta }_3$ , increases the system equivalent damping coefficient in the horizontal direction ( ${\mu }_x=\mu +\frac{{\eta }_1{\eta }_3}{2({\omega }_1^2+{\omega }^2)}$ ) without affecting the equivalent damping in the vertical direction ${\mu }_y=\mu +\frac{{\eta }_2{\eta }_4}{2({\omega }_2^2+{\omega }^2)}$ . Despite increasing of ${\eta }_1{\eta }_3$ decreases the system oscillation amplitude in $X-$ direction, it cannot eliminate the motion bifurcation, where the bistable solution region and nonlinearities still dominate the system response due to the natural nonlinear coupling between the horizontal and vertical oscillation modes. Figures 7(b) and (d), and (f) also show that activating one of the controllers can change the whirling direction of the rotor system either forward or backward depending on the spinning speed ( $\mathrm{\Omega }=\omega +\sigma )$ and the values of ${\eta }_1{\eta }_3$ .

Figure 7.

12-pole system $\sigma -$ response curve at $f=0.012, P=1.5,d=0.005,$ ${\omega }_1=1,$ and ${\eta }_2{\eta }_4={\omega }_2=0.0$ : (a, b) ${\eta }_1{\eta }_3=0.04$ , (c, d) ${\eta }_1{\eta }_3=0.08$ , and (e, f) ${\eta }_1{\eta }_3=0.16$

Numerical simulation for the temporal oscillations and whirling motion of the system according to Figure 7 at ${\eta }_1{\eta }_3=0.04, 0.08, 0.16$ , and $\sigma =0.0$ is illustrated in Figure 8 via solving equations (17)–(20) using ODE45. Figures 8(a) to (c) show the temporal oscillation and the whirling motion of the system and controller at ${\eta }_1{\eta }_3=0.04$ , Figures 8(d) to (f) illustrate the lateral vibrations and the corresponding whirling motion of the system and controller when ${\eta }_1{\eta }_3=0.08$ , and Figures 8(g) to (i) depict the temporal vibration and the corresponding whirling motion of the system and controller when ${\eta }_1{\eta }_3=0.16$ and $\sigma =0.0$ . Comparing Figures 8(a), (d), and (g), it is clear the temporal oscillations of the 12-pole rotor system decrease as ${\eta }_1{\eta }_3$ increases. In addition, Figures 8(c), (f), and (i) confirm that the system can perform elliptic forward or backward whirling motion depending on the magnitude of ${\eta }_1{\eta }_3$ when fixing $\sigma =0.0$ .

Figure 8.

The steady-state temporal oscillations of the rotor system and IRC according to Figure 7 when $\sigma =0.0$ : (a, b, c) ${\eta }_1{\eta }_3=0.04$ , (d, e, f) ${\eta }_1{\eta }_3=0.08$ , and (g, h, i) ${\eta }_1{\eta }_3=0.16$ .

The effects of the internal feedback gain ${\omega }_1$ on the oscillation amplitudes of the considered system have been illustrated in Figure 9 when ${\eta }_1{\eta }_3=0.16, {\eta }_2={\eta }_4={\omega }_2=0$ Figures 9(a) and (b) shows the system oscillation amplitudes and phase angles versus the angular velocity when ${\omega }_1=0.5$ , while the rotor system lateral oscillations and the phase angles are illustrated in Figures 9(c) and (d) when ${\omega }_1=3.0$ . In addition, the system’s steady-state lateral vibration and the phase angles have been shown in Figures 9(e) and (f) when ${\omega }_1=6.0$ . Comparing Figures 9(a) and (c), and (e), one can deduce that the increase in the internal feedback gain ${\omega }_1$ from $0.5$ to $3.0$ then to $6.0$ , increases the rotor oscillation amplitude $a$ and shifts the whole response curve to the left. In addition, Figures 9(b) and (f) illustrate that the rotor can exhibit forward whirling only regardless of the initial conditions and the system angular velocity. However, Figure 9(d) shows that the system can perform one of three whirling modes according to the rotor angular velocity, which are: (1) forward whirling only, (2) backward whirling only, and (3) forward or backward whirling.

Figure 9.

12-pole system $\sigma -$ response curve at $f=0.012, P=1.5,d=0.005,$ ${\eta }_1{\eta }_3=0.16,$ and ${\eta }_2{\eta }_4={\omega }_2=0.0$ : (a, b) ${\omega }_1=0.5$ , (c, d) ${\omega }_1=3.0$ , and (e, f) ${\omega }_1=6.0$

Numerical simulations for the system’s original equations (17)–(20) according to Figures 9(c) and (d) when $\sigma =0.015$ and $0.035$ have been shown in Figure 10. Figures 9(c) and (d) illustrate that the rotor system has a monostable solution only where ${\varphi }_1>{\varphi }_2$ at $\sigma =0.015$ , while at $\sigma =0.035$ , the system has bi-stable solutions where ${\varphi }_1>{\varphi }_2$ with one of these solutions and ${\varphi }_2>{\varphi }_1$ with the other solution. Accordingly, it is expected that the system may whirl backward with a unique solution regardless of the intimal conditions when $\sigma =0.015$ . On the other hand, it is expected that the system may whirl forward or backward at $\sigma =0.035$ depending on the initial conditions. Therefore, the two initials $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=q_x\left(0\right)=q_y\left(0\right)=0$ and $x\left(0\right)=0.3, \dot{x}\left(0\right)=0.0, y\left(0\right)=0.6, \dot{y}\left(0\right)=0.2, q_x\left(0\right)=q_y\left(0\right)=0$ are have been used to solve equations (17)–(20) numerically when $\mathrm{\Omega }=\omega +0.015$ and $\mathrm{\Omega }=\omega +0.035$ as shown in Figures 10(a) to (f), respectively. One can note from Figures 10(a) to (c) that the system executes backward whirling motion regardless of the initial conditions, while Figures 10(d) to (f) show that the 12-pole system executes either backward or forward whirling according to the initial position of the rotor.

Figure 10.

The temporal oscillations of the 12-pole rotor and IRC according to Figures 9(c) and (d) when ${\omega }_1=3.0$ and ${\omega }_2=0.0$ at the initial conditions $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=q_x\left(0\right)=q_y\left(0\right)=0$ and $x\left(0\right)=0.3, \dot{x}\left(0\right)=0.0, y\left(0\right)=0.6, \dot{y}\left(0\right)=0.2, q_x\left(0\right)=q_y\left(0\right)=0$ : (a, b, c) $\sigma =0.015$ , and (d, e, f) $\sigma =0.035$ .

The influence of the proportional gain $P$ on the control performance of the IRC has been investigated as shown in Figure 11, where the $\sigma -$ response curve for both the vibrations amplitudes and the phase angles are obtained when $P=1.5, 1.4,$ and $1.3$ . By examining Figures 11(a) and (c) and (e), one can find that the increase of the $P$ gain increases the system oscillation amplitude in $X$ direction.

Figure 11.

12-pole system $\sigma -$ response curve at $f=0.012,d=0.005,$ ${\eta }_1{\eta }_3=0.08, {\omega }_1=1.0$ , and ${\eta }_2{\eta }_4={\omega }_2=0.0$ : (a, b) $P=1.3$ , (c, d) $P=1.4$ , and (e, f) $P=1.5$

Based on Figures 9 and 11, it is clear that the rotor steady-state lateral oscillation in the horizontal direction (i.e., $a$ ) is a monotonic growing function of both the internal feedback gain ( ${\omega }_1$ ) and the proportional gain $\left(P\right)$ . This mechanism can be explained according to the introduced definition in for the equivalent damping coefficient ${\mu }_x=\mu +\frac{{\eta }_1{\eta }_3}{2({\omega }_1^2+{\omega }^2)}$ and ${\mu }_y=\mu +\frac{{\eta }_2{\eta }_4}{2({\omega }_2^2+{\omega }^2)}$ . First, remember that the system’s natural frequency $\omega$ is a monotonic growing function of the proportional gain $P$ , where $\omega =\sqrt{P(2cos\left(\alpha \right)+1)-3}$ . Therefore, increasing the system’s natural frequency (via increasing $P$ ) or/and the internal feedback gain ${\omega }_1$ , decreases the system equivalent damping coefficient ${\mu }_x$ , which ultimately increases the oscillation amplitude in $X$ direction.

The bifurcation behaviors of the 12-pole system have been explored in Figure (12) utilizing $f$ as the bifurcation parameter at two different values of $\text{σ}$ when ${\eta }_1{\eta }_3=0.016, {\omega }_1=1.0, {\eta }_2{\eta }_4={\omega }_2=0.0$ . Figures 12(a) and (c) illustrate that the rotor exhibits stable periodic oscillations as long as $0\le f\le 0.15$ . However, the system exhibits whirl backward or forward depending on the magnitude of both $f$ and $\sigma$ . Figures 12(a) and (b) illustrate that the system vibrates periodically with forward whirling motion along the eccentricity axis except for the interval ( $0.0113<f<0.0216$ ), at which the rotor may execute backward whirling. Also, Figures 12(c) and (d) illustrate that the rotor system performs backward whirling as long as $0.0225<f<0.0368$ , otherwise, the system executes forward whirling.

Figure 12.

12-pole system $f-$ response curve $P=1.5,d=0.005,$ ${\eta }_1{\eta }_3=0.16, {\omega }_1=1.0$ and ${\eta }_2{\eta }_4={\omega }_2=0.0$ : (a, b) $\sigma =0.0$ , and (c, d) $\sigma =0.05$

According to Figures 7–9, Figures 11 and 12, it is noticeable that activating one of the applied integral resonant controllers with turning off the other is resulting in the breakdown of the bifurcation symmetry of the considered system, which ultimately leads to different nonlinear phenomena such as (1) the appearance of the localized and nonlocalized oscillation and (2) the elliptic forward or/and backward whirling. However, the rotor system cannot lose its stability due to the asymmetry caused by the activation of one of the controllers only. Finally, it is important to note that if the values of the parameters ${\eta }_1, {\eta }_3,$ and ${\omega }_1$ have been interchanged with the parameters ${\eta }_2,{\eta }_4,$ and ${\omega }_2$ , respectively, one can obtain the same bifurcation diagrams that are given in Figures 7–9, Figures 11 and 12, but the bifurcation behaviors in $X$ and $Y$ directions will be interchanged.

12-pole rotor dynamics with symmetric IRC controller

The influence of the symmetric IRC on the bifurcation behaviors of the considered rotor system is explored within this section. The coupled integral resonant controllers are considered as one controller via letting the control feedback gains ${\eta }_1={\eta }_2={\eta }_c$ , the feedback gains ${\eta }_3={\eta }_4={\eta }_f$ , and the internal feedback gain ${\omega }_1={\omega }_2={\omega }_c$ . According to equations (43)–(46), one can note that the coupling of the symmetric IRC to the 12-pole system modifies the system damping $\mu$ to the equivalent coefficients ${\mu }_{Eq}=\mu +\frac{{\eta }_c{\eta }_f}{2({\omega }_c^2+{\omega }^2)}$ .

Based on the above assumption, the effects of the IRC parameters ( ${\eta }_c, {\eta }_f,$ and ${\omega }_c$ ) on the rotor system dynamics are explored through Figure 13 illustrates the system $\sigma -$ response curve at various values of ${\eta }_c\times {\eta }_f$ (i.e., when ${\eta }_c{\eta }_f=0.0, 0.04,$ and $0.16$ ). The figure depicts that the 12-pole system lateral oscillation ( $a, b$ ) are decreasing functions of ${\eta }_c{\eta }_f$ . In addition, the figure shows the symmetric bifurcation behaviors in both the $X$ and $Y$ directions. Moreover, the figure shows that increasing the gain ${\eta }_c{\eta }_f$ to $0.16$ can eliminate the system vibration close to zero. The influence of increasing the internal feedback gain ${\omega }_c$ on the rotor lateral vibration is shown in Figure 14. The figure shows the lateral amplitudes ( $a$ and $b$ ) at ${\omega }_c=1.0, 3.0,$ and $10.0$ when $f=0.012$ and ${\eta }_c{\eta }_f=0.16$ . The figure confirms that the 12-pole system behaves like the linear systems at ${\omega }_c=1.0$ and $3.0$ , but the nonlinearities may dominate the response when ${\omega }_c$ has increased to $10.0$ . Based on Figures 13 and 14, one can deduce that the oscillation amplitudes of the rotor are decreasing functions of the gain ${\eta }_c{\eta }_f$ , and monotonic growing function of the internal feedback gain ${\omega }_c$ . It is worth mentioning that the deuced results from Figures 13 and 14 agree with the introduced definition of the new damping coefficient ${\mu }_{Eq}=\mu +\frac{{\eta }_c{\eta }_f}{2({\omega }_c^2+{\omega }^2)}$ , where ${\mu }_{Eq}$ proportional to ${\eta }_c{\eta }_f$ , and inversely proportional to ${\omega }_c^2$ . Accordingly, increasing ${\eta }_c{\eta }_f$ and/or decreasing ${\omega }_c$ , increases the equivalent linear damping of the rotor system, which ultimately decreases the system lateral vibration in both $X$ and $Y$ directions.

Figure 13.

12-pole system $\sigma -$ response curve at $f=0.012,P=1.5, d=0.005,$ and ${\omega }_c=1.0$ when and ${\eta }_c{\eta }_f=0.0, 0.04$ and $0.16$ : (a) vibration amplitude in $X-$ direction, and (b) vibration amplitude in $Y-$ direction.

Figure 14.

12-pole system $\sigma -$ response curve at $f=0.012,P=1.5, d=0.005,$ and ${\eta }_c{\eta }_f=0.16$ when ${\omega }_c=1.0, 3.0$ and $10.0$ : (a) vibration amplitude in $X-$ direction, and (b) vibration amplitude in $Y-$ direction.

The control efficiency of the symmetric IRC in reducing the undesired vibrations of the 12-pole system is investigated when $P=1.2, 1.4,$ and $1.6$ as shown in Figure 15. By examining Figure 15, one can demonstrate that the increase of $P$ , increases the vibration amplitudes in both $X$ and $Y$ directions. The main cause that made the oscillation amplitudes increase with increasing the proportional gain can be explained simply with the help of equivalent damping definition ${\mu }_{Eq}=\mu +\frac{{\eta }_c{\eta }_f}{2({\omega }_c^2+{\omega }^2)}$ . It is clear that ${\mu }_{Eq}$ is inversely proportional to ${\omega }^2$ , where ${\omega }^2=P(2cos\left(\alpha \right)+1)-3$ . Therefore, increasing the proportional gain decreases the system equivalent damping coefficient, which in turn increases the lateral vibration amplitudes. Figure 16 depicts the time-responses of both the 12-poles system and the coupled IRCs according to Figure 15 at $\sigma =0.0$ (i.e., when ${\eta }_c{\eta }_f=0.16$ , ${\omega }_c=1.0$ , $f=0.012, d=0.005$ ) when decreasing the proportional gain from $P=1.6$ to $P=1.4$ at $t=300$ , and then decreasing $P$ again from $P=1.4$ to $P=1.2$ at $t=600$ . One can note from the figure that the system oscillation amplitudes decrease with the decreasing proportional gain which agrees exactly with Figure 15 at $\sigma =0.0$ . In addition, the figure illustrates the short transient time that the rotor takes to reach its steady state.

Figure 15.

12-pole system $\sigma -$ response curve at $f=0.012, d=0.005,{\eta }_c{\eta }_f=0.16$ and ${\omega }_c=1.0$ when $P=1.2, 1.4,$ and $1.6$ : (a) vibration amplitude in $X-$ direction, and (b) vibration amplitude in $Y-$ direction.

Figure 16.

The temporal vibrations of the rotor system and IRC according to Figure 15 when $\sigma =0.0$ regardless of the initial conditions when $P=1.6, 1.4,$ and $1.2$ : (a, b) the rotor temporal oscillations in $X$ and $Y$ directions, and (c, d) the temporal oscillations of the coupled integral resonant controllers.

The system eccentricity response curve before and after control are compared as shown in Figure 17 and Figure 18 when $\sigma =0.0$ and $\sigma =0.05$ , respectively. It is clear from the figure that the coupling of the symmetric IRC to the rotor system has suppressed the strong lateral vibrations, and eliminated the bifurcation behaviors, regardless of the eccentricity and angular velocity magnitudes. In addition, the nonlinear relation between the oscillation amplitudes and the eccentricity before control has been turned to the linear form when the controller is activated. The system time-response is simulated in Figure 19 according to Figure 18 when $f=0.03$ (i.e., $\sigma =0.05, {\eta }_c{\eta }_f=0.16, {\omega }_c=1.0$ ) before and after activating the integral resonant controller. The figure simulates the evolution of the instantaneous oscillations of the rotor from $t=0$ to $t=500$ while the controller is turned off. At the time instant $t=500$ , the IRC controller has turned on to mitigate the system’s instantaneous oscillation. It is clear from Figures 19(a) and (c) that the steady-state oscillation amplitudes in $X$ and $Y$ directions are $a=b\approx 0.8$ before control and $a=b\approx 0.2$ after control. Comparing these results (i.e., $a=b\approx 0.8$ before control and $a=b\approx 0.2$ after control) with Figure 18 when $f=0.03$ , one can demonstrate the excellent correspondence between analytical and numerical simulations.

Figure 17.

12-pole system $f-$ response curve $=1.5,d=0.005,$ ${\eta }_c{\eta }_f=0.16, {\omega }_c=1.0$ , and $\sigma =0.0$ : (a) vibration amplitude in $X-$ direction, (b) vibration amplitude in $Y-$ direction., and (c, d) the steady-state phase angles.

Figure 18.

12-pole system $f-$ response curve $=1.5,d=0.005,$ ${\eta }_c{\eta }_f=0.16, {\omega }_c=1.0$ , and $\sigma =0.05$ : (a) vibration amplitude in $X-$ direction, (b) vibration amplitude in $Y-$ \direction, and (c, d) the steady-state phase angles.

Figure 19.

The temporal vibrations of the rotor system before and after control according to Figure 18 when $\sigma =0.05$ and $f=0.03$ at zero initial conditions: (a, b, c) the rotor temporal oscillations in $X$ and $Y$ directions before and after turning the controller on, (b) the system phase trajectory before control, (d) the system phase trajectory after turning the controller on, (e, g) the IRC temporal oscillations before and after turning the control on, (f) the IRC phase trajectory before turning the control on, and (h) the IRC phase trajectory after turning the control on.

Conclusions

In this study, the integral resonant controller is introduced along with PD-controller as a modified control methodology to suppress the lateral vibrations of the 12-pole system. The nonlinear dynamical model of the controlled system has been derived as a two-degree-of-freedom system coupled with two first-order differential equations. The obtained nonlinear model has been transformed to its dimensionless form to generalize the analyses and discussions. Then, the nonlinear algebraic system that governs the rotor steady-state oscillation amplitudes and the phase angles has been extracted using perturbation analysis. The influence of the different control parameters on the steady-state oscillation amplitudes has been explored through many bifurcation diagrams. In addition, the effect of failing one of the applied integral resonant controllers on the rotor dynamics has been also discussed. Moreover, the transient oscillations of the controlled system have been simulated numerically. According to the introduced analysis, we can conclude the following remarks:

(1)Integrating symmetric integral resonant controller to the rotor system has modified the linear damping coefficient $\mu$ to ${\mu }_{Eq}=\mu +\frac{{\eta }_c{\eta }_f}{2({\omega }_c^2+{\omega }^2)}$ , where $\mu =\hspace{1em}\left(\cos \left(\alpha \right)+\frac{1}{2}\right)d$ and $\omega =\sqrt{2Pcos\left(\alpha \right)+P-3}$ .

(2)The system equivalent damping coefficient ${\mu }_{Eq}$ is proportional to both the derivative control gain ( $d$ ) and the multiplication of both the control and feedback gains ( ${\eta }_c{\eta }_f$ ) of the symmetric integral resonant controller.

(3)The system equivalent damping coefficient ${\mu }_{Eq}$ is inversely proportional to both the position gain ( $P$ ) and the square of the internal loop feedback gain ( ${\omega }_{\text{c}}^2$ ).

(4)According to remarks (1), (2), and (3), to get the best control performance for the combined IRC and PD controllers, the control parameters $P,d, {\eta }_c, {\eta }_f$ , and ${\omega }_c$ should be selected in a way that maximizes the function ${\mu }_{Eq}=\left(\cos \left(\alpha \right)+\frac{1}{2}\right)d+\frac{{\eta }_c{\eta }_f}{2({\omega }_c^2+2Pcos\left(\alpha \right)+P-3)}$ .

(5)The proposed symmetric resonant controller can efficiently suppress the lateral vibrations and force the nonlinear 12-pole system to perform forward circular whirling motion like the linear systems.

(6)The failure of one of the coupled integral resonant controllers may cause asymmetric oscillation for the rotor system, but the system does not lose its stability.

(7)Activating one of the applied integral resonant controllers while turning off the other one results in the breakdown of the bifurcation symmetry of the controlled system, which ultimately leads to different nonlinear phenomena such as the appearance of the localized and nonlocalized oscillation and the elliptic forward or/and backward whirling motions.