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

This paper shows a theoretical vibration analysis regarding the controller’s parameters and the gyroscopic effect, based on a simplified rotordynamic model. Combined with 600Wh 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.


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 systemfor example, mechanical unbalance and misalignment of the coupling-which may cause vibrations [14][15][16][17][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 International Journal of Rotating Machinery 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.

System Modeling
2.1.System Structure.The basic layout of a flywheel energy storage system is depicted in Figure 1, 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.
Two radial active electromagnetic bearings (RAMBs) including upper and lower RAMBs (RAMB1 and RAMB2 in Figure 1) 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 and -axes which are denoted as  1 and  1 of the upper RAMB and  2 and  2 of the lower RAMB.Moreover, one eddy-current position sensor is installed closely to TAMB to measure the axial position in the -axis which is denoted as .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 rollingelement auxiliary bearings (auxiliary bearing 1 and auxiliary bearing 2 in Figure 1) are necessary to protect AMB stators and stationary components along shaft in the even  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.

Electromagnetic Force Model.
Considering the diameter of the shaft and the loading capacity, the RAMB with 8pole legs was designed in the 600 Wh prototype system.The structural configuration of the RAMB is shown in Figure 2(a), including electromagnetic coils.And the basic magnetic bearing control loop is shown in Figure 2(b).
Taking the two pairs of poles in the -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,  0 is the nominal air gaps of the RAMB in -axes, and  1 is the deviation in -axes.So the upper available air gap is calculated as ( 0 +  1 cos ) and the lower one is calculated as ( 0 −  1 cos ). 0 is the base current and the  1 is the control current. 0 + 1 and  0 − 1 circulate the upper and the lower coils in the -axis, respectively. 1 and  2 are the total attractive electromagnetic forces in the -axis, which is given as follows:  where  0 is the permeability  0 = 4 × 10 −7 H/m,  0 is the effective cross-sectional area of one electromagnet, and  is the number of coils around the core.
Additionally, the total nonlinear attractive electromagnetic forces for the -axis can be modeled as follows: Moreover, by taking the Taylor's expansions of (1) with respect to its nominal operating position ( 1 = 0,  1 = 0), the nonlinear electromagnetic forces can be represented by the following simplified linearized electromagnetic force models: where   and   are the displacement and current stiffness parameters of the RAMB, respectively, and can be obtained from the -axis as follows: It is noted that since the coils in the -axis and the -axis are circulated by the same bias currents  0 and the nominal air gaps in the -axis and the -axis are also the same, that is the displacement and current stiffness parameters  1 and  1 obtained from the -axis are the same as the ones obtained from the -axis (where, the bearing has been assumed to be isotropic).

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 Figure 3, which consists of electromagnet, rotor, displacement sensor, PID controller, D/A digital to analog converter, and A/D analog to digital converter.
The transfer function of the PID controller can be written as where   ,   ,   , and   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: where where   and   are the displacement and current stiffness parameters of the RAMB, respectively;   and   are the amplifiable parameter and the time constant of the power amplifier, respectively;   and   are the amplifiable parameter and the time constant of the sensor, respectively.The parameters   ,   , and   , 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): By the above determination, the derivative time parameter   can be determined firstly, and then the appropriate value can be selected; similarly, the ranges of the proportional parameter   and the differential parameter   can also be gained.According to the effect of the actual control, the integral parameter   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  0 = 0.15 A. Because of the air gap of the magnetic bearing  0 = 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   of the sensor can be obtained as 1 × 10 7 unit/m.The range of the differential time parameter   determined by ( 9) is With   = 0.0001, the range of the proportional parameter determined by ( 13), (14), and ( 9) can be carried out as According to the optimal linear results, the proportional parameter   should be selected as   = 1.And the range of the differential parameter   determined by ( 11) and ( 12) can be given as Therefore, combined with (9), the controller's parameters are finally selected as follows: 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 eddycurrent loss, the saturation and hysteresis of the magnetic core material, and the coupling effects between the electromagnets.The relationship between the center of mass () of the rotor and the five-DOF AMB is shown in Figure 4.
is the body fixed coordinate system with the  axes overlapping the geometric center axis of the electronic bearing stator;  is space fixed coordinate system with the -axes overlapping the geometric center axis of the rotor;  is the center of mass. 1 and  1 are the upper RAMB (RAMB1 in Figure 4) electromagnetic force of the rotor corresponding to the and -axes;  2 and  2 are the lower RAMB (RAMB2 in Figure 4) electromagnetic force of the rotor corresponding to the and -axes;   is the TAMB electromagnetic force of the rotor corresponding to the axes;   is the centrifugal force due to the static imbalance;  1 and  2 represent the distances from the  to the upper RAMB (RAMB1 in Figure 4) and lower RAMB (RAMB2 in the Figure 4), respectively;  1 and  2 represent the distances from the  to the upper position sensor and lower position sensor, respectively;  3 represents the distances from the  to the eccentric mass point;  1 and  1 denote the displacement output of the upper position sensor;  2 and  2 denote the displacement output of the lower position sensor; , , and  denote the pitch, yaw, and spin angles displacements around the -, -, and -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 ., ,  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: Equations (19) are written in a matrix form for simple description as follows: Equation ( 21) is written as  = ,  is the sensor output signal matrix, and  is displacement of the rotor mass center.
Equations (20) are also written in a matrix form for the simple description as follows: Equation ( 22) is directly translated as {} = []{}, where {} denotes the forces vector of the mass center of the rotor as {} denotes the electromagnetic forces vector produced by the AMBs as It is noted that when the rotor is regulated perfectly ( 1 =  1 =  2 =  2 =  = 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  1 ,  2 ,  1 ,  2 , and  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  1 ,  2 ,  1 , and  2 axes and one axial subsystem according to the -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 origin  is zero;  is the mass of the rotor;  is the gravity constant;  are the coordinates of the .The potential energy of the rotor is given as  = .
The total kinetic energy of the system  consists of the translational kinetic energy   , rotational kinetic energy   , and kinetic energy of eccentric mass   .The expression for the translational kinetic energy   is

International Journal of Rotating Machinery
The expression for the rotational kinetic energy   is where   is transverse mass moments of inertia of the rotor, and   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,   is the eccentric mass and   is the eccentricity, so the expression for the kinetic energy of eccentric mass   is There is the displacement from the geometric center of the rotor to the one of the stator.  ,   represent the displacement corresponding to the and -axes, respectively, as Equation ( 28) is substituted to (27) as So the total kinetic energy of the rotor is 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.

𝑞 = [𝑥 𝑎 𝑦 𝛽]
denotes the displacement vector in the mass center of the rotor corresponding to the and -axes.
denotes the displacement vector in the position of the RAMB corresponding to the and -axes.By means of Lagrange equations, which state where  is the total kinetic energy,  is the total potential energy, and   is the generalized force for the coordinate   , a system of four equations describing the dynamics of the model is obtained, which can be written as where  1 and  1 are the upper RAMB electromagnetic force of the rotor corresponding to the and -axes;  2 and  2 are the lower RAMB electromagnetic force of the rotor corresponding to the and -axes.The function of the electromagnetic force is given as follows: where  1 ,  2 are the up and down current stiffness of the active magnetic bearing, respectively. 1 ,  2 are the up and down displacement stiffness of the active magnetic bearing, respectively.The equation is translated as matrix operations for observation as follows where [] is the mass matrix.Because of the small value, the unbalance mass is neglected to make the simple calculation where [] is the inertia matrix, and {  } is the length matrix as where {  } is unbalanced force vector and {  } is the force vector in bearing place as (37) The transfer function block diagram of the PID controller is shown in Figure 6.
Figure 6, the current output depends on the input of the referenced displacement, which states where where   = ( The differential equation of the system can be given as where  =           ,  =         −        ,  =         /  , and When   = 0, that means that the rotor is running steadily, and the differential equation can be written as where   is the differential coefficient of the PID controller and Ḟ  is the force as Then, in order to express this system in the state space form, are defined.
From the above analysis, , , and  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.

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.

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 a  = Ω diagram [21].
Positive and negative values of  correspond to forward and  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  = Ω.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 Table 1, 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  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 Figure 7, 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.
In Figure 7, 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 Figure 8, the integral parameter   is kept at 0.002, and the proportional parameter   is set as 0.5, 1, 1.5, and 2, respectively.With the increase of the proportional parameter   , the forward whirling mode curve is away from the curve  = Ω, the cross-point between the backward whirling curve and the curve  = Ω 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.
In Figure 9, the proportional parameter   is kept at 1, and the differential parameter   is set at 0.001, 0.002, 0.003, and 0.004, respectively.With the increase of the differential parameter   , the forward whirling mode curve is away from the curve  = Ω, the cross-point between the backward whirling curve and the curve  = Ω gradually declined, and the critical speed decreases as well.
The backward whirling modes change distinctly, 3200 RPM in the   = 0.001 and 2000 RPM in the   = 0.004.
The effect of differential parameters on the whirling mode is distinct in running speed 0 ∼ 6000 RPM, but as the running speed rises, the effect of the differential parameter   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 parameter   is selected as 1 and   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.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.Figure 10 is the Campbell diagram of the system with the inertia ratio of the polar moment and the transverse mass moment.
In Figure 11, the   /  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,  /  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 International Journal of Rotating Machinery 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.
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.

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  = Ω 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   /  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  = Ω.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.

Figure 1 :
Figure 1: Basic layout of a flywheel energy storage system.
1 cos   0 −  1 cos  (b) The basic magnetic bearing control loop and its elements

Figure 2 :
Figure 2: Structural configurations and the basic magnetic bearing control loop of the RAMB.

Figure 3 :
Figure 3: Block diagram of the control system.

Figure 4 :
Figure 4: Geometry relationships of rotor and AMB systems.

Figure 5 :
Figure 5: Flywheel rotor with the eccentric mass.

Figure 6 :
Figure 6: Transfer function block diagram of the PID controller.

Figure 7 :
Figure 7: Campbell diagram of the system.

Figure 8 :Figure 9 :
Figure 8: Campbell diagram of the system with the change of the proportional parameter.

Figure 10 :
Figure 10: Campbell diagram of the system with the change of the RAMB positions.

Figure 11 :
Figure 11: Campbell diagram of the system with the inertia ratio of the polar moment and the transverse mass moment.