Nonlinear dynamics and stability of wind turbine planetary gear sets under gravity effects

European Journal of Mechanics A/Solids 47 (2014) 45e57

Contents lists available at ScienceDirect

European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol

Nonlinear dynamics and stability of wind turbine planetary gear sets under gravity effects Yi Guo a, *, Jonathan Keller a, Robert G. Parker b a b

National Wind Technology Center, National Renewable Energy Laboratory, Mail Stop: 3811, 15013 Denver West Parkway, Golden, CO 80401-3305, USA Department of Mechanical Engineering, Virgina Tech, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 23 August 2012 Accepted 20 February 2014 Available online 12 March 2014

This paper investigates the dynamics of wind turbine planetary gear sets under the effect of gravity using a modified harmonic balance method that includes simultaneous excitations. This modified method along with arc-length continuation and Floquet theory is applied to a lumped-parameter planetary gear model including gravity, fluctuating mesh stiffness, bearing clearance, and nonlinear tooth contact to obtain the dynamic response of the system. The calculated dynamic responses compare well with time domain-integrated mathematical models and experimental results. Gravity is a fundamental vibration source in wind turbine planetary gear sets and plays an important role in the system dynamic response compared to excitations from tooth meshing alone. Gravity causes nonlinear effects induced by tooth wedging and bearing-raceway contacts. Tooth wedging, also known as a tight mesh, occurs when a gear tooth comes into contact on the drive-side and back-side simultaneously and it is a source of planetbearing failures. Clearance in carrier bearings decreases bearing stiffness and significantly reduces the lowest resonant frequencies of the translational modes. Gear tooth wedging can be prevented if the carrier-bearing clearance is less than the tooth backlash. Ó 2014 Elsevier Masson SAS. All rights reserved.

Keywords: Dynamics Wind turbine Planetary gear

1. Introduction The National Renewable Energy Laboratory (NREL) Gearbox Reliability Collaborative (GRC) was established by the U.S. Department of Energy in 2006. Its key goal is to understand the root causes of premature gearbox failures (Musial et al., 2007) through a combined approach of dynamometer testing, field testing, and modeling (Link et al., 2011), resulting in improved wind turbine gearbox reliability and a reduction in the cost of energy. As a part of the GRC program, this paper investigates gravity-induced dynamic behaviors of planetary gear sets in wind turbine drivetrains that could reduce gearbox life. Planetary gear sets have been used in wind turbines for decades because of their compact design and high efficiency. Despite these advantages, planetary gear sets generate considerable noise and vibration. Vibration causing high dynamic loads may result in gear tooth and bearing failures (Musial et al., 2007). Fatigue failures are a concern in long life-cycle applications. Analyzing the dynamics of wind turbine planetary gear

* Corresponding author. Tel.: þ1 303 384 7187; fax: þ1 303 384 6901. E-mail addresses: [email protected], [email protected] (Y. Guo). http://dx.doi.org/10.1016/j.euromechsol.2014.02.013 0997-7538/Ó 2014 Elsevier Masson SAS. All rights reserved.

drivetrains is important in improving gearbox life and reducing noise and vibration. The majority of wind turbines use a horizontal-axis configuration; thus, gravity becomes a periodic excitation source in the rotating carrier frame. Prior study of gravity by the authors was performed with a static analysis and focused on the effect of gravity upon bearing force and tooth wedging in a planetary spur gear set (Guo and Parker, 2010). It was found that tooth wedging, an abnormal contact situation where the tooth is in contact with both the drive-side and back-side flanks simultaneously, was caused by gravity. Tooth wedging increases planet-bearing forces and disturbs load sharing among the planets, which could lead to premature bearing failure. Significant in-plane translational gear component motions in planetary systems lead to tooth wedging. It is the combined effect of gravity and bearing clearance nonlinearity. Bearing clearance results in greater translational vibration, while gravity is the dominant excitation source causing the large motions that lead to tooth wedging. For heavy planetary gear sets, tooth wedging is likely to occur. Tooth wedging in planetary gear sets leads to unequal load sharing and excessive planet-bearing loads by disturbing the symmetry of the planet gears. This may cause bearing failure and tooth damage (Guo and Parker, 2010).

46

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57

Research on tooth separation is well-established for automotive and helicopter applications. Tooth separation was observed in spur gear pair experiments (Blankenship and Kahraman, 1996). Botman (1976) experimentally observed tooth separation in planetary gear sets. Using finite element and lumped-parameter models, Ambarisha and Parker (2007) predicted tooth separation and other nonlinear phenomena in a planetary gear set in a helicopter gearbox. Velex and Flamand (1996) investigated tooth separation at critical speeds. Bahk and Parker (2011) derived closed-form solutions for the dynamic response of planetary gear sets with tooth separation based on a purely torsional model. Nonlinear dynamics induced by bearing clearance has been studied for relatively small geared systems. Kahraman and Singh (1991) observed chaos in the dynamic response of a geared rotorbearing system with bearing clearance and backlash. Gurkan and Ozguven (2007) studied the effects of backlash and bearing clearance in a geared flexible rotor and the interactions between these two nonlinearities. Guo and Parker (2012a) investigated the nonlinear effects and instability caused by bearing clearance in helicopter planetary gear sets. Dynamic effects of bearing clearance in wind turbine planetary gear sets have not been studied in the past because the wind turbine operating speed was believed to be well below the frequency range of drivetrain dynamics. However, bearing clearance reduces some gearbox resonances significantly. The finite element program developed by Vijayakar (1991) uses a combined surface integral and finite element approach to capture tooth deformation and contact loads in geared systems. This finite element model includes bearing clearance, tooth separation, tooth wedging, fluctuating mesh stiffness, and gravity. Numerical integration is widely adopted to compute dynamic responses of mechanical systems in the time domain. Ambarisha and Parker (2007) used numerical integration to study nonlinear dynamics and the impacts of mesh phasing on vibration reduction of planetary gear sets. Velex and Flamand used numerical integration results of a planetary gear set with time-varying mesh stiffness as a benchmark to evaluate results from a Ritz method. The harmonic balance method (Thomsen, 2003) calculates nonlinear, frequency domain, steady-state response of mechanical systems. Zhu and Parker (2005) used this method to study clutch engagement loss in a belt-pulley system. Al-shyyab and Kahraman (2005a,b) investigated primary resonances, subharmonic resonances, and chaos in a multimesh gear train caused by fluctuating gear mesh stiffness. Bahk and Parker (2011) employed harmonic balance to analyze planetary gear dynamics based on a purely rotational model. Use of the harmonic balance method reduces computational time for lightly damped or physically unstable systems by avoiding the long transient decay time before a steady state is reached. Compared to numerical integration and finite element analysis, which are widely adopted approaches to compute dynamic responses, the computation time of the harmonic balance method is oneetwo orders of magnitude lower. Harmonic balance often employs arc-length continuation (Nayfeh and Balachandran, 1995) and Floquet theory (Raghothama and Narayanan, 1999; Seydel, 1994) to calculate nonlinear resonances in the dynamic response, including unstable solutions that numerical integration and finite element analysis are unable to obtain. The established harmonic balance formulation is only suitable for systems with one fundamental excitation frequency and its higher harmonics. However, wind turbine drivetrains have simultaneous internal and external excitations, including fluctuating mesh stiffness, gravity, bending-moment-induced excitations in the rotating carrier frame, wind shear, tower shadow, and other aero-induced excitations. The major objectives of this study are to: 1) develop a modified harmonic balance method to obtain the dynamic response of wind

turbine planetary gear sets considering gravity, fluctuating mesh stiffness, bearing clearance, and nonlinear tooth contact; 2) validate the proposed method by comparing the calculated results against experimental data and a numerical integration approach; and 3) investigate the gravity-induced dynamic behaviors using the developed approach, which includes tooth wedging, tooth contact loss, and bearing-raceway contacts. 2. Gearbox description This study investigates both a 750-kW wind turbine planetary gear (PG-A) used by the GRC (Link et al., 2011) and a 550-kW wind turbine planetary gear (PG-B) (Guo and Parker, 2010; Larsen et al., 2003; Rasmussen et al., 2004). These drivetrains have a main bearing that supports the main shaft and rotor weight, and two trunnion mounts that support the gearbox. These two gearboxes have similar configurations and are representative of the majority of three-point-mounted wind turbine drivetrains. PG-A is configured in a helical planetary gear arrangement with two parallel stages as shown in Fig. 1. Within the gearbox, there are two cylindrical roller bearings supporting the carrier, each with 275 mm of clearance. The planetary gear set is arranged in an inphase bridged carrier design with three equally spaced planets. The sun pinion shaft is connected to the intermediate stage through a spline joint that partially floats the sun pinion. The ring gear is bolted to the front and rear of the gearbox housing. The rated torque is 322,610 Nm and the rated speed of the input shaft is 22.2 rpm (Guo et al., 2012). PG-B is configured in a spur planetary gear arrangement with two parallel stages. Like PG-A, PG-B has three equally spaced planets. The rated torque is 180,000 Nm and the rated speed is 30 rpm. Additional key parameters of these two gearboxes are listed in Tables 1 and 2. A large ring gear mass in these designs is a result of a typical wind turbine gearbox arrangement whereby much of the gearbox housing is rigidly connected to the ring. 3. Mathematical models 3.1. Lumped-parameter model for planetary gear sets A previously developed and validated lumped-parameter model was adopted for this paper (Guo and Parker, 2010, 2012a). As depicted in Fig. 2(a), the carrier, ring, sun, and planets are rigid bodies, each having two translational and one rotational degree of freedom. The carrier rotating frame is used as the general coordinates for all of the components of planetary gear sets. This twodimensional model has 3(N þ 3) degrees of freedom, where N is the number of planets. The model includes gravity, fluctuating mesh stiffness, bearing clearance, tooth contact loss, and tooth wedging.

Fig. 1. The GRC gearbox (PG-A) configuration.

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57 Table 1 System parameters for the 750-kW Gearbox Reliability Collaborative (GRC) planetary gearbox (PG-A). The ring mass includes the gearbox housing and parallel stages. Sun Mass (kg) Moment of Inertia (kg-m2) Number of Teeth Pitch Diameter (mm) Root Diameter (mm) Average Mesh Stiffness (N/m) Bearing Stiffness (N/m) Carrier-Bearing Stiffness (N/m) Carrier-Bearing Clearance (mm) Torsional Support Stiffness (N/m)

Ring

Carrier

Planet

181.6 2633 759.9 104 3.2 144.2 59.1 3.2 21 99 e 39 215.6 1016.4 e 400.4 186.0 1047.7 e 372.9 9 9 ksp ¼ 16.9  10 , krp ¼ 19.2  10 100 102  106 5  109 6.8  109 3.2  109 0.275 45.8  106 57.4  106 0 0

Bearings are modeled using circumferentially distributed radial springs with a uniform clearance. This study focuses on the carrier bearing with clearance as shown in Fig. 2(b), which has dynamic effects on low-speed resonances (Guo and Parker, 2012a). The coordinates are shown in Fig. 2(a). Throughout this paper, the subscripts c, r, s, p denote the carrier, ring, sun, and planet; the subscripts B, g, m denote the bearing, gravity, and gear mesh; and superscripts b, d denote drive-side and back-side tooth contact. Translational displacements in the x and y directions, xw, yw, w ¼ c, r, s, are assigned to the carrier, ring, and sun, respectively, with regard to the rotating carrier frame. The origin O is at the center of the planetary gear set. The radial and tangential displacements of the jth planet are denoted by xj, hj, j ¼ 1, ., N with respect to the carrier and oriented for each planet as shown in Fig. 2(a). The rotational displacements are uv ¼ rvqv, v ¼ c, r, s, 1, ., N, where qv is the rotation in radians and rv is the base circle radius for the sun, ring, and planets and the radius to the planet center for the carrier. The masses and moments of inertias of the carrier, ring, sun, and planets are denoted by mk, Ik, k ¼ c, r, s, p. Quantities kwx, kwy, kwu denote the bearing stiffnesses of the carrier, ring, and sun supports in x, y, and u directions. The torsional stiffnesses of the carrier, ring, 2 , where k and sun supports equal kwu rw wu is the torsional stiffness with units of force/length and rw, w ¼ c, r, s is the base radius. The jth planet-bearing stiffness is kpj. The mesh stiffnesses at the jth sun-planet and ring-planet mesh are ksj, krj. The radial stiffnesses of the bearing that connects the carrier to ring gear in the x and y directions are kcrx, kcry. The nondimensionalized equations of motion of planetary gear sets are

~ 0 þ ~f d ðs; zÞ þ ~f b ðs; zÞ þ ~f ðs; zÞ ¼ Fð ~ sÞ ~ 00 þ Cz Mz B m m x ~ M ~ C ~d f dm ~b f bm ; M ¼ ; C ¼ ; fm ¼ ; fm ¼ ; 2 L M MLu MLu MLu2 ~f ¼ f B ; F ~ ¼ F B MLu2 MLu2

z ¼

(1)

where s ¼ ut and z ¼ x/L, where u is characteristic frequency, L is the characteristic length, and M is the characteristic mass. Derivations of the mass matrix M and the displacement vector x are detailed in Guo and Parker (2010). C ¼ ðU1 ÞT diagð2zn Un ÞU1 where zn are the damping ratios and Un are the natural frequencies of the linear system where all bearings are in contact and the mesh stiffnesses are averaged over a mesh cycle; U is the orthonormalized modal matrix (UTMU ¼ I). The nonlinear forces fall into three categories: the drive- and back-side tooth mesh force vectors f dm and f bm , and bearing force vector fB (Guo and Parker, 2010). The external force F(t) ¼ Fs þ Fg(t). Fs includes the static torques applied

47

to each component. Fg(t) is the gravity force vector. Gravitational force acting on the carrier, ring, sun, and planets is periodic in the rotating carrier frame, resulting in a fundamental external excitation source.

h iT y h x h x x x y x y Fg ðtÞ ¼ fcg ; fcg ; 0; frg ; frg ; 0; fsg ; fsg ; 0; f1g ; f1g ; 0; .; fNg ; fNg ; 0

(2)

h

x ; f ðw ¼ c; r; sÞ and f x ; f ðj ¼ 1; .; NÞ denote the Quantities fwg wg jg jg gravity force acting on the carrier, ring, sun, and planets 1eN. These are y

x ¼ m g sinðU tÞ fwg w c y ¼ mw g cosðUc tÞ; fwg

w ¼ c; r; s

  fjgx ¼ mj g sin Uc t þ jj   h fjg ¼ mj g cos Uc t þ jj ;

j ¼ 1; .; N

(3)

(4)

where the variable Uc ¼ um/Nr (if the ring is fixed) denotes the carrier rotation frequency. um denotes the mesh frequency. Nr denotes the number of teeth on the ring. The mesh stiffness fluctuates as the number of teeth in contact changes and is also an important internal excitation source of geared systems. These excitations are included through time-varying mesh stiffnesses. The mesh stiffnesses are calculated using the Calyx program (Nayfeh and Balachandran, 1995). Its mean amplitudes are listed in Tables 1 and 2 (Appendix). This program analyzes gear tooth contact and rolling element contact by using a combined analytical/finite element analysis detailed in Vijayakar (1991). Its results have been compared against studies of gear dynamics (Singh, 2010, 2011; Guo and Parker, 2012b). 3.2. Finite element model The two-dimensional finite element model includes 36 quadrilateral elements per tooth, four nodes per element, and each node has three degrees of freedom. It uses a combined surface integral and finite element method to capture tooth deformation and contact loads in geared systems (Kahraman and Blankenship, 1994). The software developed by Vijayakar (1991) intrinsically evaluates time-varying tooth contact forces that are specified externally with conventional simulation tools. Fluctuating mesh stiffness over a mesh cycle, tooth wedging, and tooth separation are included intrinsically in the finite element model. The finite element model also includes clearance nonlinearity at the carrier-ring bearing. Other bearings are modeled as linear stiffnesses without clearance in the finite element and analytical models. The finite element approach has been validated by experiments of geared systems (Parker et al., 2000a,b; Kahraman and Vijayakar, 2001). It is used to benchmark the established lumped-parameter model when experimental data is unavailable. 4. Extended harmonic balance method The extended harmonic balance method is used to obtain the dynamic responses of the model in Eq. (1). The formulation includes two excitation sources with excitation frequencies U1 and U2. The coupling effects between these two excitations are considered by including their side bands. Other excitation sources can be considered in a similar way. The response z is expanded into a Fourier series and assumed to include the R1, R2, R3R4R5, and R6R7R8 harmonics of excitation frequencies U1 and U2 and their

48

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57

Table 2 System parameters for the 550-kW planetary gearbox (PG-B). The ring mass includes the gearbox housing and parallel stages. Sun Mass (kg) Moment of Inertia (kg-m2) Number of Teeth Pitch Diameter (mm) Root Diameter (mm) Average Mesh Stiffness (N/m) Bearing Stiffness (N/m) Gearbox Trunnion Stiffness (N/m) Carrier-Bearing Stiffness (N/m) Carrier-Bearing Clearance (mm) Torsional Support Stiffness (N/m)

Ring

Carrier

Planet

51 4000 1330 114 61.1 2484 314.7 51.9 16 68 e 26 224 952 e 364 202 980 e 329 9 9 ksp ¼ 3.95  10 , krp ¼ 5.29  10 100 e 4  109 5.3  109 126  106 3.6  109 1 0 3  106 24.4  106 0

side bands mU1 þ lU2 and pU1  qU2. Each component zh in z then has a total of 2(R1 þ R2 þ R3R4R5 þ R6R7R8) þ 1 terms and is expressed as

2 6 G¼6 4

3

n

G1 þ G2 þ

R4 P R3 P m¼1 l¼1

G3 ðm;l;R5 Þþ

R7 P R6 P p¼1 q¼1

7

G4 ðp;q;R8 Þ 7 5 n (7)

The response derivatives are then transformed into

z00 ¼ GU21 A1 z  GU22 A2 z  G G

R7 X R6 X

R4 R3 P P m¼1 l¼1

ðmU1 þ lU2 Þ2 A3 z

ðpU1  qU2 Þ2 A4 z

p¼1 q¼1

zh ¼ zh;1 þ

R1  X

zh;2i cos iU1 s þ zh;2iþ1 sin iU1 s



z0 ¼ GU1 B1 z þ GU2 B2 z þ G

i¼1

i¼1 m¼1 l¼1

þ zh;2aði;m;lÞþ1 sin iðmU1 þ lU2 Þs þ

R7 X R6 X

þG

ðpU1  qU2 ÞB4 z

(8)

p¼1 q¼1

i¼1

R5 X R4 X R3 n X zh;2aði;m;lÞ cos iðmU1 þ lU2 Þs

ðmU1 þ lU2 ÞB3 z

m¼1 l¼1

R2 h i X þ zh;2ðiþL1 Þ cos jU2 s þ zh;2ðiþL1 Þþ1 sin iU2 s

þ

R4 X R3 X

(5)

o

where the operators A and B are also defined in Appendix A. The nonlinear force vectors, detailed in Guo and Parker (2012a), are transformed as

R8 X R7 X R6 n X zh;2bði;p;qÞ cos iðpU1  qU2 Þs i¼1 p¼1 q¼1

o þ zh;2bði;p;qÞþ1 sin iðpU1  qU2 Þs

~f d ¼ Gf d ; m m

where a(i, m, l) ¼ l þ L2 þ (m  1)R3 þ (i  1)R3R4 and b(i, p, q) ¼ q þ L3 þ (p  1)R6 þ (i  1)R6R7. L1 ¼ R1, L2 ¼ R1 þ R2,

2

z1;1 ; .; z1;2L1 þ1 ; .; z1;2L2 þ1 ; .; z1;2L3 þ1 ; .; z1;L4 ; .; |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

uN ðU2 Þ

uN ðU1 ;U2 Þ

L3 ¼ R1 þ R2 þ R3R4R5, L4 ¼ 2(R1 þ R2 þ R3R4R5 þ R6R7R8) þ 1. The time domain is discretized into n  1 evenly distributed time intervals [s1, ., sn]. Each component zh in the response z is extended into a vector in the time domain as ½zh ðs1 Þ; .; zh ðsn ÞT . The response vector transforms into z ¼ Gz by defining a function G that maps the response from the frequency domain to the time domain. Additional terms G1eG4 are defined in Appendix A.

~f ¼ Gf ; B B

~ ¼ GF F

(9)

By substituting Eqs. (8) and (9) into the equations of motion, Eq. (1) yields

3T

6 7 6 7 xc ðU1 Þ xc ðU2 Þ xc ðU1 ;U2 Þ 7 z ¼ 6 6 z3Nþ3;1 ; .; z3Nþ3;2L þ1 ; .; z3Nþ3;2L þ1 ; .; z3Nþ3;2L þ1 ; .; z3Nþ3;L 7 1 2 3 4 4 5 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} uN ðU1 Þ

~f b ¼ Gf b ; m m

(6)

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57

82 R4 X R3 < X 2 G 4  U21 MA1  U22 MA2  ðmU1 þ lU2 Þ MA3 :

and JB are described in Guo and Parker (2012a). The Jacobian matrix associated with back-side tooth loads, Jbm , is

2

m¼1 l¼1



R7 X R6 X

ðpU1  qU2 Þ MA4 þ U1 B1 þ U2 B2 2

p¼1 q¼1

þ

R4 X R3 X

ðmU1 þ lU2 ÞCB3

(10)

m¼1 l¼1

þ

R7 X R6 X

3

þ fB þ

d fm

þ

b fm

 F

K ¼  U1 MA1  U2 MA2  2

9 =

Lm;b r1;i

Lm;b r2;1

P

;

2

þ

¼ 0

R4 X R3 X

ðmU1 þ lU2 ÞCB3 þ

m¼1 l¼1

ðmU1 þ lU2 Þ MA3 2

2 P Lb

c1;i

þ Lc;B

6 6 6 6 6 vFB G¼6 H 6 6 vz 6 6 6 4 Symmetric

Lc;B Lbr þ Lc;B Lbs

Lbc2;1 . Lbc2;N 3 7 7 . 7 7 7 . 7 7 7 0 7 Lbp;1 . 7 7 5 1 Lbp;N

R7 X R6 X

(17) ðpU1  qU2 ÞCB4

p¼1 q¼1

þ KLTI (11) Eq. (10) then becomes b

Kz ¼ F  f B  f m  f m

(12)

By using these derivations, the ordinary differential equations of motion in Eq. (1) are transformed into the linear algebraic equations in Eq. (12). The global Newton method can be used to calculate solutions of Eq. (12). The Jacobian matrix J for the global Newton iteration is defined as

J ¼

  d b v Kz þ f B þ f m þ f m  F (13)

vz

Smoothing functions are employed to approximate the piecewise nonlinear forces so that the explicit formulation of J can be derived. These nonlinear forces include bearing reaction forces and tooth loads on the drive-side and back-side of the gears. The inverse Fourier transformation operator H is defined such that z ¼ Hz, which maps from the time domain to the frequency domain.

2 6 H¼6 4

3

n H1 þ H2 þ

R4 P R3 P m¼1 l¼1

H3 ðm; l; R5 Þ þ

R7 P R6 P p¼1 q¼1

7 7 5

H4 ðp; q; R8 Þ n

(14) where H1eH4 are the inverse Fourier transformation operators associated with G1eG4, respectively. The explicit formation of J is obtained as

J ¼ KþH

d vf m

/

The Jacobian matrix associated with carrier-bearing forces, JB, is

ðpU1  qU2 Þ2 MA4 þ U1 B1 þ U2 B2

d

3 7 Lm;b r2;n 7

7 m;b 7 Lm;b Lm;b s2;1 / Ls2;n 7 s1;i 7 7 7 Lm;b 7 p;1 7 7 5 1

p¼1 q¼1 R4 X R3 X

.

(16)

m¼1 l¼1



6 6 6 6 6 vf m G¼6 H 6 6 vz 6 6 6 4 Symmetric

P

Lm;b p;n

where KLTI is the time-invariant stiffness matrix defined in Appendix B. The function G is not singular, so it can be eliminated from both sides of Eq. (10). By defining

R7 X R6 X

0

ðpU1  qU2 ÞCB4 þ KLTI 5z

p¼1 q¼1



49

d Hðvf m =vzÞG

b vf m

vf B GþH GþH G vz vz vz

(15)

Individual components in Jbm are defined in Appendix B. J can be approximated numerically by employing finite difference algorithms (Turner and Walker, 1992). The numerical determination of J is time consuming in high dimensional space. Furthermore, its accuracy depends on the step size of the finite difference algorithms and can be disrupted by computational round-off errors. Using the line search technique (Grippo et al., 1986) for the Newton method improves the iteration convergence rate. 4.1. Arc-length continuation and stability analysis Arc-length continuation traces coexisting stable and unstable solutions of the dynamic response. The direction of the estimated next solution is selected to be perpendicular to the direction of the previous solution. The iteration convergence rate relies on the initial guess and the step size in the arc-length direction. Solution stability is determined using FloqueteLiapunov theory (Nayfeh and Balachandran, 1995). A perturbation Dz to the calculated periodic response z is introduced. The equations of motion are linearized for small Dz to give



Dz0 Dz00



¼

  0  JðtÞ  K

I C



Dz Dz0

 (18)

The monodromy matrix of the linear system (Eq. (18) is computed numerically by integrating the transformation matrices in multiple time intervals (Kahraman, 1992; Friedmann, 1990). Eigenvalues of the monodromy matrix are Floquet multipliers that characterize the solution stability and bifurcation. 5. Model validations Experiments on PG-A provide a benchmark for the lumpedparameter model using the extended harmonic balance method; however, the relevant available experimental data for PG-A is only at a single speed. Thus, results from the extended harmonic balance

50

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57

Fig. 2. (a) The lumped-parameter model and (b) side view of the planetary gear set.

method are also compared to the time-integrated results of the lumped-parameter model of PG-B within the speed range of 0e 200 Hz. 5.1. Comparison to experimental data The GRC project instrumented two identical 750-kW wind turbine gearboxes for the dynamometer (Fig. 3) and field testing. Internal measurements include gear tooth loads, main shaft torque and bending, internal component deflections and misalignments, and planet-bearing loads (Link et al., 2011). The dynamic response of PG-A is computed using the extended harmonic balance method at rated speed and rated torque. 5% modal damping is assumed for PG-A based on the modal damping extracted from a numerical torque impulse test of the PG-B finite element model. Torque applied to the lumped-parameter model is scaled to exclude the out-of-plane tooth loads caused by the planet gear helix angle. The translational displacement of the carrier over a carrier cycle computed using the extended harmonic balance method is compared against the time-synchronized experimental data measured in the GRC project as shown in Fig. 4. The agreement between the extended harmonic balance method and experimental data is reasonably good in capturing the maximum displacement amplitude. Differences between the analytical model and experiments are present in the frequency domain. The frequency spectrum of the carrier displacement (in Fig. 4) is shown in Fig. 5. A modal component at 6P (six times the carrier rotational frequency) is dominant in the frequency spectrum of the

measured signal. The passing frequency of the carrier planet pin bores could be the main contributor to this 6P component. The carrier is a bridged design and thus has three pin bores on both the upwind and downwind sides. These two circular sets of pin bores are misaligned because of manufacturing tolerances and torque windup in operation, which could be the cause of the 6P excitation. This 6P component might also be caused by the in-plane deformation of the carrier rim due to its varying thickness (Oyague et al., 2010). The magnitude of the 6P component is of the 1P magnitude for the calculated carrier displacement. The lumped-parameter model considers the carrier as a single lumped mass; thus, it does not consider the detailed carrier geometry. Consequently, the lumped-parameter method underestimates the magnitude of the displacement at 6P. The analytical results are dominated by 1P and 3P excitations due to the carrier rotation and planet symmetry. The experiments include shaft misalignments, manufacturing tolerances, assembly errors, which cause high vibration amplitudes compared to the analytical model that does not include these imperfections during the assembly and manufacturing. The experimental data includes more excitation sources than the analytical model, resulting in the richer frequency spectra than that of the model. In addition, the displacement of the first mesh frequency harmonic at 99P is more than an order of magnitude lower than the displacements at the low-speed excitations (1Pe6P) for both the calculated and measured results. Calculated planet-bearing forces are compared against the experimental data in Fig. 6. The mean and fluctuating amplitudes of the bearing forces computed by the extended harmonic balance

Fig. 3. A 750-kW wind turbine gearbox of the NREL GRC during a dynamometer test.

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57

Fig. 4. Carrier radial displacement with respect to the ring calculated using the extended harmonic balance method (red e) and measured by the GRC project (/). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

method largely match those of the experimental data. The planetbearing force for each individual planet bearing over a carrier cycle is asymmetric because of different pin position errors and bearing clearances for each planet. The lumped-parameter model predicts a 1.4% difference between the peak amplitudes of the planet-bearing forces over a carrier cycle because of these errors and clearances. The measured bearing force of planet 2 deviates from the other two planets, which is caused by a signal offset error. 5.2. Comparison to computational methods Dynamic responses computed by the extended harmonic balance method are compared against finite element and numerical integration analyses for PG-B. In lightly damped systems, finite element analysis and numerical integration require many time steps for the transient response to diminish so that steady-state data can be obtained. The extended harmonic balance method avoids these long duration integration simulations by calculating the steady-state dynamic response in the frequency domain. Numerical integration and finite element analysis compute only stable

51

Fig. 6. Measured (/) and calculated (red e) planet-bearing forces over a carrier cycle. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

responses. Conversely, the harmonic balance method uses arclength continuation to find unstable responses by tracing solutions as the speed changes quasi-statically. The large ring gear mass results from a typical arrangement of wind turbines whereby much of the gearbox mass is rigidly connected to the ring, which is supported on a relatively compliant foundation. The input torque is applied to the carrier, and the sun gear is the output. Gravity acts at the center of mass of all of the component gears. For the ring and sun, it is a periodically varying external excitation in the carrier reference frame in which the model is formulated. One carrier period is analyzed using the finite element model, giving 68 mesh cycles. Each mesh cycle is divided evenly into 10 intervals. The finite element solutions have single precision accuracy while the convergence tolerance is 100  1012 mm for analytical solutions. A global Newton iteration scheme with line search technique (Baker and Overman, 2000) is used to obtain accurate numerical solutions of the analytical model. The line search technique makes the predicted solution at each iteration always closer to the exact solution by adjusting the iteration step size for each step. For the numerical integration, 120 mesh cycles were simulated at each speed, of which 20 cycles are considered as transient time. Each mesh cycle is discretized into 50 intervals. For the harmonic balance method, R1 ¼ R2 ¼ 3, R3 ¼ 1, R4 ¼ 4. Each carrier period is discretized into 5291 equally spaced intervals. The solutions of the harmonic balance method have the convergence tolerance of 1  106 mm. The natural frequencies of PG-B below 400 Hz without bearing clearance and nonlinear tooth contact are listed in Table 3. Modal analysis is performed numerically on the finite element model of PG-B by applying torque and force impulses at the carrier and sun. Numerical impulse tests provide natural frequencies and their mode shapes as described in Table 3. The natural frequencies predicted by Table 3 Natural frequencies of PG-B.

Fig. 5. Frequency spectrum of the calculated (red e) and measured (/) carrier displacement. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Natural Frequency Lumped-parameter (Hz)

Natural Frequency Deviation Mode Type Finite Element (Hz) (%)

Damping (%)

40.75 69.49 157.53 301.46

40.94 69.1 163.7 312.1

5.48 5.26 7.43 5.19

0.35 0.56 3.77 3.41

Translational Rotational Translational Rotational

52

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57

Fig. 7. Root-mean-square (RMS, mean removed) dynamic response of PG-B for various mesh frequencies calculated by the extended harmonic balance method (red C) and numerical integration (/). Unstable solutions calculated by the harmonic balance method are denoted by (red B). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. (a) Dynamic carrier-bearing force of PG-B calculated using the harmonic balance (red e) and numerical integration methods (e), (b) percentage of bearing contact in a carrier cycle of PG-B with various speeds, given Dc ¼ 1 mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the lumped-parameter model match those obtained by the finite element approach. The modal damping ratios are extracted from the frequency response functions obtained from numerical impulse tests of the finite element model, with damping ratios obtained from the half-power points of the resonances listed in Table 3. Quasi-static results of the lumped-parameter model have been correlated with finite element analyses in Guo and Parker (2010) for tooth loads and planet-bearing forces. As shown in Fig. 7, translational displacements of the ring and sun are compared for the extended harmonic balance and numerical integration methods. The agreement between these two methods is very good up to 100 Hz. The nonlinear resonance at 40 Hz is caused by tooth wedging. Using speed sweeps, numerical integration predicts nonlinear discontinuities in the dynamic response. Harmonic balance calculates the unstable branches between these nonlinear jumps, denoted by the symbol (B). The agreement among the proposed method, numerical integration, and finite element analysis validates this method for predicting nonlinear dynamic responses of wind turbine planetary gear sets.

A unique nonlinear effect caused by fluctuating bearing contact is present in the dynamic response as shown in Fig. 8(a). Good agreement between the numerical integration and harmonic balance methods is evident. Two unstable branches are displayed in the crossing resonance: curves BC and DE. Increasing the speed, the response computed by numerical integration has the first discontinuity at the point B and the second at the point D. Decreasing the speed, the first discontinuity to higher amplitude occurs at the point E. Then the response has the first discontinuity to lower amplitude at the point C. The nonlinear effect displays multiple types of bifurcation points including saddle-node bifurcation points BeE and the transcritical bifurcation point A. The bearing contact in a carrier cycle under the same conditions of Fig. 8(a) is shown in Fig. 8(b). When the carrier speed increases from 16 to 20 Hz, the percentage of bearing contact in a carrier cycle decreases. The resonance bends to the low-frequency range at point B, a softening effect. When the carrier speed decreases from 22 to 20 Hz, the bearing comes into contact. The response curve bends to the high-frequency range at point E, a hardening effect.

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57

53

Fig. 9. (a) RMS dynamic response of PG-B carrier translational vibration with gravity (e) or mesh stiffness (e e) excitations. Bearing clearance is zero, infinite backlash is used, H ¼ harmonic, T ¼ translational mode, and R ¼ rotational mode; (b) RMS dynamic response of ring translational vibration when mesh stiffness is fluctuating (e) and constant (e e). The dotted lines denote the unstable branches. Dc ¼ 0.1, 0.5, 1 mm.

This resonance includes these two competitive nonlinear effects simultaneously. 6. Dynamic response of planetary gear sets with gravity Dynamic analyses of PG-A and PG-B were performed using the extended harmonic balance approach. The response included 100 harmonics of the carrier frequency, three harmonics of the mesh frequency, and eight upper and lower side bands of the mesh frequency harmonics. These parameters were selected based on the energy distribution in the frequency spectrum of the experimental data in Fig. 5. Vibration modes of two-dimensional planetary gear sets include rotational modes with distinct natural frequencies, translational modes with degenerate natural frequencies of multiplicity two, and planet modes with degenerate natural frequencies of multiplicity N  3. Translational modes include translation but no rotation of

the carrier, ring, and sun gears. The rotational modes considered have only rotation of the carrier, ring, and sun gears. The planet modes considered have only planet motions; the carrier, ring, and sun gears do not move. These vibration modes are derived from the Eigenvalue problem of the equations of motion in Eq. (1) without tooth contact and bearing clearance nonlinearities. The natural frequencies of PG-A under 800 Hz are the translational modes at 165, 311, and 796 Hz and the rotational mode at 292 Hz. The natural frequencies and mode shapes of PG-B are described in Table 3. 6.1. Gravity-induced excitation The effects of gravity and fluctuating mesh stiffness on the dynamic response of PG-B are shown in Fig. 9(a). Within the speed range of interest, the carrier response with gravity as the excitation source is one order of magnitude higher than that with fluctuating mesh stiffness. This agrees with the observation that low-frequency

Fig. 10. (a) Dynamic response of the RMS of carrier translational (upper graph) and rotational (lower graph) responses of PG-B with different backlash; (b) drive-side and back-side tooth contact at the first sun-planet mesh with one-third of nominal backlash.

54

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57

Fig. 11. RMS xr and x1 dynamic responses of PG-A with bearing clearance from 0 to 400 mm.

external excitations have a greater contribution to the dynamic response than the high-frequency internal excitation from gear meshing, as measured in the experiments shown in Fig. 5. In the dynamic response, with gravity as the excitation source, only the first and second translational modes are present. In the dynamic response, with fluctuating mesh stiffness as the excitation, the vibration amplitude is much lower (although more resonances are excited). The vibration resonances include the first two translational modes, the fourth harmonic of the third translational mode, and the second harmonic of the rotational mode. Fluctuating mesh stiffness consists of the fundamental and higher harmonics of the mesh frequency and can thus excite gearbox vibration modes beyond the frequency range studied in Fig. 9(a). Under gravity effects, the dynamic responses of ring translation with and without mesh stiffness variation are compared in Fig. 9(b) with various bearing clearance. The hardening effects shown in the figure when bearing clearance is larger than 0.5 mm is caused by bearing-raceway contacts. With increased bearing clearance, the resonances shifts to lower speed because of the reduced bearing stiffness, which is detailed in Section 6.2. When the mesh stiffness

variation is included, the vibration amplitude is higher than that with constant mesh stiffness. However, the resonance shapes are the same in nature. This indicates a mild coupling between gravity and the mesh stiffness fluctuation. 6.2. Nonlinear dynamic effects Fig. 10(a) (top) and (bottom) shows the carrier translational and rotational displacements of PG-B with different tooth backlash conditions. Backlash affects the translational mode resonance at 41 Hz as shown in Fig. 10(a) (top graph), but it does not affect the carrier rotational displacement as shown in Fig. 10(a) (bottom graph). With nominal backlash (bs ¼ br ¼ 482 mm), the shape of the resonance peak bends toward the right. This nonlinear effect is the hardening effect caused by tooth wedging at the resonance. As the backlash decreases, the shape of the resonance bends further to the right and the resonant frequency increases from 41 to 84 Hz. When the backlash is one-third of its nominal value, a closed loop caused by nonlinear tooth contact occurs at 70 Hz in the dynamic response of xc. Fig. 10(b) shows the tooth contact at both the

Fig. 12. Back-side tooth load at the first sun-planet mesh and planetary load sharing over a carrier cycle (a) with various bearing clearances and speeds; (b) with various sun stiffnesses and bearing clearances for PG-A. The carrier speed is 0.5 Hz in (b).

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57

drive-side and back-side of the first sun-planet mesh with onethird of the nominal backlash under the same condition as Fig. 10(a). In this figure, a y-axis value of “1” indicates a full-contact condition and a “0” represents an out-of-contact condition. Within the speed range of 70e85 Hz, the gear rotates significantly because of the rotational mode at 70 Hz and loses tooth contact. In the meantime, the translational mode causes the radial motion between the sun and planet 1. This radial motion exceeds the available tooth radial gap, which leads to tooth wedging. The coexistence of tooth wedging and tooth contact loss within the same speed range reflects the competition between the translational and rotational modes, resulting in the closed loop in Fig. 10(a). Fig. 11(a) and (b) shows the dynamic responses of xr and x1, with various values of bearing clearance for PG-A. The resonant frequency of the first translational mode decreases from 165 to 28 Hz when the bearing clearance increases from 0 to 400 mm as shown in Fig. 11(a). Bearing clearance affects the resonant frequencies of the translational mode shapes of the member with the bearing clearance. For large-scale wind turbine gearboxes, the resonant frequencies can be further reduced into the range of external lowspeed excitations because of structural flexibilities of the drivetrain, the blades, and the supporting components (Helsen et al., 2010). The responses without clearance and with clearance of 400 mm, which is essentially infinite clearance, correspond to two limiting systems. The mode shapes of these two limiting systems are different from each other as depicted in Fig. 11(a). In Fig. 11(a), the ring gear displacement increases 4.6 times when clearance increases from 0 to 400 mm. In the meantime, the planet translational displacement x1 decreases 25 mm when the bearing clearance increases. Clearance in the carrier bearing increases the vibration of the central members: the sun, ring, and carrier, but its effect on planet responses is much smaller. The bearing forces at the carrier support are small because they are nearly balanced by corresponding self-imposed tooth loads. Because of the absence of the

2

0 60 Gl ða; b; cÞ ¼ 6 4« 0

cosðaU1  bU2 Þs1 cosðaU1  bU2 Þs2 « cosðaU1  bU2 Þsn

sinðaU1  bU2 Þs1 sinðaU1  bU2 Þs2 « sinðaU1  bU2 Þsn

/ / / /

bearing clearance is larger than the backlash, the back-side tooth load and load sharing factor gradually increase with bearing clearance. When the speed is near the rpm that excites the resonant frequency, the tooth load and load sharing factor are higher than when the excitation speed is further away from the resonant frequency. The critical clearance when tooth wedging occurs is entirely independent of speed. The critical clearance value is reached when the bearing clearance equals the backlash. Therefore, the bearing clearance should be smaller than the backlash to avoid the potential of tooth wedging in wind turbine planetary gear sets. The sun support stiffness also affects tooth wedging. Fig. 12(b) shows the back-side tooth load over a carrier cycle at various sun support stiffness values. When the bearing clearance equals or is larger than the backlash, the tooth load rises exponentially with an increase in sun support stiffness. At zero sun support stiffness, no back-side tooth contact occurs, regardless of the value of bearing clearance. When the bearing clearance is smaller than the backlash, the back-side tooth load is zero for small sun support stiffness values. However, when the sun support stiffness has the same order of magnitude as the ring stiffness, tooth wedging occurs even with a smaller bearing clearance than the backlash. 7. Conclusions A harmonic balance method was developed to calculate the dynamic response of wind turbine planetary gear sets. Results obtained compare well with the GRC experiments of PG-A and the numerically integrated results of the lumped-parameter model of PG-B. This approach avoids protracted time integration simulations by calculating the dynamic response in the frequency domain. The time savings of this method make it suitable for parametric studies, which could be used to tune the planetary gear dynamics during the design phase to minimize vibratory response. Gravity introduces fundamental excitations in the rotating car-

sin cðaU1  bU2 Þs1 sin cðaU1  bU2 Þs2 « sin cðaU1  bU2 Þsn

nominal bearing force, the carrier-bearing clearance creates nonlinear dynamic behavior. Away from the resonance, the vibration amplitude is too low, so the bearing is not in contact. Near the resonance, the carrier-bearing rolling element displacement evolves gradually from no contact to partial contact, and eventually to full contact at the resonant frequency. The carrier-bearing rollers contacting the raceways at the resonant frequency induces the stiffening effects shown in Fig. 11(a) and (b). This stiffening effect reflects the increase in the overall bearing stiffness from the increased bearing contact (Guo and Parker, 2012a). 6.3. Tooth wedging Tooth wedging induced by gravity is a potential source for planet-bearing premature failures (Guo and Parker, 2010; Larsen et al., 2003; Rasmussen et al., 2004; Hansen et al., 2011). The influences of bearing clearance on tooth wedging and planetary load sharing are investigated in Fig. 12(a). The tooth loads are calculated for various values of the ratio between bearing clearance and backlash. When the bearing clearance is smaller than the backlash, there is no back-side tooth contact regardless of speed. When the

55

3 0 07 7 «5 0

(20)

rier frame of planetary gear sets. It plays a dominant role in the dynamics of wind turbine planetary gear sets compared to excitations caused by mesh stiffness fluctuation. Gravity causes tooth wedging and bearing contact at vibratory resonances, leading to stiffening effects in the dynamic response. A unique nonlinear effect induced by gravity illustrates the interaction of tooth wedging and tooth contact loss in the translational and rotational modes. Clearance in the carrier bearings reduces the resonant frequencies of the translational modes of the carrier, which are the lowest resonant frequencies of the examined systems. The shapes of these translational modes are affected by this bearing clearance. For wind turbine gearboxes larger than the ones studied, these frequencies could be reduced into the range of external low-speed excitations. Tooth wedging disrupts planetary load sharing and reduces the life of planet bearings. However, it can be prevented if the carrierbearing clearance is less than the tooth backlash.

56

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57

Acknowledgment

2

This work was supported by the U.S. Department of Energy under Contract No. DE-AC36-08GO28308 with the National Renewable Energy Laboratory.

3 Xðcbrj Þsin jbri sin ar Xðcbrj Þsin jbrj cos ar Xðcbrj Þsin jbrj 6 7 Lm;b ¼ 4 Xðcbrj Þcos jbrj sin ar Xðcbrj Þcos jbrj cos ar Xðcbrj Þcos jbrj 5 r2;j

Xðcbrj Þsin ar

Xðcbrj Þcos ar

Xðcbrj Þ

(27) 2

Xðcbsj Þsin jbsj cos jbsj Xðcbsj Þsin jbsj 6 b b b7 b b 2 b b 7 ¼6 Lm;b 4 Xðcsj Þsin jsj cos jsj Xðcsj Þcos jsj Xðcsj Þcos jsj 5 s1;j b b b b b Xðcsj Þsin jsj Xðcsj Þcos jsj Xðcsj Þ

Appendix A. Operators for the harmonic balance approach

2

1 cos Ui s1 sin Ui s1 cos 2Ui s1 6 1 cos Ui s2 sin Ui s2 cos 2Ui s2 Gi ¼ 6 4« « « « 1 cos U1 tn sin U1 sn cos 2U1 sn

3

3 / sin Ri Ui s1 0 / sin Ri Ui s2 0 7 7 / « «5 / sin Ri U1 sn 0

Xðcbsj Þsin2 jbsj

(28) 3

2

(19)

Xðcbsj Þsin jbsj sin as Xðcbsj Þsin jbsj cos as Xðcbsj Þsin jbsj 7 6 Lm;b ¼ 4 Xðcbsj Þcos jbsj sin as Xðcbsj Þcos jbsj cos as Xðcbsj Þcos jbsj 5 s2;j b b b Xðcsj Þcos as Xðcsj Þsin as Xðcsj Þ (29)

where i ¼ 1,2. where l ¼ 3, 4, plus or minus signs are used in the operator, respectively. When l ¼ 3, 4, c ¼ R5, R8, respectively.

2

n

Ai ¼ 4

3 A0i

Xðcbrj Þsin ar cos ar Xðcbrj Þsin ar 3 Lm;b ¼ 4 Xðcbrj Þsin ar cos ar Xðcbrj Þcos2 ar Xðcbrj Þcos ar 5 p;j 2

Xðcbrj Þsin2 ar Xðcbrj Þsin ar

5

(21)

n

B ¼ 4

n

1 2

2

3 B0i

0 6 1 6 6 6 0 Bi ¼ 6 6 6 6 4

5 n

0 2

2 0 1 0 Xi

7 7 7 7 7 7 7 7 Xi 5 0

1 2

h

i

(31)

(24)

where drj ; dsj denote the back-side mesh deflection vectors at the jth ring-planet and sun-planet meshes. Parameter s controls how close the smoothing functions are to the original piecewise nonlinear functions.

b

(32)

b



 b s L   b dbrj ¼ G zyr cos jbrj þzxr sin jbrj þzzi sin ar þzhi cos ar þzui zur  r L

dbsj ¼ G zys cos jbsj þzxs sin jbsj zzi sin as þzhi cos as zui zus 

jbsj ¼ jj þ as ; jbrj ¼ jj  ar (33)

(25)

where jj denotes the position angle of the jth planet. ar, as denote the pressure angles of the ring and sun gear teeth. br, bs denote the tooth backlashes of the ring and sun gears. References

where P ¼ diag(kl,u, ., kl,u), l ¼ c,r,s. By defining X() ¼ Hdiag()G,

2

i

cbsj ¼  kbsj ðtÞ s sec h2 ðsdbsj Þdbsj þ 1 þ tan hðsdbsj Þ

Appendix B. Jacobian matrix for the harmonic balance approach

KLTI ¼ diagð0; 0; Pc ; 0; 0; Pr ; 0; 0; Ps ; 0; .; 0Þ

h

(23) 3

1 0

Xðcbsj Þ

cbrj ¼  kbrj ðtÞ s sec h2 ðsdbrj Þdbrj þ 1 þ tan hðsdbrj Þ

where X1 ¼ R1, X2 ¼ R2, X3 ¼ R5, X4 ¼ R8.

2

Xðcbsj Þcos as

(30) (22)

i

Xðcbrj Þ

2 Xðcb Þsin2 a Xðcbsj Þsin as cos as Xðcbsj Þsin as 3 s sj Xðcbsj Þcos as 5 þ 4 Xðcbsj Þsin as cos as Xðcbsj Þcos2 as Xðcbsj Þsin as

  A0i ¼ diag 0; 12 ; 12 ; 22 ; 22 ; .; Xi2 ; Xi2 ; 0 ; i ¼ 1; 2; 3; 4

Xðcbrj Þcos ar

3

Xðcbrj Þsin2 jbrj Xðcbrj Þsin jbrj cos jbrj Xðcbrj Þsin jbrj 7 6 b b b7 b b 2 b b Lm;b ¼6 4 Xðcrj Þsin jrj cos jrj Xðcrj Þcos jrj Xðcrj Þcos jrj 5 r1;j Xðcbrj Þsin jbrj Xðcbrj Þcos jbrj Xðcbrj Þ (26)

Al-shyyab, A., Kahraman, A., 2005a. Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: period-one motions. J. Sound Vib. 284 (1e2), 151e172. Al-shyyab, A., Kahraman, A., 2005b. Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: subharmonic motions. J. Sound Vib. 279 (1e2), 417e451. Ambarisha, V.K., Parker, R.G., 2007. Nonlinear dynamics of planetary gears using analytical and finite element models. J. Sound Vib. 302, 577e595. Bahk, C.-J., Parker, R.G., 2011. Analytical solution for the nonlinear dynamics of planetary gears. ASME J. Comput. Nonlinear Dynam. 6 (2), 021007. Baker, G., Overman, E., 2000. The Art of Scientific Computing. William H. Press.

Y. Guo et al. / European Journal of Mechanics A/Solids 47 (2014) 45e57 Blankenship, G.W., Kahraman, A., 1996. Gear dynamics experiments. Part-I. Characterization of forced response. In: ASME Power Transmission and Gearing Conference, San Diego. Botman, M., 1976. Epicyclic gear vibrations. J. Eng. Ind. 97, 811e815. Friedmann, P.P., 1990. Numerical methods for the treatment of periodic systems with applications to structural dynamics and helicopter rotor dynamics. Comput. Struct. 35 (4), 329e347. Grippo, L., Lampariello, F., Lucidi, S., 1986. A nonmonotone line search technique for Newton’s method. SIAM J. Numer. Anal. 23 (4), 707e716. Guo, Y., Parker, R.G., 2010. Dynamic modeling and analysis of a spur planetary gear involving tooth wedging and bearing clearance nonlinearity. Eur. J. Mech. A Solids 29, 1022e1033. Guo, Y., Parker, R.G., 2012a. Dynamic analysis of planetary gears with bearing clearance. ASME J. Comput. Nonlinear Dynam. 7 (3), 031008. Guo, Y., Parker, R.G., 2012b. Stiffness matrix calculation of rolling element bearings using a finite element/contact mechanics model. Mech. Mach. Theor. 51, 32e45. Guo, Y., Keller, J., LaCava, W., 2012. Combined effects of input torque, non-torque load, gravity, and bearing clearance on planetary gear load share in wind turbine drivetrains. In: AGMA Fall Technical Meeting. Gurkan, N. E., Ozguven, H. N., 2007. Interactions between backlash and bearing clearance nonlinearity in geared flexible rotors. In: ASME International Design Engineering Technical Conferences DETC2007-34101. Hansen, A.M., Rasmussen, F., Larsen, T.J., 2011. Gearbox Loads caused by Double Contact Simulated with HAWC2. European Wind Energy Association. Helsen, J., Vanhollebeke, F., Coninck, F.D., Vandepitte, D., Desmet, W., 2010. Insights in wind turbine drive train dynamics gathered by validating advanced models on a newly developed 13.2 MW dynamically controlled test-rig. Mechatronics 21 (4), 737e752. Kahraman, A., 1992. On the response of a preloaded dynamic oscillator with a clearance: period-doubling and chaos. Nonlinear Dynam. 3, 183e198. Kahraman, A., Blankenship, G. W., 1994. Planet mesh phasing in epicyclic gear sets. In: International Gearing Conference, Newcastle, UK, pp. 99e104. Kahraman, A., Singh, R., 1991. Interactions between time-varying mesh stiffness and clearance non-linearities in a geared system. J. Sound Vib. 146 (1), 135e156. Kahraman, A., Vijayakar, S.M., 2001. Effect of internal gear flexibility on the quasistatic behavior of a planetary gear set. J. Mech. Des. 123 (3), 408e415. Larsen, T.J., Thomsen, K., Rasmussen, F., 2003. Dynamics of a Wind Turbine Planetary Gear Stage. Risø National Laboratory, Risø-I-2112 (EN).

57

Link, H., LaCava, W., Van Dam, J., McNiff, B., Sheng, S., Wallen, R., McDade, M., Lambert, S., Butterfield, S., Oyague, F., 2011. reportGearbox Reliability Collaborative Project Report: Findings from Phase 1 and Phase 2 Testing. National Renewable Energy Laboratory. Musial, W., Butterfield, S., McNiff, B., 2007. Improving wind turbine gearbox reliability. In: European Wind Energy Conference, Milan, Italy, NREL/CPe500e 41548. Nayfeh, A., Balachandran, B., 1995. Applied Nonlinear Dynamics. John Wiley and Sons. Oyague, F., Gorman, D., Sheng, S., 2010. NREL Gearbox Reliability Collaborative Experimental Data Overview and Analysis. NREL. Parker, R.G., Agashe, V., Vijayakar, S.M., Sep. 2000a. Dynamic response of a planetary gear system using a finite element/contact mechanics model. J. Mech. Des. 122 (3), 304e310. Parker, R.G., Vijayakar, S.M., Imajo, T., Oct. 2000b. Non-linear dynamic response of a spur gear pair: modelling and experimental comparisons. J. Sound Vib. 237 (3), 435e455. Raghothama, A., Narayanan, S., 1999. Bifurcation and chaos in geared rotor bearing system by incremental harmonic balance method. J. Sound Vib. 226 (3), 469e492. Rasmussen, F., Thomsen, K., Larsen, T., 2004. The Gearbox Problem Revisited. Risø National Laboratory, AED-RB-17 (EN). Seydel, R., 1994. Practical Bifurcation and Stability Analysis: From Equilibrium to Chaos, second ed., Springer. Singh, A., 2010. Load sharing behavior in epicyclic gears: physical explanation and generalized formulation. Mech. Mach. Theor. 45, 511e530. Singh, A., 2011. Epicyclic load sharing map e development and validation. Mech. Mach. Theor. 46, 632e646. Thomsen, J., 2003. Vibration and Stability: Advanced Theory, Analysis, and Tools, second ed., Springer. Turner, K., Walker, H.F., 1992. Efficient high accuracy solutions with gmres(m). J. Sci. Stat. Comput. 13 (3), 815e825. Velex, P., Flamand, L., Mar. 1996. Dynamic response of planetary trains to mesh parametric excitations. J. Mech. Des. 118 (1), 7e14. Vijayakar, S.M., 1991. A combined surface integral and finite element solution for a three-dimensional contact problem. Int. J. Numer. Methods Eng. 31, 524e546. Zhu, F., Parker, R.G., Jan. 2005. Non-linear dynamics of a one-way clutch in belt pulley systems. J. Sound Vib. 279 (1e2), 285e308.