Dynamic modeling of moment wheel assemblies with nonlinear rolling bearing supports

Journal of Sound and Vibration 406 (2017) 124–145

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Dynamic modeling of moment wheel assemblies with nonlinear rolling bearing supports Hong Wang a, Qinkai Han b,n, Ruizhi Luo a, Tao Qing a a b

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

a r t i c l e in f o

abstract

Article history: Received 15 September 2016 Received in revised form 12 June 2017 Accepted 12 June 2017 Handling editor: L.G. Tham

Moment wheel assemblies (MWA) have been widely used in spacecraft attitude control and large angle slewing maneuvers over the years. Understanding and controlling vibration of MWAs is a crucial factor to achieving the desired level of payload performance. Dynamic modeling of a MWA with nonlinear rolling bearing supports is conducted. An improved load distribution analysis is proposed to more accurately obtain the contact deformations and angles between the rolling balls and raceways. Then, the bearing restoring forces are then obtained through iteratively solving the load distribution equations at every time step. The effects of preload condition, surface waviness, Hertz contact and elastohydrodynamic lubrication could all be reflected in the nonlinear bearing forces. Considering the mass imbalances of the flywheel, flexibility of supporting structures and rolling bearing nonlinearity, the dynamic model of a typical MWA is established based upon the energy theorem. Dynamic tests are conducted to verify the nonlinear dynamic model. The influences of flywheel mass eccentricity and inner/outer waviness amplitudes on the dynamic responses are discussed in detail. The obtained results would be useful for the design and vibration control of the MWA system. & 2017 Elsevier Ltd All rights reserved.

Keywords: Moment wheel assemblies Nonlinear dynamic model Rolling bearing Surface waviness

1. Introduction A typical moment wheel assemblies (MWA) consists of a rotating flywheel supported by rolling ball bearings and driven by an internal brushless DC motor. The MWA can provide momentum bias and control torques without the consumption of precious fuel and can store momentum induced by very low frequency or DC torques. Thus, it has been widely used in spacecraft attitude control [1] and large angle slewing maneuvers [2] over the years. Among all devices onboard a spacecraft, the MWA is considered to be one of the largest vibrational disturbance sources [3]. The vibration forces and moments emitted by the MWA can severely degrade the performance of precision payloads in space [2,4,5]. Therefore, understanding and controlling vibration of MWAs is a crucial factor to achieving the desired level of payload performance. Appropriately modeling and analyzing vibrations of the MWA have been the focus of most current investigations. By assuming that disturbances consist of discrete harmonics of the wheel speed with amplitudes proportional to the wheel speed squared, Masterson et al. [6] first presented an empirical model to predict the vibration behavior of a wheel assembly in the Hubble Space Telescope. Because the interactions between the harmonics and structural modes of the wheel were n

Corresponding author. E-mail address: [email protected] (Q. Han).

http://dx.doi.org/10.1016/j.jsv.2017.06.019 0022-460X/& 2017 Elsevier Ltd All rights reserved.

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

125

ignored, the empirical model might under-predict the vibration disturbance in some wheel speed ranges. Elias and Miller [7] later suggested a coupled disturbance model that includes the dynamic mass matrix of each components of the MWA. Experiments performed to show that the coupled model predicts the measured dynamics more accurately than the previous model. Shigemune and Yoshiaki [8] developed a new detector for the empirical model which enables detecting the lowfrequency vibrational disturbances. Besides the empirical model mentioned above, the analytical model, established through two order differential dynamic equations by considering the mass imbalances of the flywheel and structural flexibility, has received more and more attentions in current literature. Considering the fly wheel with axial and rocking modes and gyroscopic torques on the wheel, Liu et al. [9] first established a dynamic model of the MWA for analyzing the jitter induced from both tonal and broadband disturbances. Ground tests were also performed to verify the dynamic model. The dynamic coupling problem of a MWA with its support brackets was studied by Narayan et al. [10]. Through dynamic simulations and vibration tests, it was found that the effects of gyroscopic forces arising out of rotating wheels are significant and this aspect needs to be taken into account while designing the mounting brackets. For a cantilevered MWA supported by a soft-suspension system, Zhou and Zhang [11,12] 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 [13,14]. Similar dynamic model was also developed for parameter identification of MWAs [15]. Zhang et al. [16] analyzed the coupled micro-vibrations of a reaction wheel assembly including gyroscopic effects in its accelerant state. Very recently, Peng et al. [17] proposed a complete micro-vibration dynamical model of a magnetically suspended MWA, in which the rotor dynamics, magnetic bearing control system and micro-vibration sources were considered. The proposed dynamic model was proved by comparisons with simulation and experimental results, and has been utilized to design an adaptive and robust resonant controller to eliminate the synchronous micro-vibration fluctuations [18]. A hybrid vibration isolator comprising passive and active components was proposed by Lee et al. [19], to provide an effective solution to the vibration problems caused by the MWAs. Dynamic tests on a flexible test-bed showed effective vibration isolation by the proposed vibration isolator. Addari et al. [20,21] carried out thorough experimental studies on a reaction wheel including the gyroscopic effects. Numerical simulations were also presented in order for comparisons and validations [21]. Most papers mentioned above focused on the modeling of mass imbalances of the flywheel and flexibility of supporting structures and their effects upon the vibrations of MWAs [9–21]. Being the core supporting component, the rolling element bearing plays a decisive role in the performance, operation reliability and service life of MWAs. In current dynamic models, the rolling element bearing is usually simplified by a linear spring-damper system [11–14]. Such treatment is obviously too ideal, i.e. not enough realistic. 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 MWAs. Zhou et al. [22] first introduced the bearing irregularity and nonlinear stiffness into the dynamic model of a wellbalanced MWA. 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. In their model, both the nonlinear stiffness and bearing irregularity induced by the surface waviness were given by empirical equations, which might reduce the model accuracy. It is well known that the nonlinear bearing supporting forces should be solved iteratively through load distribution analysis. Jones [23] first analyzed the modeling of an elastically constrained ball and radial roller bearings, and 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. [24,25] greatly improved Jones model by taking into account both radial clearance and Hertzian contact characteristics. Yhland [26,27] 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 [26,27], Jang and Jeong [28–31] 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 [32] proposed an improved load distribution model to study the effect of bearing waviness on the fluctuation of the cage speed. Cao and Xiao [33] 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. [34] extended the Jones model to six degrees of freedom, and found that the effect of frictional moment on the vibrations of the system was significantly dependent on the bearing load. Antoine et al. [35] pointed out that the contact forces and contact deformations of the bearing were directly dependent on the bearing contact angle through analyzing the relationship between preload, speed and contact angle with a new analytical approach. Recently, Zhang et al. [36,37] 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. From the above reviews, one can find 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, dynamic modeling of a MWA with nonlinear rolling bearing supports is conducted in this paper. The mass

126

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

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

imbalances of the flywheel, flexibility of supporting structures and nonlinearity induced by one pair of angular contact ball bearings are considered in the 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 upon a typical MWA supported by two angular contact ball bearings are conducted to verify the model. The effects of flywheel eccentricity and inner/outer waviness amplitudes on the dynamic responses are discussed in detail. Finally, some conclusions are given.

2. Dynamic modeling of the MWA A typical MWA, shown in Fig. 1, consists of a flywheel, a spindle, rolling element bearings, a DC motor, the housing and mounting brackets. The DC motor is used to drive the rotating flywheel under constant angular speed Ω . The flywheel is supported by one pair of high precision angular contact ball bearing with face-to-face installation. The spindle is stationary and non-rotating. It is fixed to the mounting brackets at one end and the other end is supported by the housing. The inner race of the bearing is fixed to the stationary spindle, and the outer race is connected to the flywheel and rotates. The axial preload for the bearing is applied through the inner and outer load sleeves, as shown in Fig. 1. Considering the mass unbalance of the flywheel, nonlinear supports from the rolling bearing and structural flexibility of the MWA, a lumped mass-spring-damper system with nonlinear stiffness couplings is established and shown in Fig. 2. The damper is not plotted in the figure, and the damping effect will be included in the dynamic model of the whole MWA system. Modeling considerations are described as follows: (1) The flywheel has complex spoke type structures. By ignoring the structural vibration, the flywheel is simplified as a lumped mass with five degrees of freedom, including three translations and two rotations. The coordinate system for the flywheel is denoted by Xfw − Yfw − Zfw , as shown in Fig. 2. (2) In operation, the unbalance excitation of the flywheel is transferred through the rolling element bearings to the spindle, the mounting brackets and then to the spacecraft. With the transmission of unbalance excitation, internal excitations from the rolling bearing element bearing would generate due to the microdisturbances, including Hertzian contacts, surface waviness, elastohydrodynamic lubrication, and so on. These internal excitations also greatly influence the

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

127

Fig. 2. Vibration model for the MWA.

vibrations of MWA. In this study, the internal excitations are taken into account by deriving the five directions of nonlinear restoring forces of the angular contact ball bearing. For the lower and upper bearings, the restoring forces are, respectively, denoted by fLbi , fUbi (i = 1, 2, 3, 4, 5). (3) Two lumped masses (ms , mf ) and a series of springs with linear stiffness are adopted in the model to consider the structural flexibilities of the spindle, housing and mounting brackets. The stationary spindle is modeled by the lumped mass ms with five degrees of freedom in the coordinate Xfw − Yfw − Zfw . Its motions are coupled with the flywheel motions through nonlinear rolling bearing forces ( fLbi , fUbi ). Besides, the transversal supports of ms are realized by the linear springs k x1, k y1, kz1 and k x2, k y2. The lumped mass mf is utilized to simulate the mounting brackets, and three translational motions are considered. The support stiffnesses are denoted by kfx , kfy, kfz , and the couplings between the ms and mf are expressed by k x1, k y1, kz1. In the following section, the above mentioned sub-systems will be, respectively, studied to gain their dynamic equations of motion. The whole system model will then be obtained by assembling the sub-systems' equations. 2.1. The flywheel The structural vibrations of the flywheel is neglected, and the rotor is simplified as a rigid mass with five degrees of freedom. Fig. 3 shows the vibration model of the flywheel. The coordinate system is Xfw − Yfw − Zfw . The rotor has three u v translations along the axes (ufw , vfw , wfw ) and two rotations around the Xfw and Yfw axes (θfw ). According equilibrium , θfw condition of forces, one can derive the equations of motion for the flywheel as follows

mfw u¨ fw = fu1 + fLb1 + fUb1 mfw v¨fw = fu2 + fLb2 + fUb2 ¨ fw = fLb3 + fUb3 − mfw g mfw w u

v

v

u

̇ = Tu1 + f Idfwθ¨fw + Ipfw Ωθfw + fUb4 Lb4 ̇ = Tu2 + f + f Idfwθ¨fw − Ipfw Ωθfw Lb5 Ub5

(1)

where mfw , Idfw , Ipfw the mass, diametral and polar moments of inertia of the flywheel, fLbi , fUbi (i = 1, 2, 3, 4, 5) the nonlinear restoring forces and moments of the lower and upper bearings, which will be solved in the following section. The fu1 , fu2 , Tu1, Tu2 in Eq. (1) are excitation forces and moments induced by the unbalance mass, which is inevitable due

128

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

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

to the non-uniform material, processing and assembly of the flywheel. Both static and dynamic imbalances are considered, and shown in Fig. 3. For the static imbalance, time-variant forces fu1 and fu2 will be induced along the Xfw and Yfw axes. Their expressions are given by

fu1 = UsΩ2cos(Ωt + φs ) fu2 = UsΩ2sin(Ωt + φs )

(2)

in which Us = mmsrms the static imbalance mass, φs the initial phase. For the dynamic imbalance, time-variant moments Tu1, Tu2 will be induced around the Xfw and Yfw axes. One can given their expressions as follows

Tu1 = UdΩ2cos(Ωt + φd ) Tu2 = UdΩ2sin(Ωt + φd )

(3)

in which Ud = mmdrmdhf the dynamic imbalance mass, φd the initial phase of the dynamic imbalance mass, and hf the axial length between the imbalance masses. 2.2. Nonlinear bearing forces The inner race of the angular contact ball bearing is fixed on the stationary spindle, and the outer race is connected to the flywheel and rotates under constant speed Ω . 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 microdisturbances. These excitations will generate vibrations and then influence the dynamic behaviors of the MWA. 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,

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

129

Fig. 4. Schematic diagram for the bearing coordinates.

load distribution analysis is conducted by taking preload condition and various microdisturbances 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 moments of the bearings are computed utilizing the Hertz contact and oil lubricant theorems.

2.2.1. Load distribution analysis Fig. 4 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 the orbital angular speed of ball ωcj. Except for the torsional rotation 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, θ u, θ v]T . Accordingly, the preloads along the five directions are fp = [fu , fv , fw , mu , mv ]T , in which fu , fv , fw are three forces and mu , mv are two moments around the X , Y axes. Usually, apart from the longitudinal force ( fw > 0), the preloads along the other four directions are all be equal to zero, i.e. fu = fv = 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 coordinate is expressed by δj = [ζj, ηj , θj ]T . Under the assumption of small deflection, the δj could be obtained from the rigid deformation of the outer race mass center u through coordinate transformations

δj = Tju + Δwj

(4)

where Tj is the transformation matrix, and is expressed as

⎤ ⎡ cosϕ sinϕ 0 0 0 cj cj ⎥ ⎢ 0 1 R osinϕcj −R ocosϕcj ⎥ Tj = ⎢ 0 ⎥ ⎢ ⎢⎣ 0 0 0 sinϕcj −cosϕcj ⎥⎦ in which ϕcj = ϕ0 + ωcjt +

2π (j Nb

(5)

− 1) the location angle of the jth ball at time t, where ϕ0 and Nb are the initial angle and

number of rolling ball. Without loss of generality, one has ϕ0 = 0. The Ro = 0.5de − (rgo − 0.5db)cosβ0 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 β0 is the nominal contact angle. Besides the rigid displacement of the outer race, the surface waviness would also induce additional displacements Δwj , as shown in Eq. (4). 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. 5 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

130

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

Fig. 5. Schematic diagram for the bearing waviness.

pij =

∑ Aincos(ninωcjt + 2π (j − 1)/Nb + ψin)

(6)

qij =

∑ Bincos(ninωcjt + 2π (j − 1)/Nb + φin)

(7)

poj =

∑ A out cos(nout (ωo − ωcj)t + 2π (j − 1)/Nb + ψout )

(8)

qoj =

∑ Boutcos(nout (ωo − ωcj)t + 2π (j − 1)/Nb + φout )

(9)

wij =

∑ Cncos(nbaωsjt + ϑba)

(10)

woj =

∑ Cncos(nbaωsjt + pi + ϑba)

(11)

where Ain , Bin, Aout , Bout , Cn , nin , nout , nba and ψin, φin , ψout , φout , ϑba are, respectively, denote the amplitudes, harmonic orders and phases of the waviness, ωsj is the self-rotation angular speed. Considering the fixed inner race and the pure rolling state, 2⎞ ⎛ dbcosβ0 0.5d one could derive the expressions of ball's self rotation speed ωsj = d e ⎜ 1 − ⎟ωo, and orbital speed of the ball de b ⎝ ⎠

(

(

ωcj = 0.5 1 +

dbcosβ0 de

)

)ω . According to the geometric relations in Figs. 4 and 5, it is easy to obtain o

⎡ poj − pij ⎤ ⎥ ⎢ Δwj = ⎢ qoj − qij ⎥ ⎥ ⎢⎣ 0 ⎦

(12)

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. 6. The contact angles for the ball with both inner and outer races are the same and equal to the nominal angle β0. 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. 6, 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. (4), the displacements of the groove curvature center of inner race are expressed by (ζj, ηj ). The displacements of ball center are unknown parameters and denoted by (vjy, vjz ). Due to the load on inner race, the inner and outer contact angles become different. The outer contact angle reduces to βoj , while the inner contact angle increases to βij . 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. 6. 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 (rgi, rgo), one has

L ij = rgi + pij − hci − (0.5db + wij )

(13)

L oj = rgo + poj − hco − (0.5db + woj )

(14)

in which hci , hco are the oil film thicknesses for the rolling ball contacting with the inner and outer raceways, respectively.

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

131

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

According to the result of Hamrock and Dowson [38], the oil film thickness is calculated by the following empirical equation

hc = 2.69R xU 0.67G 0.53W −0.067(1 − 0.61e−0.73κ ) where U =

η0uent E ′ Rx

denotes the dimensionless speed, in which

(15)

η0 is the lubricant viscosity at the atmospheric pressure and

temperature of 20 °C, Rx = db/2(1 ± db/decosβ0) represents the equivalent radius along the ball rolling direction. For the outer race, the “þ” is selected. The “  ” is chosen for the inner race. The equivalent rotating speed is uent =

dridro , 2(dri + dro)

where the

dri , dro denote the diameters of inner and outer races. The G = E′cηp is the dimensionless elastic modulus, in which cηp the viscosity-pressure coefficient. The dimensionless load is written by W =

Qj E ′ Rx2

. Obviously, the oil film thickness hc is related to

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

tanβij =

tanβoj =

lij =

loj =

L ijsinβ0 + vjz L ijcosβ0 − vjy

(16)

L ojsinβ0 + ηj − vjz L ojcosβ0 − ζj + vjy

(17)

L ijsinβ0 + vjz sinβij

(18)

L ojsinβ0 + ηj − vjz sinβoj

(19)

Eqs. (16)–(19) 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. 7 gives the forces applied on the jth rolling ball. In the figure, Fcj and Mgj denote the centrifugal and gyroscopic forces due to the rotation and revolution of the ball. Obviously, one has Fcj = 0.5mbdeωcj2 and Mgj = Ibωsjωcjsinαj , where mb,Ib the mass and moment of inertia of the ball, and αj the angle between self-rotation axis and the z axis. For the pure rolling state, one has

(

αj = tan−1

sinβ0 cosβ0 + db / de

)

.

The contact forces between the ball and inner/outer races are denoted by Qij and Qoj. According to the Hertz contact theory, one has Q ij = χij Kiδij3/2 and Q oj = χoj Koδoj3/2 , in which Ki,Ko and δij, δoj the contact stiffness coefficients and contact deflections between the ball and inner/outer races, respectively. From the above geometric analysis, one has δij = lij − L ij and δoj = loj − Loj . When δij, δoj > 0, χij = 1, χoj = 1. Otherwise, χij = 0, χoj = 0 for δij, δoj ≤ 0. The contact stiffness coefficients are

132

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

Fig. 7. Forces applied on the jth rolling ball [39].

determined by geometric dimensions of contact zone and material properties, i.e. K = (πκE′/3ξ 3ϵR/ξ ), in which E′ = E/(1 − ν 2) denote the effective elastic modulus, and E the material elastic modulus, ν the poisson's ratio, R the equivalent curvature radius, κ the elliptical ratio, ξ, ϵ the first and second kinds of elliptical integration. From the monograph of Harris [39], the expressions for these parameters are given and the contact stiffness coefficients Ki,Ko 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. 7. λ ij , λoj denote the friction coefficients between the ball and inner/outer races. Here, the values of λ ij , λoj are set to be 1. From Fig. 7, the force equilibrium equations are obtained as follows

Q ijcosβij − Q ojcosβoj +

−Q ijsinβij + Q ojsinβoj +

λ ijMgj db

sinβij −

λ ijMgj db

λ ojMgj

cosβij −

db

sinβoj + Fcj = 0

λ ojMgj db

(20)

cosβoj = 0

(21)

Eqs. (20) and (21) 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 Qoj and friction force

λojMgj db

to the outer

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

⎤ ⎡ λ M ⎢ Q ojcosβoj + oj gj sinβoj ⎥ db ⎥ ⎡Q ⎤ ⎢ ⎥ ⎢ yj ⎥ ⎢ λ ojMgj ⎢ Q zj ⎥ = ⎢ −Q ojsinβoj + cosβoj ⎥ ⎥ db ⎥ ⎢ ⎢ ⎥ ⎢⎣ Q θxj ⎥⎦ ⎢ M λ oj gj ⎥ ⎢ r go ⎥⎦ ⎢⎣ db

(22)

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

⎡Q ⎤ yj ⎥ T⎢ ⎡ ⎤ fb = ∑ ⎣ Tj⎦ ⎢ Q zj ⎥ ⎥ ⎢ j=1 ⎢⎣ Q θxj ⎥⎦ Nb

(23)

where fb = [fbu , fbv , fbw , mbu , mbv ]T 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

f p − fb = 0

(24)

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

133

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

Eqs. (20), (21) and (24) are a set of nonlinear algebraic equations, which could be solved iteratively by the Newton-Raphson method. Obviously, the dimension of the nonlinear equations is 2Nb + 5. After the unknown parameters are gained, the contact deformations (δij, δoj ) and inner/outer contact angles ( βij , βoj ) for every rolling ball could be computed by Eqs. (4) and (16)– (19). Moreover, the displacement vector uo of outer race under given values of preload could be obtained accordingly.

2.2.2. Bearing forces and moments Fig. 8 shows a pair of angular contact ball bearings (lower and upper bearings). The Xfw − Yfw − Zfw denotes the flywheel coordinate, and the XLb − YLb − ZLb and XUb − YUb − ZUb are, respectively, the analysis coordinates for lower and upper bearings. The axial distances between the two bearings and the flywheel mass center are l1, l2 . From the above descriptions, the u v T displacement vector for the flywheel is [ufw , vfw , wfw , θfw , θfw ] . The stationary spindle is represented by a lumped masses ms , u and the displacement vector is expressed as [us , vs, ws, θs , θsv]T . 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 lower and upper bearings could be gained using coordinate transformation technique

⎛ ⎡ ufw ⎤ ⎡ u ⎤⎞ ⎡ uLb ⎤ s ⎟ ⎜⎢ ⎥ ⎢v ⎥ vfw ⎥ ⎢ vs ⎥⎟ ⎜ ⎢ Lb ⎢ ⎥ ⎥ ⎢ ⎢ wLb ⎥ = T ⎜⎜ ⎢ wfw ⎥ − ⎢ ws ⎥⎟⎟ + u L ⎢ o ⎥ ⎢ u⎥ ⎢ u⎥ u ⎜ ⎢ θfw ⎥ ⎢ θs ⎥⎟ ⎢ θLb ⎥ ⎟ ⎜ ⎢⎣ θ v ⎥⎦ ⎜ ⎢⎢ θ v ⎥⎥ ⎢⎣ θsv ⎥⎦⎟ Lb ⎠ ⎝ ⎣ fw ⎦

(25)

134

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

⎛ ⎡ ufw ⎤ ⎡ u ⎤⎞ ⎡ uUb ⎤ ⎜⎢ ⎥ ⎢ s ⎥⎟ ⎢v ⎥ ⎜ ⎢ vfw ⎥ ⎢ vs ⎥⎟ ⎢ Ub ⎥ ⎢ wUb ⎥ = T ⎜⎜ ⎢ wfw ⎥ − ⎢ ws ⎥⎟⎟ + u U ⎢ o ⎥ ⎢ u⎥ ⎢ u ⎥ u ⎜ ⎢ θfw ⎥ ⎢ θs ⎥⎟ ⎢ θUb ⎥ ⎟ ⎜⎢ ⎢⎣ θ v ⎥⎦ ⎜ ⎢ θ v ⎥⎥ ⎢⎣ θsv ⎥⎦⎟ Ub ⎠ ⎝ ⎣ fw ⎦

(26)

u v T u v T where [uLb , vLb, wLb, θLb , θLb ] and [uUb , vUb, wUb, θUb , θUb ] , respectively, denote the displacement vectors of outer race of lower and upper bearings. According to the definition in Fig. 8, the transformation matrices TL, TU are derived as follows

⎡ 1 0 ⎢ 1 ⎢ 0 TL = ⎢ 0 0 ⎢ ⎢ 0 −1/l1 ⎢ ⎣ 1/l1 0

0 0 −l1⎤ ⎥ 0 l1 0 ⎥ ⎥ 1 0 0 ⎥, 0 1 0⎥ ⎥ 0 0 1⎦

⎡ 1 0 ⎢ ⎢ 0 −1 TU = ⎢ 0 0 ⎢ ⎢ 0 1/l2 ⎢ ⎣ 1/l2 0

0 0 −1 0 0

0 l2 ⎤ ⎥ l2 0 ⎥ ⎥ 0 0⎥ 1 0⎥ ⎥ 0 −1⎦

(27)

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 lower bearing is taken as an example to show the solution process for nonlinear bearing restoring forces in detail. The displacements of jth ball of lower bearing δLj could be gained by substituting the outer race displacement vector of u v T lower bearing uL = [uLb , vLb, wLb, θLb , θLb ] into Eq. (4). Similarly, the force equilibrium equations for every rolling balls of lower bearing are obtained, as shown in Eqs. (20) and (21). 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. (16)–(19) and (22), the reaction force vector from the jth ball to the outer race of lower bearing is [Q Lyj, Q Lzj, Q Lθxj ]T . Thus, in the system coordinate, the sum of reaction forces from all of the rolling ball, just the restoring forces of lower bearing, are expressed as follows

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

fLb1⎤ ⎥ fLb2 ⎥ ⎥ T fLb3 ⎥ = ⎡⎣ TL⎤⎦ ⎥ fLb4 ⎥ fLb5 ⎥⎦

⎤ ⎥ ∑ ⎡⎣ Tj⎤⎦ ⎢ Q Lzj ⎥ ⎥ ⎢ j=1 ⎢⎣ Q Lθxj ⎥⎦ Nb

⎡Q

T⎢

Lyj

(28)

Similarly, the restoring force vector of upper bearing could also be gained as follows

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

fUb1⎤ ⎥ fUb2 ⎥ ⎥ T fUb3 ⎥ = ⎡⎣ TU⎤⎦ ⎥ fUb4 ⎥ fUb5 ⎥⎦

⎡Q ⎤ Uyj ⎥ T⎢ ⎡ ⎤ ∑ ⎣ Tj⎦ ⎢ Q Uzj ⎥ ⎥ ⎢ j=1 ⎢⎣ Q Uθxj ⎥⎦ Nb

(29)

2.3. Structural coupling and flexibility The stationary spindle, housing and mounting brackets are used to support the rotating flywheel, 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. 2. The stationary spindle is simplified by a lumped mass ms with five degrees of freedom, including three translations us , vs, ws and two rotations θsu, θsv . In the coordinate Xfw − Yfw − Zfw , the ms is nonlinearly coupled with the flywheel through nonlinear bearing restoring forces and moments. Besides, the two ends of stationary spindle is supported by k x1, k y1, kz1 and k x2, k y2, 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 k x1, k y1, kz1. According to the force equilibrium condition, the equations of motion of the mass ms could be derived as follows

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

135

ms u¨ s + (k x1 + k x2)us + (k x2l2 − k x1l1)θsv − ksxzws − k x1uf = − fLb1 − fUb1 ms v¨s + (k y1 + k y2)vs + (k y1l1 − k y2l2)θsu − ksyzws − k y1vf = − fLb2 − fUb2 ¨ s + kszws − k z1wf − ksxzus − ksyzvs = − fLb3 − fUb3 − ms g ms w u Isdθ¨s + (k y1l1 − k y2l2)vs + (k y1l12 + k y2l22)θsu − k y1l1vf = − fLb4 − fUb4 v Isdθ¨s + (k x2l2 − k x1l1)us + (k x1l12 + k x2l22)θsv + k x1l1uf = − fLb5 − fUb5

(30)

in which ksxz , ksyz the coupling stiffness between the transversal and axial motions, 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

mf u¨ f + (k x1 + kfx)uf + k x1l1θsv − k x1us = 0 mf v¨f + (k y1 + kfy)vf − k y1l1θsu − k y1vs = 0 ¨ f + (k z1 + kfz )wf − k z1ws = − mf g mf w

(31)

2.4. Whole system model By assembling the equations of motion for the flywheel and its supports, the nonlinear dynamic model for the MWA could be obtained as follows

Mq¨ + (C1 + ΩG)q̇ + K1q = Fb + Fe + Fg

(32)

u v where q = [ufw , vfw , wfw , θfw , θfw , us , vs, ws, θsu, θsv, uf , vf , wf ]T the displacement vector, M , K1, C1, ΩG the mass, stiffness, damping and gyroscopic matrices of the system. The damping is assumed to be proportional, and one has C1 = ς K1, in which ς 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. According to the above derivations, the detailed expressions for these matrices and vectors are given by

M = diag ⎡⎣ mfw , mfw , mfw , Idfw, Idfw, ms , ms , ms , Isd, Isd, mf , mf , mf ⎤⎦

(33)

⎡K K fws 0 ⎤ ⎥ ⎢ fw K1 = ⎢ K sfw K s K sf ⎥ ⎥ ⎢ K fs K f ⎥⎦ ⎢⎣ 0

(34)

G = diag ⎡⎣ ΩGfw, 0⎤⎦

(35)

⎤ ⎡ f +f Ub1 ⎥ ⎢ Lb1 ⎡ U cos(Ωt + φ ) ⎤ ⎡ 0 ⎤ s s ⎥ ⎢ fLb2 + fUb2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ Ussin(Ωt + φs ) ⎥ ⎢ 0 ⎥ ⎢ fLb3 + fUb3 ⎥ ⎥ ⎢ ⎢ −mfw g ⎥ 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ fLb4 + fUb4 ⎥ ⎢ Udcos(Ωt + φd )⎥ ⎢ 0 ⎥ ⎥ ⎢ f +f ⎥ ⎢ ⎢ 0 ⎥ Ub5 ⎥ ⎢ Lb5 ⎢ Udcos(Ωt + φd )⎥ ⎢ 0 ⎥ ⎢ −fLb1 − fUb1 ⎥ ⎥ ⎢ 0 ⎥, F = ⎢ 0 ⎥ ⎥, Fb = ⎢ Fe = Ω2⎢ ⎥ ⎢ −fLb2 − fUb2 ⎥ g ⎢ 0 ⎥ ⎢ ⎢ −ms g ⎥ ⎥ ⎢ ⎥ ⎢ 0 ⎥ ⎢ ⎢ −fLb3 − fUb3 ⎥ ⎥ ⎢ ⎢ 0 ⎥ 0 ⎥ ⎢ −f ⎥ ⎢ f − ⎢ 0 ⎥ Ub4 ⎥ ⎢ Lb4 ⎥ ⎢ 0 ⎢ 0 ⎥ ⎢ −f − f ⎥ ⎥ ⎢ Lb Ub 5 5 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ 0 ⎥ 0 ⎥ ⎢ 0 ⎥ ⎢ m g − ⎢⎣ ⎦ f ⎥ ⎥ ⎢ ⎥⎦ ⎢⎣ 0 0 ⎥⎦ ⎢⎣ 0

(36)

(

(

)

)

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

136

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

Fig. 9. Solution process for the nonlinear dynamic model of MWA.

⎡ (k + k ) 0 0 −ksxz (k x2l2 − k x1l1) ⎤ x1 x2 ⎥ ⎢ 0 0 (k y1 + k y2) −ksyz (k y1l1 − k y2l2) ⎥ ⎢ ⎥ ⎢ k k k 0 0 − − Ks = ⎢ ⎥ sxz syz sz ⎥ ⎢ 2 2 k l k l k l k l 0 0 0 ( − ) ( + ) y1 1 y2 2 y1 1 y2 2 ⎥ ⎢ ⎥ ⎢ 0 0 0 (k x1l12 + k x2l22)⎦ ⎣ (k x2l2 − k x1l1)

(37)

⎡ −k 0 0 ⎤ ⎥ ⎢ x1 0 ⎥ − k y1 ⎢ 0 ⎥ ⎢ K fs = K Tsf = ⎢ 0 0 −k z1⎥ ⎢ 0 −k l 0 ⎥ y1 1 ⎥ ⎢ ⎢⎣ k x1l1 0 0 ⎥⎦

(38)

K f = diag ⎡⎣ (k x1 + kfx), (k y1 + kfy), (k z1 + kfz )⎤⎦

(

)

(39)

The sub-matrix of the gyroscopic matrix could be expressed as

Gfw

⎡0 ⎢ ⎢0 = ⎢0 ⎢0 ⎢ ⎢⎣ 0

0 0 0 0

0 0 0 0

0 0 0 0

0 0 −Ipfw

⎤ ⎥ ⎥ ⎥ Ipfw ⎥ ⎥ 0 ⎥⎦ 0 0 0

(40)

Eq. (32) is a series of two order differential equations, which represents the dynamic model of the whole WMA system. Due to the nonlinear restoring forces of the rolling bearing, these equations are nonlinearly coupled and could be solved through numerical integration method. Detailed solution process is shown in Fig. 9. For given preloads, the pre-deflections δoj and contact angles βoj of outer race of the two bearings under given preloads are gained. Through solving nonlinear algebraic equations, the restoring forces of lower and upper bearings could be obtained. Then, the ordinary differential equation solver in MATLAB (i.e. ODE45) is utilized to compute vibrational responses at the current time. At every time interval, the nonlinear iteration analysis is conducted to obtained the bearing restoring forces more preciously. Thus, compared with the current models, the nonlinear dynamic model presented in this paper would need more computing time. The calculation

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

137

will not stop until the time exceeded the end time.

3. Natural characteristics The natural characteristics should be studied to deeply understand the vibration behaviors of the WMA system. The nonlinear restoring forces induced by the rolling bearings should be linearized for natural characteristic analysis. Five linear stiffnesses (three translational and two rotational stiffnesses) for the lower and upper bearings are denoted by kLb1, kLb2, kLb3, kLb4, kLb5 and kUb1, kUb2, kUb3, kUb4, kUb5, respectively. The coupling effects between various degrees of freedom are neglected for simplification. Here, the translational stiffness of the lower bearing kLb1 is taken as an example to show the solution process of the linear bearing stiffness. According to the stiffness definition, the kLb1 could be computed by

kLb1 =

∂fu ∂u

≈

(fu + Δfu ) − fu (u + Δu) − u

(41)

where u denote the outer race deflection along the XLb axis for given preload fu. When the preload is increased from fu to fu + Δfu (usually Δfu is a small amount), the corresponding outer race deflection changes from u to u + Δu. Thus, the solution process of kLb1 is as follows: firstly, under given preload fu, the outer race deflection u could be gained through iteratively solving the nonlinear algebraic equations, as shown in Eqs. (20) and (21). Then, the preload value is increased as fu + Δfu , and the corresponding outer race deflection u + Δu are computed through similar iteration process. Finally, the linear bearing stiffness kLb1 is obtained by Eq. (41). The other bearing stiffnesses could be gained similarly. As the lower and upper bearings are a pair of contact ball bearings with same geometry, material and preload, so their linear stiffnesses are equal to each other, i.e. kLb1 = kUb1, kLb2 = kUb2, kLb3 = kUb3, kLb4 = kUb4 and kLb5 = kUb5. In the following, for convenience, the kb1, kb2, kb3, kb4, kb5 are used to represent the linear bearing stiffnesses. After the linear bearing stiffness is gained, the linearized dynamic model of the WMA could be expressed by

Mq¨ + (C1 + ΩG)q̇ + K1q = Fe + Fg

(42)

in which q, M , C1, G and Fe, Fg are the same with that of Eq. (32). Due to the linearized bearing supports, the nonlinear restoring force vector Fb does not exist, and elements of the stiffness matrix K1 in Eq. (34) should be modified to reflect the linear bearing stiffness. Detailed expressions are given as follows

K fw

⎡ ⎛ k k ⎞ ⎢ ⎜⎜ 2kb1 + b25 + b25 ⎟⎟ 0 0 ⎢ l1 l2 ⎠ ⎝ ⎢ ⎛ ⎢ k k ⎞ ⎜⎜ 2kb2 + b24 + b24 ⎟⎟ 0 0 ⎢ l1 l2 ⎠ ⎝ ⎢ =⎢ 0 0 2kb3 ⎢ ⎛ ⎢ ⎛1 1 ⎞⎞ ⎜⎜ kb2(l1 − l2) + kb4⎜ − ⎟⎟⎟ 0 ⎢ 0 l1 ⎠⎠ ⎝ l2 ⎝ ⎢ ⎢ ⎛1 1 ⎞⎞ ⎢⎛ 0 0 ⎢ ⎜⎜ kb1(l2 − l1) + kb5⎜⎝ l − l ⎟⎠⎟⎟ 1 2 ⎠ ⎣⎝ 0

⎛ ⎛1 1 ⎞⎞ ⎜⎜ kb2(l1 − l2) + kb4⎜ − ⎟⎟⎟ l1 ⎠⎠ ⎝ l2 ⎝ 0 (2kb4 + kb2l12 + kb2l22) 0

K fws = K Tsfw = − K fw

⎤ ⎛ ⎛1 1 ⎞⎞ ⎜⎜ kb1(l2 − l1) + kb5⎜ − ⎟⎟⎟⎥ l2 ⎠⎠⎥ ⎝ l1 ⎝ ⎥ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎥ 0 ⎥ ⎥ (2kb5 + kb1l12 + kb1l22) ⎦

(43)

(44)

138

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

Fig. 10. Schematic diagram for the dynamic experiment.

K s = K fw

⎡ (k + k ) 0 0 −ksxz (k x2l2 − k x1l1) ⎤ x1 x2 ⎥ ⎢ 0 0 (k y1 + k y2) −ksyz (k y1l1 − k y2l2) ⎥ ⎢ ⎥ ⎢ k z1 0 0 −ksxz −ksyz +⎢ ⎥ ⎥ ⎢ 2 2 k l k l k l k l 0 0 0 ( − ) ( + ) y1 1 y2 2 y1 1 y2 2 ⎥ ⎢ ⎥ ⎢ 0 0 0 (k x1l12 + k x2l22)⎦ ⎣ (k x2l2 − k x1l1)

(45)

The natural characteristics include the natural frequencies and corresponding mode shapes, which could be obtained by solving the eigenvalue problem of Eq. (42). Without the damping effect, the eigenvalue problem is expressed as follows

⎡ ω2M + ωΩG + K ⎤d = 0 ⎣ 1⎦

(46)

in which ω represents the eigenvalue of the MWA system and d denotes the corresponding mode shape vector. Eq. (46) shows a polynomial eigenvalue problem, which could be solved by utilizing the POLYEIG function in MATLAB. Because the existence of gyroscopic effect, the natural angular frequency would change with the flywheel speed Ω . Due to the Hertzian contacts, surface waviness and elastohydrodynamic lubrication, the bearing stiffness is nonlinearly varying with displacements. The nonlinearities of the rolling bearings mainly affect the Campbell diagram in two aspects: (1) The Hertzian contact leads to the stiffness increase with displacements, especially for the greater preload condition. This would cause the whirling frequencies in Campbell diagram increase accordingly. (2) The surface waviness in both inner and outer races induces extra displacement excitations to the system. Besides the rotating frequency excitations, the surface waviness will cause waviness frequencies excitations in the Campbell diagram. This could be found in both dynamic tests and numerical simulations of the MWA system in the following section. 2

4. Experiments and model validation Dynamic experiments are carried out to verify the nonlinear dynamic model of MWA. A high-speed flywheel is driven by a brushless DC motor. The rotating speed varies in the range of 0–5000 rpm. The test diagram is shown in Fig. 10. The tested MWA is fixed on a 3-D dynamic force measuring platform (9253B). Dynamic force and torque signals are amplified and

Fig. 11. Physical picture of the whole experimental system.

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

139

Table 1 Physical parameters of the tested MWA system. 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-transversal coupling stiffness of spindle (N/m) Axial stiffness of spindle (N/m) Mass of mounting brackets (kg) Axial stiffness of mounting brackets (N/m) Transversal stiffness of mounting brackets (N/m)

5.2 3.5e − 2 7e − 2 5e − 6 4e − 6 2e − 5 0–5000 30e − 3/30e − 3 1.5 2e − 3 1e7,1e7/6e6,6e6 1e6,1e6 1e8 1.5 1e9 1e8,1e8

Table 2 Parameters of the bearing B7004C and lubricant oil. Description

Value

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 Pa  1

12 150 31 20 42 3.1 3.1 6 2.075e11 0.3 7800 1.2e − 3 2.2e − 8

acquired by data acquisition system, and then a high performance computer is utilized for post-analysis. Fig. 11 gives the physical picture of the whole experimental system. 4.1. System parameters 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.1 μm, and the order is given in the range of 1 − Nb , where Nb the number of balls. The phase between various waviness is zero. According to Eqs. (6)–(11), the characteristic frequencies of waviness related to the inner and outer raceways are different. For the waviness of inner race, the characteristic frequencies IFn are expressed as follows

IFn = nωc =

dbcosβ0 ⎞ n⎛ ⎜1 + ⎟Ω 2⎝ de ⎠

n = 1, 2, ⋯

(47)

For the outer race, the characteristic frequencies of waviness are

OFm = m(Ω − ωc ) =

dbcosβ0 ⎞ m⎛ ⎜1 − ⎟Ω 2⎝ de ⎠

m = 1, 2, ⋯

(48)

4.2. Model verification With the rotating speed increasing from 0 rpm to 5000 rpm, the waterfall plots for the tested force in the Yfw axis and moment around the Xfw axis are shown in Figs. 12 and 13 respectively. One can see that the dominant spectra are the

140

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

Fig. 12. Waterfall plot for the tested force along the Yfw axis.

Fig. 13. Waterfall plot for the tested moment around the Xfw axis.

rotating speed Ω and its two times value 2Ω in the range of 0-400 Hz. Several orders of natural frequency of the MWA system are also found in the figures. The first two natural frequencies, as marked by ω1, ω2 , varies with rotating speed. They are identified as the backward (ω1) and forward (ω2) frequencies. Existence of whirling frequencies reflects the typical characteristics of a rotordynamic system. The third and fourth orders of natural frequency ω3, ω4 , which are found around the 250 Hz, varies slightly with the rotating speed. They are also identified as a pair of whirling frequency. The fifth natural frequency ω5, which locates around the 350 Hz, almost has no change with varying of rotating speed. Besides the rotating frequency and the natural frequencies, some frequencies related to the surface waviness are also found in the figures. According to Eqs. (47)–(48), these frequencies are identified as IF1, IF2, IF3, IF4, IF5, IF6, IF7 (related to the inner race waviness) and OF3, OF8 (related to the outer race waviness). Obviously, high orders of harmonic for the waviness

Fig. 14. First five orders of natural frequency of the linear MWA system varying with rotating speed.

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

141

Fig. 15. Numerical bearing dynamic force along the Yfw axis varying with rotating speed.

frequency could still be identified in the figure. Besides, two combination frequencies, Ω + IF2 and Ω + OF3, also exist in the dynamic response spectra. This is due to the coupling between the unbalance and waviness excited vibrations through the nonlinear bearing supports. The values of various parameters, listed in Tables 1 and 2, are substituted into the dynamic model of MWA. Then, both natural characteristics and dynamic responses could be gained through numerical analysis. The first five orders of natural frequency of the linear MWA system varying with rotating speed is shown in Fig. 14. From the figure, one can find that there are two pairs of whirling frequency, one of which is the ω1, ω2 , and the other is the ω3, ω4 . When the rotating speed is zero, the ω1, ω2 almost have the same value (about 100 Hz). The values of ω3, ω4 are about 250 Hz. With the increasing of rotating speed, the first pair of whirling frequency varies distinctly, while the change of the second pair of whirling frequency is not obvious. The fifth natural frequency ω5 maintains constant and its value is about 350 Hz. For the rotating speed at 5000 rpm, the backward frequency ω1 has dropped to 50 Hz, and the forward frequency ω2 has risen to 200 Hz. The second pair of whirling frequency ω3, ω4 varies slightly around the 250 Hz. Comparing with the test results in Figs. 12 and 13 one can see that the theoretical model basically reflects the change of the first three orders of natural frequency of the system. Both dynamic forces and moments of the bearing supports are gained by numerically solving the nonlinear dynamic model of the MWA system. The force along the Yfw axis and the moment around the Xfw axis are chosen for comparisons with the tested ones. The numerical results under various rotating speeds are given in Figs. 15 and 16, respectively. Almost all of the response spectra in the test signal could also be found in the numerical results. The simulated variations of these spectra (including the rotating frequency, natural frequencies and waviness characteristic frequencies) with the rotating speed are basically consistent with that of the test results. As the working speed of MWA is about 5000 rpm, so the response spectra for the rotating speed at 5000 rpm are selected for detailed comparisons. The results are shown in Figs. 17 and 18. One can see that the frequency values obtained by dynamic tests and numerical simulations are in good agreement, while the spectral amplitudes have more obvious differences. This is mainly due to the inevitable errors in geometric parameters and bearing surface waviness between the theoretical and tested models. It is well known in dynamic analysis that the spectral amplitudes are more sensitive than the frequency values. Thus, the nonlinear dynamic model presented in this study is still believable. To illustrate the computational efficiency, a typical convergence process for the load distribution analysis is shown in Fig. 19. The PC hardware details are: CPU INTEL CORE I5-3337U 1.8 GHz, RAM 8GB and hard disk capacity 500G. Using this PC, the convergence time is 0.5770 s. For every time steps, two load distribution computations should be done as there are

Fig. 16. Numerical bearing dynamic moment signal around the Xfw axis varying with rotating speed.

142

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

Fig. 17. Comparison between experimental and numerical dynamic forces for rotating speed equal to 5000 rpm.

Fig. 18. Comparison between experimental and numerical dynamic moments for rotating speed equal to 5000 rpm.

Fig. 19. Convergence process for the load distribution analysis.

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

143

Fig. 20. Effect of static imbalance upon the dynamic force response for rotating speed equal to 5000 rpm.

Fig. 21. Effect of dynamic imbalance upon the dynamic moment response for rotating speed equal to 5000 rpm.

two rolling bearings (lower and upper bearings). If the interest time length is 1 s, then it is equally divided into 1e4 parts. For every time step (1e  4 s), the time for solving the differential equations (i.e. ODE45) under this PC is 0.7130 s. Thus, the total computational time for dynamic response of MWA system with 1 s is (0.5770*2 þ 0.7130)*1e4 ¼ 18,760 s.

5. Discussions 5.1. Effect of unbalanced mass excitations The rotating speed is set to be 5000 rpm. Figs. 20 and 21, respectively, give the variations of dynamic response with static and dynamic mass imbalance excitations. Increasing the amount of static imbalance would increase the force response amplitude. From the response spectra in Fig. 20(b), one can find that the amplitude increase is mainly concentrated at the rotating frequency Ω and its doubling 2Ω . The other spectral amplitudes almost remain unchanged. Similar phenomenon could also be found for the case of dynamic mass imbalance, as shown in Fig. 21. Thus, in real applications, the amount of both static and dynamic imbalances should be strictly controlled to achieving the desired level of payload performance. 5.2. Effect of surface waviness Based upon the dynamic model of MWA, the effect of bearing surface waviness on the dynamic response is analyzed. Here, the rotating speed is still 5000 rpm, and the amplitudes of both inner and outer race waviness are, respectively, increased from 0.1 μm to 0.2 μm. Variations of dynamic force response are shown in Figs. 22 and 23. For the inner waviness, one can see from Fig. 22(a) that increasing the waviness amplitude greatly enlarges the dynamic force response. The amplitude increases are concentrated not only at the inner waviness related frequencies (such as: IF1, IF2, IF3, ⋯), but also at the

144

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

Fig. 22. Effect of inner waviness amplitude upon the dynamic force response for rotating speed equal to 5000 rpm.

Fig. 23. Effect of outer waviness amplitude upon the dynamic moment response for rotating speed equal to 5000 rpm.

combined frequencies Ω + IF2, Ω + IF3, as shown in Fig. 22(b). However, it is surprising to find that the amplitude of rotating frequency Ω is reduced with the increasing of inner waviness amplitude. Compared with the inner waviness, the effect of outer waviness on the dynamic force response is relatively weaker. From Fig. 23(a), one can find that the force amplitude is increased after adding the outer waviness amplitude to 0.2 μm. The increase is not significant and mainly concentrated at the outer waviness related frequencies OF3, OF8 , as shown in Fig. 23(b). Thus, the waviness amplitude should be restricted in a certain range to control the force transmitted to the spacecraft.

6. Conclusion The dynamic modeling of a MWA with nonlinear rolling bearing supports is conducted in this paper. An improved load distribution analysis is proposed to more accurately obtain the contact deformations and angles between the rolling balls and raceways. Then, the bearing restoring forces are then obtained through iteratively solving the load distribution equations at every time step. The effects of preload condition, surface waviness, Hertz contact and elastohydrodynamic lubrication could all be reflected in the nonlinear bearing forces. Considering the mass imbalances of the flywheel, flexibility of supporting structures and rolling bearing nonlinearity, the dynamic model of a typical MWA is established based upon the energy theorem. Dynamic tests upon a typical MWA supported by two angular contact ball bearings are conducted to verify the model. The influences of flywheel mass eccentricity and inner/outer waviness amplitudes on the dynamic responses are discussed in detail. The obtained results would be useful for the design and vibration control of the MWA system.

H. Wang et al. / Journal of Sound and Vibration 406 (2017) 124–145

145

Acknowledgments The research work described in the paper was supported by the National Natural Science Foundation of China under Grant No. 51405015/11472147, and the State Key Laboratory of Tribology under Grant No. SKLT2015B12.

References [1] T. Marshall, T. Gunderman, F. Mobley, Reaction wheel control of the MSX satellite, in: Proceedings of the Annual Rocky Mountain Guidance and Control Conference, AAS paper 91-038, 1991, pp. 119–138. [2] L.P. Davis, J.F. Wilson, R.E. Jewell, J.J. Roden, Hubble Space Telescope Reaction Wheel Assembly Vibration Isolation System, NASA Marshall Space Flight Center, 1986, pp. 1–37. [3] R.A. Laskin, M.S. Martin, Control/structure system design of a spaceborne optical interferometer, in: Proceedings of the AAS/AIAA Astrodynamics Specialist Conference, AAS Paper, 1989, pp. 89–424. [4] K.J. Pendergast, C.J. Schauwecker, Use of a passive reaction wheel jitter isolation system to meet the advanced X-ray astrophysics facility imaging performance requirements, Proc. SPIE, 3356, 1998, pp. 1078–1094. [5] M.D. Hasha, Reaction Wheel Mechanical Noise Variations, Space Telescope Program, Engineering Memo No. SSS 218, June 1998. [6] R.A. Masterson, D.W. Miller, R.L. Grogan, Development and validation of reaction wheel disturbance models: empirical model, J. Sound Vib. 249 (3) (2002) 575–598. [7] L.M. Elias, D.W. Miller, A coupled disturbances analysis method using dynamic mass measurement techniques, in: Proceedings of the AIAA/ASME/ ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2002, pp. 2002–1252. [8] T. Shigemune, O. Yoshiaki, Experimental and numerical analysis of reaction wheel disturbances, JSME Int. J. 46 (2) (2003) 519–526. [9] K.C. Liu, P. Maghami, C. Blaurock, Reaction wheel disturbance modeling, jitter analysis, and validation tests for solar dynamics observatory, in: Proceedings of the AIAA Guidance, Navigation and Control Conference and Exhibit, 2008, pp. 1–18. [10] 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. [11] Z. Zhang, G.S. Aglietti, W. Zhou, Microvibrations induced by a cantilevered wheel assembly with a soft-suspension system, AIAA J. 49 (5) (2011) 1067–1079. [12] 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. [13] 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. [14] 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. [15] Y. Zhao, J. Sun, H. Tian, Development of methods identifying parameters in reaction wheel assembly disturbance model, Aircr. Eng. Aerosp. Technol. 78 (4) (2006) 326–330. [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] C. Peng, J. Fang, P. Cui, Dynamics modeling and measurement of the microvibrations for a magnetically suspended flywheel, IEEE Trans. Instrum. Meas. 64 (12) (2015) 3239–3252. [18] C. Peng, Y. Fan, Z. Huang, B. Han, J. Fang, Frequency-varying synchronous micro-vibration suppression for a MSFW with application of small-gain theorem, Mech. Syst. Signal Process. 82 (2017) 432–447. [19] D.O. Lee, G. Park, J.H. Han, Hybrid isolation of micro vibrations induced by reaction wheels, J. Sound Vib. 363 (2016) 1–17. [20] 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. [21] 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. [22] W. Zhou, D. Li, Q. Lu, K. Liu, Analysis and testing of microvibrations produced by momentum wheel assemblies, Chin. J. Aeronaut. 25 (4) (2012) 640–649. [23] 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. [24] M. Tiwari, K. Gupta, O. Prakash, Dynamic response of an unbalanced rotor supported on ball bearings, J. Sound Vib. 238 (5) (2000) 757–779. [25] 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. [26] E.M. Yhland, Waviness measurement an instrument for quality control in rolling bearing industry, Proc. Inst. Mech. Eng. Conf. Proc., 182, 1967, pp. 438– 445. [27] 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. [28] 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. [29] 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. [30] 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. [31] 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. [32] C.Q. Bai, Q.Y. Xu, Dynamic model of ball bearings with internal clearance and waviness, J. Sound Vib. 294 (2006) 23–48. [33] M. Cao, J. Xiao, A comprehensive dynamic model of double-row spherical roller bearing model development and case studies on surface defects, preloads, and radial clearance, Mech. Syst. Signal Process. 22 (2008) 467–489. [34] 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. Vib. Acoust. 134 (1) (2012) 011001–011006. [35] J.F. Antoine, G. Abba, A. Molinari, A new proposal for explicit angle calculation in angular contact ball bearing, J. Mech. Des. 128 (2) (2006) 468–478. [36] 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. [37] 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. [38] B.J. Hamrock, D. Dowson, Ball Bearing Lubrication, Wiley, New York, 1981. [39] T. Harris, M. Kotzalas, Rolling Bearing Analysis, Taylor & Francis, Boca Raton, 2007.