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

Due to the elimination of wear and lubrication system and eﬀective solution to the vibration control issue, magnetic suspended dual-rotor system (MSDS) oﬀers the possibility to signiﬁcantly enhance the performance of aero-engines. However, research on MSDS vibration characteristics under ﬁxed-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 veriﬁed by ﬁnite 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 ﬁnite element method and veriﬁed 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 diﬀerent 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 ﬁxed-point rubbing occurs. With the decrease of proportional coeﬃcient and the increase of diﬀerential coeﬃcient, the maximum normal impact force in the rubbing process decreases.


Introduction
Active magnetic bearing (AMB) is a new type of highperformance 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. e 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 rubimpact tends to result in tragic accidents [2]. erefore, investigation of the MSDS rub-impact characteristics is of vital importance to ensure the stable operation of the aeroengine.
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]. e 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. e 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. e 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. e 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. e rotor dynamic responses affect the support characteristics of AMB through the mechanical system and the control system, which change the system rubbing responses. is 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][20][21][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. e 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.

Mathematical Formulation
For the MSDS structure shown in Figure 1, the inner rotor and the outer rotor are connected through two interbearings, and the outer rotor is supported by two AMBs. e 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) e interbearing is simplified as a linear elastic support. (3) e 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. e local contact stiffness is much smaller than its respective structural stiffness. (5) e influence of thermal effect and friction torque during the rubbing process is ignored. (6) e imbalances of the system are located on disk 1 and disk 3.

Electromagnetic Force Model of AMB.
e conventional used AMB nonlinear EF F xc in the horizontal direction (x direction) is expressed as [17][18][19][20] F xc � μ 0 N 2 A p cos α 4 where μ 0 is permeability of vacuum, N is the coil turn, A p is the sectional area of the stator pole, α is the angle between the central lines of the pole and pole-pair, g 0 is the normal air gap length, I b is the bias current of coils, i x is the control current, and x is the rotor displacement in the horizontal direction. e EF F y 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 Figure 2(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 F 1 and F 2 are the magnetomotive forces of the coils, and R g1 and R g2 are the air gap reluctances of the pole-pair.
According to the equivalent magnetic circuit, the total reluctance R for one pole-pair is where g 1 and g 2 are the air gap lengths between stator poles and rotor. e magnetomotive forces are erefore, the magnetic flux Φ in the air gap is e EF produced by one pole-pair is expressed as With differential drive mode, the nonlinear EF of the AMB in the horizontal direction is obtained as where g 3 and g 4 are the air gap lengths between the righthand side pole-pair and rotor. When the center of AMB rotor moves to the point (x, y), g i (i � 1, 2, 3, 4) are, respectively, given by 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. us, 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 ratio e 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 g 0 � 0.5 mm, N � 172, A p � 80.48 mm 2 , I b � 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 I x � 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. erefore, the IEMF is adopted to establish the rubbing dynamic model of the MSDS in this research.      Shock and Vibration

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. e schematic diagram of the rub-impact model is shown in Figure 5. e normal impact force F n and tangential friction force F t can be expressed as F n � k r r 1. 5 where c e is the restitution coefficient, f r 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, v 0 is the initial impact velocity during each rubbing process, and k r is the contact stiffness expressed as  where μ d , μ p , E d , E p , R d and R p are Poisson's ratio, elastic modulus, and radii of the disk and limiter, respectively. If the initial clearance between the disk and the limiter is δ 0 , the angle of the vertical direction and the line O c O F between the disk center O c (x, y) and the limiter end center us, the invading length r is given by r > 0 implies that the collision occurs. e normal impact angle θ during rubbing is If the whirl speed vector of the disk center is , and the angle α between O c B and the vertical direction is then it can be seen from Figure 3 that the normal impact velocity of the disk in the collision direction is On the basis of knowing r, v and v 0 , the impact force F n and the friction force F t during the rubbing can be obtained by equation (8), and their projections F rx and F ry in the horizontal and vertical directions are, respectively, expressed as It is obvious from equation (8) that the initial impact velocity v 0 will directly affect the magnitude of F n , 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.

Dynamic Model of MSDS.
According to the MSDS structural characteristics shown in Figure 1, the shafts are described by using the Euler-Bernoulli beam element considering moment of inertia and gyroscopic effect. e 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.
According to the rotor dynamics theory and finite element method, the equations of motion of the inner rotor are where M i , J i and K i are mass matrix, gyroscopic matrix, and stiffness matrix of the inner rotor, respectively. ω i is the inner rotor rotational speed. q 1i and q 2i are the generalized displacement vectors of the inner rotor, and they are expressed as where x k and y k (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. Q 1i and Q 2i are generalized force vectors and are expressed as where u d1 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 Q 1i and Q 2i is 19.
Similarly, the equations of motion of the outer rotor are where M o , J o and K o are mass matrix, gyroscopic matrix, and stiffness matrix of the outer rotor, respectively. ω o is the outer rotor rotational speed. q 1o and q 2o are described as Q 1o and Q 2o are expressed as where u d3 is the imbalance on the disk 3. F L mx , F L my , F R mx and F R my are EFs in the x and y directions of the AMBs at the left and right hand sides, respectively. e position numbers of the three non-zero terms in Q 1o are, respectively, 13, 27, and 35, and the same for Q 2o .
Based on the above analysis, the MSDS dynamic equations can be expressed as 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 6 Shock and Vibration In the above modeling process, the effects of the support stiffness and support damping of interbearings have not been considered yet. erefore, 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 to K(5, On this basis, the Newmark-β method is adopted to obtain the dynamic response of MSDS with fixed-point rubbing fault.

Validation of the Dynamic Model
In this section, the dynamic model of the MSDS with fixedpoint 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 where k x and k i are the displacement stiffness and current stiffness coefficients of the AMB, and they are, respectively, given as Under the PD control, i x is given by where K p and K d are the proportional coefficient and derivative coefficient, respectively. Referring to reference [29], the equivalent stiffness k e and equivalent damping c e of the AMB can be expressed as where A a and A s 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 Section 2, and the experimental results in [15] are referred to in order to verify the effectiveness of dynamic model. e 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, In the absence of rubbing, the vertical displacement responses of disk 1 are shown in Figure 8 with ω i � 252. 6   Shock and Vibration 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. e 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.
From the above analysis, it can be seen that the dynamic modeling method is effective, and the dynamic model is     Shock and Vibration carried out to investigate the rubbing characteristics of the MSDS in the next section.

Numerical Results and Discussion
e dimension parameters of MSDS structure are shown in Figure 1, 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 r s , namely, r s � ω 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.

Fixed-Point Rubbing
Response. In the case of K p � 10, K d � 2, ω i � 700 rad/s and r s � 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. e 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 f i and f o 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. is 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 f i and f o , their combined frequency components also appear in the frequency spectrum, such as f o − f i and 4f i − 3f o , which mainly stems from the combined resonance [32] of the dualrotor system under the action of the nonlinear rub-impact force. e 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 fixedpoint rubbing.
In the Lankarani-Nikravesh model, a complete rubbing process contains a compression stage and a separation stage [21]. e 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. en, 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. e 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. is 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. e 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. erefore, 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.
Comparing the dynamic responses of the inner and outer rotors in Figure 10, it is apparent that their response laws are similar. erefore, 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 K p � 10, K d � 2 and r s � 1.2. What is striking in Figure 12  frequency component f c2 /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 f o − f i is observed in Figure 12(b) when rubbing occurs, and it is submerged in the large amplitude of f c2 /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 f c2 /2 to be stimulated at a lower operating speed.
Speed ratio is a basic and important parameter during the operation process of MSDS. With r s � 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 r s � 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 f o − f i 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 r s � 1.5 is lighter than that with r s � 1.2.  Shock and Vibration Figure 15 provides the waterfall plots of the outer rotor response for both no rubbing and rubbing with K p � 10, K d � 2 and r s � 1.5. Under the condition of no rubbing, in contrast to the dynamic responses with r s � 1.2, the distinct difference is that the first appearance of frequency component f c2 /2 is at a lower rotational speed when r s � 1.5. is 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 f c2 /2 appears at a lower speed, and the fractional combined frequency component (f o − f i )/2 is excited. Since the amplitudes of the combined frequency component f o − f i 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 f c2 /2 increases significantly, the rotor moves irregularly, and instability occurs.   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 Figure 16 when r s � 1.2 and r s � 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 f c2 /2 appears and increase sharply with the speed after f c2 /2 appears. More importantly, in the absence of f c2 /2, the maximum normal impact forces are only about 1.1 N and 0.5 N for r s � 1.2 and r s � 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 r s � 1.5 is slightly lighter than that with r s � 1.2 from the perspective of the maximum normal impact force.

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 of K p � 42, K d � 2, ω i � 700 rad/s and r s � 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 K p � 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 3f i − 2f o and 2f i − f o 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 K p � 10, the maximum impact force is significantly increased, and the rubbing duration increases to 1.3 milliseconds. ere is no obvious periodicity in the changes of the rub-impact parameters, and the degree of rub-impact is aggravated.
When K p � 42, K d � 2 and r s � 1.2, the waterfall plots of the MSDS responses are shown in Figure 19. Compared with the dynamic responses of the outer rotor with K p � 10 in Figure 12, the difference is that, for the combined frequency component, the amplitudes of f o − f i increase, and 4f i − 3f o 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.
In order to have a better knowledge of the influence of the proportional coefficient on the MSDS dynamic responses, Figure 20 presents the waterfall plots of the outer rotor with K d � 2, ω i � 700 rad/s and r s � 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.   different frequency intervals when K p ≥ 40. ese 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. e 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 K p ; 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 K p ≥ 68, continuous spectrum appears in the low frequency range, the amplitude of the frequency component 2.2 Hz increases rapidly with K p to even exceed that of the inner rotor rotational frequency, and the rotor begins to move randomly. e effects of proportional coefficient on the maximum normal impact force are demonstrated in Figure 21 for the case with K p � 42, K d � 2 and r s � 1.2 and the case with K d � 2, r s � 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 K p � 42 is greater than that with K p � 10 during stable operation of the MSDS. In Figure 21(b), the maximum normal impact force increases with proportional coefficient. e influence of differential coefficient on the MSDS dynamic responses is less complicated than that of proportional coefficient. In the case of K p � 10, ω i � 700 rad/s and r s � 1.2, the steady-state responses of the MSDS with different K d are almost the same as those with K d � 2 in Figure 10 when the MSDS is capable of maintaining stability, so the responses are not repeated here. e waterfall plots of the MSDS are demonstrated in Figure 22 when K p � 10, K d � 1.2 and r s � 1.2, which are similar with the dynamic responses in Figure 12. e difference is that the frequency component f c2 /2 has a very large amplitude as soon as it appears when K d � 1.2. Besides, the inner rotor speed corresponding to the first appearance of frequency component f c2 /2 decreases with the decrease of differential coefficient. As a result, the combined frequency component 4f i − 3f o 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. e most likely cause of this phenomenon is that the AMB can produce higher electromagnetic force with increased K d 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 K d < 1.1.  e effects of differential coefficient on the maximum normal impact force are demonstrated in Figure 24 for the case with K p � 10, K d � 1.2 and r s � 1.2 and the case with K p � 10, r s � 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 K d � 2 is smaller than that with K d � 1.2. In Figure 24(b), the maximum normal impact force decreases with the increase of differential coefficient.
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.

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.
e 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 as f o − f i and 4f i − 3f o , 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. Electromagnetic forces in the x and y directions of the left and right side AMBs, respectively F n , F t :

Nomenclature
Normal impact force and tangential friction force in the rub-impact model, respectively F rx , F ry : Rubbing forces in the horizontal and vertical directions, respectively F xc , F xe : e nonlinear electromagnetic forces of the conventional model and the proposed model, respectively g 0 : e normal air gap length, m g 1 , g 2 , g 3 , g 4 : e actual air gap lengths between the poles and rotor surface i x : Control current in the horizontal direction I b : Bias current J i , J o : Gyroscopic matrices of inner and outer rotors, respectively k r : Contact stiffness between the disk and fixed limiter K p , K d : e proportional coefficient and differential coefficient, respectively K: Overall stiffness matrix of the MSDS K i , K o : Stiffness matrices of inner and outer rotors, respectively M: Overall mass matrix of the MSDS M i , M o : Mass matrices of inner and outer rotors, respectively N: Coil turn q: Overall generalized displacement vector of the MSDS q 1i , q 2i , q 1o , q 2o : Generalized displacement vectors of inner and outer rotors, respectively Q: Overall generalized force vector of the MSDS Q 1i , Q 2i , Q 1o , Q 2o : Generalized force vectors of inner and outer rotors, respectively r s : Speed ratio of the dual-rotor system R d , R p : Radii of the disk and end of the fixed limiter, respectively R g1 , R g2 : e air gap reluctances between the polepair and rotor surface x k , y k : Translational displacements of the nodes of the MSDS, k � 1, 2, 3, . . ., 35 α: e angle between the central lines of the pole and pole-pair, α � 22.5°μ 0 : Permeability of vacuum μ d , μ p : Poisson's ratio of the disk and fixed limiter, respectively δ 0 : Initial clearance between the disk and fixed limiter ω i , ω o : Rotational speed of the inner and outer rotors, respectively ϕ yk , ϕ xk : Angular displacement of the nodes of the MSDS, k � 1, 2, 3, . . ., 35 Φ: Magnetic flux in the air gap.
Data Availability e underlying data supporting the research results in the manuscript can be accessed in the website https://pan.baidu. com/s/1SnwmpVXNnv0Cd9ArzXO93Q, and the data extraction code is 4512.

Conflicts of Interest
e authors declare that they have no conflicts of interest.