Shock and Vibration
Volume 2019 (2019), Article ID 9843732, 22 pages
https://doi.org/10.1155/2019/9843732

Research Article

Dynamic Characteristics of Magnetic Suspended Dual-Rotor System by Riccati Transfer Matrix Method

Dongxiong Wang,1,2 Nianxian Wang,1,2 Kuisheng Chen,1,2 and Chun Ye1,2

¹The Key Laboratory for Metallurgical Equipment and Control of Ministry of Education, Wuhan University of Science and Technology, No. 947, Heping Venue, Qingshan District, CT, Wuhan 430081, China
²Hubei Key Laboratory of Mechanical Transmission and Manufacturing Engineering, Wuhan University of Science and Technology, No. 947, Heping Avenue, Qingshan District, CT, Wuhan 430081, China

Correspondence should be addressed to Nianxian Wang; wangnianxian1987@163.com

Received 12 March 2019; Revised 24 May 2019; Accepted 9 June 2019; Published 24 July 2019

Academic Editor: Giuseppe Petrone

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.

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)

Figure 1: (a) Structure and (b) theoretical model of the MSDS.

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 ke and equivalent damping de of the AMB can be expressed as

where KP, KI, and KD are the proportional, integral, and derivative coefficients, respectively, ω is the whirl frequency of the system, Aa the gain of the power amplifier, As is the displacement sensor, and ks and ki are the displacement stiffness coefficient and current stiffness coefficient of the AMB, and they are, respectively, given as

where μ₀ is the permeability of vacuum, α the pole face angle, Ap the sectional area of one stator pole, N the coil turn, Ib the bias current of coils, and the nominal air-gap length.

(a)

(b)

Figure 2: (a) Basic control loop and (b) control block diagram of the AMB.

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 (kxx, kxy, kyx, kyy) of which represent the stiffness properties and the other four (cxx, cxy, cyx, cyy), the damping properties, where kxy = kyx = 0 and cxy = cyx = 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 = {MyiQxiMxiQyiMyoQxoMxoQyo}T and e = {xiθyiyiθxixoθyoyoθ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)

Figure 3: (a) Fully coupled unit model and (b) its sign convention of the MSDS.

In the fully coupled unit model, it is assumed that the inner and outer rotors are divided into many elements ai, bi, ao and bo, 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 ai, an outer rotor element ao, an interbearing element cc which interlinks the inner rotor element ai and the outer rotor element ao, and two elements ci and co to only, respectively, support the inner rotor element ai and the outer rotor element ao 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 bi, an outer rotor element bo, 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 as

where and are the right-hand and left-hand 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₁₁, D₁₂, D₂₁, and D₂₂ 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 as

where and are the right-hand and left-hand 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 as

where

The partitioned matrices B₁₁, B₁₂, B₂₁, B₂₂, D₁₁, D₁₂, D₂₁, and D₂₂ are given in Appendix. Omitting the superscript L, the expansion of equation (5) gives

Introducing a Riccati transformation at station k, we get

where 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₁ = 0 and e₁ ≠ 0, the initial condition is S₁ = 0. Knowing u₁₁, u₁₂, u₂₁, and u₂₂ for each fully coupled unit, Sk + 1 can be obtained sequentially by using equation (10) repeatedly. For the right-hand side of station N, the relationship of the state vectors fN+1 and eN+1 is

With the right-hand side boundary conditions fN+1 = 0 and eN+1 ≠ 0, a nontrivial solution of equation (10) should satisfy the condition that the determinant of SN+1 is zero, namely,

However, there are many opposite-sign-infinite-type singularities in the residual quantity curve of Δ₁, 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 as

where 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 as

where

where ω (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.

Table 1: Model parameters of the AMB system.

Parameter
Description
Value

μ₀
Vacuum permeability
4π × 10⁻⁷ H/m

Fmax
Maximum bearing capacity
240 N

Nominal air-gap length
0.5 mm

N
Coil turn
160

Ib
Bias current
3.0 A

icmax
Maximum control current
3.0 A

Ap
Sectional area of stator pole
2.47 × 10⁻⁴ m²

Α
Half angle between two poles
π/8 rad

Aa
Gain of the power amplifier
0.6 A/V

As
Gain of the displacement sensor
8 × 10³ V/m

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 km = 1 × 10⁷ N/m, damping cm = 0, and control parameters KP = 8, KI = 0, KD = 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) Campbell diagram and (b) logarithmic decrements of the dual-rotor system.

Table 2: Damped critical speeds (r/min).

Order
Inner rotor excitation
Outer rotor excitation

1
8132
8132

2
9763
9356

3
31315
31203

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)

Figure 5: Mode shapes of the dual-rotor system (—, inner rotor; ---, outer rotor). Mode shapes excited by the (a) inner rotor and (b) outer rotor.

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 KD = 0.01 and KI = 0 is given in Figure 6, where the subscript number in KP1, KP2, and KP3 means the corresponding order of whirl frequency, and the same for that of KD1, KD2, and KD3 in Figure 7 and KI1, KI2, and KI3 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 KP, 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)

Figure 6: The effect of proportional coefficient on the MSDS dynamic characteristics. (a) Campbell diagram. (b) Logarithmic decrements.

(a)

(b)

Figure 7: The effect of derivative coefficient on the MSDS dynamic characteristics. (a) Campbell diagram. (b) Logarithmic decrements.

(a)

(b)

Figure 8: The effect of integral coefficient on the MSDS dynamic characteristics. (a) Campbell diagram. (b) Logarithmic decrements.

The effect of derivative coefficient with KP = 8 and KI = 0 is given in Figure 7. Increasing the derivative coefficient KD, 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 KP = 8 and KD = 0.01 is given in Figure 8. Increasing the integral coefficient KI, 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 KI = 0, KP = 5 and cm = 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)

Figure 9: The effect of derivative coefficient on DCS characteristics (cm = 0). (a) Derivative coefficient. (b) Damping exponent (Hz). (c) Derivative coefficient. (d) Derivative coefficient.

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, KD1 = 0.0047 and KD2 = 0.0123 shown in Figure 9(a), are investigated. The first three mode shapes of MSDS corresponding to KD1 = 0.0047 and KD2 = 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 KD2 = 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)

Figure 10: Mode shapes of MSDS with KD1 = 0.0047 (—, inner rotor; ---, outer rotor). Mode shapes excited by the (a) inner rotor and (b) outer rotor.

(a)

(b)

Figure 11: Mode shapes of MSDS with KD2 = 0.0123(—, inner rotor and ---, outer rotor). Mode shapes excited by the (a) inner rotor and (b) outer rotor.

When the other parameters remain the same, the effects of derivative coefficient with cm = 2,400 N·s/m and cm = 5,000 N·s/m are given in Figures 12 and 13, respectively. In the cases of cm = 2,400 N·s/m and cm = 5,000 N·s/m, the variation trends of DCS characteristics with the derivative coefficient are similar to those of the case cm = 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 cm = 0 and are, therefore, not shown.

(a)

(b)

(c)

(d)

Figure 12: The effect of derivative coefficient on DCS characteristics (cm = 2,400 N·s/m). (a) Derivative coefficient. (b) Damping exponent (Hz). (c) Derivative coefficient. (d) Derivative coefficient.

(a)

(b)

(c)

(d)

Figure 13: The effect of derivative coefficient on DCS characteristics (cm = 5,000 N·s/m). (a) Derivative coefficient. (b) Damping exponent (Hz). (c) Derivative coefficient. (d) Derivative coefficient.

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 KD1 = 0.0047 and KD2 = 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 KI = 0 and cm = 2,400 N·s/m, shown in Figures 14 and 15, respectively.

(a)

(b)

(c)

Figure 14: The effect of proportional coefficient on DCS characteristics (KD1 = 0.0047). (a) Proportional coefficient. (b) Damping exponent (Hz). (c) Proportional coefficient.

(a)

(b)

(c)

Figure 15: The effect of proportional coefficient on DCS characteristics (KD2 = 0.0123). (a) Proportional coefficient. (b) Damping exponent (Hz). (c) Proportional coefficient.

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 KD2 = 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 KD1 = 0.0047.

3.4. Effect of Integral Coefficient

Similarly, two cases corresponding to KD1 = 0.0047 and KD2 = 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 KP = 5 and cm = 2,400 N·s/m, shown in Figures 16 and 17, respectively.

(a)

(b)

(c)

Figure 16: The effect of integral coefficient on DCS characteristics (KD1 = 0.0047). (a) Integral coefficient. (b) Damping exponent (Hz). (c) Integral coefficient.

(a)

(b)

(c)

Figure 17: The effect of integral coefficient on DCS characteristics (KD2 = 0.0123). (a) Integral coefficient. (b) Damping exponent (Hz). (c) Integral coefficient.

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 KD1 = 0.0047 in Figure 16, the first two order DCSs slightly increase with integral coefficient in the case of KD2 = 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¹¹ Pa, and a weight density of 7,833.6 kg/m³. The two bearings are identical with the stiffness coefficient 3.5 × 10⁶ N/m and the damping coefficient from 0 to 3.5 × 10⁵ 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.

Figure 18: DCSs of Lund’s rotor.

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.

Figure 19: The FEM model of the MSDS.

In the case of interbearing stiffness km = 1 × 10⁷ N/m, damping cm = 0, and control parameters KP = 8, KI = 0, KD = 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⁻⁵ 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.

Table 3: Results comparison between FEM and the Riccati transfer matrix method.

Order
Damped critical speeds (Hz)
Logarithmic decrements

Riccati
FEM
Error (%)
Riccati
FEM
Error (%)

1
135.53
133.16
1.78
4.05
3.90
3.85

2
162.72
143.65
12.28
4.86
4.79
1.46

3
521.92
514.41
1.46
0.093
0.081
14.82

(a)

(b)

Figure 20: Comparison of unbalance responses between FEM and the proposed model of the MSDS. Unbalance response of (a) disk 1 and (b) disk 3.

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, KD = 0.0003), and they disappear when the derivative coefficient is large (for example, KD = 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 KP = 5, KI = 0, km = 1 × 10⁷ N/m, and cm = 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⁻⁵ 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 KD = 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 KD = 0.02, which infer that the first two order DCSs disappear. Consequently, the response amplitudes are much smaller than those of the case KD = 0.0003. Therefore, the analysis of the transient response demonstrates the validation of dynamic characteristics analysis in Section 3.2.

(a)

(b)

Figure 21: Transient response of disk 1. (a) KD = 0.0003. (b) KD = 0.02.

(a)

(b)

Figure 22: Transient response of disk 3. (a) KD = 0.0003. (b) KD = 0.02.

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 I8 × 8 represents the 8 × 8 identity matrix.

where r₁₁ = kixx + cixxs + kcxx + ccxxs, r₁₂ = kixy + cixys + kcxy + ccxys, r₂₁ = kiyx + ciyxs + kcyx + ccyxs, r₂₂ = kiyy + ciyys + kcyy + ccyys. t₁₁ = kcxx + ccxxs, t₁₂ = kcxy + ccxys, t₂₁ = kcyx + ccyxs, t₂₂ = kcyy + ccyys,  = koxx + coxxs + kcxx + ccxxs,  = koxy + coxys + kcxy + ccxys,  = koyx + coyxs + kcyx + ccyxs, and  = koyy + coyys + kcyy + ccyys.

where γi1 = li/(EIi), , , γo1 = lo/(EIo), , and .

Nomenclature

A: | Cross-sectional area of shaft, m²
B₁₁, B₁₂, B₂₁, B₂₂: | 8 × 8 partitioned matrices
cixx, cixy, ciyx, ciyy: | Damping coefficients of bearings to support inner rotor, N·s/m; 2nd and 3rd indexes are force direction and velocity direction, respectively
coxx, coxy, coyx, coyy: | Damping coefficients of bearings to support outer rotor, N·s/m, ccxx, ccxy, ccyx, ccyy = interbearing damping coefficients, N·s/m
D₁₁, D₁₂, D₂₁, D₂₂: | 8 × 8 partitioned matrices
e: | 8 × 1 displacement vector with constituents: xi, θyi, yi, θxi, xo, θyo, yo, θxo
E: | Young’s modulus, Pa
f: | 8 × 1 force vector with constituents: Myi, Qxi, Mxi, Qyi, Myo, Qxo, Mxo, Qyo
G: | Shear modulus, Pa
i: | Imaginary unit
Ii, Io: | Sectional moment of inertia of inner and outer rotors, respectively, m⁻⁴
Jp, Jd: | Polar and transverse moment of inertia, respectively, kg·m²
kixx, kixy, kiyx, kiyy: | Bearing stiffness coefficients for the inner rotor, N/m; 2nd and 3rd indexes are force direction and displacement direction, respectively
koxx, koxy, koyx, koyy: | Bearing stiffness coefficients for the outer rotor, N/m
kcxx, kcxy, kcyx, kcyy: | Interbearing stiffness coefficients, N/m
li, lo: | The inner and outer rotor lengths, respectively, m
mi, mo: | The inner and outer rotor masses, respectively, kg
Myi, Mxi, Myo, Mxo: | Bending moments, N·m, 1st index is moment direction
Qyi, Qxi, Qyo, Qxo: | Shear force, N, 1st index is shearing force direction
s: | λ + iω, the complex frequency, rad/s
t: | Time, s
u₁₁, u₁₂, u₂₁, u₂₂: | 8 × 8 partitioned matrices
xi, yi, xo, yo: | Radial shaft displacements, m
Z: | 16 × 1 state vector with state variables: Myi, Qxi, Mxi, Qyi, Myo, Qxo, Mxo, Qyo, xi, θyi, yi, θxi, xo, θyo, yo, θxo
α: | Cross-sectional shape factor for shear deformation of the shaft
Δ₁: | 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⁻¹
μ: | Poisson’s ratio
ρ: | Density of the shaft, kg/m³
ω: | 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)

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 more-electric aeroengine,” Advances in Aeronautical Science and Engineering, vol. 3, no. 4, pp. 463–467, 2012.

J. W. Lund, “Stability and damped critical speeds of a flexible rotor in fluid-film bearings,” Journal of Engineering for Industry, vol. 96, no. 2, pp. 509–517, 1974.

P. N. Bansal and R. G. Kirk, “Stability and damped critical speeds of rotor-bearing systems,” Journal of Engineering for Industry, vol. 97, no. 4, pp. 1325–1332, 1975.

N. Wang, Y. Hu, H. Wu, and R. Dong, “Research on bearing capacity model of large-air-gap hybrid magnetic bearings,” Journal of Mechanical Engineering, vol. 51, no. 1, pp. 153–160, 2015.

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.

D.-P. Wang and F.-X. Wang, “Levitation force analysis of magnetic bearing by circuit-field combination method,” Electric Machines and Control, vol. 15, no. 11, pp. 8–13, 2011.

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.

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.

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.

Y. Guojun, X. Yang, S. Zhengang, and G. Huidong, “Characteristic analysis of rotor dynamics and experiments of active magnetic bearing for HTR-10GT,” Nuclear Engineering & Design, vol. 237, no. 12-13, pp. 1363–1371, 2007.

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.

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.

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.

X. Chen and M. Liao, “Steady state characteristics of a dual-rotor 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.

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.

X. Chen and M. Liao, “Transient characteristics of a dual-rotor system with intershaft bearing subjected to mass unbalance and base motions during start-up,” in Proceedings of ASME Turbo Expo: Turbomachinery Technical Conference and Exposition, pp. 1–18, Oslo, Norway, June 2018.

P. Yu, D. Zhang, Y. Ma, and J. Hong, “Dynamic modeling and vibration characteristics analysis of the aero-engine dual-rotor system with fan blade out,” Mechanical Systems and Signal Processing, vol. 106, no. 1, pp. 158–175, 2018.

Y. Yang, D. Cao, D. Wang, and G. Jiang, “Fixed-point rubbing characteristic analysis of a dual-rotor system based on the Lankarani-Nikravesh model,” Mechanism and Machine Theory, vol. 103, no. 1, pp. 202–221, 2016.

Y. Yang, D. Cao, T. Yu, D. Wang, and C. Li, “Prediction of dynamic characteristics of a dual-rotor system with fixed point rubbing—theoretical analysis and experimental study,” International Journal of Mechanical Sciences, vol. 115-116, no. 1, pp. 253–261, 2016.

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.

G. Luo, X. Hu, and X. Yang, “Nonlinear Analysis of counter-rotating dual-rotor system,” Journal of Vibration Engineering, vol. 22, no. 3, pp. 268–273, 2009.

H.-W. D. Chiang, C.-N. Hsu, and S.-H. Tu, “Rotor-bearing analysis for turbomachinery single- and dual-rotor systems,” Journal of Propulsion and Power, vol. 20, no. 6, pp. 1096–1104, 2004.

Z.-X. Fei, S.-G. Tong, and C. Wei, “Investigation of the dynamic characteristics of a dual rotor system and its start-up simulation based on finite element method,” Journal of Zhejiang University SCIENCE A, vol. 14, no. 4, pp. 268–280, 2013.

N. Wang, D. Jiang, and H. Xu, “Dynamic characteristics analysis of a dual-rotor system with inter-shaft bearing,” Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, vol. 233, no. 3, pp. 1147–1158, 2017.

P. Bonello and P. M. Hai, “A receptance harmonic balance technique for the computation of the vibration of a whole aero-engine model with nonlinear bearings,” Journal of Sound & Vibration, vol. 324, no. 1-2, pp. 221–242, 2009.

G. Sun, A. Palazzolo, A. Provenza, C. Lawrence, and K. Carney, “Long duration blade loss simulations including thermal growths for dual-rotor gas turbine engine,” Journal of Sound & Vibration, vol. 316, no. 1-5, pp. 147–163, 2008.

P. M. Hai and P. Bonello, “An impulsive receptance technique for the time domain computation of the vibration of a whole aero-engine model with nonlinear bearings,” Journal of Sound & Vibration, vol. 318, no. 3, pp. 592–605, 2016.

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.

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.

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.

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.

K. Jiang, C. Zhu, L. Chen, and X. Qiao, “Multi-DOF rotor model based measurement of stiffness and damping for active magnetic bearing using multi-frequency excitation,” Mechanical Systems and Signal Processing, vol. 60-61, no. 1, pp. 358–374, 2015.

Z. Wang, C. Mao, and C. Zhu, “A design method of PID controllers for active magnetic bearings-rigid rotor systems,” Proceedings of the Chinese Society for Electrical Engineering, vol. 38, no. 20, pp. 1–9, 2018.

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.

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.

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.

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.