Jump and Initial-Sensitive Excessive Motion of a Class of Relative Rotation Systems and Their Control via Delayed Feedback

Abstract

Jump and excessive motion are undesirable phenomena in relative rotation systems, causing a loss of global integrity and reliability of the systems. In this work, a typical relative rotation system is considered in which jump, excessive motion, and their suppression via delayed feedback are investigated. The Method of Multiple Scales and the Melnikov method are applied to analyze critical conditions for bi-stability and initial-sensitive excessive motion, respectively. By introducing the fractal of basins of attraction and the erosion of the safe basin to depict jump and initial-sensitive excessive motion, respectively, the point mapping approach is used to present numerical simulations which are in agreement with the theoretical prediction, showing the validity of the analysis. It is found that jump between bistable attractors can be due to saddle–node bifurcation, while initial-sensitive excessive motion can be due to heteroclinic bifurcation. Under a positive coefficient of the gain, the types of delayed feedback can both be effective in reducing jump and initial-sensitive excessive motion. The results may provide some reference for the performance improvement of rotors and main bearings.

1. Introduction

Rotating machinery, widely used in engineering fields such as transmission systems of vehicles [1] and main bearings of aeroengines [2], plays an important role in the theory of nonlinear dynamics. Due to the nonlinearities of its dynamical systems, it usually undergoes complex dynamics, which have a negative impact on its operational stability and reliability [3,4]. Among these, relative rotation systems have been given much attention, as nonlinear vibration is undesirable. For instance, main shafts can be fatigue damaged by a jump [5,6] among multiple attractors, or by a chaotic response [7], and can even be fractured due to excessive motion [8]. Thus, comprehensive study of the global dynamics of relative rotation systems has become a practical issue.

So far, many researchers have studied the responses of relative rotation systems. Wang et al. [9] obtained the approximate periodic solution of a class of nonlinear systems with relative rotation by applying the Method of Multiple Scales and proved the uniqueness of the periodic response. Xiao et al. [10] verified the existence and uniqueness of periodic response in another type of relative rotation system with a time-varying stiffness. Li et al. [11] proposed a relative rotation system with nonlinear elastic force and nonlinear generalized damping force, and discussed its periodic solution problem theoretically. Shi et al. [12] numerically studied the chaotic behaviors of a relative rotation system under parametric excitation by Poincaré map and maximal Lyapunov exponent. By simulating relative torsional vibrations of an imbalanced shaft under a limited power supply, Verichev [13] found that the interaction of shaft and power supply may result in strange attractors such as classical Lorenz and Feigenbaum attractors. In a nonlinear relative rotation system with a triple-well Mathieu-Duffing oscillator, Liu et al. [14] obtained the threshold of chaos regarding Smale horseshoe commutation and exhibited the erosion process of safe basins. A class of nonlinear relative rotational systems containing two rotors was also built whose chaotic response was presented [15]. By introducing the erosion of the safe basin to describe the safe performance of the spur gear pair, Zhu et al. [16] classified the multiple meshing states and presented the transition process from safe to unsafe. To control the complex dynamics of the main transmission system of a scraper conveyor, Ju et al. [17] analyzed its local bifurcations and proposed a nonlinear state feedback controller whose effect was studied numerically. Considering the effectiveness of time-delay feedback on controlling fractal erosion of the safe basin and chaos in nonlinear dynamical systems [18,19], Zhao et al. [20] applied delayed displacement feedback in a relative torsional vibration system for reducing its response amplitude. On this basis, Shang et al. [21] discussed the effect of delayed position feedback on controlling the erosion of the safe basin and chaos. Although there has been meaningful research on the nonlinear vibration characteristics of the relative rotation system, the study of excessive motion and jump among multiple attractors is mainly carried out by numerical simulation and the mechanism of these complex dynamical behaviors and their control is still not yet clear.

To this end, we select a class of nonlinear relative rotation systems composed of two rotors and study the mechanisms behind its jump and initial-sensitive excessive motion, as well as the effect of delayed feedback on suppressing these phenomena. The paper is arranged as follows. In Section 2, the dynamical model of a relative rotation system is constructed and made dimensionless. In Section 3, the mechanism of jump and excessive motion is analyzed. In Section 4, two control strategies, namely, delay position feedback and delay velocity feedback, are applied to the original system, respectively, whose control mechanism is then discussed. Section 5 contains the discussion.

2. Dynamical Model and Unperturbed Dynamics

We consider a typical torsional vibrating system whose simplified diagram is presented in Figure 1. ${\theta }_1$ and ${\theta }_2$ are the absolute rotating angles of two rotors, respectively, and $J_1$ and $J_2$ represent their polar moments of inertia. Considering nonlinear torsional stiffness and nonlinear damping, the governing equation of this vibrating system can be derived by the momentum theorem as

$$\begin{array}{l}{J}_1{\ddot{\theta }}_1+c_{12}({\dot{\theta }}_1-{\dot{\theta }}_2)+K_1({\theta }_1-{\theta }_2)+K_3{({\theta }_1-{\theta }_2)}^3+f_{12}=T_{e1},\\ J_2{\ddot{\theta }}_2+c_{12}({\dot{\theta }}_2-{\dot{\theta }}_1)+K_1({\theta }_2-{\theta }_1)+K_3{({\theta }_2-{\theta }_1)}^3-f_{12}=T_{e2},\end{array}$$

(1)

where $K_1$ is the coefficient of linear stiffness between the two rotors; $K_3$ is the coefficient of nonlinear stiffness satisfying $K_3<0$; $c_{12}$ is the coefficient of linear damping; $T_{e1}$ and $T_{e2}$ are external rotational torques loaded on the two rotors, respectively; $f_{12}$ is a class of nonlinear stick-slip frictions [20,21] given by

$$f_{12}=c_0+c_1({\dot{\theta }}_1-{\dot{\theta }}_2)+c_2{({\dot{\theta }}_1-{\dot{\theta }}_2)}^3,$$

(2)

in which $c_0$, $c_1$ and $c_2$ are the constant of static friction, the coefficient of linear damping and the coefficient of cubic nonlinear damping, respectively. By denoting $x={\theta }_1-{\theta }_2$ in Equation (1), the relative rotation system can be obtained as below:

$$\ddot{x}+\frac{c_{12}{J}_2+c_{12}{J}_1}{J_1{J}_2}\dot{x}+\frac{K_1(J_1+J_2)}{J_1{J}_2}x+\frac{K_3(J_1+J_2)}{J_1{J}_2}{x}^3+\frac{J_1+J_2}{J_1{J}_2}(c_1\dot{x}+c_2{\dot{x}}^3)=\frac{J_2{T}_{e1}-J_1{T}_{e2}-(J_1+J_2)c_0}{J_1{J}_2}.$$

(3)

The torques $T_{e1}$ and $T_{e2}$ are usually harmonic excitations, thus can be represented as

$$\frac{J_2{T}_{e1}-J_1{T}_{e2}-(J_1+J_2)c_0}{J_1{J}_2}=F\cos (\Omega t)$$

(4)

where $F$ and $\Omega$ are the amplitude and the frequency of the excitation, respectively. Introducing the following variables,

$${\omega }_0{}^2=\frac{K_1(J_1+J_2)}{J_1{J}_2},b=-\frac{K_3}{K_1},\mu =\frac{(c_{12}+c_1)(J_1+J_2)}{J_1{J}_2{\omega }_0},g=\frac{c_2(J_1+J_2){\omega }_0}{J_1{J}_2},T={\omega }_{0}t,\omega =\frac{\Omega }{{\omega }_0},f=\frac{F}{{\omega }_0{}^2},$$

(5)

the dimensionless form of the relative rotation system (3) is obtained as

$$y=\frac{dx}{dT},\hspace{1em}\frac{dy}{dT}=-\mu y-x+b{x}^3-g{y}^3+f\cos (\omega T).$$

(6)

in which the parameters $\mu$, b, g, and $\omega$ are positive, and the position $x(T)$ and the velocity $y(T)$ represent relative rotational angular and relative angular velocity at the moment T, respectively.

Since the parameters $c_{12}$, $c_0$, $c_1$, and $c_2$ in the original system (1) are very small, the relevant terms concerned, $\mu$, g and f in Equation (6), can be considered as perturbed. The unperturbed system can be written as

$$\frac{dx}{dT}=y,\hspace{0.17em}\frac{dy}{dT}=-x+b{x}^3.$$

(7)

There will be three equilibria in the unperturbed system (7), namely, the center S₁(0,0) and two saddle points S₂($\frac{\sqrt{b}}{b}$,0) and S₃(−$\frac{\sqrt{b}}{b}$,0). The Hamiltonian of Equation (7) is

$$H(x,y)=\frac{1}{2}{y}^2+\frac{1}{2}{x}^2-\frac{1}{4}b{x}^4.$$

(8)

Accordingly, the heteroclinic orbits surrounding the center S₁(0,0) can be given by

$$x_{\pm }(t)=\pm \frac{\sqrt{b}}{b}\tanh (\frac{\sqrt{2}}{2}t),y_{\pm }(t)=\pm \frac{\sqrt{2b}}{2b}{\mathrm{sech}}^2(\frac{\sqrt{2}}{2}t).$$

(9)

Fixing b = 0.3, the unperturbed orbits are depicted in Figure 2.

3. Complex Dynamics

3.1. Multistability and Jump

The periodic solution near the equilibrium S₁(0,0) is considered. The Method of Multiple Scales (MMS) [22] is employed to obtain the approximate theoretical solution of Equation (3). Supposing a small parameter ε (0 < ε≪ 1) to rescale the nondimensional parameters as follows:

$$\mu =\epsilon \tilde{\mu },\hspace{0.17em}b=\epsilon \tilde{b},\hspace{0.17em}g=\epsilon \tilde{g},f=\epsilon \tilde{f},$$

(10)

Equation (6) becomes

$$\frac{d^{2}x}{d{T}^2}+x=\epsilon (-\frac{\tilde{\mu }dx}{dT}+\tilde{b}{x}^3-\tilde{g}{(\frac{dx}{dT})}^3+\tilde{f}\cos (\omega T)).$$

(11)

Further rescaling parameters and operators in Equation (11) yields

$$\omega =1+\epsilon \sigma ,\hspace{0.17em}{T}_i={\epsilon }^{i}T,\hspace{0.17em}x(T)=x_0(T_0,T_1,\cdots )+\epsilon x_1(T_0,T_1,\cdots )+\cdots ,\hspace{0.17em}{D}_i=\frac{\partial }{\partial T_i},\hspace{0.17em}\frac{d}{dT}={\displaystyle \sum }_{i=0}^n{\epsilon }^i{D}_i\hspace{0.17em}\left(i=0,1,2\cdots \right),$$

(12)

where $\sigma$ is the detuning parameter, and separating the coefficients of ${\epsilon }^0$ and ${\epsilon }^1$ of the system (11) one finds

$${\epsilon }^0:\hspace{1em}{D}_0{}^2{x}_0+x_0=0,$$

(13)

and

$${\epsilon }^1:\hspace{1em}{D}_0{}^2{x}_1+x_1=-2D_1{D}_0{x}_0-\tilde{\mu }{D}_0{x}_0+\tilde{b}{x}_0{}^3-\tilde{g}{(D_0{x}_0)}^3+\tilde{f}\cos (T_0+\sigma T_1).$$

(14)

The solution of Equation (13) is

$$x_0(T_0,T_1)=A(T_1)e^{i{T}_0}+cc$$

(15)

where cc denotes the complex conjugate of the preceding terms. Substituting Equation (15) into Equation (14) and separating the secular terms yields

$$D_{1}A+\frac{\tilde{\mu }A}{2}+\frac{3(\tilde{g}+i\tilde{b})A^2\overline{A}}{2}+\frac{i\tilde{f}{e}^{i\sigma T_1}}{4}=0.$$

(16)

In addition, by defining

$$A(T_1)=\frac{a(T_1)}{2}{e}^{i\beta (T_1)},\hspace{1em}\psi =\sigma T_1-\beta ,$$

(17)

separating the real and imaginary parts of Equation (16), and returning the parameters to the nondimensional parameters of the system (3), the modulation equation can be written as

$$\frac{da}{dT}=-\frac{\mu a}{2}-\frac{3g{a}^3}{8}+\frac{f\sin \psi }{2},\hspace{1em}a\frac{d\psi }{dT}=a(\omega -1)+\frac{3b{a}^3}{8}+\frac{f\cos \psi }{2}.$$

(18)

The steady-state form of Equation (18) is obtained by setting $D_{1}a=D_1\psi =0$ and returning the parameters of (18) to the nondimensional parameters as:

$$\mu a+\frac{3g{a}^3}{4}=f\sin \psi ,\hspace{1em}-2(\omega -1)a-\frac{3b{a}^3}{4}=f\cos \psi .$$

(19)

Accordingly, the approximate periodic solution can be expressed as $x_0=a\cos (\omega T-\psi )$. By eliminating $\psi$ from the relationships in Equation (19), its amplitude $a$ can be solved from the following equation

$${(\mu +\frac{3g{a}^2}{4})}^2{a}^2+{(2\omega -2+\frac{3}{4}b{a}^2)}^2{a}^2=f^2.$$

(20)

The corresponding characteristic equation is

$$(\lambda +\frac{\mu }{2}+\frac{9g{a}^2}{8})(\lambda +\frac{\mu }{2}+\frac{3g{a}^2}{8})+(\omega -1+\frac{9}{8}b{a}^2)(\omega -1+\frac{3}{8}b{a}^2)=0,$$

(21)

meaning that for

$$\frac{27(g^2+b^2)}{16}{a}^4+3(g\mu +b(2\omega -2))a^2+4{(\omega -1)}^2+{\mu }^2>0,$$

(22)

the periodic solution is unstable and becomes a saddle. Namely, when the inequation (22) becomes an equation, saddle-node bifurcation will occur in the system (6). For

$$\mu =0.03,b=0.3,g=0.03,$$

(23)

the response curve of Equation (20) is plotted in Figure 3 where the coexistence of two stable periodic branches in the same range of the excitation parameters implies bi-stability. For example, for f = 0.13 and ω within the range (0.65, 0.86), bistable periodic attractors coexist (see Figure 3a). Fixing ω = 0.80, the increase of the amplitude f can also lead to saddle-node bifurcation, and thus bi-stability. The numerical results totally match the theoretical prediction.

3.2. Initial-Sensitive Excessive Motion

Excessive motion of the system (6) can be understood as an escape or unbounded solution [21], for instance, a relative rotating angle in system (6) more than the material of the structure can bear. This may be induced by the variation of the system parameters or the initial conditions of the system. In the latter case, excessive motion coexists with bounded, meaning that a tiny disturbance of initial conditions leads to a change of dynamical behavior, such as from bounded vibration to unbounded. This is a typical initial-sensitive phenomenon, initial-sensitive excessive motion, indicating the loss of global integrity of dynamic systems [3,8]. Since it is often due to global bifurcation [15,21], we apply the Melnikov method [23] to detect the conditions for it. Substituting the heteroclinic orbits (9) into the Melnikov function yields

$$M_{\pm }(t_0)={{\displaystyle \int }}_{-\infty }^{+\infty }(-\mu y(T)-g{y}^3(T)+f\cos (\omega (T+t_0)))y(T)dT=-\frac{2\sqrt{2}\mu }{3b}-\frac{8\sqrt{2}g}{35b^2}+\sqrt{\frac{2}{b}}\pi \omega f\mathrm{csch}(\frac{\sqrt{2}\pi \omega }{2})\cos (\omega t_0).$$

(24)

For

$$\sqrt{\frac{2}{b}}\pi \omega f\mathrm{csch}(\frac{\sqrt{2}\pi \omega }{2})\ge \frac{2\sqrt{2}\mu }{3b}+\frac{8\sqrt{2}g}{35b^2},$$

(25)

Namely,

$$f\ge f_0=\frac{2(35\mu b+12g)\sinh (\frac{\pi \omega }{\sqrt{2}})}{105\pi \omega b^{\frac{3}{2}}},$$

(26)

there will be a value of $t_0$ satisfying $M_{\pm }(t_0)=0$ and $M_{\pm }{}^{\prime }(t_0)\ne 0$, meaning that the roots of Equation $M_{\pm }(t_0)=0$ are simple, enabling the existence of the transverse heteroclinic orbits, and the system (6) may undergo initial-sensitive motion. According to Equation (26), f₀ is the amplitude threshold of heteroclinic bifurcation. For parameter values given by Equation (23), the critical value for heteroclinic bifurcation is f₀ = 0.09. In contrast, it follows from Inequation (25) that the threshold of the excitation frequency ω for heteroclinic bifurcation cannot be expressed as an explicit or monotonic increasing function of other parameters. For example, for f = 0.13 and the values of parameters in Equation (3), it can be calculated from Equation (25) that for $\omega$ < 1.22, heteroclinic bifurcation will occur in the system (6).

In order to verify the criterion obtained in this subsection, numerical simulations are carried out by fixing parameter values as in Equation (23). The 4th Runge-Kutta approach and the point-mapping method [24] are employed to describe the phenomenon’s initial-sensitive excessive motion. First, some terms such as basin of attraction and safe basin are introduced briefly. A basin of attraction is defined as the set of initial conditions that can lead to the same attractor [25]. If the boundary of basin of attraction of an attractor is fractal and mixed with another, jump among multiple attractors may occur [25]. Safe basin is defined as the union of basins of attraction for all bounded attractors [26]. Fractal of the safe basin of the system (6) induces the occurrence of initial-sensitive excessive motion. In this paper, the basin of attraction is drawn in the sufficiently large space region defined as −3.0 ≤ x(0) ≤ 3.0, −3.0 ≤ y(0) ≤ 3.0 by generating a 600 × 600 array of starting conditions, for each of these starting points. The escaping set for infinite time is approximated with good accuracy by a study with 1000 excited circles. The time step is taken as 0.01. The white region represents the numerical approximation to the basins of attraction of excessive motion. The red and blue regions are the basins of attraction for attractors from lower and higher stable periodic brunches of Figure 3, respectively. Thus, the union of the red region and the blue is the so-called safe basin.

The evolution of basins of attraction with the frequency ω for f = 0.13 and with the amplitude f for ω = 0.80 can be observed in Figure 4 and Figure 5. In Figure 4, each the safe basins, i.e., the union of red region and blue, is fractal-eroded when ω ranges from 0.62 to 0.87 satisfying ω < 1.22 i.e., the condition of heteroclinic bifurcation, which shows the occurrence of initial-sensitive excessive motion, in agreement with the analytical prediction. In Figure 4a, there is only a red region, meaning that there is only one periodic attractor from the lower stable branch. When ω = 0.69 (see Figure 4b), the red region and the blue one coexists, showing bistability. Specifically, on the left side of the origin, the red region is fractal and mixed with the white region and the blue, indicating that both jump and initial-sensitive excessive motion may occur there. As ω increases to 0.81 (see Figure 4c), most areas of the red region will be eroded by the blue one, showing that most of the safe initial conditions lead to the higher-amplitude attractor. For ω = 0.87 (see Figure 4d), the red region disappears, and there is only the blue, which means that the response shifts from the lower-amplitude branch to the higher-amplitude, in agreement with the theoretical results of Figure 3a.

In Figure 5, the safe basin (the union of red and blue regions) is steadily eroded by the white region with the increase of the amplitude f, implying that initial-sensitive excessive motion can be induced by the increase of f. In Figure 5a,b, the boundary of safe basins is smooth. When f is more than 0.09, i.e., the critical threshold solved from Equation (26), the boundary of safe basin will become fractal (see Figure 5c,d). In Figure 5a–c, the red region and the blue one are mixed, indicating bi-stability as well as jump. Note that for f = 0.12, the red and blue regions are evenly mixed in the neighborhood of the origin (see Figure 5c). For f = 0.20, only the basin of attraction of the attractors in the higher stable branch of Figure 3b exists, which is seriously eroded by the white region, showing that the excessive motion is easy to trigger.

4. Time-Delay Feedback to Control Complex Dynamics

In this section, two types of linear delay feedback controllers, i.e., delayed position feedback and delayed velocity feedback, are applied to the system (1), and their effect on suppressing the phenomena jump and initial-sensitive excessive motion is discussed. The corresponding control diagram is shown in Figure 6 where Q is the coefficient of the gain, and τ time delay. $P({\theta }_1,{\theta }_2,{\dot{\theta }}_1,{\dot{\theta }}_2,{\theta }_{1\tau },{\theta }_{2\tau },{\dot{\theta }}_{1\tau },{\dot{\theta }}_{2\tau })$ represents delayed relative-angle feedback $({\theta }_1(t-\tau )-{\theta }_2(t-\tau ))-({\theta }_1(t)-{\theta }_2(t))$ or delayed relative-angular-velocity feedback $({\dot{\theta }}_1(t-\tau )-{\dot{\theta }}_2(t-\tau ))-({\dot{\theta }}_1(t)-{\dot{\theta }}_2(t))$, which is an active control transferring the motion state/velocity of the dynamical system [27]. The corresponding governing equation is

$$\begin{array}{l}{J}_1{\ddot{\theta }}_1+c_{12}({\dot{\theta }}_1-{\dot{\theta }}_2)+K_1({\theta }_1-{\theta }_2)+K_3{({\theta }_1-{\theta }_2)}^3+f_{12}-QP({\theta }_1,{\theta }_2,{\dot{\theta }}_1,{\dot{\theta }}_2,{\theta }_{1\tau },{\theta }_{2\tau },{\dot{\theta }}_{1\tau },{\dot{\theta }}_{2\tau })=T_{e1},\\ J_2{\ddot{\theta }}_2+c_{12}({\dot{\theta }}_1-{\dot{\theta }}_2)+K_1({\theta }_1-{\theta }_2)+K_3{({\theta }_1-{\theta }_2)}^3+f_{12}+QP({\theta }_1,{\theta }_2,{\dot{\theta }}_1,{\dot{\theta }}_2,{\theta }_{1\tau },{\theta }_{2\tau },{\dot{\theta }}_{1\tau },{\dot{\theta }}_{2\tau })=T_{e2}.\end{array}$$

(27)

Letting

$$G_p=\frac{J_1+J_2}{J_1{J}_2}Q,\hspace{1em}{G}_v=\frac{J_1+J_2}{J_1{J}_2}Q{\omega }_0,\hspace{1em}\tilde{\tau }={\omega }_0\tau ,$$

(28)

and based on Equations (2)–(6), we obtain the relative rotation system under delayed position feedback control, i.e.,

$$\frac{dx}{dT}=y,\hspace{1em}\frac{dy}{dT}=-\mu y-x+b{x}^3-g{y}^3+f\cos (\omega T)+G_p(x(T-\tilde{\tau })-x(T)),$$

(29)

and the following delayed-velocity-feedback controlled system

$$y=\frac{dx}{dT},\hspace{1em}\frac{dy}{dT}=-\mu y-x+b{x}^3-g{y}^3+f\cos (\omega T)+G_v(y(T-\tilde{\tau })-y).$$

(30)

In the delayed-feedback-controlled systems (29) and (30), $G_p$, $G_v$ and $\tilde{\tau }$ are independent parameters. For $G=0$ or $\tilde{\tau }=0$, the feedback terms in Equations (29) and (30) become 0, and the delayed systems (29) and (30) become the uncontrolled non-dimensional system (6). In this paper, considering the engineering application, we do not consider the periodic characteristics of $\tilde{\tau }$ but restrict that $0\le \tilde{\tau }\le 2\pi$. Since there is no signal returned to the controlled systems (29) and (30) before $T=0$, it is supposed that the initial states of the delayed system when $-\tilde{\tau }\le T<0$ satisfy $x(T)=y(T)=0$. Then, safe basins of the delayed-feedback controlled systems can also be projected into the initial-state plane $x(0)-y(0)$, similar to the uncontrolled system (6).

4.1. Delayed Position Feedback Control

4.1.1. Primary Resonant Response and Stability of Solutions

As discussed in Section 3.1, the Method of Multiple Scales is applied to obtain the periodic response of the delayed system. Rescaling the gain as $G_p=\epsilon {\tilde{G}}_p$, and following the similar procedure as in Section 3.1, one can finally obtain the slow flow equation:

$$\frac{da}{dT}=-\frac{\mu a}{2}-\frac{3g{a}^3}{8}+\frac{f\sin \psi }{2}-\frac{a{G}_p\sin \tilde{\tau }}{2},\hspace{1em}a\frac{d\psi }{dT}=(\omega -1)a+\frac{3b{a}^3}{8}+\frac{f\cos \psi }{2}-\frac{a{G}_p(1-\cos \tilde{\tau })}{2}.$$

(31)

One can get the amplitude and the phases of the response by equating the right sides of Equation (31) to zero, i.e.,

$$\mu a+\frac{3g{a}^3}{4}+a{G}_p\sin \tilde{\tau }=f\sin \psi ,\hspace{1em}-2(\omega -1)a-\frac{3b{a}^3}{4}+a{G}_p(1-\cos \tilde{\tau })=f\cos \psi .$$

(32)

Accordingly, the approximate periodic solution can be expressed as $x_0=a\cos (\omega T-\psi )$. Its amplitude $a$ can be solved from the following equation

$${(\mu +G_p\sin \tilde{\tau }+\frac{3g{a}^2}{4})}^2{a}^2+{(2\omega -2+\frac{3}{4}b{a}^2+G_p\cos \tilde{\tau }-G_p)}^2{a}^2=f^2.$$

(33)

The stability of the solutions, thus obtained, can be ascertained by computing the eigenvalues from the corresponding characteristic equation:

$$\begin{array}{l}{\lambda }^2+\lambda (\mu +G_p\sin \tilde{\tau }+\frac{3g{a}^2}{2})\\ +\frac{27(g^2+b^2)a^4}{64}+\frac{{(\mu +G_p\sin \tilde{\tau })}^2}{4}+{(\omega -1+\frac{G_p\cos \tilde{\tau }-G_p}{2})}^2+\frac{3a^2}{4}(g{G}_p\sin \tilde{\tau }+b{G}_p\cos \tilde{\tau }+b(2\omega -2-G_p)+g\mu )=0.\end{array}$$

(34)

Assuming $\lambda =iv$ in Equation (34) and separating the real and imaginary parts, we have

$$v(\mu +G_p\sin \tilde{\tau }+\frac{3g{a}^2}{2})=0$$

(35)

and

$$-v^2+\frac{27(g^2+b^2)a^4}{64}+\frac{{(\mu +G_p\sin \tilde{\tau })}^2}{4}+{(\omega -1+\frac{G_p\cos \tilde{\tau }-G_p}{2})}^2+\frac{3a^2}{4}(g{G}_p\sin \tilde{\tau }+b{G}_p\cos \tilde{\tau }+b(2\omega -2-G_p)+g\mu )=0.$$

(36)

As we know, when the delay is small, the delayed feedback can be expanded in a Taylor series so that the controlled system (29) will become

$$\frac{dx}{dT}=y,\hspace{1em}\frac{dy}{dT}=-(\mu +G_p\tilde{\tau })y-x+b{x}^3-g{y}^3+f\cos (\omega T).$$

(37)

For $G_p>0$, it is much like increasing the damping of the system (6), thus reducing the periodic vibration [20,27]. Letting $G_p>0$, one has $\mu +G_p\sin \tilde{\tau }+\frac{3g{a}^2}{2}>0$ for $0\le \tilde{\tau }\le \pi$. Thus, v = 0 in Equation (35); if

$$\frac{27(g^2+b^2)a^4}{64}+\frac{3a^2}{4}(g{G}_p\sin \tilde{\tau }+b{G}_p\cos \tilde{\tau }+b(2\omega -2-G_p)+g\mu )+{(\omega -1+\frac{G_p\cos \tilde{\tau }-G_p}{2})}^2+\frac{{(\mu +G_p\sin \tilde{\tau })}^2}{4}=0,$$

(38)

saddle-node bifurcation will occur in the periodic branch. For $G_p>0$, when $\tilde{\tau }$ is long enough to satisfy

$${(g{G}_p\sin \tilde{\tau }+b{G}_p\cos \tilde{\tau }+b(2\omega -2-G_p)+g\mu )}^2-3(g^2+b^2)({(\omega -1-\frac{G_p}{2}+\frac{G_p\cos \tilde{\tau }}{2})}^2+\frac{{(\mu +G_p\sin \tilde{\tau })}^2}{4})<0,$$

(39)

there will be no positive roots for Equation (38), meaning under a positive gain, the delayed position feedback can be used to reduce saddle-node bifurcation and thus bi-stability. For example, given G_p = 0.3 and ω = 0.8, one can calculate from Equation (39) that there is no saddle-node bifurcation for $\tilde{\tau }>0.63$. In contrast, it is hard to express the range of $\tilde{\tau }$ for saddle-node bifurcation as a function of f. Fixing the value of f and combining Equations (33), (38) and (39), one can obtain the condition for $\tilde{\tau }$ to reduce saddle-node bifurcation. For instance, given f = 0.13, one can calculate that the condition of $\tilde{\tau }$ to reduce saddle-node bifurcation is $\tilde{\tau }>$ 0.27.

Based on the above theoretical analysis, the results are verified by numerical simulation as shown in Figure 7. Figure 7a shows the frequency-response for f = 0.13, G_p = 0.3 and different values of time delay satisfying $\tilde{\tau }<<\pi$. Obviously, in Figure 7a, the amplitude of the response reduces with the increase in time delay. Note that for $\tilde{\tau }$ = 0.35 (longer than the theoretical values 0.27), there is no saddle-node bifurcation on the periodic branches, showing that bi-stability and jump will not occur. It can be seen from Figure 7a that the increase of time delay in the position feedback can reduce saddle-node bifurcation, and thus bi-stability and jump. This can also be verified by the evolution of basins of attraction with the increase of time delay (see Figure 8). It follows from Figure 8 that as $\tilde{\tau }$ increases, the basin of attraction of the higher-amplitude periodic attractor (the blue region) will be eroded by the basin of attraction of the lower-amplitude one (see the red region). Meanwhile, the safe basin (red and bule regions) becomes smooth. For ω = 0.80, one can see from Figure 7b that even though the saddle-node bifurcation still exists, increase of time delay can help to reduce the range of f where bistable attractors coexist. The corresponding sequences of basin of attraction with the increase of f and $\tilde{\tau }$ will be discussed in detail in the next subsection.

4.1.2. Heteroclinic Bifurcation

To employ the Melnikov method to analyze the critical condition for heteroclinic bifurcation in a delayed-feedback-controlled system, it is necessary to consider the delayed position feedback as a perturbed term, namely, the values of the delay and the gain should satisfy that the stability of the equilibrium of the uncontrolled linear system cannot be changed by the delayed feedback. Since the equilibrium S₁(0,0) is stable for μ > 0, the delay should not exceed the first value of the stability switch of the following linear delayed system

$$\frac{d^{2}x}{d{T}^2}+\frac{\mu dx}{dT}+x=G_p(x(T-\tilde{\tau })-x(T)).$$

(40)

Its characteristic equation is given by

$${\lambda }^2+\mu \lambda +1+G_p-G_p{e}^{-\lambda \tilde{\tau }}=0.$$

(41)

Substituting $\lambda =iv$ into (41) and separating the imaginary and real parts yields

$$-{\Omega }_p{}^2+1+G_p-G_p\cos ({\Omega }_p\tilde{\tau })=0,\hspace{1em}\mu {\Omega }_p+G_p\sin ({\Omega }_p\tilde{\tau })=0.$$

(42)

By eliminating the trigonometric functions, Equation (42) becomes

$${\Omega }_p^4-(2+2G_p-{\mu }^2){\Omega }_p^2+1+2G_p=0.$$

(43)

For

$$G_p>\max \{\frac{{\mu }^2}{2}-1,\frac{{\mu }^2}{2}+\mu \},$$

(44)

there are two different positive solutions of Equation (43) expressed as ${\Omega }_{p1}$ and ${\Omega }_{p2}$ (${\Omega }_{p1}>{\Omega }_{p2}>0$); the critical value of the delay for stability switch of S₁(0,0) is

$${\tau }_0=\frac{1}{{\Omega }_{p1}}(2\pi -\arccos (1+\frac{1-{\Omega }_{p1}^2}{G_p})).$$

(45)

For $\tilde{\tau }<{\tau }_0$, the delay position feedback can be regarded as a disturbed term.

Substituting the heteroclinic orbits (9) into the Melnikov function of the system (29) yields

$$M_{\pm }(t_0)=-\frac{2\sqrt{2}\mu }{3b}-\frac{8\sqrt{2}g}{35b^2}+\sqrt{\frac{2}{b}}\pi \omega f\mathrm{csch}(\frac{\pi \omega }{\sqrt{2}})\cos (\omega t_0)-\frac{G_p}{b}{l}_1(\tilde{\tau }),$$

(46)

where $l_1(\tilde{\tau })={\mathrm{csch}}^2(\frac{\sqrt{2}}{2}\tilde{\tau })(\sinh (\sqrt{2}\tilde{\tau })-\sqrt{2}\tilde{\tau })>0$. The critical condition for heteroclinic bifurcation can be expressed as

$$f>f^P(\tilde{\tau })=f_0+\frac{G_p\mathrm{sech}(\frac{\pi \omega }{\sqrt{2}})l_1(\tilde{\tau })}{\pi \sqrt{2b}\omega },$$

(47)

Accordingly, for $G_p$ > 0, the threshold of heteroclinic bifurcation $f^P(\tilde{\tau })$ will increase with time delay and be higher than f₀. Fixing µ = 0.02 and G_p = 0.3, it can be calculated that τ₀ = 2.55.

Figure 9 shows the variation of $f^P(\tilde{\tau })$ with the increase of time delay $\tilde{\tau }$ ($\tilde{\tau }$ << τ₀). Then the numerical values of $f^P(\tilde{\tau })$ are obtained at which the boundary of the safe basin begin to be unsmooth. Each numerical critical value of $f^P(\tilde{\tau })$ is kept at two decimal places. We make sure that if f is less than the numerical results $f^P(\tilde{\tau })$, the boundary of safe basin will be smooth. In Figure 9, the numerical results for the critical values of f are in agreement with the analytical values, demonstrating that the threshold of the amplitude f for heteroclinic bifurcation will increase monotonically with the delay for $G_p$ > 0 and $\tilde{\tau }$ < ${\tau }_0$.

Figure 10 depicts the evolution of the safe basin of the controlled system (29) with the delay $\tilde{\tau }$. In Figure 10a,e,i for $\tilde{\tau }=0$, safe basins as well as basins of attraction can be considered as those of the uncontrolled system (6), illustrating that as the amplitude f increases, safe basin will be eroded. Specifically, in Figure 10e for f = 0.30, only the basin of attraction of the higher-amplitude attractor is left, which is seriously eroded by the white region, showing the high possibility of excessive motion. For G_p = 0.3 and $\tilde{\tau }$ = 0.11, the safe basin is obviously expanded. Besides, comparing Figure 10a–d, one can observe that as $\tilde{\tau }$ increases, the safe basin changes from the union of red region and blue to the red region itself, meaning that jump and initial-sensitive excessive motion are reduced by the basin expansion of the lower-amplitude attractor. In Figure 10e–h, one can also draw this conclusion. Note that Figure 10e,f, three color regions are mixed in the neighborhood of the origin, which indicates the high probability of jump and initial-sensitive excessive motion. As $\tilde{\tau }$ increases, they can be controlled. For f = 0.30 (see Figure 10i–l), even though the safe basin is still fractal, its area becomes much larger with the increase of $\tilde{\tau }$. This shows that the possibility of excessive motion is reduced. It follows from Figure 9 that the delayed position feedback can suppress the erosion of the safe basin effectively when G_p is positive.

4.2. Delayed Velocity Feedback

4.2.1. Primary Resonant Response and Stability of Solutions

As in Section 4.1.1, by applying the Method of Multiple Scales, one can get the slow flow equation of the delayed-velocity-feedback controlled system (30) as below:

$$\frac{da}{dT}=-\frac{\mu }{2}a-\frac{3g{a}^3}{8}+\frac{f\sin \psi }{2}-\frac{G_v(1-\cos \tilde{\tau })}{2}a,\hspace{1em}a\frac{d\psi }{dT}=a(\omega -1)+\frac{3b{a}^3}{8}+\frac{f\cos \psi }{2}+\frac{G_v\sin \tilde{\tau }}{2}a.$$

(48)

The amplitude and phases of the system can be obtained via letting the right side of Equation (48) be zero, given by

$$f\sin \psi =\mu a+\frac{3g}{4}{a}^3+G_v(1-\cos \tilde{\tau })a,\hspace{1em}f\cos \psi =(2-2\omega -\frac{3b}{4}{a}^2-G_v\sin \tilde{\tau })a.$$

(49)

Thus, the amplitude of the approximated periodic solution $x_0=a\cos (\omega T-\psi )$ can be solved from the following equation:

$${(\mu +\frac{3g}{4}{a}^2+G_v(1-\cos \tilde{\tau }))}^2{a}^2+{(2-2\omega -\frac{3b}{4}{a}^2-G_v\sin \tilde{\tau })}^2{a}^2=f^2.$$

(50)

According to its characteristic equation, the stability switch will occur if there exists a real number $\tilde{v}$ satisfying

$$\begin{array}{l}v(\mu +G_v(1-\cos \tilde{\tau })+\frac{3g{a}^2}{2})=0,\\ -v^2+\frac{{(\mu +G_v(1-\cos \tilde{\tau }))}^2}{4}+\frac{{(2-2\omega -G_v\sin \tilde{\tau })}^2}{4}+\frac{\mu }{2}{G}_v(1-\cos \tilde{\tau })+\frac{3a^2}{4}(\mu g+b(2\omega -2+G_v\sin \tilde{\tau }))+\frac{27(g^2+b^2)a^4}{64}=0.\end{array}$$

(51)

For $G_v>0$, in Equation (51), one has $\mu +\frac{3g}{2}{a}^2+G_v(1-\cos \tilde{\tau })>0$, $\tilde{v}$ = 0, and

$$\frac{{(\mu +G_v(1-\cos \tilde{\tau }))}^2}{4}+\frac{{(2-2\omega -G_v\sin \tilde{\tau })}^2}{4}+\frac{\mu }{2}{G}_v(1-\cos \tilde{\tau })+\frac{3a^2}{4}(\mu g+b(2\omega -2+G_v\sin \tilde{\tau }))+\frac{27(g^2+b^2)a^4}{64}=0,$$

(52)

implying that saddle-node bifurcation will occur in the delayed-velocity-feedback controlled system. For $G_v>0$, when $\tilde{\tau }$ is long enough to satisfy

$$4{(\mu g+b(2\omega -2)+b{G}_v\sin \tilde{\tau })}^2-3(g^2+b^2)({(\mu +G_v)}^2+2\mu G_v+G_v{}^2+4{(1-\omega )}^2-2G_v(2\mu +G_v)\cos \tilde{\tau }+4(\omega -1)G_v\sin \tilde{\tau })<0,$$

(53)

there will be no positive roots for Equation (52), meaning no saddle-node bifurcation or bi-stability. For example, given G_v = 0.3 and ω = 0.8, one can calculate from Equation (53) that saddle-node bifurcation will disappear for $\tilde{\tau }>0.68$. Similarly to Section 4.1.1, although delayed velocity feedback can be used to reduce the occurrence of saddle-node bifurcation, it is hard to express the range of $\tilde{\tau }$ for SN bifurcation as a function of f. Fixing the value of f and combining Equations (50), (52) and (53), one can obtain the condition for $\tilde{\tau }$ to reduce saddle-node bifurcation. For instance, given f = 0.13, one can calculate that the condition of $\tilde{\tau }$ to reduce saddle-node bifurcation is $\tilde{\tau }>$ 0.77.

The similar as in Section 4.1.1, for G_v = 0.3 and different values of time delay in the velocity feedback, Figure 11 shows the variation of the response amplitude with the excitation. It follows from the theoretical and numerical results depicted in Figure 11 that, even though the response amplitude is not necessary to be suppressed with the increase of time delay, the velocity feedback can reduce saddle-node bifurcation, and thus bi-stability and jump. The controlling effect of velocity feedback on jump can also be verified by the evolution of basins of attraction with the increase of time delay for f = 0.13, ω = 0.78 (see Figure 12). To be different from the case of the delayed position feedback, as the delay increases, the basin of attraction of the higher-amplitude attractor is expanded, whose boundary becomes smoother and smoother. Finally, it becomes the only bounded attractor (see Figure 12d). For ω = 0.80, if can be concluded from Figure 11b that the delayed velocity feedback can suppress bi-stability. For the fixed value of ω, sequences of the basin of attraction with the increase of f and $\tilde{\tau }$ will be discussed in detail in the next subsection.

4.2.2. Heteroclinic Bifurcation under Delayed Velocity Feedback

As in Section 4.1.2, to apply the Melnikov method in the delayed-velocity-feedback controlled system (30), we also need to treat the delayed feedback as a perturbed term of this system. To this end, the delay $\tilde{\tau }$ should not exceed the first value of stability switch of S₁(0,0) in the linear delayed system whose characteristic equation can expressed as

$${\lambda }^2+\mu \lambda +1+G_v\lambda -G_v\lambda e^{-\lambda \tilde{\tau }}=0$$

(54)

.

The first value of time delay for stability switch of S₁(0,0) will occur when $\lambda =\pm {\Omega }_{v}i$. Substituting it into Equation (54) and separating the imaginary and the real parts of Equation (54) yields

$$-{\Omega }_v{}^2+1=G_v{\Omega }_v\sin ({\Omega }_v\tilde{\tau }),\hspace{1em}\mu {\Omega }_v=G_v{\Omega }_v(\cos ({\Omega }_v\tilde{\tau })-1).$$

(55)

By eliminating the trigonometric functions, the equation above becomes

$${({\Omega }_v^2-1)}^2+({\mu }^2+2\mu G_v){\Omega }_v{}^2=0.$$

(56)

Since $\mu >0$, for $G_v>0$, there will be no real roots of Equation (56), meaning that the stability of the origin will not be changed. Then the delay velocity feedback of the delayed system (30) can be regarded as a disturbed term. The Melnikov function of the controlled system (30) can be expressed as

$$\begin{array}{l}{M}_{\pm }(t_0)={{\displaystyle \int }}_{-\infty }^{+\infty }(-\mu y(T)-g{y}^3(T)+f\cos (\omega (T+t_0))-G_{v}y(T)+G_{v}y(T-\tilde{\tau }))y(T)dT\\ =-\frac{2\sqrt{2}\mu }{3b}-\frac{8\sqrt{2}g}{35b^2}+\sqrt{\frac{2}{b}}\pi \omega f\mathrm{csch}(\frac{\sqrt{2}\pi \omega }{2})\cos (\omega t_0)-\frac{2G_v}{b}{l}_2,\end{array}$$

(57)

where $l_2(\tilde{\tau })={\mathrm{csch}}^2(\frac{\sqrt{2}\tilde{\tau }}{2})(\sqrt{2}-\tilde{\tau }\coth (\frac{\sqrt{2}}{2}\tilde{\tau }))>0$. Accordingly, for

$$f>f^V(\tilde{\tau })=f_0+\frac{\sqrt{2b}{G}_v}{b\pi \omega }\sinh (\frac{\sqrt{2}\pi \omega }{2})l_2(\tilde{\tau }),$$

(58)

And heteroclinic bifurcation may occur. In Equation (58), for $G_v$ > 0, the threshold of heteroclinic bifurcation $f^V(\tilde{\tau })$ increases monopoly with the delay $\tilde{\tau }$. Given $G_v$ = 0.3, the variation of $f^V(\tilde{\tau })$ with the increase of $\tilde{\tau }$ is presented in Figure 13. Numerical and analytical results both illustrate that, for $G_v$ > 0, the delayed velocity feedback can be used to increase the threshold of the amplitude f for heteroclinic bifurcation, thus being effective in reducing the probability of initial-sensitive excessive motion.

Figure 14 shows the evolution of the safe basin of the controlled system (30) with the increase of f and $\tilde{\tau }$. For $\tilde{\tau }$ = 0, Figure 14a,e,i are the same as Figure 10a,e,i, depicting the safe basin of the uncontrolled system (6). Under a positive gain G_v = 0.3, when $\tilde{\tau }$ increases, the area of safe basin of the delayed-velocity-feedback-controlled system will expand obviously. To be different from the results of Section 4.1.2, for f = 0.11 and 0.15 where the two periodic attractors coexisting in the uncontrolled system, as $\tilde{\tau }$ increases, the basin of attraction of the higher-amplitude attractor is expanded and finally becomes the safe basin (see the blue regions in Figure 14a–h). This illustrates that, although the delayed velocity feedback does not necessary reduce the response amplitude, it reduces jump and excessive motion.

5. Discussion

In this paper, a typical relative rotation system is considered, and the phenomena of jump and initial-sensitive excessive motion as well as their suppression via delayed feedback is investigated. The Method of Multiple Scales and the Melnikov method are applied to analyze the conditions for bi-stability and initial-sensitive excessive motion, respectively. By introducing the basins of attraction and safe basin to describe the extent of jump and initial-sensitive excessive motion, respectively, the point mapping approach is employed to present the numerical results which matches the theoretical approach, verifying the validity of the analysis. Some significant results are presented as follows.

(1)

The variation of excitation may induce the coexistence of bistable periodic attractors, which can be ascribed to saddle-node bifurcation.

(2)

The increase of the excitation amplitude may cause initial-sensitive excessive motion, which can be due to heteroclinic bifurcation.

(3)

Under positive coefficients of the feedback gain, the delayed position feedback and the delayed velocity feedback can reduce saddle-node bifurcation and heteroclinic bifurcation so as to suppress jump and initial-sensitive excessive motion. Comparatively, the former can also reduce the amplitude of the response, while the latter may not; the former works well if time delay does not exceed the first stability switch of the trivial equilibrium, while the latter does not have that restriction.

This work presents a detailed analysis of jump and initial-sensitive excessive motion of a typically relative rotation system, which may be beneficial for the performance improvement of rotors and main bearings. The relevant experimental investigations will be included in our future work.