On the oscillatory behaviours and rub-impact forces of a horizontally supported asymmetric rotor system under position-velocity feedback controller

Abstract

This work is meant to investigate and control the nonlinear dynamical behaviours of a nonlinear asymmetric rotating shaft system. The studied system is modeled as a nonlinear dynamical system having both quadratic and cubic nonlinearities and excited at the primary resonance. A linear proportional-derivative controller is introduced to eliminate the nonlinear behaviours and to suppress the lateral vibrations of the rotating shaft. The influences of different control parameters on the oscillatory behaviours of the considered system are explored. The main obtained analytical results showed that the uncontrolled shaft may respond with forward and backward whirling motion. However, the proposed controller illustrated its feasibility in eliminating the nonlinear behaviours of the studied system and forcing it to behave like the linear systems. Moreover, the obtained results confirmed that the proportional control gain plays a dominant role in destabilizing the studied rotor system. Finally, the safe operation of the studied system with avoiding the rub-impact force occurrence is discussed.

Keywords:
Asymmetric rotor; forward and backward whirling motion; stability; Bi-stable solution; Tri-stable solution; Proportional-derivative controller; rub-impact forces

Graphical Abstract

1. Introduction

The nonlinear lateral vibrations are a common phenomenon in rotating machinery in general. Many sources may cause these vibrations such as the imbalances of the rotating shafts, the coupling misalignments in multi-rotor systems, the wear between the rotating shafts and the ball bearings, the propagation of the transverse cracks, and the shafts asymmetry, …etc. Lateral vibration is an undesired phenomenon that passively influences the efficiency, durability, lifetime, and reliability of the rotating machines. The existence of asymmetry in the stiffness coefficients of the rotors is a common feature of a special category of rotating machines as in the case of the two-pole rotor generator. The two-pole generator has two slots in two opposite directions in order to accommodate the electric coils. Accordingly, the asymmetric design of these machines causes the rotating shaft stiffness coefficients to be asymmetric. The asymmetry of the rotation shaft systems is resulting in completely different oscillatory behaviours than the symmetric shafts. Therefore, investigating the vibratory characteristics of the asymmetric shafts is a research subject of many researchers, where Ardayfio and Frohrib [¹] studied the nonlinear vibrations of an asymmetric rotating shaft system that is suspended utilizing two linear ball bearings. They explored the effect of the bearing stiffness coefficients on the stability of the proposed model. The authors concluded that the stiffness coefficients of the bearings have a dominant role in both stabilizing and destabilizing the studied rotor system. Iwatsubo et al. [²] explored the forced and free dynamical behaviours of an asymmetric rotor system supported using asymmetric bearings. The whole system mathematical model was derived and then investigated applying perturbation analyses. According to the obtained results, the authors reported that the asymmetry of both the rotating shaft and the bearings can destabilizing the whole system, independently. The dynamical behaviours of a two-pole asymmetric rotor generator are investigated in Ref. [³], where the main obtained result was that the two-pole rotor system may vibrate with a strong lateral oscillation when the shaft angular speed is two times the system critical speed. The nonlinear lateral-torsional oscillations of an asymmetric rotor system that is modeled as Timoshenko continuous beam is analyzed in Ref. [⁴]. The authors concluded that the system can exhibit high vibration amplitude when the shaft spinning speed becomes two times the system's critical speed. Ishida and Liu [⁵] investigated the oscillatory behaviours of an asymmetric rotor system. The authors applied a discontinuous spring to stabilize the unstable motion of the studied system. The dynamical behaviours of an asymmetrical rotating shaft that having stretching nonlinearity is investigated [⁶]. The lateral vibration of cracked asymmetric rotor systems is analyzed [⁷-¹⁰]. The authors reported that the increase of the crack size may be the main reason for stabilizing the system oscillation. Xiang et al. [¹¹] discussed numerically the dynamic behaviours of an asymmetric shaft having a transverse crack and subjected to rub and impact forces. The obtained results showed that the studied model can perform period-n and quasiperiodic motions. Przybylowicz et al. [¹²] analyzed numerically the dynamical characteristics of an asymmetric rotating shaft that is subjected to steady excitations. They illustrated the possibility of excluding the system's chaotic motion via designing stiffness and damping coefficients appropriately. Yu et al. [¹³] explored the global dynamical behaviour of a flexible asymmetric rotating shaft system. The obtained results showed that the studied model could perform a chaotic motion in the sense of Smale-Horseshoes. More detailed investigations for the dynamical behaviours of an asymmetric rotor system suspended by elastic support could be found in Boru's nice work [¹⁴]. Saeed [¹⁵, ¹⁶] introduced analytical and numerical investigations for a nonlinear asymmetric rotating shaft when supported either vertically or horizontally. Linear and nonlinear asymmetric stiffness coefficients are included in the studied model. He reported that the asymmetric system has a complicated dynamical behaviours than the symmetric one. Moreover, it could perform single-stable, bi-stable, tri-stable, and quadri-stable solutions depending on the initial conditions and the shaft angular velocity.

Besides the above-mentioned studies, many research articles have introduced different control strategies to control the oscillatory motions of rotating machines. The principle for the vibrations reduction in the rotating machines is based on one of two main strategies. The first strategy utilizes passive elements such as the vibrations absorbers and energy sinks in different configurations to get rid of the vibrational energy from the targeted systems [¹⁷, ¹⁸], while the second strategy applies active control signals via electromagnetic actuators to suppress the oscillatory motions of the rotating shafts [¹⁹-²³]. The lateral oscillations of a Jeffcott rotor system when supported vertically is discussed in [²⁴]. The authors integrated a linear dynamic absorber to the rotating shaft via four electromagnetic poles to control the system's nonlinear vibrations. Saeed and Kamel [²⁵] proposed a linear position-velocity feedback controller to eliminate the nonlinear oscillations of a rotor system supported horizontally. The authors applied the control currents via electromagnetic poles that acted as actuators. Recently, Ghasabi et al. [²⁶], and Saeed et al. [²⁷] applied two different active control algorithms to mitigate the nonlinear vibrations of an asymmetric rotating shaft. Ghasabi et al. [²⁶] employed a proportional-derivative controller with time-delay to mitigate an asymmetric shaft lateral vibration. Saeed et al. [²⁷] proposed a nonlinear position-velocity controller to suppress the nonlinear dynamical behaviours of a vertically supported asymmetric nonlinear rotating shaft. The authors have integrated the suggested control technique to the rotor utilizing an electromagnetic actuator.

As the asymmetrical rotating machinery did not get adequate research attention regarding vibration control, the motivation of this article is to design a control algorithm in order to mitigate the nonlinear oscillations of an asymmetric horizontally supported nonlinear rotor system. The studied rotating shaft is governed by two coupled nonlinear differential equations having both quadratic and cubic nonlinearities. The shaft weight, the stiffness coefficients, and the eccentricity orientation angle are included in the studied mathematical model. A linear position-velocity controller is suggested to mitigate the considered system lateral vibrations. The system mathematical model is derived and then analyzed by applying the multiple scales perturbation method. Based on the obtained analytical and numerical simulations, it is found that the system before control may suffer from bi-stable and tri-stable periodic oscillations that are result in forward or backward whirling motion. However, the applied control strategy showed its capability in forcing the studied nonlinear model to behave like a linear system with a single periodic solution. Moreover, the numerical results approved that the proportional control gain can play an important role in stabilizing and destabilizing the asymmetric rotor system.

2. Asymmetric system mathematical model

The oscillatory motion of a nonlinear rotor system supported horizontally as shown in Figs.1a and 1b is governed by the following nonlinear differential equations [²⁸, ²⁹]:

$m\ddot{u}+c_1\dot{u}+F_u=m{e}_d{\omega }^2\cos (\omega \tau +\gamma )$ (1.a)

$m\ddot{v}+c_2\dot{v}+F_v=m{e}_d{\omega }^2\sin (\omega \tau +\gamma )-W$ (1.b)

where $m$ is the rotating disk mass, $c_1$ and $c_2$ are the linear damping coefficients in $U$ and $V$directions, $F_u$ and $F_v$ are the shaft restoring forces in $U$ and $V$directions, $e_d$ represents the disk eccentricity that is the distance between the disk geometric center $G\left(u,v\right)$ and its centroid $C\left(u_e,v_e\right)$, $\omega$ is the rotating shaft spinning speed, $W=mg$ is the disk weight, and $\gamma$ denotes the orientation angle between $\overrightarrow{GC}$ and $\overrightarrow{o{z}_1}$ axes. It is assumed that the considered system performs planner vibrations in the inertial $U–V$ plane, while $z_1-z_2$ is a rotational coordinate with spinning speed equal $\omega$. For a small radial displacement $r$ of the shaft system as in Fig. 1b, the exerted restoring forces in $U$ and $V$ direction can be expressed as [²⁸-²⁹]:

$F_u=k_{1}r\cos \psi +k_2{r}^3\cos \psi =k_{1}u+k_2(u^2+v^2)u$ (2.a)

$F_v=k_{1}r\sin \psi +k_2{r}^3\sin \psi =k_{1}v+k_2(u^2+v^2)v$ (2.b)

where $k_1$ and $k_2$ are the linear and nonlinear stiffness coefficients of the elastic shaft. By inserting Eqs. (2.a) and (2.b) into Eqs. (1.a) and (1.b), we have

$m\ddot{u}+c_1\dot{u}+k_{1}u+k_2(u^2+v^2)u=m{e}_d{\omega }^2\cos (\omega \tau +\gamma )$ (3.a)

$m\ddot{v}+c_2\dot{v}+k_{1}v+k_2(u^2+v^2)v=m{e}_d{\omega }^2\sin (\omega \tau +\gamma )-W$ (3.b)

Due to the disk weight$W$, the rotor system may have small static deflection$v_s$ in $V$ direction as illustrated in Fig. 1a. Accordingly, at static equilibrium (i.e.$\ddot{u}=\dot{u}=u=\ddot{v}=\dot{v}=\omega =0$), from equation (3.b) we have

$k_1{v}_s+k_2{v}_s^3=-W$ (4)

By using the new coordinate system $X,Y$ such that $U=X$ and $V=v_s+Y$, we have $\ddot{u}=\ddot{x},\dot{u}=\dot{x},u=x,\ddot{v}=\ddot{y},\dot{v}=\dot{y},v=v_s+y$. Substituting this translation into Eqs. (3.a) and (3.b), with considering Eq. (4), we can obtain the following equations of motion:

$m\ddot{x}+c_1\dot{x}+(k_1+k_2{v}_s^2)x+2k_2{v}_{s}xy+k_2(x^2+y^2)x=m{e}_d{\omega }^2\cos (\omega \tau +\gamma )$ (5.a)

$m\ddot{y}+c_2\dot{y}+(k_1+3k_2{v}_s^2)y+k_2{v}_s(x^2+3y^2)+k_2(x^2+y^2)y=m{e}_d{\omega }^2\sin (\omega \tau +\gamma )$ (5.b)

For asymmetric rotating shafts (i.e. two-pole rotor generators ), the restoring force in $z_1$ direction may be larger than that in $z_2$ direction due to the coil slots as Fig. 1c implies. Accordingly, Eqs. (5.a) and (5.b) should be modified to:

$m\ddot{x}+c_1\dot{x}+(k_1+k_2{v}_s^2)x+2k_2{v}_{s}xy+k_2(x^2+y^2)x+\Delta F_x=m{e}_d{\omega }^2\cos (\omega \tau +\gamma )$ (6.a)

$m\ddot{y}+c_2\dot{y}+(k_1+3k_2{v}_s^2)y+k_2{v}_s(x^2+3y^2)+k_2(x^2+y^2)y+\Delta F_y=m{e}_d{\omega }^2\sin (\omega \tau +\gamma )$ (6.b)

where $\Delta F_x$ and $\Delta F_y$ are modified terms for the shaft restoring forces in $X$ and $Y$ direction because of the shaft asymmetry

Fig. 1
Asymmetric rotor system supported horizontally: (a) Jeffcott-rotor system, (b) symmetric rotating disk, (c) asymmetric rotating disk.

As Fig. 1c implies, assume that the shaft restoring force in $z_1-$ direction is greater than the restoring force in $z_2$ direction by the value $\Delta F$, where $\Delta F$ is assumed as follows [¹⁵,¹⁶, ³⁰]:

$\Delta F=\Delta k_1{z}_1+\Delta k_2{z}_1^3=\Delta k_1\left(x\cos (\omega \tau )+y\sin (\omega \tau )\right)+\Delta k_2{\left(x\cos (\omega \tau )+y\sin (\omega \tau )\right)}^3$ (7)

where $\Delta k_1$ and $\Delta k_2$ are the extra linear and nonlinear shaft stiffness coefficients in $z_1$ direction. So, the forces $\Delta F_x$ and $\Delta F_y$ in Eqs. (6.a) and (6.b) can be written as:

$\begin{array}{c}\Delta F_x=\Delta F\cos (\omega t)=\frac{1}{2}\Delta k_{1}x+\frac{3}{8}\Delta k_2{x}^3+\frac{3}{8}\Delta k_{2}x{y}^2+\frac{1}{2}\Delta k_{1}x\cos (2\omega \tau )+\frac{1}{2}\Delta k_{1}y\sin (2\omega \tau )\\ +\frac{1}{2}\Delta k_2{x}^3\cos (2\omega \tau )+\frac{3}{4}\Delta k_2{x}^{2}y\sin (2\omega \tau )+\frac{1}{4}\Delta k_2{y}^3\sin (2\omega \tau )+\frac{1}{8}\Delta k_2{x}^3\cos (4\omega \tau )\\ -\frac{1}{8}\Delta k_2{y}^3\sin (4\omega \tau )+\frac{3}{8}\Delta k_2{x}^{2}y\sin (4\omega \tau )-\frac{3}{8}\Delta k_{2}x{y}^2\cos (4\omega \tau )\end{array}$ (8.a)

$\begin{array}{c}\Delta F_y=\Delta F\sin (\omega t)=\frac{1}{2}\Delta k_{1}y+\frac{3}{8}\Delta k_2{y}^3+\frac{3}{8}\Delta k_2{x}^{2}y-\frac{1}{2}\Delta k_{1}y\cos (2\omega \tau )+\frac{1}{2}\Delta k_{1}x\sin (2\omega \tau )\\ -\frac{1}{2}\Delta k_2{y}^3\cos (2\omega \tau )+\frac{3}{4}\Delta k_{2}x{y}^2\sin (2\omega \tau )+\frac{1}{4}\Delta k_2{x}^3\sin (2\omega \tau )+\frac{1}{8}\Delta k_2{y}^3\cos (4\omega \tau )\\ +\frac{1}{8}\Delta k_2{x}^3\sin (4\omega \tau )-\frac{3}{8}\Delta k_{2}x{y}^2\sin (4\omega \tau )-\frac{3}{8}\Delta k_2{x}^{2}y\cos (4\omega \tau )\end{array}$ (8.b)

Substituting Eqs. (8.a) and (8.b) into Eqs. (6.a) and (6.b), we can obtain the following nonlinear dynamical system.

$\begin{array}{l}m\ddot{x}+c_1\dot{x}+(k_1+\frac{1}{2}\Delta k_1+k_2{v}_s^2)x+2k_2{v}_{s}xy+(k_2+\frac{3}{8}\Delta k_2)(x^2+y^2)x=m{e}_d{\omega }^2\cos (\omega \tau +\gamma )\\ -\frac{1}{2}\Delta k_{1}x\cos (2\omega \tau )-\frac{1}{2}\Delta k_{1}y\sin (2\omega \tau )-\frac{1}{2}\Delta k_2{x}^3\cos (2\omega \tau )-\frac{3}{4}\Delta k_2{x}^{2}y\sin (2\omega \tau )\\ -\frac{1}{4}\Delta k_2{y}^3\sin (2\omega \tau )-\frac{1}{8}\Delta k_2{x}^3\cos (4\omega \tau )+\frac{1}{8}\Delta k_2{y}^3\sin (4\omega \tau )-\frac{3}{8}\Delta k_2{x}^{2}y\sin (4\omega \tau )\\ +\frac{3}{8}\Delta k_{2}x{y}^2\cos (4\omega \tau )+F_{XC}\end{array}$ (9.a)

$\begin{array}{l}m\ddot{y}+c_2\dot{y}+(k_1+\frac{1}{2}\Delta k_1+3k_2{v}_s^2)y+k_2{v}_s(x^2+3y^2)+(k_2+\frac{3}{8}\Delta k_2)(x^2+y^2)y\\ =m{e}_d{\omega }^2\sin (\omega \tau +\gamma )+\frac{1}{2}\Delta k_{1}y\cos (2\omega \tau )-\frac{1}{2}\Delta k_{1}x\sin (2\omega \tau )+\frac{1}{2}\Delta k_2{y}^3\cos (2\omega \tau )\\ -\frac{3}{4}\Delta k_{2}x{y}^2\sin (2\omega \tau )-\frac{1}{4}\Delta k_2{x}^3\sin (2\omega \tau )-\frac{1}{8}\Delta k_2{y}^3\cos (4\omega \tau )-\frac{1}{8}\Delta k_2{x}^3\sin (4\omega \tau )\\ +\frac{3}{8}\Delta k_{2}x{y}^2\sin (4\omega \tau )+\frac{3}{8}\Delta k_2{x}^{2}y\cos (4\omega \tau )+F_{YC}\end{array}$ (9.b)

where $F_{XC}$ and $F_{YC}$ are additional control forces to suppress the nonlinear vibrations of the asymmetric rotor shaft governed by Eqs. (9.a) and (9.b). To apply control force on the rotating shaft without any physical contact, the suggested control forces $F_{XC}$ and $F_{YC}$ have been suggested to be applied via four electromagnetic poles as shown in Fig. 2. Poles 1 and 3 generate the electromagnetic force $F_{XC}$ to control the lateral vibration in $X$diretcion, while poles 2 and 4 generate the electromagnetic force $F_{YC}$ to control the lateral vibration in $Y$diretcion. According to the electromagnetic theory, the total force generated in both $X$ and $Y$ directions can be expressed as follows [³¹]:

$F_{XC}=F(x,i_x)=B_0\left[\frac{{\left(I_0+i_x\right)}^2}{{\left(g_0+x\right)}^2}-\frac{{\left(I_0-i_x\right)}^2}{{\left(g_0-x\right)}^2}\right]$ (10.a)

$F_{YC}=F(y,i_y)=B_0\left[\frac{{\left(I_0+i_y\right)}^2}{{\left(g_0+y\right)}^2}-\frac{{\left(I_0-i_y\right)}^2}{{\left(g_0-y\right)}^2}\right]$ (10.b)

where $B_0$ is the magnetic force constant,$\left(g_0+x\right)$ is the air-gap size between the pole-leg and the shaft in$X$direction, $\left(g_0+y\right)$ is the air-gap size between the pole-leg and the shaft in$Y$direction,$i_x$ and $i_y$ are the control currents at $X$ and $Y$ directions, respectively, and $I_0$ is the bias current.

Fig. 2
Asymmetrical horizontally supported rotor system with 4-pole active magnetic bearings

Within this article, the linear proportional-derivative control law is proposed to generate control currents that are proportional to both the vibration displacement and corresponding velocity in $X$ and $Y$ directions as follows:

$i_x={\gamma }_{1}x+{\gamma }_2\dot{x},i_y={\gamma }_{1}y+{\gamma }_2\dot{y}$ (11)

where ${\gamma }_1$ is a linear proportional control gain, ${\gamma }_2$ is a linear derivative control gain. Substituting Eq. (11) into Eqs. (10.a) and (10.b) with expanding the resulting equations using Taylor series up to the third-order approximation, and then inserting the obtained equations into Eqs. (9.a) and (9.b) with introducing the dimensionless parameters $t=\sqrt{k_1/m}\tau ,q_1=\frac{x}{g_0},q_2=\frac{y}{g_0}$, we can derive the dimensionless equations of motions that govern the controlled system dynamics as follows:

$i_x={\gamma }_{1}x+{\gamma }_2\dot{x},i_y={\gamma }_{1}y+{\gamma }_2\dot{y}$ (12.a)

$\begin{array}{l}{\ddot{q}}_2+{\mu }_2{\dot{q}}_2+{\omega }_2^2{q}_2+\eta (q_1^2+3q_2^2)+\beta (q_1^2+q_2^2)q_2=f{\Omega }^2\sin (\Omega t+\gamma )+\frac{1}{2}{\delta }_1{q}_2\cos (2\Omega t)\\ -\frac{1}{2}{\delta }_1{q}_1\sin (2\Omega t)+\frac{1}{2}{\delta }_2{q}_2^3\cos (2\Omega t)-\frac{1}{4}{\delta }_2(3q_1{q}_2^2+q_1^3)\sin (2\Omega t)+\frac{1}{8}{\delta }_2(3q_1^2{q}_2-q_2^3)\cos (4\Omega t)\\ +\frac{1}{8}{\delta }_2(3q_1{q}_2^2-q_1^3)\sin (4\Omega t)-({\alpha }_1{q}_2+{\alpha }_2{\dot{q}}_2+{\alpha }_3{q}_2^3+{\alpha }_4{q}_2^2{\dot{q}}_2+{\alpha }_5{q}_2{\dot{q}}_2^2)\end{array}$ (12.b)

where${\mu }_1=\frac{c_1}{\sqrt{k_{1}m}},{\mu }_2=\frac{c_2}{\sqrt{k_{1}m}},{\omega }_n^2=\frac{k_1}{m}=\frac{\delta I_0^2}{m{g}_0},\mathrm{\Omega }=\frac{\omega }{{\omega }_n},f=\frac{e_d}{g_0},\lambda =\frac{v_s}{g_0},\rho =\frac{g_0^2{k}_2}{k_1},\eta =\lambda \rho ,{\delta }_1=\frac{{\mathrm{\Delta }k}_1}{k_1},{\delta }_2=\frac{{g_0^2\mathrm{\Delta }k}_2}{k_1},{\omega }_1=\sqrt{1+{\lambda }^2\rho +\frac{1}{2}{\delta }_1},{\omega }_2=\sqrt{1+3{\lambda }^2\rho +\frac{1}{2}{\delta }_1},\beta =\rho +\frac{3}{8}{\delta }_2,p=\frac{g_0}{I_0}{\gamma }_1,d=\frac{g_0{\omega }_n}{I_0}{\gamma }_2,{\alpha }_1=p-1,{\alpha }_2=d,{\alpha }_3=-2+3p-p^2,{\alpha }_4=3d-2pd,{\alpha }_5=-d^2.$

3. Analytical investigations and slow flow modulating equations

To explore the oscillatory behaviour of the derived nonlinear dynamical, the multiple time scales perturbation method [³²] is utilized to obtain an approximate analytical solution. Accordingly, we can obtain a second-order approximate solution to Eqs. (12.a) and (12.b) by seeking the solution in the following form:

$q_1(t)=q_{11}(T_0,T_1,T_2)+\epsilon q_{12}(T_0,T_1,T_2)+{\epsilon }^2{q}_{13}(T_0,T_1,T_2)+O({\epsilon }^3)$ (13.a)

$q_2(t)=q_{21}(T_0,T_1,T_2)+\epsilon q_{22}(T_0,T_1,T_2)+{\epsilon }^2{q}_{23}(T_0,T_1,T_2)+O({\epsilon }^3)$ (13.b)

where$\epsilon$is the perturbation parameter that used as a book-keeping only,$T_0=t,T_1=\epsilon t$ and$T_2={\epsilon }^{2}t$ are the fast and slow time scales to capture the evolving of both the fast and slow motions of the considered system. In terms of $T_0,T_1,$ and $T_2$ the time derivatives can be expressed as follows:

$\frac{d}{dt}=D_0+\epsilon D_1+{\epsilon }^2{D}_2,\frac{d^2}{d{t}^2}=D_0^2+2\epsilon D_0{D}_1+{\epsilon }^2(D_1^2+2D_2{D}_0),D_j=\frac{\partial }{\partial T_j},j=\mathrm{0,1,2}$ (14)

A new scaling for the system parameters depending on their values is considered such that:

${\mu }_1={\epsilon }^2{\widehat{\mu }}_1,{\mu }_2={\epsilon }^2{\widehat{\mu }}_2,\beta ={\epsilon }^2\widehat{\beta },\eta =\epsilon \widehat{\eta },f=\epsilon \widehat{f},{\alpha }_n={\epsilon }^2{\widehat{\alpha }}_n(n=\mathrm{1,2,...,5})$ (15)

Inserting Eqs. (13.a) to (15) into Eqs. (18.a) and (18.b) with equating the coefficients with like power of$\epsilon$, we get the following set of linear partial differential equations.

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

$(D_0^2+{\omega }_1^2)q_{11}=0$ (16.a)

$(D_0^2+{\omega }_2^2)q_{21}=0$ (16.b)

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

$(D_0^2+{\omega }_1^2)q_{12}=-2D_0{D}_1{q}_{11}-2\widehat{\eta }{q}_{21}{q}_{11}+\widehat{f}{\Omega }^2\cos (\Omega T_0+\gamma )$ (17.a)

$(D_0^2+{\omega }_2^2)q_{22}=-2D_0{D}_1{q}_{21}-\widehat{\eta }{q}_{11}^2-3\widehat{\eta }{q}_{21}^2+\widehat{f}{\Omega }^2\sin (\Omega T_0+\gamma )$ (17.b)

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

$\begin{array}{c}(D_0^2+{\omega }_1^2)q_{13}=-2D_0{D}_1{q}_{12}-(D_1^2+2D_2{D}_0)q_{11}-{\widehat{\mu }}_1{D}_0{q}_{11}-2\widehat{\eta }{q}_{21}{q}_{12}-2\widehat{\eta }{q}_{22}{q}_{11}-\widehat{\beta }{q}_{21}^2{q}_{11}-\widehat{\beta }{q}_{11}^3\\ -\frac{1}{2}{\widehat{\delta }}_1{q}_{11}\cos (2\Omega T_0)-\frac{1}{2}{\widehat{\delta }}_1{q}_{21}\sin (2\Omega T_0)-\frac{1}{2}{\widehat{\delta }}_2{q}_{11}^3\cos (2\Omega T_0)-\frac{1}{4}{\widehat{\delta }}_2(3q_{11}^2{q}_{21}\\ +q_{21}^3)\sin (2\Omega T_0)+\frac{1}{8}{\widehat{\delta }}_2(3q_{11}{q}_{21}^2-q_{11}^3)\cos (4\Omega T_0)+\frac{1}{8}{\widehat{\delta }}_2(q_{21}^3-3q_{11}^2{q}_{21})\sin (4\Omega T_0)\\ -{\widehat{\alpha }}_1{q}_{11}-{\widehat{\alpha }}_2{D}_0{q}_{11}-{\widehat{\alpha }}_3{q}_{11}^3-{\widehat{\alpha }}_4{q}_{11}^2{D}_0{q}_{11}-{\widehat{\alpha }}_5{q}_{11}{(}^{D_0}\end{array}$ (18.a)

$\begin{array}{c}(D_0^2+{\omega }_2^2)q_{23}=-2D_0{D}_1{q}_{22}-(D_1^2+2D_2{D}_0)q_{21}-{\widehat{\mu }}_2{D}_0{q}_{21}-2\widehat{\eta }{q}_{11}{q}_{12}-6\widehat{\eta }{q}_{21}{q}_{22}-\widehat{\beta }{q}_{21}{q}_{11}^2-\widehat{\beta }{q}_{21}^3\\ -\frac{1}{2}{\widehat{\delta }}_1{q}_{11}\sin (2\Omega T_0)+\frac{1}{2}{\widehat{\delta }}_1{q}_{21}\cos (2\Omega T_0)+\frac{1}{2}{\widehat{\delta }}_2{q}_{21}^3\cos (2\Omega T_0)-\frac{1}{4}{\widehat{\delta }}_2(3q_{11}{q}_{21}^2\\ +q_{11}^3)\sin (2\Omega T_0)+\frac{1}{8}{\widehat{\delta }}_2(3q_{11}^2{q}_{21}-q_{21}^3)\cos (4\Omega T_0)+\frac{1}{8}{\widehat{\delta }}_2(3q_{11}{q}_{21}^2-q_{11}^3)\sin (4\Omega T_0)\\ -{\widehat{\alpha }}_1{q}_{21}-{\widehat{\alpha }}_2{D}_0{q}_{21}-{\widehat{\alpha }}_3{q}_{21}^3-{\widehat{\alpha }}_4{q}_{21}^2{D}_0{q}_{21}-{\widehat{\alpha }}_5{q}_{21}{(}^{D_0}\end{array}$ (18.b)

The solution of Eqs. (16.a) and (16.b), can be written as follows:

$q_{11}(T_0,T_1,T_2)=A_1(T_1,T_2)e^{i{\omega }_1{T}_0}+{\overline{A}}_1(T_1,T_2)e^{-i{\omega }_1{T}_0}$ (19.a)

$q_{21}(T_0,T_1,T_2)=A_2(T_1,T_2)e^{i{\omega }_2{T}_0}+{\overline{A}}_2(T_1,T_2)e^{-i{\omega }_2{T}_0}$ (19.b)

where $A_1\left(T_1,T_2\right)$ and$A_1\left(T_1,T_2\right)$ are unknown coefficients of$T_1$ and$T_2$ up to this stage of the analysis and they will be determined later. Inserting Eqs. (19.a) and (19.b) into Eqs. (17.a) and (17.b), we have

$(D_0^2+{\omega }_1^2)q_{12}=-2i{\omega }_1{D}_1{A}_1{e}^{i{\omega }_1{T}_0}-2\widehat{\eta }{A}_1{A}_2{e}^{i({\omega }_1+{\omega }_2)T_0}-2\widehat{\eta }{A}_1{\overline{A}}_2{e}^{i({\omega }_1-{\omega }_2)T_0}+\frac{\widehat{f}{\Omega }^2}{2}{e}^{i(\Omega T_0+\gamma )}+cc$ (20.a)

$(D_0^2+{\omega }_2^2)q_{22}=-2i{\omega }_2{D}_1{A}_2{e}^{i{\omega }_2{T}_0}-\widehat{\eta }(A_1^2{e}^{2i{\omega }_1{T}_0}+A_1{\overline{A}}_1)-3\widehat{\eta }(A_2^2{e}^{2i{\omega }_2{T}_0}+A_2{\overline{A}}_2)-\frac{i\widehat{f}{\Omega }^2}{2}{e}^{i(\Omega T_0+\gamma )}+cc$ (20.b)

where $cc$denotes the complex conjugate of the preceding terms. The simultaneous resonance ($\Omega ={\omega }_1={\omega }_2$) is considered within this work. Accordingly, the closeness of $\Omega$ to ${\omega }_1$ and ${\omega }_2$ can be expressed as follows:

$\Omega ={\omega }_1+\sigma ={\omega }_1+\epsilon \widehat{\sigma },{\omega }_2={\omega }_1+{\sigma }_1={\omega }_1+\epsilon {\widehat{\sigma }}_1$ (21)

where $\sigma$ denotes the closeness of $\Omega$ to${\omega }_1$, and${\sigma }_1$ represents the difference between${\omega }_1$and ${\omega }_2$(i.e. ${\sigma }_1={\omega }_2-{\omega }_1=\sqrt{1+3{\lambda }^2\rho +\frac{1}{2}{\delta }_1}-\sqrt{1+{\lambda }^2\rho +\frac{1}{2}{\delta }_1}$. Substituting Eq. (21) into Eqs. (20. a) and (20.b), yield

$(D_0^2+{\omega }_1^2)q_{12}=2\left(-2i{\omega }_1{D}_1{A}_1+\frac{\widehat{f}{\Omega }^2}{2}{e}^{i(\widehat{\sigma }{T}_1+\gamma )}\right)e^{i{\omega }_1{T}_0}-2\widehat{\eta }{A}_1{A}_2{e}^{i({\omega }_1+{\omega }_2)T_0}-2\widehat{\eta }{A}_1{\overline{A}}_2{e}^{i({\omega }_1-{\omega }_2)T_0}+cc$ (22.a)

$(D_0^2+{\omega }_2^2)q_{22}=\left(-2i{\omega }_2{D}_1{A}_2+\frac{\widehat{f}{\Omega }^2}{2i}{e}^{i((\widehat{\sigma }-{\widehat{\sigma }}_1)T_1+\gamma )}\right)e^{i{\omega }_2{T}_0}-\widehat{\eta }(A_1^2{e}^{2i{\omega }_1{T}_0}+A_1{\overline{A}}_1)-3\widehat{\eta }(A_2^2{e}^{2i{\omega }_2{T}_0}+A_2{\overline{A}}_2)+cc$ (22.b)

The solvability conditions for Eqs. (22.a) and (22.b) are:

$2i{\omega }_1{D}_1{A}_1=\frac{\widehat{f}{\Omega }^2}{2}{e}^{i(\widehat{\sigma }{T}_1+\gamma )}and2i{\omega }_2{D}_1{A}_2=-\frac{i\widehat{f}{\Omega }^2}{2}{e}^{i((\widehat{\sigma }-{\widehat{\sigma }}_1)T_1+\gamma )}$ (23)

According to Eq. (23), the solutions of Eqs. (22. a) and (22.b) can be obtained as:

$q_{12}(T_0,T_1,T_2)=\frac{2\widehat{\eta }}{{\omega }_2(2{\omega }_1+{\omega }_2)}{A}_1{A}_2{e}^{i({\omega }_1+{\omega }_2)T_0}+\frac{2\widehat{\eta }}{{\omega }_2({\omega }_2-2{\omega }_1)}{A}_1{\overline{A}}_2{e}^{i({\omega }_1-{\omega }_2)T_0}+cc$ (24.a)

$q_{22}(T_0,T_1,T_2)=\frac{\widehat{\eta }}{4{\omega }_1^2-{\omega }_2^2}{A}_1^2{e}^{2i{\omega }_1{T}_0}+\frac{\widehat{\eta }}{{\omega }_2^2}{A}_2^2{e}^{2i{\omega }_2{T}_0}-\frac{\widehat{\eta }{A}_1{\overline{A}}_1}{{\omega }_2^2}-\frac{3\widehat{\delta }{A}_2{\overline{A}}_2}{{\omega }_2^2}+cc$ (24.b)

Inserting Eqs. (19.a), (19.b), (24.a), and (24.b) into Eqs. (18.a) and (18.b), we get the following solvability conditions for Eqs. (18.a) and (18.b):

$\begin{array}{c}2i{\omega }_1{D}_2{A}_1=-D_1^2{A}_1-(i{\omega }_1{\widehat{\mu }}_1+{\widehat{\alpha }}_1+i{\omega }_1{\widehat{\alpha }}_2)A_1+(\frac{4{\widehat{\eta }}^2}{{\omega }_2^2}-\frac{2{\widehat{\eta }}^2}{4{\omega }_1^2-{\omega }_2^2}-3\widehat{\beta }-3{\widehat{\alpha }}_3-i{\omega }_1{\widehat{\alpha }}_4-{\omega }_1^2{\widehat{\alpha }}_5)A_1^2{\overline{A}}_1\\ +(\frac{12{\widehat{\eta }}^2}{{\omega }_2^2}-\frac{4{\widehat{\eta }}^2}{{\omega }_2({\omega }_2-2{\omega }_1)}-\frac{4{\widehat{\eta }}^2}{{\omega }_2(2{\omega }_1+{\omega }_2)}-2\widehat{\beta })A_1{A}_2{\overline{A}}_2-(\frac{2{\widehat{\eta }}^2}{{\omega }_2^2}+\frac{4{\widehat{\eta }}^2}{{\omega }_2({\omega }_2-2{\omega }_1)}+\widehat{\beta }){\overline{A}}_1{A}_2^2{e}^{2i{\widehat{\sigma }}_1{T}_1}\\ -\frac{1}{4}{\widehat{\delta }}_1{\overline{A}}_1{e}^{2i\widehat{\sigma }{T}_1}+\frac{i}{4}{\widehat{\delta }}_1{\overline{A}}_2{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}-\frac{3}{4}{\widehat{\delta }}_2{A}_1{\overline{A}}_1^2{e}^{2i\widehat{\sigma }{T}_1}-\frac{1}{4}{\widehat{\delta }}_2{A}_1^3{e}^{-2i\widehat{\sigma }{T}_1}+\frac{3i}{8}{\widehat{\delta }}_2{\overline{A}}_1^2{A}_2{e}^{i(2\widehat{\sigma }+{\widehat{\sigma }}_1)T_1}\\ -\frac{3i}{8}{\widehat{\delta }}_2{A}_1^2{A}_2{e}^{i(-2\widehat{\sigma }+{\widehat{\sigma }}_1)T_1}+\frac{3i}{4}{\widehat{\delta }}_2{A}_1{\overline{A}}_1{\overline{A}}_2{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}+\frac{3i}{8}{\widehat{\delta }}_2{A}_2{\overline{A}}_2^2{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}-\frac{i}{8}{\widehat{\delta }}_2{A}_2^3{e}^{i(-2\widehat{\sigma }+3{\widehat{\sigma }}_1)T_1}\\ +\frac{3}{16}{\widehat{\delta }}_2{\overline{A}}_1{\overline{A}}_2^2{e}^{i(4\widehat{\sigma }-2{\widehat{\sigma }}_1)T_1}-\frac{1}{16}{\widehat{\delta }}_2{\overline{A}}_1^3{e}^{4i\widehat{\sigma }{T}_1}-\frac{i}{16}{\widehat{\delta }}_2{\overline{A}}_2^3{e}^{i(4\widehat{\sigma }-3{\widehat{\sigma }}_1)T_1}+\frac{3i}{16}{\widehat{\delta }}_2{\overline{A}}_1^2{\overline{A}}_2{e}^{i(4\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}\end{array}$ (25.a)

$\begin{array}{c}2i{\omega }_2{D}_2{A}_2=-D_1^2{A}_2-(i{\omega }_2{\widehat{\mu }}_2+{\widehat{\alpha }}_1+i{\omega }_2{\widehat{\alpha }}_2)A_2+(\frac{30{\widehat{\eta }}^2}{{\omega }_2^2}-3\widehat{\beta }-3{\widehat{\alpha }}_3-i{\omega }_2{\widehat{\alpha }}_4-{\omega }_2^2{\widehat{\alpha }}_5)A_2^2{\overline{A}}_2\\ -(\frac{4{\widehat{\eta }}^2}{{\omega }_2(2{\omega }_1+{\omega }_2)}+\frac{4{\widehat{\eta }}^2}{{\omega }_2({\omega }_2-2{\omega }_1)}-\frac{12{\widehat{\eta }}^2}{{\omega }_2^2}+2\widehat{\beta })A_1{\overline{A}}_1{A}_2-(\frac{4{\widehat{\eta }}^2}{{\omega }_2({\omega }_2-2{\omega }_1)}+\frac{6{\widehat{\eta }}^2}{4{\omega }_1^2-{\omega }_2^2}\\ +\widehat{\beta })A_1^2{\overline{A}}_2{e}^{-2i{\widehat{\sigma }}_1{T}_1}+\frac{1}{4}{\widehat{\delta }}_1{\overline{A}}_2{e}^{i(2\widehat{\sigma }-2{\widehat{\sigma }}_1)T_1}+\frac{i}{4}{\widehat{\delta }}_1{\overline{A}}_1{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}+\frac{3}{4}{\widehat{\delta }}_2{A}_2{\overline{A}}_2^2{e}^{i(2\widehat{\sigma }-2{\widehat{\sigma }}_1)T_1}\\ +\frac{1}{4}{\widehat{\delta }}_2{A}_2^3{e}^{i(2{\widehat{\sigma }}_1-2\widehat{\sigma })T_1}+\frac{3i}{8}{\widehat{\delta }}_2{A}_1{\overline{A}}_2^2{e}^{i(2\widehat{\sigma }-3{\widehat{\sigma }}_1)T_1}-\frac{3i}{8}{\widehat{\delta }}_2{A}_1{A}_2^2{e}^{i({\widehat{\sigma }}_1-2\widehat{\sigma })T_1}+\frac{3i}{4}{\widehat{\delta }}_2{\overline{A}}_1{A}_2{\overline{A}}_2{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}\\ +\frac{3i}{8}{\widehat{\delta }}_2{A}_1{\overline{A}}_1^2{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}-\frac{i}{8}{\widehat{\delta }}_2{A}_1^3{e}^{i(-2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}+\frac{3}{16}{\widehat{\delta }}_2{\overline{A}}_1^2{\overline{A}}_2{e}^{i(4\widehat{\sigma }-2{\widehat{\sigma }}_1)T_1}-\frac{1}{16}{\widehat{\delta }}_2{\overline{A}}_2^3{e}^{i(4\widehat{\sigma }-4{\widehat{\sigma }}_1)T_1}\\ -\frac{3i}{16}{\widehat{\delta }}_2{\overline{A}}_1{\overline{A}}_2^2{e}^{i(4\widehat{\sigma }-3{\widehat{\sigma }}_1)T_1}+\frac{i}{16}{\widehat{\delta }}_2{\overline{A}}_1^3{e}^{i(4\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}\end{array}$ (25.b)

From equation (14), we can write

$\begin{array}{c}2i{\omega }_j\frac{d}{dt}{A}_j(T_1,T_2)=2i{\omega }_j{D}_0{A}_j(T_1,T_2)+2\epsilon i{\omega }_j{D}_1{A}_j(T_1,T_2)+2{\epsilon }^{2}i{\omega }_j{D}_2{A}_j(T_1,T_2)\\ =2\epsilon i{\omega }_j{D}_1{A}_j(T_1,T_2)+2{\epsilon }^{2}i{\omega }_j{D}_2{A}_j(T_1,T_2),j=\mathrm{1,2.}\end{array}$ (26)

Using Eq. (26), we can combine the secular terms in Eqs. (23), (25.a) and (25.b) as follows:

$\begin{array}{c}2i{\omega }_1\frac{d}{dt}{A}_1=\frac{\epsilon \widehat{f}{\Omega }^2}{2}{e}^{i(\widehat{\sigma }{T}_1+\gamma )}-\frac{{\epsilon }^2\widehat{\sigma }\widehat{f}{\Omega }^2}{4{\omega }_1}{e}^{i(\widehat{\sigma }{T}_1+\gamma )}-{\epsilon }^2(i{\omega }_1{\widehat{\mu }}_1+{\widehat{\alpha }}_1+i{\omega }_1{\widehat{\alpha }}_2)A_1+{\epsilon }^2(\frac{4{\widehat{\eta }}^2}{{\omega }_2^2}-\frac{2{\widehat{\eta }}^2}{4{\omega }_1^2-{\omega }_2^2}\\ -3\widehat{\beta }-3{\widehat{\alpha }}_3-i{\omega }_1{\widehat{\alpha }}_4-{\omega }_1^2{\widehat{\alpha }}_5)A_1^2{\overline{A}}_1+{\epsilon }^2(\frac{12{\widehat{\eta }}^2}{{\omega }_2^2}-\frac{4{\widehat{\eta }}^2}{{\omega }_2({\omega }_2-2{\omega }_1)}-\frac{4{\widehat{\eta }}^2}{{\omega }_2(2{\omega }_1+{\omega }_2)}-2\widehat{\beta })A_1{A}_2{\overline{A}}_2\\ -{\epsilon }^2(\frac{2{\widehat{\eta }}^2}{{\omega }_2^2}+\frac{4{\widehat{\eta }}^2}{{\omega }_2({\omega }_2-2{\omega }_1)}+\widehat{\beta }){\overline{A}}_1{A}_2^2{e}^{2i{\widehat{\sigma }}_1{T}_1}-\frac{1}{4}{\epsilon }^2{\widehat{\delta }}_1{\overline{A}}_1{e}^{2i\widehat{\sigma }{T}_1}+\frac{i}{4}{\epsilon }^2{\widehat{\delta }}_1{\overline{A}}_2{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}\\ -\frac{3}{4}{\epsilon }^2{\widehat{\delta }}_2{A}_1{\overline{A}}_1^2{e}^{2i\widehat{\sigma }{T}_1}-\frac{1}{4}{\epsilon }^2{\widehat{\delta }}_2{A}_1^3{e}^{-2i\widehat{\sigma }{T}_1}+\frac{3i}{8}{\epsilon }^2{\widehat{\delta }}_2{\overline{A}}_1^2{A}_2{e}^{i(2\widehat{\sigma }+{\widehat{\sigma }}_1)T_1}-\frac{3i}{8}{\epsilon }^2{\widehat{\delta }}_2{A}_1^2{A}_2{e}^{i(-2\widehat{\sigma }+{\widehat{\sigma }}_1)T_1}\\ +\frac{3i}{4}{\epsilon }^2{\widehat{\delta }}_2{A}_1{\overline{A}}_1{\overline{A}}_2{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}+\frac{3i}{8}{\epsilon }^2{\widehat{\delta }}_2{A}_2{\overline{A}}_2^2{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}-\frac{i}{8}{\epsilon }^2{\widehat{\delta }}_2{A}_2^3{e}^{i(-2\widehat{\sigma }+3{\widehat{\sigma }}_1)T_1}\\ +\frac{3}{16}{\epsilon }^2{\widehat{\delta }}_2{\overline{A}}_1{\overline{A}}_2^2{e}^{i(4\widehat{\sigma }-2{\widehat{\sigma }}_1)T_1}-\frac{1}{16}{\epsilon }^2{\widehat{\delta }}_2{\overline{A}}_1^3{e}^{4i\widehat{\sigma }{T}_1}-\frac{i}{16}{\epsilon }^2{\widehat{\delta }}_2{\overline{A}}_2^3{e}^{i(4\widehat{\sigma }-3{\widehat{\sigma }}_1)T_1}+\frac{3i}{16}{\epsilon }^2{\widehat{\delta }}_2{\overline{A}}_1^2{\overline{A}}_2{e}^{i(4\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}\end{array}$ (27.a)

$\begin{array}{c}2i{\omega }_2\frac{d}{dt}{A}_2=-\frac{i\epsilon \widehat{f}{\Omega }^2}{2}{e}^{i((\widehat{\sigma }-{\widehat{\sigma }}_1)T_1+\gamma )}-\frac{i{\epsilon }^2({\widehat{\sigma }}_1-\widehat{\sigma })\widehat{f}{\Omega }^2}{4{\omega }_2}{e}^{i((\widehat{\sigma }-{\widehat{\sigma }}_1)T_1+\gamma )}-{\epsilon }^2(i{\omega }_2{\widehat{\mu }}_2+{\widehat{\alpha }}_1+i{\omega }_2{\widehat{\alpha }}_2)A_2\\ +{\epsilon }^2(\frac{30{\widehat{\eta }}^2}{{\omega }_2^2}-3\widehat{\beta }-3{\widehat{\alpha }}_3-i{\omega }_2{\widehat{\alpha }}_4-{\omega }_2^2{\widehat{\alpha }}_5)A_2^2{\overline{A}}_2-{\epsilon }^2(\frac{4{\widehat{\eta }}^2}{{\omega }_2(2{\omega }_1+{\omega }_2)}+\frac{4{\widehat{\eta }}^2}{{\omega }_2({\omega }_2-2{\omega }_1)}-\frac{12{\widehat{\eta }}^2}{{\omega }_2^2}\\ +2\widehat{\beta })A_1{\overline{A}}_1{A}_2-{\epsilon }^2(\frac{4{\widehat{\eta }}^2}{{\omega }_2({\omega }_2-2{\omega }_1)}+\frac{6{\widehat{\eta }}^2}{4{\omega }_1^2-{\omega }_2^2}+\widehat{\beta })A_1^2{\overline{A}}_2{e}^{-2i{\widehat{\sigma }}_1{T}_1}+\frac{1}{4}{\epsilon }^2{\widehat{\delta }}_1{\overline{A}}_2{e}^{i(2\widehat{\sigma }-2{\widehat{\sigma }}_1)T_1}\\ +\frac{i}{4}{\epsilon }^2{\widehat{\delta }}_1{\overline{A}}_1{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}+\frac{3}{4}{\epsilon }^2{\widehat{\delta }}_2{A}_2\overline{A_2^2}{e}^{i(2\widehat{\sigma }-2{\widehat{\sigma }}_1)T_1}+\frac{1}{4}{\epsilon }^2{\widehat{\delta }}_2{A}_2^3{e}^{i(2{\widehat{\sigma }}_1-2\widehat{\sigma })T_1}+\frac{3i}{8}{\epsilon }^2{\widehat{\delta }}_2{A}_1{\overline{A}}_2^2{e}^{i(2\widehat{\sigma }-3{\widehat{\sigma }}_1)T_1}\\ -\frac{3i}{8}{\epsilon }^2{\widehat{\delta }}_2{A}_1{A}_2^2{e}^{i({\widehat{\sigma }}_1-2\widehat{\sigma })T_1}+\frac{3i}{4}{\epsilon }^2{\widehat{\delta }}_2{\overline{A}}_1{A}_2{\overline{A}}_2{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}+\frac{3i}{8}{\epsilon }^2{\widehat{\delta }}_2{A}_1{\overline{A}}_1^2{e}^{i(2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}\\ -\frac{i}{8}{\epsilon }^2{\widehat{\delta }}_2{A}_1^3{e}^{i(-2\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}+\frac{3}{16}{\epsilon }^2{\widehat{\delta }}_2{\overline{A}}_1^2{\overline{A}}_2{e}^{i(4\widehat{\sigma }-2{\widehat{\sigma }}_1)T_1}-\frac{1}{16}{\epsilon }^2{\widehat{\delta }}_2{\overline{A}}_2^3{e}^{i(4\widehat{\sigma }-4{\widehat{\sigma }}_1)T_1}\\ -\frac{3i}{16}{\epsilon }^2{\widehat{\delta }}_2{\overline{A}}_1{\overline{A}}_2^2{e}^{i(4\widehat{\sigma }-3{\widehat{\sigma }}_1)T_1}+\frac{i}{16}{\epsilon }^2{\widehat{\delta }}_2{\overline{A}}_1^3{e}^{i(4\widehat{\sigma }-{\widehat{\sigma }}_1)T_1}\end{array}$ (27.b)

To analyze Eqs. (27.a) and (27.b), we can express $A_1$ and$A_2$ in polar form as:

$A_1=\frac{1}{2}{a}_1(t)e^{i{\theta }_1(t)}\Rightarrow \frac{d}{dt}{A}_1=\frac{1}{2}{\dot{a}}_1(t)e^{i{\theta }_1(t)}+\frac{i}{2}{a}_1(t){\dot{\theta }}_1(t)e^{i{\theta }_1(t)}$ (28.a)

$A_2=\frac{1}{2}{a}_2(t)e^{i{\theta }_2(t)}\Rightarrow \frac{d}{dt}{A}_2=\frac{1}{2}{\dot{a}}_2(t)e^{i{\theta }_2(t)}+\frac{i}{2}{a}_2(t){\dot{\theta }}_2(t)e^{i{\theta }_2(t)}$ (28.b)

where$a_1$and$a_2$are the system oscillation amplitudes in the horizontal and vertical directions, while ${\theta }_1$and ${\theta }_2$ are the phase-angles. Inserting Eqs. (28.a) and (28.b) into Eqs. (27.a) and (27.b), and then separating the real and imaginary parts, with the coming back to the original system parameters (i.e.${\widehat{\mu }}_1=\frac{{\mu }_1}{{\epsilon }^2},{\widehat{\mu }}_2=\frac{{\mu }_2}{{\epsilon }^2},\widehat{\beta }=\frac{\beta }{{\epsilon }^2},\widehat{\eta }=\frac{\eta }{\epsilon },\widehat{f}=\frac{f}{\epsilon },{\widehat{\alpha }}_n=\frac{\alpha }{{\epsilon }^2},\left(n=\mathrm{1,2,}\dots \mathrm{,5}\right)$), we can obtain the following first-order nonlinear autonomous dynamical system.

$\begin{array}{c}{\dot{a}}_1=(1-\frac{\sigma }{2{\omega }_1})\frac{f{\Omega }^2}{2{\omega }_1}\sin ({\varphi }_1+\gamma )-\frac{1}{2}({\mu }_1+{\alpha }_2)a_1-\frac{1}{8}{\alpha }_4{a}_1^3-\frac{1}{8{\omega }_1}(\frac{2{\eta }^2}{{\omega }_2^2}+\frac{4{\eta }^2}{{\omega }_2({\omega }_2-2{\omega }_1)}\\ +\beta )a_1{a}_2^2\sin (2{\varphi }_1-2{\varphi }_2)-\frac{1}{8{\omega }_1}{\delta }_1{a}_1\sin (2{\varphi }_1)+\frac{1}{8{\omega }_1}{\delta }_1{a}_2\cos ({\varphi }_1+{\varphi }_2)-\frac{3}{32{\omega }_1}{\delta }_2{a}_1^3\sin (2{\varphi }_1)\\ +\frac{1}{32{\omega }_1}{\delta }_2{a}_1^3\sin (2{\varphi }_1)+\frac{3}{64{\omega }_1}{\delta }_2{a}_1^2{a}_2\cos (3{\varphi }_1-{\varphi }_2)-\frac{3}{64{\omega }_1}{\delta }_2{a}_1^2{a}_2\cos ({\varphi }_1+{\varphi }_2)\\ +\frac{3}{32{\omega }_1}{\delta }_2{a}_1^2{a}_2\cos ({\varphi }_1+{\varphi }_2)+\frac{3}{64{\omega }_1}{\delta }_2{a}_2^3\cos ({\varphi }_1+{\varphi }_2)-\frac{1}{64{\omega }_1}{\delta }_2{a}_2^3\cos ({\varphi }_1-3{\varphi }_2)\\ +\frac{3}{128{\omega }_1}{\delta }_2{a}_1{a}_2^2\sin (2{\varphi }_1+2{\varphi }_2)-\frac{1}{128{\omega }_1}{\delta }_2{a}_1^3\sin (4{\varphi }_1)-\frac{1}{128{\omega }_1}{\delta }_2{a}_2^3\cos ({\varphi }_1+3{\varphi }_2)\\ +\frac{3}{128{\omega }_1}{\delta }_2{a}_1^2{a}_2\cos (3{\varphi }_1+{\varphi }_2)\end{array}$ (29.a)

$\begin{array}{c}{\dot{a}}_2=(\frac{\sigma }{2{\omega }_2}-\frac{{\sigma }_1}{2{\omega }_2}-1)\frac{f{\Omega }^2}{2{\omega }_2}\cos ({\varphi }_2+\gamma )-\frac{1}{2}({\mu }_2+{\alpha }_2)a_2-\frac{1}{8}{\alpha }_4{a}_2^3-\frac{1}{8{\omega }_2}(\frac{4{\eta }^2}{{\omega }_2({\omega }_2-2{\omega }_1)}+\frac{6{\eta }^2}{4{\omega }_1^2-{\omega }_2^2}\\ +\beta )a_1^2{a}_2\sin (-2{\varphi }_1+2{\varphi }_2)+\frac{1}{8{\omega }_2}{\delta }_1{a}_2\sin (2{\varphi }_2)+\frac{1}{8{\omega }_2}{\delta }_1{a}_1\cos ({\varphi }_1+{\varphi }_2)+\frac{3}{32{\omega }_2}{\delta }_2{a}_2^3\sin (2{\varphi }_2)\\ -\frac{1}{32{\omega }_2}{\delta }_2{a}_2^3\sin (2{\varphi }_2)+\frac{3}{64{\omega }_2}{\delta }_2{a}_1{a}_2^2\cos (3{\varphi }_2-{\varphi }_1)-\frac{3}{64{\omega }_2}{\delta }_2{a}_1{a}_2^2\cos ({\varphi }_1+{\varphi }_2)\\ +\frac{3}{32{\omega }_2}{\delta }_2{a}_1{a}_2^2\cos ({\varphi }_1+{\varphi }_2)+\frac{3}{64{\omega }_2}{\delta }_2{a}_1^3\cos ({\varphi }_1+{\varphi }_2)-\frac{1}{64{\omega }_2}{\delta }_2{a}_1^3\cos ({\varphi }_2-3{\varphi }_1)\\ +\frac{3}{128{\omega }_2}{\delta }_2{a}_1^2{a}_2\sin (2{\varphi }_1+2{\varphi }_2)-\frac{1}{128{\omega }_2}{\delta }_2{a}_2^3\sin (4{\varphi }_2)-\frac{3}{128{\omega }_2}{\delta }_2{a}_1{a}_2^2\cos ({\varphi }_1+3{\varphi }_2)\\ +\frac{1}{128{\omega }_2}{\delta }_2{a}_1^3\cos (3{\varphi }_1+{\varphi }_2)\end{array}$ (29.b)

$\begin{array}{c}{\dot{\varphi }}_1=\sigma +(1-\frac{\sigma }{2{\omega }_1})\frac{f{\Omega }^2}{2{\omega }_1{a}_1}\cos ({\varphi }_1+\gamma )-\frac{{\alpha }_1}{2{\omega }_1}+\frac{1}{8{\omega }_1}(\frac{4{\eta }^2}{{\omega }_2^2}-\frac{2{\eta }^2}{4{\omega }_1^2-{\omega }_2^2}-3\beta -3{\alpha }_3-{\omega }_1^2{\alpha }_5)a_1^2\\ +\frac{1}{8{\omega }_1}(\frac{12{\eta }^2}{{\omega }_2^2}-\frac{4{\eta }^2}{{\omega }_2({\omega }_2-2{\omega }_1)}-\frac{4{\eta }^2}{{\omega }_2(2{\omega }_1+{\omega }_2)}-2\beta )a_2^2-\frac{1}{8{\omega }_1}(\frac{2{\eta }^2}{{\omega }_2^2}+\frac{4{\eta }^2}{{\omega }_2({\omega }_2-2{\omega }_1)}\\ +\beta )a_2^2\cos (2{\varphi }_1-2{\varphi }_2)-\frac{1}{8{\omega }_1}{\delta }_1\cos (2{\varphi }_1)-\frac{1}{8{\omega }_1{a}_1}{\delta }_1{a}_2\sin ({\varphi }_1+{\varphi }_2)-\frac{3}{32{\omega }_1}{\delta }_2{a}_1^2\cos (2{\varphi }_1)\\ -\frac{1}{32{\omega }_1}{\delta }_2{a}_1^2\cos (2{\varphi }_1)-\frac{3}{64{\omega }_1}{\delta }_2{a}_1{a}_2\sin (3{\varphi }_1-{\varphi }_2)-\frac{3}{64{\omega }_1}{\delta }_2{a}_1{a}_2\sin ({\varphi }_1+{\varphi }_2)\\ -\frac{3}{32{\omega }_1}{\delta }_2{a}_1{a}_2\sin ({\varphi }_1+{\varphi }_2)-\frac{3}{64{\omega }_1{a}_1}{\delta }_2{a}_2^3\sin ({\varphi }_1+{\varphi }_2)+\frac{1}{64{\omega }_1{a}_1}{\delta }_2{a}_2^3\sin ({\varphi }_1-3{\varphi }_2)\\ +\frac{3}{128{\omega }_1}{\delta }_2{a}_2^2\cos (2{\varphi }_1+2{\varphi }_2)-\frac{1}{128{\omega }_1}{\delta }_2{a}_1^2\cos (4{\varphi }_1)+\frac{1}{128{\omega }_1{a}_1}{\delta }_2{a}_2^3\sin ({\varphi }_1+3{\varphi }_2)\\ -\frac{3}{128{\omega }_1}{\delta }_2{a}_1{a}_2\sin (3{\varphi }_1+{\varphi }_2)\end{array}$ (29.c)

$\begin{array}{c}{\dot{\varphi }}_2=\sigma -{\sigma }_1-(\frac{\sigma }{2{\omega }_2}-\frac{{\sigma }_1}{2{\omega }_2}-1)\frac{f{\Omega }^2}{2{\omega }_2{a}_2}\sin ({\varphi }_2+\gamma )-\frac{{\alpha }_1}{2{\omega }_2}+\frac{1}{8{\omega }_2}(\frac{30{\eta }^2}{{\omega }_2^2}-3\beta -3{\alpha }_3-{\omega }_2^2{\alpha }_5)a_2^2\\ -\frac{1}{8{\omega }_2}(\frac{4{\eta }^2}{{\omega }_2(2{\omega }_1+{\omega }_2)}+\frac{4{\eta }^2}{{\omega }_2({\omega }_2-2{\omega }_1)}-\frac{12{\eta }^2}{{\omega }_2^2}+2\beta )a_1^2-\frac{1}{8{\omega }_2}(\frac{4{\eta }^2}{{\omega }_2({\omega }_2-2{\omega }_1)}+\frac{6{\eta }^2}{4{\omega }_1^2-{\omega }_2^2}\\ +\beta )a_1^2\cos (-2{\varphi }_1+2{\varphi }_2)+\frac{1}{8{\omega }_2}{\delta }_1\cos (2{\varphi }_2)-\frac{1}{8{\omega }_2{a}_2}{\delta }_1{a}_1\sin ({\varphi }_1+{\varphi }_2)+\frac{3}{32{\omega }_2}{\delta }_2{a}_2^2\cos (2{\varphi }_2)\\ +\frac{1}{32{\omega }_2}{\delta }_2{a}_2^2\cos (2{\varphi }_2)-\frac{3}{64{\omega }_2}{\delta }_2{a}_1{a}_2\sin (3{\varphi }_2-{\varphi }_1)-\frac{3}{64{\omega }_2}{\delta }_2{a}_1{a}_2\sin ({\varphi }_1+{\varphi }_2)\\ -\frac{3}{32{\omega }_2}{\delta }_2{a}_1{a}_2\sin ({\varphi }_1+{\varphi }_2)-\frac{3}{64{\omega }_2{a}_2}{\delta }_2{a}_1^3\sin ({\varphi }_1+{\varphi }_2)+\frac{1}{64{\omega }_2{a}_2}{\delta }_2{a}_1^3\sin ({\varphi }_2-3{\varphi }_1)\\ +\frac{3}{128{\omega }_2}{\delta }_2{a}_1^2\cos (2{\varphi }_1+2{\varphi }_2)-\frac{1}{128{\omega }_2}{\delta }_2{a}_2\cos (4{\varphi }_2)+\frac{3}{128{\omega }_2}{\delta }_2{a}_1{a}_2\sin ({\varphi }_1+3{\varphi }_2)\\ -\frac{1}{128{\omega }_2{a}_2}{\delta }_2{a}_1^3\sin (3{\varphi }_1+{\varphi }_2)\end{array}$ (29.d)

where ${\varphi }_1=\sigma t-{\theta }_1$ and ${\varphi }_2=\left(\sigma -{\sigma }_1\right)t-{\theta }_2$. Now, by inserting Eqs. (28.a) and (28.b) into Eqs. (19.a), (19.b), (24.a) and (24.b), and then inserting the obtained results into Eqs. (13.a) and (13.b), we get the solutions of the system original Eqs. (12.a) and (12.b) as follows:

$q_1(t)=a_1(t)\cos \left(\Omega t-{\varphi }_1(t)\right)$ (30.a)

$q_2(t)=a_2(t)\cos \left(\Omega t-{\varphi }_2(t)\right)$ (30.b)

Depending on Eqs. (30.a) and (30.b), $a_1\left(t\right)$ and$a_2\left(t\right)$ are the instantaneous lateral vibration amplitudes of in $X$ and $Y$ directions, while${\varphi }_1\left(t\right)$ and ${\varphi }_2\left(t\right)$ are the instantaneous phase-angles. At the steady-state motion, we have ${\dot{a}}_1={\dot{a}}_2={\dot{\varphi }}_1={\dot{\varphi }}_2=0$. Accordingly, by setting ${\dot{a}}_1={\dot{a}}_2={\dot{\varphi }}_1={\dot{\varphi }}_2=0$ into Eqs. (29.a) to (29.d), we can get four nonlinear algebraic equations that govern the steady-state vibration amplitudes and the corresponding phase-angles. So, by solving the resulting nonlinear algebraic equations simultaneously using $\sigma$ as the bifurcation parameter, one can plot the different response curves given in Sec. 4. In addition, to investigate the steady-state solution local stability, we let $a_{10},a_{20},{\varphi }_{10},and{\varphi }_{20}$ be the solution of Eqs. (29.a) to (29.d) when ${\dot{a}}_1={\dot{a}}_2={\dot{\varphi }}_1={\dot{\varphi }}_2=0$, and $a_{11},a_{21},{\varphi }_{11},{\varphi }_{21}$ is small perturbation about that solution. Accordingly, we assume

$\begin{array}{l}{a}_1=a_{10}+a_{11},a_2=a_{20}+a_{21},{\varphi }_1={\varphi }_{10}+{\varphi }_{11},{\varphi }_2={\varphi }_{20}+{\varphi }_{21}\Rightarrow \\ {\dot{a}}_1={\dot{a}}_{11},{\dot{a}}_2={\dot{a}}_{21},{\dot{\varphi }}_1={\dot{\varphi }}_{11},{\dot{\varphi }}_2={\dot{\varphi }}_{21}\end{array}\}$ (31)

Inserting Eq. (31) into Eqs. (29.a) to (29.d) and then obtaining the corresponding linear system, we have

$\left[\begin{array}{c}{\dot{a}}_{11}\\ {\dot{a}}_{21}\\ {\dot{\varphi }}_{11}\\ {\dot{\varphi }}_{21}\end{array}\right]=\left[\begin{array}{cccc}{\Gamma }_{11}& {\Gamma }_{12}& {\Gamma }_{13}& {\Gamma }_{14}\\ {\Gamma }_{21}& {\Gamma }_{22}& {\Gamma }_{23}& {\Gamma }_{24}\\ {\Gamma }_{31}& {\Gamma }_{32}& {\Gamma }_{33}& {\Gamma }_{34}\\ {\Gamma }_{41}& {\Gamma }_{42}& {\Gamma }_{43}& {\Gamma }_{44}\end{array}\right]\left[\begin{array}{c}{a}_{11}\\ a_{21}\\ {\varphi }_{11}\\ {\varphi }_{21}\end{array}\right]$ (32)

According to the Hartman-Grobman theorem, the nonlinear system given by Eqs. (29.a) to (29.d) is topologically equivalent to the linear system (32) if their equilibrium point is hyperbolic. Therefore, we can obtain the following characteristic equation:

${\lambda }^4+{\Delta }_1{\lambda }^3+{\Delta }_2{\lambda }^2+{\Delta }_3\lambda +{\Delta }_4=0$ (33)

where $\lambda$ is the eigenvalues of the linearized system (32). According to Routh-Hurwitz criterion, the necessary and sufficient conditions for the solutions of Eqs. (29.a) to (29.d) to be asymptotically stable are:

${\Delta }_1>\mathrm{0,}{\Delta }_1{\Delta }_2-{\Delta }_3>\mathrm{0,}{\Delta }_3({\Delta }_1{\Delta }_2-{\Delta }_3)-{\Delta }_1^2{\Delta }_4>\mathrm{0,}{\Delta }_4>0$ (34)

3. Oscillatory behaviours investigation

This section is intended to investigate the oscillatory behaviours of the asymmetric rotating shaft given by Eqs. (12.a) and (12.b), and to explore the performance of the proposed controller in mitigating the rotor nonlinear vibrations, which in turn eliminates the nonlinear characteristics of such systems. Based on the obtained slow-flow modulating equations (i.e. Eqs. (29.a) to (29.d)), the different system spinning-speed response-curves are obtained via solving Eqs. (29.a) to (29.d) (when${\dot{a}}_1={\dot{a}}_2={\dot{\varphi }}_1={\dot{\varphi }}_2=0$), utilizing the detuning parameter $\sigma$ as the main bifurcation parameter. Accordingly, plotting the amplitudes ($a_1$ &$a_2$) versus the detuning parameter$\sigma$, one can predict the steady-state lateral vibrations of the considered system. Besides, plotting the phase-angles (${\varphi }_1$ &${\varphi }_2$) against $\sigma$, we can determine the whirling direction either forward, or backward depending on the relation between ${\varphi }_1$ and ${\varphi }_2$. It is easy to show from Eqs. (30.a) and (30.b) that the rotating shaft can exhibit forward whirling as long as${\varphi }_2>{\varphi }_1$, while at ${\varphi }_2<{\varphi }_1$ the system will exhibit backward whirling motion. In addition, if${\varphi }_1={\varphi }_2$, the rotating disk will oscillate along a straight line.

Moreover, the local stability of the obtained solution has been checked according to the Routh-Hurwitz criterion given by Eq. (34). Besides the obtained analytical solutions (i.e. solution of Eqs. (29.a) to (29.d) when${\dot{a}}_1={\dot{a}}_2={\dot{\varphi }}_1={\dot{\varphi }}_2=0$), numerical validations for the obtained response-curves have been performed via solving the system original equations (i.e. Eqs. (12.a) to (12.d)) numerically using the standard Matlab ODE45 solver. The numerical results are plotted as small-circles when sweeping the bifurcation parameter $\sigma$ to the right, and as big-dots when sweeping the same bifurcation parameter to the left. The system parameters values are selected as in Refs. [¹⁵, ¹⁶, ²⁷-²⁹] such that $f=\frac{e_d}{g_0}=0.025,{\mu }_1=\frac{c_1}{\sqrt{k_{1}m}}=0.015,{\mu }_2=\frac{c_2}{\sqrt{k_{1}m}}=0.025,\lambda =\frac{q_{st}}{g_0}=1,\rho =\frac{g_0^2{k}_2}{k_1}=0.05,\eta =\lambda \rho ,\beta =\rho +\frac{3}{8}{\delta }_2,{\delta }_1=\frac{{\mathrm{\Delta }k}_1}{k_1}=0.025,{\delta }_2=\frac{{g_0^2\mathrm{\Delta }k}_2}{k_1}=0.025,{\omega }_1=\sqrt{1+{\lambda }^2\rho +\frac{1}{2}{\delta }_1},{\omega }_2=\sqrt{1+3{\lambda }^2\rho +\frac{1}{2}{\delta }_1},p=\frac{g_0}{I_0}{\gamma }_1=0.0,$and$d=\frac{g_0{\omega }_n}{I_0}{\gamma }_2=0.0$, unless otherwise is mentioned. Based on the definitions of the dimensionless parameters given below Eqs. (12.a) and (12.b), it clear that ${\delta }_1$ denotes the linear asymmetric stiffness coefficient, ${\delta }_2$ represents the nonlinear asymmetric stiffness coefficient, $p$ is proportional control gain, and $d$ is derivative control gain. Accordingly, if we suppose that${\delta }_1={\delta }_2=0$, this means that we investigate the dynamical behaviours of the symmetric system, otherwise, we talk about the asymmetric one. The following subsections are organized such that Sec. 4.1 is devoted to investigating the oscillatory behaviours of both the symmetric and asymmetric system before control, while Sec.4. 2 is intended to explore the effect of the control parameters ($p$ and $d$) on the nonlinear oscillations of the asymmetric system. Finally, a comparative analysis between the controlled and uncontrolled system is introduced in Sec.4.3.

4.1 Oscillatory behaviours of both symmetric and asymmetric system

Coming back to Eq. (21), it is clear that the detuning parameter $\sigma$ denotes the difference between the system natural frequency ${\omega }_1$ and the disk spinning speed$\Omega$. Accordingly, it is possible to employ the parameter $\sigma$ to characterize the system's dynamical behaviours near the primary resonance case. Within this section, the nonlinear dynamical characteristics of both the symmetric (i.e.${\delta }_1={\delta }_2=0$) and asymmetric (i.e. ${\delta }_1\ne 0$or ${\delta }_2\ne 0$) system when the applied controller is turned off (i.e. ${\alpha }_n=\mathrm{0.0,}n=\mathrm{1,2,}\dots \mathrm{,5}$) are investigated.

Fig. 3 illustrates the spinning-speed response-curve and the corresponding phase-angles of the asymmetric system at $f=0.025$ when${\delta }_1={\delta }_2={\alpha }_n=\mathrm{0.0,}n=\mathrm{1,2,}\dots \mathrm{,5}$. Fig. 3a confirms the dominance of the nonlinearities on the response-curves, where bi-stable solutions interval has appeared for$\sigma \in \left]\mathrm{0.027,}0.0586\right[$. In addition, Fig. 3b illustrates that the symmetric system may exhibit both forward and backward whirling motion at this bi-stable solutions interval, where ${\varphi }_2>{\varphi }_1$ as $\sigma$ increases and ${\varphi }_1>{\varphi }_2$ as $\sigma$ decreases. Numerical solution for the system original equations (i.e. Eqs. (12.a) and (12.b)) according to Fig. 3 at $\sigma =0.05$ is shown in Fig. 4 at two different initial conditions. Figs. 4a and 4b show the symmetric system steady-state temporal oscillations and the corresponding whirling motion at$q_1\left(0\right)=q_2\left(0\right)={\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=0.0$, while Figs. 4c and 4d illustrate the symmetric system steady-state temporal oscillations and the corresponding whirling motion at$q_1\left(0\right)=q_2\left(0\right)=\mathrm{1.5,}{\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=1.0$. By comparing Figs. 4b and 4d, it is clear the sensitivity of the symmetric system to the initial conditions at $\sigma =0.05$ as depicted in Fig. 3, where the system performs a backward whirling at the zero initial conditions, while the forward whirling is appeared at$q_1\left(0\right)=q_2\left(0\right)=\mathrm{1.5,}{\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=1.0$.

The influence of stiffness asymmetry on the oscillatory behaviours of the horizontally supported rotor system is discussed through Figs. 5, 6, and 8. Fig. 6 illustrates the system spinning-speed response-curves and the corresponding phase-angles when ${\delta }_1=0.05$ and${\delta }_2=0$. By comparing Fig. 5 with Fig. 3, it is clear from Fig. 5 that the existence of asymmetry in the linear stiffness coefficients may be resulting in increasing the system lateral vibration amplitudes. In addition, the asymmetry of the linear stiffness coefficient increases the linear natural frequencies of the system (where${\omega }_1=\sqrt{1+{\lambda }^2\rho +\frac{1}{2}{\delta }_1},{\omega }_2=\sqrt{1+3{\lambda }^2\rho +\frac{1}{2}{\delta }_1}$), which ultimately shifts the spinning-speed response-curves to the right. Moreover, the rotating shaft system that has a linear asymmetric stiffness coefficient may exhibit complex dynamical behaviours than the symmetric one. In addition, there exist a spinning-speed interval at which the system has Tri-stable periodic solutions besides the bi-stable solutions interval $\sigma \in \left]\mathrm{0.0307,}0.0793\right[$. The coexistence of Tri-stable solutions for the asymmetric system is numerically simulated as shown in Fig. 6 according to the marked points on Fig. 5 (i.e. $p_1,p_2,$ and$p_3$) at $\sigma =0.0793$. Fig. 6 shows the steady-state lateral oscillations and the corresponding whirling orbits of the asymmetric rotor via solving Eqs. (12.a) and (12.b) at the different initial conditions according to Fig. 5 (i.e.$f=\mathrm{0.025,}{\delta }_1=\mathrm{0.05,}$ and ${\delta }_2=0.0$) when$\sigma =0.0793$. It is clear from Figs. 6a, 6b, 6e, and 7f that the system performs forward whirling motion with two different whirling orbits at the same rotational speed (i.e.$\Omega ={\omega }_1+\sigma ,\sigma =0.0793$) depending only on the shaft initial position. Also, Figs. 6c and 6d illustrate that the system periodic solutions in the vertical and horizontal directions are in-phase that are result in the shaft oscillation along a straight line when the initial conditions are $q_1\left(0\right)=q_2\left(0\right)=\mathrm{1.0,}{\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=-1.6$.

The effect of asymmetric nonlinear stiffness coefficient (${\delta }_2$) on the system spinning-speed response curve is depicted in Fig. 7, where the figure shows the system response-curves when $f=\mathrm{0.025,}{\delta }_1=\mathrm{0.0,}$ and${\delta }_2=0.05$. It is clear from Fig. 7 that the asymmetric nonlinear stiffness coefficient may cause the appearance of tri-stable solutions as long as$\Omega ={\omega }_1+\sigma ,\sigma \in \left]\mathrm{0.075,}0.12\right[$ beside the bi-stable solution in the interval$\sigma \in \left]\mathrm{0.032,}0.075\right[$. The influence of both the nonlinear and linear stiffness coefficients is illustrated in Figs. 8 that is the more practical case of the asymmetric rotating machinery. The figure shows the system response-curves when$f=\mathrm{0.025,}$ and${\delta }_1={\delta }_2=0.05$. Comparing Figs. 5, and 7 with Fig. 8, we can conclude that the asymmetric rotor system with ${\delta }_1={\delta }_2=0.05$ has a different dynamical behaviours, where the system can perform forward and backward whirling motions at the narrow interval $\sigma \in \left]\mathrm{0.036,}0.0433\right[$, otherwise, the system can perform forward whirling only.

Fig. 3
Symmetric rotor system response curve at$f=0.025$ and${\delta }_1={\delta }_2=0.0$: (a) oscillation amplitudes $a_1$ and $a_2$ at the horizontal and vertical directions, and (b) the corresponding phase-Angles.

Fig. 4
Symmetric rotor system temporal oscillations according to Fig. 3 when $\sigma =0.05$ at different initial conditions: (a, b) the system steady-state temporal oscillations, and the corresponding whirling orbit at $q_1\left(0\right)=q_2\left(0\right)={\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=0.0$, and (c, d) the system steady-state temporal oscillations, and the corresponding whirling orbit at $q_1\left(0\right)=q_2\left(0\right)=\mathrm{1.5,}{\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=1.0$.

Fig. 5
Symmetric rotor system response curves at$f=0.025$ and${\delta }_1=\mathrm{0.05,}{\delta }_2=0.0$: (a) oscillation amplitude $a_1$ in the horizontal direction, (b) oscillation amplitude $a_2$ in the vertical direction, and (c) the corresponding phase-Angles.

Fig. 6
Asymmetric rotor system temporal oscillations according to Fig. 5 when $\sigma =0.0793$ at different initial conditions: (a, b) the system steady-state temporal oscillations, and the corresponding whirling orbit at $q_1\left(0\right)=q_2\left(0\right)={\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=0.0$, (c, d) the system steady-state temporal oscillations, and the corresponding whirling orbit at $q_1\left(0\right)=q_2\left(0\right)=\mathrm{1.0,}{\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=-1.6$, and (e, f) the system steady-state temporal oscillations, and the corresponding whirling orbit at $q_1\left(0\right)=q_2\left(0\right)={\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=1.5$.

4.2 Effect of the control parameters on system response-curves

The nonlinear oscillatory behaviours of the controlled asymmetric rotating shaft are explored within this section. Fig. 9 shows the system lateral vibrations at different values of the proportional gain ($p$) when${\delta }_1={\delta }_2=0.025$. The spinning-speed response-curves of the considered system are illustrated in Figs. 9a, and 9b, for the proportional gain $p=0.9$ when the derivative gain is set to be zero (i.e.$d=0.0$), while Figs. 9c and 9d show the system response-curves at $p=0.95$ and d=0.0. Besides, the asymmetric system response-curves at $p=1.0$ and d=0.0 are presented in Figs. 9e, and 9f. Generally, it is clear from Fig. 9 that the increase of the proportional control gain beyond $p=0.95$, bent the system spinning-speed response-curves to right leading to hard spring characteristics, while decreasing the proportional gain $p$ to become lower than$0.95$, bent the system response-curves to left leading to soft spring characteristics. However, the asymmetric controlled system can exhibit unbounded oscillation amplitudes if the proportional control gain$p=0.95$ and the disk spinning-speed $\Omega ={\omega }_1+\sigma$ belong to $\sigma \in \left]-\mathrm{0.03282,}-0.0089\right[$ as it is clear from Figs. 9c and 9d. To validate the possibility of destabilizing the considered system when $p=0.95$ as reported in Figs. 9c and 9d, the system original equations (i.e. Eqs. (12.a) and (12.b)) have been solved numerically according to Fig. 9 as shown in Figs. 10 and 11. Fig.10 illustrates the temporal oscillations and the corresponding phase trajectory of the asymmetric system according to Fig. 9 at $\sigma =-0.015$ when switching the proportional control gain from $p=0.9$to$p=0.95$, while Fig.11 shows the system temporal oscillations and the corresponding phase trajectory when switching the proportional control gain from $p=1.0$to$p=0.95$. It is clear from Figs. 10 and 11 that the asymmetric system can perform bounded periodic oscillations at $\sigma =-0.015$ as long as$p=0.9$ or $p=1.0$, while switching the proportional control gain from $p=1$ or $p=0.9$ to $p=0.95$ is resulting in growth unbound oscillations of the considered system that is perfectly agreed the obtained analytical solutions in Fig .9.

Fig. 7
Asymmetric rotor system response curves at$f=0.025$ and${\delta }_1=\mathrm{0.0,}{\delta }_2=0.05$: (a) oscillation amplitude $a_1$ in the horizontal direction, (b) oscillation amplitude $a_2$ in the vertical direction, and (c) the corresponding phase-Angles.

Fig. 8
Asymmetric rotor system response curves at$f=0.025$ and${\delta }_1={\delta }_2=0.05$: (a) oscillation amplitude $a_1$ in the horizontal direction, (b) oscillation amplitude $a_2$ in the vertical direction, and (c) the corresponding phase-Angles.

The effect of the derivative gain ($d$) on the asymmetric system lateral vibrations (i.e.${\delta }_1={\delta }_2=0.025$) is explored as shown in Fig .12 for two different values of the proportional control gain. Figs. 12a and 12b show the system response-curves at the different values of the derivative gain when$p=1.0$, while Figs. 12c and 12d show the same response-curves but when$p=0.8$. Generally, the figure confirms that the increasing of the derivative control gain$d$, incases the system damping coefficients, which ultimately decreases the system lateral vibrations. Moreover, increasing the derivative gain to a specific limit (i.e. $d=0.02$ or$0.03$) can force the asymmetric complex dynamical system (as seen in Figs. 3 to 8) to behave like a linear system with a unique periodic solution. In addition. Fig. 12 confirms that the vibration mitigation efficiency of the proposed proportional-derivative controller in the case of soft spring characteristics (i.e.$p=0.8$) is higher than that of hard spring characteristics (i.e.$p=1.0$).

Fig. 9
Controlled asymmetric rotor system response curves at different values of the proportional gain $p$ at$f=\mathrm{0.025,}{\delta }_1={\delta }_2=0.025$, and$d=0.0$: (a, b) oscillation amplitudes $a_1$ and $a_2$ at the horizontal and vertical directions at $p=0.9$, (c, d) oscillation amplitudes $a_1$ and $a_2$ at the horizontal and vertical directions at $p=0.95$, and (e, f) oscillation amplitudes $a_1$ and $a_2$ at the horizontal and vertical directions at $p=1.0$.

Fig. 10
Controlled asymmetric rotor system temporal oscillations according to Fig. 9 at $\sigma =-0.015$ when changing the proportional control gain $p$ from $p=0.9$ to $p=0.95$: (a, b) the system temporal oscillations in the horizontal and vertical directions at $q_1\left(0\right)=q_2\left(0\right)={\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=0.0$, and (c) and the corresponding phase plane.

Fig. 11
Controlled asymmetric rotor system temporal oscillations according to Fig. 9 at $\sigma =-0.015$ when changing the proportional control gain $p$ from $p=1.0$ to $p=0.95$: (a, b) the system temporal oscillations in the horizontal and vertical directions at $q_1\left(0\right)=q_2\left(0\right)={\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=0.0$, and (c) and the corresponding phase plane.

Fig. 12
Controlled asymmetric rotor system response curves at different values of the derivative gain $d$ at$f=\mathrm{0.025,}{\delta }_1={\delta }_2=0.025$: (a, b) oscillation amplitudes $a_1$ and $a_2$ at the horizontal and vertical directions at$p=1.0$, and (c, d) oscillation amplitudes $a_1$ and $a_2$ at the horizontal and vertical directions at $p=0.8$.

The influence of increasing ${\delta }_1$ on the vibration amplitudes of the controlled system (when $f=\mathrm{0.025,}$ ${\delta }_2=\mathrm{0.0,}$ and$d=0.05$) is investigated as shown in Fig. 13 at two different values of the proportional control gain (i.e. $p=1$ and$0.7$). It is clear from the figure that the oscillation amplitudes are a monotonic increasing function of the linear asymmetric stiffness coefficient. Moreover, at the large values of${\delta }_1$, the nonlinear behaviours dominate again the controlled system response-curves. However, the vibration suppression efficiency of the applied controller at $p=0.7$ is higher than that at$p=1.0$, where the system oscillation amplitudes in Fig. 13a and 13.b are twice the corresponding ones in Fig. 13c and 13.d.

Fig. 14 shows the influence of increasing ${\delta }_2$ on the controlled system lateral vibration amplitudes when $f=\mathrm{0.025,}$ ${\delta }_1=\mathrm{0.15,}$ and $d=0.05$ at two different values of the proportional control gain. Figs. 14a and 14b show the system nonlinear vibrations when$p=1.0$, while Figs. 14c and 14d illustrate the system oscillatory motion when$p=0.7$. It is clear from Figs. 14c and 14d that the oscillation amplitudes of the controlled system have a negligible sensitivity to increasing${\delta }_2$, while a small increase in ${\delta }_2$, causes a huge increase in the vibration amplitudes when $p=1.0$ as shown in Figs. 14a and 14b. According to Figs.12, 13, and 14, the best proportional control gain to get a higher vibration suppression efficiency for the applied control algorithm is$p\cong 0.7$.

Fig. 13
Controlled asymmetric rotor system response curves at different values of the asymmetric linear asymmetric stiffness coefficient ${\delta }_1$ at$f=\mathrm{0.025,}{\delta }_2=\mathrm{0.0,}d=0.05$: (a, b) oscillation amplitudes $a_1$ and $a_2$ at the horizontal and vertical directions at${\delta }_1=\mathrm{0.05,}\mathrm{0.1,}\mathrm{0.15,}\mathrm{0.2,}$and$p=1.0$, and (c, d) oscillation amplitudes $a_1$ and $a_2$ at the horizontal and vertical directions at${\delta }_1=\mathrm{0.05,}\mathrm{0.1,}\mathrm{0.15,}\mathrm{0.2,}$and$p=0.7$.

Fig. 14
Controlled asymmetric rotor system response curves at different values of the asymmetric nonlinear asymmetric stiffness coefficient ${\delta }_2$ at$f=\mathrm{0.025,}{\delta }_1=\mathrm{0.15,}d=0.05$: (a, b) oscillation amplitudes $a_1$ and $a_2$ at the horizontal and vertical directions at${\delta }_2=\mathrm{0.025,}\mathrm{0.05,}0.1$and$p=1.0$, and (c, d) oscillation amplitudes $a_1$ and $a_2$ at the horizontal and vertical directions at${\delta }_2=\mathrm{0.025,}\mathrm{0.05,}\mathrm{0.1,}$and$p=0.7$.

4.3 Rub-impact forces between the pole-leg and the rotating shaft

It is important to remember that the dimensionless temporal displacements $q_1\left(t\right)$ and $q_2\left(t\right)$ represent the actual displacements (i.e. $x\left(\tau \right)$ and$y\left(\tau \right)$) of the real system divided by the air-gap size $g_0$, where $q_1\left(t\right)=\frac{x\left(\tau \right)}{g_0}$ and $q_2\left(t\right)=\frac{y\left(\tau \right)}{g_0}$. In addition, the periodic solution of Eqs. (12.a) and (12.b) is proved to be $q_1\left(t\right)=a_1\left(t\right)\cos \left(\Omega t-{\varphi }_1\left(t\right)\right)$ and$q_2\left(t\right)=a_2\left(t\right)\cos \left(\Omega t-{\varphi }_2\left(t\right)\right)$. Accordingly, if we suppose that the dimensionless oscillation amplitude $a_1\ge 1$ or/and $a_2\ge 1$, this implies that the vibration amplitude of actual displacement of $x\left(t\right)$ or/and $y\left(t\right)$ is greater than the air-gap, which in turn means the appearance of rub-impact forces between the asymmetric rotor system and the pole-leg. Accordingly, the obtained spinning-speed response-curve at which $a_1\ge 1$ or/and $a_2\ge 1$, is undesired from the engineering point of view where the rub-impact forces between the rotor and the pole-leg are inevitable as in Figs. 9, 12a, 12b, 13, 14a, and 14b.

Now, the spinning-speed response-curves of the studied system before activating the proposed controller and after control are compared in Figs. 15 and 16, when $f=\mathrm{0.025,}{\delta }_1={\delta }_2=\mathrm{0.025,}d=0.05$ and$p=\mathrm{0.7,}1.0$. It is clear from the two figures that the different nonlinear characteristics of the asymmetric rotating shaft before control have been eliminated after control and the asymmetric rotor has been forced to behave like a linear system with a unique periodic solution. Besides, the backward whirling motion of the uncontrolled system has been eliminated after control as illustrated in Figs. 16c and 16d where ${\varphi }_2$ is always greater than${\varphi }_1$. Examining Fig. 16 in-depth, we can deduce that the vibration suppression efficiency of the proportional-derivative controller in at $p=0.7$ is higher than that when$p=1.0$. Moreover, the maximum vibration amplitude of the controlled system either if $p=0.7$ or $p=1.0$ is very small compared to the unity (i.e. $a_1<1$ and $a_2<1$), which guarantees the safe operation of the controlled system without the rub-impact forces between the rotating shaft and the pole-leg.

Fig. 15
Asymmetric rotor system response curves at$f=0.025$ and${\delta }_1={\delta }_2=0.025$: (a, b) oscillation amplitudes $a_1,a_2$ in the horizontal and vertical directions, and (c, d) the corresponding phase-Angles.

Fig. 16
Controlled asymmetric rotor system response curves at$f=0.025$,${\delta }_1={\delta }_2=\mathrm{0.025,}$ and$d=0.05$: (a, b) oscillation amplitudes $a_1,a_2$ in the horizontal and vertical directions at $p=\mathrm{1,}0.7$, and (c, d) the corresponding phase-Angles.

Figs. 17 and 18 illustrate a numerical validation for the temporal vibrations of the considered system before and after control according to Figs. 15 and 16 (i.e.$f={\delta }_1={\delta }_2=0.05$) when$\sigma =0.05$. Fig. 17 simulates the asymmetric system lateral vibrations before and after control when the initial conditions are$q_1\left(0\right)=q_2\left(0\right)={\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=0.0$. The figure shows evolving the uncontrolled system lateral vibrations up to $t=1000$, where at $t=1000$ the applied controller has been turned on with control gains $p=1.0$ and $d=0.05$until the time variable became $t=1250$, at that instant (i.e. $t=1250$) the proportional control gain has been decreased from $p=1.0$ to $p=0.7$. Similarly, Fig. 18 is a repetition for Fig. 17, but the system initial conditions in this case are$q_1\left(0\right)=q_2\left(0\right)=\mathrm{0.9,}{\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=0.0$. Generally, it is clear from Figs. 17 and 18 that the bi-stable solutions of the uncontrolled system have been merged into a single solution after control. Moreover, decreasing the proportional gain from $p=1.0$ to $p=0.7$ improves the vibration suppression capability of the proposed controller. Notice that the stator circumference (i.e. the electromagnetic pole-leg location) is plotted as a dashed circle to show the possibility of rub-impact forces between the rotor and the pole-leg as shown in Figs. 17d and 18d. Figs. 17d and 18d show that the controlled system can operate safely mode regardless of the system initial conditions, where the figures display the efficiency of the controller in forcing the shaft to oscillate away from the pole-leg. However, the same figures confirm the possibility of the rotor system destruction where the rub-impact forces may occur between the rotating shaft and the pole-leg if the controller is turned off.

The angle $\gamma$ denotes the angle among the eccentricity direction (i.e. $\overrightarrow{OC}$ direction) and the shaft large stiffness direction (i.e. $\overrightarrow{O{z}_1}$ direction) as shown in Fig. 1c. The effect of $\gamma$ on the nonlinear vibration of the asymmetric rotor system when $\Omega ={\omega }_1$ (i.e. at $\sigma =0$) at different values of ${\delta }_1$ and ${\delta }_2$ is explored in Figs. 19 and 20, respectively. In Fig. 19, the angle $\gamma$ is plotted against the steady-state oscillation amplitudes of the system before control at ${\delta }_1={\delta }_2=\mathrm{0.0,}\mathrm{0.025,}$ and$0.05$. The figure shows that the oscillation amplitudes are constant as $\gamma$ increases when${\delta }_1={\delta }_2=0.0$. This means that the eccentricity direction has a negligible effect on the vibration amplitudes of the symmetric rotors (i.e.${\delta }_1={\delta }_2=0.0$). However, it has a considerable effect on the vibration amplitudes of the asymmetric systems, where for ${\delta }_1={\delta }_2>0.0$ the vibration amplitudes are periodic functions of$\gamma$ with period$\pi$. Moreover, the maximum vibration amplitudes occur at$\gamma \cong {\frac{\pi }{2}}^o,{\frac{3\pi }{2}}^o$ (i.e. when $\overrightarrow{OC}$ and $\overrightarrow{O{z}_1}$ axes are perpendicular), while the minimum oscillation amplitudes occur at $\gamma \cong \mathrm{0,}{\pi }^o$(i.e. when $\overrightarrow{OC}$ and $\overrightarrow{O{z}_1}$ axes are in the same direction or the opposite direction). Fig. 20 is a repetition for Fig. 19, but when $p=0.1$ and $d=0.05$ (i.e. for the controlled system). By comparing Fig. 19 with Fig. 20, we can conclude the capability of the applied control algorithm in reducing the fluctuations of the orientation angle on steady-state lateral vibrations.

Fig. 17
Temporal oscillations, the corresponding phase plane, and whirling orbit of the asymmetric rotor system before and after control at initial conditions $q_1\left(0\right)=q_2\left(0\right)={\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=0.0$ when$\sigma =0.05$ according to Figs. 15 and 16: (a, c) the system temporal oscillations in the horizontal and vertical directions, (b) the phase plane, and (d) the steady-state whirling orbit.

Fig. 18
Temporal oscillations, the corresponding phase plane, and whirling orbit of the asymmetric rotor system before and after control at initial conditions $q_1\left(0\right)=q_2\left(0\right)=\mathrm{0.9,}{\dot{q}}_1\left(0\right)={\dot{q}}_2\left(0\right)=0.0$ when$\sigma =0.05$ according to Figs. 15 and 16: (a, c) the system temporal oscillations in the horizontal and vertical directions, (b) the phase plane, and (d) the steady-state whirling orbit.

Fig. 19
Asymmetric rotor system $\gamma$-amplitude response-curves at$f=\mathrm{0.025,}\sigma =0.0,{\delta }_1={\delta }_2=\mathrm{0.0,}\mathrm{0.025,}$and$0.05$: (a) oscillation amplitude $a_1$ in the horizontal direction, and (b) oscillation amplitude $a_2$ in the vertical direction.

Fig. 20
Controlled asymmetric rotor system $\gamma$-amplitude response-curves at$f=\mathrm{0.025,}\sigma =\mathrm{0.0,}p=\mathrm{1.0,}d=\mathrm{0.05,}{\delta }_1={\delta }_2=\mathrm{0.0,}\mathrm{0.025,}$and$0.05$: (a) oscillation amplitude $a_1$ in the horizontal direction, and (b) oscillation amplitude $a_2$ in the vertical direction.

5. Conclusion

Within this article, a linear proportional-derivative controller is applied to mitigate the nonlinear vibrations of an asymmetric nonlinear rotor system. The suggested controller is coupled to the rotor system utilizing an electromagnetic actuator consisting of four poles, where one pair of these poles is responsible for controlling the system vibrations in the horizontal direction, while the other pair is designed to control the system oscillations in the vertical direction. The four electromagnetic poles are energized by electrical currents governed by the proposed control algorithm to generate controllable magnetic forces in both the horizontal and vertical directions. Accordingly, the whole system mathematical model is derived as a nonlinear dynamical system with quadratic and cubic nonlinearities and having both external and multi-parametric excitations. Applying the perturbation analysis, four autonomous first-order nonlinear differential equations that govern the system vibration amplitudes ($a_1$&$a_2$) and the corresponding phase-angles (${\varphi }_1$&${\varphi }_2$) are obtained. Based on the obtained autonomous equations, the oscillatory behaviours of the asymmetric rotor before and after control have been explored via obtaining the different bifurcation diagrams. According to the introduced analyses, the following points can be concluded:

1. Based on the excellent agreements between the obtained analytical and numerical solutions, it is possible to predict both the whirling orbit shape and its direction for the studied complex system accurately by obtaining the spinning-speed response-curve and the corresponding phase-angles.

2. The studied nonlinear symmetric rotor system is sensitive to the initial conditions at a specific range of the spinning-speeds, where the system can execute one of two oscillatory motions (the first one is a forward whirling motion and the other is a backward motion).

3. The asymmetric nonlinear rotor system can exhibit more complicated dynamical behaviours than the symmetric one, where the asymmetric system may perform a single-stable solution, or bi-stable solutions, or tri-stable solution according to the shaft initial conditions and its angular velocity ($\Omega$).

4. The obtained results approved the feasibility of the proposed control strategy in forcing the asymmetric rotor system to respond as a linear one, where all undesired nonlinear phenomena such as the sensitivity to initial conditions, the coexistence of both forward and backward whirling motions, bi-stable solutions, and tri-stable solutions have been eliminated after control.

5. The proportional control gain $p$ should be selected within the range of [0.7, 0.8] to enhance the vibration suppression efficiency of the suggested controller.

6. The asymmetric controlled system can operate safely as long as the applied controller is working properly under the designed optimum control gains (i.e. $p=0.7$ and $d=0.0$5). However, the abrupt failure for the applied controller may be resulting in rub-impact forces occurrence between the rotor and the pole-leg, which ultimately leads to the whole system destruction.

Acknowledgments

• The authors are thankful of the Taif University. Taif University researchers supporting project number (TURSP-2020/160), Taif University, Taif, Saudi Arabia.

References

• ^[1]
  Ardayfio, D., Frohrib, D. A.: Instabilities of an asymmetric rotor with asymmetric shaft mounted on symmetric elastic supports. Journal of Engineering for Industry 98(4), 1161-1165 (1976).
• ^[2]
  Iwatsubo, T., Tsujiuchi, Y., Inouev, T.: Vibration of asymmetric rotor supported by oil film bearings. Archive of Applied Mechanics 56(1), 1-15 (1986).
• ^[3]
  Park, J.: Diagnosis of Excessive Vibration Signals of Two-Pole Generator Rotors in Balancing. KSME International Journal, 14(6), 590-596 (2000)
• ^[4]
  Hsieh, S., Chen, J., Lee, A.: A modified transfer matrix method for the coupled lateral and torsional vibrations of asymmetric rotor-bearing systems. Journal of Sound and Vibration 312 (4-5), 563-571 (2008).
• ^[5]
  Ishida, Y., Liu, J.: Elimination of unstable ranges of rotors utilizing discontinuous spring characteristics: an asymmetrical shaft system, an asymmetrical rotor system, and a rotor system with liquid. J. Vib. Acoust 132(1), 011011 (2010) (8 pages).
• ^[6]
  Shahgholi, M., Khadem, S. E.: Primary and parametric resonances of asymmetrical rotating shafts with stretching nonlinearity. Mechanism and Machine Theory 51, 131-144 (2012).
• ^[7]
  Han, Q., Chu, F.:The effect of transverse crack upon parametric instability of a rotor-bearing system with an asymmetric disk. Communications in Nonlinear Science and Numerical Simulation 17(12), 5189-5200 (2012).
• ^[8]
  Han, Q., Chu, F.: Parametric instability of a Jeffcott rotor with rotationally asymmetric inertia and transverse crack. Nonlinear Dyn 73(1-2), 827-842 (2013).
• ^[9]
  Saeed, N. A., Eissa, M.: Bifurcations of periodic motion of a horizontally supported nonlinear Jeffcott-rotor system having transversely cracked shaft. International Journal of Non-Linear Mechanics 101, 113-130‏ (2018).
• ^[10]
  Saeed, N. A., Eissa, M.: Bifurcation Analysis of a Transversely Cracked Nonlinear Jeffcott-rotor System at Different Resonance Cases. International Journal of Acoustics and Vibration 24 (2), 284-302 (2019).
• ^[11]
  Xiang, L., Gao, X., Hu, A.: Nonlinear dynamics of an asymmetric rotor-bearing system with coupling faults of crack and rub-impact under oil-film Forces. Nonlinear Dyn 86(2), 1057-1067 (2016).
• ^[12]
  Przybylowicz, P. M., Starczewski, Z., Korczak-Komorowski, P.: Sensitivity of regions of irregular and chaotic vibrations of an asymmetric rotor supported on journal bearings to structural parameters. Acta Mech 227 (11), 3101-3112 (2016).
• ^[13]
  Yu, T., Zhou, S.,Yang, X., Zhang, W.: Global dynamics of a flexible asymmetrical rotor. Nonlinear Dyn 91 (2), 1041-1060 (2018).
• ^[14]
  Boru, F.E.: Numerical and Experimental Response and Stability Investigations of Anisotropic Rotor-Bearing Systems. Kassel University Press GmbH, Kassel (2010).
• ^[15]
  Saeed, N. A.: On vibration behaviour and motion bifurcation of a nonlinear asymmetric rotating shaft. Arch Appl Mech 89, 1899-1921 (2019).
• ^[16]
  Saeed, N. A.: On the steady-state forward and backward whirling motion of asymmetric nonlinear rotor system. . European Journal of Mechanics - A/Solids 80, 103878 (2020).
• ^[17]
  Tehrani, G., G., Dardel, M.: Mitigation of nonlinear oscillations of a Jeffcott rotor System with an optimized damper and nonlinear energy sink. International Journal of Non-Linear Mechanics 98, 122-136 (2018).
• ^[18]
  Tehrani, G., G., Dardel, M., Pashaei, M. H.: Passive vibration absorbers for vibration reduction in the multi-bladed rotor with rotor and stator contact. Acta Mech 231, 597-623 (2020).
• ^[19]
  Eissa, M, Saeed, N. A, El-Ganini, W.A.: Saturation-based active controller for vibration suppression of a four-degree-of-freedom rotor-AMB system. Nonlinear Dynamics 76(1), 743-764 (2014).
• ^[20]
  Saeed, N. A., Kandil, Al.: Lateral vibration control and stabilization of the quasiperiodic oscillations for rotor-active magnetic bearings system. Nonlinear Dynamics 98 (2), 1191-1218 (2019).
• ^[21]
  Kandil, A., Sayed, M., Saeed, N. A.: On the nonlinear dynamics of constant stiffness coefficients 16-pole rotor active magnetic bearings system. European Journal of Mechanics-A/Solids 84,104051 (2020).
• ^[22]
  Saeed, N. A, Mahrous, E., Awrejcewicz, J.: Nonlinear dynamics of the six-pole rotor-AMB system under two different control configurations. Nonlinear Dynamics101 (4), 2299-2323 (2020).
• ^[23]
  Saeed, N. A; Awwad, E. M., El-Meligy, M. A, Abouel-Nasr, E.: Radial versus Cartesian control strategies to stabilize the nonlinear whirling motion of the six-pole rotor-AMBs. IEEE Access 8, 138859-138883 (2020).
• ^[24]
  Ishida, Y., Inoue, T.: Vibration Suppression of nonlinear rotor systems using a dynamic damper. Journal of Vibration and Control 13(8),1127-1143 (2007)
• ^[25]
  Saeed, N. A., Kamel, M.: Nonlinear PD-controller to suppress the nonlinear oscillations of horizontally supported Jeffcott-rotor system. International Journal of Non-Linear Mechanics 87, 109-124 (2016).
• ^[26]
  Ghasabi, S. A., Arbabtafti, M., Shahgholi, M.: Time-delayed control of a nonlinear asymmetrical rotor near the major critical speed with flexible supports. Mechanics Based Design of Structures and Machines, (2020). DOI: 10.1080/15397734.2020.1715230.
  » https://doi.org/10.1080/15397734.2020.1715230
• ^[27]
  Saeed, N. A; Awwad, E. M., El-Meligy, M. A, Abouel-Nasr, E.: Sensitivity analysis and vibration control of asymmetric nonlinear rotating shaft system utilizing 4-pole AMBs as an actuator. European Journal of Mechanics - A/Solids 86, 104145, (2021).
• ^[28]
  Yabuno, H., Kashimura, T., Inoue, T., Ishida, Y.: Nonlinear normal modes and primary resonance of horizontally supported Jeffcott-rotor. Nonlinear Dyn 66 (3), 377-387 (2011).
• ^[29]
  Saeed, N. A., El-Gohary, H. A.: On the nonlinear oscillations of a horizontally supported Jeffcott-rotor with a nonlinear restoring force. Nonlinear Dyn 88(1), 293-314 (2017).
• ^[30]
  Ishida, Y., Yamamoto, T.: Linear and Nonlinear Rotordynamics. Wiley‐VCH (2012).
• ^[31]
  Bleuler, H., Cole, M., Keogh, P., Larsonneur, R., Maslen, E., Nordmann, R., Okada, Y., Schweitzer, G., Traxler, A.: Magnetic Bearings Theory, Design, and Application to Rotating Machinery. Springer, Berlin Heidelberg (2009).
• ^[32]
  Nayfeh, A., Mook, D.: Nonlinear Oscillations. Wiley, New York (1995).

Nomenclature

• $q_1,\dot{q_1},\ddot{q_1}$  Lateral displacement, velocity, and acceleration in the horizontal direction.
• $q_2,\dot{q_2},\ddot{q_2}$  lateral displacement, velocity, and acceleration in the vertical direction.
• ${\mu }_1,{\mu }_2$  Linear damping in the horizontal and vertical directions, respectively.
• ${\omega }_1,{\omega }_2$  Natural frequencies in the horizontal and vertical directions, respectively.
• $\eta$  Quadratic nonlinearity coefficient.
• $\beta$  Cubic nonlinearity coefficient.
• ${\delta }_1$  Linear asymmetric stiffness coefficient.
• ${\delta }_2$  Nonlinear asymmetric stiffness coefficient.
• $f$  Shaft eccentricity.
• $\Omega$  Shaft spinning-speed.
• $p$  Proportional control gain.
• $d$  Derivative control gain.
• $a_1,a_2$  Steady-state vibration amplitudes in the horizontal and vertical direction, respectively.
• ${\varphi }_1,{\varphi }_2$  Steady-state phase-angles of the system motion in the horizontal and vertical direction, respectively.

Edited by

Editor:

Marcílio Alves.

Publication Dates

• Publication in this collection
  19 Mar 2021
• Date of issue
  2021

History

• Received
  21 Jan 2021
• Reviewed
  28 Jan 2021
• Accepted
  28 Jan 2021
• Published
  05 Feb 2021