Research Article  Open Access
Dongxiong Wang, Nianxian Wang, Kuisheng Chen, Chun Ye, "Dynamic Characteristics of Magnetic Suspended DualRotor System by Riccati Transfer Matrix Method", Shock and Vibration, vol. 2019, Article ID 9843732, 22 pages, 2019. https://doi.org/10.1155/2019/9843732
Dynamic Characteristics of Magnetic Suspended DualRotor System by Riccati Transfer Matrix Method
Abstract
The magnetic suspended dualrotor 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 dualrotor system in the aeroengine, namely, the magnetic suspended dualrotor 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 bearingrotor 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 bearingrotor 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 bearingrotor 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 dualrotor 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 AMBrotor 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 hightemperature gascooled 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 AMBrotor system by Lagrange equations, discussed the dynamic characteristics of an energy storage flywheel AMBrotor 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 singlerotor systems. Due to the introduction of interbearings into the dualrotor 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 dualrotor 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 dualrotor structure, namely, the core component of the MSDS, is widely used in aeroengines. Researches on the dynamic characteristics of dualrotor systems are of great significance for lowering vibration and noise and enhancing their overall performance. Considerable theoretical work has been carried out on dualrotor 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 dualrotor turbomachinery 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 dualrotor system with FEM to explore the effects of interbearing stiffness on the critical speeds of the system when corotating and counterrotating. Wang et al. [25] established the dynamic model of a dualrotor 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 dualrotor systems were merely supported by rolling bearings with negligible damping. Although damping was taken into consideration in the dualrotor 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 AMBrotor 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 dualrotor 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 eightpole 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 airgap 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 AMBrotor 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 springdamping structures. The springdamping 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 dualrotor 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 crosscoupling 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 xz and yz 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 kth (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 kth element d is written aswhere and are the righthand and lefthand side state vectors of the kth 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 kth inner and outer rotor elements b are given aswhere and are the righthand and lefthand side state vectors of the kth element b, respectively. According to the geometric compatibility and equilibrium of internal forces at the junction between the kth elements d and b, it is required that , and the same for the junction between the kth element b and (k + 1)th element d is . Furthermore, omitting the subscripts b and d, the relationship of state vectors for the kth 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 lefthand 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 righthand side of station N, the relationship of the state vectors f_{N+1} and e_{N+1} is
With the righthand 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 oppositesigninfinitetype singularities in the residual quantity curve of Δ_{1}, and this makes the oppositesign intervals and rootexistence intervals not onetoone 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 oppositesigninfinitetype singularities to samesigninfinitetype 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 dualrotor are depicted in Figure 1(a), and the design parameters of AMBs are listed in Table 1.

3.1. Basic Dynamic Characteristics
In a dualrotor system, corotation and counterrotation are two possible operation modes. The former one is with both rotors corotating 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 corotation 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 firstorder 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 firstorder 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 firstorder mode shapes of both inner and outer rotors are translational modes, the secondorder conical modes, and the thirdorder 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 firstorder 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 thirdorder 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 firstorder 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 dualrotor 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 dualrotor 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 firstorder 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:  Crosssectional 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} 
α:  Crosssectional 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 YZ plane, rad 
θ_{yi}, θ_{yo}:  Angular displacement for inner and outer rotors in the XZ 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 XZ plane 
y:  Y direction in the YZ plane Superscripts 
L:  Lefthand side of an element 
R:  Righthand 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)
References
 G. Schweitzer and E. H. Maslen, Magnetic Bearings: Theory, Design, and Application to Rotating Machinery, Springer, Berlin, Germany, 2009.
 Z. Wu, J. Li, and R. Shi, “Current research status and key technology of moreelectric aeroengine,” Advances in Aeronautical Science and Engineering, vol. 3, no. 4, pp. 463–467, 2012. View at: Google Scholar
 J. W. Lund, “Stability and damped critical speeds of a flexible rotor in fluidfilm bearings,” Journal of Engineering for Industry, vol. 96, no. 2, pp. 509–517, 1974. View at: Publisher Site  Google Scholar
 P. N. Bansal and R. G. Kirk, “Stability and damped critical speeds of rotorbearing systems,” Journal of Engineering for Industry, vol. 97, no. 4, pp. 1325–1332, 1975. View at: Publisher Site  Google Scholar
 N. Wang, Y. Hu, H. Wu, and R. Dong, “Research on bearing capacity model of largeairgap hybrid magnetic bearings,” Journal of Mechanical Engineering, vol. 51, no. 1, pp. 153–160, 2015. View at: Publisher Site  Google Scholar
 X. Ren, Y. Le, J. Sun, B. Han, and K. Wang, “Magnetic flux leakage modelling and optimisation of a CRAHMB for DC motor,” IET Electric Power Applications, vol. 11, no. 2, pp. 212–221, 2017. View at: Publisher Site  Google Scholar
 D.P. Wang and F.X. Wang, “Levitation force analysis of magnetic bearing by circuitfield combination method,” Electric Machines and Control, vol. 15, no. 11, pp. 8–13, 2011. View at: Google Scholar
 D. Wang, N. Wang, C. Ye, and K. Chen, “Research on analytical bearing capacity model of active magnetic bearings based on magnetic saturation,” IET Electric Power Applications, vol. 11, no. 9, pp. 1548–1557, 2017. View at: Publisher Site  Google Scholar
 X. Wang and J. Zhu, “Investigation on dynamic performance of active magnetic bearing rotor system,” Chinese Journal of Mechanical Engineering, vol. 37, no. 11, pp. 7–12, 2001. View at: Publisher Site  Google Scholar
 R. R. Humphris, R. D. Kelm, D. W. Lewis, and P. E. Allaire, “Effect of control algorithms on magnetic journal bearing properties,” Journal of Engineering for Gas Turbines and Power, vol. 108, no. 4, pp. 624–632, 1986. View at: Publisher Site  Google Scholar
 Y. Guojun, X. Yang, S. Zhengang, and G. Huidong, “Characteristic analysis of rotor dynamics and experiments of active magnetic bearing for HTR10GT,” Nuclear Engineering & Design, vol. 237, no. 1213, pp. 1363–1371, 2007. View at: Publisher Site  Google Scholar
 J.N. Liu, Z.Y. Ren, S.W. Wu, and Y.L. Tang, “Theoretical vibration analysis on 600 Wh energy storage flywheel rotor—active magnetic bearing system,” International Journal of Rotating Machinery, vol. 2013, no. 8, Article ID 512674, 11 pages, 2013. View at: Publisher Site  Google Scholar
 J. Tang, B. Xiang, and Y. Zhang, “Dynamic characteristics of the rotor in a magnetically suspended control moment gyroscope with active magnetic bearing and passive magnetic bearing,” ISA Transactions, vol. 53, no. 4, pp. 1357–1365, 2014. View at: Publisher Site  Google Scholar
 K. Gupta, K. D. Gupta, and K. Athre, “Unbalance response of a dual rotor system theory and experiment,” Journal of Vibration and Acoustics, vol. 115, no. 4, pp. 427–435, 1993. View at: Publisher Site  Google Scholar
 X. Chen and M. Liao, “Steady state characteristics of a dualrotor system with intershaft bearing subjected to mass unbalance and base motions,” in Proceedings of ASME Turbo Expo: Turbomachinery Technical Conference and Exposition, pp. 1–16, Oslo, Norway, June 2018. View at: Google Scholar
 K. D. Gupta, K. Gupta, and K. Athre, “Stability analysis of dual rotor system by extended transfer matrix method,” in Proceedings of the ASME International Gas Turbine and Aeroengine Congress and Exposition, New York, NY, USA, June 1989. View at: Google Scholar
 X. Chen and M. Liao, “Transient characteristics of a dualrotor system with intershaft bearing subjected to mass unbalance and base motions during startup,” in Proceedings of ASME Turbo Expo: Turbomachinery Technical Conference and Exposition, pp. 1–18, Oslo, Norway, June 2018. View at: Google Scholar
 P. Yu, D. Zhang, Y. Ma, and J. Hong, “Dynamic modeling and vibration characteristics analysis of the aeroengine dualrotor system with fan blade out,” Mechanical Systems and Signal Processing, vol. 106, no. 1, pp. 158–175, 2018. View at: Publisher Site  Google Scholar
 Y. Yang, D. Cao, D. Wang, and G. Jiang, “Fixedpoint rubbing characteristic analysis of a dualrotor system based on the LankaraniNikravesh model,” Mechanism and Machine Theory, vol. 103, no. 1, pp. 202–221, 2016. View at: Publisher Site  Google Scholar
 Y. Yang, D. Cao, T. Yu, D. Wang, and C. Li, “Prediction of dynamic characteristics of a dualrotor system with fixed point rubbing—theoretical analysis and experimental study,” International Journal of Mechanical Sciences, vol. 115116, no. 1, pp. 253–261, 2016. View at: Publisher Site  Google Scholar
 F. Wang, G.H. Luo, X.G. Yang, and H.T. Cui, “Research on modeling and dynamic characteristics of complex coaxial rotor system,” Journal of Vibroengineering, vol. 19, no. 3, pp. 1524–1545, 2017. View at: Publisher Site  Google Scholar
 G. Luo, X. Hu, and X. Yang, “Nonlinear Analysis of counterrotating dualrotor system,” Journal of Vibration Engineering, vol. 22, no. 3, pp. 268–273, 2009. View at: Google Scholar
 H.W. D. Chiang, C.N. Hsu, and S.H. Tu, “Rotorbearing analysis for turbomachinery single and dualrotor systems,” Journal of Propulsion and Power, vol. 20, no. 6, pp. 1096–1104, 2004. View at: Publisher Site  Google Scholar
 Z.X. Fei, S.G. Tong, and C. Wei, “Investigation of the dynamic characteristics of a dual rotor system and its startup simulation based on finite element method,” Journal of Zhejiang University SCIENCE A, vol. 14, no. 4, pp. 268–280, 2013. View at: Publisher Site  Google Scholar
 N. Wang, D. Jiang, and H. Xu, “Dynamic characteristics analysis of a dualrotor system with intershaft bearing,” Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, vol. 233, no. 3, pp. 1147–1158, 2017. View at: Publisher Site  Google Scholar
 P. Bonello and P. M. Hai, “A receptance harmonic balance technique for the computation of the vibration of a whole aeroengine model with nonlinear bearings,” Journal of Sound & Vibration, vol. 324, no. 12, pp. 221–242, 2009. View at: Publisher Site  Google Scholar
 G. Sun, A. Palazzolo, A. Provenza, C. Lawrence, and K. Carney, “Long duration blade loss simulations including thermal growths for dualrotor gas turbine engine,” Journal of Sound & Vibration, vol. 316, no. 15, pp. 147–163, 2008. View at: Publisher Site  Google Scholar
 P. M. Hai and P. Bonello, “An impulsive receptance technique for the time domain computation of the vibration of a whole aeroengine model with nonlinear bearings,” Journal of Sound & Vibration, vol. 318, no. 3, pp. 592–605, 2016. View at: Publisher Site  Google Scholar
 S. Zheng, B. Han, Y. Wang, and J. Zhou, “Optimization of damping compensation for a flexible rotor system with active magnetic bearing considering gyroscopic effect,” IEEE/ASME Transactions on Mechatronics, vol. 20, no. 3, pp. 1130–1137, 2015. View at: Publisher Site  Google Scholar
 H. K. Roy, A. S. Das, and J. K. Dutt, “An efficient rotor suspension with active magnetic bearings having viscoelastic control law,” Mechanism and Machine Theory, vol. 98, no. 1, pp. 48–63, 2016. View at: Publisher Site  Google Scholar
 R. Ebrahimi, M. Ghayour, and H. M. Khanlo, “Chaotic vibration analysis of a coaxial rotor system in active magnetic bearings and contact with auxiliary bearings,” Journal of Computational and Nonlinear Dynamics, vol. 12, no. 3, Article ID 031012, 2016. View at: Publisher Site  Google Scholar
 R. Ebrahimi, M. Ghayour, and H. M. Khanlo, “Effects of some design parameters on bifurcation behavior of a magnetically supported coaxial rotor in auxiliary bearings,” Engineering Computations, vol. 34, no. 7, pp. 2379–2395, 2017. View at: Publisher Site  Google Scholar
 K. Jiang, C. Zhu, L. Chen, and X. Qiao, “MultiDOF rotor model based measurement of stiffness and damping for active magnetic bearing using multifrequency excitation,” Mechanical Systems and Signal Processing, vol. 6061, no. 1, pp. 358–374, 2015. View at: Publisher Site  Google Scholar
 Z. Wang, C. Mao, and C. Zhu, “A design method of PID controllers for active magnetic bearingsrigid rotor systems,” Proceedings of the Chinese Society for Electrical Engineering, vol. 38, no. 20, pp. 1–9, 2018. View at: Google Scholar
 P. Anantachaisilp and Z. Lin, “An experimental study on PID tuning methods for active magnetic bearing systems,” International Journal of Advanced Mechatronic Systems, vol. 5, no. 2, pp. 146–154, 2013. View at: Publisher Site  Google Scholar
 P. Anantachaisilp and Z. Lin, “Fractional order PID control of rotor suspension by active magnetic bearings,” Actuators, vol. 6, no. 1, pp. 1–31, 2017. View at: Publisher Site  Google Scholar
 Z. Wang, “Singularities of the Riccati transfer matrix method and its elimination,” Journal of Vibration and Shock, vol. 2, no. 1, pp. 74–78, 1987. View at: Google Scholar
 N. Wang, D. Wang, K. Chen, and H. Wu, “Research on analytical model and design formulas of permanent magnetic bearings based on Halbach array with arbitrary segmented magnetized angle,” Journal of Magnetism and Magnetic Materials, vol. 410, no. 1, pp. 257–264, 2016. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Dongxiong Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.