THE INTERRELATED MODELLING METHOD OF THE NONLINEAR DYNAMICS OF RIGID ROTORS IN PASSIVE AND ACTIVE MAGNETIC BEARINGS

 G. Martynenko, 2016 However, an AMB allows for a rapid change of dynamic parameters such as bearing stiffness and damping properties. PMBs realise the self-aligning effect by using permanent magnets in different design versions though fail to effect full suspension merely by using them. Earnshaw [7] proved this theoretically, and Brounbeck [8] confirmed this except for cases of embedding permanent magnets in a ferromagnetic fluid [6] or for bodies made of diamagnetics or superconductors under usual atmospheric conditions [9]. In the majority of cases, such a solution is impractical and technologically unfeasible for industrial rotor machinery. Hence, based on the possibilities of practical implementation of complete magnetic bearings of rotors, this study considers the options of using either radial or axial AMBs for stabilising a rotor over all five degrees of freedom or one AMB jointly with several PMBs in different design versions. The most practical approach would be to use a combination of magnetic bearings of different types in medium-sized high-speed rotor machinery, e. g., turbo-expanders, expander-generator, and expander-compressor units [4]. They can use two radial PMBs and one axial AMB arranged in the centre or at one end of the shaft. This is due to the design features such as the presence of one or two discs arranged on the rotor cantilevers [10].


Introduction
A magnetic bearing (MB) is one of the variants of elastic-damping bearings.Its feature is the use of magnetic fields to provide stable rotor levitation.These fields create bearing force responses to rotor displacement in order to ensure automatic alignment of its bearing areas in the MB stator elements and a required level of bearing stiffness.MBs that are the most applicable from the practical viewpoint are active magnetic bearings (AMBs) [1][2][3][4][5] and passive magnetic bearings (PMBs) [6].Some MB design options are shown in Fig. 1.It shows radial and axial AMBs with electromagnets (Figs. 1, a, b) and a radial PMB with permanent annular magnets (Figs. 1, c, d), and the following notations are introduced: 1 -rotor; 2 -stators; 3 -AMB windings; 4 -AMB position sensors; 5 -comparator in AMB control system; 6 -AMB control device; 7 -amplifiers feeding control voltages to AMB windings, which are formed according to the accepted control algorithm; 8 and 9 -movable and stationary annular permanent magnets.
The stability of rotor motion in an AMB is ensured by an automatic control system and its control algorithm.The presence of this complex and costly component, and also of power sources for control electromagnets, is one of the major drawbacks of an AMB as compared to a PMB.In considering the problem of describing the dynamics of rotors in different power machinery in which magnetic bearings are used as rotor bearing assemblies, a conclusion can be made on the need to develop special approaches to mathematical modelling with account for all specific features that this kind of elastic damping bearing introduces into rotor systems.Presently, there is a large variety of studies on this subject.However, they offer a simplified approach to mathematical modelling of rotor dynamics.Therefore, building refined mathematical models will enable increasing the accuracy of numerical computation of required dynamic parameters of rotors, magnetic bearings and control systems for active magnetic bearings.This will dramatically reduce the amount of experimental investigations and increase the effectiveness of research and development efforts.

Analysis of the studies and statement of the problem
Assessing the state of problems related to analysis and synthesis allows defining one of the most topical research subjects in the field of magnetic bearings [11].This is the problem related to mathematical modelling of the dynamics of systems with magnetic bearings.All the components in the structure of MBs of different types, namely permanent magnets in a PMB or electromagnets, power amplifiers, controllers and even position detectors in an AMB are characterised by nonlinear relationships of different parameters.Hence, any system with an MB is nonlinear in essence.In such a system, nonlinear phenomena of rotor dynamics can occur, as this was shown in [12] for rotors in nonlinear bearings of other types.
As follows from analysis of studies [1,2,4,5], nonlinearity is neglected, as a rule, when modelling the dynamic behaviour of rotor systems with MBs, and linear mathematical models are used to simplify analysis.The main reason of this is the challenges in implementing nonlinear approaches.At the same time, failing to account for nonlinear properties inherent to the system can result during modelling to substan-tial misrepresentation of results, both quantitative and qualitative ones (e.g., failing to represent certain phenomena).
Linear mathematical models are built by linearizing in the "zero point", i. e. in the rotor's central position, the nonlinear electromagnetic (in an AMB) or magnetic (in a PMB) forces.The linearized dependencies of these forces on gaps (in an AMB and a PMB) and currents (in an AMB) allow proceeding to the known linear mathematical model by replacing an MB with linear springs with damping [13].However, a linear approximation of a magnetic force can be acceptable only in a small vicinity, i. e.only for small currents and small rotor deflections from the position of equilibrium.For instance, if rotor deflection exceeds one-half of the gap, the total magnetic force created by the MB can differ by 40 % from the value obtained from the linearized relationship [14].Actually, depending on the MB features, this difference can be yet greater, resulting in incorrectness of modelling rotor dynamics and preventing to use the entire potential of a magnetic bearing when designing based on such models.Therefore, numerical analysis of the effect of nonlinearity on the dynamic behaviour of the "rotor in an MB" system with application of refined mathematical models is a highly topical problem.Its solution will allow, first, to obtain fundamental scientific data on system performance for different operating conditions and, second, to use parametric modelling to find optimal design solutions to ensure dynamic stability and improve reliability and fail-safe operation.In doing so, one should take into account that the quality of the active magnetic suspension depends directly on the control law and an algorithm implemented by the control system controller.Its development directly depends on in-depth knowledge of the system's dynamic characteristics.If the control unit shall be developed based on erroneous modelling, this can lead to operational failures.
One can distinguish the following most essential causes of nonlinearity of the "rotor in MB" system, which limit the possibilities of linearising mathematical models and using linear approaches to analysing phenomena occurring in the system [11]: (1) nonlinear dependence of magnetic forces on gaps between moving and stationary parts in a PMB and an AMB, and on currents in the windings of AMB electromagnets; (2) current delay in the windings of AMB electromagnets, i. e. nonlinearity related to coil inductance; (3) geometric links between the electromagnets of one AMB and links between all AMBs in one rotor, which result in coupling of processes in orthogonal directions; (4) delays in the control system controller; (5) presence of eddy currents and dissipation fluxes among other factors.Investigation of all these factors has shown that each of them can result in origination of subharmonic and superharmonic vibrations, nonlinear resonant response, amplitude failures and jumps, chaotic motion, and also loss of rotor motion stability, whereas the overall action of nonlinear elements can only aggravate the situation.
In spite of the variety of studies on modelling the dynamics of rotors in an MB with account for nonlinearities, the majority of them have been done while considering systems with one or two degrees of freedom and only with account for some of the above-mentioned causes of nonlinearity.For instance, in [15], this is the nonlinear dependence of magnetic force presented by a simplified expression in the form of a piecewise linear characteristic with account for a linear control system with feedback and possible saturation of its force elements.A downside of the approach to modelling dynamics presented in the studies [16,17] accounts for only nonlinear effects caused by the geometric interrelation of two pairs of electromagnet poles that create forces in mutually perpendicular directions and the dependence of forces on displacements and currents in AMB coils.The study [18] considers only nonlinearities occurring due to the effect of a thrust AMB on radial vibrations.
Models with two degrees of freedom can partially account for certain geometric interrelations.In the study [19], such a case is mathematically described by equations of motion in the horizontal and vertical directions with quadratic and cubic terms and excitation parameters.Based on this model, the study [20] has shown that in a system with variable bearing stiffness jumps from mode to mode can occur.Paper [21] offers a series of studies conducted by using numerical integration of equations of motion with account for both soft and rigid force characteristics of radial AMBs.These studies have proved the likelihood of occurrence of a large variety of nonlinear phenomena in the system.However, these mathematical models accounted for but several nonlinear dependencies.
The present study is an attempt to eliminate the drawbacks of the existing models and to build new mathematical models describing rotor dynamics in MBs of different kinds with account for as much as possible nonlinear effects inherent in such systems.

Research objective and tasks
The objective of the research is to develop a common method for mathematical description of the dynamic behaviour of rotors in systems with magnetic bearings of different kinds.Its distinctive features are approach generality and completeness of accounting for the nonlinear interrelations of process occurring in such a system -electric, magnetic, and mechanical.
To achieve the stated research objective, the study specifies and solves the following tasks: -to develop a common approach to building analytical models for a mathematical description of the dynamics of rigid rotors in complete magnetic suspensions that are implemented by magnetic bearings of different kinds, with the models accounting for the nonlinear interrelations of processes associated with the majority of practically significant physical, design and functional features of the system; -to adapt the analytical method of analysing electromagnetic circuits of active magnetic bearings to find expressions for magnetic energy and forces with account for the control laws to introduce the dynamics of rotors in magnetic bearings into mathematical models; -to develop an imitation numerical model and use it for numerical research to prove the adequacy of the approach suggested by comparing numerical and experimental data.

Description of rotors in magnetic bearings as dynamic electromagnetic mechanical systems
One of the methods of describing the dynamics of electromechanical systems is using Lagrange-Maxwell equations, having the equations of mechanics structure: where Т and U are kinetic and potential energy respectively; W¢ is magnetic field co-energy; V is electric field energy; R is dissipation function; q r are generalised mechanical coordinates; Q r are nonpotential generalised forces; M is the number of generalised mechanical coordinates; D is an electric dissipation function; E k is an algebraic sum of extraneous electromotive forces; i k are circuit currents; χ k are charges of capacitors; and N is the number of closed-loop circuits, whereas: where B k are magnetic inductions of magnetic circuits; L ks are induction coefficients of electromagnetic circuits; ε is a dielectric constant; C k are capacitances of capacitors; and r ks are active resistances.
The study considers a complete rotor magnetic suspension in two radial and one axial magnetic bearings.Such a variant is a most common one in small and medium-sized rotor machinery for different purposes (centrifugal compressors and pumps, turbo-expanders, and expander-compressor units).However, the approach suggested for mathematical modelling can be used for describing the dynamics of the "rotor in an MB" system with any number of magnetic bearings of any type.The design schematic for a rigid rotor in an MB is shown in Fig. 2.
equal to d r1 , d r2 , and d a , respectively.The rotor rotates with a constant angular speed ω.To determine rotor positions, the study suggests use of generalised coordinates q={x 1 , y 1 , x 2 , y 2 , z 3 }, i. e. the coordinates of the centres of MB thrust areas, with their total number being M=5.In this case, when conductance currents are closed and there are no capacitors in electric branches, electromechanical systems, apart from Lagrange-Maxwell equations in the form (1), can be described by equations identical to Routh equations in mechanics [22]: where W is energy of the MB magnetic field (W+W¢=Ψi); N is the number of closed-loop circuits with currents i k ; Ψ k is an induction flux (magnetic flux linkage); r C ks are active resistances of electric circuits; and E k is an algebraic sum of extraneous electromotive forces in the k-th circuit, with: where L ks =L ks (x 1 , y 1 , x 2 , y 2 , z 3 ) are self and mutual inductance coefficients of circuits.The mathematical model ( 5) for the case of an AMB in the system contains differential equations for the flux linkage Ψ k , corresponding to Kirchhoff's second law for magnetic circuits, and they are a form of notation of Ampere's law for each k-th electric circuit in the system (AMB windings).Such an approach accounts for resistance forces and delays in such an electromechanical system as an AMB.
The kinetic energy of a rigid rotor has the form: where m is the rotor mass; v 0 is the pole O vector velocity; Ω is an angular speed vector; R′ C is a radius vector of OC; and Θ 0 is an inertia tensor in pole O.
In applying Euler's angles, the projections of vectors of pole velocity v 0 and of angular speed Ω on the system of coordinate axes shall be written down using directional cosines.By using the expansion of trigonometric functions into exponential series, expressing these angles through generalised coordinates, substituting the expression for T into the Lagrange equation and retaining only the first, second and third-order terms for generalised coordinates, their derivatives and unbalance parameters, we obtain nonlinear motion equations for a rigid rotor in an MB.Such a structure of representing nonlinear inertia forces corresponds to the structure of potential forces in an MB, and is sufficient for a correct description of nonlinear vibrations.
Terms -∂W/∂q r in equations ( 5) are ponderomotive forces, i.e. generalised forces caused by the mechanical action of a magnetic or electromagnetic field.These forces are suggested to be determined by differentiating for generalised mechanical coordinates of the magnetic energy function obtained by using equivalent circuits for the AMB electromagnetic circuits.The practical application of the suggested approach is further demonstrated whilst analysing the dynamics of one of the variants of a complete rotor magnetic suspension.

Mathematical model of rotor dynamics in a laboratory setup
The dynamic behaviour of a rotor in an MB was analysed for a laboratory setup of a rotor in a complete combined passive-active magnetic suspension, which was a prototype of the magnetic suspension for a rotor in an expander compressor unit (ECU).A schematic diagram of a complete magnetic suspension of a rotor, including two radial PMB1.2 and one axial AMB3, is shown in Fig. 3. Here, a PMB (Fig. 1, с) is used as PMB1 and PMB2, and an AMB is used as AMB3 (Fig. 1 If the energy of the AMB3 magnetic field is W=W(x 1 , y 1 , x 2 , y 2 , z 3 , Ψ c1 , Ψ c2 ), then the currents in the windings of its coils i c1 and i c2 are linked to total magnetic fluxes through the circuits of coils Ψ c1 , Ψ c2 (flux linkage in windings of respective AMB3 electromagnets) by the following expressions: We also assume that all generalised coordinates -displacements x 1 ,..., z 3 , unbalance parameters e 1 , e 2 and γ 1 , γ 2 , and gaps in the MB -d r1 , d r2 , and d a -have the same order of smallness.Then, with account for this assumption on smallness of generalised coordinates and their derivatives, the nonlinear addends of equations of motion can be considered small as compared to the linear terms.By excluding from consideration the addends of the equations of motion, with an order of smallness higher than three, we derive a completely coupled system of seven nonlinear differential equations describing the dynamics of this electromagnetic mechanical system: where f′′ qr (x 1 ,...,z 3 ) and f′′′ qr (x 1 ,...,z 3 ) are nonlinear terms of the equations of motion due to inertia forces and the second and third-order potential field; b x1,...,z3 are viscosity coefficients; r c 1,...,N are active resistances in winding circuits; u c 1,…,N are control voltages applied across AMB windings whose magnitude is formed according to the accepted control law depending on the rotor current position ; m ks , j are inertial and gyroscopic coefficients with the following values: Addenda -F Mqr are potential forces that depend only on the generalised coordinates.In this case, these are the magnetic forces in PMB1.2 (Fig. 3, b).The magnetic forces dependencies were obtained using the Maxwell tension tensor by solving a series of magnetic statics problems in the finite element statement for a fixed number of rotor magnet positions corresponding to certain discrete values of its displacement, though they can be described by analytical expressions as in [23].Terms - ¶W/ ¶q r are ponderomotive forces, i. e. the electromagnetic responses of AMB3.Their dependence on the generalised coordinates and currents in windings is suggested in this study to be obtained analytically by considering magnetic circuits with the use of equivalent circuits and the loop fluxes method.Forces Q qr are other generalised forces, in particular, the force of gravity, and H qr (t) are external time-dependent exciting forces and moments, in particular, caused by unbalance, whereas: For second-order nonlinear terms f′′ qr (x 1 ,...,z 3 ), required to account for the nonlinear features of the electromechanical system being considered, the following relations hold: Third-order nonlinear terms f′′′ qr (x 1 ,...,z 3 ) are not shown here due to their cumbersome notation; however, it is they that demonstrate the full co-relation between all generalised coordinates with the help of terms having no dependence on unbalance parameters.In characterising analytical model (9), it is necessary to note that, along with second-order differential equations of motion for generalised coordinates, it contains also first-order differential equations for flux linkages.However, the last two equations in (9) are present in the mathematical model because AMB takes part in the magnetic suspension of the rotor.The number of these equations equals that of the windings in all AMBs.The time constant of current change in the windings of electromagnets in these equations is accounted for by self and mutual induction coefficients (self-induction and mutual induction of circuits in electromagnets).The values that are inverted to these parameters are the coefficients with Ψ 2 c k /2 and Ψ c k Ψ c s /2 in the expression for magnetic energy W. The AMB damping properties are determined by the values of active resistances r c k and control system parameters accounted for in u c k .The expression for magnetic energy W has the magnetic resistances of sections of AMB magnetic core circuits and the resistances of air gaps whose formulae also include constants defined by their geometry and the values of the gaps proper, and hence, the generalised coordinates as well, with the dependence on them being nonlinear.Besides, magnetic energy has a square-law dependence on all flux linkages.Such an expression for magnetic energy couples equations of motion and equations for flux linkages into a unique system and determines the correlation of electric and magnetic dynamic phenomena with mechanical ones.
The expressions for magnetic energy W and the links between flux linkages Ψ c 1,2 and currents i c 1,2 (Fig. 3) can be found by considering equivalent circuits for the AMB3 electromagnetic circuit and schematisation of magnetic flux paths.The suggested method allows obtaining an expression for magnetic energy for any arbitrary complex and branched circuit of any radial or axial AMB.In the general case, such an approach accounts for both resistances of air gaps between stator poles and the rotor, and separate sections of magnetic cores and sections carrying dissipation fluxes [24].
Let us consider the schematic construction of AMB3 (Fig. 1, b) used in the laboratory setup and the equivalent circuit for magnetic circuits with account for dissipation x ,y ,y ,z ); r fluxes through slots filled with copper windings.They are shown in Fig. 4. The following notations are accepted: Φ k are magnetic fluxes in circuit sections; Φ Ck are loop fluxes; R k are magnetic resistances of circuit sections, where R pk are for the poles, R gk are for the air gaps between the poles and the rotor, R sk are for the sections of the stator yoke between two poles, R lk are for dissipation between the poles (slots filled with copper winding), and R ak are for the disk; E k =i c k w k are electromotive forces of coils.
The magnetic circuit is analysed using the loop fluxes method, which is an analog of the loop circuit one.The system of algebraic equations for loop fluxes is as follows: where The expressions for fluxes in circuit sections in terms of loop fluxes (the fluxes in nonadjacent branches are equal to loop fluxes if their directions coincide, and they are equal to loop fluxes with the opposite sign if they do not coincide; the fluxes in adjacent branches are determined in the same manner): The magnetic energy of the axial AMB (the energy of the entire magnetic circuit is the sum of energies of sections of this circuit): To determine the loop fluxes Φ C1 , Φ C2 , Φ C3 , Φ C4 as a function of currents, it is necessary to solve a system of all four equations (13).Next, transferring them into expressions for magnetic fluxes in the circuits (15), and these, in turn, into the expression for magnetic energy (16), yields the formula for W as a function of currents i c1 and i c2 and resistances of sections of the magnetic circuit (gaps), depending on generalised coordinates R g1 =R g1 (q r ), R g2 =R g2 (q r ) and R g3 =R g3 (q r ), R g4 =R g4 (q r ), where q r =(x 1 ,...,z 3 ).The system of the last two equations in (13) allows finding the dependencies of Φ C3 and Φ C4 on Φ C1 and Φ C2 because it is written namely for loop fluxes Φ C3 and Φ C4 .By using these expressions in (15), we obtain relationships: where, with account for ( 14), the coefficients depend on the generalised coordinates and are equal to: The magnetic fluxes in circuit sections with coils can be replaced with flux linkages if we assume that the magnetic flux pertinent to each turn of the coil is the same.Then the total or full magnetic flux through the coil circuit shall exceed this flux by the number of times equal to the number of coil turns: .
Accounting for this link of magnetic fluxes Φ 1 and Φ 2 in circuit sections with coils with flux linkages Ψ={Ψ с1 , Ψ с2 }, we obtain a final expression for magnetic energy in terms of flux linkage and generalised mechanical coordinates q r = ={x 1 , y 1 , x 2 , y 2 , z 3 }: Fig. 4. Design chart and the equivalent chart for substitution of magnetic circuits in an axial AMB where, with account for expressions ( 14) and ( 18), the coefficients with flux linkages Ψ depend on generalised coordinates q r and are equal to the following: In coefficients (21), the magnetic resistances of magnetic core sections R pk R sk , R lk and of the disk R ak are constant values depending on the geometric and physical parameters of AMB, and they can be found by using schematisation of magnetic flux paths [24].The dependence of coefficients (21) in the magnetic energy expression (20) on generalised coordinates q r is determined by the dependence thereon of the resistances of air gaps between the poles and the rotor R gk (x 1 , y 1 , x 2 , y 2 , z 3 ) based on the formulae ( 14) and (18).These resistances can be found through using expressions for magnetic conductances of air gaps under four poles of the axial AMB with account for the height of gaps depending on the rotor spatial position, i.e. on generalised mechanical coordinates: where g gk are conductances of air gaps under AMB poles, h z3± is gap height under the elementary area of poles located on the side of the positive and negative directions of axis z, respectively [24].
Analysis of results shows that, close to the rotor central position (±δ a /2), the values of gap conductance obtained by using simplified expressions (depending only on z 3 ) as compared to refined formulae (depending on x 1 , y 1 , x 2 , y 2 , z 3 ) differ by less than 0.5 %, and can be considered coinciding [24].However, the error of calculating these conductances as the real gap tends to zero can be as much as 20-40 % in some cases [24].This can result in falsehood in determining magnetic resistances and fluxes and, hence, the magnetic field energy and electromagnetic responses of the suspension.Such an approach introduces an error into the mathematical model (9), which cannot be neglected in some cases.Thus, for instance, failing to account for this fact when determining or checking the parameters of the control system and control algorithm with the help of numerical experiments can narrow the actual range of stability of rotor motion within AMB gaps [25].

Numerical and experimental research
Research was conducted for a rotor with a mass of 2.5 kg in a complete magnetic-electromagnetic suspension (Fig. 3), in which, for AMB3, a unique control method and an algorithm were used, i.e. formation of u c k in (9) [26]: u c2,1 = =(u max -2u min )z 3 2 /(2d a 2 )±u max z 3 /(2d a )+u min .The basic parameters have the following values: m=2.69 kg, l 1 =0.118 m, l 2 =0.166 m, J 1 =0.00997 kg×m 2 , J 3 =0.00347kg×m 2 , d r =5.5× ×10 -3 m, d r =3×10 -3 m, e=6×10 -5 m, g=0.003 rad, u max =24 V, and Q Rqr =b qr ¶q r / ¶t, where b qr =2.325 kg/s.A laboratory setup with such parameters was developed as a prototype of a complete magnetic suspension for an ECU rotor.It was used for experimental studying of possible nonlinear dynamic phenomena in the system when the angular rotational speed changes within 0 to 3'000 rpm.
The result of a series of experiments was the amplitude-frequency response (AFR) shown in Fig. 5.It allows evaluating the presence of resonant modes in the area being investigated and the kind of rotor motion corresponding to different rotational speeds.Thus, the following was found: -bifurcation of the first (~10.5 and ~12 Hz) and the second (~22.5 and ~33 Hz) resonances due to different PMB stiffness in the horizontal and vertical directions (anisotropy of bearings) due to different static equilibrium positions (x 1st =x 2st =0, y 1st and y 2st ≠0) with respect to centres of bearings that occur owing to the force of gravity; -direct (~10.5 Hz) and reverse (~12 Hz) cylindrical precessions as well as direct (~22.5 Hz) and reverse (~35 Hz) conical precessions (Fig. 5 shows vibration modes corresponding to these motions); -loss of vibrations with transition from one stable mode to another stable mode (the dashed area in Fig. 5).
Besides, our analysis of the results detected in the system concerned the following: harmonic vibrations with an excitation (rotational) frequency, subharmonic and superharmonic vibrations, multiple sub and super resonances, and a link between radial and axial vibrations.A detailed description of the results is given in the following, in comparison with the results of numerical modelling.During numerical modelling, the system of equations ( 9) was solved with the 4 th -5 th -order Runge-Kutta method for discrete angular speed values.The many-valuedness of the solution was checked and excluded by multiple computations for each frequency and different initial conditions.In doing so, stationary areas were searched for, whereas time intervals corresponding to transient processes were excluded from consideration.Hence, the results of the numerical analysis of forced vibrations are solutions for stationary areas and generalised coordinates x 1 , y 1 , x 2 , and y 2 in the angular speed range of 0-100π rad/s.They are shown in Fig. 6 as harmonics amplitudes A obtained by using the fast Fourier transform versus the driving force angular frequency ω 0 caused by the rotor's own unbalance.This frequency relates to the rotor angular speed as ω 0 =ω.The following notations are used in Fig. 6: A (1) is the first harmonic amplitude (Fig. 6), A (1/n) is the subharmonic amplitude (Fig. 6, a, c) and A (n) is the superharmonic amplitude (Fig. 6, b, d), where the number in parentheses is the multiplicity of the harmonic frequency with respect to the fundamental frequency ω 0 , with the dashed lines showing the skeleton curves.The natural frequencies of a nonrotating rotor that were calculated by using a linearized system of equations without account for damping [2p rad/s] are as follows: p 1x =10.39; p 1y =11.72; p 2x =25.03; p 2y =27.66.The curves in Fig. 6 show the dynamic behaviour of a rotor in the investigated range and, in essence, are the projections of 3-dimensional spectra on the coordinate planes OωA.Thereat, the dependence of the amplitude of the first harmonic of forced vibrations А (1) on the frequency of the harmonic driving force is the amplitude-frequency response, and the graphic representation of this dependence (Fig. 6) is the resonance curve.Analysis of these responses has shown that the first resonant mode (ω 1z ) corresponds to axial vibrations (Fig. 6).
Next, analysis of the results (Fig. 6) has demonstrated the following: superharmonic vibrations in the area of the second resonant mode (I); bifurcation of the second resonance due to anisotropy of PMB rigidity in the horizontal and vertical directions when at ω<ω 1x and ω>ω 1y the rotor's motion is of the direct cylindrical precession type, and in the range between these critical speeds ω 1x <ω<ω 1y the motion is of the reverse cylindrical precession type (II); super-resonant vibrations ω 2x (2) , which coincide also with the inner resonance ω 2x(2) =ω 1y (III); bifurcation of the third resonance due to anisotropy of PMB rigidity when at ω<ω 2x and ω>ω 2y rotor motion is of the direct cylindrical precession type, and in the range between these critical speeds (ω 2x <ω<ω 2y ) when the motion is of the reverse conical precession type (IV); subresonant vibrations ω 1y(1/4) , which are enhanced by inner resonance ω 1y(1/4) =4ω 1y =4ω 2x(2) (Fig. 6, a, b), with these subharmonic vibrations occurring at relatively high excitation frequencies, and their amplitudes significantly exceeding the amplitudes of the first harmonic (V); the form of resonance curves in the area of the third resonant mode (ω 2x and ω 2y ) is specific to systems with rigid characteristics of the restoring force, which is true for PMB (VI); the third resonant mode is more dangerous than the second one because it is accompanied by a significant amplitude increase as during motion of the conical precession type (angular vibrations) the flatness of gaps in the axial AMB is disturbed, resulting in a moment coinciding in the direction with the angular deflection of the rotor (VII); in the area of frequencies, wherein two stable forced vibration modes with two different amplitudes are possible, a failure of vibrations is observed (VIII); the fundamental and the superharmonic resonant vibrations in the axial direction are excited by a load acting in the radial direction (by the rotor's own unbalance), with the peaks of super-resonant axial vibrations coinciding with the peaks of the fundamental radial vibrations (Fig. 6, a-d), which is the result of accounting for the interrelation between radial and axial generalised coordinates with nonlinear terms in the equations of motion (9) (IX).The same resonant modes and phenomena were found in the system also during experimental research.The adequacy of the mathematical model representing a system of nonlinear completely mutually coupled by generalised mechanical coordinates x 1 ,...,z 3 and flux linkages Ψ c1 and Ψ c2 equations (with account for the control law, i. e. voltages u 1,2 also dependent on x 1 ,...,z 3 ) can be judged by the results of comparing the calculated data (Fig. 6) with experimentally obtained amplitude-frequency responses (Fig. 5) and dependencies of harmonic amplitudes that differ from the fundamental one by the driving force frequency.Thus, comparative analysis of the results has shown an identity for both qualitative representation of processes in the system and quantitative determination of their parameters: for the amplitude, the difference is within 2-3 %, and for the values of resonant frequencies, the difference is within 0.2-0.5 %.

Discussion of the findings on the dynamics of a rotor in magnetic bearings
The result of this research, which is a follow-up of a large variety of studies, is the development and practical implementation of the method of mathematical description of linear and nonlinear rotor dynamics phenomena in systems with magnetic bearings of different types, which affect the vibration activity of power rotor machinery.The advantage of the approach suggested is that the interrelationship of electric, magnetic and mechanical stationary and nonstationary processes can be accounted for as shown by the example of a laboratory setup.Applying this approach to modelling the dynamics of rotor systems with magnetic bearings improves the dynamic parameters of a whole class of rotor machinery due to a more correct description of dynamic processes and phenomena occurring therein.In turn, this has the effect of cutting the cost of development activities at the design and commissioning stages, and reducing operational and power resources costs.

Conclusions
1.It is known that analysis based on linearized models allows judging only the stability of equilibrium states with small deflections.The negligible nonlinear equation terms in this case, when investigating motion with increasing deflection, allow expanding the information content of the mathematical model about nonlinear effects occurring in the system.Therefore, the study has suggested an analytical model that accounts for the nonlinear interrelation of mechanical and electromagnetic processes in the "rotor in magnetic bearings" system.Its adequacy and feasibility for studying the dynamics of rigid rotors in magnetic bearings of different types has been proved experimentally.
2. The main distinctive feature of the model implemented, among other things, by using the analytical method of analysing electromagnetic circuits of active magnetic bearings, is being able to account for the following: -a nonlinear dependence of magnetic forces on gaps between movable and stationary parts in PMBs and AMBs, and on currents in the windings of AMB electromagnets; -a current delay in the windings of AMB electromagnets, i. e. nonlinearity linked to the inductance of the coils; -geometric links between the electromagnets of one AMB and links between all AMB of one rotor, which results in coupling of processes in orthogonal directions; -practically any AMB control law; -limitations on the control current caused by physical constraints in the control system controller; -delays in the control system controller; -dissipation fluxes as well as magnetic resistances of AMB magnetic core sections, making the mathematical model insensitive as regards origination of "nonzero" gaps and currents.
3. The developed imitation of the mathematical model built around the suggested mathematical modelling method has been used for numerical research into a complete magnetic suspension in a laboratory setup.The study has shown how the model can be used for investigating the mechanisms of excitation of spatial vibrations in rotating rigid rotors in an MB, finding out the conditions of existence of different resonant modes (including super, sub and inner resonances), super-and subharmonic, and combined vibrations, and also testing control algorithms and selecting optimal suspension parameters.

Fig. 2 .+e 2 2 )
Fig. 2. Schematic diagram of a complete rotor suspension in an MB Here, a spatial fixed orthogonal system of coordinates O * xyz has been introduced.Its axis O * z passes through the centres of radial MBs.Point C is centre of the mass, m is the rotor mass, J x , J y , and J z are the main central moments of inertia; e x =e 1 and e y =e 2 (e 2 =e 1 2 +e 2 2 ) are linear whereas γ x =γ 1 and γ y =γ 2 are angular unbalance parameters.Points O, O 1 , O 2 , and O 3 are on the rotor axis of stiffness, with O being in one plane with C (perpendicular to the axis); O 1 and O 2 are centres of the thrust areas of radial MBs (with l 1 +l 2 =l); and O 3 is that of the axial MB.Gaps in MB1, MB2, and MB3

Fig. 3 .
Fig. 3. Complete suspension for a laboratory setup rotor in an MB: a -schematic diagram; b -restoring magnetic forces in PMB1 and 2 vs. rotor displacements