# Research on Vibration Characteristic Analysis of Magnetic Suspended Dual-Rotor System with Fixed-Point Rubbing

Authors: Dongxiong Wang; Tao Yi; Xuan Qin; Ju Zhang; Changcheng Yin; Baohua Wang
Journal: Shock and Vibration
 (2022)

---

## Abstract

Due to the elimination of wear and lubrication system and effective solution to the vibration control issue, magnetic suspended dual-rotor system (MSDS) offers the possibility to significantly enhance the performance of aero-engines. However, research on MSDS vibration characteristics under fixed-point rubbing with the consideration of nonlinear support characteristics of active magnetic bearing (AMB) has rarely been addressed. In this research, an improved EF model of AMB is proposed by equivalent magnetic circuit method and verified by finite element model. Then, the Lankarani-Nikravesh model is employed to describe the impact force during rubbing process. On this basis, the rubbing dynamic model is established by finite element method and verified by comparing with the experimental results in literature. Finally, the MSDS vibration characteristics are explored, and some interesting phenomena are found. Fixed-point rubbing will cause abundant frequency components equally spaced around the combined frequency components of rotor rotational speeds at different frequency intervals. Due to the stimulation of the AMB nonlinear EF, the fractional frequency component of MSDS second-order critical speed emerges in a certain speed range and implies potential instability of the system, and the rotational speed corresponding to the initial appearance of this frequency component will decrease when fixed-point rubbing occurs. With the decrease of proportional coefficient and the increase of differential coefficient, the maximum normal impact force in the rubbing process decreases.


---

## Body

## 1. Introduction

Active magnetic bearing (AMB) is a new type of high-performance bearing, which can achieve stable rotor suspension through controllable electromagnetic force (EF). Due to the advantages of noncontact, high speed, high reliability, and adjustable support characteristics over traditional mechanical bearings [1], the application of AMB in high-speed rotating machinery has received high attention in recent years. In order to obtain high efficiency and high thrust weight ratio, the dual-rotor structure is widely employed as a core rotor component in aero-engine. The employment of AMB instead of mechanical bearings to support the dual-rotor system, namely, magnetic suspended dual-rotor system (MSDS), can significantly reduce the system complexity and weight, optimize the structure, obtain higher reliability and maintainability, and greatly improve the overall performance of the dual-rotor system.With the continuous pursuit of high thrust weight ratio and high power density in modern aero-engine, the reserved gaps between rotating and static components are forced to be tighter, causing greater susceptibility for rubbing faults. As a secondary failure caused by failures such as imbalance, misalignment, loose bearings, thermal expansion, oil whirl, and cracks, rub-impact is of great detriment. Weak rub-impact will cause aggravated structural vibration, seal wear, and increased noise, thus shortening the service life of machinery, while serious rub-impact tends to result in tragic accidents [2]. Therefore, investigation of the MSDS rub-impact characteristics is of vital importance to ensure the stable operation of the aero-engine.Rub-impact is a very common and typical type of fault in aero-engine. According to the contact area, rubbing can be divided into fixed-point rubbing [3], partial rubbing [4], and full annular rubbing [5]. The dynamic characteristics of single rotor system with different rubbing forms in the aspects of theory [6, 7], simulation [8, 9], and experiment [10, 11] have been studied thoroughly, while the dynamic characteristics of single rotor systems are significantly different from those of dual-rotor systems. The introduction of interbearings becomes a cross-excitation source in the dual-rotor system, which intensifies the vibration coupling between the inner and outer rotors, rendering the dynamic characteristics of the dual-rotor system more complex. Taking gyroscopic effects and the coupling effects among inner and outer rotors and casing into consideration, Ma et al. [12] established the dynamic model of a dual-rotor-casing system and studied the effects of rotational speed and speed ratio on rubbing vibration characteristics of dual-rotor-casing system. Considering the elastic deformation and motion of the casing, Haiqin et al. [13] established a dynamic model of a dual-rotor system with rubbing failure and analyzed the rubbing characteristics of the system with parameters such as speed and unbalance and contact stiffness. The results showed that the doubling frequency and combined frequency of the rotors would appear in the spectrum when the rub-impact fault occurred. Based on the coupled bending and torsional dynamic equations of a dual-rotor system with rub-impact, Zhang et al. [14] analyzed the spectral characteristics and bifurcation characteristics of the bending vibration and torsion vibration and observed that the torsional vibration exhibits characteristic frequencies similar to bending vibration, including multiple frequency, fractional frequency, and the combined frequency of working frequency and multiple or fractional frequency. With the softening characteristics of coatings applied on the surface of disks and fixed limiter taken into consideration, a dual-rotor dynamic model with fixed-point rubbing was established by Yang et al. [15, 16], where the Lankarani-Nikravesh model was adopted to describe the impact force. On this basis, the effects of parameters such as speed ratio, initial clearance, and radius of curvature of the limiter on the rubbing response were inspected. Ma et al. [17] investigated the nonlinear dynamic behavior of a dual-rotor system supported by rolling bearings and squeeze film damper (SFD), explored the steady-state vibration characteristics with different support parameters, and found that the increase of ball bearing clearance leads to the strengthening of resonance hysteresis phenomenon strengthening, and small clearances of ball bearing and oil film contribute to the improvement of system stability. Wang et al. [18] built a finite element model of a dual-rotor blade-casing system and investigated the effects of speed ratio and unbalance and contact stiffness on rubbing characteristics based on the casing acceleration response. The results showed that the impact frequency was equal to the product of the rotational frequency and the numbers of blades, and the fraction frequency components of rotating speed difference were excited on both sides of impact frequency and its multiple frequency components.From above literature, it can be seen that researches on rubbing dynamics of dual-rotor system supported by mechanical bearings are mainly focused on the dynamic modeling, rubbing characteristics analysis, stability, and nonlinear dynamic behavior, and the rubbing responses of dual-rotor system are complicated and significantly affected by the support characteristics of bearings. In the MSDS, the AMB support characteristics are not only related to the structural parameters, but also determined by the control system and rotor vibration conditions. The rotor dynamic responses affect the support characteristics of AMB through the mechanical system and the control system, which change the system rubbing responses. This coupling effect will make the dynamic responses of the dual-rotor system more complicated. Considering the contact effect between auxiliary bearing and inner rotor, Ebrahimi et al. [19, 20] built a ten-degree-of-freedom dynamic model of MSDS control based on Lagrange equation and discussed the nonlinear dynamic behaviors with different speed ratio, gravity, inner rotor stiffness, and interbearing stiffness. However, the possible rubbing fault was neglected between dual-rotor system and casing. Besides, in the dynamic modeling of magnetic suspended rotor system, the conventional nonlinear EF model of AMB adopted in present researches only considers the influence of the rotor displacement in a certain direction on the EF in the same direction [19–22], while it ignores the influence of rotor displacement component in the orthogonal direction of EF, which reduces the accuracy of dynamic model. Moreover, there are few publications that report the rubbing responses of MSDS. Nevertheless, researches on the rubbing dynamics are the foundation of reasonable design of AMB control system, and they are of considerable importance for the stable operation and fault diagnosis of MSDS in aero-engine.Various coatings are usually applied on the surface of key components in practical aero-engine to extend the service life [23]. For rub-impact between rotating and stationary components with soften characteristics of coatings, the impact stiffness depends on the local contact stiffness of the coatings rather than the structural stiffness. In this case, Cao et al. [24] proposed that the Lankarani-Nikravesh model [25] was more suitable to describe rub-impact than the piecewise linear model [26]. In this research, aiming at large modeling error of the conventional EF due to the ignorance of rotor displacement in the orthogonal direction, an improved AMB EF model is proposed by equivalent magnetic circuit method and verified by finite element model. Considering the softening characteristics of coatings printed on the rubbing components, the Lankarani-Nikravesh model and the Coulomb model are employed to, respectively, describe the impact force and friction force during the rubbing process. The fixed-point rubbing dynamic model of the MSDS is established by finite element method and verified by experimental results in literature. On this basis, the effects of AMB control parameters on the MSDS dynamic responses are investigated.


## 2. Mathematical Formulation

For the MSDS structure shown in Figure1, the inner rotor and the outer rotor are connected through two interbearings, and the outer rotor is supported by two AMBs. The left end of the outer rotor and the right end of the inner rotor are, respectively, connected to two high-speed motors through elastic diaphragm couplings, allowing a certain axial offset and angular deflection, so that the axial freedom of the inner and outer rotors can be adjusted. For the convenience of analysis, the following assumptions are made about the dynamic modeling process of the MSDS with fixed-point rubbing fault.(1)
Only the radial vibration is considered, and the effects of torsional and axial vibrations are ignored.(2)
The interbearing is simplified as a linear elastic support.(3)
The AMB rotor is thermally sleeved on the outer rotor shaft, treated as a whole with the rotating shaft.(4)
Rubbing occurs between the disk 3 on the outer rotor and the fixed limiter, both of which are covered by coatings. The local contact stiffness is much smaller than its respective structural stiffness.(5)
The influence of thermal effect and friction torque during the rubbing process is ignored.(6)
The imbalances of the system are located on disk 1 and disk 3.Figure 1
The MSDS structure. ①⑧ Interbearing. ②⑦ AMB rotor. ③ Disk 2. ④ Outer rotor. ⑤ Inner rotor. ⑥ Disk 3. ⑨Disk 1.


### 2.1. Electromagnetic Force Model of AMB

The conventional used AMB nonlinear EFFxc in the horizontal direction (x direction) is expressed as [17–20](1)Fxc=μ0N2Apcosα4Ib+ixg0−x2−Ib−ixg0+x2,where μ0 is permeability of vacuum, N is the coil turn, Ap is the sectional area of the stator pole, α is the angle between the central lines of the pole and pole-pair, g0 is the normal air gap length, Ib is the bias current of coils, ix is the control current, and x is the rotor displacement in the horizontal direction. The EF Fy in the vertical direction can be easily obtained by replacing x with y in (1). It is apparent that the rotor displacement in x direction only affects the EF in this direction.For a typical eight-pole AMB, the configuration in the horizontal direction is depicted in Figure2(a). In the case of ignoring flux leakage, magnetic saturation, and magnetic coupling, the equivalent magnetic circuit for the left-hand side pole-pair in Figure 2(a) is shown in Figure 2(b), where F1 and F2 are the magnetomotive forces of the coils, and Rg1 and Rg2 are the air gap reluctances of the pole-pair.Figure 2
AMB configuration and equivalent magnetic circuit for one pole-pair. (a) AMB configuration in horizontal direction. (b) Equivalent magnetic circuit.
(a)(b)According to the equivalent magnetic circuit, the total reluctanceR for one pole-pair is(2)R=Rg1+Rg2=g1+g2μ0Ap,where g1 and g2 are the air gap lengths between stator poles and rotor. The magnetomotive forces are(3)F1=F2=NIb+ix2.Therefore, the magnetic flux Φ in the air gap is(4)Φ=F1+F2R.The EF produced by one pole-pair is expressed as(5)F=Φ2cosαμ0Ap.With differential drive mode, the nonlinear EF of the AMB in the horizontal direction is obtained as(6)Fxe=μ0N2ApcosαIb+ixg3+g42−Ib−ixg1+g22,where g3 and g4 are the air gap lengths between the right-hand side pole-pair and rotor. When the center of AMB rotor moves to the point (x, y), gi (i = 1, 2, 3, 4) are, respectively, given by(7)g1=g0+xcosα+ysinα,g2=g0+xcosα−ysinα,g3=g0−xcosα−ysinα,g4=g0−xcosα+ysinα.By the same argument, the AMB EF in the vertical direction can be obtained. Comparing the conventional EF model (CEFM) in equation (1) and the improved EF model (IEFM) in equation (6), it is obvious that the conventional EF is only a function of displacement x, while the improved EF is a nonlinear function of both displacements x and y. Thus, the effect of components of rotor arbitrary displacement in the horizontal and vertical directions on the EF is taken into consideration simultaneously in the IEFM.In order to evaluate the effectiveness of IEMF, the finite element model (FEM) with high precision is taken as the reference, and the results of IEFM and CEFM are compared with FEM. For the facilitation of analysis, the rotor eccentric ratioe is defined as the ratio of the rotor vibration amplitude to the nominal air gap length, and the eccentric angle β is defined as β = arctan(x/y). With g0 = 0.5 mm, N = 172, Ap = 80.48 mm2, Ib = 3A, α = 22.5°, e = 10% and β = 60°, Figure 3 shows the EF curves and relative error curves at different control current. With other parameters invariant, Figure 4 depicts the EF curves and relative error curves at different eccentric angle when Ix = 1A. Due to the consideration of the effect of rotor displacement components in both horizontal and vertical directions, it is apparent from Figures 3 and 4 that the calculation accuracy of IEMF is obviously higher than that of CEMF at different control current and eccentric angle. Therefore, the IEMF is adopted to establish the rubbing dynamic model of the MSDS in this research.Figure 3
Result comparison between IEF and CEF with different control current. (a) Electromagnetic force. (b) Relative error.
(a)(b)Figure 4
Result comparison between IEF and CEF with different eccentric angles. (a) Electromagnetic force. (b) Relative error.
(a)(b)


### 2.2. Fixed-Point Rubbing Model

Considering the softening characteristics of the coatings applied on the surface of both disk and fixed limiter, the Lankarani-Nikravesh model suitable for describing the local deformation between the interfaces is employed to depict the impact force. Meanwhile, the Coulomb model is used to describe the friction force during the contact process. The schematic diagram of the rub-impact model is shown in Figure5. The normal impact force Fn and tangential friction force Ft can be expressed as(8)Fn=krr1.51+31−ce2v4v0,Ft=frFn,where ce is the restitution coefficient, fr is the friction coefficient, r is the distance between the disk and limiter end in the collision direction, v is the impact velocity of the disk, v0 is the initial impact velocity during each rubbing process, and kr is the contact stiffness expressed as(9)kr=431−μd2/Ed+1−μp2/EdRdRpRd+Rp0.5,where μd, μp, Ed, Ep, Rd and Rp are Poisson’s ratio, elastic modulus, and radii of the disk and limiter, respectively.Figure 5
Schematic diagram of rub-impact force between the disk and fixed limiter.If the initial clearance between the disk and the limiter isδ0, the angle of the vertical direction and the line OcOF between the disk center Oc(x, y) and the limiter end center OF(0, Rd + Rp + δ0) is θ, and then OcOF can be expressed as(10)OcOF=Rd+Rp+δ0−y2+x2.Thus, the invading lengthr is given by(11)r=Rd+Rp−OcOF=Rd+Rp−Rd+Rp+δ0−y2+x2,r > 0 implies that the collision occurs. The normal impact angle θ during rubbing is(12)θ=arctanAOcAOF=arctanxRd+Rp+δ0−y.If the whirl speed vector of the disk center isOcB⟶=vx,vy, and the angle α between OcB and the vertical direction is(13)α=arctanvxvy.then it can be seen from Figure3 that the normal impact velocity of the disk in the collision direction is(14)v=vx2+vy2cosθ+α.On the basis of knowingr, v and v0, the impact force Fn and the friction force Ft during the rubbing can be obtained by equation (8), and their projections Frx and Fry in the horizontal and vertical directions are, respectively, expressed as(15)Frx=Fnsinθ+Ftcosθ,Fry=−Fncosθ+Ftsinθ.It is obvious from equation (8) that the initial impact velocity v0 will directly affect the magnitude of Fn, so it is very important to accurately determine the initial moment of rubbing occurrence and precisely calculate the value of the initial impact velocity. In this work, the linear interpolation method [27] is adopted to modify the integration step to improve the calculation accuracy.


### 2.3. Dynamic Model of MSDS

According to the MSDS structural characteristics shown in Figure1, the shafts are described by using the Euler-Bernoulli beam element considering moment of inertia and gyroscopic effect. The inner rotor is discretized into 11 beam elements and 12 nodes, and the outer rotor is discretized into 22 beam elements and 23 nodes. For convenient modeling, the disks are regarded as rigid. As a result, the finite element model of MSDS is depicted in Figure 6.Figure 6
Finite element model of the MSDS.According to the rotor dynamics theory and finite element method, the equations of motion of the inner rotor are(16)Miq¨1i+ωiJiq˙2i+Kiq1i=Q1i,Miq¨2i−ωiJiq˙1i+Kiq2i=Q2i,where Mi, Ji and Ki are mass matrix, gyroscopic matrix, and stiffness matrix of the inner rotor, respectively. ωi is the inner rotor rotational speed. q1i and q2i are the generalized displacement vectors of the inner rotor, and they are expressed as(17)q1i=x1ϕy1x2ϕy2...x12ϕy12T,q2i=y1−ϕx1y2−ϕx2...y12−ϕx12T,where xk and yk (k = 1, 2, …, 12) represent the translational displacement of the nodes 1–12, and ϕyk and ϕxk (k = 1, 2, …, 12) represent corresponding angular displacement. Q1i and Q2i are generalized force vectors and are expressed as(18)Q1i=000...ud1ωi2cosωit00000T,Q2i=000...ud1ωi2sinωit00000T,where ud1 is the imbalance on the disk 1, and t is the time. It should be noted that the position number value of non-zero elements in Q1i and Q2i is 19.Similarly, the equations of motion of the outer rotor are(19)Moq¨1o+ωoJoq˙2o+Koq1o=Q1o,Moq¨2o−ωoJoq˙1o+Koq2o=Q2o,where Mo, Jo and Ko are mass matrix, gyroscopic matrix, and stiffness matrix of the outer rotor, respectively. ωo is the outer rotor rotational speed. q1o and q2o are described as(20)q1o=x13ϕy13x14ϕy14...x35ϕy35T,q2o=y13−ϕx13y14−ϕx14...y35−ϕx35T,Q1o and Q2o are expressed as(21)Q1o=00...FmxL0...ud3ωo2cosωot+Frx0...FmxR0...00T,Q2o=00...FmyL0...ud3ωo2sinωot+Fry0...FmyR0...00T,where ud3 is the imbalance on the disk 3. FmxL, FmyL, FmxR and FmyR are EFs in the x and y directions of the AMBs at the left and right hand sides, respectively. The position numbers of the three non-zero terms in Q1o are, respectively, 13, 27, and 35, and the same for Q2o.Based on the above analysis, the MSDS dynamic equations can be expressed as(22)Mq¨+Cq˙+Kq=Q,where q, M, C, and K are the overall generalized displacement vector, overall mass matrix, overall damping matrix, and overall stiffness matrix of the MSDS, and they are, respectively, expressed as(23)q=q1iq1oq2iq2o,(24)M=diagMiMoMiMo,(25)C=0diagωiJiωoJo−diagωiJiωoJo0,(26)K=diagKiKoKiKo.In the above modeling process, the effects of the support stiffness and support damping of interbearings have not been considered yet. Therefore, the motion equations of MSDS (equations (22)–(26)) need to be modified subsequently. According to the node numbers of interbearings shown in Figure 6, the stiffness matrix and damping matrix need to be separately corrected, and the specific correction process of the motion equation is as follows.For the interbearing between node 3 of the inner rotor and node 14 of the outer rotor, its stiffness is added toK(5, 5), K(27, 27), K(75, 75), and K(97, 97), and its damping is added to C(5, 5), C(27, 27), C(75, 75), and C(97, 97). In addition, its stiffness should be subtracted from K(5, 27), K(27, 5), K(75, 97), and K(97, 75), and its damping should be subtracted from C(5, 27), C(27, 5), C(75, 97), and C(97, 75). Similarly, the stiffness is added to K(15, 15), K(67, 67), K(85, 85), and K(137, 137), and the damping is added to C(15, 15), C(67, 67), C(85, 85), and C(137, 137) for the interbearing between node 8 of the inner rotor and node 34 of the outer rotor. Besides, its stiffness should be subtracted from K(15, 67), K(67, 15), K(85, 137), and K(137, 85), and its damping should be subtracted from C(15, 67), C(67, 15), C(85, 137), and C(137, 85).On this basis, the Newmark-β method is adopted to obtain the dynamic response of MSDS with fixed-point rubbing fault.


## 2.1. Electromagnetic Force Model of AMB

The conventional used AMB nonlinear EFFxc in the horizontal direction (x direction) is expressed as [17–20](1)Fxc=μ0N2Apcosα4Ib+ixg0−x2−Ib−ixg0+x2,where μ0 is permeability of vacuum, N is the coil turn, Ap is the sectional area of the stator pole, α is the angle between the central lines of the pole and pole-pair, g0 is the normal air gap length, Ib is the bias current of coils, ix is the control current, and x is the rotor displacement in the horizontal direction. The EF Fy in the vertical direction can be easily obtained by replacing x with y in (1). It is apparent that the rotor displacement in x direction only affects the EF in this direction.For a typical eight-pole AMB, the configuration in the horizontal direction is depicted in Figure2(a). In the case of ignoring flux leakage, magnetic saturation, and magnetic coupling, the equivalent magnetic circuit for the left-hand side pole-pair in Figure 2(a) is shown in Figure 2(b), where F1 and F2 are the magnetomotive forces of the coils, and Rg1 and Rg2 are the air gap reluctances of the pole-pair.Figure 2
AMB configuration and equivalent magnetic circuit for one pole-pair. (a) AMB configuration in horizontal direction. (b) Equivalent magnetic circuit.
(a)(b)According to the equivalent magnetic circuit, the total reluctanceR for one pole-pair is(2)R=Rg1+Rg2=g1+g2μ0Ap,where g1 and g2 are the air gap lengths between stator poles and rotor. The magnetomotive forces are(3)F1=F2=NIb+ix2.Therefore, the magnetic flux Φ in the air gap is(4)Φ=F1+F2R.The EF produced by one pole-pair is expressed as(5)F=Φ2cosαμ0Ap.With differential drive mode, the nonlinear EF of the AMB in the horizontal direction is obtained as(6)Fxe=μ0N2ApcosαIb+ixg3+g42−Ib−ixg1+g22,where g3 and g4 are the air gap lengths between the right-hand side pole-pair and rotor. When the center of AMB rotor moves to the point (x, y), gi (i = 1, 2, 3, 4) are, respectively, given by(7)g1=g0+xcosα+ysinα,g2=g0+xcosα−ysinα,g3=g0−xcosα−ysinα,g4=g0−xcosα+ysinα.By the same argument, the AMB EF in the vertical direction can be obtained. Comparing the conventional EF model (CEFM) in equation (1) and the improved EF model (IEFM) in equation (6), it is obvious that the conventional EF is only a function of displacement x, while the improved EF is a nonlinear function of both displacements x and y. Thus, the effect of components of rotor arbitrary displacement in the horizontal and vertical directions on the EF is taken into consideration simultaneously in the IEFM.In order to evaluate the effectiveness of IEMF, the finite element model (FEM) with high precision is taken as the reference, and the results of IEFM and CEFM are compared with FEM. For the facilitation of analysis, the rotor eccentric ratioe is defined as the ratio of the rotor vibration amplitude to the nominal air gap length, and the eccentric angle β is defined as β = arctan(x/y). With g0 = 0.5 mm, N = 172, Ap = 80.48 mm2, Ib = 3A, α = 22.5°, e = 10% and β = 60°, Figure 3 shows the EF curves and relative error curves at different control current. With other parameters invariant, Figure 4 depicts the EF curves and relative error curves at different eccentric angle when Ix = 1A. Due to the consideration of the effect of rotor displacement components in both horizontal and vertical directions, it is apparent from Figures 3 and 4 that the calculation accuracy of IEMF is obviously higher than that of CEMF at different control current and eccentric angle. Therefore, the IEMF is adopted to establish the rubbing dynamic model of the MSDS in this research.Figure 3
Result comparison between IEF and CEF with different control current. (a) Electromagnetic force. (b) Relative error.
(a)(b)Figure 4
Result comparison between IEF and CEF with different eccentric angles. (a) Electromagnetic force. (b) Relative error.
(a)(b)


## 2.2. Fixed-Point Rubbing Model

Considering the softening characteristics of the coatings applied on the surface of both disk and fixed limiter, the Lankarani-Nikravesh model suitable for describing the local deformation between the interfaces is employed to depict the impact force. Meanwhile, the Coulomb model is used to describe the friction force during the contact process. The schematic diagram of the rub-impact model is shown in Figure5. The normal impact force Fn and tangential friction force Ft can be expressed as(8)Fn=krr1.51+31−ce2v4v0,Ft=frFn,where ce is the restitution coefficient, fr is the friction coefficient, r is the distance between the disk and limiter end in the collision direction, v is the impact velocity of the disk, v0 is the initial impact velocity during each rubbing process, and kr is the contact stiffness expressed as(9)kr=431−μd2/Ed+1−μp2/EdRdRpRd+Rp0.5,where μd, μp, Ed, Ep, Rd and Rp are Poisson’s ratio, elastic modulus, and radii of the disk and limiter, respectively.Figure 5
Schematic diagram of rub-impact force between the disk and fixed limiter.If the initial clearance between the disk and the limiter isδ0, the angle of the vertical direction and the line OcOF between the disk center Oc(x, y) and the limiter end center OF(0, Rd + Rp + δ0) is θ, and then OcOF can be expressed as(10)OcOF=Rd+Rp+δ0−y2+x2.Thus, the invading lengthr is given by(11)r=Rd+Rp−OcOF=Rd+Rp−Rd+Rp+δ0−y2+x2,r > 0 implies that the collision occurs. The normal impact angle θ during rubbing is(12)θ=arctanAOcAOF=arctanxRd+Rp+δ0−y.If the whirl speed vector of the disk center isOcB⟶=vx,vy, and the angle α between OcB and the vertical direction is(13)α=arctanvxvy.then it can be seen from Figure3 that the normal impact velocity of the disk in the collision direction is(14)v=vx2+vy2cosθ+α.On the basis of knowingr, v and v0, the impact force Fn and the friction force Ft during the rubbing can be obtained by equation (8), and their projections Frx and Fry in the horizontal and vertical directions are, respectively, expressed as(15)Frx=Fnsinθ+Ftcosθ,Fry=−Fncosθ+Ftsinθ.It is obvious from equation (8) that the initial impact velocity v0 will directly affect the magnitude of Fn, so it is very important to accurately determine the initial moment of rubbing occurrence and precisely calculate the value of the initial impact velocity. In this work, the linear interpolation method [27] is adopted to modify the integration step to improve the calculation accuracy.


## 2.3. Dynamic Model of MSDS

According to the MSDS structural characteristics shown in Figure1, the shafts are described by using the Euler-Bernoulli beam element considering moment of inertia and gyroscopic effect. The inner rotor is discretized into 11 beam elements and 12 nodes, and the outer rotor is discretized into 22 beam elements and 23 nodes. For convenient modeling, the disks are regarded as rigid. As a result, the finite element model of MSDS is depicted in Figure 6.Figure 6
Finite element model of the MSDS.According to the rotor dynamics theory and finite element method, the equations of motion of the inner rotor are(16)Miq¨1i+ωiJiq˙2i+Kiq1i=Q1i,Miq¨2i−ωiJiq˙1i+Kiq2i=Q2i,where Mi, Ji and Ki are mass matrix, gyroscopic matrix, and stiffness matrix of the inner rotor, respectively. ωi is the inner rotor rotational speed. q1i and q2i are the generalized displacement vectors of the inner rotor, and they are expressed as(17)q1i=x1ϕy1x2ϕy2...x12ϕy12T,q2i=y1−ϕx1y2−ϕx2...y12−ϕx12T,where xk and yk (k = 1, 2, …, 12) represent the translational displacement of the nodes 1–12, and ϕyk and ϕxk (k = 1, 2, …, 12) represent corresponding angular displacement. Q1i and Q2i are generalized force vectors and are expressed as(18)Q1i=000...ud1ωi2cosωit00000T,Q2i=000...ud1ωi2sinωit00000T,where ud1 is the imbalance on the disk 1, and t is the time. It should be noted that the position number value of non-zero elements in Q1i and Q2i is 19.Similarly, the equations of motion of the outer rotor are(19)Moq¨1o+ωoJoq˙2o+Koq1o=Q1o,Moq¨2o−ωoJoq˙1o+Koq2o=Q2o,where Mo, Jo and Ko are mass matrix, gyroscopic matrix, and stiffness matrix of the outer rotor, respectively. ωo is the outer rotor rotational speed. q1o and q2o are described as(20)q1o=x13ϕy13x14ϕy14...x35ϕy35T,q2o=y13−ϕx13y14−ϕx14...y35−ϕx35T,Q1o and Q2o are expressed as(21)Q1o=00...FmxL0...ud3ωo2cosωot+Frx0...FmxR0...00T,Q2o=00...FmyL0...ud3ωo2sinωot+Fry0...FmyR0...00T,where ud3 is the imbalance on the disk 3. FmxL, FmyL, FmxR and FmyR are EFs in the x and y directions of the AMBs at the left and right hand sides, respectively. The position numbers of the three non-zero terms in Q1o are, respectively, 13, 27, and 35, and the same for Q2o.Based on the above analysis, the MSDS dynamic equations can be expressed as(22)Mq¨+Cq˙+Kq=Q,where q, M, C, and K are the overall generalized displacement vector, overall mass matrix, overall damping matrix, and overall stiffness matrix of the MSDS, and they are, respectively, expressed as(23)q=q1iq1oq2iq2o,(24)M=diagMiMoMiMo,(25)C=0diagωiJiωoJo−diagωiJiωoJo0,(26)K=diagKiKoKiKo.In the above modeling process, the effects of the support stiffness and support damping of interbearings have not been considered yet. Therefore, the motion equations of MSDS (equations (22)–(26)) need to be modified subsequently. According to the node numbers of interbearings shown in Figure 6, the stiffness matrix and damping matrix need to be separately corrected, and the specific correction process of the motion equation is as follows.For the interbearing between node 3 of the inner rotor and node 14 of the outer rotor, its stiffness is added toK(5, 5), K(27, 27), K(75, 75), and K(97, 97), and its damping is added to C(5, 5), C(27, 27), C(75, 75), and C(97, 97). In addition, its stiffness should be subtracted from K(5, 27), K(27, 5), K(75, 97), and K(97, 75), and its damping should be subtracted from C(5, 27), C(27, 5), C(75, 97), and C(97, 75). Similarly, the stiffness is added to K(15, 15), K(67, 67), K(85, 85), and K(137, 137), and the damping is added to C(15, 15), C(67, 67), C(85, 85), and C(137, 137) for the interbearing between node 8 of the inner rotor and node 34 of the outer rotor. Besides, its stiffness should be subtracted from K(15, 67), K(67, 15), K(85, 137), and K(137, 85), and its damping should be subtracted from C(15, 67), C(67, 15), C(85, 137), and C(137, 85).On this basis, the Newmark-β method is adopted to obtain the dynamic response of MSDS with fixed-point rubbing fault.


## 3. Validation of the Dynamic Model

In this section, the dynamic model of the MSDS with fixed-point rubbing is verified by the experimental results in literature. When the external load disturbance is small, the rotor only vibrates in a small range near the equilibrium position. Under this situation, the nonlinear support characteristics of AMB can be linearized at the equilibrium position and modeled as the equivalent stiffness and equivalent damping model [28, 29]. For the nonlinear EF model of the AMB in Section 2.1, ignoring the quadratic and higher order terms, the Taylor expansion of equation (6) near the equilibrium position can be expressed as(27)Fxe=kxx+kiix,where kx and ki are the displacement stiffness and current stiffness coefficients of the AMB, and they are, respectively, given as(28)kx=μ0ApN2Ib2cos2αg03,ki=−μ0ApN2Ibcosαg02.Under the PD control,ix is given by(29)ix=−Kpx−Kdx˙,where Kp and Kd are the proportional coefficient and derivative coefficient, respectively. Referring to reference [29], the equivalent stiffness ke and equivalent damping ce of the AMB can be expressed as(30)ke=kiAaAsKp−kx,ce=kiAaAsKd,where Aa and As are the gains of the power amplifier and the displacement sensor, respectively.With the equivalent stiffness and equivalent damping of AMBs, respectively, equal to the stiffness and damping of the mechanical bearings, the support characteristics of the system remain unchanged after the substitution of AMBs for mechanical bearings. On this basis, the dynamic responses of the dual-rotor system supported by mechanical bearings are calculated by the dynamic modeling method mentioned in Section2, and the experimental results in [15] are referred to in order to verify the effectiveness of dynamic model. The dual-rotor system and its main parameters are depicted in Figure 7 and Table 1, respectively. In the following calculation, the AMB system parameters are set as α = π/8, Ap = 2.47 × 10−4 m2, N = 160, Ib = 3A, g0 = 5 × 10–4 m, Aa = 0.6 A/V, As = 8 × 103 V/m, Kp = 8.25, and Kd = 0.05.Figure 7
Finite element model of the dual-rotor system [15].Table 1
The main parameters of the dual-rotor system [15].
Physical parameterValuePhysical parameterValueInner rotor lengthli (m)0.706Outer radius of disk 4Rd4 (m)0.127Outer rotor lengthlo (m)0.5011Inner radius of disk 4rd4 (m)0.02Outer radius of inner rotorRi (m)0.0125Thickness of diskshd (m)0.0273Inner radius of inner rotorri (m)0.0075Density of disksρd (kg/m3)7928.56Outer radius of outer rotorRo (m)0.02Elastic modulus of disksEd (Pa)2.1 × 1011Inner radius of outer rotorro (m)0.015Elastic modulus of limiterEp (Pa)2.1 × 1011Density of rotorsρ (kg/m3)7850Poisson ratio of disksμd0.3Elastic modulus of rotorsE (Pa)2.1 × 1011Poisson ratio of limiterμp0.3Poisson ratio of rotorsν0.3Radius of limiter endRP (m)1.25 × 10–4Outer radius of disk 1Rd1 (m)0.125Eccentricities of disk 2 and 4e (m)2 × 10–5Inner radius of disk 1rd1 (m)0.0125Friction coefficientfr0.1Outer radius of disk 2Rd2 (m)0.127Restitution coefficientce0.9Inner radius of disk 2rd2 (m)0.0125Initial rubbing clearanceδ0 (m)1 × 10–5Outer radius of disk 3Rd3 (m)0.125Stiffness of bearingsks (N/m)3 × 106Inner radius of disk 3rd3 (m)0.02Damping of bearingscs (N·s/m)2100In the absence of rubbing, the vertical displacement responses of disk 1 are shown in Figure8 with ωi = 252.6 rad/s and ωo = 301.2 rad/s. It can be seen that the time domain response is performed as beat vibration. Obviously, the simulation results are in consistent with the experimental results. In the case of ωi = 184.1 rad/s and ωo = 276.1 rad/s, the response spectra diagrams of disk 1 are depicted in Figure 9. It is apparent that the response spectrums are characterized by only rotational frequencies ωi and ωo of inner and outer rotors without rubbing, and featured by combined frequency components like ωo + ωi in addition to ωi and ωo when the rubbing occurs. The distinction between the simulation and experimental results is the amplitude of each frequency component, which is mainly caused by the difference in the imbalance position and the initial conditions of the system.Figure 8
Dynamic responses of dual-rotor system without rubbing. (a) Simulation result. (b) Experimental result [15].
(a)(b)Figure 9
Response spectra diagrams of the dual-rotor system. (a) Simulation result without rubbing. (b) Simulation result with rubbing. (c) Experimental result without rubbing [15]. (d) Experimental result with rubbing [15].
(a)(b)(c)(d)From the above analysis, it can be seen that the dynamic modeling method is effective, and the dynamic model is carried out to investigate the rubbing characteristics of the MSDS in the next section.


## 4. Numerical Results and Discussion

The dimension parameters of MSDS structure are shown in Figure1, and the other key parameters are listed in Table 2. For the convenience of analysis, the ratio of the outer rotor speed ωo to the inner rotor speed ωi is defined as the speed ratio rs, namely, rs = ωo/ωi. Numerical simulation is performed to explore the vibration characteristics of the MSDS with fixed-point rubbing. In the simulation, it is assumed that the MSDS is accelerated from zero to the specified speed with a constant acceleration.Table 2
The main parameters of the MSDS.
Physical parameterValuePhysical parameterValueVacuum permeabilityμ0 (H/m)4π × 107Radius of disk 3Rd (m)0.046Half angle between two polesα (°)22.5Radius of limiter endRP (m)1.25 × 10–4Nominal air gap lengthg0 (mm)0.5Poisson ratio of disk 3μd0.3Sectional area of stator poleAp (mm2)247Poisson ratio of limiterμp0.3Coil turnN160Elastic modulus of disk 3Ed (Pa)2.09 × 1011Bias currentIb (A)3Elastic modulus of limiterEp (Pa)2.09 × 1011Maximum control currenticmax (A)3Restitution coefficientce0.9Initial rubbing clearanceδ0 (m)1 × 10–5Friction coefficientfr0.3Interbearing stiffnesskc (N/m)1 × 107Imbalance of disk 1ud1 (kg·m)5 × 10–5Interbearing dampingcc (N·s/m)0Imbalance of disk 3ud3 (kg·m)5 × 10–5


### 4.1. Fixed-Point Rubbing Response

In the case ofKp = 10, Kd = 2, ωi = 700 rad/s and rs = 1.2, time-domain waveforms, dynamic trajectories, and spectrum graphs of disk 1 on the inner rotor and disk 3 on the outer rotor without and with rub-impact are, respectively, depicted in Figure 10. The position of the dotted line in Figure 10(d) represents the initial rubbing clearance. Since the rotational speeds of the inner and outer rotors are relatively close, the time-domain waveforms are exhibited as beat vibration phenomena. Such vibration characteristic is consistent with the research results in many literature [15, 30, 31]. In Figures 10(a) and 10(b) without rub-impact, the rotational frequencies fi and fo concurrently exist in the frequency spectrums of the inner and outer rotor responses, which, respectively, correspond to the unbalance excitation frequencies of the inner and outer rotors and occupy an absolute advantage in the system responses. This reveals the cross-excitation phenomenon in a dual-rotor system due to the coupling effect of interbearings. When the rubbing occurs in Figures 10(c) and 10(d), in addition to the rotational frequency components fi and fo, their combined frequency components also appear in the frequency spectrum, such as fo − fi and 4fi − 3fo, which mainly stems from the combined resonance [32] of the dual-rotor system under the action of the nonlinear rub-impact force. The comparison of the dynamic trajectories shows that the entire steady-state trajectories of the rotors significantly translate downward along the vertical direction for a certain distance and slightly translate to the right along the horizontal direction. It indicates that the position of the entire steady-state trajectory of the system can be changed by fixed-point rubbing.Figure 10
Dynamic responses of the MSDS withKp = 10, Kd = 2 and rs = 1.2. (a) Inner rotor response without rub-impact. (b) Outer rotor response without rub-impact. (c) Inner rotor response with rub-impact. (d) Outer rotor response with rub-impact.
(a)(b)(c)(d)In the Lankarani-Nikravesh model, a complete rubbing process contains a compression stage and a separation stage [21]. The trajectory in Figure 10(d) shows the compression and separation stage during a rubbing process in one motion cycle, which are represented by bold solid and bold dotted lines, respectively. Figure 11 illustrates the change process of the normal impact force, invading length, normal impact angle (the angle is set to zero when the rubbing disappears), and the normal impact velocity of the disk during a motion cycle. In the compression stage, the impact force gradually increases from zero with the increase of the invading length and reaches the maximum at the same time as the invading length, while the impact velocity decreases to almost zero. Then, the rubbing begins to enter the separation stage. In this stage, the impact force and invading length both decrease to zero, while the impact velocity increases. The entire rubbing process is extremely short, lasting only 0.8 millisecond. It can be observed from Figure 11(c) that disk 3 begins to collide the limiter from the bottom right, and the duration of disk 3 located at the bottom right of the limiter is 0.52 milliseconds, which is more than half of the entire rubbing time. This explains why the entire dynamic trajectory of the rotor translates to the lower right after rub-impact. Generally, since part of the energy is lost during collision due to the damping effect of the material, the value of the impact velocity at the beginning of the collision is greater than the value at the end when two objects collide directly along the normal direction of their outer surface. Nevertheless, the phenomenon shown in Figure 11(d) is inconsistent with the conventional case. The main reason is that the impact velocity not only is related to the impact angle, but also depends on the magnitude and direction of whirl speed of disk 3, which can be seen from equation (15). In the rubbing process, the direction of the whirl speed is constantly changing. Therefore, the initial impact velocity may be slightly greater than the final impact velocity. Furthermore, the maximum normal impact force is only 1.1 N, which indicates that the severity of rubbing is slight.Figure 11
Rubbing process of the MSDS withKp = 10, Kd = 2 and rs = 1.2. (a) Normal impact force. (b) Invading length. (c) Normal impact angle. (d) Normal impact velocity
(a)(b)(c)(d)Comparing the dynamic responses of the inner and outer rotors in Figure10, it is apparent that their response laws are similar. Therefore, only the dynamic responses of the outer rotor are analyzed as follows. Figure 12 demonstrates the waterfall plots of the outer rotor response in the case of no rubbing and rubbing with Kp = 10, Kd = 2 and rs = 1.2. What is striking in Figure 12(a) is the appearance of the constant frequency component fc2/2 = 35 Hz in addition to the rotational frequencies of inner and outer rotors, where fc2 equals the second-order critical speed of the linear MSDS (where the AMBs are modeled as equivalent stiffness and equivalent damping model) corresponding to the nonlinear MSDS. The possible reason for fc2 is that the unbalance force increases with the rotational speed, which results in the aggravated rotor vibration. In order to pull the rotor back to the equilibrium position, the AMBs must provide larger electromagnetic forces. At this time, the nonlinearity degree of the electromagnetic force acting on the rotor is intensified, which excites the response frequency component fc2/2. The frequency component fc2/2 appears when the rotational speed of inner rotor exceeds 1280 rad/s, and its amplitude increases sharply. At this moment, the rotor starts to travel irregularly, and the trajectory becomes complicated. In contrast to the waterfall plot without rubbing, the combined frequency component fo − fi is observed in Figure 12(b) when rubbing occurs, and it is submerged in the large amplitude of fc2/2 when the inner rotor speed exceeds 1210 rad/s, while it still exists rather than disappearing. Moreover, the inner rotor speed corresponding to the initial appearance of the constant frequency component is reduced from 1280 rad/s without rubbing to 1170 rad/s, which implies that the existence of rub-impact will cause the frequency component fc2/2 to be stimulated at a lower operating speed.Figure 12
Waterfall plots of the MSDS withKp = 10, Kd = 2 and rs = 1.2. (a) No rubbing. (b) Rubbing.
(a)(b)Speed ratio is a basic and important parameter during the operation process of MSDS. Withrs = 1.5 and other parameters being the same as those of the above case, the dynamic responses and rubbing process of the MSDS are, respectively, presented in Figures 13 and 14. Compared with the rubbing response with rs = 1.2, it can be found in Figure 13 that the number of rub-impacts significantly increases in the same time interval, the dynamic trajectory is relatively simple, and only one combined frequency component fo − fi is observed other than the rotational frequency components. In the rubbing process, the maximum normal impact force is 0.52 N, the rubbing duration is 0.56 milliseconds, the impact angle spans about 1.45 × 10–4 rad, and the maximum normal impact velocity is about 3×10–4 m/s. It is obvious that, under the same conditions, the severity of the rubbing response with rs = 1.5 is lighter than that with rs = 1.2.Figure 13
Dynamic responses of the MSDS withKp = 10, Kd = 2 and rs = 1.5. (a) Outer rotor response without rub-impact. (b) Outer rotor response with rub-impact.
(a)(b)Figure 14
Rubbing process of the MSDS withKp = 10, Kd = 2 and rs = 1.5. (a) Normal impact force. (b) Invading length. (c) Normal impact angle. (d) Normal impact velocity.
(a)(b)(c)(d)Figure15 provides the waterfall plots of the outer rotor response for both no rubbing and rubbing with Kp = 10, Kd = 2 and rs = 1.5. Under the condition of no rubbing, in contrast to the dynamic responses with rs = 1.2, the distinct difference is that the first appearance of frequency component fc2/2 is at a lower rotational speed when rs = 1.5. This may be mainly caused by the sharp increase in the unbalanced force provided by the outer rotor with higher speed ratio. Because the outer rotor speed increases with the speed ratio when the inner rotor speed keeps constant, and the unbalanced force is proportional to the square of the speed, when the rubbing occurs, due to the existence of impact force and the friction force, the germination of frequency component fc2/2 appears at a lower speed, and the fractional combined frequency component (fo − fi)/2 is excited. Since the amplitudes of the combined frequency component fo − fi are two orders of magnitude smaller than those of the rotational frequency components, they are hardly observed in the waterfall plot represented by Cartesian coordinate system. When the inner rotor speed exceeds 960 rad/s, locally continuous spectrum appears at low frequency band, the amplitude of frequency component fc2/2 increases significantly, the rotor moves irregularly, and instability occurs.Figure 15
Waterfall plots of the MSDS withKp = 10, Kd = 2 and rs = 1.5. (a) No rubbing. (b) Rubbing.
(a)(b)In order to further investigate the severity of rubbing, the change curves of maximum normal impact force in the rubbing process with inner rotor speed are illustrated in Figure16 when rs = 1.2 and rs = 1.5. Combining Figure 14 and the waterfall plots in the case of rubbing, it is apparent that the maximum normal impact forces rarely change with the inner rotor speed before the frequency component fc2/2 appears and increase sharply with the speed after fc2/2 appears. More importantly, in the absence of fc2/2, the maximum normal impact forces are only about 1.1 N and 0.5 N for rs = 1.2 and rs = 1.5, respectively. However, the maximum impact force of the dual-rotor system during stable operation can be as high as 100–700 N in reference [14], which is much greater than that of the MSDS. Obviously, MSDS can effectively reduce the severity of fixed-point rubbing faults and significantly attenuate the wear of rotating and stationary components caused by rubbing, thus improving the reliability of components and extending the service life of the aero-engine. Such a huge advantage is attributed to the unique nonlinear support characteristics of the AMB. In addition, the severity of rubbing responses with rs = 1.5 is slightly lighter than that with rs = 1.2 from the perspective of the maximum normal impact force.Figure 16
Maximum normal impact forces of the MSDS withKp = 10 and Kp = 2. (a) rs = 1.2. (b) rs = 1.5.
(a)(b)


### 4.2. Influence of AMB Control Parameters

Generally, support characteristics of bearings usually have significant influence on the vibration characteristics of rotor systems. For the AMB with structural parameters determined, its nonlinear electromagnetic force mainly depends on the control strategy and the rotor motion state. Under the classic PD control, the effects of proportional coefficient and differential coefficient on the MSDS dynamic response are conducted as follows.In the case ofKp = 42, Kd = 2, ωi = 700 rad/s and rs = 1.2, the dynamic responses and the rubbing process of the MSDS are shown in Figures 17 and 18, respectively. Under the condition of no rubbing, the dynamic responses in Figure 17(a) are almost the same as those of Kp = 10. When rubbing occurs, the trajectory of the rotor becomes complicated, which is manifested by the aggregation and interweaving of many similar curves. As a result, the frequency spectrum contains more frequency components. In addition to the rotational frequency components and the combined frequency components found earlier, new combined frequency components 3fi − 2fo and 2fi − fo are seen in Figure 17(b). Besides, one interesting finding is that there is a sideband family whose frequency interval is the 1/7 of rotational frequency difference of the inner and outer rotors (3.2 Hz = 1/7 × (133.7–111.4) Hz) on both sides of these combined frequency components. Compared with the rubbing response with Kp = 10, the maximum impact force is significantly increased, and the rubbing duration increases to 1.3 milliseconds. There is no obvious periodicity in the changes of the rub-impact parameters, and the degree of rub-impact is aggravated.Figure 17
Dynamic responses of the outer rotor withKp = 42, Kd = 2 and rs = 1.2. (a) Outer rotor response without rub-impact. (b) Outer rotor response with rub-impact.
(a)(b)Figure 18
Rubbing process of the MSDS withKp = 42, Kd = 2 and rs = 1.2. (a) Normal impact force. (b) Invading length. (c) Normal impact angle.
(a)(b)(c)(d)WhenKp = 42, Kd = 2 and rs = 1.2, the waterfall plots of the MSDS responses are shown in Figure 19. Compared with the dynamic responses of the outer rotor with Kp = 10 in Figure 12, the difference is that, for the combined frequency component, the amplitudes of fo − fi increase, and 4fi − 3fo appears in a wider speed range, while rubbing occurs. Moreover, the initial instability speed of the inner rotor is reduced from 1220 rad/s to 1130 rad/s. It can be seen that the proportional coefficient has a significant influence on the MSDS rubbing responses.Figure 19
Waterfall plots of the MSDS withKp = 42, Kd = 2 and rs = 1.2. (a) No rubbing. (b) Rubbing.
(a)(b)In order to have a better knowledge of the influence of the proportional coefficient on the MSDS dynamic responses, Figure20 presents the waterfall plots of the outer rotor with Kd = 2, ωi = 700 rad/s and rs = 1.2. In the case of no rubbing, there are only rotational frequency components. When the rubbing occurs, the spectrum contains considerable frequency components other than the rotational frequency components. The combined frequency components fo − fi, 4fi − 3fo and 3fi − 2fo are observed, around which many frequency components are equally spaced with different frequency intervals when Kp ≥ 40. These frequency intervals are the fractions of the rotational frequency difference of the inner and outer rotors. Specifically, these frequency intervals are 3.2 Hz, 2.8 Hz, and 2.2 Hz for 40 ≤ Kp ≤ 54, 55 ≤ Kp ≤ 57, and 58 ≤ Kp ≤ 67, which are about 1/7, 1/8, and 1/10 of the rotational frequency difference 22.3 Hz, respectively. The reason for this phenomena may be that, according to equations (6) and (7), when the rotor has a unit displacement, the AMB will generate greater electromagnetic force with increased Kp; namely, the support stiffness of the MSDS is enhanced, which will lead to the alteration of the system critical speeds. As a result, the critical speed of a certain order of the system increases and gradually approaches the operating speed of the inner rotor, resulting in a more severe rubbing response and causing fractional frequency components of the rotational frequency difference. When Kp ≥ 68, continuous spectrum appears in the low frequency range, the amplitude of the frequency component 2.2 Hz increases rapidly with Kp to even exceed that of the inner rotor rotational frequency, and the rotor begins to move randomly.Figure 20
Waterfall plots of the MSDS withKd = 2, rs = 1.2 and ωi = 700 rad/s. (a) No rubbing. (b) Rubbing.
(a)(b)The effects of proportional coefficient on the maximum normal impact force are demonstrated in Figure21 for the case with Kp = 42, Kd = 2 and rs = 1.2 and the case with Kd = 2, rs = 1.2 and ωi = 700 rad/s. Comparing Figures 16(a) and 21(a), it can be seen that the maximum normal impact force at each speed with Kp = 42 is greater than that with Kp = 10 during stable operation of the MSDS. In Figure 21(b), the maximum normal impact force increases with proportional coefficient.Figure 21
Maximum normal impact forces of the MSDS withKd = 2 and rs = 1.2. (a) Kp = 42. (b) ωi = 700 rad/s.
(a)(b)The influence of differential coefficient on the MSDS dynamic responses is less complicated than that of proportional coefficient. In the case ofKp = 10, ωi = 700 rad/s and rs = 1.2, the steady-state responses of the MSDS with different Kd are almost the same as those with Kd = 2 in Figure 10 when the MSDS is capable of maintaining stability, so the responses are not repeated here. The waterfall plots of the MSDS are demonstrated in Figure 22 when Kp = 10, Kd = 1.2 and rs = 1.2, which are similar with the dynamic responses in Figure 12. The difference is that the frequency component fc2/2 has a very large amplitude as soon as it appears when Kd = 1.2. Besides, the inner rotor speed corresponding to the first appearance of frequency component fc2/2 decreases with the decrease of differential coefficient. As a result, the combined frequency component 4fi − 3fo is also excited. What is interesting about the influence of differential coefficient in Figure 23 is that the duration required for the system to accelerate from zero to the steady-state response decreases with the increase of the differential coefficient. The most likely cause of this phenomenon is that the AMB can produce higher electromagnetic force with increased Kd according to equations (6) and (7); namely, the system will possess better damping characteristics with greater differential coefficient. Consequently, transient large-scale vibration can be attenuated faster, and the system reaches a stable state rapidly. Nevertheless, the instability occurs due to poor damping capacity when Kd < 1.1.Figure 22
Waterfall plots of the MSDS withKp = 10, Kd = 1.2 and rs = 1.2. (a) No rubbing. (b) Rubbing.
(a)(b)Figure 23
Waterfall plots of the MSDS withKp = 10, rs = 1.2 and ωi = 700 rad/s. (a) No rubbing. (b) Rubbing.
(a)(b)The effects of differential coefficient on the maximum normal impact force are demonstrated in Figure24 for the case with Kp = 10, Kd = 1.2 and rs = 1.2 and the case with Kp = 10, rs = 1.2 and ωi = 700 rad/s. Comparing Figures 16(a) and 24(a), it is obvious that the maximum normal impact force at each speed with Kd = 2 is smaller than that with Kd = 1.2. In Figure 24(b), the maximum normal impact force decreases with the increase of differential coefficient.Figure 24
Maximum normal impact forces of the MSDS withKp = 10 and rs = 1.2. (a) Kd = 1.2. (b) ωi = 700 rad/s.
(a)(b)From the above analysis, it can be seen that the AMB control parameters have substantial influence on the dynamic responses of the MSDS. Under the circumstance of grasping the laws of the influence of AMB control parameters on the fixed-point rubbing characteristics of the MSDS, the severity of rubbing can be further attenuated by reasonably selecting the control parameters, which is just the unparalleled advantage of the active controllability of AMB support characteristics over mechanical bearings.


## 4.1. Fixed-Point Rubbing Response

In the case ofKp = 10, Kd = 2, ωi = 700 rad/s and rs = 1.2, time-domain waveforms, dynamic trajectories, and spectrum graphs of disk 1 on the inner rotor and disk 3 on the outer rotor without and with rub-impact are, respectively, depicted in Figure 10. The position of the dotted line in Figure 10(d) represents the initial rubbing clearance. Since the rotational speeds of the inner and outer rotors are relatively close, the time-domain waveforms are exhibited as beat vibration phenomena. Such vibration characteristic is consistent with the research results in many literature [15, 30, 31]. In Figures 10(a) and 10(b) without rub-impact, the rotational frequencies fi and fo concurrently exist in the frequency spectrums of the inner and outer rotor responses, which, respectively, correspond to the unbalance excitation frequencies of the inner and outer rotors and occupy an absolute advantage in the system responses. This reveals the cross-excitation phenomenon in a dual-rotor system due to the coupling effect of interbearings. When the rubbing occurs in Figures 10(c) and 10(d), in addition to the rotational frequency components fi and fo, their combined frequency components also appear in the frequency spectrum, such as fo − fi and 4fi − 3fo, which mainly stems from the combined resonance [32] of the dual-rotor system under the action of the nonlinear rub-impact force. The comparison of the dynamic trajectories shows that the entire steady-state trajectories of the rotors significantly translate downward along the vertical direction for a certain distance and slightly translate to the right along the horizontal direction. It indicates that the position of the entire steady-state trajectory of the system can be changed by fixed-point rubbing.Figure 10
Dynamic responses of the MSDS withKp = 10, Kd = 2 and rs = 1.2. (a) Inner rotor response without rub-impact. (b) Outer rotor response without rub-impact. (c) Inner rotor response with rub-impact. (d) Outer rotor response with rub-impact.
(a)(b)(c)(d)In the Lankarani-Nikravesh model, a complete rubbing process contains a compression stage and a separation stage [21]. The trajectory in Figure 10(d) shows the compression and separation stage during a rubbing process in one motion cycle, which are represented by bold solid and bold dotted lines, respectively. Figure 11 illustrates the change process of the normal impact force, invading length, normal impact angle (the angle is set to zero when the rubbing disappears), and the normal impact velocity of the disk during a motion cycle. In the compression stage, the impact force gradually increases from zero with the increase of the invading length and reaches the maximum at the same time as the invading length, while the impact velocity decreases to almost zero. Then, the rubbing begins to enter the separation stage. In this stage, the impact force and invading length both decrease to zero, while the impact velocity increases. The entire rubbing process is extremely short, lasting only 0.8 millisecond. It can be observed from Figure 11(c) that disk 3 begins to collide the limiter from the bottom right, and the duration of disk 3 located at the bottom right of the limiter is 0.52 milliseconds, which is more than half of the entire rubbing time. This explains why the entire dynamic trajectory of the rotor translates to the lower right after rub-impact. Generally, since part of the energy is lost during collision due to the damping effect of the material, the value of the impact velocity at the beginning of the collision is greater than the value at the end when two objects collide directly along the normal direction of their outer surface. Nevertheless, the phenomenon shown in Figure 11(d) is inconsistent with the conventional case. The main reason is that the impact velocity not only is related to the impact angle, but also depends on the magnitude and direction of whirl speed of disk 3, which can be seen from equation (15). In the rubbing process, the direction of the whirl speed is constantly changing. Therefore, the initial impact velocity may be slightly greater than the final impact velocity. Furthermore, the maximum normal impact force is only 1.1 N, which indicates that the severity of rubbing is slight.Figure 11
Rubbing process of the MSDS withKp = 10, Kd = 2 and rs = 1.2. (a) Normal impact force. (b) Invading length. (c) Normal impact angle. (d) Normal impact velocity
(a)(b)(c)(d)Comparing the dynamic responses of the inner and outer rotors in Figure10, it is apparent that their response laws are similar. Therefore, only the dynamic responses of the outer rotor are analyzed as follows. Figure 12 demonstrates the waterfall plots of the outer rotor response in the case of no rubbing and rubbing with Kp = 10, Kd = 2 and rs = 1.2. What is striking in Figure 12(a) is the appearance of the constant frequency component fc2/2 = 35 Hz in addition to the rotational frequencies of inner and outer rotors, where fc2 equals the second-order critical speed of the linear MSDS (where the AMBs are modeled as equivalent stiffness and equivalent damping model) corresponding to the nonlinear MSDS. The possible reason for fc2 is that the unbalance force increases with the rotational speed, which results in the aggravated rotor vibration. In order to pull the rotor back to the equilibrium position, the AMBs must provide larger electromagnetic forces. At this time, the nonlinearity degree of the electromagnetic force acting on the rotor is intensified, which excites the response frequency component fc2/2. The frequency component fc2/2 appears when the rotational speed of inner rotor exceeds 1280 rad/s, and its amplitude increases sharply. At this moment, the rotor starts to travel irregularly, and the trajectory becomes complicated. In contrast to the waterfall plot without rubbing, the combined frequency component fo − fi is observed in Figure 12(b) when rubbing occurs, and it is submerged in the large amplitude of fc2/2 when the inner rotor speed exceeds 1210 rad/s, while it still exists rather than disappearing. Moreover, the inner rotor speed corresponding to the initial appearance of the constant frequency component is reduced from 1280 rad/s without rubbing to 1170 rad/s, which implies that the existence of rub-impact will cause the frequency component fc2/2 to be stimulated at a lower operating speed.Figure 12
Waterfall plots of the MSDS withKp = 10, Kd = 2 and rs = 1.2. (a) No rubbing. (b) Rubbing.
(a)(b)Speed ratio is a basic and important parameter during the operation process of MSDS. Withrs = 1.5 and other parameters being the same as those of the above case, the dynamic responses and rubbing process of the MSDS are, respectively, presented in Figures 13 and 14. Compared with the rubbing response with rs = 1.2, it can be found in Figure 13 that the number of rub-impacts significantly increases in the same time interval, the dynamic trajectory is relatively simple, and only one combined frequency component fo − fi is observed other than the rotational frequency components. In the rubbing process, the maximum normal impact force is 0.52 N, the rubbing duration is 0.56 milliseconds, the impact angle spans about 1.45 × 10–4 rad, and the maximum normal impact velocity is about 3×10–4 m/s. It is obvious that, under the same conditions, the severity of the rubbing response with rs = 1.5 is lighter than that with rs = 1.2.Figure 13
Dynamic responses of the MSDS withKp = 10, Kd = 2 and rs = 1.5. (a) Outer rotor response without rub-impact. (b) Outer rotor response with rub-impact.
(a)(b)Figure 14
Rubbing process of the MSDS withKp = 10, Kd = 2 and rs = 1.5. (a) Normal impact force. (b) Invading length. (c) Normal impact angle. (d) Normal impact velocity.
(a)(b)(c)(d)Figure15 provides the waterfall plots of the outer rotor response for both no rubbing and rubbing with Kp = 10, Kd = 2 and rs = 1.5. Under the condition of no rubbing, in contrast to the dynamic responses with rs = 1.2, the distinct difference is that the first appearance of frequency component fc2/2 is at a lower rotational speed when rs = 1.5. This may be mainly caused by the sharp increase in the unbalanced force provided by the outer rotor with higher speed ratio. Because the outer rotor speed increases with the speed ratio when the inner rotor speed keeps constant, and the unbalanced force is proportional to the square of the speed, when the rubbing occurs, due to the existence of impact force and the friction force, the germination of frequency component fc2/2 appears at a lower speed, and the fractional combined frequency component (fo − fi)/2 is excited. Since the amplitudes of the combined frequency component fo − fi are two orders of magnitude smaller than those of the rotational frequency components, they are hardly observed in the waterfall plot represented by Cartesian coordinate system. When the inner rotor speed exceeds 960 rad/s, locally continuous spectrum appears at low frequency band, the amplitude of frequency component fc2/2 increases significantly, the rotor moves irregularly, and instability occurs.Figure 15
Waterfall plots of the MSDS withKp = 10, Kd = 2 and rs = 1.5. (a) No rubbing. (b) Rubbing.
(a)(b)In order to further investigate the severity of rubbing, the change curves of maximum normal impact force in the rubbing process with inner rotor speed are illustrated in Figure16 when rs = 1.2 and rs = 1.5. Combining Figure 14 and the waterfall plots in the case of rubbing, it is apparent that the maximum normal impact forces rarely change with the inner rotor speed before the frequency component fc2/2 appears and increase sharply with the speed after fc2/2 appears. More importantly, in the absence of fc2/2, the maximum normal impact forces are only about 1.1 N and 0.5 N for rs = 1.2 and rs = 1.5, respectively. However, the maximum impact force of the dual-rotor system during stable operation can be as high as 100–700 N in reference [14], which is much greater than that of the MSDS. Obviously, MSDS can effectively reduce the severity of fixed-point rubbing faults and significantly attenuate the wear of rotating and stationary components caused by rubbing, thus improving the reliability of components and extending the service life of the aero-engine. Such a huge advantage is attributed to the unique nonlinear support characteristics of the AMB. In addition, the severity of rubbing responses with rs = 1.5 is slightly lighter than that with rs = 1.2 from the perspective of the maximum normal impact force.Figure 16
Maximum normal impact forces of the MSDS withKp = 10 and Kp = 2. (a) rs = 1.2. (b) rs = 1.5.
(a)(b)


## 4.2. Influence of AMB Control Parameters

Generally, support characteristics of bearings usually have significant influence on the vibration characteristics of rotor systems. For the AMB with structural parameters determined, its nonlinear electromagnetic force mainly depends on the control strategy and the rotor motion state. Under the classic PD control, the effects of proportional coefficient and differential coefficient on the MSDS dynamic response are conducted as follows.In the case ofKp = 42, Kd = 2, ωi = 700 rad/s and rs = 1.2, the dynamic responses and the rubbing process of the MSDS are shown in Figures 17 and 18, respectively. Under the condition of no rubbing, the dynamic responses in Figure 17(a) are almost the same as those of Kp = 10. When rubbing occurs, the trajectory of the rotor becomes complicated, which is manifested by the aggregation and interweaving of many similar curves. As a result, the frequency spectrum contains more frequency components. In addition to the rotational frequency components and the combined frequency components found earlier, new combined frequency components 3fi − 2fo and 2fi − fo are seen in Figure 17(b). Besides, one interesting finding is that there is a sideband family whose frequency interval is the 1/7 of rotational frequency difference of the inner and outer rotors (3.2 Hz = 1/7 × (133.7–111.4) Hz) on both sides of these combined frequency components. Compared with the rubbing response with Kp = 10, the maximum impact force is significantly increased, and the rubbing duration increases to 1.3 milliseconds. There is no obvious periodicity in the changes of the rub-impact parameters, and the degree of rub-impact is aggravated.Figure 17
Dynamic responses of the outer rotor withKp = 42, Kd = 2 and rs = 1.2. (a) Outer rotor response without rub-impact. (b) Outer rotor response with rub-impact.
(a)(b)Figure 18
Rubbing process of the MSDS withKp = 42, Kd = 2 and rs = 1.2. (a) Normal impact force. (b) Invading length. (c) Normal impact angle.
(a)(b)(c)(d)WhenKp = 42, Kd = 2 and rs = 1.2, the waterfall plots of the MSDS responses are shown in Figure 19. Compared with the dynamic responses of the outer rotor with Kp = 10 in Figure 12, the difference is that, for the combined frequency component, the amplitudes of fo − fi increase, and 4fi − 3fo appears in a wider speed range, while rubbing occurs. Moreover, the initial instability speed of the inner rotor is reduced from 1220 rad/s to 1130 rad/s. It can be seen that the proportional coefficient has a significant influence on the MSDS rubbing responses.Figure 19
Waterfall plots of the MSDS withKp = 42, Kd = 2 and rs = 1.2. (a) No rubbing. (b) Rubbing.
(a)(b)In order to have a better knowledge of the influence of the proportional coefficient on the MSDS dynamic responses, Figure20 presents the waterfall plots of the outer rotor with Kd = 2, ωi = 700 rad/s and rs = 1.2. In the case of no rubbing, there are only rotational frequency components. When the rubbing occurs, the spectrum contains considerable frequency components other than the rotational frequency components. The combined frequency components fo − fi, 4fi − 3fo and 3fi − 2fo are observed, around which many frequency components are equally spaced with different frequency intervals when Kp ≥ 40. These frequency intervals are the fractions of the rotational frequency difference of the inner and outer rotors. Specifically, these frequency intervals are 3.2 Hz, 2.8 Hz, and 2.2 Hz for 40 ≤ Kp ≤ 54, 55 ≤ Kp ≤ 57, and 58 ≤ Kp ≤ 67, which are about 1/7, 1/8, and 1/10 of the rotational frequency difference 22.3 Hz, respectively. The reason for this phenomena may be that, according to equations (6) and (7), when the rotor has a unit displacement, the AMB will generate greater electromagnetic force with increased Kp; namely, the support stiffness of the MSDS is enhanced, which will lead to the alteration of the system critical speeds. As a result, the critical speed of a certain order of the system increases and gradually approaches the operating speed of the inner rotor, resulting in a more severe rubbing response and causing fractional frequency components of the rotational frequency difference. When Kp ≥ 68, continuous spectrum appears in the low frequency range, the amplitude of the frequency component 2.2 Hz increases rapidly with Kp to even exceed that of the inner rotor rotational frequency, and the rotor begins to move randomly.Figure 20
Waterfall plots of the MSDS withKd = 2, rs = 1.2 and ωi = 700 rad/s. (a) No rubbing. (b) Rubbing.
(a)(b)The effects of proportional coefficient on the maximum normal impact force are demonstrated in Figure21 for the case with Kp = 42, Kd = 2 and rs = 1.2 and the case with Kd = 2, rs = 1.2 and ωi = 700 rad/s. Comparing Figures 16(a) and 21(a), it can be seen that the maximum normal impact force at each speed with Kp = 42 is greater than that with Kp = 10 during stable operation of the MSDS. In Figure 21(b), the maximum normal impact force increases with proportional coefficient.Figure 21
Maximum normal impact forces of the MSDS withKd = 2 and rs = 1.2. (a) Kp = 42. (b) ωi = 700 rad/s.
(a)(b)The influence of differential coefficient on the MSDS dynamic responses is less complicated than that of proportional coefficient. In the case ofKp = 10, ωi = 700 rad/s and rs = 1.2, the steady-state responses of the MSDS with different Kd are almost the same as those with Kd = 2 in Figure 10 when the MSDS is capable of maintaining stability, so the responses are not repeated here. The waterfall plots of the MSDS are demonstrated in Figure 22 when Kp = 10, Kd = 1.2 and rs = 1.2, which are similar with the dynamic responses in Figure 12. The difference is that the frequency component fc2/2 has a very large amplitude as soon as it appears when Kd = 1.2. Besides, the inner rotor speed corresponding to the first appearance of frequency component fc2/2 decreases with the decrease of differential coefficient. As a result, the combined frequency component 4fi − 3fo is also excited. What is interesting about the influence of differential coefficient in Figure 23 is that the duration required for the system to accelerate from zero to the steady-state response decreases with the increase of the differential coefficient. The most likely cause of this phenomenon is that the AMB can produce higher electromagnetic force with increased Kd according to equations (6) and (7); namely, the system will possess better damping characteristics with greater differential coefficient. Consequently, transient large-scale vibration can be attenuated faster, and the system reaches a stable state rapidly. Nevertheless, the instability occurs due to poor damping capacity when Kd < 1.1.Figure 22
Waterfall plots of the MSDS withKp = 10, Kd = 1.2 and rs = 1.2. (a) No rubbing. (b) Rubbing.
(a)(b)Figure 23
Waterfall plots of the MSDS withKp = 10, rs = 1.2 and ωi = 700 rad/s. (a) No rubbing. (b) Rubbing.
(a)(b)The effects of differential coefficient on the maximum normal impact force are demonstrated in Figure24 for the case with Kp = 10, Kd = 1.2 and rs = 1.2 and the case with Kp = 10, rs = 1.2 and ωi = 700 rad/s. Comparing Figures 16(a) and 24(a), it is obvious that the maximum normal impact force at each speed with Kd = 2 is smaller than that with Kd = 1.2. In Figure 24(b), the maximum normal impact force decreases with the increase of differential coefficient.Figure 24
Maximum normal impact forces of the MSDS withKp = 10 and rs = 1.2. (a) Kd = 1.2. (b) ωi = 700 rad/s.
(a)(b)From the above analysis, it can be seen that the AMB control parameters have substantial influence on the dynamic responses of the MSDS. Under the circumstance of grasping the laws of the influence of AMB control parameters on the fixed-point rubbing characteristics of the MSDS, the severity of rubbing can be further attenuated by reasonably selecting the control parameters, which is just the unparalleled advantage of the active controllability of AMB support characteristics over mechanical bearings.


## 5. Conclusions

In this research, an improved nonlinear EF model of AMB is developed by the equivalent magnetic circuit method. With the combination of AMB model and Lankarani-Nikravesh model describing the impact force, the fixed-point rubbing dynamic model of the MSDS is established using finite element method and verified by experimental results in literature. On this basis, the influence of AMB control parameters on the MSDS vibration characteristics is investigated. The main conclusions can be summarized as follows.Due to the coupling effect of interbearings, strong cross-excitation phenomenon exists in the MSDS, and the rotational frequency components of the inner and outer rotors concurrently appear in the spectrum. Combination resonance is excited by the fixed-point rubbing fault, resulting in the occurrence of combined frequency components, such asfo − fi and 4fi − 3fo, as well as a cluster of frequency components equally spaced at different frequency intervals around the combined frequency components. Under the nonlinear EF of AMB, the fractional frequency component of the second-order critical speed of the MSDS can be stimulated in a certain speed range, which indicates the potential instability of the system. Moreover, fixed-point rubbing will reduce the system rotational speed at which this frequency component appears for the first time, which implies the potential instability of the system at a lower speed. As the proportional coefficient decreases, and the differential coefficient increases, the maximum normal impact force in the rubbing process decreases.



---

*Source: 9254013-2022-03-28.xml*