On the Performance of a Nonlinear Position-Velocity Controller to Stabilize Rotor-Active Magnetic-Bearings System

Abstract

The performance of a nonlinear position-velocity controller in stabilising the lateral vibrations of a rotor-active magnetic-bearings system (RAMBS) is investigated. Cubic nonlinear position-velocity and linear position-velocity controllers are introduced to stabilise RAMBS lateral oscillations. According to the proposed control law, the nonlinear system model is established and then investigated with perturbation analysis. Nonlinear algebraic equations that govern the steady-state oscillation amplitudes and the corresponding phases are derived. Depending on the obtained algebraic equations, the different frequency response curves and bifurcation diagrams are plotted for the studied model. Sensitivity analysis for the linear and nonlinear controllers’ gains is explored. Obtained analytical results demonstrated that the studied model had symmetric bifurcation behaviours in both the horizontal and vertical directions. In addition, the integration of the cubic position controller made the control algorithm more flexible to reshape system dynamical behaviours from the hardening spring characteristic to the softening spring characteristic (or vice versa) to avoid resonance conditions. Moreover, the optimal design of the cubic position gain and/or cubic velocity gain could stabilise the unstable motion and eliminate the nonlinear effects of the system even at large disc eccentricities. Lastly, numerical validations for all acquired results are performed, where the presented simulations show accurate correspondence between numerical and analytical investigations.

1. Introduction

An active magnetic bearings system (AMBS) is a bearings system that supports rotating shafts without physical contacts via electromagnetic attractive force. The attractive electromagnetic forces of AMBS can be balanced by using a pair of electromagnetic poles at each opposite directions. These balanced forces allow for the rotor to rotate within the air gap without any frictional losses. The frictionless operation of the rotor AMBS (RAMBS) excludes the necessity for the lubrication of bearings components. Accordingly, AMBS technology provides a clean environment, requires less maintenance, and a long working time for the rotating machines. AMBS is the main part of many high-speed rotating machines such as in aerospace engines, the automotive industry, machine tools, military industries, and autonomous power engineering. Due to the many advantages of the AMBS compared to conventional bearings, vibration analysis and control of the RAMBS are the main global research subject of scientists and engineers, and many configurations with different control strategies of the RAMBS have been investigated. Ji et al. [1] explored the bifurcation behaviours of a 4-pole RAMBS controlled via a linear position-velocity controller. The authors investigated their model as an autonomous and nonautonomous system, and the autonomous system could exhibit saddle-node, Hopf, and saddle-connection bifurcations, while the nonautonomous system could have a homoclinic orbit transversal intersection. Saeed et al. [2,3] discussed the nonlinear dynamical characteristics of a 6-pole RAMBS, and applied the linear position-velocity controller in both Cartesian control configuration [2] and the radial control technique [3]. The authors concluded that the Cartesian control method is more effective than the radial control strategy is in suppressing system vibrations. However, the radial control technique is more stable than the Cartesian control method, especially in the case of large disc eccentricities.

The nonlinear dynamics of an 8-pole RAMBS having a constant stiffness coefficient was extensively investigated [4,5,6,7,8]. Ji and Hansen [4], and Saeed et al. [5] explored the lateral vibration control of the 8-pole RAMBS when the rotor’s angular speed was closed to the system natural frequency utilising the linear position-velocity controller. The authors reported many nonlinear phenomena, including the existence of bistable or tristable solutions at the same angular speed. In addition, they concluded that the proposed control law may lose its stability at large excitation forces. Ji and Leung [6] discussed the lateral oscillations of the 8-pole RAMBS at a superharmonic resonance condition applying the linear position-velocity controller where both bistable and tristable solutions were reported. Yang et al. [7] investigated the modal motion of the 8-pole RAMBS controlled by a linear position-velocity controller. They utilised phase-difference and energy-ratio methods along with the multiple scales perturbation technique to explore the nature of the system orbital motion. The authors reported that the 8-pole RAMBS could perform either elliptic or quasiperiodic motions. The effects of both the impact and rub forces between rotor and stator (i.e., pole legs) on the oscillatory behaviours of the 8-pole RAMBS were investigated by Saeed et al. [8]. The introduced analyses demonstrated that the RAMBS could exhibit either fully annular rub motion or partial rub-impact motion depending on the dynamic friction and impact stiffness coefficients. Zhang and co-workers [9,10,11,12,13,14] studied the 8-pole RAMBS that had a time-varying stiffness coefficient. The authors concluded the system had Shilnikov multipulse chaotic behaviours. In addition, numerical analysis showed that the time-varied RAMBS could perform periodic-$\mathrm{n}$, quasiperiodic, and chaotic vibrations.

The 16-pole RAMBS was investigated with either a constant [15,16] or time-varying [17,18,19,20,21] stiffness coefficient. Saeed et al. [15,16] explored the nonlinear dynamics of the constant stiffness coefficient 16-pole RAMBS. They investigated system oscillatory behaviours in the case of both the Cartesian and radial control configurations utilising the linear position-velocity controller as the main control algorithm. The authors reached the same conclusion as that in [2,3]. Zhang et al. [17,18,19,20,21] investigated 16-pole RAMBS having time-varying stiffness in detail, where many nonlinear phenomena were reported. Saeed et al. [22,23] suggested two advanced control algorithms along with the linear position-velocity controller to enhance the dynamical characteristics of the 8-pole RAMBS. Saeed et al. [22] applied a combination of both linear position-velocity controller and nonlinear saturation controller (NSC) to suppress RAMBS lateral oscillations, while Saeed and Kandil [23] integrated the positive position feedback controller (PPFC) with the linear position-velocity controller to control the 8-pole RAMBS. The authors reported that both the PPFC and NSC were highly efficient in mitigating the nonlinear oscillations of the RAMBS. The contactless feature of the AMBS attracted many researchers to utilise the active magnetic bearings system as an actuator, through which the advanced control algorithms could be applied to modify or to reshape the dynamical characteristics of a wide range of rotating machinery [24,25,26,27,28,29,30,31]. Ishida and Inoue [24] developed a nonlinear vibration absorber utilising 4-pole AMBS to control the nonlinear oscillations of a vertically suspended rotor system. Saeed et al. [25,26,27,28] applied different control schemes of the position and velocity controllers to suppress the lateral oscillations of the rotating shafts using the 4-pole AMBS.

In all previous studies regarding the RAMBS, the Cartesian control methodology was more effective than the radial control strategy was in suppressing system vibrations. However, one of the main drawbacks of Cartesian control is that the system may lose its stability at a large disc eccentricity (i.e., at large excitation force) [2,3,4,5,8,15,16]. Accordingly, the instability problem of the Cartesian control mechanism at a large disc eccentricity was resolved here via integrating a nonlinear position-velocity controller along with the conventional linear position-velocity controller, where the feasibility of the nonlinear position-velocity controller in improving the stability of different dynamical systems was validated [29,30,31,32,33,34,35,36,37,38,39]. Therefore, the system mathematical model was developed and analysed by utilising the multiple-scale perturbation technique according to the suggested control algorithm. The influence of the different control gains was explored. The obtained analytical and numerical investigations showed that the proposed nonlinear controller could stabilise the 8-pole RAMBS at large disc eccentricities and suppress corresponding nonlinear phenomena. In addition, the efficiency of the system in the case of the nonlinear control law was improved compared to that with only the linear position-velocity controller [4,5,6,7,8]. Moreover, the plotted response curves revealed that the considered system had symmetric bifurcation behaviours in both the horizontal and vertical directions.

2. Mathematical Modelling

A schematic diagram of the rotor system is illustrated in Figure 1. The rotating disc was assumed to be a two-degree-of-freedom system supported by controlled magnetic attractive forces $f_i$ ($i=1,2,3,\dots ,8$) generated via eight identical electromagnetic poles. Accordingly, the system equations of motion can be written as follows [40,41].

$$m\ddot{x}(t)=me{\psi }^2\cos (\psi t)+f_X$$

(1)

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

(2)

where $\ddot{x}\left(t\right)$ and $\ddot{y}\left(t\right)$ are the instantaneous lateral acceleration of the rotor in the $X$ and $Y$ directions, respectively; $f_X$ and $f_Y$ are the nonlinear electromagnetic attractive forces that are generated via the magnetic poles in the $X$ and $Y$ directions, respectively; $m$ is the mass of the rotating disc; $e$ is disc eccentricity; and $\psi$ is disc angular velocity. All poles were assumed to be identical. Therefore, the attractive electromagnetic force of each pole can be written as follows [41].

$$f_i=\frac{1}{4}{\mu }_0{N}^{2}A\frac{I_i^2}{{\delta }_i^2}\cos (\phi ),\hspace{1em}\hspace{1em}i=1,2,\cdots ,8$$

(3)

where ${\mu }_0$ is the air-gap permeability, $N$ is the turn number of each electrical coil suspended on each pole, $A\cos \left(\phi \right)$ is the effective cross-sectional area, $I_i$ is the electrical current in the $i$th pole, and ${\delta }_i$ is the instantaneous radial clearance between the rotating disc and the $i$th pole. According to the studied system geometry, for displacements $x$ and $y$ of the rotating disc away from origin $O,$ as shown in Figure 1b, radial clearance ${\delta }_i$ can be expressed as

$$\begin{array}{l}{\delta }_i=c_0\pm x\sin (\alpha )\mp y\cos (\alpha ),\hspace{1em}\hspace{1em}i=1,5\\ {\delta }_i=c_0\pm x\sin (\alpha )\pm y\cos (\alpha ),\hspace{1em}\hspace{1em}i=4,8\\ {\delta }_i=c_0\pm x\cos (\alpha )\mp y\sin (\alpha ),\hspace{1em}\hspace{1em}i=2,6\\ {\delta }_i=c_0\pm x\cos (\alpha )\pm y\sin (\alpha ),\hspace{1em}\hspace{1em}i=3,7\hspace{1em}\hspace{1em}\end{array}\}$$

(4)

where $c_0$ is the nominal air-gap size as shown in Figure 1a and Figure 2a is the angle between any two successive poles.

Electrical currents ($I_1, I_2, \dots , I_8$) in the eight electrical coils were designed such that

$$I_1=I_8=I_0-i_y,\hspace{1em}{I}_2=I_3=I_0+i_x,\hspace{1em}{I}_4=I_5=I_0+i_y,\hspace{1em}{I}_6=I_7=I_0-i_x.$$

(5)

where $i_x$ and $i_y$ are control currents that are responsible for restoring the rotating disc to its nominal position when any deviation occurs due to disc eccentricity, while $I_0$ is a constant electrical current known as premagnetising current. In most work regarding RAMBS, the linear position-velocity controller was only applied to suppress the system’s nonlinear oscillations [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]. However, both the cubic-position and cubic-velocity controllers proved their feasibility and applicability in controlling the dynamical behaviours of a wide range of nonlinear systems [29,30,31,32,33,34,35,36,37,38,39]. Accordingly, a combination of the linear and cubic position-velocity controllers is suggested here to control the nonlinear vibrations of the considered system. Therefore, control currents $i_x$ and $i_y$ are proposed, such that

$$i_x=k_{1}x+k_2{x}^3+k_3\dot{x}+k_4{\dot{x}}^3,\hspace{1em}\hspace{1em}\hspace{1em}{i}_y=k_{1}y+k_2{y}^3+k_3\dot{y}+k_4{\dot{y}}^3.$$

(6)

where $k_1$ and $k_2$ denote linear and cubic position control gains, while $k_3$ and $k_4$ represent linear and cubic velocity control gains, respectively. By substituting Equations (4)–(6) into Equation (3) by expanding the resulting equations using the Maclaurin series up to the third order, we can obtain forces $f_i, \left(i=1,2,\dots ,8\right)$ as given in Appendix A. According to the system geometry in Figure 1, the resultant attractive forces ($f_X$ and $f_Y$) of the eight magnetic poles can be expressed, such that

$$f_X=(f_6+f_7-f_2-f_3)\cos (\alpha )+(f_5+f_8-f_1-f_4)\sin (\alpha )$$

(7)

$$f_Y=(f_1+f_8-f_4-f_5)\cos (\alpha )+(f_2+f_7-f_3-f_6)\sin (\alpha )$$

(8)

By substituting Equations (7) and (8) into Equations (1) and (2), taking into account the forces $f_i, \left(i=1,2,\dots ,8\right)$ as given in Appendix A, we obtain

$$\begin{array}{l}m\ddot{x}(t)-\frac{1}{4c_0^5}{\mu }_0{N}^{2}A\cos (\phi )[\left(8I_0^2{c}_0^2-8I_0{c}_0\cos (\alpha )\right)x-\left(8I_0{c}_0^3{k}_3\cos (\alpha )\right)\dot{x}+(16I_0^2-32I_0^2{\cos }^2(\alpha )\\ -16I_0{c}_0^3{k}_2{\cos }^2(\alpha )-24I_0{c}_0{k}_1{\cos }^3(\alpha )+64I_0^2{\cos }^4(\alpha )+8c_0^2{k}_1^2{\cos }^2(\alpha ))x^3+(8c_0^2{k}_1^2-8c_0^2{k}_1^2{\cos }^2(\alpha )\\ +96I_0^2{\cos }^2(\alpha )-96I_0^2{\cos }^4(\alpha )-72I_0{c}_0{k}_1\cos (\alpha )+72I_0{c}_0{k}_1{\cos }^3(\alpha ))x{y}^2+(16k_1{k}_2{c}_0^2{\cos }^2(\alpha )\\ -24I_0{k}_3{c}_0{\cos }^3(\alpha ))x^2\dot{x}+\left(8k_3^2{c}_0^2{\cos }^2(\alpha )\right)x{\dot{x}}^2+(24I_0{k}_3{c}_0{\cos }^3(\alpha )-24I_0{k}_3{c}_0\cos (\alpha ))\dot{x}{y}^2\\ +\left(8k_3^2{c}_0^2-8k_3^2{c}_0^2{\cos }^2(\alpha )\right)x{\dot{y}}^2+(16k_1{k}_3{c}_0^2-16k_1{k}_3{c}_0^2{\cos }^2(\alpha )-48I_0{k}_3{c}_0\cos (\alpha )\\ +48I_0{k}_3{c}_0{\cos }^3(\alpha ))xy\dot{y}-\left(8I_0{k}_4{c}_0^3\cos (\alpha )\right){\dot{x}}^3]=me{\psi }^2\cos (\psi t)\end{array}$$

(9)

$$\begin{array}{l}m\ddot{y}(t)-\frac{1}{4c_0^5}{\mu }_0{N}^{2}A\cos (\phi )[\left(8I_0^2{c}_0^2-8I_0{c}_0\cos (\alpha )\right)y+(-8I_0{c}_0^3{k}_3\cos (\alpha ))\dot{y}+(16I_0^2-32I_0^2{\cos }^2(\alpha )\\ -16I_0{c}_0^3{k}_2{\cos }^2(\alpha )-24I_0{c}_0{k}_1{\cos }^3(\alpha )+32I_0^2{\cos }^4(\alpha )+8c_0^2{k}_1^2{\cos }^2(\alpha ))y^3+{(8c}_0^2{k}_1^2-8c_0^2{k}_1^2{\cos }^2(\alpha )\\ +96I_0^2{\cos }^2(\alpha )-96I_0^2{\cos }^4(\alpha )-72I_0{c}_0{k}_1\cos (\alpha )+72I_0{c}_0{k}_1{\cos }^3(\alpha ))y{x}^2+(16k_1{k}_2{c}_0^2{\cos }^2(\alpha )\\ -24I_0{k}_3{c}_0{\cos }^3(\alpha ))y^2\dot{y}+\left(8k_3^2{c}_0^2{\cos }^2(\alpha )\right)y{\dot{y}}^2+\left(24I_0{k}_3{c}_0{\cos }^3(\alpha )-24I_0{k}_3{c}_0\cos (\alpha )\right)\dot{x}{y}^2\\ +\left(8k_3^2{c}_0^2-8k_3^2{c}_0^2{\cos }^2(\alpha )\right)y{\dot{x}}^2+(16k_1{k}_3{c}_0^2-16k_1{k}_3{c}_0^2{\cos }^2(\alpha )-48I_0{k}_3{c}_0\cos (\alpha )\\ +48I_0{k}_3{c}_0{\cos }^3(\alpha ))yx\dot{x}+\left(-8I_0{k}_4{c}_0^3\cos (\alpha )\right){\dot{y}}^3]=me{\psi }^2\sin (\psi t)\end{array}$$

(10)

Introducing dimensionless variables and parameters $x^{*}=\frac{x}{c_0}, y^{*}=\frac{y}{c_0}$, $t^{*}={\omega }_{n}t, {\dot{x}}^{*}=\dot{\frac{x}{{\omega }_n{c}_0}}$, ${\dot{y}}^{*}=\dot{\frac{y}{{\omega }_n{c}_0}}, {\ddot{x}}^{*}=\frac{\ddot{x}}{{\omega }_n^2{c}_0}, {\ddot{y}}^{*}=\frac{\ddot{y}}{{\omega }_n^2{c}_0}, p=\frac{c_0}{I_0}{k}_1, d=\frac{c_0{\omega }_n}{I_0}{k}_3, {\lambda }_1=\frac{c_0^3}{I_0}{k}_2, {\lambda }_2=\frac{c_0{\omega }_n^3}{I_0}{k}_4 , f=\frac{e}{c_0}, \mathsf{\Omega }=\frac{\omega }{{\omega }_n}$, ${\omega }_n=\sqrt{\frac{{\mu }_0{I}_0^2{N}^{2}A\cos \left(\phi \right)}{4m{c}_0^3}}$ into Equations (9) and (10) while omitting the asterisks for brevity, we obtain the following normalised equations of motion

$$\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}+{\alpha }_8{\dot{x}}^3)={\mathrm{\Omega }}^{2}f\cos (\mathrm{\Omega }t)$$

(11)

$$\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}+{\alpha }_8{\dot{y}}^3)={\mathrm{\Omega }}^{2}f\sin (\mathrm{\Omega }t)$$

(12)

where $\mu =4d\cos (\alpha ),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\omega }^2=8p\cos (\alpha )-8$,                                        $\begin{array}{l}{\alpha }_1=16-32{\cos }^2(\alpha )+32{\cos }^4(\alpha )-24p{\cos }^3(\alpha )+8p^2{\cos }^2(\alpha )-8{\lambda }_1\cos (\alpha ),\\ {\alpha }_2=8p^2-8p^2{\cos }^2(\alpha )+96{\cos }^2(\alpha ){\sin }^2(\alpha )-72p\cos (\alpha ){\sin }^2(\alpha ),\\ {\alpha }_3=16pd{\cos }^2(\alpha )-24d{\cos }^3(\alpha ),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} {\alpha }_4=-24d\cos (\alpha ){\sin }^2(\alpha ),\\ {\alpha }_5=8d^2{\sin }^2(\alpha ),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\alpha }_6=8d^2{\cos }^2(\alpha ),\\ {\alpha }_7=16pd{\sin }^2(\alpha )-48d\cos (\alpha ){\sin }^2(\alpha ),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\alpha }_8=8{\lambda }_2\cos (\alpha ).\end{array}$

3. System Periodic Solution and Slow-Flow Modulating Equations

This section is devoted to obtaining an analytical solution for the system equations of motion (i.e., Equations (11) and (12)). As the system model is governed by two coupled nonlinear differential equations and only comprises cubic nonlinearity, the first-order multiple-time scale solution was assumed to be as follows [42,43].

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

(13)

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

(14)

where $\epsilon$ is a book-keeping parameter only [41], and $T_0=t$ and $T_1=\epsilon t$ are the fast and slow time scales, respectively. In terms of $T_0$ and $T_1$, ordinary derivatives $\frac{d}{dt}$ and $\frac{d^2}{d{t}^2}$ can be expressed 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{1em}{D}_j=\frac{\partial }{\partial T_j},\hspace{1em}\hspace{1em}j=0,1$$

(15)

To perform the multiple-time-scale procedure, system parameters were scaled, such that

$$\mu =\epsilon \widehat{\mu },\hspace{1em}{\alpha }_j=\epsilon {\widehat{\alpha }}_j,\hspace{1em}f=\epsilon \widehat{f},\hspace{1em}j=1,\dots ,8$$

(16)

By substituting Equations (13)–(16) into Equations (11) and (12), and equating the coefficients of the same power of $\epsilon$, we have

O (${\epsilon }^0$)

$${(D}_0^2+\omega {}_{}^2)x_0=0$$

(17)

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

(18)

O ($\epsilon$)

$$\begin{array}{cc}\hfill \left(D_0^2+{\omega }^2\right)x_1& =-2D_0{D}_1{x}_0-2\widehat{\mu }{D}_0{x}_0+{\widehat{\alpha }}_1{x}_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{\left(D_0{y}_0\right)}^2\hfill \\ & +{\widehat{\alpha }}_6{x}_0{\left(D_0{x}_0\right)}^2+{\widehat{\alpha }}_7{x}_0{y}_0{D}_0{y}_0+{\widehat{\alpha }}_8{\left(D_0{x}_0\right)}^3+\widehat{f}{\mathrm{\Omega }}^2\cos (\mathrm{\Omega }{T}_0)\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \left(D_0^2+{\omega }^2\right)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{\left(D_0{x}_0\right)}^2\hfill \\ & +{\widehat{\alpha }}_6{y}_0{\left(D_0{y}_0\right)}^2+{\widehat{\alpha }}_7{y}_0{x}_0{D}_0{x}_0+{\widehat{\alpha }}_8{\left(D_0{y}_0\right)}^3+\widehat{f}{\mathrm{\Omega }}^2\sin (\mathrm{\Omega }{T}_0)\hfill \end{array}$$

(20)

The solutions of Equations (17) and (18) are

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

(21)

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

(22)

where $A\left(T_1\right)$ and $B\left(T_1\right)$ are unknown functions, while $\overline{A}\left(T_1\right)$ and $\overline{B}\left(T_1\right)$ are the complex conjugate of $A\left(T_1\right)$ and $B\left(T_1\right)$, respectively. Substituting Equations (21) and (22) into Equations (19) and (20), we have

$$\begin{array}{cc}\hfill \left(D_0^2+{\omega }^2\right)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}\hfill \\ & +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}\hfill \\ & +i{\widehat{\alpha }}_7\omega \overline{A}{B}^2+3i{\widehat{\alpha }}_8{\omega }^3{A}^2\overline{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\hfill \\ & -{\widehat{\alpha }}_5{\omega }^{2}A{B}^2-{\widehat{\alpha }}_6{\omega }^2{A}^3+i{\widehat{\alpha }}_7\omega A{B}^2-i{\widehat{\alpha }}_8{\omega }^3{A}^3]e^{3i\omega T_0}+\frac{1}{2}{\mathrm{\Omega }}^2\widehat{f}{e}^{i\mathrm{\Omega }{T}_0}\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \left(D_0^2+{\omega }^2\right)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}\hfill \\ & +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{1em}+{\widehat{\alpha }}_6{\omega }^2{B}^2\overline{B}\hfill \\ & +i{\widehat{\alpha }}_7\omega \overline{B}{A}^2+3i{\widehat{\alpha }}_8{\omega }^3{B}^2\overline{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\hfill \\ & -{\widehat{\alpha }}_5{\omega }^{2}B{A}^2-{\widehat{\alpha }}_6{\omega }^2{B}^3+i{\widehat{\alpha }}_7\omega B{A}^2-i{\widehat{\alpha }}_8{\omega }^3{B}^3]e^{3i\omega T_0}-\frac{1}{2}i{\mathrm{\Omega }}^2\widehat{f}{e}^{i\mathrm{\Omega }{T}_0}\hfill \end{array}$$

(24)

Equations (23) and (24) show that the possible resonance cases are the primary resonance (i.e., $\mathsf{\Omega }=\omega$) and subharmonic resonance cases (i.e., $\mathsf{\Omega }=3\omega$. The efficiency of the applied nonlinear controller is examined at the primary resonance here. Accordingly, the closeness of rotor angular speed $\mathsf{\Omega }$ to RAMBS natural frequency $\omega$ can be expressed via introducing parameter $\sigma$, such that

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

(25)

By substituting Equation (25) into the small devisor and secular terms of Equations (23) and (24), we have

$$\begin{array}{cc}\hfill (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}\hfill \\ & +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}\hfill \\ & +i{\widehat{\alpha }}_7\omega \overline{A}{B}^2+3i{\widehat{\alpha }}_8{\omega }^3{A}^2\overline{A}+\frac{1}{2}{(\omega +\sigma )}^2\widehat{f}{e}^{i\epsilon \widehat{\sigma }{T}_0}]e^{i\omega T_0}+[{\widehat{\alpha }}_1{A}^3+{\widehat{\alpha }}_{2}A{B}^2\hfill \\ & +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-i{\widehat{\alpha }}_8{\omega }^3{A}^3]e^{3i\omega T_0}\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill (D_0^2+{\omega }^2)y_1& =[-2i\omega (D_{1}B)-2i\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}\hfill \\ & +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{1em}+{\widehat{\alpha }}_6{\omega }^2{B}^2\overline{B}\hfill \\ & +i{\widehat{\alpha }}_7\omega \overline{B}{A}^2+3i{\widehat{\alpha }}_8{\omega }^3{B}^2\overline{B}-\frac{1}{2}i{(\omega +\sigma )}^2\widehat{f}{e}^{i\epsilon \widehat{\sigma }{T}_0}]e^{i\omega T_0}+[{\widehat{\alpha }}_1{B}^3+{\widehat{\alpha }}_{2}B{A}^2\hfill \\ & +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-i{\widehat{\alpha }}_8{\omega }^3{B}^3]e^{3i\omega T_0}\hfill \end{array}$$

(27)

The solvability conditions of Equations (26) and (27) are the vanishing of the coefficients of $e^{i\omega T_0}$. Therefore, we have

$$\begin{array}{l}-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+3i{\widehat{\alpha }}_8{\omega }^3{A}^2\overline{A}+\frac{1}{2}{(\omega +\sigma )}^2\widehat{f}{e}^{i\epsilon \widehat{\sigma }{T}_0}=0\end{array}$$

(28)

$$\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{\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+{\widehat{\alpha }}_6{\omega }^2{B}^2\overline{B}+i{\widehat{\alpha }}_7\omega \overline{B}{A}^2+3i{\widehat{\alpha }}_8{\omega }^3{B}^2\overline{B}-\frac{1}{2}i{(\omega +\sigma )}^2\widehat{f}{e}^{i\epsilon \widehat{\sigma }{T}_0}=0\end{array}$$

(29)

To obtain the solution of Equations (28) and (29), it is convenient to express $A\left(T_1\right)$ and $B\left(T_1\right)$ in their polar form, such that

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

(30)

where $a\left(T_1\right)$ and $b\left(T_1\right)$ are the steady-state oscillation amplitudes of RAMBS in the horizontal and vertical directions, respectively, while ${\theta }_1$ and ${\theta }_2$ are the phases. Substituting Equation (30) into Equations (28) and (29) by separating the real and imaginary parts, and restoring the scaled parameters to the original form (i.e., $T_1=\epsilon t,\sigma =\epsilon \widehat{\sigma },\mu =\epsilon \widehat{\mu },f=\epsilon \widehat{f},{\alpha }_j=\epsilon {\widehat{\alpha }}_j, j=1,\dots ,8$), we obtain

$$\begin{array}{cc}\hfill \frac{da}{dt}=H_1(a,b,{\phi }_1,{\phi }_2)& =-\mu a+\frac{1}{8}{a}^3\left({\alpha }_3+3{\alpha }_8{\omega }^2\right)+\frac{1}{4}{\alpha }_{4}a{b}^2+\frac{1}{8}a{b}^2\left(-{\alpha }_4+{\alpha }_7\right)\cos (2{\phi }_2-2{\phi }_1)\hfill \\ & \frac{1}{8\omega }a{b}^2\left({\alpha }_2-{\alpha }_5{\omega }^2\right)\sin (2{\phi }_2-2{\phi }_1)+\frac{1}{2\omega }{(\omega +\sigma )}^{2}f\sin {\phi }_1\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \frac{db}{dt}=H_2(a,b,{\phi }_1,{\phi }_2)& =-\mu b+\frac{1}{8}({\alpha }_3+3{\alpha }_8{\omega }^2)b^3+\frac{1}{4}{\alpha }_{4}b{a}^2+\frac{1}{8}(-{\alpha }_4+{\alpha }_7)b{a}^2(\cos (2{\phi }_2-2{\phi }_1)\hfill \\ & +\frac{1}{8}(\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega )b{a}^2\sin (2{\phi }_2-2{\phi }_1)-\frac{1}{2\omega }{(\omega +\sigma )}^{2}f\cos {\phi }_2\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \frac{d{\phi }_1}{dt}=H_3(a,b,{\phi }_1,{\phi }_2)& =\sigma +\frac{1}{8\omega }{a}^2\left(3{\alpha }_1+{\alpha }_6\omega \right)+\frac{1}{4\omega }{b}^2\left({\alpha }_2+{\alpha }_5{\omega }^2\right)+\frac{1}{8\omega }{b}^2\left({\alpha }_2-{\alpha }_5{\omega }^2\right)\cos (2{\phi }_2-2{\phi }_1)\hfill \\ & +\frac{1}{8}{b}^2\left(-{\alpha }_4+{\alpha }_7\right)\sin (2{\phi }_2-2{\phi }_1)+\frac{1}{2a\omega }f{\left(\omega +\sigma \right)}^2\cos {\phi }_1\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \frac{d{\phi }_2}{dt}=H_4(a,b,{\phi }_1,{\phi }_2)& =\sigma +\frac{1}{8}(\frac{3{\alpha }_1}{\omega }+{\alpha }_6\omega )b^2+\frac{1}{4}(\frac{{\alpha }_2}{\omega }+{\alpha }_5\omega )a^2+\frac{1}{8}(\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega )a^2(\cos (2{\phi }_2-2{\phi }_1)\hfill \\ & \frac{1}{8}(-{\alpha }_4+{\alpha }_7)a^2\sin (2{\phi }_2-2{\phi }_1)+\frac{1}{2b\omega }{(\omega +\sigma )}^{2}f\sin {\phi }_2\hfill \end{array}$$

(34)

where ${\phi }_1=\sigma t-{\theta }_1, {\phi }_2=\sigma t-{\theta }_2$. Substituting Equations (21), (22), and (30) into Equations (13) and (14), we have

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

(35)

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

(36)

where Equations (35) and (36) are the solution of Equations (11) and (12). At steady-state conditions, we have $\frac{da}{dt}=\frac{db}{dt}=\frac{d{\phi }_1}{dt}=\frac{d{\phi }_2}{dt}=0$. Substituting this condition into Equations (31)–(34), we obtain the following nonlinear system of algebraic equations:

$$\begin{array}{l}-\mu a+\frac{1}{8}({\alpha }_3+3{\alpha }_8{\omega }^2)a^3+\frac{1}{4}{\alpha }_{4}a{b}^2+\frac{1}{8}(-{\alpha }_4+{\alpha }_7)a{b}^2\cos (2{\phi }_2-2{\phi }_1)\\ -\frac{1}{8}(\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega )a{b}^2\sin (2{\phi }_2-2{\phi }_1)+\frac{1}{2\omega }{(\omega +\sigma )}^{2}f\sin ({\phi }_1)=0\end{array}$$

(37)

$$\begin{array}{l}-\mu b+\frac{1}{8}({\alpha }_3+3{\alpha }_8{\omega }^2)b^3+\frac{1}{4}{\alpha }_{4}b{a}^2+\frac{1}{8}(-{\alpha }_4+{\alpha }_7)b{a}^2\cos (2{\phi }_2-2{\phi }_1)\\ +\frac{1}{8}(\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega )b{a}^2\sin (2{\phi }_2-2{\phi }_1)-\frac{1}{2\omega }{(\omega +\sigma )}^{2}f\cos ({\phi }_2)=0\end{array}$$

(38)

$$\begin{array}{l}-\sigma -\frac{1}{8}(\frac{3{\alpha }_1}{\omega }+{\alpha }_6\omega )a^2-\frac{1}{4}(\frac{{\alpha }_2}{\omega }+{\alpha }_5\omega )b^2-\frac{1}{8}(\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega )b^2\cos (2{\phi }_2-2{\phi }_1)\\ -\frac{1}{8}(-{\alpha }_4+{\alpha }_7)b^2\sin (2{\phi }_2-2{\phi }_1)-\frac{1}{2a\omega }{(\omega +\sigma )}^{2}f\cos ({\phi }_1)=0\end{array}$$

(39)

$$\begin{array}{l}-\sigma -\frac{1}{8}(\frac{3{\alpha }_1}{\omega }+{\alpha }_6\omega )b^2-\frac{1}{4}(\frac{{\alpha }_2}{\omega }+{\alpha }_5\omega )a^2-\frac{1}{8}(\frac{{\alpha }_2}{\omega }-{\alpha }_5\omega )a^2\cos (2{\phi }_2-2{\phi }_1)\\ +\frac{1}{8}(-{\alpha }_4+{\alpha }_7)a^2\sin (2{\phi }_2-2{\phi }_1)-\frac{1}{2b\omega }{(\omega +\sigma )}^{2}f\sin ({\phi }_2)=0\end{array}$$

(40)

By solving the above algebraic equations in terms of the system and control parameters (i.e., $\sigma , f, p, d, {\lambda }_1, {\lambda }_2$), we can obtain the different bifurcation diagrams as given in Section 4. In addition, the stability of the obtained solution (i.e., the solution of Equations (37)–(40)) can be investigated by checking the eigenvalues of the Jacobian matrix of the dynamical system given by Equations (31)–(34). Accordingly, let the solution of Equations (37)–(40) be ($a_0, b_0,{\phi }_{10},{\phi }_{20}$), and suppose ($a_1, b_1,{\phi }_{11},{\phi }_{21}$) is a small perturbation about this solution. So, we can write the following:

$$\begin{array}{r}a=a_0+a_1, b=b_0+b_1, {\phi }_1={\phi }_{10}+{\phi }_{11}, {\phi }_2={\phi }_{20}+{\phi }_{21},\\ \frac{da}{dt}=\frac{d{a}_1}{dt},\hspace{1em}\frac{db}{dt}=\frac{d{b}_1}{dt},\hspace{1em}\frac{d{\phi }_1}{dt}=\frac{d{\phi }_{11}}{dt},\hspace{1em}\frac{d{\phi }_2}{dt}=\frac{d{\phi }_{21}}{dt}\end{array}\}$$

(41)

By inserting Equation (41) into Equations (31)–(34), we can obtain the linearised model of the system given by Equations (31)–(34) as

$$\left(\begin{array}{c}\frac{d{a}_1}{dt}\\ \frac{d{b}_1}{dt}\\ \frac{d{\phi }_{11}}{dt}\\ \frac{d{\phi }_{21}}{dt}\end{array}\right)=\left(\begin{array}{cccc}\frac{\partial H_1}{\partial a_1}& \frac{\partial H_1}{\partial b_1}& \frac{\partial H_1}{\partial {\phi }_{11}}& \frac{\partial H_1}{\partial {\phi }_{21}}\\ \frac{\partial H_2}{\partial a_1}& \frac{\partial H_2}{\partial b_1}& \frac{\partial H_2}{\partial {\phi }_{11}}& \frac{\partial H_2}{\partial {\phi }_{21}}\\ \frac{\partial H_3}{\partial a_1}& \frac{\partial H_3}{\partial b_1}& \frac{\partial H_3}{\partial {\phi }_{11}}& \frac{\partial H_3}{\partial {\phi }_{21}}\\ \frac{\partial H_4}{\partial a_1}& \frac{\partial H_4}{\partial b_1}& \frac{\partial H_4}{\partial {\phi }_{11}}& \frac{\partial H_4}{\partial {\phi }_{21}}\end{array}\right)\left(\begin{array}{c}{a}_1\\ b_1\\ {\phi }_{11}\\ {\phi }_{21}\end{array}\right)$$

(42)

The coefficients of the Jacobian matrix given by Equation (42) are given in Appendix B. According to the Hartman–Grobman theorem [44], nonlinear autonomous system (31)–(34) is topologically equivalent to linear system (42) at the hyperbolic equilibrium point ($a_0,b_0, {\varphi }_{10},{\varphi }_{20},$). Hence, the solution of the nonlinear system given by Equations (31)–(34) is asymptotically stable if and only if the eigenvalues of the Jacobian matrix in (42) have a real negative part.

4. Sensitivity Investigations

Within this section, the different response curves of the RAMBS are obtained via solving the nonlinear algebraic Equations (37)–(40) numerically applying the Newton–Raphson algorithm with a continuation method, using parameters $\sigma$, $f$, ${\lambda }_1,$ and ${\lambda }_2$ as bifurcation control parameters [45,46]. The sensitivity of the system vibration amplitudes to the change in control parameters $p, d, {\lambda }_1$ and ${\lambda }_2$ was investigated. The obtained bifurcation diagrams are shown as a solid line for stable solutions, and a dotted line for unstable solutions. Moreover, numerical confirmations for the plotted response curves were introduced by solving system temporal Equations (11) and (12), utilising the ODE45 MATLAB solver. Numerical results are plotted as a small circle during the increment of the bifurcation parameter, and as a large dot during the decrement of the bifurcation parameter. Simulation results were established using the following system parameters: $p=1.22, d=0.005, \alpha ={22.5}^{°}$, ${\lambda }_1={\lambda }_2=0.0,$ $f=0.015$, and $\mathsf{\Omega }=\omega +\sigma$ unless otherwise mentioned [4,5,6,7,8]. Dimensionless parameters $p, d, {\lambda }_1$, and ${\lambda }_2$ are defined such that $p=\frac{c_0}{I_0}{k}_1, d=\frac{c_0{\omega }_n}{I_0}{k}_3, {\lambda }_1=\frac{c_0^3}{I_0}{k}_2, {\lambda }_2=\frac{{\mathrm{c}}_0{\mathsf{\omega }}_{\mathrm{n}}^3}{I_0}{k}_4,$ as given in Equation (10). Accordingly, $p$ and $d$ denote the dimensionless linear-position and linear-velocity control gains, respectively. In addition, ${\lambda }_1$ and ${\lambda }_2$ represent the dimensionless cubic-position and cubic-velocity control gains, respectively (Equation (6)). In the following subsections, the efficiency of the linear position-velocity and cubic position-velocity controllers in controlling the oscillation amplitudes ($a$ and $b$) of the RAMBS is explored by solving Equations (37)–(40) in terms of control gains ($p, d, {\lambda }_1, {\lambda }_2$), disc eccentricity ($f$), and disc spinning speed ($\mathsf{\Omega }=\omega +\sigma$).

4.1. Sensitivity Analysis of Linear Position-Velocity Controller ($p$ and $d$)

The performance of the linear position-velocity controller only (i.e., ${\lambda }_1={\lambda }_2=0$) in eliminating the vibrations of the RAMBS is investigated here. According to Equation (25), if $\sigma =0$, the system works at perfect primary resonance (i.e., $\mathsf{\Omega }=\omega$); $\sigma >0$ implies that the disc spinning speed ($\mathsf{\Omega }$) is greater than the system natural frequency ($\omega$) by value $\sigma$. In addition, $\sigma <0$ implies that the disc spinning speed ($\mathsf{\Omega }$) is smaller than the system natural frequency ($\omega$) by value $\sigma$. In other words, parameter $\sigma$ represents the closeness of RAMBS natural frequency ($\omega$) to disc angular speed ($\mathsf{\Omega }$). Therefore, we employed $\sigma$ as a bifurcation control parameter denoting the angular-speed response curve.

Accordingly, the angular-speed response curves for RAMBS at two different magnitudes of disc eccentricity (i.e., $f=0.05$ and $0.015$) are presented in Figure 2. Figure 2a,b show that the RAMBS responded as a linear system, exhibiting single-valued vibration amplitudes for every excitation frequency $\mathsf{\Omega }$ as long as the disc eccentricity had small values (i.e., $f=0.005$). However, Figure 2c,d show that nonlinear behaviours, such as the coexistence of bi- or tri-stable solutions and sensitivity to the initial conditions, dominated the spinning-speed response curves when the rotating disc had large eccentricity.

The sensitivity of RAMBS oscillation amplitudes ($a, b$) to the position control gain ($p$) when the other control parameters were fixed constants (i.e., $d=0.005, {\lambda }_1={\lambda }_2=0.0$) is illustrated in Figure 3. Figure 3a,b show the RAMBS response curves at $p=1.35$, while Figure 3c,d illustrate the response curves at $p=1.45$. Comparing Figure 3a,b with Figure 3c,d shows that the increase in position gain from $p=1.35$ to $p=1.45$ increased the system natural frequency (i.e., $\omega =\sqrt{8p\cos \left(\alpha \right)-8}$), which ultimately shifted the response curves to the right. Figure 3 shows that increasing the proportional gain enhanced RAMBS vibrational behaviours at the perfect tuning (i.e., when $\sigma =0$), where the system exhibited bistable vibration amplitudes at $p=1.35$ that were merged into a single solution when $p=1.45$. On the other hand, Figure 3 confirms that increasing $p$ may increase the system lateral vibrations when the RAMBS operates at a spinning speed ($\mathsf{\Omega }$) higher than the natural frequency ($\omega$) (i.e., when $\sigma =\mathsf{\Omega }-\omega >0$).

According to Figure 3 at $\sigma =0.1$, numerical simulations for the temporal lateral oscillations ($x\left(t\right)$ & $y\left(t\right)$) and the corresponding orbital motion of the RAMBS are illustrated in Figure 4, Figure 5 and Figure 6 at three different values of the initial conditions when $p=1.35$ and $1.45$. Figure 4 shows the RAMBS vibration amplitudes in $X$ and $Y$, and the corresponding orbital motion according to Figure 3 when $\sigma =0.1$ (i.e., when $\mathsf{\Omega }=\omega +0.1$) via numerically solving Equations (11) and (12) using the ODE45 MATLAB algorithm at initial conditions $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=0$. Figure 4c shows that the system executed a forward whirling motion, where the whirling amplitude at $p=1.45$ was larger than that at $p=1.35$. Figure 5 is a repetition of Figure 4 at initial conditions $x\left(0\right)=\dot{x}\left(0\right)=0.5, y\left(0\right)=\dot{y}\left(0\right)=0$, while Figure 6 is a repetition of Figure 4 at initial conditions $x\left(0\right)=\dot{x}\left(0\right)=0, y\left(0\right)=\dot{y}\left(0\right)=-1.9$. Figure 5c and Figure 6c show that the RAMBS could execute a forward whirling motion when position gain $p=1.35$, while increasing position gain to $p=1.45$ resulted in backward whirling oscillation. Figure 4, Figure 5 and Figure 6 generally illustrate that the RAMBS had tristable solutions at the same spinning speed (i.e., at $\mathsf{\Omega }=\omega +0.1$), where the system exhibited each of their solutions depending on the initial conditions. In addition, Figure 4c, Figure 5c, and Figure 6c confirmed that the system oscillation amplitudes when $p=1.35$ were always smaller than the oscillation amplitudes in the case of $p=1.45$.

The sensitivity of the RAMBS oscillation amplitudes ($a, b$) to linear velocity control gain ($d$) is explored through Figure 7. Figure 7a,b show the system spinning speed response curves at $d=0.01$, and Figure 7c,d illustrate the same response curves at $d=0.015$. Comparing Figure 7a,b with Figure 7c,d shows that the increase in linear velocity gain from $d=0.01$ to $d=0.015$ decreased system lateral vibrations and eliminated different nonlinear phenomena. such as the coexistence of multisolution and sensitivity to initial conditions, as Figure 7c,d show. The decreasing RAMBS oscillation amplitudes as velocity control gain increased can be explained depending on the definition of the system parameters in Equations (11) and (12), where the magnitude of linear damping coefficient $\mu =4d\cos \left(\alpha \right)$ was directly dependent on velocity gain $d$. Therefore, increasing the linear velocity control gain increased the whole system’s linear damping coefficient, which ultimately decreased the sensitivity of the RAMBS to the large excitation of disc eccentricity $f$.

4.2. Sensitivity Analysis of Nonlinear Position-Velocity Controller (${\lambda }_1$ and ${\lambda }_2$)

This section explores the performance of the nonlinear cubic position-velocity controller (i.e., ${\lambda }_1$ and ${\lambda }_2$) in mitigating the RAMBS lateral oscillations when the other control parameters were fixed constant (i.e., $p=1.22$ and $d=0.005$). First, the influence of the cubic position gain (${\lambda }_1$) on system dynamical behaviours is discussed through Figure 8 and Figure 9. Figure 8 shows the system spinning speed response curves when cubic position gain ${\lambda }_1=0.05$ and ${\lambda }_1=0.1$, while Figure 9 shows the same response curve at ${\lambda }_1=-0.1$ and $-0.3$. Figure 8 shows that the increase in nonlinear position gain (${\lambda }_1$) from ${\lambda }_1=0.05$ to ${\lambda }_1=0.1$ bent the RAMBS response curve to the right and simultaneously increased the system oscillation amplitudes. Figure 9 illustrates that the increase in negative values of the cubic position gain from ${\lambda }_1=-0.1$ (Figure 9a,b) to ${\lambda }_1=-0.3$ (Figure 9c,d) bent the system frequency response curve to the left, leading to softening spring behaviours. In addition, Figure 9 shows that the increase in negative values of ${\lambda }_1$ from ${\lambda }_1=-0.1$ to ${\lambda }_1=-0.3$ minimised the system oscillation amplitudes, and the system response curves became simpler, where intervals at which the system could respond with a tristable solution (as in Figure 8 and Figure 9a,b) were eliminated at ${\lambda }_1=-0.3$, as shown in Figure 9c,d.

According to Figure 8 and Figure 9, the bifurcation behaviours of the RAMBS were explored by utilising the cubic position gain as the main bifurcation parameter at two different values of disc spinning speed $\mathsf{\Omega }$ ($\mathsf{\Omega }=\omega +\sigma , \sigma =0$ and $0.05$), as shown in Figure 10. Comparing Figure 10a,b with Figure 10c,d shows that the cubic position gain interval at which the RAMBS could exhibit multiple solutions shifted to the right as the detuning parameter ($\sigma$) increased. In addition, the figure confirms that the RAMBS could exhibit a small oscillation amplitude with a single periodic solution only if ${\lambda }_1<-0.2$, regardless of the magnitude of $\sigma$. According to Figure 8, Figure 9 and Figure 10, to avoid the high oscillation amplitudes and sensitivity to the initial conditions (i.e., avoiding the multistable solution interval), the cubic position control gain should be selected to be negative (i.e., ${\lambda }_1<-0.2$) if the system operates at rotational speed ($\mathsf{\Omega }$) higher than or equal to its natural frequency ($\omega$), as Figure 9a,b and Figure 10 show. On the other hand, if the RAMBS operates at a lower rotational speed ($\mathsf{\Omega }$) than its natural frequency (i.e., $\sigma <0$), the cubic position gain (${\lambda }_1$) should be designed to be positive (i.e., ${\lambda }_1>0$), as Figure 8c,d and Figure 10 show.

The effect of nonlinear velocity gain (${\lambda }_2$) on the vibrational motion of the RAMBS is explored through Figure 11 and Figure 12. Figure 11 illustrates the RAMBS spinning speed response curve at two different positive values of ${\lambda }_2$, while Figure 12 shows the effect of the negative cubic velocity gain on the oscillatory behaviours of the RAMBS. Figure 11 shows that the cubic velocity controller works as a damping element as in the case of the linear velocity controller, where increasing ${\lambda }_2$ from $0.05$ (as in Figure 11a,b) to $0.15$ (as in Figure 11c,d) suppressed the system nonlinearity, and the bistable solution interval disappeared. However, Figure 12 shows that the negative gain (i.e., ${\lambda }_2<0$) of the cubic velocity controller is unacceptable, where the RAMBS may lose its stability to respond with unbounded oscillation amplitudes at a specific range of the disc spinning speed (i.e., lose its stability when $-0.03<\sigma <0.05$). Accordingly, the optimal working condition for the cubic velocity controller is designing its gain to be a positive value (i.e., ${\lambda }_2$ must be greater than zero).

According to Figure 12 at $\sigma =0.0$, the temporal oscillation and the corresponding orbital motions of the RAMBS are illustrated in Figure 13 by numerically solving Equations (11) and (12) using MATLAB ODE45 solver at the zero initial conditions. Comparing Figure 12 and Figure 13 at $\sigma =0.0$ confirms the accurate correspondence between the two figures, where numerical results illustrate the unstable motion of the RAMBS at $\sigma =0.0$.

According to Figure 11 and Figure 12, the cubic velocity control gain (${\lambda }_2$) was employed as the main bifurcation parameter at two different values of the detuning parameter ($\sigma$), as shown in Figure 14. Figure 14a,b show ${\lambda }_2$-response curves when $\sigma =0.0$ (i.e., at perfect tuning $\mathsf{\Omega }=\omega$), while Figure 14c,d illustrate ${\lambda }_2$-response curves at $\sigma =0.05$. Figure 14 shows hat as the cubic velocity feedback gain increased, the RAMBS oscillation amplitudes ($a$ and $b$) dramatically decreased, and system nonlinear characteristics disappeared (i.e., the tristable solution interval was merged into a single solution at the large values of ${\lambda }_2$).

The temporal oscillations and the corresponding orbital motion of the RAMBS are illustrated in Figure 15 and Figure 16 according to Figure 14c,d when ${\lambda }_2=0$ and $0.1$, respectively. Figure 14c,d show that the RAMBS had a tristable solution at ${\lambda }_2=0$, where the system responded to one of them only depending on the initial conditions. Therefore, the system equations of motion (i.e., Equations (11) and (12)) were numerically solved at the three initial conditions $x\left(0\right)=\dot{x}\left(0\right)=y\left(0\right)=\dot{y}\left(0\right)=0.0$, $x\left(0\right)=\dot{x}\left(0\right)=0.3, y\left(0\right)=\dot{y}\left(0\right)=0.0$, and $x\left(0\right)=\dot{x}\left(0\right)=0.0, y\left(0\right)=\dot{y}\left(0\right)=0.3$ as depicted in Figure 15 when cubic velocity control gain ${\lambda }_2=0.0$. Figure 15 confirms the existence of three stable solutions of the RAMBS at ${\lambda }_2=0.0$. In addition, numerical simulations for the RAMBS according to Figure 14c,d at ${\lambda }_2=0.1$ are illustrated in Figure 16 at the same initial conditions that were utilised to obtain Figure 15. Figure 16 shows that the three stable solutions reported in Figure 15 were merged into a single solution.

RAMBS eccentricity response curves (i.e., Figure 17a,b) and corresponding bifurcation diagrams (i.e., Figure 17c,d) were explored when the nonlinear control gains are set to be zero (i.e., ${\lambda }_1={\lambda }_2=0.0$). Figure 17a,b show that the RAMBS responded as a linear system with a single periodic attractor as long as disc eccentricity $f<0.005$. However, the nonlinear behaviours dominated the system response when the rotor eccentricity was greater than $0.05$. Figure 17a,b indicate that the system may have bistable solutions as long as $0.005<f<0.026$. Moreover, the figures show that the system may lose its stability when the disc eccentricity exceeded a critical value (i.e., when $f>0.026$). According to Figure 17a,b, the RAMBS bifurcation diagrams are illustrated in Figure 17c,d to investigate the nature of the system’s unstable motions when disc $f>0.026$. Figure 17c,d were developed by numerically solving Equations (11) and (12) utilising the ODE45 MATLAB solver, where the steady-state Poincaré map was obtained and then plotted by employing the disc eccentricity as the bifurcation parameter. Comparing Figure 17a,b with Figure 17c,d confirmed that the RAMBS could periodically oscillate as long as $f<0.026$, but if the disc eccentricity exceeded the critical value (i.e., if $f>0.026$), the system exhibited aperiodic motion.

A numerical simulation for the RAMBS temporal oscillations according to Figure 17 when $f=0.04$ at zero initial condition is depicted in Figure 18, where Figure 18a,b show the system temporal oscillations in the $X$ and $Y$ directions, Figure 18c illustrates the corresponding orbital motion, Figure 18d,e represent the Poincaré map, and Figure 18f,g are the corresponding frequency spectrum. The Poincaré map and the frequency spectrum showed that the system performed a quasiperiodic oscillation at $f=0.04$ that accurately agreed with the obtained results in Figure 17.

The influence of cubic position control gain (${\lambda }_1$) on the RAMBS eccentricity response curves is explored in Figure 19. Figure 19a,b show the system response curves at ${\lambda }_1=-0.1$, and Figure 19c illustrates the system response curves at ${\lambda }_1=-0.2$ when $\sigma =0.0$. Comparing Figure 19 with Figure 17a,b shows that the unstable system motion at the large disc eccentricity (i.e., when $f>0.026$) that is shown in Figure 17a,b was stabilised, as shown in Figure 19, when the cubic position gain was activated (i.e., when ${\lambda }_1=-0.1$ or $-0.2$). Moreover, comparing Figure 19a,b with Figure 19c,d confirmed that the increase in cubic position gain from ${\lambda }_1=-0.1$ to ${\lambda }_1=-0.2$ eliminated the system nonlinearities, where the bistable solutions interval in Figure 19a,b disappeared, as shown in Figure 19c,d, and the system responded as a linear system with a single periodic attractor along the eccentricity axis. Accordingly, the cubic position controller with negative control gain (i.e., ${\lambda }_1<0$) could stabilise the unstable motion of the RAMBS and force the system to respond as a linear one.

Numerical validations to Figure 19 when $f=0.04$ at ${\lambda }_1=-0.1$ and $-0.2$ are presented in Figure 20 and Figure 21, respectively. Figure 20 simulates the RAMBS oscillation amplitudes and the corresponding orbital motion according to Figure 19a,b when $f=0.04$ at two different initial conditions, while Figure 21 simulates the RAMBS oscillation amplitudes and the corresponding orbital motion according to Figure 19c,d when $f=0.04$ at the same initial conditions. Comparing Figure 18, Figure 20, and Figure 21 shows that the unstable motion shown in Figure 18 was stabilised by applying the cubic position feedback as illustrated in Figure 20 and Figure 21. In addition, the increase in cubic position from ${\lambda }_1=-0.1$ as in Figure 20 to ${\lambda }_1=-0.2$ as in Figure 21 eliminated the sensitivity of the system to the initial conditions. The effect of the cubic velocity control gain (${\lambda }_2$) on the RAMBS eccentricity response curves is investigated through Figure 22. where Figure 22a,b show the eccentricity response curves at ${\lambda }_2=0.05$, and Figure 22c,d depict the system response curves at ${\lambda }_2=0.15$ when $\sigma =0.0$. Comparing Figure 22 and Figure 17a,b confirmed that the unstable system motion at the large disc eccentricity (i.e., when $f>0.026$) in Figure 17a,b was stabilised, as shown in Figure 22, either when ${\lambda }_2=0.05$ or $0.15$. In addition, comparing Figure 22a,b and Figure 22c,d shows that the increase in cubic velocity gain from ${\lambda }_2=0.05$ to ${\lambda }_2=0.15$ eliminated system nonlinearities at a wide range of the disc eccentricity, where the bistable solution interval in Figure 22a,b was narrowed, as shown in Figure 22c,d. Therefore, the cubic velocity controller with positive control gain (i.e., ${\lambda }_2>0$) could stabilise the unstable motion of the RAMBS and force the system to respond as a linear one.

5. Conclusions

A cubic position-velocity feedback controller was proposed to enhance the control performance of a rotor-active magnetic-bearings system. The suggested nonlinear controller was along with a conventional linear position-velocity controller into an 8-pole RAMBS. According to the introduced control law, the system dynamical model was established and then analysed utilising perturbation techniques. Slow-flow autonomous differential equations that govern system vibration amplitudes and the modified phases were derived. The influence of both the linear and nonlinear control gains on the system dynamics were explored through different response curves and bifurcation diagrams. The acquired analytical solutions and corresponding numerical simulations confirmed that the nonlinear controller could enhance the dynamical characteristics of the studied system by adding many important features to the 8-pole system, summarised as follows:

• Optimal linear position gain $p$ should be as small as possible; however, it should be greater than $\frac{1}{\cos \left(\alpha \right)}$ (i.e., gain $p>\frac{1}{\cos \left(\alpha \right)}$) to guarantee system stability by making system natural frequency $\omega =\sqrt{8\left(p\cos \left(\alpha \right)-1\right)}$ always have a positive value.
• Integrating the cubic position controller (${\lambda }_1$) into the linear controller makes the control algorithm more flexible to changing the system dynamical behaviours from the hardening spring characteristic to the softening spring characteristic (or vice versa) by designing the suitable values of ${\lambda }_1$ without any constraints to avoid the resonance conditions.
• Selecting the cubic position gain (${\lambda }_1$) with large negative values can simplify the system dynamical behaviours and mitigate system oscillations, even at resonance conditions.
• The good design of the cubic position gain (i.e., ${\lambda }_1<0$) can stabilise the unstable motion and eliminate the nonlinear effects of the system at large disc eccentricities.
• Integrating the cubic velocity controller (${\lambda }_2$) to the linear controller added a nonlinear damping term to the controlled system that improved system stability or destabilised its motion, depending on the control gain sign.
• The optimal design of the cubic velocity gain (i.e., ${\lambda }_2>0$) could stabilise the unstable motion and eliminate the nonlinear effects of the system at large disc eccentricities