Dynamic modeling of a non-ideal gyroscopic rotor system with nonlinear damping and nonlinear rigidity of an elastic support

Abstract

The article considers the effect of joint linear and nonlinear cubic damping on dynamics of a gyroscopic rigid rotor interacting with an electric motor with a rectilinear characteristic, taking into account the nonlinear rigidity of the support material. The method of regulating the control parameter (voltage on the motor), the amplitude of vibration, and the angular velocity of the shaft in the frequency equation, depending on the value of the coefficient of nonlinear cubic damping of the support, offers the most effective options for controlling resonant oscillations of large amplitudes. It is shown that the greater the value of the coefficient of nonlinear cubic damping, the easier it is to control these oscillations. Moreover, it is proved that the Sommerfeld effect (of the first kind) can also be weakened and eliminate with the help of joint linear and nonlinear damping. To do this, in the case of a rigid characteristic of the nonlinear elasticity of the support material, in a rotor system with a nonideal energy source to eliminate the bistability region, that is, jumping effects, more nonlinear damping of support or energy from a nonideal energy source will be required than in the case of an ideal rotor system.

Introduction

The functioning of many important and responsible vibration machines is based on the use of energy sources of various types and characteristics. There are cases when the power consumed by the oscillatory system is comparable in magnitude, with the power of the exciter-energy source called “nonideal.” The movement of the oscillatory system under the action of such sources is accompanied by the mutual influence of the energy source and the oscillatory system on each other. The interaction of an oscillatory system with an energy source should manifest itself in both stationary and nonstationary modes of motion, especially in such practically important cases as the case of the passage of oscillatory systems through resonances.

The process of passing the rotor system through critical speeds, considering its interaction with the energy source, was considered in many studies.^1–11 In particular, it was found that under conditions when the power reserve of the energy source is small, the course of the process significantly depends on the characteristics of the energy source. Changes in the rotor speed turn out to be closely related to the change in vibration amplitude. As the rotor approaches the lower limit of any unstable speed range, the rotor speed is locked at that speed with increasing swirl amplitudes and the rotor speed does not respond to increased motor power unless there is enough excess power to accelerate the rotor through the corresponding unstable speed range. With such excess power, the transition from a lower limiting rate of instability to a much higher stable rate occurs as a nonlinear surge phenomenon.^1–4 These symptoms, commonly referred to as the Sommerfeld effect, arise from the internal energy coupling between the drive and the driven system and are important design considerations in the design of various rotating machines with flexible supports (bearings).

Dasgupta and Rajan ⁵ studied the Sommerfeld effect using the steady-state amplitude obtained by the instantaneous power balance method and additionally verified by numerical simulation. It also shows the transient characteristics of the dimensionless amplitude and the speed of the shaft with time change through the resonance of the first mode. A rotor system with asymmetric bending flexibility and support bearings has many stable and unstable speed ranges. In contrast to the usual Sommerfeld effect (of the first kind), where the power deficit at resonance is the cause of velocity capture, ideally, there is no need for a residual rotor unbalance in the Sommerfeld effect of the second kind. In both cases of speed capture, excess engine power is spent on increasing swirl amplitudes. The Sommerfeld effect of the first kind refers to the resonance at a synchronous rotor vortex (critical speeds), while the second kind effect refers to the instability of the rotor vortex. ⁴

Due to the anisotropy in the rigid shaft supports, both forward and reverse vortex motions of the rotor are excited with an unbalanced disk with a nonideal energy source – a direct current motor (DC). ⁶ This causes the Sommerfeld multieffect, when two nonlinear surges appear at the first forward and reverse critical speeds during the rotor up and down coast. When the two critical speeds are close to each other, the complex phenomenon of speed lock and surge occurs. Accordingly, transitions through both resonant states are studied. First, analytical solutions are obtained under the assumption of a stationary state; then, transient simulations are performed using a link graph model (BG).

In study of Cveticanin et al., ⁷ two particular cases of motion depending on the frequency properties of the system were assessed. When the frequency properties in both orthogonal directions are equal, only one resonance occurs. If the frequency in one direction is twice as high as in the other, two different resonances occur: one in the x-direction and the other in the y-direction. The conditions for the surge phenomena and the Sommerfeld effect are presented.

The decomposition of the equations of motion, proposed in the article of Awrejcewicz et al., ⁸ allows separating the vibration of the rotor from its rotations. The presented approach can be used to separate vibration from rotation in many other mechanical and mechatronic systems. The behavior of the considered nonideal system near two simultaneously occurring resonances is studied using the Krylov–Bogolyubov averaging method. An analysis of the stability of the resonant response is also carried out.

Kafi and Hosseini ⁹ are studying nonlinear nonstationary vibrations of a rotating composite shaft during the passage of a critical velocity, excited by a nonideal energy source. Geometric nonlinearity, gyroscopic effect, rotational inertia, and coupling due to material anisotropy are taken into account, while shear deformation is not taken into account. Using the Hamilton principle, the equations of motion of axial, bending, torsional, and rotational motion for a composite shaft with a variable speed of rotation are obtained. The influence of external damping, eccentricity, fiber angle, the moment of resistance, and nonlinear terms and the influence of tensile-torsional coupling on the occurrence of the Sommerfeld effect are studied. It is shown that nonlinearity and coupling have a significant effect on the prediction of the Sommerfeld effect. Another article Bisoi et al. ¹⁰ studied the Sommerfeld effect in a strongly gyroscopic rotary dynamic system. The dynamics of an overhanging rotor system near the Sommerfeld effect modes is studied using a discrete and continuous shaft-rotor model in combination with a nonideal motor drive model.

A rigid rotor system driven through a universal joint (U-joint) and mounted on weakly damped flexible supports demonstrates parametric instability of the main and combined types at certain combinations of torque and speed. ¹¹ It is shown that the dynamic coupling between a DC motor and a rotating rotor in a nonideal system can stabilize the system in speed ranges where the corresponding ideal system is unstable. This stabilization depends on constant torque/speed; that is, a motor with a higher constant characteristic tends to stabilize the system better. Input-limited and output-limited stabilization occur through repeated transitions between stable and unstable velocity regions.

Thus, if there is not enough power to overcome the resonance, then the rotor speed may stop at this resonance, or it will take a long time (or additional power) to get out of resonance; thereby damaging the system: causing failure of the rotor shaft, bearings, and other parts of the structure. Therefore, in order to stabilize the movement, it is very important to control the exit from the instability area with the Sommerfeld effect and weaken the effect with the help of vibration isolation.

Moreover, Cveticanin et al. ¹² have researched the dynamics of a nonideal mechanical system that contains an engine, a nonideal source of energy, and an oscillator with a mass slowly changing in time. The system is modeled by two coupled second-order equations with time-varying parameters, in which the engine torque is assumed to be a linear function of the angular velocity. Based on the approximate averaging equations, the amplitude-frequency ratios are determined. The results of numerical solutions show that as the mass increases or decreases, the number of almost stationary positions changes. Based on the results obtained, it is proposed to develop a method for controlling motion in an oscillatory system with a nonideal mass. We would like to highlight the interesting work of Awrejcewicz and Dzyubak,^13,14 in which the nonlinear dynamics of a 2-DOF rotor supported by a magnetohydrodynamic bearing was investigated using perturbation analysis. Non-resonant and various resonant cases (primary resonance with and without internal resonance) are considered. In the case of hard magnetic materials, in order to characterize the closed-loop hysteresis motion control of the rotor, the regions of instability, as well as the contours of the amplitude level of the vertical and horizontal oscillations of the rotor, are obtained in a form suitable for technical applications, that is, in the parametric plane “frequency-amplitude.”

Sghaier et al. ¹⁵ proposed a new dynamic model for unbalanced high-speed rotors with less stringent assumptions. The finite element model takes into account bending, torsion, and tension-compression, resulting in six degrees of freedom at each node. A nonideal energy source is considered and the rotor operates under nonstationary operating conditions and crosses supercritical speeds. The angular displacement is defined in such a way that it combines both its own nominal rotation and torsional deformation. The ability of the proposed model to take into account the mutual influence of transverse and torsional behavior is emphasized and demonstrated using time-frequency analysis.

The Sommerfeld effect causes instability in high-speed rotors. In the work of Goncalves et al., ¹⁶ it was found that the Sommerfeld effect depends on some system parameters and engine operation procedures. These parameters are investigated in order to avoid resonant trapping in the Sommerfeld effect. Thus, reducing this effect is extremely important for the smooth operation of the rotors at high speeds. The purpose of the article of Jha and Dasgupta ¹⁷ is to reduce the Sommerfeld effect of a system with internal damping of an unbalanced flexible shaft and disk using linearized active magnetic bearings. The characteristic equation of the fifth-order polynomial in rotor speed is obtained by balancing the energy of the input power and the mechanical power dissipated in steady state. Using MATLAB simulations, frequency responses are obtained close to system resonance for several active magnetic bearing bias currents. Thus, the Sommerfeld effect is weakened as the bias current is gradually increased. The complete disappearance of the Sommerfeld effect is also reported when the bias current reaches a certain value under certain conditions.

Recently, there has been a growing interest in nonlinear approaches in the design of vibration isolation systems.

The Sommerfeld effect can also be weakened using linear damping and nonlinear cubic damping of viscoelastic flexible rubber supports. ¹⁸

Iskakov ¹⁹ investigated the efficiency of an elastic support with linear damping and nonlinear cubic damping in resonant and nonresonant vibration regions of a rigid gyroscopic rotor with linear stiffness. It is proved that, unlike linear damping, nonlinear damping not only significantly suppresses the maximum resonant vibration amplitude but also preserves vibration isolation in the region beyond the resonant shaft rotation speed. In his next study ²⁰ additionally takes into account the influence of the cubic nonlinear stiffness of the elastic support material on the performance of the insulator. It is shown that nonlinear viscous damping not only reduces the amplitude of high-speed vibration but also eliminates the surge effect. It turns out that if the linear damping narrows the instability area only near the resonant velocity, then the nonlinear cubic damping of the support narrows this area from all sides. It also presents the results of experimental studies of the effect of combined linear and nonlinear cubic damping of rubber support on the frequency response of the rotor, which are in good agreement with the results of analytical studies. In another article, ²¹ the effect of nonlinear cubic damping of elastic support on the response of an ideal rigid gyroscopic rotor system during a nonstationary resonant transition is studied and effects similar to those found for stationary oscillations are confirmed. He proposed a methodology is proposed for determining and identifying linear damping and nonlinear cubic damping of elastic support of a rigid gyroscopic rotor, in which surge effects in the frequency characteristics of a weakly nonlinear oscillating system are eliminated. ²² Al-Solihat and Behdinan ²³ numerically investigated the nonlinear dynamic characteristics and characteristics of the transmitted force of a flexible rotor shaft-disk system supported by a suspension system with nonlinear stiffness and damping. The harmonic balance (HB) method, in combination with the continuation circuit, is used to determine the nonlinear frequency response and force transfer curves due to the disk unbalance force. The influence of bearing stiffness nonlinearity on the system dynamics under the support of hard and soft linear bearing stiffness and the influence of bearing damping nonlinearity are studied and compared with the influence of linear damping. The effects of linear and nonlinear vibration damping of the rubber material of the elastic support are experimentally confirmed. Therefore, in and beyond the resonant region, only nonlinear cubic damping can maintain the performance of the vibration isolator.

Review ²⁴ has discussed the applications of nonlinearity in passive vibration control devices to provide an understanding of how nonlinearity is applied and useful in an implemented system.

In work of Joubert et al. ²⁵ the Rayleigh scattering function is generalized so that anisotropic nonlinear damping can be introduced into the equations of motion. Using a mixture of numerical and symbolic analysis of ODEs of movement of a vibrating gyroscope (VG) for nonlinear damping of anisotropic light, we demonstrate (up to an approximate mean) that Bryan’s law is affected by any form of such damping, causing the pattern to shift, reducing the accuracy of VG.

In the article of Huang et al., ²⁶ the study is related to the isolation characteristics of nonlinear velocity-displacement damping (VDD), which has an arbitrary nonnegative velocity and displacement exponent. Stability analysis is performed for both integer power-law damping and rational degrees of nonlinearities less than unity. The proposed damping can not only suppress the response at resonance but also improve isolation at high frequencies, if the speed and displacement parameters, as well as the damping factor, satisfy certain conditions. The VDD condition is obtained, which does not depend on the excitation amplitude. The measurement results from experimental studies are in good agreement with theoretical predictions, thereby confirming the effects of VDD during base excitation.

Lv and Yao ²⁷ investigated a vibration isolator of systems with one degree of freedom, having nonlinear viscous damping, under force excitation. The stability of the steady-state periodic response was discussed. The effect of damping coefficients on the force transmission and displacement transmission has been studied. The results show that the performance of the nonlinear isolator has some positive effects compared to the linear isolator within a certain range.

Dong et al. ²⁸ suggested in order to improve the efficiency of low-frequency vibration isolation of an insulator with high static and low dynamic stiffness (HSLDS), a new geometric nonlinear damping (GND) design containing semiactive electromagnetic shunt damping. GND is dependent on displacement and vibration speed, which can cause the HSLDS isolator to have different damping characteristics in different frequency ranges. The results show that the HSLDS isolator assembled with GND can meet the requirements of the isolation system under both base and force excitation, providing broadband vibration isolation and low resonant peak while maintaining high-frequency damping. Moreover, GND is superior to linear damping regardless of whether base excitation or force excitation is applied. For base excitation, GND exhibits some desirable properties that cubic nonlinear damping does not have at high frequencies.

For the first time, nonlinear damping is accurately derived from a fractional viscoelastic standard solid model by introducing a geometric nonlinearity into it. ²⁹ The resulting damping model is nonlinear and its frequency response can be tuned with a fractional derivative according to the behavior of the material. Experiments show a strong increase in damping with the oscillation amplitude for nonlinear oscillations of plates and shells. Of the two different continuous structural elements, the free-edge plate is interesting in that it represents the case where energy does not escape through the boundary. ³⁰

The method proposed Lisitano and Bonisoli ³¹ is theoretically derived to identify nonlinear damping forces depending on the degrees of displacement and velocity and is applied to identify linear and nonlinear damping matrices of a system with several degrees of freedom with a localized nonlinear magnetic damper. The coefficients of the nonlinear magnetic damping force are determined for two configurations of the magnetic damper.

The nonlinear vibration isolation system is promising for providing high-performance broadband isolation. Zhang et al. ³² created a generalized vibration isolation system with nonlinear stiffness, nonlinear viscous damping, and Bouc–Wen (BW) hysteresis damping. To evaluate the damping effect, a generalized equivalent damping coefficient with varying stiffness characteristics is determined. It has been found that the damping factor of the linear damping is related to the stiffness change characteristics, while the damping factors of the two types of nonlinear damping are related to the response amplitudes. Linear damping, hysteresis damping, and nonlinear viscous damping are suitable for small amplitude, medium amplitude, and large amplitude conditions, respectively. Hysteresis damping has the added benefit of broadband isolation.

In the article of Balasubramanian et al., ³³  three different models with the same degree of freedom were compared with the same experimental data; each model has a different damping description. In particular, the models are based on modified Duffing generators with linear, quadratic, and cubic stiffness: (i) linear viscous damping; (ii) nonlinear viscoelastic dissipation described by a loss factor; (iii) standard linear solid viscoelastic model with nonlinear springs. The dissipation revealed by various models is discussed, and the main nonlinear character of damping depending on the oscillation amplitude is confirmed.

The single-degree-of-freedom linear oscillator (LO) under the influence of harmonic action is connected to a light mount, which acts, in fact, as a nonlinear energy sink (NES). ³⁴ The strong rigidity and nonlinearity of damping in this system are provided by connecting LO and NES with two inclined linear spring-damping elements. The slope of the connecting elements during movement introduces highly nonlinear geometric effects into the forced dynamics. The analysis results show that increasing the initial tilt angle makes it possible to reduce and even completely eliminate unwanted high amplitude steady-state LO responses that coexist with desirable low amplitude-induced responses in certain frequency ranges. The presented results and associated analytical modeling can be used to design and improve the performance of nonlinear vibration absorbers as vibration reduction devices.

Mofidian and Bardaweel ³⁵ studied the influence of nonlinear cubic viscous damping in a vibration isolation system consisting of a magnetic spring with a positive nonlinear stiffness and a mechanically inclined spring with a geometric nonlinear negative stiffness is studied. The results show that the introduction of nonlinear damping improves the performance of the vibration isolation system, eliminating the undesirable frequency surge phenomena traditionally encountered in vibration isolation systems with quasizero stiffness. Compared to a competing linear vibration isolation system, the described nonlinear system transmits fewer vibrations around the resonant peak.

In work of Lu et al., ³⁶  a method was proposed based on measuring both the frequency of surges and the amplitudes of surges of a system subject to sinusoidal excitations, in terms of frequency and amplitude, respectively. The experimental data give frequency response curves at a fixed excitation amplitude and amplitude response curves at a fixed frequency. Based on the measured surge frequencies and displacement surge amplitudes, the system parameters can be determined by the harmonic balance method. The evaluation results show that the proposed method allows one to estimate the stiffness and damping parameters of a system with strong nonlinearities. Therefore, using knowledge about the phenomena of a nonlinear surge is a promising approach to parameter estimation.

In the article of Luongo et al. ³⁷ the influence of nonlinear damping on the postcritical behavior of the Ziegler column is discussed. The classical double Ziegler pendulum is considered in the mode of finite displacements, in which, in addition, nonlinear damping of the van der Pol type is introduced on the hinges. The nature of the Hopf bifurcation, namely, supercritical or subcritical, and the occurrence of the phenomenon of “hard loss of stability” are investigated. Finally, the influence of nonlinear damping on the limit cycle amplitude for various linearly damped columns is studied. In study of Luongo et al., ³⁸ the effect of linear damping on the postcritical behavior of a Ziegler column is discussed. For this purpose, the well-known double pendulum loaded at the free end with a repeating force, first introduced by Ziegler, is considered in the regime of finite displacements. It is shown that the resulting bifurcation equations are useful for providing qualitative information about the nonlinear mechanical response of the column in the entire damping plane. The verification of the asymptotic solution, carried out using numerical analysis of the exact equations of motion, indicates the effectiveness of the proposed analysis also from the quantitative point of view.

The above review clearly shows the insufficient study of the effect of joint linear and nonlinear cubic damping of support on the dynamics of a rotor system with a nonideal excitation source in stationary and nonstationary modes of motion.

Iskakov et al. ³⁹ considered the effect of joint linear and nonlinear cubic damping on the dynamics of a gyroscopic rigid unbalanced rotor with a nonlinear rigidity of elastic support and a nonideal energy source is considered. It is proven that the combined linear and nonlinear cubic damping of elastic support can more effectively suppress the maximum amplitude and amplitude after resonant oscillations, eliminate abrupt phenomena, significantly reduce the amplitude of the oscillation frequency, and significantly narrow the unstable area until it is completely eliminated.

The purposes of this work, in contrast to the previous, are the effective control of the exit from the region of instability with the Sommerfeld effect, the weakening and the elimination of the Sommerfeld effect using nonlinear cubic support damping and the study of the effect of nonlinear cubic damping on the response dynamics of a gyroscopic rigid rotor with nonlinear rigidity of the support when interacting with a nonideal energy source with a rectilinear characteristic, both stationary and nonstationary modes of the resonant transition, including in the case of an anisotropic linear stiffness of the elastic support.

Theory/calculation

Equations of motion

A rotor system is considered, consisting of a disk and a shaft (Figure 1). A disk with mass $m$ , moment inertia about its axis of rotation $J_P$ , and the moment of inertia about its diametrical axis $J_T$ , fixed without distortion at the free end of a rigid shaft. The shaft with a length $L$ , was installed vertically with the help of the lower hinge and spaced from it at a distance $l_0$ upper elastic support. With this arrangement of the shaft with the disk relative to the supports and at a sufficiently high speed of rotation of the shaft $\stackrel{·}{\phi }$ , the rotor can be considered as a gyroscope. The position of the geometric center of the disk $S$ with respect to the fixed coordinate system $\mathrm{Oxyz}$ is determined by coordinates $x,y,z$ , the position of the shaft in space by the Euler angles $\alpha$ , $\beta$ , and angle of rotation $\phi$ . One of the coordinates will be taken equal and excluded from consideration $z=L=\mathrm{const}$ . Next, we denote the coordinates of the center of mass $m$ driven through $x_m$ and $y_m$ . We also assume that the linear eccentricity $e$ lies along the N - axis of the ONKZ coordinate system. We restrict ourselves to small deviations of the rotor axis. Then $\sin \alpha \approx \alpha ,\sin \beta \approx \beta ,\cos \alpha \approx 1,\cos \beta \approx 1$ .

Figure 1.

Structural diagram of the rotor.

Under these assumptions, the projections of the angular velocity on the coordinate axes of the ONKZ system can be written as follows:

$${\omega }_N\approx -\stackrel{·}{\beta },{\omega }_K\approx \stackrel{·}{\alpha },{\omega }_Z\approx \stackrel{·}{\phi }+\stackrel{·}{\alpha }·\beta .$$

(1)

Relationship between coordinates $x$ and $x_m$ , $y$ and $y_m$ is expressed by the following relations:

$$\begin{array}{c}{x}_m=x+e\cos \phi =L\alpha +e\cos \phi ,\hfill \\ y_m=y+e\sin \phi =L\beta +e\sin \phi .\hfill \end{array}$$

(2)

Thus, the position of the rigid shaft and the rotor as a whole in space is determined by the angles $\alpha ,\beta$ , and $\phi$ .

The equations of motion of the rotor system are derived by the Lagrange method (Appendix A) from the expressions for kinetic and potential energy, the Rayleigh function, the moments of gravity, the force of inertia of the mass imbalance, and the dynamic moment of the engine:

$$\begin{array}{c}(J_T+m{L}^2)\stackrel{··}{\alpha }+J_P\stackrel{·}{\phi }\stackrel{·}{\beta }+{\mu }_{d1}\stackrel{·}{\alpha }+{\mu }_{d3}{\stackrel{·}{\alpha }}^3+(c_1{l}_0^2-\mathrm{mgL})\alpha \hfill \\ \multicolumn{1}{c}{+c_3{l}_0^4{\alpha }^3=\mathrm{me}{\stackrel{·}{\phi }}^{2}L\cos \phi ,}\\ (J_T+m{L}^2)\stackrel{··}{\beta }-J_P\stackrel{·}{\phi }\stackrel{·}{\alpha }+{\mu }_{d1}\stackrel{·}{\beta }+{\mu }_{d3}{\stackrel{·}{\beta }}^3+(c_1{l}_0^2-\mathrm{mgL})\beta \hfill \\ \multicolumn{1}{c}{+c_3{l}_0^4{\beta }^3=\mathrm{me}{\stackrel{·}{\phi }}^{2}L\sin \phi ,}\\ J_P\stackrel{··}{\phi }=\mathrm{meL}(\stackrel{··}{\alpha }\sin \phi -\stackrel{··}{\beta }\cos \phi )-J_P(\stackrel{··}{\alpha }\beta +\stackrel{·}{\alpha }\stackrel{·}{\beta })\hfill \\ \multicolumn{1}{c}{+C_M\mathrm{\Phi }U/R-(C_M{C}_E{\mathrm{\Phi }}^2/R+q_m)\stackrel{·}{\phi }.}\end{array}$$

(3)

Let us introduce the following dimensionless parameters:

$$\begin{array}{c}l=l_0/L;\overline{t}=t{\omega }_0;{\overline{J}}_p=J_p/(m{L}^2);{\overline{J}}_T=J_T/(m{L}^2);\hfill \\ \multicolumn{1}{c}{\overline{G}=g/\left(L{\omega }_0^2\right);{\overline{C}}_1=c_1/\left(m{\omega }_0^2\right);e_r=e/\left[L(1+{\overline{J}}_T)\right];}\\ C_3=c_3{l}_0^4/\left[m{L}^2{\omega }_0^2(1+{\overline{J}}_T)\right];{\mu }_1={\mu }_{d1}/\left[m{L}^2{\omega }_0(1+{\overline{J}}_T)\right];\hfill \\ \multicolumn{1}{c}{{\mu }_3={\mu }_{d3}{\omega }_0/\left[m{L}^2(1+{\overline{J}}_T)\right];}\\ u_1=C_M\mathrm{\Phi }{\omega }_{0}U/\left[\mathrm{Rm}{L}^2(1+{\overline{J}}_T)\right];\hfill \\ \multicolumn{1}{c}{u_2=(C_M{C}_E{\mathrm{\Phi }}^2/R+q_m)/\left[m{L}^2{\omega }_0(1+{\overline{J}}_T)\right];}\\ \multicolumn{1}{c}{J_{p1}={\overline{J}}_p/(1+{\overline{J}}_T),}\end{array}$$

(4)

where

$${\omega }_0=\sqrt{\frac{c_1{l}_0^2-\mathrm{mgL}}{m{L}^2-(J_p-J_T)}}$$

(5)

is the natural frequency of the rotor system (3). Using (4), we give the equations of motion (3) a compact dimensionless form:

$$\begin{array}{c}\alpha ″+J_{p1}\phi \prime \beta \prime +{\mu }_1\alpha \prime +{\mu }_3{\alpha \prime }^3+{\omega }_n^2\alpha \hfill \\ \multicolumn{1}{c}{+C_3{\alpha }^3=e_r{\phi \prime }^2\cos \phi ,}\\ \beta ″-J_{p1}\phi \prime \alpha \prime +{\mu }_1\beta \prime +{\mu }_3{\beta \prime }^3+{\omega }_n^2\beta +C_3{\beta }^3\hfill \\ \multicolumn{1}{c}{=e_r{\phi \prime }^2\sin \phi ,}\\ \phi ″=[e_r(\alpha ″\sin \phi -\beta ″\cos \phi )-J_{p1}(\alpha ″\beta +\alpha \prime \beta \prime )+u_1-u_2\phi \prime ]/J_{p1},\hfill \end{array}$$

(6)

where

$${\omega }_n=\sqrt{\frac{{\overline{C}}_1{l}^2-\overline{G}}{1+{\overline{J}}_T}}$$

(7)

is the dimensionless natural frequency of the rotor system (6), where it is accepted ${\overline{J}}_T>>{\overline{J}}_P$ . In the third equation of equations (6),

$$u_1-u_2\phi \prime =M$$

(8)

is the dimensionless dynamic torque of the engine. Here, $u_1$ is the control parameter, depending on the voltage on the motor; $u_2$ is the parameter depending on the type of energy source. ⁴⁰

Here, the prime denotes the derivative with respect to the dimensionless time $\overline{t}$ .

We are interested in the influence of nonlinear cubic damping on the dynamics of a weakly nonlinear rotor system, taking into account the nonlinear nature of the support rigidity and the change in voltage on the motor with a characteristic taken as a straight line. ⁴⁰

For the analytical solution of the equations of motion of the rotor by the Bogolyubov method,^40–42 the following restrictions are accepted. Projections of moments of damping forces ${\mu }_1\alpha \prime ,{\mu }_1\beta \prime$ and ${\mu }_3\alpha {\prime }^3,{\mu }_3\beta {\prime }^3$ and the moment of the nonlinear restoring force $C_3{\alpha }^3,C_3{\beta }^3$ , moments of inertia forces $e_r\phi {\prime }^2\cos \phi ,e_r\phi {\prime }^2\sin \phi ,$ $e_r\alpha ″\sin \phi ,-e_r\beta ″\cos \phi$ , and passive gyroscopic moments $-J_{p1}\phi \prime \alpha \prime$ , $J_{p1}\phi \prime \beta \prime$ , $-J_{P1}(\alpha ″\beta +\alpha \prime \beta \prime )$ are assumed to be small compared to other moments of forces acting in the system. We confine ourselves to consideration of regimes close to stationary motions; therefore, the acceleration $\phi ″$ will be a small quantity, also considering motion in the resonance region, where the frequency of free oscillations ${\omega }_n$ is close to the frequency of forced oscillations $\mathrm{\Xi }$ ; that is, ${\omega }_n-\mathrm{\Xi }=\epsilon \tilde{\sigma },\epsilon <<1.$

After introducing a small parameter $\epsilon$ in accordance with the accepted restrictions of the problem and the new notation for quantities,

$$\begin{array}{c}{\mu }_1=\epsilon {\tilde{\mu }}_1,{\mu }_3=\epsilon {\tilde{\mu }}_3,C_3=\epsilon {\tilde{C}}_3,J_{P1}=\epsilon {\tilde{J}}_P,\hfill \\ e_r=\epsilon \tilde{e},u_1=\epsilon {\tilde{u}}_1,u_2=\epsilon {\tilde{u}}_2.\hfill \end{array}$$

(9)

The system of equations (6) will take the following form:

$$\begin{array}{c}\alpha ″+{\omega }_n^2\alpha =\epsilon (\tilde{e}{\phi \prime }^2\cos \phi -{\tilde{J}}_P\phi \prime \beta \prime -{\tilde{\mu }}_1\alpha \prime -{\tilde{\mu }}_3{\alpha \prime }^3-{\tilde{C}}_3{\alpha }^3),\hfill \\ \beta ″+{\omega }_n^2\beta =\epsilon (\tilde{e}{\phi \prime }^2\sin \phi +{\tilde{J}}_P\phi \prime \alpha \prime -{\tilde{\mu }}_1\beta \prime -{\tilde{\mu }}_3{\beta \prime }^3-{\tilde{C}}_3{\beta }^3),\hfill \\ \phi ″=\epsilon [\tilde{e}(\alpha ″\sin \phi -\beta ″\cos \phi )-{\tilde{J}}_P(\alpha ″\beta +\alpha \prime \beta \prime )+{\tilde{u}}_1-{\tilde{u}}_2\phi \prime ]/J_{P1}.\hfill \end{array}$$

(10)

Equations (10) are a system of nonlinear ordinary differential equations of the second order with respect to $\alpha ,\beta$ , and $\phi$ with a small parameter $\epsilon$ .

Solutions of the equations of motion

Equations (10) for $\epsilon =0$ describe harmonic oscillations and rotation at a constant frequency $\frac{d\phi }{d\overline{t}}=\mathrm{const}$ , at $\epsilon \ne 0$ , oscillations close to harmonic oscillations with a frequency $\frac{d\phi }{d\overline{t}}$ , close to constant, that is, slowly changing.

Solutions to the system of equations (10) are sought in the following form:

$$\alpha =A\cos (\phi +\gamma ),$$

(11)

$$\beta =A\sin (\phi +\gamma ).$$

(12)

Here, the variables $A,\gamma ,\mathrm{\Xi }$ will be slowly varying functions of time $\overline{t}$ . They represent the most essential movement parameters: $A$ is the oscillation amplitude; $\gamma$ is the phase shift angle between the angular coordinate $\alpha$ or $\beta$ and disturbing moment.

Following Kononenko, ⁴² we accept new variables corresponding to the above representations:

$$\frac{d\alpha }{d\overline{t}}=-A{\omega }_n\sin (\phi +\gamma ),$$

(13)

$$\frac{d\beta }{d\overline{t}}=A{\omega }_n\cos (\phi +\gamma ),$$

(14)

$$\frac{d\phi }{d\overline{t}}=\mathrm{\Xi },$$

(15)

where $\mathrm{\Xi }$ is the frequency of the disturbing moment or the angular speed of rotation of the motor shaft.

Differentiating expression (11) and taking into account expression (15), we have the following:

$$\frac{d\alpha }{d\overline{t}}=\frac{\mathrm{dA}}{d\overline{t}}\cos (\phi +\gamma )-(\mathrm{\Xi }+\frac{d\gamma }{d\overline{t}})A\sin (\phi +\gamma )$$

and equating it with expression (13), we obtain

$$\frac{\mathrm{dA}}{d\overline{t}}\cos (\phi +\gamma )-\frac{d\gamma }{d\overline{t}}A\sin (\phi +\gamma )=(\mathrm{\Xi }-{\omega }_n)A\sin (\phi +\gamma ).$$

(16)

Differentiating expression (13) and taking into account expression (15), we obtain

$$\frac{d^2\alpha }{d{\overline{t}}^2}=-\frac{\mathrm{dA}}{d\overline{t}}{\omega }_n\sin (\phi +\gamma )-(\mathrm{\Xi }+\frac{d\gamma }{d\overline{t}})A{\omega }_n\cos (\phi +\gamma ).$$

(17)

Expression (17) together with (11), (13), and (14) will be substituted into the first and third equations of system (10):

$$\begin{array}{c}\frac{\mathrm{dA}}{d\overline{t}}{\omega }_n\sin (\phi +\gamma )+\frac{d\gamma }{d\overline{t}}A{\omega }_n\cos (\phi +\gamma )\hfill \\ \multicolumn{1}{c}{-A{\omega }_n({\omega }_n-\mathrm{\Xi })\cos (\phi +\gamma )=}\\ -\epsilon \left[\begin{array}{c}\tilde{e}{\mathrm{\Xi }}^2\cos \phi +{\tilde{\mu }}_1{\omega }_{n}A\sin (\phi +\gamma )+{\tilde{\mu }}_3{\left({\omega }_{n}A\right)}^3{\sin }^3(\phi +\gamma )\hfill \\ \multicolumn{1}{c}{-{\tilde{C}}_3{A}^3{\cos }^3(\phi +\gamma )}\end{array}\right],\hfill \\ \frac{d\mathrm{\Xi }}{d\overline{t}}=\frac{1}{J_{P1}}\epsilon \left\{\begin{array}{c}\tilde{e}\left[A{\omega }_n\mathrm{\Xi }(-\sin \phi \cos (\phi +\gamma )+\cos \phi \sin (\phi +\gamma ))\right]\hfill \\ \multicolumn{1}{c}{+{\tilde{J}}_P{A}^2{\omega }_n({\omega }_n+\mathrm{\Xi })\cos (\phi +\gamma )\sin (\phi +\gamma )+{\tilde{u}}_1-{\tilde{u}}_2\phi \prime }\end{array}\right\}.\hfill \end{array}$$

(18)

We solve equations (18) and equation (16) with respect to derivatives $\frac{\mathrm{dA}}{d\overline{t}}$ and $\frac{d\gamma }{d\overline{t}}$ , and we get

$$\begin{array}{c}\frac{\mathrm{dA}}{d\overline{t}}=-\frac{\epsilon }{{\omega }_n}\left[\begin{array}{c}\tilde{e}{\mathrm{\Xi }}^2\cos \phi +{\tilde{\mu }}_1{\omega }_{n}A\sin (\phi +\gamma )\hfill \\ \multicolumn{1}{c}{+{\tilde{\mu }}_3{\left({\omega }_{n}A\right)}^3{\sin }^3(\phi +\gamma )-{\tilde{C}}_3{A}^3{\cos }^3(\phi +\gamma )}\end{array}\right]\hfill \\ \multicolumn{1}{c}{\sin (\phi +\gamma ),\frac{d\gamma }{d\overline{t}}={\omega }_n-\mathrm{\Xi }-\frac{\epsilon }{A{\omega }_n}}\\ \multicolumn{1}{c}{\left[\begin{array}{c}\tilde{e}{\mathrm{\Xi }}^2\cos \phi +{\tilde{\mu }}_1{\omega }_{n}A\sin (\phi +\gamma )\hfill \\ \multicolumn{1}{c}{+{\tilde{\mu }}_3{\left({\omega }_{n}A\right)}^3{\sin }^3(\phi +\gamma )-{\tilde{C}}_3{A}^3{\cos }^3(\phi +\gamma )}\end{array}\right]\cos (\phi +\gamma ).}\end{array}$$

(19)

The system of equations (19) is equivalent to system (10) so that the system of equations (10) is transformed as follows:

$$\begin{array}{l}\frac{d\Xi }{d\overline{t}}=\frac{1}{J_{P1}}\\ \epsilon \left\{\begin{array}{l}{\tilde{u}}_1-{\tilde{u}}_2\Xi +\tilde{e}A{\omega }_n\Xi \left[-\sin \phi \cos \left(\phi +\gamma \right)+\cos \phi \sin \left(\phi +\gamma \right)\right]\\ +{\tilde{J}}_P{A}^2{\omega }_n\left({\omega }_n+\Xi \right)\cos \left(\phi +\gamma \right)\sin \left(\phi +\gamma \right)\end{array}\right\},\\ \frac{dA}{d\overline{t}}=-\frac{\epsilon }{{\omega }_n}\left[\begin{array}{l}\tilde{e}{\Xi }^2\cos \phi +{\tilde{\mu }}_1{\omega }_{n}A\sin \left(\phi +\gamma \right)\\ +{\tilde{\mu }}_3{\left({\omega }_{n}A\right)}^3{\sin }^3\left(\phi +\gamma \right)-{\tilde{C}}_3{A}^3{\cos }^3\left(\phi +\gamma \right)\end{array}\right]\\ \sin \left(\phi +\gamma \right)+{\epsilon }^2\dots ,\\ \frac{d\gamma }{d\overline{t}}=\\ \epsilon \left\{\tilde{\sigma }-\frac{1}{A{\omega }_n}\left[\begin{array}{l}\tilde{e}{\Xi }^2\cos \phi +{\tilde{\mu }}_1{\omega }_{n}A\sin \left(\phi +\gamma \right)\\ +{\tilde{\mu }}_3{\left({\omega }_{n}A\right)}^3{\sin }^3\left(\phi +\gamma \right)-{\tilde{C}}_3{A}^3{\cos }^3\left(\phi +\gamma \right)\end{array}\right]\cos \left(\phi +\gamma \right)\right\}+{\epsilon }^2\dots \end{array}$$

(20)

Following the Bogolyubov method,^40–42 approximate solutions of equations (20) can be sought in the following form:

$$\mathrm{\Xi }=\mathrm{\Omega }+\epsilon E_1(\overline{t},\mathrm{\Omega },a,\chi ),A=a+\epsilon E_2(\overline{t},\mathrm{\Omega },a,\chi ),\gamma =\chi +\epsilon E_3(\overline{t},\mathrm{\Omega },a,\chi ),$$

(21)

where $\epsilon E_1(\overline{t},\mathrm{\Omega },a,\chi ),\epsilon E_2(\overline{t},\mathrm{\Omega },a,\chi ),\epsilon E_3(\overline{t},\mathrm{\Omega },a,\chi )$ is small periodic functions of time $\overline{t}$ .

The equations of the first approximation will be obtained by averaging the right-hand sides of (20) over an explicitly specified time $\overline{t}$ . It is easy to see that the same result will be obtained if the averaging takes place directly over the coordinate $\phi$ because here $\phi$ is the quasicyclic coordinate. Proceeding in this way, we write the equations of the first approximation to determine $\mathrm{\Omega },a,\chi$ in the following form:

$$\begin{array}{c}\frac{d\mathrm{\Omega }}{d\overline{t}}=\hfill \\ \multicolumn{1}{c}{\frac{\epsilon }{2\pi }\underset{0}{\overset{2\pi }{\int }}\frac{1}{J_{P1}}\left[\begin{array}{c}{\tilde{u}}_1-{\tilde{u}}_2\mathrm{\Omega }+\tilde{e}a{\omega }_n\mathrm{\Omega }(-\sin \phi \cos \psi +\cos \phi \sin \psi )\hfill \\ \multicolumn{1}{c}{+{\tilde{J}}_P{a}^2{\omega }_n({\omega }_n+\mathrm{\Omega })\cos (\phi +\chi )\sin (\phi +\chi )}\end{array}\right]d\phi ,}\\ \multicolumn{1}{c}{\psi =\phi +\chi ,}\\ \frac{\mathrm{da}}{d\overline{t}}=-\frac{\epsilon }{2\pi {\omega }_n}\underset{0}{\overset{2\pi }{\int }}\left[\begin{array}{c}\tilde{e}{\mathrm{\Omega }}^2\cos \phi +{\tilde{\mu }}_1{\omega }_{n}a\sin \psi \hfill \\ \multicolumn{1}{c}{+{\tilde{\mu }}_3{\left({\omega }_{n}a\right)}^3{\sin }^3\psi -{\tilde{C}}_3{a}^3{\cos }^3\psi }\end{array}\right]\sin \psi d\phi ,\hfill \\ \frac{d\chi }{d\overline{t}}=\epsilon \{\sigma -\frac{1}{2\pi {\omega }_n}\underset{0}{\overset{2\pi }{\int }}\frac{1}{a}\left[\begin{array}{c}\tilde{e}{\mathrm{\Omega }}^2\cos \phi +{\tilde{\mu }}_1{\omega }_{n}a\sin \psi \hfill \\ \multicolumn{1}{c}{+{\tilde{\mu }}_3{\left({\omega }_{n}a\right)}^3{\sin }^3\psi -{\tilde{C}}_3{a}^3{\cos }^3\psi }\end{array}\right]\cos \psi d\phi \},\hfill \\ \epsilon \sigma ={\omega }_n-\mathrm{\Omega }.\hfill \end{array}$$

(22)

After performing the averaging operation and eliminating the small parameter $\epsilon$ in accordance with the formulas for the transition of values (9), we obtain

$$\begin{array}{c}\frac{d\mathrm{\Omega }}{d\overline{t}}=\frac{1}{J_{P1}}(u_1-u_2\mathrm{\Omega }+e_{r}a{\omega }_n\mathrm{\Omega }\sin \chi ),\hfill \\ \frac{\mathrm{da}}{d\overline{t}}=-\frac{1}{2}(e_r\frac{{\mathrm{\Omega }}^2}{{\omega }_n}\sin \chi +{\mu }_{1}a+\frac{3}{4}{\mu }_3{\omega }_n^2{a}^3),\hfill \\ \frac{d\chi }{d\overline{t}}={\omega }_n-\mathrm{\Omega }-\frac{e_r}{2{\omega }_{n}a}{\mathrm{\Omega }}^2\cos \chi +\frac{3C_3{a}^2}{8{\omega }_n},\hfill \\ \frac{d\phi }{d\overline{t}}=\mathrm{\Omega }.\hfill \end{array}$$

(23)

Equations (23) can also be obtained using relations (12)–(15) and the Bogolyubov method^40–42 as applied to the second equation of system (10).

Under conditions

$$\frac{d\mathrm{\Omega }}{d\overline{t}}=0,\frac{\mathrm{da}}{d\overline{t}}=0,\frac{d\chi }{d\overline{t}}=0$$

from (23), we obtain equations that determine the stationary modes of motion:

$$\begin{array}{c}{u}_1-u_2\mathrm{\Omega }+e_{r}a{\omega }_n\mathrm{\Omega }\sin \chi =0,\hfill \\ e_r{\mathrm{\Omega }}^2\sin \chi /{\omega }_n+{\mu }_{1}a+(3/4){\mu }_3{\omega }_n^2{a}^3=0,\hfill \\ {\omega }_n-\mathrm{\Omega }-e_r{\mathrm{\Omega }}^2\cos \chi /\left(2{\omega }_{n}a\right)+3C_3{a}^2/\left(8{\omega }_n\right)=0.\hfill \end{array}$$

(24)

From the system of equations (24), we find the equation for determining the frequency:

$$u_1-u_2\mathrm{\Omega }-({\mu }_1+0.75{\mu }_3{\omega }_n^2{a}^2){\omega }_n^2{a}^2/\mathrm{\Omega }=0.$$

(25)

The expression for determining the amplitude of oscillations $a$ is as follows:

$$f(\mathrm{\Omega },a)=\{{({\mu }_1{\omega }_n+0.75{\mu }_3{\omega }_n^3{a}^2)}^2+{[0.75C_3{a}^2+2{\omega }_n({\omega }_n-\mathrm{\Omega })]}^2\}{a}^2-{\left(e_r{\mathrm{\Omega }}^2\right)}^2=0.$$

(26)

Expression for determining the phase of oscillations $\chi$ is as follows:

$$\tan \chi =-\frac{{\mu }_1{\omega }_n+0.75{\mu }_3{\omega }_n^3{a}^2}{0.75C_3{a}^2+2{\omega }_n({\omega }_n-\mathrm{\Omega })}.$$

(27)

In the absence of nonlinear damping, formulas (25)–(27) completely coincide with similar results of Kononenko’s work ⁴² for an oscillatory system with a nonlinear restoring force.

Stability criterion

Analysis of the stability of periodic solutions is performed using approximate first-order differential equations (23), which can be written in an abbreviated form:

$$\frac{d\mathrm{\Omega }}{d\overline{t}}=F_1(a,\mathrm{\Omega },\chi ),\frac{\mathrm{da}}{d\overline{t}}=F_2(a,\mathrm{\Omega },\chi ),\frac{d\chi }{d\overline{t}}=F_3(a,\mathrm{\Omega },\chi ).$$

(28)

In the stationary state, equations (28) are equal to zero:

$$F_1(a,\mathrm{\Omega },\chi )=0,F_2(a,\mathrm{\Omega },\chi )=0,F_3(a,\mathrm{\Omega },\chi )=0.$$

Taking into account the difference between the perturbed and unperturbed equations, the differential equations in variations are written as follows:

$$\begin{array}{c}\frac{d{\delta }_{\mathrm{\Omega }}}{d\overline{t}}={\left(\frac{\partial F_1}{\partial \mathrm{\Omega }}\right)}_0{\delta }_{\mathrm{\Omega }}+{\left(\frac{\partial F_1}{\partial a}\right)}_0{\delta }_a+{\left(\frac{\partial F_1}{\partial \chi }\right)}_0{\delta }_{\chi },\hfill \\ \frac{d{\delta }_a}{d\overline{t}}={\left(\frac{\partial F_2}{\partial \mathrm{\Omega }}\right)}_0{\delta }_{\mathrm{\Omega }}+{\left(\frac{\partial F_2}{\partial a}\right)}_0{\delta }_a+{\left(\frac{\partial F_2}{\partial \chi }\right)}_0{\delta }_{\chi },\hfill \\ \frac{d{\delta }_{\chi }}{d\overline{t}}={\left(\frac{\partial F_3}{\partial \mathrm{\Omega }}\right)}_0{\delta }_{\mathrm{\Omega }}+{\left(\frac{\partial F_3}{\partial a}\right)}_0{\delta }_a+{\left(\frac{\partial F_3}{\partial \chi }\right)}_0{\delta }_{\chi }.\hfill \end{array}$$

(29)

The characteristic determinant of equation (29) has the following form:

$$\left|\begin{array}{ccc}{\left(\frac{\partial F_1}{\partial \mathrm{\Omega }}\right)}_0{\delta }_{\mathrm{\Omega }}-\rho \hfill & {\left(\frac{\partial F_1}{\partial a}\right)}_0{\delta }_a\hfill & {\left(\frac{\partial F_1}{\partial \chi }\right)}_0{\delta }_{\chi }\hfill \\ \multicolumn{1}{c}{{\left(\frac{\partial F_2}{\partial \mathrm{\Omega }}\right)}_0{\delta }_{\mathrm{\Omega }}}& {\left(\frac{\partial F_2}{\partial a}\right)}_0{\delta }_a-\rho \hfill & {\left(\frac{\partial F_2}{\partial \chi }\right)}_0{\delta }_{\chi }\hfill \\ \multicolumn{1}{c}{{\left(\frac{\partial F_3}{\partial \mathrm{\Omega }}\right)}_0{\delta }_{\mathrm{\Omega }}}& {\left(\frac{\partial F_3}{\partial a}\right)}_0{\delta }_a\hfill & {\left(\frac{\partial F_3}{\partial \chi }\right)}_0{\delta }_{\chi }-\rho \hfill \end{array}\right|=0.$$

(30)

Index “0” means that partial derivatives of F are functions at the equilibrium point. The stability of approximate solutions (23) depends on the roots of the characteristic equation (30). The solutions are stable if the eigenvalues of the characteristic equation have negative real parts.

In equations (30), the partial derivatives have the following expressions:

$$\begin{array}{c}{\left(\frac{\partial F_1}{\partial \mathrm{\Omega }}\right)}_0=\frac{1}{J_{P1}}(-u_2+e_{r}a{\omega }_n\sin \chi ),\hfill \\ \multicolumn{1}{c}{{\left(\frac{\partial F_1}{\partial a}\right)}_0=\frac{1}{J_{P1}}\left(e_r{\omega }_n\mathrm{\Omega }\sin \chi \right),{\left(\frac{\partial F_1}{\partial \chi }\right)}_0=\frac{1}{J_{P1}}\left(e_{r}a{\omega }_n\mathrm{\Omega }\cos \chi \right),}\\ {\left(\frac{\partial F_2}{\partial \mathrm{\Omega }}\right)}_0=-e_r\frac{\mathrm{\Omega }}{{\omega }_n}\sin \chi ,{\left(\frac{\partial F_2}{\partial a}\right)}_0=-\frac{1}{2}({\mu }_1+\frac{9}{4}{\mu }_3{\omega }_n^2{a}^2),\hfill \\ \multicolumn{1}{c}{{\left(\frac{\partial F_2}{\partial \chi }\right)}_0=-e_r\frac{{\mathrm{\Omega }}^2}{2{\omega }_n}\cos \chi ,}\\ {\left(\frac{\partial F_3}{\partial \mathrm{\Omega }}\right)}_0=-(1+\frac{e_r}{{\omega }_{n}a}\mathrm{\Omega }\cos \chi ),\hfill \\ \multicolumn{1}{c}{{\left(\frac{\partial F_3}{\partial a}\right)}_0=\frac{e_r}{2{\omega }_n{a}^2}{\mathrm{\Omega }}^2\cos \chi +\frac{3C_{3}a}{4{\omega }_n},}\\ \multicolumn{1}{c}{{\left(\frac{\partial F_3}{\partial \chi }\right)}_0=\frac{e_r}{2{\omega }_{n}a}{\mathrm{\Omega }}^2\sin \chi ,}\end{array}$$

(31)

where

$$\sin \chi =-\frac{1}{e_{r}a{\omega }_n\mathrm{\Omega }}(u_1-u_2\mathrm{\Omega }),\cos \chi =\sqrt{1-{\sin }^2\chi }=\sqrt{1-{\left[\frac{1}{e_{r}a{\omega }_n\mathrm{\Omega }}(u_1-u_2\mathrm{\Omega })\right]}^2}.$$

(32)

Equations of motion in the case of anisotropy of support flexibility

In the case of an anisotropic linear stiffness of elastic support in two mutually perpendicular directions, under the assumption that the stiffness and damping nonlinearities are small and isotropic, the nonlinear differential equations of rotor motion are derived by the Lagrange method (Appendix B) and take the following form:

$$\begin{array}{c}(J_T+m{L}^2)\stackrel{··}{\alpha }+J_P\stackrel{·}{\phi }\stackrel{·}{\beta }+{\mu }_{d1}\stackrel{·}{\alpha }+{\mu }_{d3}{\stackrel{·}{\alpha }}^3+\hfill \\ \multicolumn{1}{c}{(c_{11}{l}_0^2-\mathrm{mgL})\alpha +c_3{l}_0^4{\alpha }^3=\mathrm{me}{\stackrel{·}{\phi }}^{2}L\cos \phi ,}\\ (J_T+m{L}^2)\stackrel{··}{\beta }-J_P\stackrel{·}{\phi }\stackrel{·}{\alpha }+{\mu }_{d1}\stackrel{·}{\beta }+{\mu }_{d3}{\stackrel{·}{\beta }}^3+\hfill \\ \multicolumn{1}{c}{(c_{12}{l}_0^2-\mathrm{mgL})\beta +c_3{l}_0^4{\beta }^3=\mathrm{me}{\stackrel{·}{\phi }}^{2}L\sin \phi ,}\\ J_P\stackrel{··}{\phi }=\mathrm{meL}(\stackrel{··}{\alpha }\sin \phi -\stackrel{··}{\beta }\cos \phi )-J_P(\stackrel{··}{\alpha }\beta +\stackrel{·}{\alpha }\stackrel{·}{\beta })+\hfill \\ \multicolumn{1}{c}{C_M\mathrm{\Phi }U/R-(C_M{C}_E{\mathrm{\Phi }}^2/R+q_m)\stackrel{·}{\phi }.}\end{array}$$

(33)

Natural frequencies (critical speeds) of a dampened rotor system (33) are as follows:

$${\omega }_{1,2}=\sqrt{\tilde{b}/2\mp \sqrt{{\tilde{b}}^2/4-\tilde{c},}}$$

(34)

where

$$\tilde{b}=\frac{(J_T+m{L}^2)[(c_{11}{l}_0^2-\mathrm{mgL})+(c_{12}{l}_0^2-\mathrm{mgL})]}{{(J_T+m{L}^2)}^2-J_P^2},\tilde{c}=\frac{(c_{11}{l}_0^2-\mathrm{mgL})(c_{12}{l}_0^2-\mathrm{mgL})}{{(J_T+m{L}^2)}^2-J_P^2},$$

(35)

when $J_T>>$ $J_P$ : $\tilde{b}=(c_{11}{l}_0^2-\mathrm{mgL})/(J_T+m{L}^2)+(c_{12}{l}_0^2-\mathrm{mgL})/(J_T+m{L}^2),$ $\tilde{c}=[(c_{11}{l}_0^2-\mathrm{mgL})/(J_T+m{L}^2)][(c_{12}{l}_0^2-\mathrm{mgL})/(J_T+m{L}^2)].$

Let us introduce the following dimensionless parameters:

$$\begin{array}{c}l=l_0/L;\overline{t}=t{\omega }_1;{\overline{J}}_p=J_p/(m{L}^2);{\overline{J}}_T=J_T/(m{L}^2);\hfill \\ \multicolumn{1}{c}{\overline{G}=g/\left(L{\omega }_1^2\right);{\overline{C}}_{11}=c_{11}/\left(m{\omega }_1^2\right);}\\ {\overline{C}}_{12}=c_{12}/\left(m{\omega }_1^2\right);e_r=e/\left[L(1+{\overline{J}}_T)\right];\hfill \\ \multicolumn{1}{c}{C_3=c_3{l}_0^4/\left[m{L}^2{\omega }_1^2(1+{\overline{J}}_T)\right];}\\ {\mu }_1={\mu }_{d1}/\left[m{L}^2{\omega }_1(1+{\overline{J}}_T)\right];\hfill \\ \multicolumn{1}{c}{{\mu }_3={\mu }_{d3}{\omega }_1/\left[m{L}^2(1+{\overline{J}}_T)\right];J_{p1}={\overline{J}}_p/(1+{\overline{J}}_T);}\\ u_1=C_M\mathrm{\Phi }{\omega }_{1}U/\left[\mathrm{Rm}{L}^2(1+{\overline{J}}_T)\right];\hfill \\ \multicolumn{1}{c}{u_2=(C_M{C}_E{\mathrm{\Phi }}^2/R+q_m)/\left[m{L}^2{\omega }_1(1+{\overline{J}}_T)\right].}\end{array}$$

(36)

Using (36), the equation of motion (33) can be represented in a dimensionless form:

$$\begin{array}{c}\alpha ″+J_{p1}\phi \prime \beta \prime +{\mu }_1\alpha \prime +{\mu }_3{\alpha \prime }^3+{\omega }_{n1}^2\alpha +C_3{\alpha }^3=e_r{\phi \prime }^2\cos \phi ,\hfill \\ \beta ″-J_{p1}\phi \prime \alpha \prime +{\mu }_1\beta \prime +{\mu }_3{\beta \prime }^3+{\omega }_{n2}^2\beta +C_3{\beta }^3=e_r{\phi \prime }^2\sin \phi ,\hfill \\ \phi ″=\frac{1}{J_{p1}}[e_r(\alpha ″\sin \phi -\beta ″\cos \phi )-J_{p1}(\alpha ″\beta +\alpha \prime \beta \prime )+u_1-u_2\phi \prime ],\hfill \end{array}$$

(37)

where

$${\omega }_{n1}=\sqrt{\frac{{\overline{C}}_{11}{l}^2-\overline{G}}{1+{\overline{J}}_T}}$$

(38)

is the dimensionless natural frequency of the rotor system (37) at ${\overline{C}}_1={\overline{C}}_{11}$ ,

$${\omega }_{n2}=\sqrt{\frac{{\overline{C}}_{12}{l}^2-\overline{G}}{1+{\overline{J}}_T}}$$

(39)

is the dimensionless natural frequency of the rotor system (37) at ${\overline{C}}_1={\overline{C}}_{12}$ , where ${\overline{J}}_T>>$ ${\overline{J}}_P$ is accepted.

In the case of resonant transitions close to steady-state modes from the equations of motion of the rotor, in accordance with (37), can obtain an energy relation that determines the power spent on damping the oscillatory motion

$$N=M\phi \prime =(u_1-u_2\phi \prime )\phi \prime ={\mu }_1({\alpha \prime }^2+{\beta \prime }^2)+{\mu }_3({\alpha \prime }^4+{\beta \prime }^4).$$

(40)

Results and discussion

Regular fluctuations and their analysis

The numerical results were obtained using the MATLAB package (R2021a (9.10.0.1602886) 64-bit (win 61) 17 February 2021). The calculations were carried out for the system parameters with the following values: e _r  = 0.0346, ω_n ≈ 1, J_P1 = 0.021 ( ${\overline{J}}_P=0.026,{\overline{J}}_T=0.213$ ) and the characteristic of the energy source was taken as a straight line (8). These system parameters were chosen in accordance with various design parameters of the centrifuge used experimental studies.^20,39 The values of the parameters C₃ = 0.1, μ₁, μ₃, u₁, u₂ are chosen in the course of numerical experimentation, as well as the values of μ₁ and μ₃, taking into account those values at which hopping effects disappear, in accordance with the rest of the known design parameters necessary to create effective vibration isolation for the gyroscopic rotor. In calculations using equations (25) and (26), the control parameters were taken $u_1\in (1.1,1.5)$ and depending on the type of engine $u_2=1.245$ from correct physical considerations. Solving the approximate equations (25) and (26), numerically finding the roots of the characteristic equation (30), and taking into account (31) and (32), the dependencies of the oscillation amplitude $a$ and angular speed of rotation of the shaft $\mathrm{\Omega }$ from the control parameter $u_1$ for different values of the coefficient of nonlinear cubic damping μ₃ = 0.01, 0.02 are plotted, in Figures 2 and 3. Stable and unstable analytical solutions are indicated by solid and dotted lines, respectively. The boundaries of the region of instability are found by jointly solving the stability criterion:

$$\frac{d\mathrm{\Omega }}{\mathrm{da}}\approx \frac{\partial f}{\partial a}=0$$

(41)

Figure 2.

Vibration near the main forced resonance at the coefficient of nonlinear cubic damping 0.01: (a) amplitudes and (b) angular velocity.

Figure 3.

Vibration near the main forced resonance at the coefficient of nonlinear cubic damping 0.02: (a) amplitudes and (b) angular velocity.

or

$$\begin{array}{c}[{\mu }_1{\omega }_n+0.75{\mu }_3{\omega }_n^3{a}^2][{\mu }_1{\omega }_n+2.25{\mu }_3{\omega }_n^3{a}^2]+\hfill \\ [0.75C_3{a}^2+2{\omega }_n({\omega }_n-\mathrm{\Omega })][2.25C_3{a}^2+2{\omega }_n({\omega }_n-\mathrm{\Omega })]=0\hfill \end{array}$$

(42)

and the frequency equation (25), and they are marked with small two triangles connected by a line (Figures 2(a) and 3(a)). Increasing the parameter $u_1$ , we observe a periodic movement, which is obtained as a result of the interaction of a forced oscillation and an energy source – a DC motor. In the area of $u_1\in (1.100,1.409)$ for μ₃ = 0.01 and $u_1\in (1.100,1.357)$ for μ₃ = 0.02, a synchronization effect is observed. There is an amplitude maximum resulting from the strong interaction between the limited excitation energy source and the vibrating rotor, with parameters at μ₃ = 0.01: $u_{1max}=1.409,a_{max}=1.507$ , at μ₃ = 0.02: $u_{1max}=1.359,a_{max}=1.207$ (Figures 2(a) and 3(a)). At C₃ = 0.1, as the value of μ₃ increases, the value of the parameter $u_{max}$ shifts in the direction of decrease. Leaving the region of synchronization, the system can pass from periodic motion to quasiperiodic motion and, finally, quasiperiodic oscillates. When the parameter $u_1$ decreases, the transition from quasiperiodic to periodic motion can be carried out through the region of instability, at μ₃ = 0.01 in the regions $u_1=1.37$ and $u_1=1.34$ and at μ₃ = 0.02 in the region $u_1=1.34$ . These transitions in Figures 2 and 3 are indicated by two dots connected by a line. With the same parameter value $u_1=1.34$ , the height of the transition through the region of instability for the case of μ₃ = 0.02 is much less than in the case of μ₃ = 0.01 (Figures 2(a) and 3(a)).

At the points of these transitions, the roots of the characteristic equation (30) have negative real parts. One root has no imaginary part and the remaining two have identical real parts in absolute value so small compared to unity; they can be approximately considered equal to zero.

It should be noted that natural transitions between periodic modes of motion under hopping effects also occur through the region of instability: in the case of μ₃ = 0.01 at $u_1=1.409,1.330$ ; in the case of μ₃ = 0.02 at $u_1=1.357,1.306$ . If the transition point from quasiperiodic to periodic motion is taken below the lower boundary of the instability region, then the upward transition can be carried out without crossing the instability region. In this case, in the case with μ₃ = 0.02, the transition point with coordinate parameters $u_1=1.260$ and $a=0.440$ is higher than a similar point with coordinate parameters $u_1=1.254$ and $a=0.398$ for the case with μ₃ = 0.01.

Thus, with an increase in the coefficient of nonlinear cubic damping μ₃ from a value of 0.01 to a value of 0.02 with constant linear damping with a coefficient μ₁ = 0.01, the “distances” of transitions from a periodic to a quasiperiodic motion mode and vice versa decrease. This is positively reflected in the safety of these transitions and in control of oscillations of large amplitudes with the regulation of the voltage on the motor and vibration amplitude.

In Figures 2(b) and 3(b), a deviation of the dependency graph is noticed $\mathrm{\Omega }=\mathrm{\Omega }\left(u_1\right)$ from a straight line, the size of this region in the case of μ₃ = 0.02 is smaller than the size of the indicated region in the case of μ₃ = 0.01. The size of the bistable region also decreases with an increase in the value of the coefficient of nonlinear cubic damping from 0.01 to 0.02. At the boundaries of the bistable region, which correspond to jumps tangent to the resonance curves vertically, that is,

$$\frac{\partial a}{\partial \delta }=\infty ,$$

(43)

where $\delta =\mathrm{\Omega }-{\omega }_n$ is the frequency detuning.

Differentiating (26) with respect to $\delta$ and taking into account (43), we obtain

$$({\mu }_1{\omega }_n+0.75{\mu }_3{\omega }_n^3{a}^2)({\mu }_1{\omega }_n+2.25{\mu }_3{\omega }_n^3{a}^2)+(0.75C_3{a}^2-2{\omega }_n\delta )(2.25C_3{a}^2-2{\omega }_n\delta )=0.$$

(44)

By resolving (26) and (44) together, we can find the values of the amplitude of oscillations and the speed of rotation of the shaft (or the detuning of the frequency of oscillations): $a_1$ , $a_2$ , and ${\mathrm{\Omega }}_1$ (or ${\delta }_1$ ), ${\mathrm{\Omega }}_2$ (or ${\delta }_2$ ) corresponding to jumping effects.

The extreme values of the oscillation amplitude can be determined from the following condition:

$$\frac{\partial a}{\partial \delta }=0$$

or

$$2e_r^2{(\delta +{\omega }_n)}^3+({\mu }_1+0.75{\mu }_3{\omega }_n^2{a}^2)({\mu }_1+2.25{\mu }_3{\omega }_n^2{a}^2){\omega }_n{a}^2+2(0.75C_3{a}^2-2{\omega }_n\delta )({\omega }_n-\delta )a^2=0.$$

(45)

By resolving the system from equations (26) and (45), it is possible to find the maximum resonant oscillation amplitude.

In order to determine the degree and nature of the influence of the energy source on the dynamics of an oscillating system in comparison with the case with an ideal energy source, we represent the system of equations of motion of an ideal gyroscopic rotor system in the following form:

$$\begin{array}{c}\alpha ″+{\omega }_n^2\alpha =e_r{\mathrm{\Omega }}^2\cos \mathrm{\Omega }\overline{t}-J_{P1}\mathrm{\Omega }\beta \prime -{\mu }_1\alpha \prime -{\mu }_3{\alpha \prime }^3-C_3{\alpha }^3,\hfill \\ \beta ″+{\omega }_n^2\beta =e_r{\mathrm{\Omega }}^2\sin \mathrm{\Omega }\overline{t}+J_{P1}\mathrm{\Omega }\alpha \prime -{\mu }_1\beta \prime -{\mu }_3{\beta \prime }^3-C_3{\beta }^3.\hfill \end{array}$$

(46)

Having solved equations (46) by the method of varying amplitude, we obtain the averaged equations of nonstationary oscillations of the rotor ²⁰ in the following form:

$$\begin{array}{c}a\prime =-0.5e_r{\mathrm{\Omega }}^2\sin \xi -0.5{\mu }_1\mathrm{\Omega }a-0.375{\mu }_3{\mathrm{\Omega }}^3{a}^3,\hfill \\ a\xi \prime =-0.5e_r{\mathrm{\Omega }}^2\cos \xi -(\mathrm{\Omega }-{\omega }_n-0.5J_{P1}\mathrm{\Omega })\mathrm{\Omega }a+0.375C_3{a}^3.\hfill \end{array}$$

(47)

Under conditions $a\prime =0$ and $\xi \prime =0$ from (47), we obtain an expression for determining the amplitude of stationary oscillations:

$$F(a,\mathrm{\Omega })=({\mu }_1\mathrm{\Omega }+0.75{\mu }_3{\mathrm{\Omega }}^3{a}^2{)}^2+{[0.75C_3{a}^2-2(\delta -0.5J_{P1}\mathrm{\Omega })\mathrm{\Omega }]}^2-{(e_r{\mathrm{\Omega }}^2/a)}^2=0.$$

(48)

Acting with function (48) in the same way as in the case of a rotor system with a nonideal energy source, we obtain the expressions necessary to determine the boundaries of the bistability region and the maximum resonant oscillation amplitude, respectively:

$$({\mu }_1\mathrm{\Omega }+\frac{3}{4}{\mu }_3{\mathrm{\Omega }}^3{a}^2)({\mu }_1\mathrm{\Omega }+\frac{9}{4}{\mu }_3{\mathrm{\Omega }}^3{a}^2)+[\frac{3}{4}{C}_3{a}^2-2(\delta -\frac{1}{2}{J}_{P1}\mathrm{\Omega })\mathrm{\Omega }][\frac{9}{4}{C}_3{a}^2-2(\delta -\frac{1}{2}{J}_{P1}\mathrm{\Omega })\mathrm{\Omega }]=0$$

(49)

and

$$\frac{3}{4}{C}_3{a}^2-2(\delta -\frac{1}{2}{J}_{P1}\mathrm{\Omega })\mathrm{\Omega }=0.$$

(50)

From the system of equations (48) and (49), we find the coordinates of the boundaries of the bistability region in terms of the oscillation amplitude a and the shaft rotation speed Ω; from the system of equations (48) and (50), we find the maximum resonant amplitude a_m and the corresponding shaft rotation speed Ωm.

The dependencies of the coordinates of the boundaries of the bistability region on the coefficient of nonlinear cubic damping for the positive and negative nonlinear coefficients of support stiffness are shown in Figures 4 and 5. Of these, it is noticeable that the case of a rigid nonlinear characteristic of the elastic support material for a rotor system with a nonideal energy source has a wider bistability area than for a rotor system with an ideal energy source, which contrasts with the case of a soft nonlinear characteristic of the elastic support material. In both variants of the nonlinear elastic characteristic of the support material, as the coefficient of nonlinear cubic damping μ₃ increases, the area of bistability narrows and when certain values of μ₃ are reached, it will disappear. The bistability region of the rotor system with a negative nonlinear support stiffness coefficient with the value C₃ = –0.1 is located in the speed range below the natural oscillation frequency ω_n, and the rotor system bistability region with a positive nonlinear support stiffness coefficient C₃ = 0.1 is located in the speed range above the natural frequency oscillations ω_n of system (6).

Figure 4.

Dependencies of oscillation amplitude (a) and shaft rotation speed (b), corresponding to the boundaries of the bistability region, on the coefficient of nonlinear cubic damping at a non-linear rigid characteristic.

Figure 5.

Dependencies of the oscillation amplitude (a) and shaft rotation speed (b), corresponding to the boundaries of the bistability region, on the coefficient of nonlinear cubic damping at a non-linear soft characteristic.

It should be noted that in vibrating nonideal systems, the most common case is with a positive value of the coefficient of nonlinear stiffness C₃ > 0 of the support material, that is, with a bistability region located behind the natural oscillation frequency ω_n. In this case, a rotor system with a nonideal energy source, to eliminate the bistability region, that is, jumping effects, will require more nonlinear cubic damping of support or energy from a nonideal energy source than in the case of an ideal rotor system.

The dependencies of the maximum amplitude of the frequency response of the rotor system on the coefficient of nonlinear cubic damping for the cases of hard and soft characteristics of the nonlinear elasticity of the support are shown in Figures 6 and 7. It is easy to see that with an increase in the value of the coefficient of nonlinear cubic damping μ₃, a significant decrease in the maximum amplitude a_m of the resonance curves and the approximation of the corresponding shaft rotation speed Ω_m to the natural frequency of the system ω_n ≈ 1, in the case of a rigid characteristic of the nonlinear elasticity of the support material from the side of the larger from ω_n the speed of rotation of the shaft, in the case of a soft characteristic of the nonlinear elasticity of the support material from the side of the shaft rotation speed Ω_m that is less from ω_n.

Figure 6.

Dependencies of the maximum resonant oscillation amplitude: (a) and the corresponding shaft rotation speed and (b) on the coefficient of nonlinear cubic damping with a nonlinear rigid characteristic.

Figure 7.

Dependencies of the maximum resonant oscillation amplitude: (a) and the corresponding shaft rotation speed and (b) on the coefficient of nonlinear cubic damping with a nonlinear soft characteristic.

In a nonideal rotor system with a support material with C₃ > 0, due to the strong interaction of the system with the energy source, the acquisition of additional energy by the system from the source for a steep increase in the response amplitude at a practically unchanged shaft rotation speed and the maximum amplitude of the frequency response is greater than in an ideal rotary system, the same value of C₃ support material.

Nonstationary resonant transition and their analysis

Now consider the nonstationary transition of the rotor through the resonant region with a “slowly” time-varying control parameter $u_1$ with regularity $u_1=u_{10}+\nu \epsilon \overline{\tau }=u_{10}+\nu \overline{t}$ . ⁴⁰ Here $\overline{t}=\epsilon \overline{\tau }$ , where $\overline{\tau }$ is the “slow” time and $\epsilon <<1$ is the small parameter.

The approximate equations of rotor motion (23) were modeled in the MATLAB -Simulink package (R2021a (9.10.0.1602886) 64-bit (win 61) 17 February 2021).

Damping parameters have the following values: μ₁ = 0.01; μ₃ = 0.01, 0.02. The control parameter $u_1$ changed slowly. The abscissa axis has a parameter scale $u_1$ and the corresponding time scale $\overline{t}$ .

The values of the parameters for the initial conditions are borrowed from the frequency characteristics of the stationary oscillation for Ω₀< ω_n and Ω₀ > ω_n.

For C₃ = 0.1 (rigid nonlinear characteristic of the elasticity of the support), we choose the initial conditions for the case ν = 0.00025: $\overline{t}$  = 0, Ω₀ = 0.88, u₁₀ = 1.096, a₀ = 0.11, ξ₀ = –0.0420 and for the case ν = − 0.00025: $\overline{t}$ =0, Ω₀ = 1.2, u₁₀ = 1.50, a₀ = 0.12, ξ₀ = 0.0253.

With C₃ = –0.1 (soft nonlinear characteristic of the elasticity of the support), we accept the following initial conditions for the case: ν = 0.00025: $\overline{t}$  = 0, Ω₀ = 0.88, u₁₀ = 1.096, a₀ = 0.112, ξ₀ = –0.0424 and for the case ν = –0.00025: $\overline{t}$  = 0, Ω₀ = 1.2, u₁₀ = 1.494, a₀ = 0.124, ξ₀ = 0.0253.

Dependency graphs $a=a\left(u_1\right)$ and $\mathrm{\Omega }=\mathrm{\Omega }\left(u_1\right)$ for nonstationary oscillations of the rotor, built on the results of modeling equations (23), are shown in Figures 8 to 11.

Figure 8.

Transition through resonance: (a) amplitude and (b) angular velocity, in the case of a nonlinear rigid characteristic with a slowly increasing value of the control parameter.

Figure 9.

Transition through resonance: (a) amplitude and (b) angular velocity, in the case of a nonlinear rigid characteristic with a slowly decreasing value of the control parameter.

Figure 10.

Transition through resonance: (a) amplitude and (b) angular velocity, in the case of a nonlinear soft characteristic with slowly increasing control value parameter.

Figure 11.

Transition through resonance: (a) amplitude and (b) angular velocity, in the case of a nonlinear soft characteristic with a slowly decreasing value of the control parameter.

From these figures, the damping effect of the parameter is primarily evident ${\mu }_3$ to the maximum $a_m$ and on the quasiperiodic variation of the oscillation amplitude behind the amplitude peak and the influence ${\mu }_3$ on the values of the control parameter $u_m$ and shaft speed ${\mathrm{\Omega }}_m$ corresponding to the maximum amplitude. Nonlinear cubic damping shifts the control parameter $u_m$ and shaft rotation speed ${\mathrm{\Omega }}_m$ corresponding to the maximum amplitude at $C_3\succ 0$ (Figure 8(a)) downwards (Figure 8(b)), whereas $C_3\prec 0$ (Figure 11(a)) shifts upwards (Figure 11(b)), that is, in both cases, the characteristics of the nonlinear stiffness of the support ${\mathrm{\Omega }}_m$ approaching ${\omega }_n$ .

The values of the maximum amplitude and the corresponding control parameter and shaft rotation speed in the resonance curves during acceleration (Figures 8 and 10) and run-out (Figures 9 and 11) of the rotary machine approximately determine the positions of the jumping effects. As the coefficient of nonlinear cubic damping increases, the distance between these positions decreases (the displacements of these positions in Figures 8 and 11 are much larger than in Figures 9 and 10), and its further increase can completely eliminate jumping phenomena.

Thus, with an increase in the nonlinear cubic damping, the Sommerfeld effect with a nonlinear jump can be significantly weakened up to its complete elimination. ⁴³

A change in the characteristic of the nonlinear stiffness of elastic support significantly affects the description of the dependency graphs $a=a\left(u_1\right)$ and $\mathrm{\Omega }=\mathrm{\Omega }\left(u_1\right)$ . The difference between the values of the maximum amplitude, the values of the corresponding control parameter, and the shaft rotation speed with increasing (Figures 8 and 10) and decreasing (Figures 9 and 11) control parameter (shaft rotation speed) is explained by hopping transitions with different values of these parameters during the run-up and run-out of the rotary machine. In the case of $C_3\succ 0$ and at $\nu \succ 0$ , the jump is carried out from a large amplitude to a smaller one (Figure 8), with $\nu \prec 0$ from a smaller amplitude to a larger one (Figure 9), and in the case of $C_3\prec 0$ vice versa (Figures 10 and 11, respectively). This is usually observed in experimental studies during acceleration and deceleration of the machine. ²⁰ . With a rigid nonlinear characteristic of the elasticity of the support $C_3\succ 0$ , jumps will be located in the region of shaft rotation speed, where $\mathrm{\Omega }\succ {\omega }_n$ (Figures 8 and 9), with a soft nonlinear elastic characteristic of the support $C_3\prec 0$ – in the range of shaft rotation speed, where $\mathrm{\Omega }\prec {\omega }_n$ (Figures 10 and 11). It must be taken into account that the jump transition from an unstable rotation speed with a larger amplitude to a stable shaft rotation speed with a smaller amplitude is unsafe.

After the transition of the amplitude peak with a slow change in the control parameter and the shaft rotation speed, the amplitude of the system response, having made damped oscillations, tends to certain values, regardless of the values of the coefficient of nonlinear cubic damping.

To verify the reliability of the process of transition through resonance, from equation (26), assuming to be equal to zero the expression $e{\mathrm{\Omega }}^2$ , due to the moment of inertia of the mass imbalance, and the damping coefficients ${\mu }_1$ and ${\mu }_3$ , we obtain the equation of the backbone curve:

$$\mathrm{\Omega }={\omega }_n+\frac{3C_3{a}^2}{8{\omega }_n},$$

(51)

where $\mathrm{\Omega }\succ {\omega }_n$ at $C_3\succ 0$ ; $\mathrm{\Omega }\prec {\omega }_n$ at $C_3\prec 0$ .

The maximum amplitudes and the corresponding shaft angular velocities and the control parameters of the resonance curves approximately satisfy equations (51) and (25). With system parameters ω_n ≈ 1, J_P1 = 0.021, μ₁ = 0.01, ν = 0.00025, correspondences between numerical (graphical) values of the maximum amplitude $a$ and resonant values $\mathrm{\Omega }$ and $u_1$ for given values of the linear damping coefficient μ₁ and different values of the coefficient of nonlinear stiffness $C_3$ and coefficient of nonlinear cubic damping μ₃ (Figures 8(a) and 10(a)) according to formulas (50) and (25) are given in Table 1.

Table 1.

Correspondence between oscillatory characteristics: numerical (graphical) values maximum amplitude, resonant values of rotation speed, and control parameter according to equations (51) and (25).

$C_3$	${\mu }_3$	a	$\mathrm{\Omega }$	$u_1$	$C_3$	${\mu }_3$	a	$\mathrm{\Omega }$	$u_1$
0.1	0.010	1.510	1.086	1.41	−0.1	0.010	1.239	0.942	1.21
	0.020	1.210	1.055	1.36		0.020	1.084	0.956	1.22

The values of oscillation parameters given in Table 1 completely coincide with the analytical results of solving equations (26) and (25).

For comparative analysis, we present the solutions of the system of equations (47) with previously known initial parameters of the system in the form of dependency graphs $a=a\left(\mathrm{\Omega }\right)$ at different values ${\mu }_3=0.01,0.02$ in cases with $C_3=0.1$ and $C_3=-0.1$ as shown in Figure 12 and the equation of the backbone curve

$$\mathrm{\Omega }=\frac{{\omega }_n}{2-J_{P1}}+\sqrt{{\left(\frac{{\omega }_n}{2-J_{P1}}\right)}^2+\frac{3C_3{a}^2}{4(2-J_{P1})}}$$

(52)

Figure 12.

Transition through resonance with a slow increase in the rotation speed in cases of non-linear (a) hard and (b) soft characteristic.

and a similar Table 2 for a rotor with an ideal energy source.

Table 2.

Correspondence between oscillatory characteristics: numerical (graphical) values of maximum amplitude, resonant values of rotation speed according to the backbone curve equation (52).

C 3	μ 3	a	Ω	C 3	μ 3	A	Ω
0.1	0.010	1.360	1.075	−0.1	0.010	1.163	0.960
	0.020	1.125	1.060		0.020	1.050	0.967

In formula (52), $\mathrm{\Omega }\succ {\omega }_n$ if $C_3\succ 0$ and $\mathrm{\Omega }\prec {\omega }_n$ if $C_3\prec 0$ .

By comparing the tabular values of the maximum amplitude and the corresponding shaft rotation speed at different values of the coefficient of nonlinear cubic damping for cases with a nonideal source and an ideal energy source, one can verify the reliability of the previously given conclusions.

Expressions (44) and (49) without frequency response equations (26) and (48) determine the boundaries of the instability regions for the cases of a rotor with nonideal and ideal energy sources. The boundaries of the instability region constructed according to the stability criteria (44) and (49) at the value of the linear damping coefficient μ₁ = 0.01 and different values of the coefficient of nonlinear cubic damping μ₃ are shown in Figure 13.

Figure 13.

The boundaries of the region of instability at the coefficient of nonlinear cubic damping: (a) 0.01 and (b) 0.02: (1) the case of a nonideal energy source and (2) the case of an ideal energy source.

In both cases, as the value of μ₃ increased from 0.01 to 0.02, the instability regions narrowed. It is clearly seen that in the case of an ideal energy source, the dimensions of the instability region along the Ω and a axes are noticeably smaller than those in the case of a nonideal energy source; this is more distinct at μ₃ = 0.02.

With a further increase in the value of μ₃, the dimensions of the instability region decreased: in the case with a nonideal energy source along the a axis, along the Ω axis from the side of the origin, and in the case with an ideal energy source along both axes comprehensively. The regions of instability disappear: in the case of an ideal energy source at μ₃ = 0.03788; in the case of a nonideal energy source at μ₃ = 0.0565.

The power consumed for damping the oscillatory motion can be obtained from equation (25) by multiplying it by Ω:

$$N=({\mu }_1+0.75{\mu }_3{\omega }_n^2{a}^2){\omega }_n^2{a}^2.$$

(53)

On the other hand, the power of the energy source and the power expended to overcome the forces of resistance to the rotation of the motor rotor are as follows:

$$N\left(\mathrm{\Omega }\right)=(u_1-u_2\mathrm{\Omega })\mathrm{\Omega }.$$

(54)

From previous studies, it was found that the greater the probability of interaction with an energy source in the resonant region, with a slow increase in the speed of rotation of the shaft, the amplitude of oscillations increasing steeply, and in order to suppress oscillations with large amplitudes in accordance with (47) and (48), the more the damping and power of the energy source ((53) and (54)).

To confirm the analytical studies, equations (10) were solved directly numerically using the MATLAB-Simulink package (R2021a (9.10.0.1602886) 64-bit (win 61) 17 February 2021). Figure 14 shows the obtained numerical results in the form of dependency graphs $\alpha =\alpha \left(u_1\right)$ , to pass through the resonance with a slowly increasing control parameter $u_1$ in cases with hard and soft characteristics of nonlinear elasticity of the support $C_3$ .

Figure 14.

Transition through resonance with ν = 0.00025 according to the results of the numerical solution of equations (10) for non-linear: (a) hard and (b) soft characteristic.

An analysis of these oscillograms shows that with increasing the coefficient of nonlinear cubic damping, the maximum and beyond the maximum oscillation amplitudes decrease. Oscillations in the region behind the resonance have the character of damped beats, which are more pronounced for small values of the coefficient of nonlinear cubic damping than for large values. The difference between Figures 8(a), 10(a), and 14 lies in the magnitude of the maximum amplitude and its corresponding value of the control parameter, but the basic behavior of the transient through resonance remains the same.

Now, for a comparative analysis of the results of direct modeling of the equations of motion of a nonideal and ideal system, we present in Figure 15 results of numerical solution of equations (46).

Figure 15.

Transition through resonance with ν = 0.00025 according to the results of the numerical solution of equations (43) for non-linear: (a) hard and (b) soft characteristic.

From the oscillograms in Figures 14 and 15, it can be seen that the response amplitudes (including maximum amplitude) of a nonideal system in the transition region are much greater than the vibration amplitudes of an ideal system, which is the result of a strong interaction of the vibrating rotary machine with the engine in the shaft speed region $\mathrm{\Omega }\succ {\omega }_n$ . The moments of passage with the maximum amplitude are different, but the behavior of the transient process in a nonideal and ideal system is generally the same.

Damped beats observed in Figures 8 to 12, 14, and 15 occur due to the superposition of forced nonstationary oscillations and damped natural oscillations with frequencies that closely coincide with the vicinity of the resonance. ²²

According to equations (10), the direct modeling results can also be provided in the form of phase trajectories $\alpha \prime =\alpha \prime \left(\alpha \right)$ and the dependence graph $\alpha =\alpha \left(\beta \right)$ of oscillating systems with different characteristics of nonlinear elasticity and different values of the coefficient of nonlinear cubic damping of the support, in Figures 16 and 17. Figure 16 shows that the value ${\mu }_3$ dampens not only angular misalignment $\alpha$ but also the angular velocity $\alpha \prime$ ; the sign change $C_3$ affects the size and density of phase trajectories.

Figure 16.

Phase trajectories: (a) C₃ = 0.1, μ₃ = 0.01, (b) C₃ = 0.1, μ₃ = 0.02, (c) C₃ = −0.1, μ₃ = 0.01, and (d) C₃ = −0.1, μ₃ = 0.02

Figure 17.

The configuration space (beta, alpha): (a) C₃ = 0.1, μ₃ = 0.01, (b) C₃ = 0.1, μ₃ = 0.02, (c) C₃ = −0.1, μ₃ = 0.01, (d) C₃ = −0.1, μ₃ = 0.02.

Increasing value ${\mu }_3$ from 0.01 to 0.02 reduces the range of fluctuations in angular coordinates $\alpha$ and $\beta$ and changes the size of the shape of the dependence graph $\alpha =\alpha \left(\beta \right)$ , and the sign change $C_3$ changes the type of this graphical dependence (Figure 17).

For numerical calculations by (37), some geometric and dynamic parameters of the gyroscopic rotor system are taken from the experimental setup: L = 0.46 m, l₀ = 0.33 m, m = 2 kg, c₁₁ = 2 10⁴ N/m, c₁₂ = 2.4 10⁴ N/m, J_T = 0.090 kgm², J_P = 0.011 kgm², e = 0.0193 m used in study. ²⁰

Now, according to formula (34), taking into account (35), it is possible to calculate the natural frequencies (critical speeds) of the gyroscopic rotor system (33): ω₁ = 64.94 s⁻¹, ω₂ = 71.34 s⁻¹.

Dependencies ω₁ and ω₂ from c₁₂ at constant c₁₁ are shown in Figure 18. Hence, a noticeably weak dependence of ω₁ on c₁₂ and ω₂ decreases as c₁₂ decreases, approaching the value of ω₁.

Figure 18.

Dependences of critical speeds on the linear stiffness of the support material.

The angular speed of rotation of shaft ω₁ in the interval [64.95, 64.48] can be approximately considered constant and equal to ω₁ = 65 s⁻¹. This allows dimensionless parameters $\overline{G},{\overline{C}}_{11},{\overline{C}}_{12},C_3,{\mu }_1,{\mu }_3,u_1$ and $u_2$ independent of changes in parameter ω₁.

Calculation by formulas (36), (38), and (39) will provide the following values of dimensionless parameters: $l=0.72;$ ${\overline{J}}_p=0.026;$ ${\overline{J}}_T=0.213;$ $\overline{G}=0.0051;$ ${\overline{C}}_{11}=2.37;$ ${\overline{C}}_{12}=2.84;$ $e_r=0.0346;$ $J_{p1}=0.021;$ $C_3=0.1;$ ${\mu }_1=0.01;$ ${\mu }_3=0.01,0.02;$ $u_2=1.263,$ ${\omega }_{n1}=1,$ ${\omega }_{n2}=1.1$ .

Figures 19 and 20 show waveforms α = α(u₁) and β = β(u₁) at ${\mu }_1=0.01,{\mu }_3=0.01;0.02$ , $\nu =\pm 0.00025$ and for different values ${\overline{C}}_{12}$ .

Figure 19.

Oscillograms of angular coordinates α and β at $\nu =0.00025$ and different values ${\overline{C}}_{12}$ , at: (a) ${\mu }_3=0.01$ and (b) ${\mu }_3=0.02$ .

Figure 20.

Oscillograms of angular coordinates α and β at $\nu =-0.00025$ and different values ${\overline{C}}_{12}$ , at: (a) ${\mu }_3=0.01$ and (b) ${\mu }_3=0.02$ .

Control parameter $u_1$ depending on time $\overline{t}$ changes linearly $u_1=u_{10}+\nu \overline{t}$ .

To simulate a direct resonant transition, the initial conditions are taken for sufficiently small angular coordinates α and β of stationary vibrations, under which the anisotropy properties of the support material rigidity practically do not manifest themselves: $\overline{t}$  = 0, $\alpha$  = 0.055, $\alpha \prime$  = 0.000019, $\beta$  = –0.000024, $\beta \prime$  = 0.044, $\phi$  = 0, $\phi \prime$  = 0.8, $u_{10}=0.996$ .

At a slowly increasing parameter $u_1$ the value of $\nu$ equals 0.00025.

In the oscillograms α = α(u₁) of Figure 19, two resonance regions are well noticed, first resonance region in terms of manifestation area and maximum amplitude of oscillations prevails in them, and the second small resonance region is observed in the beating region. In the oscillograms β = β(u₁) of Figure 19 in terms of manifestation area and maximum oscillation range the second resonance region prevails. When the value ${\overline{C}}_{12}$ approaches the value ${\overline{C}}_{11}$ the difference between the values of critical speeds is reduced the second resonant region shifts toward the first, at ${\overline{C}}_{12}=2.84$ a practically separate region of the second critical velocity of the oscillogram β = β(u₁) increases to occupy the regions of the first critical velocity of the oscillogram α = α(u₁), until it is fully occupied at ${\overline{C}}_{12}=2.37.$ Full fusion of these regions is clearly visible at ${\mu }_3=0.02$ .

A nonlinear jump from a large amplitude of oscillations to a smaller amplitude is observed in the first resonance region of the oscillogram α = α(u₁) and the second resonance region of the oscillogram β = β(u₁). If runout in oscillograms α = α(u₁) is observed after the jump effect, in oscillograms β = β(u₁) – in the region with the control parameter $u_1\succ 1.371$ .

A comparative analysis of the corresponding oscillograms in Figure 19(a) and (b) shows that nonlinear cubic damping significantly suppresses the amplitudes of oscillations in the resonance regions and in the runout region, especially in the oscillograms α = α(u₁) in the interval ${\overline{C}}_{12}\in (2.84,2.38)$ and in all oscillograms β = β(u₁). When ${\mu }_3=0.02$ runout is observed in the interval ${\overline{C}}_{12}\in (2.84,2.38)$ , approximately in the region with the control parameter $u_1\succ 1.371$ , with ${\overline{C}}_{12}=2.37$ in the region with the control parameter $u_1\succ 1.434$ .

When choosing the initial conditions to simulate the reverse resonant transition the equations of motion (51), we took into account the results of solving the same equations for a direct resonant transition. Thus, the initial conditions for μ₃ = 0.01: $\overline{t}$  = 0, $\alpha$  = 0.02902, $\alpha \prime$  = 0.1149, $\beta$  = −0.1221, $\beta \prime$  = 0.0648, $\phi$  = 3111.973, $\phi \prime$  = 1.3636, u₁₀= 1.746, and at μ₃= 0.02: $\overline{t}$  = 0, $\alpha$  = 0.04828, $\alpha \prime$  = 0.0877, $\beta$  = −0.0812, $\beta \prime$  = 0.07385, $\phi$  = 3246.607, $\phi \prime$  = 1.3824, u₁₀= 1.746.

At a slowly decreasing parameter $u_1$ the value of $\nu$ equals −0.00025.

From the oscillograms α = α(u₁) in Figure 20, two resonant positions are clearly visible with the rotor machine runout and at the values ${\overline{C}}_{12}=2.84;2.50$ . Both positions correspond to a jump from a smaller amplitude of oscillations to a larger amplitude in resonant regions. As the value decreases, ${\overline{C}}_{12}$ the first resonant region approaches the second. If the case with ${\overline{C}}_{12}=2.50$ corresponds to the maximum approximation of these areas, then the case with ${\overline{C}}_{12}=2.38$ – to the beginning of the merger, and the case with ${\overline{C}}_{12}=2.37$ - to their complete merger. When the value of ${\overline{C}}_{12}=2.37$ the oscillogram β = β(u₁) completely merges with the oscillogram α = α(u₁). Comparison of the oscillograms of Figure 20(a) and (b) with the corresponding values of ${\overline{C}}_{12}$ shows the damping effect of the magnitude ${\mu }_3$ .

Graphic dependencies of the angular speed of rotation of the shaft on the control parameter: $\phi \prime =\phi \prime \left(u_1\right)$ , the dependence of the dynamic torque of the engine on the control parameter according to the formula $M=u_1-u_2\phi \prime$ , and the dependence of the power expended on the damping of the oscillatory motion on the control parameter according to the formula $N=(u_1-u_2\phi \prime )\phi \prime$ at values μ₃ = 0.02 and at $\nu =\pm 0.00025$ are shown in Figures 21 and 22. From these graphs, first of all, you can observe the deviations of the dependencies $\phi \prime =\phi \prime \left(u_1\mathrm{or}\overline{t}\right)$ from a straight line and change $\phi \prime$ depending on $\overline{t}$ or $u_1$ , similar dependencies $M(u_1\mathrm{or}\overline{t})$ and $N(u_1\mathrm{or}\overline{t})$ , and hopping effects in resonant regions. ⁴⁴ . When two critical speeds are close to each other, for example, when ${\overline{C}}_{12}=2.50,2.38$ (Figure 21) exit from resonance at a lower speed, it can lead to capture at another resonance at a higher speed. The average power spent on fluctuations in angular displacements and on increasing the speed of rotation of the rotor shaft in the second resonant region is greater than that in the first, respectively, and possibly, the harm from the Sommerfeld effect in the region of the second critical velocity is greater than in the region of the first critical velocity ( ${\overline{C}}_{12}=2.84,2.50$ ). With a further decrease ${\overline{C}}_{12}$ and reaching a value of 2.38, the severity of the Sommerfeld effect on each of the resonant regions becomes comparable, and at ${\overline{C}}_{12}=2.37$ , two resonant regions merge into one. An increasing beat is observed in the second resonant region. Then, after the hopping effect, the oscillatory process passes into damped oscillations.

Figure 21.

Dependences of the rotation speed, dynamic torque and power of the motor shaft on the increasing control parameter at the coefficient of nonlinear cubic damping 0.02.

Figure 22.

Dependences of the rotation speed, dynamic torque and power of the motor shaft on the decreasing control parameter at the coefficient of nonlinear cubic damping 0.02.

In Figure 22, the positions and magnitudes of the maximum amplitude in the resonant regions correspond to nonlinear jumps from a smaller amplitude to a larger one with decreasing values ${\overline{C}}_{12}$ and decreasing parameter u₁. At ${\overline{C}}_{12}=2.84,2.50$ in the resonant areas, large oscillatory processes in the values of M and N are observed in significant ranges $u_1$ and $\phi \prime$ , which is undesirable.

Thus, during the acceleration of the rotary machine, due to the anisotropy of the linear stiffness of the support material, the two critical speeds that appear are so close to each other that the exit from the first resonance at a lower speed can lead to capture at the second resonance at a higher speed, the power acquired by the system in the second resonant region is several times greater than that in the first, at least comparable. In a significant range of values of the control parameter or shaft rotation speed, the severity of the Sommerfeld effects is significant. In both resonant regions of speed capture, excess engine power is spent on increasing the amplitude of oscillations of angular (linear) displacements. Therefore, the isotropy of the linear stiffness of the support material is a desirable property for designing a vibration isolator of the machine.

Subsequent to the results of present and previous studies, and on the basis of the model developed in work of Iskakov and Bissembayev ²⁰ 3D, a prototype of a centrifuge installation based on a gyroscopic rigid rotor (Figure 23) and the structure of its elastic support (Figure 24) were created for which a patent request for invention was filed. The installation consists of the main parts: the rotor bowl, the shaft, the motor, the lower hinged and the upper elastic support and the outer housing attached to the platform. The rotor is made in the form of a cylindrical bowl to satisfy the condition $J_T>>J_p$ and is mounted on the upper end of the shaft resting on the lower cardan joint and the upper elastic support. The electric motor rigidly placed inside the cylindrical casing is located closer to the lower hinged support. An encoder for feedback to the control unit is attached to the lower part of the electric motor. The elastic support is made in the form of rubber plates with the shape of a rectangular parallelepiped, tightly installed in six grooves adjacent to the structure coupling. The shaft rotates freely on the bearing located in the coupling. Such a design ensures easy changeability of materials with different elastic and damping characteristics. By changing the number and location of rubber plates, a combination of plates with different elastic and damping properties it is possible to adjust the rigidity and damping of the support. In the case of determining the range of operating speeds for the critical frequency(s) of rotation of the shaft, the rubber pads are selected so that the damping of the specially manufactured rubber or elastic material helped to avoid the jumping with Sommerfeld effects in the resonance curve(s) resulting from the action of the nonlinear component of the elastic force. This in turn enables the rotor to pass safely through the critical speed(s).

Figure 23.

Centrifuge based on vertical rigid gyroscopic rotor: 1 – rotor bowl, 2 – electric motor, and 3 – control panel.

Figure 24.

Elastic support structure: 1 – rubber plate, 2 - coupling, 3 – rotor shaft, and 4 – bearing.

Conclusion

Differential equations of motion of a gyroscopic rigid unbalanced rotor with nonlinear cubic damping and nonlinear rigidity of elastic support are constructed, taking into account the interaction with a nonideal energy source. These equations were solved by Bogolyubov’s asymptotic method and numerical methods.

In the constructed graphic dependencies of the amplitude of oscillations and the speed of rotation of the shaft on the control parameter, possible options for controlling resonant oscillations of large amplitudes are given, including the most effective, by the method of controlling the control parameter (voltage on the motor), the vibration amplitude, and the angular velocity of the shaft in the frequency equation.

It is proved that in the case of a rigid characteristic of the nonlinear elasticity of the support material, in a rotor system with a nonideal energy source, in order to eliminate the bistability region, that is, jumping effects and weakening of the Sommerfeld effect, more nonlinear cubic damping of the support or energy from a nonideal energy source is required than in the case of an ideal rotor system.

It is shown that during the resonant transition, the nonlinear cubic damping of the support significantly suppresses not only the vibration with the maximum amplitude but also the resonant damped beating of such vibrations. It shifts the control parameter and the angular velocity of the shaft corresponding to the maximum amplitude with a rigid nonlinear elastic characteristic of the support material toward a decrease and with a soft nonlinear elastic characteristic of the support material toward an increase.

The effect of the nature of the nonlinear stiffness and the anisotropy of the linear stiffness of the support material on the maximum oscillation amplitudes with Sommerfeld effects and the corresponding control parameters (shaft rotation velocities) with increasing and decreasing control parameter is studied.

There is an agreement between the results of analytical solutions and numerical solutions of the rotor motion equations for transient processes.

The research results can be used in scientific research and design calculations to create elastic supports with the best damping parameters for a vibrating gyroscopic rotary machine.