Output torque modeling of control moment gyros considering rolling element bearing induced disturbances

Mechanical Systems and Signal Processing 115 (2019) 188–212

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Output torque modeling of control moment gyros considering rolling element bearing induced disturbances Hong Wang a, Qinkai Han b,⇑, Daning Zhou a a b

Beijing Institute of Control Engineering, Beijing 100000, China The State Key Laboratory of Tribology, Tsinghua University, Beijing 100084, China

a r t i c l e

i n f o

Article history: Received 18 January 2018 Received in revised form 6 April 2018 Accepted 25 May 2018

Keywords: Control moment gyros Output torque Rolling element bearing Disturbances

a b s t r a c t Control moment gyros (CMG) have been widely used in spacecraft attitude control and large angle slewing maneuvers over the years. Understanding and suppressing highfrequency disturbances in CMG’s output torques is a crucial factor to achieving the desired level of payload performance. Output torque modeling of a single gimbal CMG (SGCMG) with nonlinear rolling bearing supports is conducted in this paper. Taking the installation errors and micro-vibrations of the flywheel into account, three axis output torques of a SGCMG are derived based on Newton-Euler approach and theorem of moment of momentum. Dynamic model is then constructed to obtain the micro-vibration responses of the rotary flywheel. Mass imbalances of the flywheel, flexibility of supporting structures and nonlinearity induced by one pair of angular contact ball bearings are considered in the dynamic model. Especially for the rolling bearing, an improved load distribution analysis is proposed to more accurately obtain the contact deformations and angles between the rolling balls and raceways. Various factors, including the preload condition, surface waviness, Hertz contact and elastohydrodynamic lubrication, are included in the analysis. The bearing restoring forces are then obtained through iteratively solving the load distribution equations at every time step. Dynamic tests on a typical SGCMG supported by angular contact ball bearings are conducted to verify the output torque model. The effects of flywheel dynamic/static eccentricities, inner/outer surface waviness amplitudes, bearing axial preload and installation skew angles on the dynamic output torques are discussed. The obtained results would be useful for the optimal design and vibration control of the SGCMG system. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction A control moment gyro (CMG) consists of a spinning rotor (flywheel) and one (or more) motorized gimbals that tilt the rotor’s angular momentum. As the rotor tilts, the changing angular momentum causes an output torque (also the gyroscopic torque) that rotates the spacecraft. Because of some superior properties, such as large torque amplification [1] and high power efficiency [2], CMG based attitude control systems have been widely used in various space applications [3–5]. When CMGs provide output torque for controlling the attitude of the spacecraft, undesirable high frequency torques with low amplitude appear during these operation due to various micro-vibration disturbances. These high-frequency output torques

⇑ Corresponding author. E-mail address: [email protected] (Q. Han). https://doi.org/10.1016/j.ymssp.2018.05.042 0888-3270/Ó 2018 Elsevier Ltd. All rights reserved.

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

189

can seriously degrade the performance of instruments with high pointing precision and stability [6]. Therefore, understanding and controlling high-frequency output torques of CMGs becomes a significant problem for the modern high-performance spacecraft [7]. As the key component of the CMG, the rotary flywheel produces much of the disturbance, and has dominant effects on the high-frequency output torques of the CMG [8]. Most current research focused on appropriately modeling and analyzing vibrations of the rotary flywheel. Masterson et al. [9] proposed an empirical model to predict the vibration behavior of a wheel assembly in the Hubble Space Telescope. Liu, Maghami and Blaurock [10] considered the axial/rocking modes of the flywheel and gyroscopic torques on the wheel, and then established a dynamic model of the flywheel to analyze the jitter induced from both tonal and broadband disturbances. Narayan, Nair and Ghosal [11] discussed the dynamic coupling problem of a flywheel with its supporting brackets. For a cantilevered flywheel supported by a soft-suspension system, Zhou and Zhang [12,13] presented a comprehensive linear dynamic model to study the micro-vibrations induced by the mass imbalance of wheel and broadband noises. Later, the linear dynamic model has been successfully used for the design of both passive and intelligent flywheel vibration isolation system [14,15]. Zhang, Aglietti and Ren [16] analyzed the coupled microvibrations of a reaction wheel assembly including gyroscopic effects in its accelerant state. Addari, Aglietti and Remedia [17,18] carried out thorough experimental studies on a reaction wheel including the gyroscopic effects. Numerical simulations were also presented in order for comparisons and validations [18]. The above mentioned studies aimed to the flywheels in the reaction and moment wheel assemblies (R/MWAs). Although R/MWAs and CMGs both achieve reaction torque via high speed rotary flywheels, their operating principles are different [8]. R/MWAs provide control torque for spacecraft by adjusting he rotational acceleration of the flywheel, while CMGs generate output gyroscopic torques by gimballing the angular momentum vector of the flywheel. The rotation speed keeps constant in CMGs and it is much higher than that of R/MWAs (as high as 6000–10,000 rev/min) [8]. Thus, some researcher also paid a great deal of effort to study dynamic behaviors and vibration controls of the rotary flywheel in CMGs. Luo, Li and Jiang [19] conducted coupled dynamic analysis of a single gimbal CMG (SGCMG) cluster integrated with an isolation system. Their results showed that the gyroscopic effects produced by the rotary flywheel will stiffen or soften several of the structural modes of the coupled system. Recently, they also presented an innovative work on the optimum design of microvibration isolation for multiple flywheel system of spacecraft [20]. After analyzing the micro-vibration characteristics of the flywheel, a passive isolator for SGCMG using the viscoelastic material was put forward by Shi, Li and Luo [21]. Numerical simulations indicated that disturbances with frequency above 40 Hz can be isolated effectively by the viscoelastic material. Taking the static and dynamic imbalances and installation errors of the flywheel into consideration, Zhang and Zhang [22] conducted the disturbance characteristic analysis of CMGs. The constraint conditions of the CMG performance indexes were obtained to meet the requirements of attitude precision and attitude stabilization. Based on the dynamic model of the CMGs in a pyramid configuration considered static and dynamic imbalances of the flywheel, Zhang et al. [23] presented an activepassive integrated strategy for vibration control of CMGs. Through detailed simulations in both time and frequency domains, it was shown that remarkable attitude control improvement could be achieved by applying the proposed control strategy. Most papers mentioned above considered the modeling of mass imbalances of the flywheel and flexibility of supporting structures and their effects upon the micro-vibrations of flywheel assemblies in CMGs or R/MWAs [10–23]. Rolling element bearings are the core supporting component, and play a decisive role in the performance, operation reliability and service life of CMGs. Currently, a linear spring-damper model has been widely adopted to describe the internal compliances of rolling element bearings [12–16,19–23]. Such treatment is obviously too ideal and could not reflect the imperfections and disturbances accurately. Some key factors, including the Hertzian contacts, bearing radial clearance, surface waviness, preload condition and so on, would significantly affect the dynamic behaviors of rolling element bearings and then the output torques of CMGs. Zhou et al. [24] first introduced the bearing irregularity and nonlinear stiffness into the dynamic model of a wellbalanced flywheel. Simulation results indicated that the dynamic amplification is even greater than the mass imbalance when the bearing noise intersects with the translation mode at high rotational speed. Considering the nonlinear compliance of the internal bearing systems and the coupled motions of the rotating flywheel and the gimbal, Luo et al. [25] employed the energy method to develop a theoretical output torque model of a SGCMG. High frequency disturbances were found in the time history of output torques predicted by the theoretical model. In their model [24,25], both the nonlinear stiffness and bearing irregularity induced by the surface waviness were given by empirical equations, which might reduce the model accuracy. Currently, most researchers believed that load distribution analysis should be utilized to iteratively solve the nonlinear bearing supporting forces. By applying elastic/geometric constraints between a rolling ball and radial roller bearings, Jones [26] first proposed a general load distribution model (Jones model) with five degrees of freedom considering gyroscopic moments and centrifugal forces of the balls. Tiwari et al. [27,28] greatly improved Jones’ model by taking into account both radial clearance and Hertzian contact characteristics. Yhland [29,30] first pointed out that the waviness is the global imperfection on the surfaces of the bearing components, and conducted experimental measurements of waviness. Based on the work of Yhland [29,30], Jang and Jeong [31–34] introduced the bearing waviness into the Jones’ model and carried out a series of work to analyze the effects of bearing waviness on the vibration forces and frequencies and on the stability of the rotating system. Later, Bai and Xu [35] proposed an improved load distribution model to study the effect of bearing waviness on the fluctuation of the cage speed. Cao and Xiao [36] developed a comprehensive model of the double-row spherical roller bearing. In their model, several excitations such as the waviness, clearance, surface defects were included. Babu et al. [37] extended the Jones’ model to six degrees of freedom, and found that the effect of frictional moment on the vibrations of

190

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

the system was significantly dependent on the bearing load. Cao et al. [38] and Niu et al. [39] proposed a new dynamic model of ball-bearing rotor systems based on rigid body element method, and conducted detailed investigations on the occurrence of stable cage whirl motions in ball bearings. Zhang et al. [40–42] considered the rotor imbalance excitation as a type of bearing preloads, and showed the variations of dynamic stability and response behaviors due to the rotor imbalance. The authors [43,44] presented a nonlinear dynamic model for angular contact ball bearings, and then applied the model for microvibration analysis of MWAs. Very recently, Cao et al. [45] conducted a comprehensive review on the mechanical model development of rolling bearing-rotor systems, and pointed out that the bearing dynamic model becomes comparatively mature after nearly sixty years development. Various effects, including the Hertzian contacts, bearing radial clearance, surface waviness and preload condition, could be considered in detail. Currently, most bearing dynamic models utilize the displacement coordinate relations to approximately obtain the contact deformations between the rolling element and raceways, and then the nonlinear restoring forces of the rolling bearing could be calculated accordingly. Although the calculation efficiency is relatively higher, the accuracy is lower as the actual contact deformations should be solved through iterative analysis. Thus, output torque modeling of a SGCMG with nonlinear rolling bearing supports is conducted in this paper. Taking the installation errors and micro-vibrations of the flywheel into account, the output torque of a SGCMG is derived according to the Newton-Euler approach and theorem of moment of momentum. Dynamic model is then constructed to obtain the microvibration responses of the rotary flywheel. Mass imbalances of the flywheel, flexibility of supporting structures and nonlinearity induced by one pair of angular contact ball bearings are considered in the dynamic model. Especially for the rolling bearing, an improved load distribution analysis is proposed to more accurately obtain the contact deformations and angles between the rolling balls and raceways. Various factors, including the preload condition, surface waviness, Hertz contact and elastohydrodynamic lubrication, are included in the analysis. The bearing restoring forces are then obtained through iteratively solving the load distribution equations at every time step. Dynamic tests on a typical SGCMG supported by angular contact ball bearings are conducted to verify the output torque model. The effects of flywheel eccentricity and inner/outer waviness amplitudes on the dynamic output torques are discussed in detail. Finally, some conclusions are given. 2. Output torque modeling of a SGCMG The structure of a small SGCMG, as shown in Fig. 1, usually includes the high-speed flywheel, bearing assemblies, housing and mounting shell, low-speed gimbal and its supporting structures. A pair of high-precision angular contact ball bearings (see the left and right bearings in the figure) with face-to-face installation is fixed on the stationary spindle. The high-speed flywheel is fixed with the bearing’s out rings through the mounting shell, and driven by a DC motor in the working state. One end of the stationary spindle is connected with the mounting brackets, and the other end is supported by the mounting shell. A four-contact-point ball bearing is utilized to support the low-speed gimbal, which is connected with the high-speed flywheel through the mounting brackets. The gimbal support structure is fixed with the satellite. Under the drive of the motor, the gimbal rotates at low speed and then produces output torques to control the attitude of satellites. In order to derive the output torque of the SGCMG, a series of coordinates should be defined first, as shown in Fig. 2. The center of f g ðX g ; Y g ; Z g Þ is at the mass center of gimbal, Z g axis coincides with the rotating axis of the gimbal, X g axis is perpendicular to the Z g axis, but it does not coincides with the rotating axis of the flywheel. The coordinates f w ðX w ; Y w ; Z w Þ and f f ðX f ; Y f ; Z f Þ are defined for the rotary flywheel, and their centers are all at the geometric center of flywheel. The X f axis and X w axis are all coincides with the rotating axis of flywheel. The f w ðX w ; Y w ; Z w Þ rotates with the flywheel, while the f f ðX f ; Y f ; Z f Þ remains stationary. In the initial state, the f w ðX w ; Y w ; Z w Þ coincides with f I ðX I ; Y I ; Z I Þ, which represents the inertial coordinate of the flywheel. Its center is at the mass center of flywheel, and X I ; Y I denote the directions of transversal moments of inertia. In the f I coordinate, the vector of moment of inertia of the flywheel could be written as I wI ¼ diagð½Ipfw ; Idfw ; Idfw Þ, where Ipfw ; Idfw denote the polar and transversal moments of inertia of the flywheel. Considering the micro-vibrations of flywheel, the vector of moment of inertia changes in the f w coordinate. Applying the coordinate transformation, it is expressed as follows

I ww ¼ ATIw I wI AIw

ð1Þ

in which AIw is the transformation matrix from f w to f I , the superscript ‘‘T” denotes the matrix transposition. Three Euler angles, nfw ; gfw ; /fw , are utilized to describe the micro-vibrations and self-rotation of the flywheel. Three continuous rotations are needed to make the f w parallel to the f I . The specific rotation sequence is shown in Fig. 3(a). First, it is rotated by angle gfw around axis Z w , then by angle nfw around the new axis Y w1 and last by angle /fw around the new axis X w2 . Thus, AIw could be gained by the product of coordinate transformation matrices

2

1

0

0

32

cosðnfw Þ 0  sinðnfw Þ

6 76 0 AIw ¼ Ax ð/fw ÞAy ðnfw ÞAz ðgfw Þ ¼ 4 0 cosð/fw Þ  sinð/fw Þ 54 0 sinð/fw Þ cosð/fw Þ sinðnfw Þ

1 0

32

cosðgfw Þ  sinðgfw Þ 0

76 54 sinðgfw Þ cosðnfw Þ 0 0

cosðgfw Þ 0

3

7 05

1 ð2Þ

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

191

Fig. 1. Schematic diagram for a SGCMG with angular contact ball bearing supports.

As the micro-vibration of flywheel always has low amplitudes, so the Euler angles nfw ; gfw are assumed to be small quantities. For the flywheel with constant rotational speed X, one has /fw ¼ Xt. Substituting Eq. (2) to Eq. (1) and ignoring the second-order small quantities, the expression of I ww is gained as

2

I ww

6 Ipfw 6 6 6 4 Sym:



gfw cos Xt

 ðIpfw  Idfw Þ

nfw sin Xt Idfw

nfw cos Xt

!

3

ðIdfw  Ipfw Þ 7 7 þgfw sin Xt 7 7 0 5 Idfw

ð3Þ

In the actual applications, due to the manufacturing and installation errors, the gimbal axis does not coincide with the flywheel axis. Three Euler angles, agf ; bgf ; vgf , are defined to describe the angular skewness between the gimbal and flywheel. Fig. 3(b) shows that three continuous rotations should be done to make the f g parallel to the f f . First, it is rotated by angle agf around axis Z g , then by angle bgf around the new axis Y g1 and last by angle vgf around the new axis X g2 . As the skewness angles due to installation errors are also small quantities, so the transformation matrix from f g to f f could be approximately expressed by

2

1

6 Agf ¼ Ax ðvgf ÞAy ðbgf ÞAz ðagf Þ  4 agf bgf

agf 1 vgf

bgf

3

vgf 7 5 1

ð4Þ

192

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

Fig. 2. Coordinates used in the analysis.

Fig. 3. Specific rotation sequence: (a) from f w to f I ; (b) from f g to f f .

Considering the micro-vibrations and installation errors, the absolute angular momentum of the flywheel is obtained through the above coordinate transformations

Hgw ¼ I ww Agf X

ð5Þ

in which X ¼ ½X; 0; 0T denotes the angular speed vector of the flywheel. According to theorem of moment of momentum, disturbance torque vector of the flywheel in the gimbal coordinate f g is expressed as follows

Td ¼ Agg

d ðHgw Þ dt

ð6Þ

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

193

where Agg ¼ ½0; 0; 1; 0; 1; 0; 1; 0; 0. For given gimbal rotational speed xg , under the ideal condition (without disturbances induced by micro-vibrations of flywheel and installation errors, etc.), the output torque vector for the rotary flywheel could be expressed as ½0; Ipfw Xxg ; 0T . In our dynamic tests, the output torque is acquired in the fixed frame coordinate. Thus, in the derivation of total output torques of the SGCMG, the influence of gimbal rotation should also be taken into account. In the fixed frame coordinate for the gimbal, considering the disturbance torques induced by the micro-vibrations and installation errors, the total output torque vector for the SGCMG could be written by

2

Tout

3

2 cos xg t d 6 7 6 ¼ Td þ Az ðxg tÞ4 Ipfw Xxg 5 ¼ Agg ðHgw Þ þ 4 sin xg t dt 0 0 0

 sin xg t cos xg t 0

0

32

0

3

7 76 0 54 Ipfw Xxg 5 0 1

ð7Þ

_ ¼ 0. Substituting Eqs. (5) and (6) into Eq. (7) and simplifying, the The rotation speed of the flywheel is constant, i.e. X expressions for the output torque components of the SGCMG are written by

T out;x ¼ KX2 þ Ipfw Xxg sin xg t

ð8Þ

T out;y ¼ NX2 þ Ipfw Xxg cos xg t

ð9Þ

T out;z ¼ ðagf N þ bgf KÞX2

ð10Þ

where N; K are determined by the micro-vibrations of the flywheel

N ¼ ðgfw sin Xt  nfw cos XtÞðIpfw  Idfw Þ

ð11Þ

K ¼ ðnfw sin Xt þ gfw cos XtÞðIdfw  Ipfw Þ

ð12Þ

From Eqs. (11) and (12), one can find that N ¼ K ¼ 0 if the micro-vibrations and installation errors are ignored (gfw ¼ nfw ¼ 0 and agf ¼ bgf ¼ 0). In this case, the output torque components of the SGCMG could be rewritten as

Tb out;x ¼ Ipfw Xxg sin xg t

ð13Þ

Tb out;y ¼ Ipfw Xxg cos xg t

ð14Þ

Tb out;z ¼ 0

ð15Þ

Obviously, Eqs. (13)–(15) describe the output torques of the SGCMG under the ideal condition. However, due to the influences of micro-vibrations and installation skewness, disturbances appears in the output torques, as shown in Eqs. (8)–(10). In order to compute the total output torques, the skew angles abf ; bbf ; vbf and micro-vibration responses nfw ; gfw should be gained first. The former belongs to the input parameter, while the latter is time-variable and closely related to the mass imbalances and bearing induced disturbances. The values of nfw ; gfw could be gained by numerically solving dynamic differential equations of the rotary flywheel. In the following section, the dynamic behaviors of high-speed rotary flywheel will be analyzed to determine the time-variable responses nfw ; gfw . 3. Dynamic analysis of high-speed flywheels Under the working state, the rotational speed of flywheel in SGCMG is higher than 5000 rev/min, while the maximum speed of the gimbal is no more than 10 rev/min. As the difference of working speed between the flywheel and gimbal is distinct, so the dynamic coupling between them is weak from the view point of vibration. The low-speed gimbal mainly supports the high-speed flywheel. Therefore, in the analysis of output torque of SGCMG, the low-speed gimbal is assumed to be rigid, and the high-speed flywheel is regarded as a rigid-flexible coupling dynamic system. Such system is usually described by a series of second order differential equations, which will be constructed in this section. Then, the influence of the microvibrations on the output torque of the SGCMG could be analyzed accordingly. In the modeling of micro-vibrations of high-speed flywheel, three subsystems, including the rotor, rolling bearing assembly, stationary spindle and supporting structures, are considered to established their nonlinear mass-spring lumped parameter models, as shown in Fig. 4. The detailed information for these subsystems is described as follows: (1) The actual flywheel rotor has complex wheel-spoke type of structures. In our analysis, the structural vibrations of the rotor itself are ignored. A lumped mass with five degrees of freedom in the coordinate f I ðX I ; Y I ; Z I Þ, including three translations and two rotations around the Y I ; Z I axes, is utilized to describe the three-dimensional movements of the rotor.

194

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

Fig. 4. Vibration model for the high-speed flywheel.

(2) Under the working condition, the rotor’s imbalance disturbance force is transmitted to the spindle through the rolling element bearings, and then it is transmitted to the satellite by the mounting bracket and the low-speed gimbal. At the same time, the internal excitation forces can also be generated due to various micro factors, such as the Hertz contact, elastohydrodynamic lubrication and surface waviness in the rolling bearings. In this study, the internal excitations and disturbances induced by the rolling bearings are considered by establishing the nonlinear restoring forces/torques in five directions, i.e. disturbing forces at three translational directions and disturbing torques along two rotational directions (except for the torques along the bearing’s rotational axis). Fig. 4 also shows these nonlinear forces and torques. For the left and right bearings, the forces and torques are denoted by f bLi ; f bRi , in which i ¼ 1; 2; 3; 4; 5. (3) Two lumped masses ms and mf are used to simulate the inertial properties of the stationary spindle and mounting bracket, respectively. Similar to the flywheel rotor, under the coordinate f I , the lumped mass ms also has three translational motions and two rotational motions. These motions are coupled with that of the flywheel rotor through bearing’s nonlinear disturbing forces/torques f bLi ; f bRi . In order to simplify the analysis, the lumped mass mf is considered with only three translational degrees of freedom. It is coupled with the ms through three linear springs, whose stiffnesses are denoted by kx1 ; ky1 ; kz1 . The lumped mass mf is connected to the satellite through the gimbal and its supporting structures. The flexibility of the supporting structures can be simplified as three linear springs, whose stiffnesses are expressed by kfx ; kfy ; kfz . In addition, the right end of the stationary spindle is usually supported by the mounting shell, and the supporting stiffnesses along two translational directions are denoted by ky2 ; kz2 . In the following section, the equations of motion for these subsystems will be, respectively, derived and then the dynamic model of high-speed flywheel will be gained through assembling subsystems’ equations of motion.

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

195

3.1. Rotor vibration model The structural vibrations of the rotor are neglected. It is simplified as a rigid mass with five degrees of freedom. Fig. 5 shows the vibration model of the rotor. The coordinate system is f I ðX I ; Y I ; Z I Þ. The rotor has three translations along the axes (ufw ; v fw ; wfw ) and two rotations around the Y I and Z I axes (nfw ; gfw ). According to the equilibrium condition of forces, one can derive the equations of motion for the rotor as follows

€ fw ¼ f bL1 þ f bR1 mfw u mfw v€ fw ¼ f u1 þ f bL2 þ f bR2  mfw g € fw ¼ f u2 þ f bL3 þ f bR3 mfw w Idfw €nfw þ Ipfw Xg_ fw ¼ T u1 þ f bL4 þ f bR4 € fw  Ipfw Xn_ fw ¼ T u2 þ f bL5 þ f bR5 Idfw g

ð16Þ

where mfw , Idfw , Ipfw the mass, diametral and polar moments of inertia of the rotor, f bLi , f bRi (i ¼ 1; 2; 3; 4; 5) the nonlinear restoring forces and torques of the left and right bearings, which will be solved in the following section. The f u1 ; f u2 ; T u1 ; T u2 in Eq. (16) are forces and torques induced by the unbalance mass, which is inevitable due to the nonuniform material, processing and assembly of the rotor. Both static and dynamic imbalances are considered, and shown in Fig. 5. For the static imbalance, time-variant forces f u1 and f u2 will be induced along the Y I and Z I axes. Their expressions are given by

f u1 ¼U s X2 cosðXt þ us Þ f u2 ¼U s X2 sinðXt þ us Þ

ð17Þ

in which U s ¼ mms rms the static imbalance mass, us the initial phase. For the dynamic imbalance, time-variant torques T u1 ; T u2 will be induced around the Y I and Z I axes. One can get their expressions as follows

Fig. 5. Vibration model for the rotor: (a) coordinate system; (b) static imbalance; (c) dynamic imbalance.

196

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

T u1 ¼U d X2 cosðXt þ ud Þ

ð18Þ

T u2 ¼U d X2 sinðXt þ ud Þ

in which U d ¼ mmd rmd hf the dynamic imbalance mass, ud the initial phase of the dynamic imbalance mass, and hf the axial length between the imbalance masses. 3.2. Nonlinear restoring forces of rolling bearings The inner race of the angular contact ball bearing is fixed on the stationary spindle, and the outer race is connected to the rotor and rotates under constant speed X. The mass unbalance excitations are transmitted through rolling bearings to the spacecraft. With the transmission of unbalanced excitations, internal excitations of rolling element bearings are induced by the Hertzian contact, surface waviness, elastohydrodynamic lubrication and other micro disturbances. These excitations will generate vibrations and then influence the dynamic behaviors of the flywheel. In the study, the internal excitations are considered by modeling the nonlinear restoring forces and moments of the angular contact ball bearing. Besides, appropriate value of preload should be applied to maintain the effective contact between the rolling element and the raceways, and then improve the bearing rotation accuracy and bearing stiffness. Thus, load distribution analysis is conducted by taking preload condition and various micro disturbances into account. Through nonlinear iterative solutions, both the contact deformations and angles between the rolling ball and inner/outer raceways are obtained. Finally, the nonlinear restoring forces and torques of the bearings are computed utilizing the Hertz contact and oil lubricant theorems. 3.2.1. Load distribution analysis Fig. 6 gives the global coordinate of the bearing (X  Y  Z) and the local coordinate of the jth ball (xj  yj  zj ), of which the former is fixed at the rotating axis center of the bearing and the later is rotating around the Z axis under ball’s orbital speed xcj . Except for the torsional motion around the Z axis, the other five degrees of freedom of the outer race, including two translations, two rotations and one axial motion, are considered in the model. The vector of degree of freedom is expressed by u ¼ ½u; v ; w; hu ; hv  . Accordingly, the preloads along the five directions are f p ¼ ½f u ; f v ; f w ; mu ; mv  , in which f u ; f v ; f w are three forces and mu ; mv are two torques around the X; Y axes. Usually, apart from the axial force (f w > 0), the preloads along the other four directions are all be equal to zero, i.e. f u ¼ f v ¼ 0; mu ¼ mv ¼ 0. The structural flexibility of outer race is neglected. The external forces cause the outer race deform as a rigid body. For the groove curvature center of the outer race corresponding to the jth ball, its rigid body deformation vector in the local coorT

T

dinate is expressed by dj ¼ ½fj ; gj ; hj T . Under the assumption of small deflection, the dj could be obtained from the rigid deformation of the outer race mass center u through coordinate transformations

dj ¼ Tj u þ Dwj

ð19Þ

Fig. 6. Schematic diagram for the bearing coordinates.

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

197

where Tj is the transformation matrix, and is expressed as

2

6 Tj ¼ 4

cos /cj

sin /cj

0

0

0

0

1 Ro sin /cj

0

0

0

sin /cj

0

3

Ro cos /cj 7 5

ð20Þ

 cos /cj

in which /cj ¼ /0 þ xcj t þ 2Np ðj  1Þ the location angle of the jth ball at time t, where /0 and N b are the initial angle and numb

ber of rolling ball. Without loss of generality, one has /0 ¼ 0. The Ro ¼ 0:5de  ðr go  0:5db Þ cos b0 denotes distance from the rotating axis center of outer race to the groove curvature center, de is the pitch diameter of the bearing, rgo is the groove curvature radius, db is the ball diameter and b0 is the nominal contact angle. Besides the rigid displacement of the outer race, the surface waviness would also induce additional displacements Dwj , as shown in Eq. (19). Waviness is the global sinusoidally shaped imperfection on the surfaces of the bearing components. In general, these imperfections are due to the irregularities during the grinding and honing process of the bearing components. Nowadays, the amplitude of the waviness in rolling element bearings is at the micrometer level. Despite that, waviness can still produce significant vibrations. Fig. 7 shows the schematic diagram of the waviness appeared in the inner/outer raceways and rolling balls. In the figure, pij ; poj and qij ; qoj denote circumferential and axial waviness of inner and outer raceways when contact with the jth ball. The wij ; woj represent the rolling ball’s waviness. These values all periodically vary with time, and could be expressed using the sum of a series of harmonic functions

pij ¼ qij ¼

X

X

poj ¼ qoj ¼ wij ¼ woj ¼

Ain cosðnin xcj t þ 2pðj  1Þ=Nb þ win Þ

ð21Þ

Bin cosðnin xcj t þ 2pðj  1Þ=Nb þ uin Þ

ð22Þ

X X

Aout cosðnout ðxo  xcj Þt þ 2pðj  1Þ=N b þ wout Þ

ð23Þ

Bout cosðnout ðxo  xcj Þt þ 2pðj  1Þ=Nb þ uout Þ

ð24Þ

C n cosðnba xsj t þ #ba Þ

ð25Þ

C n cosðnba xsj t þ pi þ #ba Þ

ð26Þ

X X

where Ain ; Bin ; Aout ; Bout ; C n ; nin ; nout ; nba and win ; uin ; wout ; uout ; #ba are, respectively, denote the amplitudes, harmonic orders and phases of the waviness, xsj is the self-rotation angular speed. Considering the fixed inner race and the pure rolling state, one   2  b0 e 1  db cos xo , and orbital speed of the ball could derive the expressions of ball’s self rotation speed xsj ¼ 0:5d db de   db cos b0 xcj ¼ 0:5 1 þ de xo . According to the geometric relations in Figs. 6 and 7, it is easy to obtain

2

poj  pij

3

6 7 Dwj ¼ 4 qoj  qij 5

ð27Þ

0

Fig. 7. Schematic diagram for the bearing waviness.

198

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

Without load, the ball just contacts with the inner and outer races. As the deflection is zero, so the ball center and groove curvature centers of inner and outer races lie in a straight line, which is shown as the dash line in Fig. 8. The contact angles for the ball with both inner and outer races are the same and equal to the nominal contact angle b0 . As long as the load is applied, the outer race will first move and then the ball also moves due to the contact deflections. Equilibrium state, shown by the solid line in Fig. 8, could be achieved between the ball and inner/outer races. As the inner race is fixed, so its groove curvature center keeps unchanged. From Eq. (19), the displacements of the groove curvature center of inner race are expressed by (fj ; gj ). The displacements of ball center are unknown parameters and denoted by (v jy ; v jz ). Due to the load on inner race, the inner and outer contact angles become different. The outer contact angle reduces to boj , while the inner contact angle increases to bij . Before deflection, the distance from the groove curvature center of inner race to the ball mass center is Lij . After deflection, the distance becomes lij , as shown in Fig. 8. Similarly, the distance from the groove curvature center of outer race to the ball mass center changes from Loj (before deflection) to loj (after deflection). For the given groove curvature radii of the inner/ outer races (r gi ; rgo ), one has

Lij ¼ rgi þ pij  hci  ð0:5db þ wij Þ

ð28Þ

Loj ¼ r go þ poj  hco  ð0:5db þ woj Þ

ð29Þ

in which hci ; hco are the oil film thicknesses for the rolling ball contacting with the inner and outer raceways, respectively. According to the result of Hamrock and Dowson [46], the oil film thickness is calculated by the following empirical equation

hc ¼ 2:69Rx U 0:67 G0:53 W 0:067 ð1  0:61e0:73j Þ where U ¼

g0 uent E0 Rx

ð30Þ

denotes the dimensionless speed, in which g0 is the lubricant viscosity at the atmospheric pressure and tem-

perature of 200 C; Rx ¼ db =2ð1  db =de cos b0 Þ represents the equivalent radius along the ball rolling direction. For the outer dro race, the ‘‘þ” is selected. The ‘‘” is chosen for the inner race. The equivalent rotating speed is uent ¼ 2ðddriþd , where the ro Þ ri

dri ; dro denote the diameters of inner and outer races. The G ¼ E0 cgp is the dimensionless elastic modulus, in which cgp the Q

viscosity-pressure coefficient. The dimensionless load is written by W ¼ E0 Rj2 . Obviously, the oil film thickness hc is related x

to the contact force Q j . Thus, the value of hc varies with time due to the variation of Q j . In the following dynamic analysis, such variation will be considered in the iteration. According to the geometric relations shown in Fig. 8, one can obtain

tan bij ¼

Lij sin b0 þ v jz Lij cos b0  v jy

ð31Þ

tan boj ¼

Loj sin b0 þ gj  v jz Loj cos b0  fj þ v jy

ð32Þ

lij ¼

Lij sin b0 þ v jz sin bij

ð33Þ

Fig. 8. Geometric relations for the bearing before and after deflection.

199

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

loj ¼

Loj sin b0 þ gj  v jz sin boj

ð34Þ

Eqs. (31)–(34) are geometric equations, which should be satisfied in the load distribution analysis. Besides the geometric equations, the force equilibrium equations for the rolling ball and outer race should also be established. Fig. 9 gives the forces applied on the jth rolling ball. In the figure, F cj and M gj denote the centrifugal and gyroscopic forces due to the rotation and revolution of the ball. Obviously, one has F cj ¼ 0:5mb de x2cj and M gj ¼ Ib xsj xcj sin aj , where mb ; Ib the mass and moment of inertia of the ball, and aj the angle between self-rotation axis and the z axis. For the pure rolling state, one has   aj ¼ tan1 cos bsin0 þdb0 =de . b

The contact forces between the ball and inner/outer races are denoted by Q ij and Q oj . According to the Hertz contact theand Q oj ¼ voj K o d3=2 ory, one has Q ij ¼ vij K i d3=2 ij oj , in which K i ; K o and dij ; doj the contact stiffness coefficients and contact deflections between the ball and innner/outer races, respectively. From the above geometric analysis, one has dij ¼ lij  Lij and doj ¼ loj  Loj . When dij ; doj > 0; vij ¼ 1; voj ¼ 1. Otherwise, vij ¼ 0; voj ¼ 0 for dij ; doj  0. The contact stiffness coefficients are pffiffiffiffiffiffiffiffiffiffiffiffiffiffi determined by geometric dimensions of contact zone and material properties, i.e. K ¼ ðpjE0 =3n 3R=nÞ, in which E0 ¼ E=ð1  m2 Þ denote the effective elastic modulus, and E the material elastic modulus, m the poisson’s ratio, R the equivalent curvature radius, j the elliptical ratio, n;  the first and second kinds of elliptical integration. From the monograph of Harris [47], the expressions for these parameters are given and the contact stiffness coefficients K i ; K o could then be gained accordingly. In the load distribution analysis, the friction force out of the plane is ignored. The gyroscopic moment just equals to the moment induced by the friction force in the y  z plane, as shown in Fig. 9. kij ; koj denote the friction coefficients between the ball and inner/outer races. Here, the values of kij ; koj are set to be 1. From Fig. 9, the force equilibrium equations are obtained as follows

Q ij cos bij  Q oj cos boj þ

kij Mgj koj Mgj sin bij  sin boj þ F cj ¼ 0 db db

Q ij sin bij þ Q oj sin boj þ

ð35Þ

kij M gj koj M gj cos bij  cos boj ¼ 0 db db

ð36Þ

Eqs. (35) and (36) show the equilibrium state for the forces applied on the jth rolling ball. For the other rolling balls, similar force equilibrium equations could be expressed. Not only the rolling ball is in equilibrium state, but also the outer race in analyzing the load distribution of rolling bearings. The jth ball provides contact force Q oj and friction force

koj Mgj db

to the outer

race, and the directions are opposite with the directions shown in Fig. 9. By transforming these forces to the groove curvature center of the outer race, one has

3 k M Q oj cos boj þ ojd gj sin boj b 7 6 k M 7 6Q 7 6 4 zj 5 ¼ 6 Q oj sin boj þ ojdb gj cos boj 7 5 4 koj M gj Q hxj r go d 2

Q yj

3

2

b

Fig. 9. Forces applied on the jth rolling ball [47].

ð37Þ

200

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

Thus, the force equilibrium equations of the inner race could be expressed as

2 3 Q yj Nb X  T 6 7 fb ¼ Tj 4 Q zj 5 j¼1 Q hxj

ð38Þ T

where f b ¼ ½f bu ; f bv ; f bw ; mbu ; mbv  denotes the resultant force vector of the outer race due to all of the contacting balls, which is just the restoring force vector of the rolling bearing. Obviously, the following equilibrium equation exists for the outer race

fp  fb ¼ 0

ð39Þ

Eqs. (35), (36) and (39) are a set of nonlinear algebraic equations, which could be solved iteratively by the NewtonRaphson method. Obviously, the dimension of the nonlinear equations is 2N b þ 5. After the unknown parameters are gained, the contact deformations (dij ; doj ) and inner/outer contact angles (bij ; boj ) for every rolling ball could be computed by Eqs. (19)–(31), (20)–(34). Moreover, the displacement vector uo of outer race under given values of preload could be obtained accordingly. 3.2.2. Bearing forces and moments Fig. 10 shows a pair of angular contact ball bearings (left and right bearings). The X I  Y I  Z I denotes the flywheel rotor coordinate, and the X bL  Y bL  Z bL and X bR  Y bR  Z bR are, respectively, the analysis coordinates for left and right bearings. The axial distances between the two bearings and the flywheel mass center are l1 ; l2 . From the above descriptions, the displacement vector for the rotor is ½ufw ; v fw ; wfw ; nfw ; gfw T . The stationary spindle is represented by a lumped masses ms , and the

displacement vector is expressed as ½us ; v s ; ws ; hvs ; hw s  . As the inner races of the rolling bearings are rigidly connected to the spindle, so the displacements of the inner race could be determined through coordinate transformations. Considering the spindle’s deflection and pre-deflection of outer race under given preloads, the vibrational displacements of outer race of left and right bearings could be gained using coordinate transformation technique T

02 3 3 2 31 ufw us ubL B6 v 7 6 7 C 7 6 B6 fw 7 6 v s 7C 6 v bL 7 B6 7 7 6 7C 6 B6 7 7 6 7C 6 6 wbL 7 ¼ TL B6 wfw 7  6 ws 7C þ uo B6 7 7 6 7C 6 B6 n 7 6 hv 7C 6 hv 7 @4 fw 5 4 s 5A 4 bL 5 2

hw bL 2

ubR

hw s

gfw 3

02

ufw

ð40Þ

3

2

us

31

B6 v 7 6 7C 7 6 B6 fw 7 6 v s 7C 6 v bR 7 B6 7 7 6 7C 6 B6 7 7 6 7C 6 w ¼ T 6 bR 7 U B6 wfw 7  6 ws 7C þ uo B6 7 7 6 7C 6 B6 n 7 6 hv 7C 6 hv 7 @4 fw 5 4 s 5A 4 bR 5 hw bR

ð41Þ

hw s

gfw

v w where ½ubL ; v bL ; wbL ; hvbL ; hw bL  and ½ubR ; v bR ; wbR ; hbR ; hbR  , respectively, denote the displacement vectors of outer race of left and right bearings. According to the definition in Fig. 10, the transformation matrices TL ; TU are derived as follows T

T

Fig. 10. Bearing coordinates for computing the nonlinear forces and moments.

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

2

1

6 60 6 TL ¼ 6 60 6 40

0

0

0

1

0

0

0

1

l1

0

1=l1

1

0

0

0 1=l1

0

3

7 l1 7 7 0 7 7; 7 0 5 1

2

1

6 6 0 6 TU ¼ 6 6 0 6 4 0 0

0

0

0

1

0

0

0

1

l2

0

1=l2

1

1=l2

0

0

0

201

3

7 l2 7 7 0 7 7 7 0 5 1

ð42Þ

After the outer race displacements for the two bearings are obtained, the corresponding restoring forces could be gained through load distribution analysis presented in the previous section. In the following, the left bearing is taken as an example to show the solution process for nonlinear bearing restoring forces in detail. The displacements of jth ball of left bearing dLj could be gained by substituting the outer race displacement vector of left bearing uL ¼ ½ubL ; v bL ; wbL ; hvbL ; hw bL  into Eq. (19). Similarly, the force equilibrium equations for every rolling balls of left bearing are obtained, as shown in Eqs. (35) and (36). Utilizing the Newton-Raphson method, such nonlinear algebraic equations are solved iteratively, and the displacements of every ball in the local coordinate are gained. According to Eqs. (31)–(34) and T

Eq. (37), the reaction force vector from the jth ball to the outer race of left bearing is ½Q Lyj ; Q Lzj ; Q Lhxj T . Thus, in the system coordinate, the sum of reaction forces from all of the rolling ball, just the restoring forces of left bearing, are expressed as follows

3 f bL1 3 2 7 6 Q Lyj 6 f bL2 7 Nb 7 6 X 7 6   T 7 6 T 7 Tj 6 6 f bL3 7 ¼ ½TL  4 Q Lzj 5 7 6 j¼1 7 6 Q Lhxj 4 f bL4 5 2

ð43Þ

f bL5 Similarly, the restoring force vector of right bearing could also be gained as follows

2

f bR1

3

3 2 7 6 Q Ryj 6 f bR2 7 Nb 7 6 X 7  T 6 7 6 T 7 Tj 6 6 f bR3 7 ¼ ½TU  4 Q Rzj 5 7 6 j¼1 7 6 Q Rhxj 4 f bR4 5

ð44Þ

f bR5 3.3. Stationary spindle and support structures The stationary spindle, housing and mounting brackets are used to support the flywheel rotor, as shown in Fig. 1. In order to consider the coupling effect and structural flexibilities of these supports, two lumped masses (ms ; mf ) and a series of springs with linear stiffness are adopted in the vibration model, as shown in Fig. 4. The stationary spindle is simplified by a lumped mass ms with five degrees of freedom, including three translations us ; v s ; ws and two rotations hvs ; hw s . In the coordinate X I  Y I  Z I , the ms is nonlinearly coupled with the flywheel rotor through nonlinear bearing restoring forces and torques. Besides, two ends of stationary spindle are supported by kx1 ; ky1 ; kz1 and kx2 ; ky2 , respectively. The mounting brackets is modeled by mf with three translational degrees of freedom. The support stiffnesses of mf are kfx ; kfy ; kfz , and the coupling stiffnesses with the ms are denoted by kx1 ; ky1 ; kz1 . According to the force equilibrium condition, the equations of motion of the mass ms could be derived as

€ s þ kx1 us  kx1 uf ¼ f bL1  f bR1 ms u ms v€ s þ ðky1 þ ky2 Þv s þ ðky2 l2  ky1 l1 Þhvs  ky1 v f ¼ f bL2  f bR2  ms g € s þ ðkz1 þ kz2 Þws þ ðkz1 l1  kz2 l2 Þhw ms w s  kz1 wf ¼ f bL3  f bR3

ð45Þ

2 2 Isd €hvs þ ðky2 l2  ky1 l1 Þv s þ ðky1 l1 þ ky2 l2 Þhvs  kz1 l1 wf ¼ f bL4  f bR4 2 2 w Isd €hw s þ ðkz1 l1  kz2 l2 Þws þ ðkz1 l1 þ kz2 l2 Þhs þ ky1 l1 v f ¼ f bL5  f bR5

in which Isd the moment of inertia of the spindle. For the lumped mass mf , its equation of axial motion could also be derived based on the force equilibrium condition

€ f þ ðkx1 þ kfx Þuf  kx1 us ¼ 0 mf u mf v€ f þ ðky1 þ kfy Þv f þ ky1 l1 hw s  ky1 v s ¼ mf g v

€ f þ ðkz1 þ kfz Þwf  kz1 l1 hs  kz1 ws ¼ 0 mf w

ð46Þ

202

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

3.4. Dynamic model of high-speed flywheels By assembling the equations of motion for the rotor and its supports, the nonlinear dynamic model for the high-speed flywheel could be obtained as follows

€ þ ðC1 þ XGÞq_ þ K1 q ¼ Fb þ Fe þ Fg Mq

ð47Þ

where q ¼ ½ufw ; v fw ; wfw ; nfw ; gfw ; us ; v s ; ws ; hvs ; hw s ; uf ; v f ; wf  the displacement vector, M; K1 ; C1 ; XG the mass, stiffness, damping and gyroscopic matrices of the system. The damping is assumed to be proportional, and one has C1 ¼ 1K1 , in which 1 the damping coefficient. The unbalanced excitations, nonlinear resorting forces of the rolling bearing and self weights of various components are reflected in Fe ; Fb ; Fg , respectively. Expressions for these matrices and vectors are given by T

M ¼ diagð½mfw ; mfw ; mfw ; Idfw ; Idfw ; ms ; ms ; ms ; Isd ; Isd ; mf ; mf ; mf Þ 2

Kfw 6 K1 ¼ 6 4 Ksfw 0

Kfws

0

ð48Þ

3

7 Ksf 7 5 Kf

Ks Kfs

ð49Þ

G ¼ diagð½XGfw ; 0Þ

ð50Þ

3 0 7 6 6 U s cosðXt þ us Þ 7 7 6 7 6 6 U s sinðXt þ us Þ 7 7 6 7 6 6 U d cosðXt þ ud Þ 7 7 6 7 6 6 U d cosðXt þ ud Þ 7 7 6 7 6 0 7 6 7 6 7 6 Fe ¼ X2 6 0 7; 7 6 7 6 0 7 6 7 6 7 6 0 7 6 7 6 7 6 0 7 6 7 6 7 6 0 7 6 7 6 7 6 0 5 4 0 2

3 f bL1 þ f bR1 7 6 6 f bL2 þ f bR2 7 7 6 7 6 6 f bL3 þ f bR3 7 7 6 7 6 6 f bL4 þ f bR4 7 7 6 7 6 6 f bL5 þ f bR5 7 7 6 7 6 6 f bL1  f bR1 7 7 6 7 6 Fb ¼ 6 f bL2  f bR2 7; 7 6 7 6 6 f bL3  f bR3 7 7 6 6 f  f 7 6 bL4 bR4 7 7 6 6 f  f 7 6 bL5 bR5 7 7 6 7 6 0 7 6 7 6 7 6 0 5 4 0 2

3 0 7 6 6 mfw g 7 7 6 7 6 6 0 7 7 6 7 6 6 0 7 7 6 7 6 6 0 7 7 6 7 6 6 0 7 7 6 7 6 Fg ¼ 6 ms g 7 7 6 6 0 7 7 6 7 6 6 0 7 7 6 7 6 6 0 7 7 6 7 6 6 0 7 7 6 7 6 6 mf g 7 5 4 0 2

ð51Þ

The motions of the flywheel and spindle are coupled through the nonlinear restoring forces of the two bearings. Thus, some of the sub-matrices in the stiffness matrix K1 are set to be zero matrices, i.e. Kfw ¼ 0; Kfws ¼ KTsfw ¼ 0. The other submatrices are expressed as follows

2

kx1 6 6 0 6 6 0 Ks ¼ 6 6 6 6 0 4

0

0

ðky1 þ ky2 Þ

0

0

ðkz1 þ kz2 Þ

ðky2 l2  ky1 l1 Þ

0

0

ðkz1 l1  kz2 l2 Þ

0 2 Kfs ¼

KTsf

kx1

6 6 0 6 ¼6 6 0 6 4 0 0

0 ky1 0 0 ky1 l1

0

0

þ ky2 l2 Þ

7 7 7 7 ðkz1 l1  kz2 l2 Þ 7 7 7 7 0 5

0

ðkz1 l1 þ kz2 l2 Þ

ðky2 l2  ky1 l1 Þ 0 2 ðky1 l1

3

0

2

0

2

ð52Þ

2

3

7 7 7 kz1 7 7 7 kz1 l1 5 0 0

Kf ¼ diagð½ðkx1 þ kfx Þ; ðky1 þ kfy Þ; ðkz1 þ kfz ÞÞ The sub-matrix of the gyroscopic matrix could be expressed as

ð53Þ

ð54Þ

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

2

Gfw

0 60 6 6 ¼6 60 6 40

0 0 0 0

0 0

0 0

0

0 0

0

0 0 0 Ipfw

3 0 0 7 7 7 0 7 7 7 Ipfw 5 0

203

ð55Þ

Eq. (45) is a series of two order differential equations, which represent the dynamic model of the high-speed flywheel system. Due to the nonlinear restoring forces and torques of the rolling bearing, these equations are nonlinearly coupled. Micro-vibration responses nfw ; gfw of the flywheel could be gained by numerically solving Eq. (45). Then, substituting the values of nfw ; gfw into Eqs. (8)–(10), the time-variable output torques of the SGCMG with various sources of micro-vibration disturbances are determined accordingly. Detailed solution process is shown in Fig. 11. For given preloads, the pre-deflections doj and contact angles boj of outer race of the two bearings under given preloads are gained. Through solving nonlinear algebraic equations, the restoring forces of left and right bearings could be obtained. Then, the ordinary differential equation solver in MATLAB (i.e. ODE45) is utilized to solve Eq. (45) to obtain the micro-vibration responses nfw ; gfw of the flywheel at current time. From Eqs. (8)–(10), output tor-

Fig. 11. Solution process for the output torques of SGCMG.

204

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

que values of the SGCMG at current time step are acquired accordingly. At every time interval, the nonlinear iteration analysis is conducted to obtained the bearing restoring forces more preciously. The calculation will not stop until the time exceeded the end time. Time history of output torques is finally obtained for post-processing analysis. 4. Model verification by dynamic tests 4.1. Experimental apparatus For a typical SGCMG, dynamic tests are carried out to verify the output torque model proposed by this study. Fig. 12 shows a schematic diagram of the output torque test system. The tested SGCMG is supported by the 3-axis air bearing table and placed on the Kistler-9253B dynamometer platform. The signal testing system mainly includes a charge amplifier, an OR35 data acquisition instrument and a personal computer. In dynamic tests, the gimbal coordinate is defined to be the test coordinate. The Z g axis is the rotation axis, which is also perpendicular to the direction of the dynamometer. The X g ; Y g axes are perpendicular to the gimbal rotational axis, and also the horizontal width and length directions of the dynamometer, as shown in Fig. 12. In the test, the SGCMG flywheel rotates around the gimbal axis. The rotational vector direction of the flywheel rotor is initially coincided with X g axis. The dynamometer will measure the real time dynamic torques around the three axes of the gimbal coordinate when the SGCMG is in working state. The physical pictures of tested SGCMG and output torque test system are, respectively, shown in Fig. 13(a) and (b). 4.2. Parameters of tested SGCMG Table 1 gives the values of system parameters, including the flywheel, supporting system and unbalanced mass excitation. The parameters of bearing B7004C and lubricant oil are listed in Table 2. Due to the limitation of experimental conditions, the surface waviness of raceways and rolling balls could not be accurately measured. In this study, the waviness amplitudes are set to be 0.2 lm, and the order is given in the range of 1  N b , where N b the number of balls. The phase between various waviness is assumed to be zero. According to Eqs. (21)–(26), the characteristic frequencies of waviness related to the inner and outer raceways are different. For the waviness of inner race, the characteristic frequencies IF n are expressed as follows

IF n ¼ nxc ¼

  n db cos b0 X n ¼ 1; 2;    1þ 2 de

ð56Þ

For the outer race, the characteristic frequencies of waviness are

OF m ¼ mðX  xc Þ ¼

  m db cos b0 1 X m ¼ 1; 2;    2 de

Fig. 12. Schematic diagram for dynamic output torque test of a SGCMG.

ð57Þ

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

205

Fig. 13. Physical pictures of (a) tested SGCMG and (b) output torque test system.

Table 1 Physical parameters of the tested SGCMG. Description

Value

Flywheel mass (kg) Diametral moment of inertia (kg m2 ) Polar moment of inertia (kg m2 ) Static mass eccentricity (kg m) Dynamic mass eccentricity (kg m2 ) Proportional damping coefficient Rotating speed (rpm) Bearing distances (m) Mass of spindle (kg) Moment of inertia of spindle (kg m2 ) Transversal stiffness of spindle (N/m) Axial stiffness of spindle (N/m) Mass of mounting brackets (kg) Axial support stiffness of mounting brackets (N/m) Transversal support stiffness of mounting brackets (N/m)

5.2 1.99e2 3.98e2 1e6 5e8 2e5 6000 20e3/30e3 1.5 2e-2 1.2e7, 1.2e7/1e7, 1e7 1e8 1.5 1.8e9 1e8, 1e8

Table 2 Parameters of the bearing B7004C and lubricant oil. Description Number of balls Nominal contact angle Pitch diameter (mm) Diameter of inner raceway (mm) Diameter of outer raceway (mm) Curvature radius of inner race (mm) Curvature radius of outer race (mm) Ball diameter (mm) Elastic modulus (N=m2 ) Poisson’s ratio Ball density (kg=m3 ) Lubricating oil viscosity Pa s Viscosity-pressure coefficient Pa1 Preload N

Value 12 150 31 20 42 3.1 3.1 6 2.075e11 0.3 7800 1.2e3 2.2e8 40

4.3. Comparisons and verification In the experiment, time domain waveforms of three axis output torques are, respectively, measured for the gimbal speed equal to 20 =s and 10 =s. The proposed theoretical model is also utilized to numerically obtain the time domain waveforms

206

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

of output torques. Comparisons between numerical and theoretical results are shown in Figs. 14 and 15. When the gimbal speed is 20 =s, the rotation period of the gimbal is 18s. From the time domain waveform of output torques, except for the torque T out;z around the Z g axis, the other two axis torques change with the gimbal rotation period, and there are 5.5 cycles in the 100s test time, as shown in Fig. 14. This is mainly because the output torques are modulated by the rotation of gimbal, which is consistent with the theoretical expectation (see Eqs. (13)–(15)). For the gimbal speed of 10 =s, as shown in Fig. 15, the amplitude of the time domain waveform is reduced to half and the period is increased to 36 s compared with that of the gimbal speed 20 =s. Compared with the experimental results, it can be seen that the time domain waveform of the theoretical model is basically consistent with the experimental waveform. Good agreements are found in the period and amplitude of the waveform. Both experimental and numerical results all show that besides the fluctuation with the gimbal period, there are highfrequency disturbances in the output torque waveform, which is mainly caused by the micro-vibrations of the high-speed flywheel. Fast Fourier Transformation (FFT) spectra of dynamic output torques T out;x ; T out;y for gimbal speed equal to 20 =s and 10 =s are, respectively, given in Figs. 16–19. For gimbal speed of 20 =s, besides the gimbal rotational speed xg , the flywheel rotational speed and its multiples (X; 2X; 3X) appear in the FFT spectra. This is caused by the static and dynamic mass imbalances of flywheel rotor. Bearing inner waviness frequencies (IF 1 ; IF 2 ; IF 3 ; IF 4 ; IF 5 ; IF 6 ) and outer waviness frequencies (OF 2 ; OF 3 ; OF 4 ; OF 5 ; OF 6 ; OF 7 ; OF 8 ; OF 9 ; OF 10 ) also have distinct amplitudes, indicating that the rolling bearing waviness has significant effect on the output torques of SGCMG. In addition, one can also find combination frequencies between flywheel rotation frequency and outer waviness frequency, i.e. X þ OF 1 ; X þ OF 3 ; X þ OF 4 ; X þ OF 6 ; X þ OF 7 , which shows the nonlinear coupling effect between the rotor imbalance vibration and waviness induced vibration. Comparing the results of T out;x and T out;y , the difference in terms of frequency components and amplitudes is not obvious. For gimbal speed of 10 =s, as shown in Figs. 18, 19, except for the gimbal rotation frequency xg , there is no obvious change in the frequency values and amplitudes of other high-frequency components comparing with the case of the gimbal speed equal to 20 =s. This verifies the hypothesis in theoretical derivation, which is the weak coupling assumption between the flywheel and gimbal vibrations. This is true as the difference between flywheel rotor speed and gimbal speed is too distinct. From the numerical predicted FFT spectra of the output torques (see Figs. 16(b), 17(b), 18(b) and 19(b)), one can find that numerical results also reflect the modulation effect on the output torque due to the micro-vibrations of the flywheel rotor. Various frequencies, including the flywheel rotational frequency and its multiples (X; 2X; 3X), bearing inner waviness frequencies (IF 1 ; IF 2 ; IF 3 ; IF 4 ; IF 5 ; IF 6 ), bearing outer waviness frequencies (OF 2 ; OF 3 ; OF 4 ; OF 5 ; OF 6 , OF 7 ; OF 8 ; OF 9 ; OF 10 ) and combination frequencies (X þ OF 1 ; X þ OF 3 ; X þ OF 4 ; X þ OF 6 ; X þ OF 7 ), which have been identified in the experimental results, could also be found in the numerical results. The numerical predicted frequency values are in good agreement with tested

Fig. 14. Waveforms of three axis output torques for gimbal speed equal to 20 =s: (a) experimental results and (b) numerical results.

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

Fig. 15. Waveforms of three axis output torques for gimbal speed equal to 10 =s: (a) experimental results and (b) numerical results.

Fig. 16. FFT spectra of T out;x for gimbal speed equal to 20 =s: (a) experimental results and (b) numerical results.

207

208

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

Fig. 17. FFT spectra of T out;y for gimbal speed equal to 20 =s: (a) experimental results and (b) numerical results.

Fig. 18. FFT spectra of T out;x for gimbal speed equal to 10 =s: (a) experimental results and (b) numerical results.

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

Fig. 19. FFT spectra of T out;y for gimbal speed equal to 10 =s: (a) experimental results and (b) numerical results.

Fig. 20. Effect of static (a) and dynamic (b) eccentricities on the value of DT.

209

210

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

frequency values, and maximum relative error is less than 5%. This verifies the accuracy of the theoretical model. However, the spectral amplitudes predicted by the theoretical model have distinct differences with that of tested results, especially for the spectral amplitudes of both inner/outer waviness frequencies and the combination frequencies. This is due to the limitation of the experimental conditions, which cause the measurement of actual shapes and amplitudes of rolling bearing inner/outer waviness difficult and impossible. In addition, the theoretical results also indicate that changing the gimbal speed has little effect on the high-frequency components of the output torque, which is also consistent with the experimental results. 5. Discussions Based on the SGCMG output torque model, various parameters, including rotor mass imbalances, bearing waviness, bearing axial preload and installation skew angles, are discussed for their effects on the high-frequency disturbances in the output torque. From Eqs. (13)–(15), one can see that the ideal output torques of SGCMG should be harmonic forms. However, due to the micro-vibration of flywheel and installation errors, there are always high-frequency disturbances in the output torques, as shown in Eqs. (8)–(10). The output torques along the X g ; Y g axes are similar with 90 degrees of phase difference. The root mean square (RMS) value of the high-frequency disturbances along the X g (or Y g ) axis could be defined as follows

DT ¼ RMSðT out;x  Tb out;x Þ

ð58Þ

Obviously, the DT could be utilized to measure the amplitude of the high-frequency disturbances in the output torque. Fig. 20 gives the variation of DT with static and dynamic mass imbalance excitations. Increasing the amount of static (or dynamic) eccentricity all increases the value of DT, indicating that the rotor mass imbalances have dominant influence on the

Fig. 21. Effect of bearing outer race waviness on the value of DT.

Fig. 22. Effect of bearing axial preload on the value of DT.

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

211

Fig. 23. Effect of installation skew angles abf ; bbf on the value of DT.

high-frequency disturbances in the output torque. In real applications, the amount of both static and dynamic imbalances should be strictly controlled to achieving the desired level of payload performance. The effect of bearing surface waviness on the value of DT is shown in Fig. 21. The amplitude of outer race waviness is increased from 0.2 lm to 0.4 lm. Increasing the waviness amplitude greatly enlarges the high-frequency disturbance amplitude DT. The effect of inner race waviness on the value of DT is similar with that of the outer waviness, and the result is excluded from the paper. Thus, the waviness amplitude should be restricted in a certain range to control the highfrequency disturbance torques transmitted to the spacecraft. Fig. 22 also discusses the effect of bearing axial preload on the value of DT. One can see that the bearing preload greater than 40 N almost has no effect on the high-frequency disturbances of the output torque. This means that the increase of bearing preload can guarantee the normal work of the bearing, but it does not play a significant role in reducing the high-frequency components of the output torque. From Eqs. (8)–(10), one can find that the installation skew angle vbf is not reflected in the expressions, and thus has no effect on the output torque. Here, the values of abf ; bbf are all increased from 1e-4 rad to 4e-4 rad, and the variations of DT is given in Fig. 23. With the increasing of abf ; bbf , the value of DT is also increased significantly, which means that the installation skew angles have great influence on the output torque disturbances. In the assembling of SGCMG, the skew angles should be strictly controlled to reduce the high-frequency disturbances in the output torque. 6. Conclusions The output torque modeling of a SGCMG with nonlinear rolling bearing supports is conducted in this paper. Taking the installation errors and micro-vibrations of the flywheel into account, three axis output torques of a SGCMG are derived based on Newton-Euler approach and theorem of moment of momentum. Dynamic model is then constructed to obtain the microvibration responses of the rotary flywheel. Mass imbalances of the flywheel, flexibility of supporting structures and nonlinearity induced by one pair of angular contact ball bearings are considered in the dynamic model. Especially for the rolling bearing, an improved load distribution analysis is proposed to more accurately obtain the contact deformations and angles between the rolling balls and raceways. Various factors, including the preload condition, surface waviness, Hertz contact and elastohydrodynamic lubrication, are included in the analysis. The bearing restoring forces are then obtained through iteratively solving the load distribution equations at every time step. Dynamic tests on a typical SGCMG supported by angular contact ball bearings are conducted to verify the output torque model. The effects of flywheel dynamic/static eccentricities, inner/outer waviness amplitudes, bearing axial preload and installation skew angles on the dynamic output torques are discussed. The obtained results would be useful for the optimal design and vibration control of the SGCMG system. Acknowledgments The research work described in the paper was supported by the National Science Foundation of China under Grant No. 11472147/51405015, and the State Key Laboratory of Tribology under Grant No. SKLT2015B12. References [1] S.P. Bhat, P.K. Tiwari, Controllability of spacecraft attitude using control moment gyroscopes, IEEE Trans. Autom. Control 54 (3) (2009) 585–590. [2] R. Votel, D. Sinclair, Comparison of control moment gyros and reaction wheels for small Earth-observing satellites, in: 26th Annual AIAA/USU Conference on Small Satellites, 2012, SSC12-X-1.

212

H. Wang et al. / Mechanical Systems and Signal Processing 115 (2019) 188–212

[3] V. Lappas, P. Oosthuizen, P. Madle, L. Cowie, G. Yuksel, D. Fertin, Micro CMGs for agile small satellites: design and in-orbit tests, in: Proceeding of the 6th International ESA Conference on Guidance, Navigation and Control Systems, Loutraki, Greece, 2006, pp. 1–9. [4] J. Gersh, M. Peck, Violet: a high-agility nanosatellite for demonstrating small control-moment gyroscope prototypes and steering laws, in: AIAA Atmospheric Flight Mechanics Conference, Chicago, IL, 2009. [5] B. Girouart, I. Sebbag, J.M. Lachiver, Performance of the Pleiades-HR agile attitude control system, in: Proceedings of the 5th ESA International Conference on Spacecraft Guidance Navigation and Control Systems, Frascati, Italy, 2002, pp. 497–500. [6] R.G. Cobb, J.M. Sullivan, A. Das, L.P. Davis, T.T. Hyde, T. Davis, Z.H. Rahman, J.T. Spanos, Vibration isolation and suppression system for precision payloads in space, Smart Mater. Struct. 8 (1999) 798–812. [7] C. Liu, X. Jing, S. Daley, F. Li, Recent advances in micro-vibration isolation, Mech. Syst. Signal Process. 56–57 (2015) 55–80. [8] C.E. Eyerman, A system engineering approach to disturbance minimization for spacecraft utilizing controlled structures technology (Master Degree thesis), MIT, 1990. [9] R.A. Masterson, D.W. Mille, R.L. Grogan, Development and validation of reaction wheel disturbance models: empirical model, J. Sound Vib. 249 (3) (2002) 575–598. [10] K.C. Liu, P. Maghami, C. Blaurock, Reaction wheel disturbance modeling, jitter analysis, and validation tests for solar dynamics observatory, in: AIAA Guidance, Navigation and Control Conference and Exhibit, 2008, pp. 1–18. [11] S.S. Narayan, P.S. Nair, A. Ghosal, Dynamic interaction of rotating momentum wheels with spacecraft element, J. Sound Vib. 315 (4–5) (2008) 970–984. [12] Z. Zhang, G.S. Aglietti, W. Zhou, Micro-vibrations induced by a cantilevered wheel assembly with a soft-suspension system, AIAA J. 49 (5) (2011) 1067– 1079. [13] W.Y. Zhou, G.S. Aglietti, Z. Zhang, Modelling and testing of a soft suspension design for a reaction/momentum wheel assembly, J. Sound Vib. 330 (18– 19) (2011) 4596–4610. [14] W.Y. Zhou, D.X. Li, Design and analysis of an intelligent vibration isolation platform for reaction/momentum wheel assemblies, J. Sound Vib. 331 (13) (2012) 2984–3005. [15] Z. Wei, D. Li, Q. Luo, J. Jiang, Modeling and analysis of a flywheel microvibration isolation system for spacecrafts, Adv. Space Res. 55 (2) (2015) 761–777. [16] Z. Zhang, G.S. Aglietti, W. Ren, Coupled microvibration analysis of a reaction wheel assembly including gyroscopic effects in its accelerance, J. Sound Vib. 332 (22) (2013) 5748–5765. [17] D. Addari, G.S. Aglietti, M. Remedia, Dynamic mass of a reaction wheel including gyroscopic effects: an experimental approach, AIAA J. 55 (1) (2016) 274–285. [18] D. Addari, G.S. Aglietti, M. Remedia, Experimental and numerical investigation of coupled microvibration dynamics for satellite reaction wheels, J. Sound Vib. 386 (2017) 225–241. [19] Q. Luo, D. Li, J. Jiang, Coupled dynamic analysis of a single gimbal control moment gyro cluster integrated with an isolation system, J. Sound Vib. 333 (2014) 345–363. [20] Q. Luo, D. Li, J. Jiang, Analysis and optimization of micro-vibration isolation for multiple flywheel systems of spacecraft, AIAA J. 54 (5) (2016) 1719– 1731. [21] S. Shi, D. Li, Q. Luo, Design and dynamic analysis of micro-vibration isolator for Single Gimbal Control Moment Gyro, Proc. Eng. 99 (2015) 551–559. [22] Y. Zhang, J. Zhang, Disturbance characteristics analysis of CMG due to imbalances and installation errors, IEEE Trans. Aerosp. Electron. Syst. 50 (2) (2014) 1017–1026. [23] Y. Zhang, Y. Zang, M. Li, Y. Wang, W. Li, Active-passive integrated vibration control for control moment gyros and its application to satellites, J. Sound Vib. 394 (2017) 1–14. [24] W. Zhou, D. Li, Q. Luo, K. Liu, Analysis and testing of microvibrations produced by momentum wheel assemblies, Chin. J. Aeronaut. 25 (4) (2012) 640– 649. [25] Q. Luo, D. Li, W. Zhou, J. Jiang, G. Yang, X. Wei, Dynamic modelling and observation of micro-vibrations generated by a Single Gimbal Control Moment Gyro, J. Sound Vib. 332 (2013) 4496–4516. [26] A. Jones, A general theory for elastically constrained ball and radial roller bearings under arbitrary load and speed conditions, J. Basic Eng. 82 (2) (1960) 309–320. [27] M. Tiwari, K. Gupta, O. Prakash, Dynamic response of an unbalanced rotor supported on ball bearings, J. Sound Vib. 238 (5) (2000) 757–779. [28] M. Tiwari, K. Gupta, O. Prakash, Effect of radial internal clearance of a ball bearing on the dynamics of a balanced horizontal rotor, J. Sound Vib. 238 (5) (2000) 723–756. [29] E.M. Yhland, Waviness measurement an instrument for quality control in rolling bearing industry, Proc. Inst. Mech. Eng. Conf. Proc. 182 (1967) 438– 445. [30] E.M. Yhland, A linear theory of vibration caused by ball bearings with form errors operating at moderate speeds, Trans. ASME J. Tribol. 114 (1992) 348– 359. [31] G.H. Jang, S.W. Jeong, Nonlinear excitation model of ball bearing waviness in a rigid rotor supported by two or more ball bearings considering five degrees of freedom, Trans. ASME J. Tribol. 124 (2002) 82–90. [32] G.H. Jang, S.W. Jeong, Stability analysis of a rotating system due to the effect of ball bearing waviness, Trans. ASME J. Tribol. 125 (2003) 91–101. [33] G.H. Jang, S.W. Jeong, Analysis of a ball bearing with waviness considering the centrifugal force and gyroscopic moment of the ball, J. Tribol. 125 (3) (2003) 487–498. [34] G.H. Jang, S.W. Jeong, Vibration analysis of a rotating system due to the effect of ball bearing waviness, J. Sound Vib. 269 (2004) 709–726. [35] C.Q. Bai, Q.Y. Xu, Dynamic model of ball bearings with internal clearance and waviness, J. Sound Vib. 294 (2006) 23–48. [36] M. Cao, J. Xiao, A comprehensive dynamic model of double-row spherical roller bearingModel development and case studies on surface defects, preloads, and radial clearance, Mech. Syst. Signal Process. 22 (2008) 467–489. [37] C.K. Babu, N. Tandon, R.K. Pandey, Vibration modeling of a rigid rotor supported on the lubricated angular contact ball bearings considering six degrees of freedom and waviness on balls and races, Trans. ASME J. Vibr. Acoust. 134 (1) (2012) 011001–011006. [38] H. Cao, Y. Li, X. Chen, A new dynamic model of ball-bearing rotor systems based on rigid body element, J. Manuf. Sci. Eng. 138 (7) (2016) 071007. [39] L. Niu, H. Cao, Z. He, Y. Li, An investigation on the occurrence of stable cage whirl motions in ball bearings based on dynamic simulations, Tribol. Int. 103 (2016) 12–24. [40] X. Zhang, Q. Han, Z. Peng, A new nonlinear dynamic model of the rotor-bearing system considering preload and varying contact angle of the bearing, Commun. Nonlinear Sci. Numer. Simul. 22 (1) (2015) 821–841. [41] X. Zhang, Q. Han, Z. Peng, F. Chu, Stability analysis of a rotor-bearing system with time-varying bearing stiffness due to finite number of balls and unbalanced force, J. Sound Vib. 332 (25) (2013) 6768–6784. [42] X. Zhang, Q. Han, Z. Peng, F. Chu, A comprehensive dynamic model to investigate the stability problems of the rotor-bearing system due to multiple excitations, Mech. Syst. Signal Process. 70–71 (2016) 1171–1192. [43] H. Wang, Q. Han, D. Zhou, Nonlinear dynamic modeling of rotor system supported by angular contact ball bearings, Mech. Syst. Signal Process. 85 (2017) 16–40. [44] H. Wang, Q. Han, R. Luo, T. Qing, Dynamic modeling of moment wheel assemblies with nonlinear rolling bearing supports, J. Sound Vib. 406 (2017) 124–145. [45] H. Cao, L. Niu, S. Xi, X. Chen, Mechanical model development of rolling bearing-rotor systems: a review, Mech. Syst. Signal Process. 102 (2018) 37–58. [46] B.J. Hamrock, D. Dowson, Ball Bearing Lubrication, Wiley, New York, 1981. [47] T. Harris, M. Kotzalas, Rolling Bearing Analysis, Taylor & Francis, Boca Raton, 2007.