Integral Resonant Controller to Suppress the Nonlinear Oscillations of a Two-Degree-of-Freedom Rotor Active Magnetic Bearing System

Abstract

Within this article, the nonlinear vibration control of the rotor active magnetic bearings system is tackled utilizing the integral resonant controller for the first time. Two integral resonant controllers are proposed to mitigate the system lateral oscillations in the horizontal and vertical directions. Based on the suggested control technique, the whole system dynamical model is derived as a two-degree-of-freedom nonlinear system (i.e., rotor system) coupled linearly to two first-order filters (i.e., the integral resonant controllers). The nonlinear autonomous system that governs the oscillation amplitudes of the controlled system as a function of the control parameters is extracted by applying perturbation analysis. The obtained autonomous system showed that the linear damping coefficients of the rotor system are functions of the control gains, feedback gains, and internal loop feedback gains of the coupled controller. Accordingly, the sensitivity of the rotor oscillation amplitudes to the different control parameters is explored. The stability margins and the optimal control gains are reported via plotting the different stability charts in two-dimensional space. The main acquired results demonstrated that the vibration suppression efficiency of the proposed controller is proportional to the product of both the control and feedback signal gains, and inversely proportional to the square of the internal loop feedback gains. In addition, the analytical investigations confirmed that the proposed integral resonant control method can force the rotor system to respond as a linear one with a single periodic attractor when the control parameters are designed properly. Finally, numerical simulations are performed that have illustrated the excellent correspondence with the obtained analytical results.

1. Introduction

The nonlinear vibration control of the Active Magnetic Bearings System (AMBS) is the main subject of researchers and engineers worldwide due to its huge advantages over the conventional bearings system. The working principle of the magnetic bearings system relies on generating controllable magnetic forces to suspend the rotating shaft in its hovering position without physical contact with the stator. This working mechanism has earned the active magnetic bearings system many benefits over the conventional one such as (1) frictionless operation, (2) no necessity for lubrication, (3) less maintenance, (4) high-speed operation, (5) long working time, (6) clean environment, etc. Therefore, many researchers have investigated the dynamical behaviors of this system, as well as proposed different control strategies to enhance its dynamical stability. Ji et al. [1] studied the nonlinear behaviors of the 4-pole AMBS utilizing the conventional position-velocity feedback control technique. The author investigated the system dynamics in two oscillation modes. Firstly, he explored the system as an autonomous one via neglection of the rotor eccentricity, where the obtained results showed that the AMBS can exhibit saddle-node, saddle-connection, and Hopf bifurcations. Secondly, the author included the effect of the rotor eccentricity in the studied model, where the analysis illustrated that the nonautonomous system may have homoclinic orbit and transversal intersection. Saeed et al. [2,3] investigated the 6-pole AMBS with two different control methods, where they explored the vibration suppression efficiency of a proportional-derivative controller designed based on the Cartesian displacements and Cartesian velocities of the rotor system in [2]. In addition, the performance of a proportional-derivative controller implemented based on the radial displacements and radial velocities has been investigated in [3]. Based on the analysis introduced, they reported that the Cartesian control method has vibration suppression efficiency higher than the radial one. On the other hand, the radial control technique exhibited more stability features than the Cartesian method. The vibration control of the 8-pole AMBS with constant stiffness is investigated in [4,5,6,7,8]. Ji et al. [4,5] studied the position-velocity controller efficiency in suppressing the oscillatory behaviors of an 8-pole system at both the primary and superharmonic resonances. Yang et al. [6] discussed the oscillation modes of the 8-pole AMBS using the energy and phase-difference technique, where the obtained analysis demonstrated that the rotor system can perform either elliptic or quasiperiodic oscillations. Saeed et al. [7] introduced a combination of both the linear and nonlinear position-velocity controllers to enhance the dynamical behaviors of the 8-pole AMBS. They concluded that the integration of a cubic position controller can modify the system oscillatory behavior to respond as a hardening or softening spring oscillator. Saeed et al. [8] investigated the effect of rub-impact force on motion bifurcation of the 8-pole AMBS. The authors explored the nature of system motion utilizing both the impact stiffness and the dynamic friction coefficients as the bifurcation parameters. They reported that the system can oscillate with one of three vibration modes, which are (a) full annular rub, (b) period-n rub-impact, and (c) and quasiperiodic rub-impact. The efficiency of a time-varying proportional-derivative controller in mitigating the nonlinear vibrations of the 8-pole AMBS is introduced by Zhang et al. [9,10,11,12,13,14]. Based on the investigations introduced, they reported that the system exhibits Shilnikov multipulse chaotic motion. Furthermore, the system can oscillate with periodic-n, quasiperiodic, and chaotic motion. El-Shourbagy et al. [15] explored the oscillatory motion of the 12-pole AMBS for the first time utilizing the conventional PD-controller, where the authors reported that the proportional gain can play an important role in reshaping the system dynamics. Different control methodologies to control the 16-pole AMBS have been introduced [16,17,18,19,20,21]. Saeed and Kandil [16] utilized the PD-controller with two configuration methods to suppress the nonlinear vibrations of the 16-pole AMBS. The authors implemented the proposed control based on both the Cartesian and radial displacements and velocities of the rotor system. They found that the Cartesian control can suppress vibrations better than the radial control technique, while the radial control has higher stability characteristics than the Cartesian one. Zhang et al. [17,18,19,20,21] investigated the bifurcation behaviors of the 16-pole system utilizing time time-varying PD-controller. The reliability and durability features of the active magnetic bearing system made many of the researchers use it not only as a replacement for the conventional bearings system but also as an active actuator with advanced control algorithms to control the lateral vibrations in different rotating machinery [22,23,24,25,26,27].

The Integral Resonant Controller (IRC) proved its capability in suppressing the nonlinear oscillations and eliminating the catastrophic bifurcation behaviors for different categories of the dynamical systems [28,29,30,31,32,33,34]. Diaz et al. [28] control the vibrations of a lightweight structure utilizing IRC. Al-Mamun et al. [29] investigated the resonant vibrations control of a piezoelectric micro-actuator utilized the integral resonant controller. The authors concluded that the IRC is the best compared to the conventional control method and the use of a notch filter. Omidi and Mahmoodi [30], and MacLean and Sumeet [31] introduced modified versions of the IRC to control the transversal oscillation of a smart structure, where the obtained results demonstrated the efficiency of the applied controller. Omidi and Mahmoodi [32,33] integrated the IRC along with the Positive Position Feedback Controller (PPFC) to mitigate the primary resonance of a nonlinear flexible structure. The authors concluded that the IRC has high efficiency in suppressing the response peaks due to the PPFC. Recently, Saeed et al. [34] employed the nonlinear integral resonant controller with time delays to eliminate the principal parametric oscillation of a vertically suspended cantilever beam system. The authors showed that the introduced control method can mitigate the system oscillation amplitude to zero. Despite the high efficiency of the IRC in suppressing the nonlinear oscillations of different dynamical systems [28,29,30,31,32,33,34], it is not applied before to control the lateral vibrations of the magnetic bearing system [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27].

Within this work, the integral resonant controller has been introduced within this article as a new control strategy to mitigate the nonlinear oscillations of the 8-pole rotor active magnetic bearings system for the first time. Two integral resonant controllers are coupled to the 8-pole rotor system to suppress its horizontal and vertical oscillations. According to the proposed control method, the system equations of motion are derived relying on the principle of classical mechanics. The perturbation analysis is utilized to extract the controlled system amplitude-phase modulating equations. Then, the different response curves and stability charts are plotted. The influence of the different control gains on the system lateral vibrations is explored, and the system stability margins are reported. The main obtained results demonstrated that the linear damping coefficients of the rotor system are proportional to the product of both the control and feedback signal gains (i.e., ${\eta }_1 {\eta }_3$ and ${\eta }_2 {\eta }_4$), and inversely proportional to the square of the internal loop feedback gains (i.e., ${\lambda }_1^2$ and ${\lambda }_2^2$) of the coupled IRC controllers. In addition, it is reported that the IRC can force the rotor system to vibrate as a linear one with negligible oscillation amplitudes if the control parameters have been tuned properly.

2. Mathematical Modelling

This section is dedicated to obtaining the mathematical model of the eight-pole system that is controlled by both proportional-derivative and integral resonant controllers. The rotor system is modeled as shown in Figure 1 as a two-degree-of-freedom system that performs instantaneous displacements $x\left(t\right)$ and $y\left(t\right)$ in $X$ and $Y$ directions when spinning with angular speed $\upsilon$ about its axes. In addition, eight identical electromagnetic poles are utilized as an actuator to control the system oscillations via producing control force based on a predefined control law. Accordingly, the system equations of motions can be written as follows [35,36]:

$$m\ddot{x}(t)=me{\upsilon }^2\cos (\upsilon t)+F_{XC}$$

(1)

$$m\ddot{y}(t)=me{\upsilon }^2\sin (\upsilon t)+F_{YC}$$

(2)

where $m$ is the system mass, $x\left(t\right)$ and $y\left(t\right)$ are the temporal oscillations of the system in $X$ and $Y$ directions, $\upsilon$ denotes the system angular speed, $e$ is the eccentricity of the rotor, $F_{XC}$, $F_{YC}$ represent the net control forces that are generated via the eight electromagnetic poles to mitigate the rotor oscillations in $X$ and $Y$ directions. According to the electromagnetic theory [36], the electromagnetic attractive forces $F_j, \left(j=1,2,3,\dots ,8\right)$ of the eight poles shown in Figure 1B can be expressed as follows:

$$F_j=\Delta \frac{I_j^2}{h_j^2}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,2,\cdots ,8$$

(3)

where $\mathsf{\Delta }$ is constant such that $\mathsf{\Delta }=\frac{1}{4}{\mu }_0{N}^{2}A \cos \left(\theta \right)$ (${\mu }_0$ denotes the air-gap permeability, $N$ is the coil turn number of each pole, $A \cos \left(\theta \right)$ is the effective cross-sectional area), $I_j$ denotes the electrical current in ampere of the $j$^th pole, $h_j$ represent the instantaneous air gap size as shown in Figure 1B. Let $s_0$ be the nominal air gap of the rotor system, and $2\alpha$ be the angle between every two consecutive poles as shown in Figure 1B. So, for the small deviations $x\left(t\right)$ and $y\left(t\right)$ of the rotor, the instantaneous air gap size $h_j$ can be expressed as follows:

$$\begin{array}{l}{h}_j=s_0\pm x\sin (\alpha )\mp y\cos (\alpha ),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,5\\ h_j=s_0\pm x\sin (\alpha )\pm y\cos (\alpha ),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=4,8\\ h_j=s_0\pm x\cos (\alpha )\mp y\sin (\alpha ),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=2,6\\ h_j=s_0\pm x\cos (\alpha )\pm y\sin (\alpha ),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=3,7\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\end{array}\}$$

(4)

According to the schematic diagram shown in Figure 1C, the eight poles current ($I_1, I_2, \dots , I_8$) can be proposed such that:

$$\begin{array}{l}{I}_1(t)=I_8(t)=I_0-i_y(t),\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}{I}_2(t)=I_3(t)=I_0+i_x(t),\\ I_4(t)=I_5(t)=I_0+i_y(t),\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}{I}_6(t)=I_7(t)=I_0-i_x(t).\hspace{0.17em}\end{array}\}$$

(5)

where $I_0$ is constant current, and $i_x\left(t\right)$ and $i_y\left(t\right)$ are the control currents that depend on the proposed control law. Within this work, a combination of both the linear proportional-derivative controller and the integral resonant controller (i.e., $\mathrm{PD}+\mathrm{IRC}$ controllers) are suggested to control the nonlinear vibrations of the rotor AMB system for the first time. Therefore, the currents $i_x\left(t\right)$ and $i_y\left(t\right)$ are designed such that:

$$i_x(t)=k_{1}x(t)+k_2\dot{x}(t)+k_{3}u(t),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{i}_y(t)=k_{1}y(t)+k_2\dot{y}(t)+k_{4}v(t)$$

(6)

where $k_1$ is the linear proportional gain, $k_2$ denotes the linear derivative gain, $k_3$ represents the control gain of the integral resonant controller that is coupled to the horizontal vibration mode $x\left(t\right)$, and $k_4$ is the control gain of the integral resonant controller that is coupled to the vertical vibration mode $y\left(t\right)$.

Accordingly, the proposed integral resonant controllers are modeled as linear first-order filters as follows [28,29,30,31,32,33,34]:

$$\dot{u}(t)+{\rho }_{1}u(t)=k_{5}x(t)$$

(7)

$$\dot{v}(t)+{\rho }_{2}v(t)=k_{6}y(t)$$

(8)

where $u\left(t\right)$ and $\dot{u}\left(t\right)$ denote the displacement and velocity of the IRC that is connected to the horizontal vibration mode. $v\left(t\right)$ and $\dot{v}\left(t\right)$ denote the displacement and velocity of the IRC that connected to the vertical vibration mode. ${\rho }_1, {\rho }_2$ are constants, $k_5, k_6$ are the feedback control gains. Figure 2 shows the engineering implementation of the proposed control strategy in more detail.

Now, by substituting Equations (4)–(6) into Equation (3), we have

$$\begin{array}{l}{F}_1=\Delta {\left(\frac{I_0-k_{1}y-k_2\dot{y}-k_{4}v}{s_0+x\sin (\alpha )-y\cos (\alpha )}\right)}^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}{F}_2=\Delta {\left(\frac{I_0+k_{1}x+k_2\dot{x}+k_{3}u}{s_0+x\cos (\alpha )-y\sin (\alpha )}\right)}^2,\\ F_3=\Delta {\left(\frac{I_0+k_{1}x+k_2\dot{x}+k_{3}u}{s_0+x\cos (\alpha )+y\sin (\alpha )}\right)}^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}{F}_4=\Delta {\left(\frac{I_0+k_{1}y+k_2\dot{y}+k_{4}v}{s_0+x\sin (\alpha )+y\cos (\alpha )}\right)}^2,\\ F_5=\Delta {\left(\frac{I_0+k_{1}y+k_2\dot{y}+k_{4}v}{s_0-x\sin (\alpha )+y\cos (\alpha )}\right)}^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}{F}_6=\Delta {\left(\frac{I_0-k_{1}x-k_2\dot{x}-k_{3}u}{s_0-x\cos (\alpha )+y\sin (\alpha )}\right)}^2,\\ F_7=\Delta {\left(\frac{I_0-k_{1}x-k_2\dot{x}-k_{3}u}{s_0-x\cos (\alpha )-y\sin (\alpha )}\right)}^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}{F}_8=\Delta {\left(\frac{I_0-k_{1}y-k_2\dot{y}-k_{4}v}{s_0-x\sin (\alpha )-y\cos (\alpha )}\right)}^2\end{array}\}$$

(9)

Based on the system geometry shown in Figure 1B, the net control forces $F_{XC}$ and $F_{YC}$ in $X$ and $Y$ directions, can be expressed such that:

$$F_{XC}=(F_6+F_7-F_2-F_3)\cos (\alpha )+(F_5+F_8-F_1-F_4)\sin (\alpha )$$

(10)

$$F_{YC}=(F_1+F_8-F_4-F_5)\cos (\alpha )+(F_2+F_7-F_3-F_6)\sin (\alpha )$$

(11)

Expanding Equations (10) and (11) using Taylor series up to third-order approximation as given in Appendix A, and then substituting the resulting equations into Equations (1) and (2) with introducing the dimensionless variables and parameter $t^{*}={\omega }_{n}t, x^{*}=\frac{x}{s_0}, y^{*}=\frac{y}{s_0}, u^{*}=\frac{u}{s_0}, v^{*}=\frac{v}{s_0}, p=\frac{s_0}{I_0}{k}_1, d=\frac{s_0{\omega }_n}{I_0}{k}_2, {\eta }_1=-8\cos \left(\alpha \right)\frac{s_0}{I_0}{k}_3,$ ${\eta }_2=-8\cos \left(\alpha \right)\frac{s_0}{I_0}{k}_4, {\eta }_3=\frac{k_5}{{\omega }_n}, {\eta }_4=\frac{k_6}{{\omega }_n}, {\lambda }_1=\frac{{\rho }_1}{{\omega }_n}, {\lambda }_2=\frac{{\rho }_2}{{\omega }_n}, f=\frac{e}{s_0}, \mathsf{\Omega }=\frac{\upsilon }{{\omega }_n}, {\omega }_n=\sqrt{\mathsf{\Delta }/m{s}_0^3}$, we get the equations of motion (after neglecting the asterisks for brevity) as follows:

$$\begin{array}{l}\ddot{x}+2\mu \dot{x}+{\omega }^{2}x-({\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}\\ +{\beta }_1{x}^{2}u+{\beta }_{2}x\dot{x}u+{\beta }_{3}x{u}^2+{\beta }_{4}xyv+{\beta }_{5}x{v}^2+{\beta }_{6}x\dot{y}v+{\beta }_{7}u{y}^2)={\Omega }^{2}f\cos (\Omega t)+{\eta }_{1}u\end{array}$$

(12)

$$\begin{array}{l}\ddot{y}+2\mu \dot{y}+{\omega }^{2}y-({\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}\\ +{\gamma }_1{y}^{2}v+{\gamma }_{2}y\dot{y}v+{\gamma }_{3}y{v}^2+{\gamma }_{4}yxu+{\gamma }_{5}y{u}^2+{\gamma }_{6}y\dot{x}u+{\gamma }_{7}v{x}^2)={\Omega }^{2}f\sin (\Omega t)+{\eta }_{2}v\end{array}$$

(13)

$$\dot{u}+{\lambda }_{1}u={\eta }_{3}x$$

(14)

$$\dot{v}+{\lambda }_{2}v={\eta }_{4}y$$

(15)

where the coefficients $\mu , \omega , {\alpha }_j, {\beta }_j, {\gamma }_j, j=1,2,\dots ,7$ are given in Appendix B. Equations (12)–(15) represent the whole governing equations of motion of the controlled rotor system, where Equations (12) and (13) represent the equations of motion of the controlled rotor system, and Equations (14) and (15) are the equations of motion of the connected integral resonant controllers.

3. Analytical Investigations

To investigate the performance of the proposed control law (i.e., $\mathrm{IRC}$ controllers), the multiple scales perturbation method is applied within this section to obtain the periodic solutions of the nonlinear dynamical system given by Equations (12)–(15). Accordingly, the first-order approximate solution to Equations (12)–(15) is proposed as follows [37,38]:

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

(16)

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

(17)

$$u(t,\epsilon )=\epsilon u_0(T_0,T_1)+{\epsilon }^2{u}_1(T_0,T_1)$$

(18)

$$v(t,\epsilon )=\epsilon v_0(T_0,T_1)+{\epsilon }^2{v}_1(T_0,T_1)$$

(19)

where $\epsilon$ is an artificial parameter used as a book-keeping only [38], $T_0=t$, and $T_1=\epsilon t$. Accordingly, the derivatives $\frac{d}{dt}$ and $\frac{d^2}{d{t}^2}$ can be expressed in terms of $T_0$ and $T_1$ as follows:

$$\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},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=0,1$$

(20)

To obtain the system solution using the multiple scales method, the system parameters should be scaled such that:

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

(21)

Substituting Equations (16)–(21) into Equations (12) and (15) with equating the coefficients that have the same power of $\epsilon$, we get

O (${\epsilon }^0$):

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

(22)

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

(23)

O ($\epsilon$):

$$\left(D_0+{\lambda }_1\right)u_0={\widehat{\eta }}_3{x}_0$$

(24)

$$\left(D_0+{\lambda }_2\right)v_0={\widehat{\eta }}_4{y}_0$$

(25)

$$\begin{array}{ll}(D_0^2+{\omega }^2)x_1& =-2D_0{D}_1{x}_0-2\widehat{\mu }{D}_0{x}_0+{\widehat{\alpha }}_1{u}_0^3+{\widehat{\alpha }}_2{x}_0{y}_0^2+{\widehat{\alpha }}_3{x}_0^2{D}_0{x}_0+{\widehat{\alpha }}_4{y}_0^2{D}_0{x}_0\\ & +{\widehat{\alpha }}_5{x}_0{(D_0{y}_0)}^2+{\widehat{\alpha }}_6{x}_0{(D_0{x}_0)}^2+{\widehat{\alpha }}_7{x}_0{y}_0{D}_0{y}_0+{\beta }_1{x}_0^2{u}_0+{\beta }_2{x}_0{D}_0{x}_0{u}_0\\ & +{\beta }_3{x}_0{u}_0^2+{\beta }_4{x}_0{y}_0{v}_0+{\beta }_5{x}_0{v}_0^2+{\beta }_6{x}_0{D}_0{y}_0{v}_0+{\beta }_7{y}_0^2{u}_0+{\Omega }^2\widehat{f}\cos (\Omega t)\\ & +{\eta }_1{u}_0\end{array}$$

(26)

$$\begin{array}{ll}(D_0^2+{\omega }^2)y_1& =-2D_0{D}_1{y}_0-2\widehat{\mu }{D}_0{y}_0+{\widehat{\alpha }}_1{y}_0^3+{\widehat{\alpha }}_2{y}_0{x}_0^2+{\widehat{\alpha }}_3{y}_0^2{D}_0{y}_0+{\widehat{\alpha }}_4{x}_0^2{D}_0{y}_0\\ & +{\widehat{\alpha }}_5{y}_0{(D_0{x}_0)}^2+{\widehat{\alpha }}_6{y}_0{(D_0{y}_0)}^2+{\widehat{\alpha }}_7{y}_0{x}_0{D}_0{x}_0+{\gamma }_1{y}_0^2{v}_0+{\gamma }_2{y}_0{D}_0{y}_0{v}_0\\ & +{\gamma }_3{y}_0{v}_0^2+{\gamma }_4{y}_0{x}_0{u}_0+{\gamma }_5{y}_0{u}_0^2+{\gamma }_6{y}_0{D}_0{x}_0{u}_0+{\gamma }_7{x}_0^2{v}_0+{\Omega }^2\widehat{f}\sin (\Omega t)\\ & +{\eta }_2{v}_0\end{array}$$

(27)

The steady-state solution of Equations (22)–(25) can be expressed as follows:

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

(28)

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

(29)

$$u_0(T_0,T_1)={\delta }_{1}A(T_1)e^{i\omega T_0}+{\overline{\delta }}_1\overline{A}(T_1)e^{-i\omega T_0}$$

(30)

$$v_0(T_0,T_1)={\delta }_{2}B(T_1)e^{i\omega T_0}+{\overline{\delta }}_2\overline{B}(T_1)e^{-i\omega T_0}$$

(31)

where $i=\sqrt{-1}, {\delta }_1=\frac{{\lambda }_1-i\omega }{{\lambda }_1^2+{\omega }^2}{\widehat{\eta }}_3, {\delta }_2=\frac{{\lambda }_2-i\omega }{{\lambda }_2^2+{\omega }^2}{\widehat{\eta }}_4$, and $A\left(T_1\right)$, $B\left(T_1\right)$ are unknown functions that will be determined later. By substituting Equations (28)–(31) into Equations (26) and (27), we have

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

(32)

$$\begin{array}{ll}(D_0^2+{\omega }^2)y_1& =[-2i\omega (D_{1}B)-2i\widehat{\mu }\omega B+3{\widehat{\alpha }}_1{B}^2\overline{B}+2{\widehat{\alpha }}_{2}BA\overline{A}+{\widehat{\alpha }}_2\overline{B}{A}^2+i{\widehat{\alpha }}_3\omega B^2\overline{B}\\ & +2i{\widehat{\alpha }}_4\omega BA\overline{A}-i{\widehat{\alpha }}_4\omega \overline{B}{A}^2+2{\widehat{\alpha }}_5{\omega }^{2}BA\overline{A}-{\widehat{\alpha }}_5{\omega }^2\overline{B}{A}^2\hspace{0.17em}+{\widehat{\alpha }}_6{\omega }^2{B}^2\overline{B}\\ & +i{\widehat{\alpha }}_7\omega \overline{B}{A}^2+2{\gamma }_1{\delta }_1{B}^2\overline{B}+{\gamma }_3{\delta }_2^2{B}^2\overline{B}+{\gamma }_4({\delta }_{1}BA\overline{A}+{\delta }_1\overline{B}{A}^2)+{\gamma }_5{\delta }_1^2\overline{B}{A}^2\\ & +{\gamma }_6(-i\omega {\delta }_{1}BA\overline{A}+i\omega {\delta }_1\overline{B}{A}^2)+2{\gamma }_7{\delta }_{1}BA\overline{A}+{\eta }_2{\delta }_{1}B]e^{i\omega T_0}+[{\widehat{\alpha }}_1{B}^3\\ & +{\widehat{\alpha }}_{2}B{A}^2+i{\widehat{\alpha }}_3\omega B^3+i{\widehat{\alpha }}_4\omega B{A}^2-{\widehat{\alpha }}_5{\omega }^{2}B{A}^2-{\widehat{\alpha }}_6{\omega }^2{B}^3+i{\widehat{\alpha }}_7\omega B{A}^2\\ & +{\gamma }_1{\delta }_1{B}^3+i{\gamma }_2\omega N{B}^3{e}^{3i\omega T_0}+{\gamma }_3{\delta }_2^2{B}^3+{\gamma }_4{\delta }_{1}B{A}^2{e}^{3i\omega T_0}+{\gamma }_5{\delta }_1^{2}B{A}^2\\ & +i{\gamma }_6\omega {\delta }_{1}B{A}^2+{\gamma }_7{\delta }_{1}B{A}^2]e^{3i\omega T_0}-\frac{1}{2}i{\Omega }^2\widehat{f}{e}^{i\Omega T_0}+cc\end{array}$$

(33)

where $cc$ denotes the complex conjugate term. To obtain the solvability conditions of Equations (32) and (33), the closeness of the rotor spinning speed ($\mathsf{\Omega }$) to the system natural frequency ($\omega$) is described via introducing the detuning parameter $\sigma$ as follows:

$$\Omega =\omega +\sigma =\omega +\epsilon \widehat{\sigma }$$

(34)

Substituting Equation (34) into Equations (32) and (33), one can obtain the following solvability conditions:

$$\begin{array}{l}-2i\omega (D_{1}A)-2i\widehat{\mu }\omega A+3\hspace{0.17em}{\widehat{\alpha }}_1{A}^2\overline{A}+2\hspace{0.17em}{\widehat{\alpha }}_{2}AB\overline{B}+{\widehat{\alpha }}_2\overline{A}{B}^2+i\omega \hspace{0.17em}{\widehat{\alpha }}_3{A}^2\overline{A}+2i\omega \hspace{0.17em}{\widehat{\alpha }}_{4}AB\overline{B}\\ -i\omega \hspace{0.17em}{\widehat{\alpha }}_4\overline{A}{B}^2+2{\omega }^2\hspace{0.17em}{\widehat{\alpha }}_{5}AB\overline{B}-{\omega }^2\hspace{0.17em}{\widehat{\alpha }}_5\overline{A}{B}^2+{\omega }^2\hspace{0.17em}{\widehat{\alpha }}_6{A}^2\overline{A}+i\omega \hspace{0.17em}{\widehat{\alpha }}_7\overline{A}{B}^2+2{\beta }_1{\delta }_1{A}^2\overline{A}\\ +{\beta }_3{\delta }_1^2{A}^2\overline{A}+{\beta }_4{\delta }_{2}AB\overline{B}+{\beta }_4{\delta }_2\overline{A}{B}^2+{\beta }_5{\delta }_2^2{N}^2\overline{A}{B}^2-i\omega {\beta }_6{\delta }_{2}AB\overline{B}+i\omega {\beta }_6{\delta }_2\overline{A}{B}^2\\ +2{\beta }_7{\delta }_{1}AB\overline{B}+{\eta }_1{\delta }_{1}A+\frac{1}{2}{(\omega +\sigma )}^2\widehat{f}{e}^{i\epsilon \widehat{\sigma }{T}_0}=0\end{array}$$

(35)

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

(36)

To obtain the amplitude-phase equations of the controlled system, the functions $A\left(T_1\right)$ and $B\left(T_1\right)$ can be expressed in the polar form as follows:

$$A(T_1)=\frac{1}{2}a(T_1)e^{i{\theta }_1(T_1)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}B(T_1)=\frac{1}{2}b(T_1)e^{i{\theta }_2(T_1)}\hspace{0.17em}\hspace{0.17em}$$

(37)

Substituting Equation (37) into Equations (35) and (36), and then separating the real and imaginary parts, we have the following amplitude-phase modulating equations:

$$\begin{array}{ll}\dot{a}& =G_1(a,b,{\phi }_1,{\phi }_2)=-\left(\mu +\frac{{\eta }_3{\eta }_1}{2({\lambda }_1^2+{\omega }^2)}\right)a+\frac{1}{8}\left({\alpha }_3-\frac{2{\beta }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2}-\frac{2{\beta }_3{\lambda }_1{\eta }_3^2}{{({\lambda }_1^2+{\omega }^2)}^2}\right)a^3\\ & +\frac{1}{8}(2{\alpha }_4-\frac{{\beta }_4{\eta }_4}{{\lambda }_2^2+{\omega }^2}-\frac{{\beta }_6{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2}-\frac{2{\beta }_7{\eta }_3}{{\lambda }_1^2+{\omega }^2})a{b}^2+\frac{1}{8}(-{\alpha }_4+{\alpha }_7-\frac{{\beta }_4{\eta }_4}{{\lambda }_2^2+{\omega }^2}\\ & +\frac{2{\beta }_5{\lambda }_2{\eta }_4^2}{{({\lambda }_2^2+{\omega }^2)}^2}+\frac{{\beta }_6{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2})a{b}^2\cos (2{\phi }_2-2{\phi }_1)-\frac{1}{8}(\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega +\frac{{\beta }_4{\eta }_4{\lambda }_2}{\omega ({\lambda }_2^2+{\omega }^2)}\\ & +\frac{{\beta }_5{\eta }_4^2({\lambda }_2^2-{\omega }^2)}{\omega {({\lambda }_2^2+{\omega }^2)}^2}-\frac{{\beta }_6{\eta }_4\omega }{{\lambda }_2^2+{\omega }^2})a{b}^2\sin (2{\phi }_2-2{\phi }_1)+\frac{f}{2\omega }{\left(\omega +\sigma \right)}^2\sin ({\phi }_1)\end{array}$$

(38)

$$\begin{array}{ll}\dot{b}& =G_2(a,b,{\phi }_1,{\phi }_2)=-\left(\mu +\frac{{\eta }_2{\eta }_4}{2({\lambda }_2^2+{\omega }^2)}\right)b+\frac{1}{8}\left({\alpha }_3-\frac{2{\gamma }_1{\eta }_4}{{\lambda }_2^2+{\omega }^2}-\frac{2{\gamma }_3{\lambda }_2{\eta }_4^2}{{({\lambda }_2^2+{\omega }^2)}^2}\right)b^3\\ & +\frac{1}{8}(2{\alpha }_4-\frac{{\gamma }_4{\eta }_3}{{\lambda }_1^2+{\omega }^2}-\frac{{\gamma }_6{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2}-\frac{2{\gamma }_7{\eta }_4}{{\lambda }_2^2+{\omega }^2})a^{2}b+\frac{1}{8}(-{\alpha }_4+{\alpha }_7-\frac{{\gamma }_4{\eta }_3}{{\lambda }_1^2+{\omega }^2}\\ & -\frac{2{\gamma }_5{\lambda }_1{\eta }_3^2}{{({\lambda }_1^2+{\omega }^2)}^2}+\frac{{\gamma }_6{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2})a^{2}b\cos (2{\phi }_2-2{\phi }_1)+\frac{1}{8}(\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega +\frac{{\gamma }_4{\lambda }_1{\eta }_3}{\omega ({\lambda }_1^2+{\omega }^2)}\\ & +\frac{{\gamma }_5{\eta }_3^2({\lambda }_1^2-{\omega }^2)}{\omega {({\lambda }_1^2+{\omega }^2)}^2}+\frac{{\gamma }_6{\eta }_3\omega }{{\lambda }_1^2+{\omega }^2})a^{2}b\sin (2{\phi }_2-2{\phi }_1)-\frac{f}{2\omega }{\left(\omega +\sigma \right)}^2\cos ({\phi }_2)\end{array}$$

(39)

$$\begin{array}{ll}{\dot{\phi }}_1& =G_3(a,b,{\phi }_1,{\phi }_2)=\left(\sigma +\frac{{\lambda }_1{\eta }_3{\eta }_1}{2\omega ({\lambda }_1^2+{\omega }^2)}\right)+\frac{1}{8\omega }(3{\alpha }_1+{\alpha }_6{\omega }^2+\frac{2{\beta }_1{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2}\\ & +\frac{{\beta }_3{\eta }_3^2({\lambda }_1^2-{\omega }^2)}{{({\lambda }_1^2+{\omega }^2)}^2})a^2+\frac{1}{8\omega }\left(2{\alpha }_2+2{\alpha }_5{\omega }^2+\frac{{\beta }_4{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2}-\frac{{\beta }_6{\eta }_4{\omega }^2}{{\lambda }_2^2+{\omega }^2}+\frac{2{\beta }_7{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2}\right)b^2\\ & +\frac{1}{8\omega }({\alpha }_2-{\alpha }_5{\omega }^2+\frac{{\beta }_4{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2}+\frac{{\beta }_5{\eta }_4^2({\lambda }_2^2-{\omega }^2)}{{({\lambda }_2^2+{\omega }^2)}^2}+\frac{{\beta }_6{\eta }_4{\omega }^2}{{\lambda }_2^2+{\omega }^2})b^2\cos (2{\phi }_2-2{\phi }_1)\\ & +\frac{1}{8}(-{\alpha }_4+{\alpha }_7-\frac{{\beta }_4{\eta }_4}{{\lambda }_2^2+{\omega }^2}-\frac{2{\beta }_5{\lambda }_2{\eta }_4^2}{{({\lambda }_2^2+{\omega }^2)}^2}+\frac{{\beta }_6{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2})b^2\sin (2{\phi }_2-2{\phi }_1)\\ & +\frac{f}{2a\omega }{\left(\omega +\sigma \right)}^2\cos ({\phi }_1)\end{array}$$

(40)

$$\begin{array}{ll}{\dot{\phi }}_2& =G_4(a,b,{\phi }_1,{\phi }_2)=\left(\sigma +\frac{{\lambda }_2{\eta }_2{\eta }_4}{2\omega ({\lambda }_2^2+{\omega }^2)}\right)+\frac{1}{8\omega }(3{\alpha }_1+{\alpha }_6{\omega }^2+\frac{2{\gamma }_1{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2}\\ & +\frac{{\gamma }_3{\eta }_4^2({\lambda }_2^2-{\omega }^2)}{{({\lambda }_2^2+{\omega }^2)}^2})b^2+\frac{1}{8\omega }\left(2{\alpha }_2+2{\alpha }_5{\omega }^2+\frac{{\gamma }_4{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2}-\frac{{\gamma }_6{\eta }_3{\omega }^2}{{\lambda }_1^2+{\omega }^2}+\frac{2{\gamma }_7{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2}\right)a^2\\ & +\frac{1}{8\omega }({\alpha }_2-{\alpha }_5{\omega }^2+\frac{{\gamma }_4{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2}+\frac{{\gamma }_5{\eta }_3^2({\lambda }_1^2-{\omega }^2)}{{({\lambda }_1^2+{\omega }^2)}^2}+\frac{{\gamma }_6{\eta }_3{\omega }^2}{{\lambda }_1^2+{\omega }^2})a^2\cos (2{\phi }_2-2{\phi }_1)\\ & -\frac{1}{8}(-{\alpha }_4+{\alpha }_7+\frac{{\gamma }_4{\eta }_3}{{\lambda }_1^2+{\omega }^2}-\frac{2{\gamma }_5{\lambda }_1{\eta }_3^2}{{({\lambda }_1^2+{\omega }^2)}^2}-\frac{{\gamma }_6{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2})a^2\sin (2{\phi }_2-2{\phi }_1)\\ & +\frac{f}{2b\omega }{\left(\omega +\sigma \right)}^2\sin ({\phi }_2)\end{array}$$

(41)

Inserting Equations (28)–(31) and (37) into Equations (16)–(19) with considering Equation (37), one can obtain the analytical solution of the coupled system given by Equations (12)–(15) as follows:

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

(42)

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

(43)

$$u(t)=\frac{{\eta }_3}{({\lambda }_1^2+{\omega }^2)}a(t)[{\lambda }_1\cos (\Omega t-{\phi }_1(t))+\omega \sin (\Omega t-{\phi }_1(t))]$$

(44)

$$v(t)=\frac{{\eta }_4}{({\lambda }_2^2+{\omega }^2)}b(t)[{\lambda }_2\cos (\Omega t-{\phi }_2(t))+\omega \sin (\Omega t-{\phi }_2(t))]$$

(45)

where $a\left(t\right)$ and $b\left(t\right)$ are the system oscillation amplitudes in $X$ and $Y$ directions, respectively, and ${\phi }_1\left(t\right)=\sigma t-{\theta }_1, {\phi }_2\left(t\right)=\sigma t-{\theta }_2$ are the corresponding phase angles. At steady-state oscillations, we have $\dot{a}=\dot{b}={\dot{\phi }}_1={\dot{\phi }}_2=0$. Substituting $\dot{a}=\dot{b}={\dot{\phi }}_1={\dot{\phi }}_2=0$ into Equations (38)–(41), we get

$$\begin{array}{ll}{f}_1(a,b,{\phi }_1,{\phi }_2)& =-\left(\mu +\frac{{\eta }_3{\eta }_1}{2({\lambda }_1^2+{\omega }^2)}\right)a+\frac{1}{8}\left({\alpha }_3-\frac{2{\beta }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2}-\frac{2{\beta }_3{\lambda }_1{\eta }_3^2}{{({\lambda }_1^2+{\omega }^2)}^2}\right)a^3\\ & +\frac{1}{8}(2{\alpha }_4-\frac{{\beta }_4{\eta }_4}{{\lambda }_2^2+{\omega }^2}-\frac{{\beta }_6{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2}-\frac{2{\beta }_7{\eta }_3}{{\lambda }_1^2+{\omega }^2})a{b}^2+\frac{1}{8}(-{\alpha }_4+{\alpha }_7\\ & -\frac{{\beta }_4{\eta }_4}{{\lambda }_2^2+{\omega }^2}+\frac{2{\beta }_5{\lambda }_2{\eta }_4^2}{{({\lambda }_2^2+{\omega }^2)}^2}+\frac{{\beta }_6{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2})a{b}^2\cos (2{\phi }_2-2{\phi }_1)\\ & -\frac{1}{8}(\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega +\frac{{\beta }_4{\eta }_4{\lambda }_2}{\omega ({\lambda }_2^2+{\omega }^2)}+\frac{{\beta }_5{\eta }_4^2({\lambda }_2^2-{\omega }^2)}{\omega {({\lambda }_2^2+{\omega }^2)}^2}\\ & -\frac{{\beta }_6{\eta }_4\omega }{{\lambda }_2^2+{\omega }^2})a{b}^2\sin (2{\phi }_2-2{\phi }_1)+\frac{f}{2\omega }{\left(\omega +\sigma \right)}^2\sin ({\phi }_1)=0\end{array}$$

(46)

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

(47)

$$\begin{array}{ll}{f}_3(a,b,{\phi }_1,{\phi }_2)& =\left(\sigma +\frac{{\lambda }_1{\eta }_3{\eta }_1}{2\omega ({\lambda }_1^2+{\omega }^2)}\right)+\frac{1}{8\omega }(3{\alpha }_1+{\alpha }_6{\omega }^2+\frac{2{\beta }_1{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2}\\ & +\frac{{\beta }_3{\eta }_3^2({\lambda }_1^2-{\omega }^2)}{{({\lambda }_1^2+{\omega }^2)}^2})a^2+\frac{1}{8\omega }(2{\alpha }_2+2{\alpha }_5{\omega }^2+\frac{{\beta }_4{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2}-\frac{{\beta }_6{\eta }_4{\omega }^2}{{\lambda }_2^2+{\omega }^2}\\ & +\frac{2{\beta }_7{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2})b^2+\frac{1}{8\omega }({\alpha }_2-{\alpha }_5{\omega }^2+\frac{{\beta }_4{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2}+\frac{{\beta }_5{\eta }_4^2({\lambda }_2^2-{\omega }^2)}{{({\lambda }_2^2+{\omega }^2)}^2}\\ & +\frac{{\beta }_6{\eta }_4{\omega }^2}{{\lambda }_2^2+{\omega }^2})b^2\cos (2{\phi }_2-2{\phi }_1)+\frac{1}{8}(-{\alpha }_4+{\alpha }_7-\frac{{\beta }_4{\eta }_4}{{\lambda }_2^2+{\omega }^2}\\ & -\frac{2{\beta }_5{\lambda }_2{\eta }_4^2}{{({\lambda }_2^2+{\omega }^2)}^2}+\frac{{\beta }_6{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2})b^2\sin (2{\phi }_2-2{\phi }_1)+\frac{f}{2a\omega }{\left(\omega +\sigma \right)}^2\cos ({\phi }_1)=0\end{array}$$

(48)

$$\begin{array}{ll}{f}_4(a,b,{\phi }_1,{\phi }_2)& =\left(\sigma +\frac{{\lambda }_2{\eta }_2{\eta }_4}{2\omega ({\lambda }_2^2+{\omega }^2)}\right)+\frac{1}{8\omega }\left(3{\alpha }_1+{\alpha }_6{\omega }^2+\frac{2{\gamma }_1{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2}+\frac{{\gamma }_3{\eta }_4^2({\lambda }_2^2-{\omega }^2)}{{({\lambda }_2^2+{\omega }^2)}^2}\right)b^2\\ & +\frac{1}{8\omega }(2{\alpha }_2+2{\alpha }_5{\omega }^2+\frac{{\gamma }_4{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2}-\frac{{\gamma }_6{\eta }_3{\omega }^2}{{\lambda }_1^2+{\omega }^2}+\frac{2{\gamma }_7{\lambda }_2{\eta }_4}{{\lambda }_2^2+{\omega }^2})a^2\\ & +\frac{1}{8\omega }({\alpha }_2-{\alpha }_5{\omega }^2+\frac{{\gamma }_4{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2}+\frac{{\gamma }_5{\eta }_3^2({\lambda }_1^2-{\omega }^2)}{{({\lambda }_1^2+{\omega }^2)}^2}+\frac{{\gamma }_6{\eta }_3{\omega }^2}{{\lambda }_1^2+{\omega }^2})a^2\cos (2{\phi }_2-2{\phi }_1)\\ & -\frac{1}{8}(-{\alpha }_4+{\alpha }_7+\frac{{\gamma }_4{\eta }_3}{{\lambda }_1^2+{\omega }^2}-\frac{2{\gamma }_5{\lambda }_1{\eta }_3^2}{{({\lambda }_1^2+{\omega }^2)}^2}-\frac{{\gamma }_6{\lambda }_1{\eta }_3}{{\lambda }_1^2+{\omega }^2})a^2\sin (2{\phi }_2-2{\phi }_1)\\ & +\frac{f}{2b\omega }{\left(\omega +\sigma \right)}^2\sin ({\phi }_2)=0\end{array}$$

(49)

The nonlinear algebraic system $f_j\left(a,b,{\phi }_1,{\phi }_2\right)=0, j=1,2,3,4$ given by Equations (46)–(49) governs the steady-state oscillation amplitudes ($a$ and $b$) and phase angles (${\phi }_1$ and ${\phi }_2$) of the controlled rotor system as a function of IRC control parameters (${\eta }_1, {\eta }_2, {\eta }_3, {\eta }_4, {\lambda }_1, {\lambda }_2$). Therefore, one can investigate the performance of the applied IRC in suppressing the steady-state oscillation amplitudes of the controlled system via solving the nonlinear system $f_j\left(a,b,{\phi }_1,{\phi }_2\right)=0, j=1,2,3,4$ utilizing $f$ or $\sigma$ as a bifurcation parameter as illustrated in Section 4. In addition, the solution stability of Equations (46)–(49) can be explored via linearizing the nonlinear system given by Equations (38)–(41). Therefore, if we assume ($a_0, b_0,{\phi }_{10},{\phi }_{20}$) is the solution of $f_j\left(a,b,{\phi }_1,{\phi }_2\right)=0, j=1,2,3,4$, and ($a_1, b_1,{\phi }_{11},{\phi }_{21}$) is small deviations about that solution [39]. Accordingly, we have

$$\begin{array}{l}a=a_0+a_1,\hspace{0.17em}\hspace{0.17em}b=b_0+b_1,\hspace{0.17em}\hspace{0.17em}{\phi }_1={\phi }_{10}+{\phi }_{11},\hspace{0.17em}\hspace{0.17em}{\phi }_2={\phi }_{20}+{\phi }_{21},\\ \dot{a}={\dot{a}}_1,\hspace{0.17em}\hspace{0.17em}\dot{b}={\dot{b}}_1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\dot{\phi }}_1={\dot{\phi }}_{11},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\dot{\phi }}_2={\dot{\phi }}_{21}\end{array}\}$$

(50)

Substituting Equation (50) into Equations (38)–(41), one can obtain the following linear autonomous system:

$$\left(\begin{array}{c}{\dot{a}}_1\\ {\dot{b}}_1\\ {\dot{\phi }}_{11}\\ {\dot{\phi }}_{21}\end{array}\right)=\left(\begin{array}{cccc}\frac{\partial G_1}{\partial a_1}& \frac{\partial G_1}{\partial b_1}& \frac{\partial G_1}{\partial {\phi }_{11}}& \frac{\partial G_1}{\partial {\phi }_{21}}\\ \frac{\partial G_2}{\partial a_1}& \frac{\partial G_2}{\partial b_1}& \frac{\partial G_2}{\partial {\phi }_{11}}& \frac{\partial G_2}{\partial {\phi }_{21}}\\ \frac{\partial G_3}{\partial a_1}& \frac{\partial G_3}{\partial b_1}& \frac{\partial G_3}{\partial {\phi }_{11}}& \frac{\partial G_3}{\partial {\phi }_{21}}\\ \frac{\partial G_4}{\partial a_1}& \frac{\partial G_4}{\partial b_1}& \frac{\partial G_4}{\partial {\phi }_{11}}& \frac{\partial G_4}{\partial {\phi }_{21}}\end{array}\right)\left(\begin{array}{c}{a}_1\\ b_1\\ {\phi }_{11}\\ {\phi }_{21}\end{array}\right)$$

(51)

where the coefficient of the above square Jacobian matrix is given in Appendix C. The linear dynamical system given by Equation (51) is topologically equivalent to the nonlinear system given by Equations (38)–(41). Accordingly, one can check the solution stability of Equations (46)–(49) via exploring the eigenvalues of Equation (51). In the following section, the performance of the IRC is explored via solving Equations (46)–(49) at different values of the control parameters (${\eta }_1, {\eta }_2, {\eta }_3, {\eta }_4, {\lambda }_1, {\lambda }_2$) utilizing $f$ or $\sigma$ as the bifurcation parameters. In addition, the stability of the obtained solution (i.e., solution of Equations (46)–(49)) is investigated via checking the eigenvalues of Equation (51).

4. Bifurcation Analysis and Control Performance

The sensitivity analysis and control performance of the applied integral resonant controller in suppressing the nonlinear oscillations of the rotor system are investigated within this section. The different bifurcation diagrams are obtained via solving the nonlinear system $f_j\left(a,b,{\phi }_1,{\phi }_2\right)=0, j=1,2,3,4$ utilizing one of the system or controller parameters ($\sigma , f, {\eta }_1, {\eta }_2,{\eta }_3, {\eta }_4$) [40]. In addition, the solution stability of Equations (46)–(49) has been studied via checking the eigenvalues of the linearized model given by Equation (51). During plotting the different bifurcation diagrams via solving $f_j\left(a,b,{\phi }_1,{\phi }_2\right)=0, j=1,2,3,4$, the proposed MATLAB code is designed to plot the stable solution as a solid line and the unstable solution as a dotted line. Moreover, to guarantee the accuracy of the obtained bifurcation diagrams, numerical validations are performed via solving numerically Equations (12)–(15) using ODE45 (MATLAB function). The numerical solution is plotted as a circle when sweeping the bifurcation parameter forward and as a big dot when sweeping the bifurcation parameter backward. The system and controllers parameters are selected such that: $f=0.015, p=1.22, d=0.005, \alpha ={22.5}^{°}$,${\lambda }_1={\lambda }_2=1.0,{\eta }_1={\eta }_2={\eta }_3={\eta }_4=0.0$, and $\mathsf{\Omega }=\omega +\sigma$ [7,8]. Before proceeding further, it is important to remember that the dimensionless control parameters ${\eta }_1, {\eta }_2, {\eta }_3, {\eta }_4, {\lambda }_1,$ and ${\lambda }_2$ are defined before Equation (12) such that ${\eta }_1=-\frac{8s_0\cos \left(\alpha \right)}{I_0}{k}_3, {\eta }_2=-\frac{8\cos \left(\alpha \right)s_0}{I_0}{k}_4,{\eta }_3=\frac{k_5}{{\omega }_n}, {\eta }_4=\frac{k_6}{{\omega }_n}, {\lambda }_1=\frac{{\rho }_1}{{\omega }_n},$ and ${\lambda }_2=\frac{{\rho }_2}{{\omega }_n}$. Therefore, one can easily deduce that ${\eta }_1$ and ${\eta }_2$ represent the dimensionless control signal gains of the IRCs, ${\eta }_3$ and ${\eta }_4$ denote the dimensionless feedback signal gains of the IRCs, and ${\lambda }_1$ and ${\lambda }_2$ are the internal feedback gains of the IRCs (i.e., see Figure 2). Accordingly, the influence of the IRCs parameters (${\eta }_1, {\eta }_2, {\eta }_3, {\eta }_4, {\lambda }_1,$ and ${\lambda }_2$) on the dynamical behaviors of the rotor system is investigated within the following subsections considering the proportional and derivative control gain are fixed constant (i.e., $p=1.22, d=0.005$).

4.1. The Rotor System Dynamics without IRC

The steady-state nonlinear vibration of the rotor system is explored within this section when turning the integral resonant controller off (i.e., when ${\eta }_1={\eta }_2={\eta }_3={\eta }_4={\lambda }_1={\lambda }_2=0.0$), while the proportional gain $p=1.22$ and the derivative gain $d=0.005$. Remember from Equation (34) that $\sigma$ represents the closeness of the rotor angular speed $\mathsf{\Omega }$ to its natural frequency $\omega$. Accordingly, $\sigma$ is utilized within this article to describe the system oscillatory behaviors close to the resonance conditions (i.e., when $\mathsf{\Omega }\to \mathsf{\omega }$). In Figure 3, the rotor system steady-state vibration amplitudes ($a$ and $b$) are plotted versus $\sigma$ at two different values of the rotor eccentricity $f=0.01$ and $0.02$. Figure 3A,B show that the rotor system exhibits a bistable periodic attractor if $-0.018<\sigma <0.033$, otherwise, the rotor responds as a linear system with a single periodic attractor when $f=0.01$. In addition, Figure 3C,D illustrate that the rotor system has a complex $\sigma -$response curve when increasing the eccentricity to $f=0.02$. Figure 3C,D illustrates that the rotor system can oscillate with one of three oscillations modes, which are: single periodic solutions, bistable periodic solutions, and tri-stable periodic solutions.

To visualize the bifurcation behaviors of the rotor system for a wide range of the eccentricity, $f$ is utilized as the main bifurcation parameter as shown in Figure 4 at two different values of the rotor spinning speed $\mathsf{\Omega }=\omega +\sigma , \sigma =0.0, 0.05$. It is clear from the figure that the rotor system may lose its stability at strong eccentricity beside the three oscillations modes reported in Figure 3. Due to the complex dynamical behaviors of the considered system that are reported in Figure 3 and Figure 4, the integral resonant controller is introduced as a novel control technique to mitigate the nonlinear oscillation and to suppress the motion bifurcations of such system as discussed in Section 4.2.

4.2. The Rotor System Dynamics with IRC

This section is dedicated to exploring the influence of one of the connected integral resonant controllers on the oscillatory behaviors of the rotor system when the other one is turned off. Before investigating the effects of the IRC using the different bifurcation diagrams, let us first examine the derived amplitude-phase equations (i.e., Equations (38)–(41)). It is clear from Equations (38)–(41) that the coupling of the IRC to the rotor system has modified the rotor linear damping coefficient ($\mu$) to the equivalent damping coefficients ${\mu }_x=\mu +\frac{{\eta }_1{\eta }_3}{2\left({\lambda }_1^2+{\omega }^2\right)}$ for the horizontal oscillation mode, and ${\mu }_y=\mu +\frac{{\eta }_2{\eta }_4}{2\left({\lambda }_2^2+{\omega }^2\right)}$ for the vertical oscillation mode. In addition, the detuning parameter ($\sigma$) has been modified to ${\sigma }_x=\sigma +\frac{{\lambda }_1{\eta }_1{\eta }_3}{2\omega \left({\lambda }_1^2+{\omega }^2\right)}$ for the horizontal oscillation mode, and ${\sigma }_y=\sigma +\frac{{\lambda }_2{\eta }_2{\eta }_4}{2\omega \left({\lambda }_2^2+{\omega }^2\right)}$ for the vertical oscillation mode. Accordingly, the IRC adds the linear damping magnitude $\frac{{\eta }_1{\eta }_3}{2\left({\lambda }_1^2+{\omega }^2\right)}$ to the horizontal oscillation mode and the linear damping magnitude $\frac{{\eta }_2{\eta }_4}{2\left({\lambda }_2^2+{\omega }^2\right)}$ to the vertical oscillation mode. It is worthy to note that ${\mu }_x$ is proportional to the product of both the control signal gain (${\eta }_1$) and feedback signal gain (${\eta }_3$), and inversely proportional to the square of the internal loop feedback gain (${\lambda }_1$) of the IRC that is connected to the horizontal oscillation mode. In addition, ${\mu }_y$ is proportional to the product of both the control signal gain (${\eta }_2$) and feedback signal gain (${\eta }_4$), and inversely proportional to the square of the internal loop feedback gain (${\lambda }_2$) of the IRC that is connected to the vertical oscillation mode. Moreover, it should be noted that ${\eta }_1, {\eta }_2, {\eta }_3, {\eta }_4, {\lambda }_1,$ and ${\lambda }_2$ are the dimensionless forms of the original parameters $k_3, k_4, k_5, k_6, {\rho }_1,$ and ${\rho }_2$, respectively. Therefore, we can deduce that the IRC that connected to the horizontal oscillation mode adds linear damping magnitude proportional to the product of the actual feedback and control signal gains ($k_3{k}_5$), while the IRC that is connected to the vertical oscillation mode adds linear damping magnitude proportional to the product of the actual feedback and control signal gains ($k_4{k}_6$) as shown in Figure 2.

For simplicity, let us denote the IRC that is connected to the horizontal oscillation mode as U-IRC and the IRC that is connected to the vertical oscillation mode as V-IRC. Based on the previous discussion, setting ${\eta }_1$ or ${\eta }_3$ or both to be equal to zero means the U-IRC is turned off. Also, setting ${\eta }_2$ or ${\eta }_4$ or both to be equal to zero means the V-IRC is turned off. Accordingly, the influence of the proposed IRC on the response curves of the rotor system has been investigated in two stages of analysis. Firstly, it is considered that the U-IRC is in action, while the V-IRC is deactivated via setting ${\eta }_2=0$ or ${\eta }_4=0$ or ${\eta }_2={\eta }_4=0$. Secondly, the effect of both U-IRC and V-IRC on the rotor system is investigated simultaneously as one symmetric IRC.

• The influence of asymmetric IRC on the rotor system dynamics

The effect of the U-IRC on the rotor system dynamics has been explored within this section considering the V-IRC controller is turned off via setting ${\eta }_2=0$ or ${\eta }_4=0$ or ${\eta }_2={\eta }_4=0$. Figure 5 shows the controlled rotor system $\sigma -$response curve at three different values of the control signal gain ${\eta }_1$ (i.e., ${\eta }_1=0.05, 0.2$, and $0.4$) when ${\eta }_2=0, {\eta }_3={\eta }_4=0.2$. It is clear from the figure that the increase of the control signal gain ${\eta }_1$, decreases the system oscillations in the horizontal direction only, and eliminates the complex bifurcation behaviors of the whole system close to the resonance condition (i.e., when $\sigma$ close to or equal to zero). The obtained results in Figure 5 can be explained based on the modified linear damping coefficients ${\mu }_x=\mu +\frac{{\eta }_1{\eta }_3}{2\left({\lambda }_1^2+{\omega }^2\right)}$ and ${\mu }_y=\mu +\frac{{\eta }_2{\eta }_4}{2\left({\lambda }_2^2+{\omega }^2\right)}$, where increasing ${\eta }_1$ with letting ${\eta }_2=0, {\eta }_3={\eta }_4=0.2$ means increasing ${\mu }_x$ with fixing ${\mu }_y=\mu$. Therefore, the system oscillation at the horizontal direction decreases monotonically with increasing ${\eta }_1$.

Numerical simulation for the controlled system according to Figure 5 when $\sigma =0.0$ is illustrated in Figure 6 and Figure 7 via solving Equations (12)–(15) using ODE45 at the two different initial conditions $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=u\left(0\right)=v\left(0\right)=0$ and $x\left(0\right)=0.2, \dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=u\left(0\right)=v\left(0\right)=0$. Figure 6 shows the steady-state time-response and the corresponding orbital motions of the rotor system and the connected IRCs when $\sigma =0.0, {\eta }_1=0.05, {\eta }_2=0, {\eta }_3={\eta }_4=0.2$. It is clear from the figure that the rotor system is sensitive to the initial conditions, where the system has two bistable periodic solutions as reported in Figure 5A,B. Moreover, the rotor whirling motion may be forward or backward depending on the initial conditions.

Figure 7 illustrates the time-response and corresponding periodic orbits of the rotor system and IRCs according to Figure 5E,F when $\sigma =0.0, {\eta }_1=0.4, {\eta }_2=0, {\eta }_3={\eta }_4=0.2$ using the same initial conditions that were used to obtain Figure 6. The figure confirms that the rotor system became insensitive to the initial condition when increasing the control signal gain ${\eta }_1$ to $0.4$, where the rotor exhibits backward whirling motion only regardless of the initial conditions. It is clear from Figure 7 that the system oscillation amplitude at the $X-$direction is about half of the system oscillation amplitude in $Y-$direction that agrees with the system response curve given in Figure 5E,F.

According to Equations (38) and (41), the equivalent linear damping coefficient ${\mu }_x=\mu +\frac{{\eta }_1{\eta }_3}{2\left({\lambda }_1^2+{\omega }^2\right)}$ is proportional to the product of the control signal gain ${\eta }_1$ and the feedback signal gain ${\eta }_3$ regardless of the magnitude of each one individually. Also, the equivalent linear damping coefficient ${\mu }_y=\mu +\frac{{\eta }_2{\eta }_4}{2\left({\lambda }_2^2+{\omega }^2\right)}$ is proportional to the product of the control gain ${\eta }_2$ and the feedback signal gain ${\eta }_4$ regardless of the magnitude of each one individually. The influence of product of the control and feedback signal gains of the U-IRC (i.e., ${\eta }_1{\eta }_3$) on the rotor system $\sigma -$response curve when ${\eta }_2{\eta }_4=0.0$ is illustrated in Figure 8. The figure shows the system response curves at ${\eta }_1{\eta }_3=0.01, 0.04,$ and $0.08$. Comparing Figure 8A,B with Figure 5A,B, we can note that the two figures are the same, where the product of ${\eta }_1{\eta }_3=0.01$ and ${\eta }_2{\eta }_4=0.0$. Also, Figure 5C,D and Figure 8C,D are identical, where ${\eta }_1{\eta }_3=0.04$ and ${\eta }_2{\eta }_4=0.0$. In addition, Figure 5E,F and Figure 8E,F are the same because ${\eta }_1{\eta }_3=0.08$ and ${\eta }_2{\eta }_4=0.0$. According to Figure 5 and Figure 8, we can confirm that the coupling of an IRC to any dynamical system can mitigate the system’s nonlinear oscillations via adding a linear damping coefficient that is proportional to the product of the feedback and control signal gains of the coupled controller.

The influence of increasing the eccentricity ($f$) on the rotor system oscillation amplitude when activating the U-IRC such that ${\eta }_1{\eta }_3=0.04$ with turning off the V-IRC (i.e., ${\eta }_2{\eta }_4=0.0$) is illustrated in Figure 9. The figure shows the system $f$—response curve at two different values of the rotor spinning speed ($\mathsf{\Omega }=\omega +\sigma , \sigma =0.0, 0.05$). It is clear from the figure that the rotor system responds as a linear system with a single periodic attractor at the small eccentricity magnitudes ($f$). However, the nonlinear behaviors may dominate the system response when $f$ is increased beyond a critical value. Figure 9A,B show that the system may respond with a single periodic solution as long as $f<0.01199$, but the rotor may respond with one of two stable periodic solutions depending on the initial conditions if $0.01199<f<0.02779$. The figure also illustrates that the rotor system can exhibit either a periodic or aperiodic solution depending on the initial conditions if $0.02779<f<0.0366$, but the system performs aperiodic oscillations if the eccentricity is increased beyond $0.0366$. Figure 9C,D illustrate the rotor system $f-$response curve when $\sigma =0.05$. It is clear from the figure the complex dynamics of the rotor system when increasing $f$ from $zero$ to $0.05$, where the system responds with a single periodic solution as long as $f\in [0,0.01245[\cup ]0.01652, 0.02058]$, but the system has bistable periodic solutions as long as $f\in [0.01245, 0.01652[\cup ] 0.02058, 0.044]$. In addition, the rotor system exhibits periodic or aperiodic solution depending on the initial conditions if $f>0.044$.

The stability charts of the 8-pole rotor system are plotted in $f-{\eta }_1{\eta }_3$ plane to show the stable and unstable solution regions when $\sigma =0.0, 0.05$ as illustrated in Figure 10. It is clear from the figure that the increase of ${\eta }_1{\eta }_3$, increases the stable solution region which guarantees the stable oscillations of the system at the large magnitude of the rotor eccentricity. Also, the figure confirms that the rotor system has three oscillation modes when activating the U-IRC controller only. These three modes are (1) stable periodic solution only regardless of the initial conditions, (2) stable or unstable periodic solution depending on the initial conditions, and (3) unstable periodic solution only regardless of the initial conditions. To validate the accuracy of the obtained stability chart in Figure 10, numerical simulations for the system equations of motion (i.e., Equations (12)–(15)) are performed as shown in Figure 11, Figure 12 and Figure 13 according to the three points (${\mathrm{p}}_1,{ \mathrm{p}}_2,{ \mathrm{p}}_3$) that marked in Figure 10B. Figure 11 shows the steady-state temporal oscillations of the rotor system and the connected IRCs according to the point ${\mathrm{p}}_1$ that marked in Figure 10B (i.e., when $f=0.075$, ${\eta }_1{\eta }_3=0.05, {\eta }_2=0.0, {\eta }_4=1$). The figure shows the temporal oscillations and the corresponding Poincaré-map at the initial conditions $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=u\left(0\right)=v\left(0\right)=0$ and $x\left(0\right)=0.3, \dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=u\left(0\right)=v\left(0\right)=0$. It is clear from Figure 11 that both the system and the controllers oscillate with aperiodic motion regardless of the initial conditions that agree with the point ${\mathrm{p}}_1$ that marked in Figure 10B. Figure 12 illustrates the temporal oscillations and the corresponding Poincaré-map of both the rotor system and the IRCs according to the point ${\mathrm{p}}_2$ that showed in Figure 10B (i.e., when $f=0.075$, ${\eta }_1{\eta }_3=0.15, {\eta }_2=0.0, {\eta }_4=1$). It is clear from the figure that the rotor system vibrates with aperiodic motion at zero initial conditions (i.e., $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=u\left(0\right)=v\left(0\right)=0$), but modifying the initial conditions so that $x\left(0\right)=0.3, \dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=u\left(0\right)=v\left(0\right)=0$, has forced the rotor system to oscillate periodically. Accordingly, Figure 12 confirms that the rotor system may respond with a periodic or aperiodic motion depending on the initial conditions that agree with the point ${\mathrm{p}}_2$ that marked in Figure 10B. The nonlinear oscillations of the rotor system and the connected controller are simulated as shown in Figure 13 according to the point ${\mathrm{p}}_3$ that illustrated in Figure 10B (i.e., when $f=0.075$, ${\eta }_1{\eta }_3=0.25, {\eta }_2=0.0, {\eta }_4=1$). The figure is obtained using the two different initial conditions that are used to obtain Figure 11 and Figure 12. It is clear from Figure 13 that the system is insensitive to the initial conditions, where the system responds periodically with the same oscillation amplitude regardless of the initial conditions that agree with the point ${\mathrm{p}}_3$ that is marked in Figure 10B.

b.

The influence of symmetric IRC on the rotor system dynamics

The dynamical behavior of the rotor system when activating both the U-IRC and V-IRC controllers simultaneously is discussed in this section. Figure 14 shows the steady-state oscillation amplitudes of the rotor system ($a \& b$) versus the parameter $\sigma$ at three different values of the product of the control and feedback gains (i.e., when ${\eta }_1{\eta }_3={\eta }_2{\eta }_4=0.01, 0.04, 0.08$). It is clear from the figure that the increase of ${\eta }_1{\eta }_3={\eta }_2{\eta }_4$, increases the system equivalent linear damping coefficients ${\mu }_x$ and ${\mu }_y$ as discussed before, which ultimately decreases the system oscillation amplitudes and eliminates the motion bifurcations as illustrated in Figure 14E,F. Comparing Figure 14 with Figure 8, it is clear from Figure 8E,F that the rotor system responds as a nonlinear system with a bistable periodic attractor at a specific range of the rotor spinning speed when ${\eta }_1{\eta }_3=0.08 \& {\eta }_2{\eta }_4=0.0$, while Figure 14E,F show that the system responds as a linear one with a single periodic attractor regardless of the rotor spinning speed when ${\eta }_1{\eta }_3={\eta }_2{\eta }_4=0.08$. Accordingly, we can deduce that the dynamical behavior of the rotor system when activating one of the IRCs is completely different than the system dynamic when both the U-IRC and V-IRC are activated simultaneously. It is clear from Figure 14 that the worst resonance occurs at $\sigma =0.0$, therefore ${\eta }_1{\eta }_3={\eta }_2{\eta }_4$ are used as a bifurcation parameter to explore the influence of increasing ${\eta }_1{\eta }_3={\eta }_2{\eta }_4$ on the system oscillation amplitudes when $\sigma =0.0$ as illustrated in Figure 15. It is clear from the figure that the system oscillations amplitudes are a monotonic decreasing function of ${\eta }_1{\eta }_3={\eta }_2{\eta }_4$. In addition, the figure confirms the system bistable oscillation amplitudes have been merged into a single periodic attractor when ${\eta }_1{\eta }_3={\eta }_2{\eta }_4$ exceeded a specific value.

It has been concluded before that the coupling of the IRCs to the rotor system modified the rotor linear damping coefficient ($\mu$) to the equivalent damping coefficients ${\mu }_x=\mu +\frac{{\eta }_1{\eta }_3}{2\left({\lambda }_1^2+{\omega }^2\right)}$ and ${\mu }_y=\mu +\frac{{\eta }_2{\eta }_4}{2\left({\lambda }_2^2+{\omega }^2\right)}$. In addition, the detuning parameter ($\sigma$) has been modified to ${\sigma }_x=\sigma +\frac{{\lambda }_1{\eta }_1{\eta }_3}{2\omega \left({\lambda }_1^2+{\omega }^2\right)}$ and ${\sigma }_y=\sigma +\frac{{\lambda }_2{\eta }_2{\eta }_4}{2\omega \left({\lambda }_2^2+{\omega }^2\right)}$. Accordingly, the influence of internal loop feedback gains (${\lambda }_1$ and ${\lambda }_2$) of the IRCs on the oscillation amplitudes ($a \& b$) is illustrated in Figure 16, where Figure 16A,B show the rotor system oscillation amplitudes at ${\lambda }_1={\lambda }_2=0.1, 1.0$ when ${\eta }_1{\eta }_3={\eta }_2{\eta }_4=0.08$. It is clear from Figure 16A,B that the rotor system vibration amplitudes ($a$ and $b)$ when ${\lambda }_1={\lambda }_2=0.1$ is smaller than its oscillation amplitudes when ${\lambda }_1={\lambda }_2=1.0$. It is worthy to mention that the obtained results in Figure 16A,B agree with the derived equations of the equivalent linear damping coefficients ${\mu }_x$ and ${\mu }_y$, where ${\mu }_x$ and ${\mu }_y$ inversely proportional to $2{\lambda }_1^2$ and to $2{\lambda }_2^2$, respectively. Moreover, Figure 16A,B indicate that internal loop feedback gains (${\lambda }_1$ and ${\lambda }_2$) play an important role in shifting the system $\sigma -$response curve to the right or the left via modifying the detuning parameter $\sigma$ to ${\sigma }_x=\sigma +\frac{{\lambda }_1{\eta }_1{\eta }_3}{2\omega \left({\lambda }_1^2+{\omega }^2\right)}$ and ${\sigma }_y=\sigma +\frac{{\lambda }_2{\eta }_2{\eta }_4}{2\omega \left({\lambda }_2^2+{\omega }^2\right)}$. Figure 16C,D shows the rotor system oscillation amplitude at ${\eta }_1{\eta }_3={\eta }_2{\eta }_4=0.08$ when the internal loop feedback gains are increased to ${\lambda }_1={\lambda }_2=2.0$. Comparing Figure 16C,D with Figure 16A,B, we can confirm that the increase of the internal loop feedback gains from ${\lambda }_1={\lambda }_2=1.0$ to ${\lambda }_1={\lambda }_2=2.0$, has decreased equivalent damping coefficients ${\mu }_x$ and ${\mu }_y$, which ultimately forced the rotor system to respond as a nonlinear one with high vibration amplitudes.

The $f-$response curve of the rotor system when ${\eta }_1{\eta }_3={\eta }_2{\eta }_4=0.08$ are illustrated in Figure 17 at two different values of the rotor spinning speed $\mathsf{\Omega }=\omega +\sigma , \sigma =0.0, 0.05$. It is clear from Figure 17 that the controlled rotor system will respond with periodic motion as long as $0.0\le f\le 0.05$. However, the system may respond as a linear system with a single periodic attractor or as a nonlinear one with a bistable periodic attractor depending on the eccentricity magnitude. Comparing Figure 17 with Figure 4, we can confirm that the activation of both the U-IRC and V-IRC simultaneously has eliminated the unstable oscillation of the uncontrolled rotor system that showed in Figure 4. In addition, the four oscillation modes that are noticed in Figure 9 when activating the U-IRC only, have been reduced to two desirable oscillation modes only as illustrated in Figure 17, which are single and bistable periodic motion.

A full picture for the stability boundary in $f-{\eta }_1{\eta }_3={\eta }_2{\eta }_4$ is visualized as shown in Figure 18 when $\sigma =0.0$ and $0.05$. It is clear from the figure the increase of the controller gain (${\eta }_1{\eta }_3={\eta }_2{\eta }_4$), stabilizes the system motion for a wide range of the rotor eccentricity, where the width of the stable solution regions along the $f$ axis is increased exponentially as a function of the control gain ${\eta }_1{\eta }_3={\eta }_2{\eta }_4$. Comparing Figure 18 with Figure 10, one can deduce that the controlled system has simple bifurcation behaviors when activating both the U-IRC and V-IRC simultaneously. In addition, the system can oscillate periodically without losing its stability for a wide range of the excitation force $f$.

To validate the accuracy of the obtained stability charts in Figure 18, numerical simulations for the system temporal equations (i.e., Equations (12)–(15)) are illustrated in Figure 19 and Figure 20 according to the points ${\mathrm{p}}_1$ and ${\mathrm{p}}_2$ that illustrated in Figure 18A, respectively. The temporal oscillation, frequency spectrum, and the corresponding Poincaré-map of the controlled rotor system are illustrated in Figure 19 according to the points ${\mathrm{p}}_1\left(f, {\eta }_1{\eta }_3={\eta }_2{\eta }_4\right)=\left(0.15, 0.1\right)$ that is marked in Figure 18A within the unstable solution region. It is clear from Figure 19 that the rotor system exhibits chaotic oscillations (i.e., aperiodic motion) as predicted in the stability chart given in Figure 18A. In addition, Figure 20 simulates the temporal oscillations, frequency spectrum, and the corresponding Poincaré-map of the rotor system according to the points ${\mathrm{p}}_2\left(f, {\eta }_1{\eta }_3={\eta }_2{\eta }_4\right)=\left(0.15, 0.25\right)$ illustrated in Figure 18A within the stable solution region. By examining Figure 20, one can conclude that the rotor system performs a stable periodic motion that agrees with the stability chart given in Figure 18A with excellent accuracy.

Finally, the controlled system stability chart in $\sigma -{\eta }_1{\eta }_3={\eta }_2{\eta }_4$ has been established as shown in Figure 21 at three different magnitudes of the eccentricity $f$ (i.e., $f=0.05, 0.1,$ and $0.15$). It is clear from the figure that the unstable solutions area lies close to the primary resonance (i.e., close to $\sigma =0.0$). In addition. The figure shows that the increase of the excitation force $f$, increases the unstable solution region. However, one can avoid the instability of the rotor system even at the large eccentricity (i.e., $f=0.15$) if the control gains are adjusted so that ${\eta }_1{\eta }_3={\eta }_2{\eta }_4\ge 0.15$ regardless of the rotor spinning speed $\mathsf{\Omega }=\omega +\sigma$ (See Figure 21C). Numerical validations for the accuracy of the obtained stability charts that are given in Figure 21 have been illustrated in Figure 22 and Figure 23. The nonlinear vibrations of the rotor system according to the point ${\mathrm{p}}_1\left(\sigma , {\eta }_1{\eta }_3={\eta }_2{\eta }_4\right)={\mathrm{p}}_1\left(0, 0.1\right)$ that lies within the unstable solution region in Figure 21B, have been simulated numerically as shown in Figure 22 via solving Equations (14)–(17) using ODE45. By exploring Figure 22, we can find that the rotor system performs an unstable quasiperiodic motion that agrees with the point ${\mathrm{p}}_1$. On the other hand, the temporal oscillations of the considered system according to the point ${\mathrm{p}}_2\left(\sigma , {\eta }_1{\eta }_3={\eta }_2{\eta }_4\right)={\mathrm{p}}_2\left(0, 0.15\right)$ that lies within the stable solution region in Figure 21B, have been illustrated in Figure 23, where the system has responded with periodic oscillation as expected from Figure 21B.

5. Case Study

Within this section, a practical case study for the 8-pole rotor system with the proposed control method (i.e., PD+IRC controller) has been investigated to show how to simulate the dynamics of a physical system based on the obtained dimensionless results given in Section 4. Consider the dimensionless system parameters that are used to obtain Figure 23 and corresponding physical parameters of the eight-pole system given in

Table 1. Actual system parameters and the corresponding dimensionless parameters.

Dimensionless System Parameters for Figure 23		The Corresponding Physical System Parameters That Used to Obtain Figure 24
$f=\frac{e}{s_0}$	$0.1$	Disk radius	$R=0.15 \mathrm{m}$
$p=\frac{s_0}{I_0}{k}_1$	$1.22$	Disk thickness	$D=0.015 \mathrm{m}$
$d=\frac{s_0{\omega }_n}{I_0}{k}_2$	$0.005$	Disk-mass	$m=8 \mathrm{kg}$
${\eta }_1=-8\cos \left(\alpha \right)\frac{s_0}{I_0}{k}_3$	$\sqrt{0.15}$	Disk eccentricity	$e=5\times {10}^{-4} \mathrm{m}$
${\eta }_2=-8\cos \left(\alpha \right)\frac{s_0}{I_0}{k}_4$	$\sqrt{0.15}$	The angle between each two poles	$\alpha ={22.5}^{°}$
${\eta }_3=\frac{k_5}{{\omega }_n}$	$\sqrt{0.15}$	Air-gap size	$s_0=5\times {10}^{-3} \mathrm{m}$
${\eta }_4=\frac{k_6}{{\omega }_n}$	$\sqrt{0.15}$	The magnetic-pole cross-sectional area	$A\cos \left(\theta \right)=7.44\times {10}^{-4} {\mathrm{m}}^2$
${\lambda }_1=\frac{{\rho }_1}{{\omega }_n}$	$1.0$	Coil turn-numbers	$N=1000$
${\lambda }_2=\frac{{\rho }_2}{{\omega }_n}$	$1.0$	Bias current	$I_0=2 \mathrm{A}$
$\mu =4d\cos \left(\alpha \right)$	$0.0185$	Magnetic permeability	${\mu }_0=4\pi \times {10}^{-7} \mathrm{H}/\mathrm{m}$
$\omega =\sqrt{8p\cos \left(\alpha \right)-8}$	$1.008$	Magnetic force constant	$\mathsf{\Delta }=\frac{1}{4}{\mu }_0{N}^{2}A\cos \left(\theta \right)=2.34\times {10}^{-4} \mathrm{Hm}$
$\sigma$	$0$	Normalized natural frequency	${\omega }_n=\sqrt{\mathsf{\Delta }/m{s}_0^3}=15.3 {\mathrm{s}}^{-1}$
$\mathsf{\Omega }$	$\mathsf{\omega }+\sigma$	Proportional control gain	$k_1=488 \mathrm{A}/\mathrm{m}$
		Derivative control gain	$K_2=0.1307 \mathrm{A}.\mathrm{s}/\mathrm{m}$
		IRC control gains	$K_3=K_4=-20.96 \mathrm{A}/\mathrm{m}$
		IRC feedback gains	$K_5=K_6=5.92{ \mathrm{s}}^{-1}$
		IRC internal feedback gains	${\rho }_1={\rho }_2=15.3{ \mathrm{s}}^{-1}$

[41,42]. The actual temporal oscillations (i.e., $x\left(t\right)=s_0{x}^{*}\left(t\right)$ & $y\left(t\right)=s_0{y}^{*}\left(t\right)$) and the corresponding actual control currents ($i_x\left(t\right)=k_{1}x\left(t\right)+k_2\dot{x}\left(t\right)+k_{3}u\left(t\right) \& i_y\left(t\right)=k_{1}y\left(t\right)+k_2\dot{y}\left(t\right)+k_{4}v\left(t\right)$ of the rotor system are simulated in Figure 24, where the actual control current can be obtained based on Equation (6) as follows: $i_x\left(t\right)=k_1{s}_0{x}^{*}\left(t\right)+k_2{s}_0{\omega }_n{\dot{x}}^{*}\left(t\right)+k_3{s}_0{u}^{*}\left(t\right)$ and $i_y\left(t\right)=k_1{s}_0{y}^{*}\left(t\right)+k_2{s}_0{\omega }_n{\dot{y}}^{*}\left(t\right)+k_4{s}_{0}v\left(t\right)$.

Figure 24 shows the actual temporal oscillations of the 8-pole rotor system in meter at $X$ and $Y$ directions, and the corresponding actual control currents in ampere in each electromagnetic pole according to the system physical parameters given in Table 1. It is clear from Figure 24A,B that the rotor lateral oscillations are always smaller than the air gap along the time axis, which confirm the rotational motions of the shaft without any contact with the legs of the electromagnetic poles. Also, Figure 24C,D show that the maximum instantaneous control currents in each pole do not exceed two amperes, which confirm the feasibility of the applied control algorithm in mitigating the rotor lateral oscillation with small control currents.

6. Conclusions

The nonlinear vibrations control of the 8-pole rotor system has been tackled within this article. The integral resonant controller (IRC) has been proposed to control the system lateral vibrations for the first time. The U-IRC and V-IRC controllers have been coupled to the horizontal and vertical oscillation modes of the considered system, respectively. The whole system equations of motion are derived as two second-order nonlinear differential equations (i.e., rotor system) coupled to two first-order linear differential equations (i.e., IRC controllers). Then, the asymptotic analysis is applied to obtain an approximate solution for the derived nonlinear mathematical model. The different bifurcation diagrams have been plotted and the sensitivity analysis for all controller parameters has been explored. In addition, many stability charts in the two-dimensional space have been established to investigate the stability margins of the different control parameters. Moreover, numerical validations for the different response curves have been performed. According to the above discussions, the following important remarks can be summarized:

• The coupling of the U-IRC to the horizontal oscillation mode of the rotor system has modified the system linear damping coefficient $\mu$ to ${\mu }_x=\mu +\frac{{\eta }_1 {\eta }_3}{2\left({\lambda }_1^2+{\omega }^2\right)}.$
• The coupling of the V-IRC to the vertical oscillation mode of the rotor system has modified the system linear damping coefficient $\mu$ to ${\mu }_y=\mu +\frac{{\eta }_2 {\eta }_4}{2\left({\lambda }_2^2+{\omega }^2\right)}.$
• According to the concluded points (1) and (2), the vibration suppression efficiency of the IRC relies on adjusting the linear damping coefficient for the targeted system via designing the optimum values of the control gains (${\eta }_1, {\eta }_2$), feedback gains (${\eta }_3, {\eta }_4$), and internal feedback gains (${\lambda }_1, {\lambda }_2$).
• The coupling of an IRC to a nonlinear oscillatory system is resulting in modifying its linear damping coefficient, where the equivalent linear damping of the controlled system is proportional to the product of the control and feedback gains of the IRC, and inversely proportional to the square of internal loop feedback gain.
• The optimum vibration suppression efficiency of the IRC controller could be achieved via designing its control and feedback gains so that their product is at the maximum possible value, as well as its internal feedback gain should be at the smallest possible value.
• The proper selection of the IRC control parameters can eliminate the catastrophic static bifurcation behaviors of the rotor system and force it to oscillate with negligible vibration amplitudes.