#### Abstract

The magnetic suspended dual-rotor system (MSDS) can effectively increase the thrust weight ratio of aeroengines. However, the MSDS dynamic characteristics have rarely been investigated. In this research, a MSDS with the outer rotor supported by two active magnetic bearings (AMBs) is designed, and the PID control is employed. The Riccati transfer matrix method using complex variables is adopted to establish the MSDS dynamic model. Subsequently, the influences of AMBs’ control parameters on the MSDS dynamic characteristics are explored. According to the analysis, two rigid mode shapes remain unchanged with the variation of the relationship between their corresponding damped critical speeds (DCSs). Moreover, the rigid DCSs disappear with large derivative coefficient. Eventually, the validity of the dynamic model and the appearance of rigid DCSs are verified.

#### 1. Introduction

Considering many advantages of active magnetic bearings (AMBs) over conventional bearings such as contactless operation, low power consumption, high reliability, adjustable dynamic characteristics, and extended life [1], the substitution of AMBs for mechanical bearings to support the dual-rotor system in the aeroengine, namely, the magnetic suspended dual-rotor system (MSDS), can greatly enhance the system functionality, reduce the system complexity and the engine weight, save the cost, optimize the engine structure, and obtain higher reliability and maintainability [2]. Therefore, MSDS provides the possibility to significantly improve aeroengine performance.

It is already known that the prediction of critical speed of rotating systems is indispensable to the rotor design. With negligible damping in the bearing-rotor system, undamped resonant frequencies are calculated as the conventional critical speeds. Generally, damping is the basic parameter of rotor dynamic analysis, and it possesses a great effect on the dynamic characteristics of various types of bearing-rotor systems. In the general case, with the consideration of damping, it is the damped natural frequencies that are to be determined, namely, the damped critical speeds (DCSs), which establish the actual critical speeds of bearing-rotor systems, including the stiffening effect of the damping in the bearings. These results give a more realistic base than the conventional critical speed calculation for assessing any potential trouble from passing through or operating close to a critical speed [3, 4]. For the MSDS, one of the outstanding features is that the stiffness and damping of the AMBs are controllable and can be adjusted in a wide range with the control algorithm, which requests the exploration of the critical speeds with variable damping. To achieve excellent dynamic characteristics, considerable damping is necessarily introduced by the AMBs to control the vibration of the dual-rotor system, whose ignorance will seriously reduce the accuracy of dynamic analysis results. Therefore, it is of vital importance to investigate the DCSs of the MSDS, which is the main interest in this paper.

As one of the key components of MSDS, in the last few decades, AMBs have been under heavy development and found applications in rotating machinery where conventional bearings might not be feasible. Researches on the bearing characteristics of AMBs have drawn an increasing attention. In addition to the bearing capacity models of AMBs taken flux leakage [5, 6] or magnetic saturation [7, 8] into account, much more concern has been focused on the rotor dynamics analysis of the AMB-rotor system, including the calculation of the critical speeds, which are based on the equivalent stiffness coefficients and equivalent damping coefficients of the AMBs [9, 10]. Guojun et al. [11] adopted AMBs to support the rotor of a 10 MW high-temperature gas-cooled reactor and analyzed the effects of AMB stiffness on the critical speeds of the rotor system. However, they ignored the effects of AMB damping characteristics on the system critical speeds. Liu et al. [12] established a rigid dynamic model of the AMB-rotor system by Lagrange equations, discussed the dynamic characteristics of an energy storage flywheel AMB-rotor system using PID control, and analyzed the effects of control parameters and the position of AMBs on the system critical speeds. Tang et al. [13] investigated the relationship between modal frequency and stiffness and damping in a magnetic suspended control moment gyroscope with radial AMB and axial passive magnetic bearing. Even though these publications considered the AMB equivalent damping, their emphases are not on the DCSs. In addition, the application of AMBs in the above work is limited to the single-rotor systems. Due to the introduction of interbearings into the dual-rotor system, the AMBs in the MSDS will need not only to resist the external load but also to counteract the load caused by the coupling of the internal motion of the dual-rotor system. As a result, the influence of AMB bearing characteristics on the MSDS dynamic characteristics is more complex. For the purpose of better control system design of the MSDS in the future, it is essential to achieve a more thorough understanding of the influence of AMB bearing characteristics, especially for the effects of AMB equivalent damping on the critical speeds, which are the base of further rotor dynamic analysis.

On the other hand, the dual-rotor structure, namely, the core component of the MSDS, is widely used in aeroengines. Researches on the dynamic characteristics of dual-rotor systems are of great significance for lowering vibration and noise and enhancing their overall performance. Considerable theoretical work has been carried out on dual-rotor systems supported by mechanical bearings, covering various aspects like critical speeds, unbalance response [14, 15], stability analysis [16], transient response [17, 18], rubbing characteristic [19, 20], and nonlinear dynamics [21, 22]. On the research of critical speeds, Chiang et al. [23] developed the dynamic model of a dual-rotor turbo-machinery system based on finite element model (FEM) to predict the critical speeds and investigated the effects of the speed ratio and the interbearing stiffness on the critical speeds of the system. Fei et al. [24] built the vibration equation of a dual-rotor system with FEM to explore the effects of interbearing stiffness on the critical speeds of the system when co-rotating and counterrotating. Wang et al. [25] established the dynamic model of a dual-rotor system by FEM to determine the system critical speeds of the first two orders and verified the calculation by experimental results. It should be noted that all of these papers calculated the undamped critical speeds due to the reason that their dual-rotor systems were merely supported by rolling bearings with negligible damping. Although damping was taken into consideration in the dual-rotor system supported by squeeze film dampers by some researchers [26–28], the influence of damping on the dynamic characteristics was not considered owing to the fact that damping in mechanical bearings can hardly be precisely controlled. However, for better vibration abatement and stable operation, an appropriate amount of controllable damping will be inevitably introduced into the AMB-rotor system [29, 30], including the MSDS. Thus, it is the DCSs that need to be determined. Relatively, there are only a small number of publications performed on the dual-rotor systems supported by AMB, i.e., MSDS. Ebrahimi et al. established the nonlinear motion equations of a MSDS by Lagrange’s equations to investigate the chaotic vibration [31] and bifurcation behavior [32], where the MSDS contains two AMBs and an interbearing. The current work is mainly focused on the nonlinear dynamic analysis of MSDS, and it is insufficient to solve the problem of dynamic characteristics description for MSDS during the practical operation. Furthermore, the basic dynamic analysis including the DCSs has not been addressed, which is the foundation of the design of control system of MSDS. Therefore, it is necessary to predict the DCSs for the safety design and stable operation of the MSDS in aeroengines.

In this research, a structure of MSDS is designed, whose inner and outer rotors are interlinked by two rolling bearings as interbearings and the outer rotor is supported by two AMBs with PID control. Then, the Riccati transfer matrix method using complex variables and the fully coupled unit model are employed to develop the dynamic equation of the system, which is solved by Muller’s quadratic interpolation technique. Subsequently, the influences of the AMBs’ control parameters on the DCSs and stability are carried out. Finally, the dynamic model and some of the conclusions are verified.

#### 2. System Modelling

The structure and theoretical model of the MSDS are depicted in Figure 1, where the inner and outer rotors are interlinked through two rolling bearings acting as intershaft bearings, and the outer rotor is supported by two AMBs. In the MSDS, the left end of the outer rotor and the right end of the inner rotor are, respectively, connected to a motor through an elastic diaphragm coupling, which allows a certain axial deviation and deflection angle. Thus the axial degree of freedoms (DOFs) of the inner and outer rotors can be regulated. Therefore, the main focus is on the translational and rotational DOFs of the MSDS in radial direction in this research.

**(a)**

**(b)**

##### 2.1. AMB Modelling

The bearing characteristics are the foundation of the rotor dynamics analysis for either AMBs or traditional mechanical bearings. In mechanical bearings, the stiffness and the damping are the commonly used parametric representations for bearing characteristics, and they are also used for AMBs, namely, the equivalent stiffness and the equivalent damping [13, 33, 34]. Compared with mechanical bearing system, the AMB system is open loop unstable and feedback control is needed for rotor levitation. PID control is adopted to model the AMB bearing characteristics herein [35, 36].

For an eight-pole AMB with differential drive mode, the basic control loop and control block diagram in the *y* axis are shown in Figure 2. In the case of ignoring flux leakage and magnetic saturation, according to [13, 33], the equivalent stiffness *k*_{e} and equivalent damping *d*_{e} of the AMB can be expressed aswhere *K*_{P}, *K*_{I}, and *K*_{D} are the proportional, integral, and derivative coefficients, respectively, *ω* is the whirl frequency of the system, *A*_{a} the gain of the power amplifier, *A*_{s} is the displacement sensor, and *k*_{s} and *k*_{i} are the displacement stiffness coefficient and current stiffness coefficient of the AMB, and they are, respectively, given aswhere *μ*_{0} is the permeability of vacuum, *α* the pole face angle, *A*_{p} the sectional area of one stator pole, *N* the coil turn, *I*_{b} the bias current of coils, and the nominal air-gap length.

**(a)**

**(b)**

Equation (1) indicates that the equivalent stiffness and equivalent damping are determined by not only the AMB’s size and structure but also the control system, which reveals the fact that the equivalent stiffness and equivalent damping of AMB can be adjusted by the control system when the AMB’s size and structure have been determined. It must be noted that the equivalent stiffness and equivalent damping of AMBs are linear, which is based on the small displacement hypothesis, namely, the AMB-rotor operates near the equilibrium position. In fact, the assumption is usually satisfied during the practical operating process of AMBs. Considering the weak coupling effects of the AMB in *x* and *y* axes, the equivalent cross stiffness and equivalent cross damping are neglected. In this way, a rotor is supported by two AMBs. Each AMB support can be equivalent to two mutually perpendicular spring-damping structures. The spring-damping structure of each direction is completely independent. Therefore, the decentralized PID control is adopted in this research.

##### 2.2. Dynamic Model of MSDS

The mathematical model of the MSDS is shown in Figure 1(b). Lumped masses *a* connected by massless beam elements *b* represent the dual-rotor of the MSDS, where bearing element *c* expressed by eight linearized coefficients is modelled as forces acting on the masses at the appropriate axial locations. For simplicity, the following assumptions are made to build the dynamic model:(1)The position where the magnetic force is applied is considered as the place where the displacement is tested by the sensor(2)The dynamic properties of bearings are represented by 8 linearized coefficients, four (*k*_{xx}, *k*_{xy}, *k*_{yx}, *k*_{yy}) of which represent the stiffness properties and the other four (*c*_{xx}, *c*_{xy}, *c*_{yx}, *c*_{yy}), the damping properties, where *k*_{xy} = *k*_{yx} = 0 and *c*_{xy} = *c*_{yx} = 0(3)Flexibilities and cross-coupling force effects of the foundation are ignored

In the MSDS, there are sixteen state variables, two deflections, two slopes, two bending moments, and two shear force components for inner and outer rotors, respectively. The Riccati transfer matrix method divides the sixteen state variables *Z* into two groups, namely, *Z* = {*f e*}, where *f* = {*M*_{yi}*Q*_{xi}*M*_{xi}*Q*_{yi}*M*_{yo}*Q*_{xo}*M*_{xo}*Q*_{yo}}^{T} and *e* = {*x*_{i}*θ*_{yi}*y*_{i}*θ*_{xi}*x*_{o}*θ*_{yo}*y*_{o}*θ*_{xo}}^{T}, and these parameters are described in nomenclature. For the simplification of interbearing modelling and being easy to program, the fully coupled unit model is adopted, and it is shown in Figure 3(a). The sign convention for the constituents of the state vector in the *x*-*z* and *y*-*z* planes is described in Figure 3(b).

**(a)**

**(b)**

In the fully coupled unit model, it is assumed that the inner and outer rotors are divided into many elements *a*_{i}, *b*_{i}, *a*_{o} and *b*_{o}, respectively, and the total element number *N* of the inner rotor should be equal to that of the outer rotor for logical alignment and overall transfer of state variables (for example, in Figure 1(b)). Then, an inner rotor element *a*_{i}, an outer rotor element *a*_{o}, an interbearing element *c*_{c} which interlinks the inner rotor element *a*_{i} and the outer rotor element *a*_{o}, and two elements *c*_{i} and *c*_{o} to only, respectively, support the inner rotor element *a*_{i} and the outer rotor element *a*_{o} are organized to assemble the element *d* shown in Figure 3(a). Finally, an overall transfer element of MSDS, namely, the fully coupled unit model, is constructed by an inner rotor element *b*_{i}, an outer rotor element *b*_{o}, and an element *d*. To achieve this mathematically, premultiply the element matrix *b* by the element matrix *d*. Therefore, the MSDS can be represented by *N* fully coupled units whose constituents are completely identical. When there are no bearings at the *k*th (1 ≤ *k* ≤ *N*) element, their stiffness and damping coefficients are just set to be zero.

Considering the MSDS represented by *N* sectional matrices, the relationship of state vectors for the *k*th element *d* is written aswhere and are the right-hand and left-hand side state vectors of the *k*th element *d*, respectively. The superscript “T” in represents the state vectors transposed. It should be noted that if these constituents of partitioned matrices *D*_{11}, *D*_{12}, *D*_{21}, and *D*_{22} involving the support characteristics of interbearings are all set equal to be zero, equation (3) reduces to that of two uncoupled straight rotors. The state vectors for the *k*th inner and outer rotor elements *b* are given aswhere and are the right-hand and left-hand side state vectors of the *k*th element *b*, respectively. According to the geometric compatibility and equilibrium of internal forces at the junction between the *k*th elements *d* and *b*, it is required that , and the same for the junction between the *k*th element *b* and (*k *+* *1)th element *d* is . Furthermore, omitting the subscripts *b* and *d*, the relationship of state vectors for the *k*th fully coupled unit is expressed aswhere

The partitioned matrices *B*_{11}, *B*_{12}, *B*_{21}, *B*_{22}, *D*_{11}, *D*_{12}, *D*_{21}, and *D*_{22} are given in Appendix. Omitting the superscript *L*, the expansion of equation (5) gives

Introducing a Riccati transformation at station *k*, we getwhere the 8 × 8 matrix *S* is called the Riccati transfer matrix, and it is a function of the complex variable *s*. Equation (7) is further described as

According to equations (8) and (9), the recurrence formula of the Riccati transfer matrix is

Since the left-hand side boundary conditions are *f*_{1} = 0 and *e*_{1} ≠ 0, the initial condition is *S*_{1} = 0. Knowing *u*_{11}, *u*_{12}, *u*_{21}, and *u*_{22} for each fully coupled unit, *S*_{k + 1} can be obtained sequentially by using equation (10) repeatedly. For the right-hand side of station *N*, the relationship of the state vectors *f*_{N+1} and *e*_{N+1} is

With the right-hand side boundary conditions *f*_{N+1} = 0 and *e*_{N+1} ≠ 0, a nontrivial solution of equation (10) should satisfy the condition that the determinant of *S*_{N+1} is zero, namely,

However, there are many opposite-sign-infinite-type singularities in the residual quantity curve of Δ_{1}, and this makes the opposite-sign intervals and root-existence intervals not one-to-one correspondences and results in wrong roots and leaking roots. In order to avoid wrong roots and leaking roots, the technique described in [37] is employed to change opposite-sign-infinite-type singularities to same-sign-infinite-type singularities. In this way, the nonsingularity frequency equation is transformed aswhere the “sign” and “∏” represent the sign function and continued multiplication, respectively. Substituting the equivalent stiffness and equivalent damping of the AMB in equation (1) for the bearing characteristics to support the outer rotor, the dynamic model of MSDS is obtained. Herein, Muller’s quadratic interpolation technique [4] is employed to solve equation (13).

For a given operating condition, the solution to the elliptical motion is assumed aswherewhere *ω* (rad/s) is the damped natural frequency of the system and *λ* (rad/s) is corresponding damping exponent, which represents the decay or growth of the vibration. If *λ* < 0, the vibration dies out exponentially with time and the system is stable. Conversely, *λ* > 0 represents an unstable system. Therefore, the DCSs and stability of the MSDS can be determined simultaneously from the eigenvalue *s*.

##### 2.3. Stability Criterion

As mentioned above, *λ* in equation (15) can be used to evaluate the stability of the system. For further quantify representation, the logarithmic decrement is adopted to conveniently expressed the stability margin of the rotor system, and it is defined as

Obviously, *δ* > 0 and *δ* < 0 correspond to a stable and unstable rotor system, respectively. When the value of *δ* exceeds 1, that particular mode is well damped [3].

#### 3. Dynamic Characteristics Analysis of MSDS

To demonstrate the application of the above analysis, simulation is carried out to investigate the dynamic characteristics of MSDS. The basic dimensions (unit, mm) of the dual-rotor are depicted in Figure 1(a), and the design parameters of AMBs are listed in Table 1.

##### 3.1. Basic Dynamic Characteristics

In a dual-rotor system, co-rotation and counterrotation are two possible operation modes. The former one is with both rotors co-rotating and the latter one is with the rotors counterrotating with respect to each other. The main difference between the two operation modes is caused by gyroscopic moment. Due to the fact that gyroscopic moment has a minor effect on the amplitudes of DCS characteristics and it has little effect on the trends of DCS characteristics with the variations of AMBs’ control parameters, the typical case of co-rotation is considered in this research. The speed ratio *r* of the system is defined as the ratio of outer rotor speed *ω*_{o} to inner rotor speed *ω*_{i}, i.e., *r* = *ω*_{o}/*ω*_{i}. In the case of interbearing stiffness *k*_{m} = 1 × 10^{7} N/m, damping *c*_{m} = 0, and control parameters *K*_{P} = 8, *K*_{I} = 0, *K*_{D} = 0.01, and *r* = 1.5, the Campbell diagram and corresponding logarithmic decrements of the MSDS are plotted in Figure 4. It should be noted that the intersections (marked with “o”) of the hyphen lines and the frequency lines in Figure 4(a) represent the DCSs, which are listed in Table 2.

**(a)**

**(b)**

Figure 4(a) shows the variation of the system whirl frequencies of the first three orders with the rotating speed. Among these, the first order frequency curves of both forward and backward whirl overlap and they are almost horizontal straight lines. In Figure 4(b), the logarithmic decrements of the first two orders are greater than 1, which indicates that the first two order modes are well damped, and the vibrations of the system can be effectively attenuated when passing through or operating close to the first two order DCSs. However, the values of logarithmic decrements for the third order are slightly greater than 0, and the corresponding stability margin is very small.

The first three mode shapes of the MSDS are shown in Figure 5, where the first order mode shapes of both inner and outer rotors are the translational modes, the second order the conical modes, and the third order the first order bending modes. Therefore, combined with Figures 4 and 5, the MSDS can safely pass through the rigid DCSs, while it may encounter instability problems experiencing the first-order bending DCS.

**(a)**

**(b)**

As stated above, the equivalent stiffness and equivalent damping are dependent on the controller with the AMB’s size and structure determined, and they generally have substantial effects on the DCSs, vibration modes, and system stability. Therefore, it is necessary to investigate the influence of PID control parameters on the MSDS dynamic characteristics.

The effect of proportional coefficient on the MSDS dynamic characteristics with *K*_{D} = 0.01 and *K*_{I} = 0 is given in Figure 6, where the subscript number in *K*_{P1}, *K*_{P2}, and *K*_{P3} means the corresponding order of whirl frequency, and the same for that of *K*_{D1}, *K*_{D2}, and *K*_{D3} in Figure 7 and *K*_{I1}, *K*_{I2}, and *K*_{I3} in Figure 8. The ordinate values of the intersections of the lines (*ω*_{i} = *ω*, *ω*_{o} = 1.5*ω*) and forward whirl frequency curves are the DCSs. With the increase of proportional coefficient *K*_{P}, the first two order DCSs increase, while the corresponding logarithmic decrements decrease. The main reason is that increasing the proportional coefficient leads to increased equivalent stiffness, which provides the tendency to oscillate. Therefore, the proportional coefficient plays an important role in the DCSs and stability of the MSDS rigid modes.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

The effect of derivative coefficient with *K*_{P} = 8 and *K*_{I} = 0 is given in Figure 7. Increasing the derivative coefficient *K*_{D}, the first two order DCSs decrease, while the corresponding logarithmic decrements increase. This is due to that the increased derivative coefficient leads to the increased equivalent damping, which reduces the oscillation of the system. Therefore, the derivative coefficient is also a key parameter to determining the DCSs and stability of the rigid modes.

The effect of integral coefficient with *K*_{P} = 8 and *K*_{D} = 0.01 is given in Figure 8. Increasing the integral coefficient *K*_{I}, the first two order DCSs almost do not vary, and the corresponding logarithmic decrements slightly decrease. This may be due to the fact that the integral coefficient has no effect on equivalent stiffness and has some effect on equivalent damping of the system, which can be concluded from equation (1).

From the above analysis, the proportional and derivative coefficients possess significant influence on the dynamic characteristics of the MSDS rigid modes, while the integral coefficient has a minor effect on the rigid modes. However, these PID control parameters have little influence on the dynamic characteristics of the first-order bending mode. In order to further explore the influence law of these parameters on the dynamic characteristics in detail, more cases are performed from the perspective of DCS characteristics.

##### 3.2. Effect of Derivative Coefficient

The effect of derivative coefficient with *K*_{I} = 0, *K*_{P} = 5 and *c*_{m} = 0 is given in Figure 9. In Figure 9(a), with the derivative coefficient increasing, the DCSs of the first two orders decrease to zero and that of the third order slightly decreases to a certain value. In Figure 9(b), the arrow direction represents the direction in which the derivative coefficient increases. It can be observed that the absolute values of damping exponents of the first two orders increase, and along with the process, the corresponding DCSs reduce to zero. For the third order, the absolute value of damping exponent firstly increases and then decreases. Thus, according to equation (16), the logarithmic decrements *δ* of the first two orders increase dramatically, while that of the third order firstly increases and then it decreases, as shown in Figures 9(c) and 9(d). These observations reveal the process how the rigid modes gradually disappear with the derivative coefficient. When the rigid DCSs decrease, the corresponding absolute values of damping exponent increase and their logarithmic decrements approach infinity during the process.

**(a)**

**(b)**

**(c)**

**(d)**

In Figure 9(a), there is an intersection between the first two order DCS curves. In order to illustrate the variation of dynamic characteristics with the derivative coefficient, the mode shapes corresponding to two derivative coefficients before and after the intersection, for example, *K*_{D1} = 0.0047 and *K*_{D2} = 0.0123 shown in Figure 9(a), are investigated. The first three mode shapes of MSDS corresponding to *K*_{D1} = 0.0047 and *K*_{D2} = 0.0123 are depicted in Figures 10 and 11, respectively. In Figure 10, the first-order mode shapes of both inner and outer rotors are translational modes, the second-order conical modes, and the third-order bending modes, which are the same as those shown in Figure 5. In Figure 11, it is obvious that the first three order mode shapes of MSDS are almost the same as those in Figure 10, even though the first-order DCS is greater than that of the second order when *K*_{D2} = 0.0123. It indicates an interesting phenomenon that the status of the rigid mode shapes remains invariant. The reason may be that the MSDS has the instinct of preserving its inherent characteristics. Moreover, it may infer that the mode shapes are more essential characteristics of the system than the DCSs.

**(a)**

**(b)**

**(a)**

**(b)**

When the other parameters remain the same, the effects of derivative coefficient with *c*_{m} = 2,400 N·s/m and *c*_{m} = 5,000 N·s/m are given in Figures 12 and 13, respectively. In the cases of *c*_{m} = 2,400 N·s/m and *c*_{m} = 5,000 N·s/m, the variation trends of DCS characteristics with the derivative coefficient are similar to those of the case *c*_{m} = 0. It can be observed that the interbearing damping mainly affects the dynamic characteristics of the third order. More specifically, the interbearing damping increases the third-order DCS, and there is an appropriate interbearing damping for the maximum of the logarithmic decrement curve. Moreover, intersections also exist between the DCS curves of the first two orders in Figures 12(a) and 13(a), and their corresponding mode shapes remain unchanged with the derivative coefficient, which are similar to those of the case *c*_{m} = 0 and are, therefore, not shown.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

Another interesting observation in Figures 9, 12, and 13 is that the first two order DCSs disappear when the derivative coefficient is greater than 0.015. It is due to the fact that the increased equivalent damping resulted from the derivative coefficient critically damps the rigid modal mode. It means that the MSDS rigid modes can be eliminated by increasing the equivalent damping of the AMBs, which makes the system response more stable when the MSDS operates under the first-order bending DCS.

##### 3.3. Effect of Proportional Coefficient

In this part, two cases corresponding to *K*_{D1} = 0.0047 and *K*_{D2} = 0.0123 before and after the intersection in Figure 9(a) are selected to investigate the effects of proportional coefficient on the dynamic characteristics with *K*_{I} = 0 and *c*_{m} = 2,400 N·s/m, shown in Figures 14 and 15, respectively.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

In Figure 14, with the increase of the proportional coefficient, the first two order DCSs increase while that of the third order barely varies, and the absolute values of damping exponents of the first two orders decrease while that of the third order slightly increases. Thus, the logarithmic decrements of the first two orders decrease and that of third order is nearly the same according to equation (16). In Figure 15 with *K*_{D2} = 0.0123, there are two differences from Figure 14. One is that there is also an intersection between the DCS curves of the first two orders in Figure 15(a). Likewise, the mode shape characteristics corresponding to the points before and after the intersection remain the same with the variation of the proportional coefficient, which are also similar to those of the cases above and are, therefore, not shown. The other is that the damping exponent of the third order decreases in Figure 15(b), but its amplitude variation is very small. The variation trends of DCSs and logarithmic decrements with proportional coefficient are, respectively, similar to those of the case *K*_{D1} = 0.0047.

##### 3.4. Effect of Integral Coefficient

Similarly, two cases corresponding to *K*_{D1} = 0.0047 and *K*_{D2} = 0.0123 before and after the intersection in Figure 9(a) are selected to illustrate the effects of integral coefficient on the dynamic characteristics with *K*_{P} = 5 and *c*_{m} = 2,400 N·s/m, shown in Figures 16 and 17, respectively.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

In Figure 16, with the increase of integral coefficient, the first three order DCSs almost remain invariant, and the absolute values of damping exponents of the first two orders decrease while that of the third order keeps unchanged and it is represented by the blue “o” in Figure 16(b). Thus, based on equation (16), the logarithmic decrements of the first two orders decrease, and that of third order stays the same. Comparing with *K*_{D1} = 0.0047 in Figure 16, the first two order DCSs slightly increase with integral coefficient in the case of *K*_{D2} = 0.0123 in Figure 17. In addition, there is also an intersection between the DCS curves of the first two orders in Figure 17(a), and the mode shape characteristics corresponding to the points before and after the intersection remain the same with different integral coefficients.

#### 4. Verification of Results

In order to verify the above analysis, the validity of the model and the characteristic analysis are carried out. The model verification is performed by the DCSs of Lund’s rotor [3] and the comparison between the MSDS’s calculation results of the FEM and the proposed model, and the verification of the characteristic analysis is demonstrated by the simulation results of the transient responses of the MSDS.

##### 4.1. Damped Critical Speeds of Lund’s Rotor

As stated above, a dual-rotor system gets degenerated into two straight rotors by setting all parameters related to interbearings equal to zero. The model is validated by the investigation of DCSs of Lund’s rotor [3], a uniform shaft supported at the ends in two damped bearings. The shaft has a length of 1.27 m, a diameter of 0.1016 m, Young’s modulus of 2.06 × 10^{11} Pa, and a weight density of 7,833.6 kg/m^{3}. The two bearings are identical with the stiffness coefficient 3.5 × 10^{6} N/m and the damping coefficient from 0 to 3.5 × 10^{5} N·s/m. The system is modelled as 41 lumped masses and 40 massless shaft elements. The DCS curves are depicted in Figure 18, which are in agreement with Lund’s results. Notice that the first and second modes become critically damped with the increase of bearing damping.

##### 4.2. Comparison of Results of the Proposed Model and FEM Model

It is necessary to verify the dynamic model of the MSDS as the modelling of the dual-rotor structure is simplified. The model is verified by FEM due to its high precision [8, 38]. The FEM model is constructed with elements BEAM 188 and COMBI 214 in ANSYS 15.0 code package. The BEAM 188 elements are used to model the shafts, and their sections and node positions are obtained according to structure size of the MSDS in Figure 3(a). The COMBI 214 elements are adopted to model the AMBs and the interbearings, and the stiffness and damping values of the element COMBI 214 correspond to the stiffness and damping values of the AMBs and the interbearings. For the AMBs, one end of the element COMBI 214 is linked with the node of the outer rotor and the other is completely constrained. For the interbearings, one end of the element COMBI 214 is connected with the node of the inner rotor and the other is connected with the node of the outer rotor. As a result, there are 337 nodes and 337 elements in total, and the FEM model of the MSDS is illustrated in Figure 19.

In the case of interbearing stiffness *k*_{m} = 1 × 10^{7} N/m, damping *c*_{m} = 0, and control parameters *K*_{P} = 8, *K*_{I} = 0, *K*_{D} = 0.01, and *r* = 1.5, the first three order DCSs and corresponding logarithmic decrements of the MSDS are listed in Table 3, and the unbalanced responses of disk 1 and disk 3 are given in Figure 20 with disk 1 being applied an unbalance of 1 × 10^{−5} kg·m. It can be seen from Figure 20 and Table 3 that the local response curves near the first two order DCSs are smooth, while that of the third order is steep, because the logarithmic decrements of the first two orders are much larger than 1, while that of the third order is much smaller. These results show that agreements between the results of Riccati transfer matrix method and FEM are achieved, which illustrates the validity of the MSDS dynamic model.

**(a)**

**(b)**

##### 4.3. Transient Response of the MSDS

It is known from Figure 9(a) that the first two order DCSs exist when the derivative coefficient is small (for example, *K*_{D} = 0.0003), and they disappear when the derivative coefficient is large (for example, *K*_{D} = 0.02). In order to examine the dynamic characteristics, the transient response of the MSDS is simulated by the combination of transfer matrix method and Newmark-*β* method.

With *K*_{P} = 5, *K*_{I} = 0, *k*_{m} = 1 × 10^{7} N/m, and *c*_{m} = 0 and *r* = 1.5, when the inner rotor of the MSDS moves uniformly from 0 to 10000 r/min by the exponential law *ω*_{i}(*t*) = 10000(1 − *e*^{−1.05t}) and an unbalance of 1 × 10^{−5} kg·m is applied on disk 1, the transient responses of disk 1 and disk 3 are given in Figures 21 and 22, respectively. In the case of *K*_{D} = 0.0003, it can be observed from Figures 21(a) and 22(a) that disk 1 and disk 3 simultaneously achieve two apparent peak amplitudes at specific speeds. Actually, the two speeds corresponding to the peaks in Figure 21(a) are the inner rotor speeds corresponding to the first two order DCSs of MSDS excited by the outer rotor. Similarly, the two distinct peaks in Figure 22(a) correspond to the first two order DCSs of MSDS excited by the outer rotor. In contrast, there are no peaks in both Figures 21(b) and 22(b) in the case of *K*_{D} = 0.02, which infer that the first two order DCSs disappear. Consequently, the response amplitudes are much smaller than those of the case *K*_{D} = 0.0003. Therefore, the analysis of the transient response demonstrates the validation of dynamic characteristics analysis in Section 3.2.

**(a)**

**(b)**

**(a)**

**(b)**

#### 5. Conclusion

In this research, a MSDS consisting of inner and outer rotors, two interbearings to connect the two rotors, and two AMBs to suspend the outer rotor is designed. With PID control adopted for the AMBs, the dynamic model of the MSDS is developed by the Riccati transfer matrix method using complex variables and the fully coupled unit model. On this basis, the dynamic characteristics of the MSDS are investigated.

According to the analysis, the dynamic characteristics of the MSDS rigid modes are more sensitive to the proportional and derivative coefficients and less sensitive to the integral coefficient. In contrast, the dynamic characteristics of first-order bending mode are insensitive to the PID control parameters. Moreover, two interesting phenomena are captured. The first one is that the mode shapes corresponding to the first two order DCSs still remain invariant even when the value of the first order DCS is greater than that of the second order. The second one is that the first two order DCSs disappear when the derivative coefficient is greater than a certain value, which indicates that the vibration response due to rigid modes can be eliminated or reduced by increasing the equivalent damping of AMBs. Finally, the dynamic model is verified by FEM and the second phenomenon is tested by the transient response simulation of the MSDS.

#### Appendix

In Section 2.2, the partitioned matrices used to formulate the overall frequency equation are described herein. Due to the limitation of the space, the mathematical details have been omitted. The related parameters of formulae are described in nomenclature:where *I*_{8 × 8} represents the 8 × 8 identity matrix.where *r*_{11} = *k*_{ixx} + *c*_{ixx}*s* + *k*_{cxx} + *c*_{cxx}*s*, *r*_{12} = *k*_{ixy} + *c*_{ixy}*s* + *k*_{cxy} + *c*_{cxy}*s*, *r*_{21} = *k*_{iyx} + *c*_{iyx}*s* + *k*_{cyx} + *c*_{cyx}*s*, *r*_{22} = *k*_{iyy} + *c*_{iyy}*s* + *k*_{cyy} + *c*_{cyy}*s*. *t*_{11} = *k*_{cxx} + *c*_{cxx}*s*, *t*_{12} = *k*_{cxy} + *c*_{cxy}*s*, *t*_{21} = *k*_{cyx} + *c*_{cyx}*s*, *t*_{22} = *k*_{cyy} + *c*_{cyy}*s*, = *k*_{oxx} + *c*_{oxx}*s* + *k*_{cxx} + *c*_{cxx}*s*, = *k*_{oxy} + *c*_{oxy}*s* + *k*_{cxy} + *c*_{cxy}*s*, = *k*_{oyx} + *c*_{oyx}*s* + *k*_{cyx} + *c*_{cyx}*s*, and = *k*_{oyy} + *c*_{oyy}*s* + *k*_{cyy} + *c*_{cyy}*s*.where *γ*_{i1} = *l*_{i}/(*EI*_{i}), , , *γ*_{o1} = *l*_{o}/(*EI*_{o}), , and .

#### Nomenclature

A: | Cross-sectional area of shaft, m^{2} |

B_{11}, B_{12}, B_{21}, B_{22}: | 8 × 8 partitioned matrices |

c_{ixx}, c_{ixy}, c_{iyx}, c_{iyy}: | Damping coefficients of bearings to support inner rotor, N·s/m; 2nd and 3rd indexes are force direction and velocity direction, respectively |

c_{oxx}, c_{oxy}, c_{oyx}, c_{oyy}: | Damping coefficients of bearings to support outer rotor, N·s/m, c_{cxx}, c_{cxy}, c_{cyx}, c_{cyy} = interbearing damping coefficients, N·s/m |

D_{11}, D_{12}, D_{21}, D_{22}: | 8 × 8 partitioned matrices |

e: | 8 × 1 displacement vector with constituents: x_{i}, θ_{yi}, y_{i}, θ_{xi}, x_{o}, θ_{yo}, y_{o}, θ_{xo} |

E: | Young’s modulus, Pa |

f: | 8 × 1 force vector with constituents: M_{yi}, Q_{xi}, M_{xi}, Q_{yi}, M_{yo}, Q_{xo}, M_{xo}, Q_{yo} |

G: | Shear modulus, Pa |

i: | Imaginary unit |

I_{i}, I_{o}: | Sectional moment of inertia of inner and outer rotors, respectively, m^{−4} |

J_{p}, J_{d}: | Polar and transverse moment of inertia, respectively, kg·m^{2} |

k_{ixx}, k_{ixy}, k_{iyx}, k_{iyy}: | Bearing stiffness coefficients for the inner rotor, N/m; 2nd and 3rd indexes are force direction and displacement direction, respectively |

k_{oxx}, k_{oxy}, k_{oyx}, k_{oyy}: | Bearing stiffness coefficients for the outer rotor, N/m |

k_{cxx}, k_{cxy}, k_{cyx}, k_{cyy}: | Interbearing stiffness coefficients, N/m |

l_{i}, l_{o}: | The inner and outer rotor lengths, respectively, m |

m_{i}, m_{o}: | The inner and outer rotor masses, respectively, kg |

M_{yi}, M_{xi}, M_{yo}, M_{xo}: | Bending moments, N·m, 1st index is moment direction |

Q_{yi}, Q_{xi}, Q_{yo}, Q_{xo}: | Shear force, N, 1st index is shearing force direction |

s: | λ + iω, the complex frequency, rad/s |

t: | Time, s |

u_{11}, u_{12}, u_{21}, u_{22}: | 8 × 8 partitioned matrices |

x_{i}, y_{i}, x_{o}, y_{o}: | Radial shaft displacements, m |

Z: | 16 × 1 state vector with state variables: M_{yi}, Q_{xi}, M_{xi}, Q_{yi}, M_{yo}, Q_{xo}, M_{xo}, Q_{yo}, x_{i}, θ_{yi}, y_{i}, θ_{xi}, x_{o}, θ_{yo}, y_{o}, θ_{xo} |

α: | Cross-sectional shape factor for shear deformation of the shaft |

Δ_{1}: | The frequency equation with singularity |

Δ_{p}: | The nonsingularity frequency equation |

θ_{xi}, θ_{xo}: | Angular displacement for inner and outer rotors in the Y-Z plane, rad |

θ_{yi}, θ_{yo}: | Angular displacement for inner and outer rotors in the X-Z plane, rad |

λ: | Real part of complex frequency s or damping exponent of vibrations of the shaft, s^{−1} |

μ: | Poisson’s ratio |

ρ: | Density of the shaft, kg/m^{3} |

ω: | Imaginary part of complex frequency s or damped natural frequency, rad/s |

Ω: | Angular speed of the shaft, rad/s Subscripts |

b: | A massless elastic beam |

c: | Interbearing |

d: | An element defined in Section 2 which consists of two lumped mass elements and three bearing elements |

i: | Inner rotor |

k: | Rotor station number |

o: | Outer rotor |

x: | X direction in the X-Z plane |

y: | Y direction in the Y-Z plane Superscripts |

L: | Left-hand side of an element |

R: | Right-hand side of an element. |

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (51405351), major projects of Natural Science Foundation of Hubei Province (2014CFA013), China Scholarship Council (no. 201808420108), and Hubei Youth Science and Technology Morning Light Project (2018B22).

#### Supplementary Materials

Figure 1: the schematic residual quantity curve of classic RTMM. Figure 2: the schematic residual quantity curve of the improved RTMM by the authors.* (Supplementary Materials)*