# Theoretical Vibration Analysis on 600Wh Energy Storage Flywheel Rotor—Active Magnetic Bearing System

Authors: Jing-na Liu; Zheng-yi Ren; Shan-wei Wu; Yin-long Tang
Journal: International Journal of Rotating Machinery
 (2013)

---

## Abstract

This paper shows a theoretical vibration analysis regarding the controller’s parameters and the gyroscopic effect, based on a simplified rotordynamic model. Combined with 600 Wh energy storage flywheel rotor system mathematical model, the Campbell diagram of the rotor system was obtained by the calculation of the whirl frequency under different parameters of the controller in MATLAB to analyze the effect of the controller parameter on the whirl frequency and to limit the operating speed and acceleration or deceleration of the rotor. The result of the analysis can be used to set the support position of the rotor system, limit the ratio of transverse moment of inertia and the polar moment of inertia, and direct the flywheel prototype future design. The presented simplified rotordynamic model can also be applied to rotating machines.


---

## Body

## 1. Introduction

Later in the 1970s, flywheel energy storage was proposed as a primary objective for electric vehicles and stationary power backup [1]. With the improvements in materials, magnetic bearing technology, and power electronics, flywheel energy storage technology has large developments. Compared with traditional battery energy storage system, flywheel energy storage system has many advantages such as higher energy storage density, higher specific power and power density, lower risk of overcharge or overdischarge, wide range of operation temperature, very long life cycle, and environmental friendliness [2].Many problems appear as the development of flywheel energy storage, and one of them is the bearing. Besides, the active magnetic bearing (AMB) implies that bearing forces are actively controlled by means of electromagnets, a well-designed closed control loop, and other components such as position sensors and power amplifiers [3]. Therefore, the rotor of the AMB can be suspended to the predefined positions by the controlled electromagnetic forces without mechanical contact and friction between the magnetic bearing and the rotor [4]. Based on the noncontact and frictionless characteristics, the magnetic suspension of AMB offers many practical and promising advantages over conventional bearings such as longer life, lower rotating frictional losses, higher rotational speed, and elimination of the lubrication [5]. Hence, AMBs have been successfully and widely implemented in various high-performance applications including the rotating devices such as turbine engines [6], flywheel energy and storage devices [7], bearingless motor [8], and vacuum pump [9] and the nonrotating devices such as motion control stage [10], biomedical applications [11], and manufacturing equipment [12]. Furthermore, the adjustable stiffness and damping characteristics make the AMB suitable for elimination of vibration [13].Many different kinds of excitations exist in rotor system—for example, mechanical unbalance and misalignment of the coupling—which may cause vibrations [14–18]. Besides these typical excitations, also specific excitations associated with the type of the rotating machine occur. This paper focuses on vibrations of rotor system supported by active electromagnet bearing; therefore, also electromagnetic forces have to be considered, which may cause vibrations by an eccentricity. The eccentricity can be divided into static eccentricity and dynamic eccentricity. Static eccentricity is, for example, caused by production tolerances regarding concentricity and fitting tolerance between stator housing, bearing housing, and so on. Dynamic eccentricity is caused if the rotor is bent or if the rotor core is eccentrically positioned on the rotor shaft and so on.The aim of this paper is to show a theoretical vibration analysis of the rotor system of 600 Wh flywheel energy storage system, based on a simplified rotordynamic model. The result of the analysis can be used to set the support position of the rotor system, limit the ratio of transverse moment of inertia and the polar moment of inertia, and direct the flywheel prototype future design.


## 2. System Modeling

### 2.1. System Structure

The basic layout of a flywheel energy storage system is depicted in Figure1, that is, the 600 Wh prototype system, a rigid vertical rotor shaft with a rigid motor and flywheel connected to shaft. Flywheel energy storage system is a complex construction where energy is stored mechanically and transferred to and from the flywheel by an integrated motor/generator.Figure 1
Basic layout of a flywheel energy storage system.Two radial active electromagnetic bearings (RAMBs) including upper and lower RAMBs (RAMB1 and RAMB2 in Figure1) and one thrust magnetic bearing (TAMB) are fixed on the platform to suspend and regulate the rotor in the radial and axial DOF, respectively. The positions of the rotor in five axes are defined as the displacements deviated from the nominal air gaps. Therefore, two pairs of the perpendicular eddy-current position sensors are installed closely to respective RAMBs to measure the respective radial positions in x- and y-axes which are denoted as x1 and y1 of the upper RAMB and x2 and y2 of the lower RAMB. Moreover, one eddy-current position sensor is installed closely to TAMB to measure the axial position in the z-axis which is denoted as z. After the rotor positions in five axes are all measured and sent to the control core through the position signals line, the required electromagnet currents are generated from the drive system and they circulate the coils through the power line. Therefore, the rotor can be regulated and stabilized in the centers of the apertures of the two RAMBs and the thrust disc can be centered in the middle of the air gap of the TAMB, respectively.One 70 kg flywheel is mounted on the rotor whose position can be adjusted to modify the rotational characteristics. Furthermore, one motor/generator is equipped between the lower RAMB and flywheel and is responsible for the torque generation or electricity generation based on the rotational speed requirement. However, in the rotor system, the rolling-element auxiliary bearings (auxiliary bearing 1 and auxiliary bearing 2 in Figure1) are necessary to protect AMB stators and stationary components along shaft in the even t of AMB failure or high transient loads. Under normal operation of AMBs, the rotor maintains a positive clearance with auxiliary bearings, which is less than a clearance with AMBs.


### 2.2. AMB Modeling

#### 2.2.1. Electromagnetic Force Model

Considering the diameter of the shaft and the loading capacity, the RAMB with 8-pole legs was designed in the 600 Wh prototype system. The structural configuration of the RAMB is shown in Figure2(a), including electromagnetic coils. And the basic magnetic bearing control loop is shown in Figure 2(b).Structural configurations and the basic magnetic bearing control loop of the RAMB.
(a)
Structural configuration of the RAMB with 8-pole legs
(b)
The basic magnetic bearing control loop and its elementsTaking the two pairs of poles in they-axes of one RAMB, for example, the electromagnetic force model will be established [19]. As shown in Figure 2(b), the basic magnetic bearing control loop includes differential driving, position sensor, controller, and power amplifier. The included differential driving is adopted in this study to obtain maximum range of the force dynamic and good linearity of the control dynamic. α is the half of the angle between the two-pole legs, y0 is the nominal air gaps of the RAMB in y-axes, and y1 is the deviation in y-axes. So the upper available air gap is calculated as (y0+y1cosα) and the lower one is calculated as (y0-y1cosα).i0 is the base current and the iy1 is the control current. i0+iy1 and i0-iy1 circulate the upper and the lower coils in the y-axis, respectively. f1 and f2 are the total attractive electromagnetic forces in the y-axis, which is given as follows:
(1)f1=μ0A0N24(i0+iy1)2(y0+y1)2cosα(2)f2=μ0A0N24(i0+iy1)2(y0-y1)2cosα,
where µ0 is the permeability μ0=4π×10-7 H/m, A0 is the effective cross-sectional area of one electromagnet, and N is the number of coils around the core.Additionally, the total nonlinear attractive electromagnetic forces for they-axis can be modeled as follows:
(3)fy=f1-f2=u0A0N24[(i0+iy1y0+y1)2-(i0-iy1y0-y1)2]cosα.
Moreover, by taking the Taylor’s expansions of (1) with respect to its nominal operating position (y1=0,iy1=0), the nonlinear electromagnetic forces can be represented by the following simplified linearized electromagnetic force models:
(4)fy1≅kii1+kyy1,
where ky and ki are the displacement and current stiffness parameters of the RAMB, respectively, and can be obtained from the y-axis as follows:
(5)ky≡∂fy1(y1,iy1)∂y1|y1=0,iy1=0=-4μ0A0N2cos2αi02y03ki≡∂fy1(y1,iy1)∂iy1|y1=0,iy1=0=4μ0A0N2cosαi0y02.
It is noted that since the coils in the x-axis and the y-axis are circulated by the same bias currents i0 and the nominal air gaps in the x-axis and the y-axis are also the same, that is the displacement and current stiffness parameters kx1 and ki1 obtained from the y-axis are the same as the ones obtained from the x-axis (where, the bearing has been assumed to be isotropic).


#### 2.2.2. The Controller of the AMB

The development of PID control has been around for 90 years and is still popular for industries and academies nowadays [20]. Due to the complexity of the rotor vibration problem, most engineers have turned to primitive control technologies such as the proportional-integral-derivative (PID) controller, which presently accounts for over 95% of all industrial control applications. This is not because practicing engineers are unaware of recently developed control methods, but because they find the advanced controllers difficult to tune and that it requires years of training. The overwhelming advantage in selecting PID over more advanced controllers is its ease of use, with only three tuning parameters and applicability to a vast range of plant models. It is very important to reduce the controller complexity whenever possible and thus the overall complexity of the closed loop system.The block diagram of the total control system of the AMB is shown in Figure3, which consists of electromagnet, rotor, displacement sensor, PID controller, D/A digital to analog converter, and A/D analog to digital converter.Figure 3
Block diagram of the control system.The transfer function of the PID controller can be written as(6)G(s)=U(s)E(s)=Kp+Kis+Kds1+Tds,
where Kp,Ki,Kd, and Td are the proportional parameter, the integral parameter, the derivative parameter, and the time constant of the PID controller, respectively.Combined with the gains of the power amplifier and the sensor, the characteristic equation of the system can be obtained as follows:(7)b5s5+b4s4+b3s3+b2s2+b1s1=0,
where
(8)b1=AaAsKpki-kxb2=ki(AaAsKd+AaAsKpTd)-kx(Ta+Td)b3=m-kxTdTab4=m(Ta+Td)b5=mTaTd,
where kx and ki are the displacement and current stiffness parameters of the RAMB, respectively; Aa and Ta are the amplifiable parameter and the time constant of the power amplifier, respectively; As and Ts are the amplifiable parameter and the time constant of the sensor, respectively.The parametersKp,Kd, and Td, of the controller can be determined by the Routh stability criterion, and the results are shown as follows (the detailed derivation is not listed in this section):
(9)0<TaTd<mkxx(10)(Kd+KaTd)<m(Ta+Td)AaAsTaTdkix(11)Kd>(s2-s22-4s1s32AaAss1kix+kxx(Ta+Td)AaAskix-KpTd)(12)Kd<(s2+s22-4s1s32AaAss1kix+kxx(Ta+Td)AaAskix-KpTd)(13)Kp≥(kixAaAskxx+(m-kxxTaTd)4mkixAaAsTaTd)(14)Kp>kxxAaAskix.
By the above determination, the derivative time parameter Td can be determined firstly, and then the appropriate value can be selected; similarly, the ranges of the proportional parameter Kp and the differential parameter Kd can also be gained. According to the effect of the actual control, the integral parameter Ki can be added to eliminate the lag error.The parameters of the bearing structure are set asµ0=4π×10-7 (H/m), and the initial bias current is selected as I0=0.15 A. Because of the air gap of the magnetic bearing S0=3×10-4 m, the scope of the actual changes of the rotor is in the range of 0 to 0.6 mm, the change of the A/D is in the range of 0 to 6,000 unit, and the equivalent gain As of the sensor can be obtained as 1×107 unit/m.The range of the differential time parameterKp determined by (9) is
(15)0<Td<0.8323.
With Td=0.0001, the range of the proportional parameter determined by (13), (14), and (9) can be carried out as
(16)0.154<Kp<4.578.
According to the optimal linear results, the proportional parameter Kp should be selected as Kp=1. And the range of the differential parameter Kd determined by (11) and (12) can be given as
(17)5.009×10-5<Kd<1.921×10-1.
Therefore, combined with (9), the controller’s parameters are finally selected as follows:
(18)Kp=1,Kd=2×10-3,Td=1×10-4.


### 2.3. Rotor System Modeling

In this study, the rotor is assumed to be a rigid and symmetric body. It is assumed that all magnets have identical structure. For simplicity, we neglect the magnetic flux leakage, the fringe magnetic flux, the eddy-current loss, the saturation and hysteresis of the magnetic core material, and the coupling effects between the electromagnets. The relationship between the center of mass(O) of the rotor and the five-DOF AMB is shown in Figure 4.Figure 4
Geometry relationships of rotor and AMB systems.o
′
x
′
y
′
z
′ is the body fixed coordinate system with the z′-axes overlapping the geometric center axis of the electronic bearing stator;Oxyz is space fixed coordinate system with the z-axes overlapping the geometric center axis of the rotor; O is the center of mass. fx1 and fy1 are the upper RAMB (RAMB1 in Figure 4) electromagnetic force of the rotor corresponding to the x- and y-axes; fx2 and fy2 are the lower RAMB (RAMB2 in Figure 4) electromagnetic force of the rotor corresponding to the x- and y-axes; fz is the TAMB electromagnetic force of the rotor corresponding to the z-axes; fb is the centrifugal force due to the static imbalance; lz1 and lz2 represent the distances from the O to the upper RAMB (RAMB1 in Figure 4) and lower RAMB (RAMB2 in the Figure 4), respectively; lc1 and lc2 represent the distances from the O to the upper position sensor and lower position sensor, respectively; l3 represents the distances from the O to the eccentric mass point; x1 and y1 denote the displacement output of the upper position sensor; x2 and y2 denote the displacement output of the lower position sensor; α,β, and φ denote the pitch, yaw, and spin angles displacements around the x-, y-, and z-axes of the rotor; θ is the initial phase of the centrifugal force due to the static imbalance; the rotational speed of the rotor can be denoted by ω.x
,
y
,
z is the displacement of the rotor mass center referring to the absolute coordinate system. The relationship between the angular and translational motions of the mass center and rotation, stress, and bearing is given as follows:
(19)x1=x+lc1αy1=y-lc1βx2=x-lc2αy2=y+lc2βz=zωt=φ(20)fx=fx1+fx2fy=fy1+fy2fz=fzTα=fx1lz1-fx2lz2Tβ=fy1lz2-fy1lz1Tφ=Td.Equations (19) are written in a matrix form for simple description as follows:
(21){x1y1x2y2zωt}=[1lc10000001-lc1001-lc20000001lc200000010000001]{xαyβzφ}.Equation (21) is written as Y=AX, Y is the sensor output signal matrix, and X is displacement of the rotor mass center.Equations (20) are also written in a matrix form for the simple description as follows:
(22){fxTαfyTβfzTφ}=[101000lz10-lz20000101000-lz10lz200000010000001]{fx1fy1fx2fy2fzTd}.Equation (22) is directly translated as {G}=[B]{F}, where {G} denotes the forces vector of the mass center of the rotor as
(23){G}=[fxTafyTβfzTφ]T.{
F
} denotes the electromagnetic forces vector produced by the AMBs as
(24){F}=[fx1fy1fx2fy2fzTd]T.It is noted that when the rotor is regulated perfectly(x1=y1=x2=y2=z=0), the rotational speed of the rotor can be denoted by ω. As shown in Figure 3, the control characteristics of the five-DOF AMB system are highly nonlinear and time varying because of the system parameter variations, external disturbances, and inherent nonlinearity such as the coupling effects among five axes and gyroscopic effects of rotation. Therefore, the dynamic states that is, the rotor positions x1,x2,y1,y2, and z of the five-DOF AMB system, are coupled and affected by more than one force. But, observing from (22), the couple of the radial and axial motion of the system studied in this paper is so small that the dynamic model of the five-DOF AMB system is decoupled as five independent subsystems including four radial subsystems according to x1,x2,y1, and y2 axes and one axial subsystem according to the z-axis to achieve the decentralized control. So this study emphasizes the dynamics behavior analysis of the four radial subsystems.When the rotor is suspended steadily, the relative coordinate system of the flywheel rotor coincides with the absolute coordinate system of the AMB’s stator. The potential energy of the stator in the originO is zero; m is the mass of the rotor; g is the gravity constant; z are the coordinates of the O. The potential energy of the rotor is given as U=mgz.The total kinetic energy of the systemT consists of the translational kinetic energy Tp, rotational kinetic energy Tz, and kinetic energy of eccentric mass Tb.The expression for the translational kinetic energyTp is
(25)Tp=m2[(dxdt)2+(dydt)2]=m2(x˙2+y˙2).
The expression for the rotational kinetic energy Tz is
(26)Tz=12Jd(α˙2+β˙2)+12Jpω(αβ˙-α˙β),
where Jd is transverse mass moments of inertia of the rotor, and Jp is the polar moment of inertia of the rotor. Due to the actual machining error and the assembly error, the rotor has the eccentric mass. If the rotor has a residual static unbalance following the ISO quality grade G6.3, it will produce the kinetic energy when the rotor is spinning. As shown in Figure 5, mb is the eccentric mass and rb is the eccentricity, so the expression for the kinetic energy of eccentric mass Tb is
(27)Tb=12mb(b˙x2+b˙y2+ω2rb2)+mb(b˙xωrbcosωt-b˙yωrbsinωt).
There is the displacement from the geometric center of the rotor to the one of the stator. bx, by represent the displacement corresponding to the x- and y-axes, respectively, as
(28)bx=x+l3αby=y-l3β.
Equation (28) is substituted to (27) as
(29)Tb=12mb[(x˙+l3α˙)2+(y˙-l3β˙)2+ω2rb2]+mb(y˙-l3β˙)ωrbcosωt-mb(x˙+l3α˙)ωrbsinωt.
So the total kinetic energy of the rotor is
(30)T=Tp+Tz+Tb=12m(x˙2+y˙2)+12Jd(α˙2+β˙2)+12Jpω(αβ˙-α˙β)+12mb[(x˙+l3α˙)2+(y˙-l3β˙)2+ω2rb2]+mb(y˙-l3β˙)ωrbcosωt-mb(x˙+l3α˙)ωrbsinωt.
According to the above analysis, the axial motion and the radial one are approximately orthogonal, and the influence of the radial bearing on the rotor only is analyzed as follows.Figure 5
Flywheel rotor with the eccentric mass.q
=
[
x
a
y
β
]
T denotes the displacement vector in the mass center of the rotor corresponding to the x- and y-axes.q
z
=
[
x
z
1
y
z
1
x
z
2
y
z
2
]
T denotes the displacement vector in the position of the RAMB corresponding to the x- and y-axes.By means of Lagrange equations, which state(31)ddt(∂T∂q˙i)-∂T∂qi+∂U∂qi=Qi,
where T is the total kinetic energy, U is the total potential energy, and Qi is the generalized force for the coordinate qi, a system of four equations describing the dynamics of the model is obtained, which can be written as
(32)(m+mb)x¨+mbl3α¨=fx1+fx2+ω2rbsinωtmbl3x¨+(Jd+mbl32)α¨-Jpβ˙ω=fx1lz1-fx2lz2+ω2rbl3sinωt(m+mb)y¨-mbl3β¨=fy1+fy2+ω2rbcosωt-mbl3y¨+(Jd+mbl32)β¨+Jpα˙ω=fy2lz2-fy1lz1+ω2rblbsinωt,
where fx1 and fy1 are the upper RAMB electromagnetic force of the rotor corresponding to the x- and y-axes; fx2 and fy2 are the lower RAMB electromagnetic force of the rotor corresponding to the x- and y-axes. The function of the electromagnetic force is given as follows:
(33)fx1=ki1i1+kx1x1fy1=ki1i1+kx1y1fx2=ki2i2+kx2x2fy2=ki2i2+kx2y2,
where ki1,ki2 are the up and down current stiffness of the active magnetic bearing, respectively. kx1,kx2 are the up and down displacement stiffness of the active magnetic bearing, respectively. The equation is translated as matrix operations for observation as follows
(34)[M]{q¨}+ω[G]{q˙}=[Bz]{Fz}+{Fb},
where [M] is the mass matrix. Because of the small value, the unbalance mass is neglected to make the simple calculation
(35)M=[m+mbmblb00mblbJd+mblb20000m+mb-mblb00-mblbJd+mblb2]≈[m0000Jd0000m0000Jd],
where [G] is the inertia matrix, and {Bz} is the length matrix as
(36)G=[0000000-Jp00000Jp00],Bz=[1010lz10-lz2001010-lz10lz2],
where {Fz} is unbalanced force vector and {Fb} is the force vector in bearing place as
(37)FZ=(fx1fy1fx2fy2)=(ki1i1+kx1x1ki1i1+kx1y1ki2i2+kx2x2ki2i2+kx2y2),Fb=(ω2rbsinωtω2rbl3sinωtω2rbcosωt-ω2rbl3cosωt).
The transfer function block diagram of the PID controller is shown in Figure 6.Figure 6
Transfer function block diagram of the PID controller.In Figure6, the current output depends on the input of the referenced displacement, which states
(38)I(s)=Kp(1+Tds)[-Y(s)+Yb(s)],
where i=Kpyz+KpTdy˙z-Kpy-KpTdy˙,y=[x1y1x2y2]T(39)BzFz=BzkiKp[(-Czq+yz)+Td(-Czq˙+y˙z)]+Bzkδqb,
where Cz=(1lc100001-lc11-lc100001lc1),qb=BzTq.The differential equation of the system can be given as(40)Mq¨+(ωG+D)q¨+Kq˙+Tq=F˙b+F1,
where D=kiKpTdBzCz,K=kiKpBzCz-kδBzBzT,T=KpkiBzCz/Ti, and F1=kiKpBzy˙z+kiKpTdBzy˙z.Whenyz=0, that means that the rotor is running steadily, and the differential equation can be written as
(41)Mq¨+(ωG+D)q¨+Kq˙+Tq=F˙b,
where Ti is the differential coefficient of the PID controller and F˙b is the force as
(42)F˙b=(rbω3cosωtrblbω3cosωt-rbω3sinωtrblbω3sinωt).
Then, in order to express this system in the state space form,
(43)x˙=Axx=[qq˙q¨]T,A=[0I000I-M-1Z-M-1K-M-1(ωG+D)]
are defined.From the above analysis,D, K, and T are directly involved with the controller’s parameters. The model will be variant according to the variation of the controller’s parameters. The dynamics characteristic of the system is given as follows.


## 2.1. System Structure

The basic layout of a flywheel energy storage system is depicted in Figure1, that is, the 600 Wh prototype system, a rigid vertical rotor shaft with a rigid motor and flywheel connected to shaft. Flywheel energy storage system is a complex construction where energy is stored mechanically and transferred to and from the flywheel by an integrated motor/generator.Figure 1
Basic layout of a flywheel energy storage system.Two radial active electromagnetic bearings (RAMBs) including upper and lower RAMBs (RAMB1 and RAMB2 in Figure1) and one thrust magnetic bearing (TAMB) are fixed on the platform to suspend and regulate the rotor in the radial and axial DOF, respectively. The positions of the rotor in five axes are defined as the displacements deviated from the nominal air gaps. Therefore, two pairs of the perpendicular eddy-current position sensors are installed closely to respective RAMBs to measure the respective radial positions in x- and y-axes which are denoted as x1 and y1 of the upper RAMB and x2 and y2 of the lower RAMB. Moreover, one eddy-current position sensor is installed closely to TAMB to measure the axial position in the z-axis which is denoted as z. After the rotor positions in five axes are all measured and sent to the control core through the position signals line, the required electromagnet currents are generated from the drive system and they circulate the coils through the power line. Therefore, the rotor can be regulated and stabilized in the centers of the apertures of the two RAMBs and the thrust disc can be centered in the middle of the air gap of the TAMB, respectively.One 70 kg flywheel is mounted on the rotor whose position can be adjusted to modify the rotational characteristics. Furthermore, one motor/generator is equipped between the lower RAMB and flywheel and is responsible for the torque generation or electricity generation based on the rotational speed requirement. However, in the rotor system, the rolling-element auxiliary bearings (auxiliary bearing 1 and auxiliary bearing 2 in Figure1) are necessary to protect AMB stators and stationary components along shaft in the even t of AMB failure or high transient loads. Under normal operation of AMBs, the rotor maintains a positive clearance with auxiliary bearings, which is less than a clearance with AMBs.


## 2.2. AMB Modeling

### 2.2.1. Electromagnetic Force Model

Considering the diameter of the shaft and the loading capacity, the RAMB with 8-pole legs was designed in the 600 Wh prototype system. The structural configuration of the RAMB is shown in Figure2(a), including electromagnetic coils. And the basic magnetic bearing control loop is shown in Figure 2(b).Structural configurations and the basic magnetic bearing control loop of the RAMB.
(a)
Structural configuration of the RAMB with 8-pole legs
(b)
The basic magnetic bearing control loop and its elementsTaking the two pairs of poles in they-axes of one RAMB, for example, the electromagnetic force model will be established [19]. As shown in Figure 2(b), the basic magnetic bearing control loop includes differential driving, position sensor, controller, and power amplifier. The included differential driving is adopted in this study to obtain maximum range of the force dynamic and good linearity of the control dynamic. α is the half of the angle between the two-pole legs, y0 is the nominal air gaps of the RAMB in y-axes, and y1 is the deviation in y-axes. So the upper available air gap is calculated as (y0+y1cosα) and the lower one is calculated as (y0-y1cosα).i0 is the base current and the iy1 is the control current. i0+iy1 and i0-iy1 circulate the upper and the lower coils in the y-axis, respectively. f1 and f2 are the total attractive electromagnetic forces in the y-axis, which is given as follows:
(1)f1=μ0A0N24(i0+iy1)2(y0+y1)2cosα(2)f2=μ0A0N24(i0+iy1)2(y0-y1)2cosα,
where µ0 is the permeability μ0=4π×10-7 H/m, A0 is the effective cross-sectional area of one electromagnet, and N is the number of coils around the core.Additionally, the total nonlinear attractive electromagnetic forces for they-axis can be modeled as follows:
(3)fy=f1-f2=u0A0N24[(i0+iy1y0+y1)2-(i0-iy1y0-y1)2]cosα.
Moreover, by taking the Taylor’s expansions of (1) with respect to its nominal operating position (y1=0,iy1=0), the nonlinear electromagnetic forces can be represented by the following simplified linearized electromagnetic force models:
(4)fy1≅kii1+kyy1,
where ky and ki are the displacement and current stiffness parameters of the RAMB, respectively, and can be obtained from the y-axis as follows:
(5)ky≡∂fy1(y1,iy1)∂y1|y1=0,iy1=0=-4μ0A0N2cos2αi02y03ki≡∂fy1(y1,iy1)∂iy1|y1=0,iy1=0=4μ0A0N2cosαi0y02.
It is noted that since the coils in the x-axis and the y-axis are circulated by the same bias currents i0 and the nominal air gaps in the x-axis and the y-axis are also the same, that is the displacement and current stiffness parameters kx1 and ki1 obtained from the y-axis are the same as the ones obtained from the x-axis (where, the bearing has been assumed to be isotropic).


### 2.2.2. The Controller of the AMB

The development of PID control has been around for 90 years and is still popular for industries and academies nowadays [20]. Due to the complexity of the rotor vibration problem, most engineers have turned to primitive control technologies such as the proportional-integral-derivative (PID) controller, which presently accounts for over 95% of all industrial control applications. This is not because practicing engineers are unaware of recently developed control methods, but because they find the advanced controllers difficult to tune and that it requires years of training. The overwhelming advantage in selecting PID over more advanced controllers is its ease of use, with only three tuning parameters and applicability to a vast range of plant models. It is very important to reduce the controller complexity whenever possible and thus the overall complexity of the closed loop system.The block diagram of the total control system of the AMB is shown in Figure3, which consists of electromagnet, rotor, displacement sensor, PID controller, D/A digital to analog converter, and A/D analog to digital converter.Figure 3
Block diagram of the control system.The transfer function of the PID controller can be written as(6)G(s)=U(s)E(s)=Kp+Kis+Kds1+Tds,
where Kp,Ki,Kd, and Td are the proportional parameter, the integral parameter, the derivative parameter, and the time constant of the PID controller, respectively.Combined with the gains of the power amplifier and the sensor, the characteristic equation of the system can be obtained as follows:(7)b5s5+b4s4+b3s3+b2s2+b1s1=0,
where
(8)b1=AaAsKpki-kxb2=ki(AaAsKd+AaAsKpTd)-kx(Ta+Td)b3=m-kxTdTab4=m(Ta+Td)b5=mTaTd,
where kx and ki are the displacement and current stiffness parameters of the RAMB, respectively; Aa and Ta are the amplifiable parameter and the time constant of the power amplifier, respectively; As and Ts are the amplifiable parameter and the time constant of the sensor, respectively.The parametersKp,Kd, and Td, of the controller can be determined by the Routh stability criterion, and the results are shown as follows (the detailed derivation is not listed in this section):
(9)0<TaTd<mkxx(10)(Kd+KaTd)<m(Ta+Td)AaAsTaTdkix(11)Kd>(s2-s22-4s1s32AaAss1kix+kxx(Ta+Td)AaAskix-KpTd)(12)Kd<(s2+s22-4s1s32AaAss1kix+kxx(Ta+Td)AaAskix-KpTd)(13)Kp≥(kixAaAskxx+(m-kxxTaTd)4mkixAaAsTaTd)(14)Kp>kxxAaAskix.
By the above determination, the derivative time parameter Td can be determined firstly, and then the appropriate value can be selected; similarly, the ranges of the proportional parameter Kp and the differential parameter Kd can also be gained. According to the effect of the actual control, the integral parameter Ki can be added to eliminate the lag error.The parameters of the bearing structure are set asµ0=4π×10-7 (H/m), and the initial bias current is selected as I0=0.15 A. Because of the air gap of the magnetic bearing S0=3×10-4 m, the scope of the actual changes of the rotor is in the range of 0 to 0.6 mm, the change of the A/D is in the range of 0 to 6,000 unit, and the equivalent gain As of the sensor can be obtained as 1×107 unit/m.The range of the differential time parameterKp determined by (9) is
(15)0<Td<0.8323.
With Td=0.0001, the range of the proportional parameter determined by (13), (14), and (9) can be carried out as
(16)0.154<Kp<4.578.
According to the optimal linear results, the proportional parameter Kp should be selected as Kp=1. And the range of the differential parameter Kd determined by (11) and (12) can be given as
(17)5.009×10-5<Kd<1.921×10-1.
Therefore, combined with (9), the controller’s parameters are finally selected as follows:
(18)Kp=1,Kd=2×10-3,Td=1×10-4.


## 2.2.1. Electromagnetic Force Model

Considering the diameter of the shaft and the loading capacity, the RAMB with 8-pole legs was designed in the 600 Wh prototype system. The structural configuration of the RAMB is shown in Figure2(a), including electromagnetic coils. And the basic magnetic bearing control loop is shown in Figure 2(b).Structural configurations and the basic magnetic bearing control loop of the RAMB.
(a)
Structural configuration of the RAMB with 8-pole legs
(b)
The basic magnetic bearing control loop and its elementsTaking the two pairs of poles in they-axes of one RAMB, for example, the electromagnetic force model will be established [19]. As shown in Figure 2(b), the basic magnetic bearing control loop includes differential driving, position sensor, controller, and power amplifier. The included differential driving is adopted in this study to obtain maximum range of the force dynamic and good linearity of the control dynamic. α is the half of the angle between the two-pole legs, y0 is the nominal air gaps of the RAMB in y-axes, and y1 is the deviation in y-axes. So the upper available air gap is calculated as (y0+y1cosα) and the lower one is calculated as (y0-y1cosα).i0 is the base current and the iy1 is the control current. i0+iy1 and i0-iy1 circulate the upper and the lower coils in the y-axis, respectively. f1 and f2 are the total attractive electromagnetic forces in the y-axis, which is given as follows:
(1)f1=μ0A0N24(i0+iy1)2(y0+y1)2cosα(2)f2=μ0A0N24(i0+iy1)2(y0-y1)2cosα,
where µ0 is the permeability μ0=4π×10-7 H/m, A0 is the effective cross-sectional area of one electromagnet, and N is the number of coils around the core.Additionally, the total nonlinear attractive electromagnetic forces for they-axis can be modeled as follows:
(3)fy=f1-f2=u0A0N24[(i0+iy1y0+y1)2-(i0-iy1y0-y1)2]cosα.
Moreover, by taking the Taylor’s expansions of (1) with respect to its nominal operating position (y1=0,iy1=0), the nonlinear electromagnetic forces can be represented by the following simplified linearized electromagnetic force models:
(4)fy1≅kii1+kyy1,
where ky and ki are the displacement and current stiffness parameters of the RAMB, respectively, and can be obtained from the y-axis as follows:
(5)ky≡∂fy1(y1,iy1)∂y1|y1=0,iy1=0=-4μ0A0N2cos2αi02y03ki≡∂fy1(y1,iy1)∂iy1|y1=0,iy1=0=4μ0A0N2cosαi0y02.
It is noted that since the coils in the x-axis and the y-axis are circulated by the same bias currents i0 and the nominal air gaps in the x-axis and the y-axis are also the same, that is the displacement and current stiffness parameters kx1 and ki1 obtained from the y-axis are the same as the ones obtained from the x-axis (where, the bearing has been assumed to be isotropic).


## 2.2.2. The Controller of the AMB

The development of PID control has been around for 90 years and is still popular for industries and academies nowadays [20]. Due to the complexity of the rotor vibration problem, most engineers have turned to primitive control technologies such as the proportional-integral-derivative (PID) controller, which presently accounts for over 95% of all industrial control applications. This is not because practicing engineers are unaware of recently developed control methods, but because they find the advanced controllers difficult to tune and that it requires years of training. The overwhelming advantage in selecting PID over more advanced controllers is its ease of use, with only three tuning parameters and applicability to a vast range of plant models. It is very important to reduce the controller complexity whenever possible and thus the overall complexity of the closed loop system.The block diagram of the total control system of the AMB is shown in Figure3, which consists of electromagnet, rotor, displacement sensor, PID controller, D/A digital to analog converter, and A/D analog to digital converter.Figure 3
Block diagram of the control system.The transfer function of the PID controller can be written as(6)G(s)=U(s)E(s)=Kp+Kis+Kds1+Tds,
where Kp,Ki,Kd, and Td are the proportional parameter, the integral parameter, the derivative parameter, and the time constant of the PID controller, respectively.Combined with the gains of the power amplifier and the sensor, the characteristic equation of the system can be obtained as follows:(7)b5s5+b4s4+b3s3+b2s2+b1s1=0,
where
(8)b1=AaAsKpki-kxb2=ki(AaAsKd+AaAsKpTd)-kx(Ta+Td)b3=m-kxTdTab4=m(Ta+Td)b5=mTaTd,
where kx and ki are the displacement and current stiffness parameters of the RAMB, respectively; Aa and Ta are the amplifiable parameter and the time constant of the power amplifier, respectively; As and Ts are the amplifiable parameter and the time constant of the sensor, respectively.The parametersKp,Kd, and Td, of the controller can be determined by the Routh stability criterion, and the results are shown as follows (the detailed derivation is not listed in this section):
(9)0<TaTd<mkxx(10)(Kd+KaTd)<m(Ta+Td)AaAsTaTdkix(11)Kd>(s2-s22-4s1s32AaAss1kix+kxx(Ta+Td)AaAskix-KpTd)(12)Kd<(s2+s22-4s1s32AaAss1kix+kxx(Ta+Td)AaAskix-KpTd)(13)Kp≥(kixAaAskxx+(m-kxxTaTd)4mkixAaAsTaTd)(14)Kp>kxxAaAskix.
By the above determination, the derivative time parameter Td can be determined firstly, and then the appropriate value can be selected; similarly, the ranges of the proportional parameter Kp and the differential parameter Kd can also be gained. According to the effect of the actual control, the integral parameter Ki can be added to eliminate the lag error.The parameters of the bearing structure are set asµ0=4π×10-7 (H/m), and the initial bias current is selected as I0=0.15 A. Because of the air gap of the magnetic bearing S0=3×10-4 m, the scope of the actual changes of the rotor is in the range of 0 to 0.6 mm, the change of the A/D is in the range of 0 to 6,000 unit, and the equivalent gain As of the sensor can be obtained as 1×107 unit/m.The range of the differential time parameterKp determined by (9) is
(15)0<Td<0.8323.
With Td=0.0001, the range of the proportional parameter determined by (13), (14), and (9) can be carried out as
(16)0.154<Kp<4.578.
According to the optimal linear results, the proportional parameter Kp should be selected as Kp=1. And the range of the differential parameter Kd determined by (11) and (12) can be given as
(17)5.009×10-5<Kd<1.921×10-1.
Therefore, combined with (9), the controller’s parameters are finally selected as follows:
(18)Kp=1,Kd=2×10-3,Td=1×10-4.


## 2.3. Rotor System Modeling

In this study, the rotor is assumed to be a rigid and symmetric body. It is assumed that all magnets have identical structure. For simplicity, we neglect the magnetic flux leakage, the fringe magnetic flux, the eddy-current loss, the saturation and hysteresis of the magnetic core material, and the coupling effects between the electromagnets. The relationship between the center of mass(O) of the rotor and the five-DOF AMB is shown in Figure 4.Figure 4
Geometry relationships of rotor and AMB systems.o
′
x
′
y
′
z
′ is the body fixed coordinate system with the z′-axes overlapping the geometric center axis of the electronic bearing stator;Oxyz is space fixed coordinate system with the z-axes overlapping the geometric center axis of the rotor; O is the center of mass. fx1 and fy1 are the upper RAMB (RAMB1 in Figure 4) electromagnetic force of the rotor corresponding to the x- and y-axes; fx2 and fy2 are the lower RAMB (RAMB2 in Figure 4) electromagnetic force of the rotor corresponding to the x- and y-axes; fz is the TAMB electromagnetic force of the rotor corresponding to the z-axes; fb is the centrifugal force due to the static imbalance; lz1 and lz2 represent the distances from the O to the upper RAMB (RAMB1 in Figure 4) and lower RAMB (RAMB2 in the Figure 4), respectively; lc1 and lc2 represent the distances from the O to the upper position sensor and lower position sensor, respectively; l3 represents the distances from the O to the eccentric mass point; x1 and y1 denote the displacement output of the upper position sensor; x2 and y2 denote the displacement output of the lower position sensor; α,β, and φ denote the pitch, yaw, and spin angles displacements around the x-, y-, and z-axes of the rotor; θ is the initial phase of the centrifugal force due to the static imbalance; the rotational speed of the rotor can be denoted by ω.x
,
y
,
z is the displacement of the rotor mass center referring to the absolute coordinate system. The relationship between the angular and translational motions of the mass center and rotation, stress, and bearing is given as follows:
(19)x1=x+lc1αy1=y-lc1βx2=x-lc2αy2=y+lc2βz=zωt=φ(20)fx=fx1+fx2fy=fy1+fy2fz=fzTα=fx1lz1-fx2lz2Tβ=fy1lz2-fy1lz1Tφ=Td.Equations (19) are written in a matrix form for simple description as follows:
(21){x1y1x2y2zωt}=[1lc10000001-lc1001-lc20000001lc200000010000001]{xαyβzφ}.Equation (21) is written as Y=AX, Y is the sensor output signal matrix, and X is displacement of the rotor mass center.Equations (20) are also written in a matrix form for the simple description as follows:
(22){fxTαfyTβfzTφ}=[101000lz10-lz20000101000-lz10lz200000010000001]{fx1fy1fx2fy2fzTd}.Equation (22) is directly translated as {G}=[B]{F}, where {G} denotes the forces vector of the mass center of the rotor as
(23){G}=[fxTafyTβfzTφ]T.{
F
} denotes the electromagnetic forces vector produced by the AMBs as
(24){F}=[fx1fy1fx2fy2fzTd]T.It is noted that when the rotor is regulated perfectly(x1=y1=x2=y2=z=0), the rotational speed of the rotor can be denoted by ω. As shown in Figure 3, the control characteristics of the five-DOF AMB system are highly nonlinear and time varying because of the system parameter variations, external disturbances, and inherent nonlinearity such as the coupling effects among five axes and gyroscopic effects of rotation. Therefore, the dynamic states that is, the rotor positions x1,x2,y1,y2, and z of the five-DOF AMB system, are coupled and affected by more than one force. But, observing from (22), the couple of the radial and axial motion of the system studied in this paper is so small that the dynamic model of the five-DOF AMB system is decoupled as five independent subsystems including four radial subsystems according to x1,x2,y1, and y2 axes and one axial subsystem according to the z-axis to achieve the decentralized control. So this study emphasizes the dynamics behavior analysis of the four radial subsystems.When the rotor is suspended steadily, the relative coordinate system of the flywheel rotor coincides with the absolute coordinate system of the AMB’s stator. The potential energy of the stator in the originO is zero; m is the mass of the rotor; g is the gravity constant; z are the coordinates of the O. The potential energy of the rotor is given as U=mgz.The total kinetic energy of the systemT consists of the translational kinetic energy Tp, rotational kinetic energy Tz, and kinetic energy of eccentric mass Tb.The expression for the translational kinetic energyTp is
(25)Tp=m2[(dxdt)2+(dydt)2]=m2(x˙2+y˙2).
The expression for the rotational kinetic energy Tz is
(26)Tz=12Jd(α˙2+β˙2)+12Jpω(αβ˙-α˙β),
where Jd is transverse mass moments of inertia of the rotor, and Jp is the polar moment of inertia of the rotor. Due to the actual machining error and the assembly error, the rotor has the eccentric mass. If the rotor has a residual static unbalance following the ISO quality grade G6.3, it will produce the kinetic energy when the rotor is spinning. As shown in Figure 5, mb is the eccentric mass and rb is the eccentricity, so the expression for the kinetic energy of eccentric mass Tb is
(27)Tb=12mb(b˙x2+b˙y2+ω2rb2)+mb(b˙xωrbcosωt-b˙yωrbsinωt).
There is the displacement from the geometric center of the rotor to the one of the stator. bx, by represent the displacement corresponding to the x- and y-axes, respectively, as
(28)bx=x+l3αby=y-l3β.
Equation (28) is substituted to (27) as
(29)Tb=12mb[(x˙+l3α˙)2+(y˙-l3β˙)2+ω2rb2]+mb(y˙-l3β˙)ωrbcosωt-mb(x˙+l3α˙)ωrbsinωt.
So the total kinetic energy of the rotor is
(30)T=Tp+Tz+Tb=12m(x˙2+y˙2)+12Jd(α˙2+β˙2)+12Jpω(αβ˙-α˙β)+12mb[(x˙+l3α˙)2+(y˙-l3β˙)2+ω2rb2]+mb(y˙-l3β˙)ωrbcosωt-mb(x˙+l3α˙)ωrbsinωt.
According to the above analysis, the axial motion and the radial one are approximately orthogonal, and the influence of the radial bearing on the rotor only is analyzed as follows.Figure 5
Flywheel rotor with the eccentric mass.q
=
[
x
a
y
β
]
T denotes the displacement vector in the mass center of the rotor corresponding to the x- and y-axes.q
z
=
[
x
z
1
y
z
1
x
z
2
y
z
2
]
T denotes the displacement vector in the position of the RAMB corresponding to the x- and y-axes.By means of Lagrange equations, which state(31)ddt(∂T∂q˙i)-∂T∂qi+∂U∂qi=Qi,
where T is the total kinetic energy, U is the total potential energy, and Qi is the generalized force for the coordinate qi, a system of four equations describing the dynamics of the model is obtained, which can be written as
(32)(m+mb)x¨+mbl3α¨=fx1+fx2+ω2rbsinωtmbl3x¨+(Jd+mbl32)α¨-Jpβ˙ω=fx1lz1-fx2lz2+ω2rbl3sinωt(m+mb)y¨-mbl3β¨=fy1+fy2+ω2rbcosωt-mbl3y¨+(Jd+mbl32)β¨+Jpα˙ω=fy2lz2-fy1lz1+ω2rblbsinωt,
where fx1 and fy1 are the upper RAMB electromagnetic force of the rotor corresponding to the x- and y-axes; fx2 and fy2 are the lower RAMB electromagnetic force of the rotor corresponding to the x- and y-axes. The function of the electromagnetic force is given as follows:
(33)fx1=ki1i1+kx1x1fy1=ki1i1+kx1y1fx2=ki2i2+kx2x2fy2=ki2i2+kx2y2,
where ki1,ki2 are the up and down current stiffness of the active magnetic bearing, respectively. kx1,kx2 are the up and down displacement stiffness of the active magnetic bearing, respectively. The equation is translated as matrix operations for observation as follows
(34)[M]{q¨}+ω[G]{q˙}=[Bz]{Fz}+{Fb},
where [M] is the mass matrix. Because of the small value, the unbalance mass is neglected to make the simple calculation
(35)M=[m+mbmblb00mblbJd+mblb20000m+mb-mblb00-mblbJd+mblb2]≈[m0000Jd0000m0000Jd],
where [G] is the inertia matrix, and {Bz} is the length matrix as
(36)G=[0000000-Jp00000Jp00],Bz=[1010lz10-lz2001010-lz10lz2],
where {Fz} is unbalanced force vector and {Fb} is the force vector in bearing place as
(37)FZ=(fx1fy1fx2fy2)=(ki1i1+kx1x1ki1i1+kx1y1ki2i2+kx2x2ki2i2+kx2y2),Fb=(ω2rbsinωtω2rbl3sinωtω2rbcosωt-ω2rbl3cosωt).
The transfer function block diagram of the PID controller is shown in Figure 6.Figure 6
Transfer function block diagram of the PID controller.In Figure6, the current output depends on the input of the referenced displacement, which states
(38)I(s)=Kp(1+Tds)[-Y(s)+Yb(s)],
where i=Kpyz+KpTdy˙z-Kpy-KpTdy˙,y=[x1y1x2y2]T(39)BzFz=BzkiKp[(-Czq+yz)+Td(-Czq˙+y˙z)]+Bzkδqb,
where Cz=(1lc100001-lc11-lc100001lc1),qb=BzTq.The differential equation of the system can be given as(40)Mq¨+(ωG+D)q¨+Kq˙+Tq=F˙b+F1,
where D=kiKpTdBzCz,K=kiKpBzCz-kδBzBzT,T=KpkiBzCz/Ti, and F1=kiKpBzy˙z+kiKpTdBzy˙z.Whenyz=0, that means that the rotor is running steadily, and the differential equation can be written as
(41)Mq¨+(ωG+D)q¨+Kq˙+Tq=F˙b,
where Ti is the differential coefficient of the PID controller and F˙b is the force as
(42)F˙b=(rbω3cosωtrblbω3cosωt-rbω3sinωtrblbω3sinωt).
Then, in order to express this system in the state space form,
(43)x˙=Axx=[qq˙q¨]T,A=[0I000I-M-1Z-M-1K-M-1(ωG+D)]
are defined.From the above analysis,D, K, and T are directly involved with the controller’s parameters. The model will be variant according to the variation of the controller’s parameters. The dynamics characteristic of the system is given as follows.


## 3. Simulation Result Analysis

### 3.1. Systems Parameters

Suitable values for the parameters involved in the model of (41) were determined. Some of these parameters will be varied, with the goal of obtaining useful information for the optimization of designs. The parameters are shown in Table 1. However, if not otherwise stated, the following values will be the ones used for analysis.Table 1
Model parameters.
Parameter
Description
Value
Unit
m
Mass of rotor
70
kg
L
Length of rotor
0.505
m
L
z
1
Distance betweenO and upper RAMB
0.3047
m
L
z
2
Distance betweenO and lower RAMB
0.2463
m
L
c
1
Distance betweenO and upper position sensor
0.3112
m
L
c
2
Distance betweenO and lower position sensor
0.2518
m
J
d
Transverse mass moments of inertia of rotor about thex- and y-axes
1.6426
Kg
·
m
2
J
p
Polar mass moment of inertia of rotor about thez-axis
1.9235
Kg
·
m
2
r
b
Distance between CG and external disturbances
0.263
m
S
0
Nominal air gaps in thex- and y-axes of RAMB
0
.
3
mm
μ
0
Electromagnet parameter
4
π
×
10
-
7
H/m
i
0
Bias currents
0.15
A
z
0
Nominal air gap in thez-axis of TAMB
0
.
3
mm
ω
max
Design maximum speed
20000
rpm
ω
Rated speed
5000~16000
rpm


### 3.2. Simulation Result Analysis

A planer motion of a rotor is called a whirling motion or a whirl. And a circular whirl in the same direction as the shaft rotation is called a forward whirl, and that in the opposite direction is termed a backward whirl. Plot of these natural angular frequencies versus the rotational speed is called a natural angular frequency diagrams (or shortly natural frequency diagram) or ap=Ω diagram [21]. Positive and negative values of p correspond to forward and backward whirls, respectively. Diagram of natural frequencies versus the rotational speed is also used. Its natural frequency diagram is illustrated in Figure 6 where natural angular frequencies become horizontal straight lines in the diagram. Resonance mentioned later occurs at the angular frequency given by the cross-point of this straight line and the line p=Ω. A forced oscillation becomes large in the neighborhood of this resonance frequency. Sometimes, radii of these whirling motions are represented by relative sizes of circles in the diagram as shown in Figure 6. These diagrams are called Campbell diagram, which is used to analyze the transient dynamics characteristic of the rotor system.According to Table1, the differential coefficient of the AMB’s controller is set as 0.002, and the proportional coefficient is set as 1. Extracting the eigenvalues of A for the parameters given above with the running speed, ω, varying from 0 to 20,000 RPM, the Campbell diagram of the rotor is obtained.Examination of the results of the simulation using the model of the rotor allows us to argument that rotational modes of vibration have natural frequencies that are considerably larger than those of translational modes. The forward whirling modes that represent critical speeds, that is, the ones that intersect with the line describing the frequencies equal to the running speeds in Figure7, are all translational. The rotational modes do not intersect with that line but with their backward directions (negative slopes). These natural frequencies increase rapidly with the running speed and do not represent critical speeds, and thus stability problems are not expected to occur in these modes.Figure 7
Campbell diagram of the system.In Figure7, we can know that the critical speed of the rotor is about 2800 RPM, the backward whirling frequency is about 48 Hz, and the translational mode is kept as a constant value 42 Hz. So the stability problems are not expected to occur in working running speeds. It is needed to accelerate or decelerate quickly over 2800 RPM, so as to keep the system in stability state.In Figure8, the integral parameter Kd is kept at 0.002, and the proportional parameter Kp is set as 0.5, 1, 1.5, and 2, respectively. With the increase of the proportional parameter Kp, the forward whirling mode curve is away from the curve P=Ω, the cross-point between the backward whirling curve and the curve P=Ω gradually rises, and the critical speed increases as well. And the translational mode also increases, which distinctly showed that proportional parameter affects the stiffness of the rotor system; that is, the proportional parameter plays a very important role in the stiffness of the rotor system.Figure 8
Campbell diagram of the system with the change of the proportional parameter.In Figure9, the proportional parameter Kp is kept at 1, and the differential parameter Kd is set at 0.001, 0.002, 0.003, and 0.004, respectively. With the increase of the differential parameter Kd, the forward whirling mode curve is away from the curve P=Ω, the cross-point between the backward whirling curve and the curve P=Ω gradually declined, and the critical speed decreases as well.Figure 9
Campbell diagram of the system with the change of the differential parameter.The backward whirling modes change distinctly, 3200 RPM in theKd=0.001 and 2000 RPM in the Kd=0.004.The effect of differential parameters on the whirling mode is distinct in running speed0~6000RPM, but as the running speed rises, the effect of the differential parameter Kd becomes more and more smaller. And the translational mode also increases, which showed that the differential parameter also affects the stiffness of the rotor system, but the influence is not very distinctly compared with the influence of the proportional parameter.The following analysis shows how the dynamics of the system would change if the distance between the two RAMBs changed in the axial direction. The controller parameterKp is selected as 1 and Kd is set to 0.002; the Campbell diagram of the system with the change of the RAMB positions is shown in Figure 10. Curve 1 shows that the positions of RAMBs are set close to the center 50 mm, respectively. Curve 2 shows that the positions of RAMBs are not be changed. Curve 3 shows that the positions of RAMBs are set far away from the center of the system 50 mm, respectively. Compared to the controller’s parameters, the change of the RAMB positions does not change the system's translational mode but just changes the whirling mode of the system, which can verify that the change of the RAMB positions does not change the overall stiffness of the system. With the positions of RAMBs far away from the center of the rotor, the critical speed of the backward whirling mode increases from 2400 rpm to 3800 rpm. The result of the analysis, which indicates that the shorter distance between the locations of the two RAMBs bearing can effectively reduce the whirling frequency of the rotor, can be used to direct the design of the rotor system.Figure 10
Campbell diagram of the system with the change of the RAMB positions.The purpose of the flywheel storage machinery is to store energy as much as possible. The polar mass moment of inertia of rotor determines the energy storage of the flywheel storage machinery. So the effect of the whirling modes of the system is gained as follows. The polar mass moment and transverse mass moments of inertia of rotor depend on the geometry of the rotor. Figure10 is the Campbell diagram of the system with the inertia ratio of the polar moment and the transverse mass moment.In Figure11, the Jp/Jd of the curve 1 is 1/3, that of curve 2 is 1/2, that of curve 3 is 1.17 (the designed one of the rotor system), and that of curve 4 is 2. What we can see from this figure is that with the ratio of the moments of inertia increasing,Jp/Jd from almost 1/3 to about 2, and the forward rotational mode has increased; at the same time the backward rotational mode has become smaller. The main consequence of this occurrence is that this mode now intersects (the curve 1) with the line describing the frequencies equal to the running speeds and thus becomes a critical speed, with the potential of becoming an unstable mode for a certain running speed. Compared to the controller’s parameters, the change of the ratio of the moments of inertia also does not change the system's translational mode but just changes the whirling mode of the system, which can verify that the change of the RAMB positions does not change the overall stiffness of the system.Figure 11
Campbell diagram of the system with the inertia ratio of the polar moment and the transverse mass moment.With a constant magnetic force and a constant polar mass moment of inertia of rotor, the transverse mass moments of inertia of rotor can be changed just by the change of the distance between the local center to the overall center of the rotor to avoid the work speed in the strong whirling mode, which can greatly improve the stability of the system.


## 3.1. Systems Parameters

Suitable values for the parameters involved in the model of (41) were determined. Some of these parameters will be varied, with the goal of obtaining useful information for the optimization of designs. The parameters are shown in Table 1. However, if not otherwise stated, the following values will be the ones used for analysis.Table 1
Model parameters.
Parameter
Description
Value
Unit
m
Mass of rotor
70
kg
L
Length of rotor
0.505
m
L
z
1
Distance betweenO and upper RAMB
0.3047
m
L
z
2
Distance betweenO and lower RAMB
0.2463
m
L
c
1
Distance betweenO and upper position sensor
0.3112
m
L
c
2
Distance betweenO and lower position sensor
0.2518
m
J
d
Transverse mass moments of inertia of rotor about thex- and y-axes
1.6426
Kg
·
m
2
J
p
Polar mass moment of inertia of rotor about thez-axis
1.9235
Kg
·
m
2
r
b
Distance between CG and external disturbances
0.263
m
S
0
Nominal air gaps in thex- and y-axes of RAMB
0
.
3
mm
μ
0
Electromagnet parameter
4
π
×
10
-
7
H/m
i
0
Bias currents
0.15
A
z
0
Nominal air gap in thez-axis of TAMB
0
.
3
mm
ω
max
Design maximum speed
20000
rpm
ω
Rated speed
5000~16000
rpm


## 3.2. Simulation Result Analysis

A planer motion of a rotor is called a whirling motion or a whirl. And a circular whirl in the same direction as the shaft rotation is called a forward whirl, and that in the opposite direction is termed a backward whirl. Plot of these natural angular frequencies versus the rotational speed is called a natural angular frequency diagrams (or shortly natural frequency diagram) or ap=Ω diagram [21]. Positive and negative values of p correspond to forward and backward whirls, respectively. Diagram of natural frequencies versus the rotational speed is also used. Its natural frequency diagram is illustrated in Figure 6 where natural angular frequencies become horizontal straight lines in the diagram. Resonance mentioned later occurs at the angular frequency given by the cross-point of this straight line and the line p=Ω. A forced oscillation becomes large in the neighborhood of this resonance frequency. Sometimes, radii of these whirling motions are represented by relative sizes of circles in the diagram as shown in Figure 6. These diagrams are called Campbell diagram, which is used to analyze the transient dynamics characteristic of the rotor system.According to Table1, the differential coefficient of the AMB’s controller is set as 0.002, and the proportional coefficient is set as 1. Extracting the eigenvalues of A for the parameters given above with the running speed, ω, varying from 0 to 20,000 RPM, the Campbell diagram of the rotor is obtained.Examination of the results of the simulation using the model of the rotor allows us to argument that rotational modes of vibration have natural frequencies that are considerably larger than those of translational modes. The forward whirling modes that represent critical speeds, that is, the ones that intersect with the line describing the frequencies equal to the running speeds in Figure7, are all translational. The rotational modes do not intersect with that line but with their backward directions (negative slopes). These natural frequencies increase rapidly with the running speed and do not represent critical speeds, and thus stability problems are not expected to occur in these modes.Figure 7
Campbell diagram of the system.In Figure7, we can know that the critical speed of the rotor is about 2800 RPM, the backward whirling frequency is about 48 Hz, and the translational mode is kept as a constant value 42 Hz. So the stability problems are not expected to occur in working running speeds. It is needed to accelerate or decelerate quickly over 2800 RPM, so as to keep the system in stability state.In Figure8, the integral parameter Kd is kept at 0.002, and the proportional parameter Kp is set as 0.5, 1, 1.5, and 2, respectively. With the increase of the proportional parameter Kp, the forward whirling mode curve is away from the curve P=Ω, the cross-point between the backward whirling curve and the curve P=Ω gradually rises, and the critical speed increases as well. And the translational mode also increases, which distinctly showed that proportional parameter affects the stiffness of the rotor system; that is, the proportional parameter plays a very important role in the stiffness of the rotor system.Figure 8
Campbell diagram of the system with the change of the proportional parameter.In Figure9, the proportional parameter Kp is kept at 1, and the differential parameter Kd is set at 0.001, 0.002, 0.003, and 0.004, respectively. With the increase of the differential parameter Kd, the forward whirling mode curve is away from the curve P=Ω, the cross-point between the backward whirling curve and the curve P=Ω gradually declined, and the critical speed decreases as well.Figure 9
Campbell diagram of the system with the change of the differential parameter.The backward whirling modes change distinctly, 3200 RPM in theKd=0.001 and 2000 RPM in the Kd=0.004.The effect of differential parameters on the whirling mode is distinct in running speed0~6000RPM, but as the running speed rises, the effect of the differential parameter Kd becomes more and more smaller. And the translational mode also increases, which showed that the differential parameter also affects the stiffness of the rotor system, but the influence is not very distinctly compared with the influence of the proportional parameter.The following analysis shows how the dynamics of the system would change if the distance between the two RAMBs changed in the axial direction. The controller parameterKp is selected as 1 and Kd is set to 0.002; the Campbell diagram of the system with the change of the RAMB positions is shown in Figure 10. Curve 1 shows that the positions of RAMBs are set close to the center 50 mm, respectively. Curve 2 shows that the positions of RAMBs are not be changed. Curve 3 shows that the positions of RAMBs are set far away from the center of the system 50 mm, respectively. Compared to the controller’s parameters, the change of the RAMB positions does not change the system's translational mode but just changes the whirling mode of the system, which can verify that the change of the RAMB positions does not change the overall stiffness of the system. With the positions of RAMBs far away from the center of the rotor, the critical speed of the backward whirling mode increases from 2400 rpm to 3800 rpm. The result of the analysis, which indicates that the shorter distance between the locations of the two RAMBs bearing can effectively reduce the whirling frequency of the rotor, can be used to direct the design of the rotor system.Figure 10
Campbell diagram of the system with the change of the RAMB positions.The purpose of the flywheel storage machinery is to store energy as much as possible. The polar mass moment of inertia of rotor determines the energy storage of the flywheel storage machinery. So the effect of the whirling modes of the system is gained as follows. The polar mass moment and transverse mass moments of inertia of rotor depend on the geometry of the rotor. Figure10 is the Campbell diagram of the system with the inertia ratio of the polar moment and the transverse mass moment.In Figure11, the Jp/Jd of the curve 1 is 1/3, that of curve 2 is 1/2, that of curve 3 is 1.17 (the designed one of the rotor system), and that of curve 4 is 2. What we can see from this figure is that with the ratio of the moments of inertia increasing,Jp/Jd from almost 1/3 to about 2, and the forward rotational mode has increased; at the same time the backward rotational mode has become smaller. The main consequence of this occurrence is that this mode now intersects (the curve 1) with the line describing the frequencies equal to the running speeds and thus becomes a critical speed, with the potential of becoming an unstable mode for a certain running speed. Compared to the controller’s parameters, the change of the ratio of the moments of inertia also does not change the system's translational mode but just changes the whirling mode of the system, which can verify that the change of the RAMB positions does not change the overall stiffness of the system.Figure 11
Campbell diagram of the system with the inertia ratio of the polar moment and the transverse mass moment.With a constant magnetic force and a constant polar mass moment of inertia of rotor, the transverse mass moments of inertia of rotor can be changed just by the change of the distance between the local center to the overall center of the rotor to avoid the work speed in the strong whirling mode, which can greatly improve the stability of the system.


## 4. Conclusions

Variation studies were conducted to assess the influence of controller and the rotor geometry, with running speed on rotor dynamic stability. Each point on the lines (or surface) represents the threshold running speed above which the rotor becomes unstable for a certain parameter configuration. This means that there is a distinctive operating speed below which the system is always stable for a given parametric configuration.Considering the influence of controller, the dynamic model of the rigidity flywheel rotor supporting by AMBs was established to analyze the dynamics characteristic. The influence analysis of the controller’s parameters to the dynamics characteristic was obtained.(
1
) The rotational mode and critical speed will reduce as the proportional parameter increases. And the translational mode also increases, which distinctly showed that proportional parameter affects the stiffness of the rotor system; that is, the proportional parameter plays a very important role to the stiffness of the rotor system.The rotational mode will increase as the differential parameter increases. And the translational mode also increases, which showed that the differential parameter also affects the stiffness of the rotor system, but the influence is not very distinctly compared with the influence of the proportional parameter.(
2
) The change of the distance between the mass center of the rotor to the position of the RAMB will change the cross-point of this straight line and the line P=Ω but will not change the trend of the total whirling frequency curve. Compared to the controller’s parameters, the change of the RAMB positions does not change the system's translational mode but just changes the whirling mode of the system, which can verify that the change of the RAMB positions does not change the overall stiffness of the system.(
3
) The ratio Jp/Jd is 1.17 and is close to critical value 1. There is not any cross-point of the curve of the forward whirling mode and the line P=Ω. These natural frequencies increase rapidly with the running speed and do not represent critical speeds, and thus stability problems are not expected to occur in these modes. Compared to the controller’s parameters, the change of the ratio of the moments of inertia also does not change the system's translational mode but just changes the whirling mode of the system, which can verify that the change of the RAMB positions does not change the overall stiffness of the system.The result of the analysis can be used to set the support position of the rotor system, limit the ratio of transverse moment of inertia and the polar moment of inertia, and direct the flywheel prototype future design. The presented simplified rotordynamic model can also be applied to rotating machines.



---

*Source: 512674-2013-05-20.xml*