Advantages and drawbacks of applying periodic time-variant modal analysis to spur gear dynamics

ARTICLE IN PRESS Mechanical Systems and Signal Processing 24 (2010) 1495–1508

Contents lists available at ScienceDirect

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

Advantages and drawbacks of applying periodic time-variant modal analysis to spur gear dynamics Rune Pedersen a, Ilmar F. Santos b,, Ivan A. Hede a a b

Siemens Wind Power A/S, Borupvej 16, DK-7330 Brande, Denmark Technical University of Denmark, Department of Mechanical Engineering, Nils Koppels Alle´, Bygning 403, DK-2800 Kgs. Lyngby, Denmark

a r t i c l e in fo

abstract

Article history: Received 18 August 2009 Received in revised form 18 November 2009 Accepted 22 December 2009 Available online 7 January 2010

A simplified torsional model with a reduced number of degrees-of-freedom is used in order to investigate the potential of the technique. A time-dependent gear mesh stiffness function is introduced and expanded in a Fourier series. The necessary number of Fourier terms is determined in order to ensure sufficient accuracy of the results. The method of time-variant modal analysis is applied, and the changes in the fundamental and the parametric resonance frequencies as a function of the rotational speed of the gears, are found. By obtaining the stationary and parametric parts of the timedependent modes shapes, the importance of the time-varying component relative to the stationary component is investigated and quantified. The method used for calculation and subsequent sorting of the left and right eigenvectors based on a first order Taylor expansion is explained. The advantages and drawbacks of applying the methodology to wind turbine gearboxes are addressed and elucidated. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Modal analysis Gear dynamics Parametric vibration Time-variant systems

1. Introduction Modal analysis method is very frequently used with the aim of solving time-invariant linear equations of motion. The eigenvalues and eigenvectors obtained from the solution of eigenvalue problems allow engineers to interpret and visualize the dynamic behavior of different mechanical systems. For time-variant equations of motion the use of modal analysis in its well-known form is not possible. In this case a stability analysis normally relies on Floquet theory. It does not deliver the complete homogenous solution, but gives only information about the stability of the system represented by equations of motion with time-variant coefficients [1]. In many cases as flexible rotating blades, flexible rotating discs, rotating shafts with non-symmetrical cross section, and gear dynamics, the coefficients of the equations of motion vary in a periodic way. In Xu and Gasch [2] the complete homogenous solution of periodic time-variant linear equations of motions is presented, based on Hill’s approximation. The periodic time-variant matrix systems are expanded in Fourier series. Assuming that the solution also can be expanded in Fourier series, it is possible to obtain a general homogenous solution by solving a hyper-eigenvalue problem, i.e. when the number of Fourier coefficients is not infinite. The solution of the hyper-eigenvalue problem leads also to eigenvalues and eigenvectors, nevertheless, the eigenvectors become also periodic time-variant and the eigenvalues become dependent on the periodicity of the parameter variation. Different contributions to the problem of periodic time-variant modal analysis are presented in the literature, i.e. non-symmetric rotors [1,3], flexible rotating discs [4,5], rotor-blade dynamics [2,6] and

 Corresponding author. Tel.: + 45 4525 6269; fax: + 45 4593 1577.

E-mail address: [email protected] (I.F. Santos). 0888-3270/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2009.12.009

ARTICLE IN PRESS 1496

R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

later active control of rotor-blade dynamics [7–10]. The mathematical foundations of modal analysis for time-varying linear systems is clearly and very nicely presented by Irretier in [11], and Bucher and Ewins in [12]. In [1] the dynamics of a simple flexible shaft with non-symmetric cross section, supported by anisotropic bearings, is theoretically investigated using Hill’s approach. In [3] such an investigation is carried out theoretically as well as experimentally. Flexible rotating discs are also an example of period time-variant structure. Their dynamics are carefully investigated using periodic-time variant modal analysis in [4,5]. The theoretical work presented in [4], and the experimental validation in [5], illustrates the continuation of Irretier’s work [13]. Rotating flexible blades are also an example of a periodic time-variant system. Their dynamics are also carefully investigated using periodic time-variant modal analysis, as it can be seen in [2,14,6,15]. In [15] a contribution to the experimental validation of linear and non-linear dynamic models for representing rotor-blades parametric coupled vibrations is given. The rotor-blade dynamics is described by using three models with different levels of complexity followed by experimental validation of such models. A deeper physical understanding of the dynamic coupling and the behavior of the parametric vibrations are achieved. Such an understanding is of fundamental importance while developing active control strategies. In [7] the design of time-variant modal controllers is in focus. Time-variant modal analysis and modern control theory are integrated in an elegant way allowing the development of new control strategies. The modal controllability and observability of bladed discs are strongly dependent on the angular velocity, a detailed analysis of such a dependency is presented in [8]. To control rotor and blade vibration using only shaft actuation is a very challenging problem. In [9] such a problem is investigated theoretically as well as experimentally using different control strategies. The electromagnetic actuators are used to control a horizontal rotor-blade system (blades periodically excited by the gravity). In [10] the same problem is theoretically as well as experimentally investigated and new strategies are developed to control vertical rotor-blade systems. The existing dynamic gear models can be classified based on the excitation source: transmission error (TE)-excited models and parametrically excited models. As shown by Blankenship and Singh [16], the TE and the mesh stiffness depend on each other, making the dynamic gear mesh modeling very complicated. To simplify the calculations, it is common to use the TE as the only external excitation source, and/or use the varying stiffness as a parametric excitation. The validity of this procedure is examined by Velex and Ajmi [17], who mention the problem of defining TE for helical gears as a limitation. Also the assumption of quasi-stationarity often used in TE-excited models is a disadvantage, as it excludes the possibility of using the model to correctly predict dynamic behavior in a resonant region. This problem is then partly solved by Kahraman and Singh in [18–21], who use the methods of non-linear dynamics for solving the equations of motion in a resonant region of a gear pair with a clearance-type non-linearity. The effect of the non-linearity is amplified by introducing the varying mesh stiffness. In the non-TE-excited models, no knowledge about the TE is needed before the calculation is performed. Peeken et al. [22–25] present models which are parametrically excited. It is shown that the equation of motion for the torsional 1-DOF (degree-of-freedom) system reduces to the Mathieu equation, when the gear mesh stiffness is replaced by a cosine function. More complicated mesh stiffness functions are introduced via the lowest terms of their Fourier series. In [26], Velex and Berthe solve the equations of motion after splitting these in a stationary part (results from the mean external load) and a dynamic part (from the varying part of the external load). The mesh stiffness is described by a Fourier series in the time domain. Velex and Maatar [27] and Ajmi and Velex [28] extend the model by further investigations concerning the mesh stiffness. The goal of this work is to apply the theory of time-variant modal analysis to spur gear dynamics. The advantages and drawbacks of using the technique for analyzing vibrations in spur gears are investigated.

2. Mathematical modeling 2.1. Gear mesh stiffness Finding the gear mesh stiffness is an iterative process, which involves calculation of the load distribution across the tooth and the load sharing between the contacting tooth pairs. The gear mesh model is shown schematically—and slightly simplified—in Fig. 1(a). In this example, three tooth pairs are taken into account, marked by 1, 2, and 3. Each tooth pair is divided into three sections in the longitudinal direction of the tooth, as shown in Fig. 1(b). The elastic coupling between sections within one tooth is handled by applying beam theory to simplify the plate deflection problem, as shown by Schmidt [32]. The goal for the stiffness model is to find the stiffness matrices KHertz for the Hertzian deformation, and Ktooth for the bending and shear in the gear tooth. For a low load P as shown in Fig. 1(a), only one section is in contact. With increasing load, more sections will be in contact with the mating gear flank, and thus contribute to the gear mesh stiffness. The stiffness model is based on the true involute tooth form, which is calculated from basic gear data using the model described by Padieth [29]. The stiffness calculations are performed using the methods of Weber and Banaschek [30], Ziegler [31], and Schmidt [32]. Following these references, the gear body deformation is included in Ktooth by regarding the gear body as an elastic, semi-infinite space. When applying a force to a tooth, the resulting stresses are transferred to the underlying gear body, which is assumed to deform in a zone that stretches 2–3 tooth widths from the loaded tooth.

ARTICLE IN PRESS R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

1497

Fig. 1. (a) Gear mesh stiffness model; (b) tooth pair discretization.

The inclusion of grinding corrections is an important aspect of the contact modeling. In Fig. 1(a), the topology to the right shows the grinding corrections. These will affect the number of contacting sections for a given load, which will then affect the gear mesh stiffness. The gear mesh stiffness is calculated for several relative positions of the meshing gears. This gives a load- and position dependent stiffness. Under the high loads considered in the following dynamic analysis, it is assumed that tooth separation will not occur. After completing the stiffness calculation, a constant rotational speed is assumed. By these assumptions, the gear mesh stiffness becomes a periodically varying parameter in the time domain, that is, the non-linearity is removed. It must also be noted that all teeth are considered identical, i.e. pitch errors are not included. The resulting system can be investigated using the theory for modal analysis of linear time-variant systems. 2.2. Modal analysis of time-variant system When the only time-variant parameter in the system is the gear mesh stiffness, the equation of motion can be written as Mq€ þCq_ þ KðtÞq ¼ f

ð1Þ

{ "

# " 0 q_  M1 KðtÞ q€

#"

I M

1

C

q q_

#

 ¼

0 0

0 M1

  0 f

ð2Þ

or "

_ zAðtÞz ¼ p;

q z¼ _ q

# ð3Þ

The first part of the analysis of Eq. (3) is to find the solutions to the homogeneous part of the equation, that is, solve _ zAðtÞz ¼ 0 for the unknown z. This is done by expanding the theory of modal analysis of time-invariant systems [33] to include periodically time-varying parameters, as shown in [6,2]. By the substitutions zðtÞ ¼ rðtÞelt and p ¼ 0, the homogeneous part of Eq. (3) can be written as the eigenvalue problem r_ ðtÞ þ ðlIAðtÞÞrðtÞ ¼ 0

ð4Þ

Under the assumption of constant rotational speed, the gear mesh frequency in units rad/s is called O. The stiffness matrix KðtÞ and therefore the state matrix AðtÞ are both periodic with period T ¼ 2p=O and can be expanded into infinite, complex Fourier series.pAlso ffiffiffiffiffiffiffi the time-variant eigenvector rj ðtÞ, belonging to the j th eigenvalue, is assumed to be periodic. In all equations, i ¼ 1: KðtÞ ¼

1 X

Kk eikOt

ð5Þ

Aa eiaOt

ð6Þ

k ¼ 1

AðtÞ ¼

1 X a ¼ 1

rj ðtÞ ¼

1 X r ¼ 1

rj;r eirOt

for j ¼ 1 . . . 2N

ð7Þ

ARTICLE IN PRESS 1498

R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

It can be shown that the components of the eigenvectors rj ðtÞ can be found by re-writing Eq. (4): 3 3 2 2   7 6 0 7 6r 7 6 j;2 7 6 7 7 6 6 6 rj;1 7 6 0 7 7 7 6 6 7 7 6 ^ 6 ^ AÞ ðlj I 6 rj;0 7 ¼ 6 0 7 7 7 6 6 6 rj;1 7 6 0 7 7 7 6 6 7 7 6 6 4 rj;2 5 4 0 5  

ð8Þ

where 2

& 6  6 6 6  6 6 A^ ¼ 6    6 6  6 6  4 

3













2iOI þA0 A1

A1 iOI þ A0

A2 A1

A3 A2

A4 A3

A2

A1

A0

A1

A2

A3

A2

A1

iOI þ A0

A1

 7 7 7  7 7  7 7 7  7 7  7 5 &

A4

A3

A2

A1

2iOI þA0











ð9Þ

In an exact solution to Eq. (8), A^ is infinitely large. However, the magnitude of the Fourier components of AðtÞ depend directly on the Fourier components of KðtÞ. Therefore, if the gear mesh stiffness function included in KðtÞ is smooth, only a few of the Fourier components of AðtÞ will have a magnitude that will significantly influence the resulting displacements and velocities in the vector zðtÞ. In a later section, it will be shown how to calculate the necessary number of Fourier components to be ^ in order to obtain accurate results. If the number of included Fourier components is n, in the sense that included in A, KðtÞ ¼

n X

Kk eikOt

ð10Þ

k ¼ n

and the number of degrees of freedom in the model is N, the size of matrix A^ will be 2Nð2n þ 1Þ  2Nð2n þ 1Þ. There is a great amount of redundant information in the solution of Eq. (8). Only the basis eigenvalues and basis eigenvectors, which are identified as described in a later section, are needed in the further analysis. The 2N basis eigenvalues are stored in a diagonal matrix K, and the Fourier components of the basis eigenvectors rj;0 are stored in the three-dimensional array R. This array can be visualized as ð2n þ1Þ ‘‘layers’’ of matrices of size 2N  2N. The k th layer of R will have the structure Rk ¼ ½r1;k r2;k    r2N1;k r2N;k  for k ¼ n; . . . ; n To find the solution q of Eq. (1), also the left eigenvectors LðtÞ, which are the solution to the equation RðtÞLðtÞ ¼ I are needed. These can be found by solving the matrix equation: RðtÞLðtÞ ¼ I { 1 X

Rr eirOt

r ¼ 1

{ 2

 6 6 6 6 6 6 6 6 6 6 6 4 

ð11Þ 1 X

Ll eilOt ¼ I

ð12Þ

l ¼ 1











R0

R1

R2

R3

R4

R1

R0

R1

R2

R3

R2

R1

R0

R1

R2

R3

R2

R1

R0

R1

R4 

R3 

R2 

R1 

R0 

3 3 2   6 7 7 6    76 L2 7 6 0 7 7 7 7 6 76 6 L1 7 6 0 7 7 7 7 6 76 7 7 6 76    76 L 0 7 ¼ 6 I 7 7 7 6 76 7 7 6 6 7 76 L 1 7 6 0 7 7 7 6 76    54 L 2 5 4 0 5    

32

ð13Þ

In most practical cases, the magnitudes of the submatrices Rn will decrease by increasing n. Therefore, a solution based on the central rows and columns of the infinitely large system of Eqs. (13) will be sufficiently accurate. For n ¼ 2, Eq. (13) can be approximated by 32 3 2 3 2 L2 R0 R1 R2 0 0 0 7 7 607 6R 6 R R R 0 L 0 1 2 76 1 7 6 7 6 1 76 7 6 7 6 7 6 7 6 R 2 R1 6 R0 R1 R2 7 ð14Þ 76 L 0 7 ¼ 6 I 7 6 7 7 6 7 6 0 6 R R R R L 0 4 5 4 2 1 0 1 54 1 5 0 0 R2 R1 R0 L2 0

ARTICLE IN PRESS R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

1499

With Ll found for l ¼ n    n, the solutions to the homogeneous and the inhomogeneous equation of motion, Eq. (1), can now be found [6,2]. After introducing the modal coordinates yðtÞ as _ _ ¼ RðtÞyðtÞ _ zðtÞ ¼ RðtÞyðtÞ ) zðtÞ þ RðtÞyðtÞ

ð15Þ

and inserting into Eq. (3), the terms are multiplied from the left with LðtÞ and rearranged. The result is the decoupled modal equation of motion 2 3 &   0 6 7 lk _ ð16Þ yðtÞ 4 5yðtÞ ¼ LðtÞ 1 fðtÞ M & P  By inserting LðtÞ ¼ nl¼ n Ll eilOt and the oscillating force fðtÞ ¼ p þ eiO t into the right hand side of (16), the equation becomes 2 3 &   n X 0 6 7 þ iðO þ lOÞt lk _ Ll ð17Þ yðtÞ 4 5yðtÞ ¼ 1 p e M l ¼ n & 

Now a modal solution yðtÞ ¼ yeiðO þ lOÞt is assumed and inserted into (17). The transformation back to the zðtÞ coordinates is Pn performed using Eq. (15), and RðtÞ is substituted by its equivalent sum of r ¼ n Rr . After rearranging the terms, the resulting equation shows the forced response: 3 2 & n n 7  0  6 X X 1 7 6 þ iðO þ Oðr þ lÞÞt Rr 6 ð18Þ zinhom ðtÞ ¼ 7Ll  1 p e O þl O Þ l ið 5 4 M k r ¼ n l ¼ n & When using a finite n and the initial condition zðt0 Þ ¼ z0 , the solution to the homogeneous equation of motion becomes 2 3 & 6 7 elk ðtt0 Þ ð19Þ zhom ðtÞ ¼ RðtÞ4 5Lðt0 Þz0 &

2.3. Sorting of eigenvalues and eigenvectors To find the basis eigenvectors and basis eigenvalues needed to set up RðtÞ and lk in Eq. (18), the eigenvalues and corresponding vectors must be sorted. This sorting is also necessary to correctly calculate the eigenvector normalization factors, as described in a later section. The problem can be described in Fig. 2, in which the imaginary part of the eigenvalues is plotted versus the gear mesh rotational frequency O. Both figures are zoomed to show only two ‘‘families’’ of eigenfrequencies. In Fig. 2(a), n ¼ 2 Fourier components are included in the analysis. Here, the sorting is easy: the 2n þ 1 lower frequencies belong to the low-frequency family centered around f1 . In Fig. 2(b), n ¼ 4 and the situation is more complicated—it is no longer a trivial task to decide which eigenvalues belong to which family. A typical eigenvalue distribution in the complex plane is shown in Fig. 3(a). It is clear, how the different families of eigenvalues can be identified by their real part, which is nearly identical for all members of the family [2]. This method will solve the problem visualized in Fig. 2(b). However, because of the transformation into state space from Eqs. (1) to (3), the eigenvalues are found in complex conjugate pairs, as depicted in Figs. 3(a) and (b). Obviously, these cannot be separated based on their real part alone. Instead, a sorting algorithm based on a first order Taylor expansion of the eigenvalues as a function of O is used. For a given Ok , where k is an integer index, the expected value for l is given by

lk;expected ¼ lk1 þ

dl l l DO  lk1 þ k1 k2 ðOk Ok1 Þ dO Ok1 Ok2

ð20Þ

The eigenvalues lk;expected predicted by Eq. (20) are then compared to those actually obtained by solving the hypereigenvalue problem, lk . The absolute difference between lk;expected and lk is calculated, and the eigenvalues in lk are then identified with the eigenvalue in lk;expected , which shows the minimum difference. A few requirements must be fulfilled in order to use the Taylor expansion sorting method: 1. The hyper-eigenvalue problem must be solved for a monotonically increasing value of O. 2. At least two solutions, for k2 and for k1 must be computed without the Taylor expansion sorting. 3. The values of O should not increase too much in each step in order to get a good estimation of lk from Eq. (20), i.e., a certain smoothness of the function lðOÞ is required.

ARTICLE IN PRESS 1500

R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

f2

Im (λ)

Im (λ)

f2

f1

f1

Angular frequency Ω [rad/sec]

Angular frequency Ω [rad/sec]

Fig. 2. Solutions to hyper-eigenvalue problem: (a) low n or low O; (b) high n or high O.

72 eigenvalues 8 basis eigenvalues

Im (λ)

Im (λ)

f1 0

0 −f1

0

Angular frequency Ω [rad/sec]

Re (λ) Fig. 3. Complex conjugate eigenvalues: (a) eigenvalue distribution, (b) high n or high O.

When the eigenvalues have been sorted within each family, the families are sorted relative to each other based on the mean value of imaginary part of the eigenvalues in the family. The basis eigenvalues are now defined as the central member of the family, based on the imaginary part. 3. Results 3.1. Presentation of model The gear pair, the attached shafts and the inertias J1 J4 are shown in Fig. 4. The degrees of freedom (DOF) of the model are the rotational angles of the four inertias and will be referred to as q ¼ ½q1 q2 q3 q4 T . The model has intentionally been kept very simple in order to keep the results easy to interpret. The basic data of the gears used in this work are presented in Table 1. The data are loosely based on the intermediate stage of a 1 MW wind turbine gear box, except a zero degree helix angle and zero profile shift on both gear and pinion is used. The calculated gear mesh stiffness function is shown in Fig. 5(a). For a constant rotational speed O, the periodic stiffness can be described using its complex Fourier expansion. When t is the time, the Fourier series is defined as f ðtÞ ¼

1 X j ¼ 1

cj eijt

ð21Þ

ARTICLE IN PRESS R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

J2

1501

J1

J4

J3

Fig. 4. System model.

Table 1 Basic gear data. Driving gear z

95

an

Unit

Description

22

– deg deg mm mm – mm Nm

Number of teeth Normal pressure angle Helix angle Center distance Normal module Profile shift Tooth width Torque

20 0 486 8

b

0 215 74875.84

0 225 17339.67

30

3

25

2.5

A = 2 abs (c) [N (mm μm)−1]

Gear mesh stiffness [N (mm μm)−1]

a mn x b Mt

Driven gear

20 15 10 Original Using 2 Fourier coeffs. Using 5 Fourier coeffs. Using 20 Fourier coeffs.

5 0 0

π

2π

3π

4π

Gear mesh position [rad]

5π

2 1.5 1 0.5 0

6π

1

5

10

15

20

25

30

Fourier component number

Fig. 5. (a) Gear mesh stiffness, (b) Fourier expansion of gear mesh stiffness.

In Fig. 5(b), the first 30 Fourier components of the signals are shown, expressed as Aj ¼ 2  absðcj Þ. Because of the symmetry in the complex Fourier components, where cj ¼ conjðcj Þ, only cj for j ¼ 1; . . . ; 30 are shown. The mean value of the gear mesh stiffness c0 ¼ 26:66 N ðmm mmÞ1 is not shown in order to illustrate the other components more clearly.

ARTICLE IN PRESS 1502

R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

3.2. Stiffness model validation The stiffness model described in Section 2.1 must be validated. This is done by comparing to stiffness results obtained using the commercial software KissSoft. The comparison is shown in Fig. 6. In Fig. 6(a) is shown the gear mesh stiffness, when deformations of the gear teeth are not included in the kinematic analysis to determine tooth contact positions. In this type of calculation, the total contact ratio for this spur gear pair equals the theoretical ea of 1.67. A good correlation between the models is seen. In Fig. 6(b), the elastic deformation in the gear mesh is included. In this case, the deformations cause the tooth to touch before the theoretical start of contact, and to stay in contact a longer period of time. This leads to an increase in contact ratio beyond the theoretical value. Also in this case a good correlation with the KissSoft results are shown. 3.3. Eigenfrequencies and mode shapes of the system To identify each of the 2N eigenfrequencies and the associated eigenvectors, a time-invariant (n ¼ 0) calculation has been performed. The 2N modes are numbered according to Table 2. The mode shapes rj ðtÞ can be evaluated at a given time t using the formula (7), using n and n as the limits for the sum. The result is shown in Fig. 7 for n ¼ 10, where the displacement parts of the eigenvectors are plotted versus the gear mesh position y in the interval ½0; 2p, corresponding to the time interval ½0; T at a constant gear mesh rotational frequency O. Modes 1 and 2—the rigid body modes—directly show the gear ratio 95/22. Modes 3 through 8 are elastic modes showing torsional deformations in the shafts and the gear mesh. It can be seen how modes 3–6 are strongly affected by the varying stiffness, while modes 1, 2, 7, and 8 are almost constant, r1 ðtÞ  r1 , r2 ðtÞ  r2 , r7 ðtÞ  r7 , and r8 ðtÞ  r8 . 3.4. ‘‘Quasi-static’’ eigenfrequency calculation As a preliminary analysis, the eigenfrequencies for the system can be calculated ‘‘quasi-statically’’ for a given time t ¼ t1 . To do this, AðtÞ in Eq. (4) is evaluated at t ¼ t1 using the first 2n þ 1 terms in Eq. (6), where n is the number of Fourier

30 Gear mesh stiffness [N (mm μm)−1]

Gear mesh stiffness [N (mm μm)−1]

30 25 20 15 c cγ

10

cth KissSoft

5

cth,γ KissSoft 0

25 20 15 10

c cγ

5

c KissSoft cγ KissSoft

0 0

π

2π

3π

4π

0

Gear mesh position [rad]

π

2π

3π

4π

Gear mesh position [rad]

Fig. 6. Validation of stiffness c and mean stiffness cg by comparing to program KissSoft: (a) no load-induced increase in contact ratio, (b) including loadinduced increase in contact ratio.

Table 2 Fundamental eigenfrequencies for n ¼ 0. Mode no.

Frequency (Hz)

1 2 3 4 5 6 7 8

0 0 399  399 585  585 3216  3216

ARTICLE IN PRESS R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

Mode 1, 2

1503

Mode 3, 4

1.5

1

1 q1 q2

0.5

0.5

q3 q4

0

0

−0.5

−0.5 0

π/2

π

3/2π

0

2π

Mode 5, 6

π/2

π

3/2π

2π

π/2 π 3/2π Gear mesh position [rad]

2π

Mode 7, 8

1.5

1.5

1

1

0.5

0.5

0

0

−0.5

−0.5 0

π/2 π 3/2π Gear mesh position [rad]

2π

0

Fig. 7. Mode shapes as a function of gear mesh position y ¼ Ot, calculated for n ¼ 10.

components included in the analysis: Aðt1 Þ ¼

n X

Aa eiaOt1

ð22Þ

a ¼ n

Aðt1 Þ is then inserted into Eq. (4), and the eigenvalue problem is solved. The resulting (positive) eigenfrequencies are shown in Fig. 8 for two different values of n. From the figure, it can be seen that the highest eigenfrequency at roughly 3200 Hz is much more sensitive to the variations in the gear mesh stiffness than the other eigenfrequencies. Also, n has an influence on the extent of the frequency interval, in which the highest eigenfrequency is located. As a conclusion to this ‘‘quasi-static’’ analysis, n can be expected to play an important role when determining the free and the forced response of the time-variant system. 3.5. Necessary number of Fourier components The accuracy of the time-variant modal analysis depends on the number of Fourier components, n, included in the analysis. When the applied force f in Eq. (1) is constant or zero, only the changing stiffness in KðtÞ can excite the system. From Eq. (10), it can be expected, that only frequencies up to approximately f ¼ nO will be excited. From the time-invariant modal analysis of the system it is known that the highest eigenfrequency is around 3215 Hz. When the gear mesh frequency is O ¼ 1357:59 rad=s  216 Hz, this eigenfrequency can be expected to be excited only when n Z3215=216  15. The accuracy of the method can be determined by comparing the solution in the time domain to the corresponding solution obtained by numerical integration. For this purpose, q3 has been chosen, as this DOF shows the largest difference between the two calculation methods, as seen from Fig. 9(a). It might be interesting to represent the time variations of the dynamic mesh force (as opposed to q3 in Fig. 9) since it also gives an indication on the presence of non-linearity (contact losses when the mesh force is negative). Since this loss of contact is not considered in the present gear mesh model, the dynamic mesh force is calculated as the difference between the displacements of the two gears multiplied by the mesh stiffness, FðtÞ ¼ cðtÞðq3 ðtÞr3 q2 ðtÞr2 Þ, where r2 and r3 are the radii of the two gears. Comparison based on a single DOF is the simplest way, since FðtÞ offers no new information on the dynamic behavior of the system. It is clear from Fig. 9(a) how the

ARTICLE IN PRESS R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

3500

3500

3000

3000 Eigenfrequencies [Hz]

Eigenfrequencies [Hz]

1504

2500 2000 1500 1000 500

2500 2000 1500 1000 500

0

0 π/2

0

π

3/2π

2π

0

Gear mesh position [rad]

π/2

π

3/2π

2π

Gear mesh position [rad]

Fig. 8. ‘‘Quasi-static’’ eigenfrequency analysis: (a) n ¼ 1; (b) n ¼ 18.

x 10−4

Numerical integration Time−variant modal analysis 10−4

2

q2 Difference [rad]

1 q3 (t) [rad]

q1

0 −1

q3

10−5

q4

10−6 −2

0.1

0.102 0.104 0.106 0.108 Time [s]

0.11

0

5 10 15 20 Number of Fourier coefficients, n

Fig. 9. Comparison of solutions: (a) maximum difference between numerical integration and time-variant modal analysis, (b) comparison of timedomain solutions (zoom). Only q3 ðtÞ is shown.

accuracy of the modal solution increases significantly when n reaches 15. An increase in n beyond 18 does not increase accuracy much. A zoom of the two solutions is shown in Fig. 9(b). 3.6. Change in fundamental frequencies For each number of n, the eigenvalue problem, Eq. (8), changes. It can therefore be expected that both the fundamental eigenfrequencies (for k ¼ 0) and the higher order parametric eigenfrequencies (for k 2 ½n; n; ka0) will change as a function of n. In the following analysis, the change in the fundamental eigenfrequencies, evaluated at the nominal speed O ¼ 1357:59 rad=s, is investigated. The results are shown in Fig. 10, where the change (in per cent) of the eigenfrequencies belonging to the three elastic modes, relative to a time-invariant modal analysis (n ¼ 0) are plotted. It can be seen that there are no significant changes in the fundamental eigenfrequencies. At no point, the frequencies change more than 0.16 per cent. 3.7. Normalization factors For each of the 2N eigenvectors, the importance of the k th harmonic parametric vibration mode can be calculated. This is done by finding the relative magnitude of the harmonic components of the eigenvector. For the k th harmonic

ARTICLE IN PRESS R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

1505

0 Modes 7 and 8 Modes 5 and 6 Modes 3 and 4

Change [per cent]

−0.05

−0.1

−0.15

−0.2 0

5

10

15

20

Number of Fourier coefficients, n Fig. 10. Change in fundamental frequencies.

component of the j th eigenvector, rj;k , the normalization factor NFj;k is defined as 2

1 NFj;k ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi T rj;k rj;k

rj;n

3

6  7 7 6 7 6 6 rj;1 7 7 6 7 6 with rj ¼ 6 rj;0 7 for j ¼ 1; . . . ; 2N 7 6 6 rj;1 7 7 6 6  7 5 4 rj;n

ð23Þ

In Figs. 11 and 12, normalization factors for all eight modes are plotted versus the gear mesh frequency O (in rad/s), for n ¼ 3 and 10, respectively. The NF-axes in the plots are logarithmic and show the interval NF 2 ½102 ; 101 . NF for the fundamental harmonic component, k ¼ 0, is set to NFj;0 ¼ 1 in all cases. With these definitions, the normalization factors show how the fundamental part of the mode shape vector r0 is related to the time-varying parts rk for k ¼ n; . . . ; n depending on the operational speed O. The k th harmonic component showing the smallest NFj;k will be predominant in the j th mode shape at a particular gear mesh frequency O. A number of observations can be made from the normalization factor plots, Figs. 11 and 12: 1. For the rigid-body modes 1 and 2, the fundamental harmonic k ¼ 0 is predominant for all gear mesh frequencies O 4 0. 2. For all elastic modes, there is a symmetry between the positive (left column) and the negative frequencies (right column in the plots): For the j th positive eigenfrequency, the normalization factor NFj;k is equal to NFj;k for the corresponding negative eigenfrequency. 3. For modes 3–6, there exist certain mesh frequency ranges, for which the time-varying part of the mode shape (parametric vibrations) will be more important than the stationary. These frequency ranges are independent of the number of Fourier components included in the analysis. For instance, for mode 3 in Fig. 11, NF for the þ2 harmonic component shows a local minimum around O ¼ 600 rad=s. This ‘‘potential þ 2 harmonic resonance’’ is found for all n Z 2. 4. Modes 7 and 8 at roughly f ¼ 7 3216 Hz behave fundamentally different from the other elastic modes (modes 3–6). No n-independent frequencies with low NF are observed. 5. For modes 7 and 8, frequency intervals exist in which NF for one or more harmonic components of the mode shape are smaller than 1. In these frequency intervals, which depend on the number of Fourier components included but are generally located in O 2 ½0; 600 rad=s, the content of the time-varying part of the mode in the overall mode shape is larger than the fundamental component. These modes are strongly dependent on the relative angular movement between the two gears, and are strongly affected by the time-variant tooth stiffness.

ARTICLE IN PRESS 1506

R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

Mode 1

−3

Mode 2

−1

100

0

NF

NF

100

−2

10−2

1 2

10−2 0

500

1000

1500

0

Mode 3

500

1000

1500

3

Mode 4

NF

100

NF

100

10−2

10−2 0

500

1000

1500

0

Mode 5

500

1000

1500

Mode 6

NF

100

NF

100

10−2

10−2 0

500

1000

1500

0

Mode 7

500

1000

1500

Mode 8

NF

100

NF

100

10−2

10−2 0

500 1000 Ω [rad/sec]

1500

0

500 1000 Ω [rad/sec]

1500

Fig. 11. Normalization factors, n ¼ 3.

4. Conclusion The theory of modal analysis of time-variant systems has been applied to a simple spur gear pair with a periodically time-varying gear mesh stiffness. It has been made clear that a large number of terms in the Fourier expansion of the system matrices is necessary in order to yield results of sufficient accuracy. This is a direct result of the jumps in the gear mesh stiffness function. In the cases of flexible rotors with non-symmetrical cross section, flexible rotating discs, and flexible rotating blades, a very reduced number of Fourier components is needed, normally n ¼ 2, as the time-variant coefficients are normally sine and cosine functions. It has been shown that there exist regions, in which the higher order parametric contributions to the overall mode shapes will be significant. For the case studied here, these effects were mainly observed at low gear mesh frequencies. The elastic mode with the highest frequency behaved differently from the other elastic modes. While the two lower modes showed parametric resonance frequencies that were largely independent of the number of Fourier components included in the analysis n, the parametric resonance areas of the high frequency mode strongly depended on n. Overall, the vibrations related to the higher frequency mode seemed to be more sensitive to the time-variant nature of the gear mesh stiffness. The system studied in this work consisted of a single gear stage. A typical modern wind turbine gearbox consists of one or two planetary stages followed by one or two parallel gear stages, with a total of 8–15 gear meshes. When applying the theory of modal analysis of time-variant systems to such a complex system, great care must be taken when interpreting the results. The larger the number of Fourier components needed to expand the periodic time-variant coefficients, the larger the hyper-eigenvalue problem becomes. It means also, it becomes more complicated and complex to physically interpret the basic and parametric mode shapes. Nevertheless, by exploring the definition of NF, the importance of the time-varying part of the mode shapes (parametric modes) can be investigated and quantified as a function of the gear mesh frequency. Compared to time-step integration schemes, the main advantage of the time-variant modal analysis is

ARTICLE IN PRESS R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

Mode 1

Mode 2 −10

100

−9

NF

NF

100

10−2

−8 −7

10−2 0

500 1000 Mode 3

1500

0

500 1000 Mode 4

1500

−1 0

10−2 500 1000 Mode 5

1500

0

500 1000 Mode 6

1500

1 2 3 4

100

5

NF

NF

100

−5

−2

NF 10−2

−6 −4 −3

100

NF

100

0

1507

6 −2

−2

10

10

0

500 1000 Mode 7

1500

0

500 1000 Mode 8

1500

7 8 9 10

NF

100

NF

100

10−2

10−2 0

500 1000 Ω [rad/sec]

1500

0

500 1000 Ω [rad/sec]

1500

Fig. 12. Normalization factors, n ¼ 10.

that it offers an analytical solution to the vibration problem. Therefore a modal truncation is possible, which removes the need for the very short time step used for a numerical integration of a system with high eigenfrequencies. Also the method allows to expand the analysis to the concepts of observability and controllability [8], which offer a quantification of the parametric vibrations.

References [1] M. Ertz, A. Reister, R. Nordmann, Zur Berechnung der Eigenschwingungen von Strukturen mit periodisch zeitvarianten Bewegungsgleichungen, in: H. Irretier, R. Nordmann (Eds.), Schwingungen in rotierenden Maschinen, vol. III, Springer, Vieweg, Braunschweig, Germany, 1995, pp. 288–296. [2] J. Xu, R. Gasch, Modale Behandlung linearer periodisch zeitvarianter Bewegungsgleichungen, Archive of Applied Mechanics—Ingenieur Archiv 65 (3) (1995) 178. [3] F.E. Boru, H. Irretier, Numerical and experimental dynamic analysis of a rotor with non-circular shaft mounted in anisotropic bearings, in: SIRM 2009—8th International Conference on Vibrations in Rotating Machines, Vienna, Austria, 2009, pp. 1–10. ¨ [4] H. Irretier, F. Reuter, Frequenzgange rotierender periodisch zeitvarianter Systeme, in: H. Irretier, R. Nordmann (Eds.), Schwingungen in rotierenden Maschinen, vol. IV, Springer, Vieweg, Braunschweig, Germany, 1997, pp. 113–121. [5] F. Reuter, Coupling of elastic and gyroscopic modes of rotating disc structures, in: H. Irretier, R. Nordmann (Eds.), Fifth International Conference on Rotor Dynamics, Darmstadt, Germany, 1998, pp. 443–455. [6] I.F. Santos, C.M. Saracho, Modal analysis in periodic, time-varying systems with emphasis to the coupling between flexible rotating beams and nonrotating flexible structures, in: Proceedings of the Xth International Symposium on Dynamic Problems of Mechanics, Sa~ o Paulo, Brazil, 2003, pp. 399–404. [7] R. Christensen, I. Santos, Design of active controlled rotor-blade systems based on time-variant modal analysis, Journal of Sound and Vibration 280 (3–5) (2005) 863–882. [8] R. Christensen, I. Santos, Modal controllability and observability of bladed disks and their dependency on the angular velocity, Journal of Vibration and Control 11 (6) (2005) 801–828. [9] R. Christensen, I. Santos, Active rotor-blade vibration control using shaft-based electromagnetic actuation, Transactions of the ASME. Journal of Engineering for Gas Turbines and Power 128 (3) (2006) 644–652. [10] R.H. Christensen, I.F. Santos, Control of rotor-blade coupled vibrations using shaft-based actuation, Shock and Vibration 13 (4–5) (2006) 255–271. [11] H. Irretier, Mathematical foundations of experimental modal analysis in rotor dynamics, Mechanical Systems and Signal Processing 13 (2) (1999) 183–191. [12] I. Bucher, D.J. Ewins, Modal analysis and testing of rotating structures, Philosophical Transactions of the Royal Society of London A 359 (2001) 61–96.

ARTICLE IN PRESS 1508

R. Pedersen et al. / Mechanical Systems and Signal Processing 24 (2010) 1495–1508

[13] H. Irretier, F. Reuter, Experimentelle Modalanalyse an einer rotierenden Scheibe, in: H. Irretier, R. Nordmann (Eds.), Schwingungen in rotierenden Maschinen, vol. I, Springer, Vieweg, Braunschweig, Germany, 1991, pp. 66–77. [14] J. Bienert, Anwendung der Strukturmodifikation zur Vorhersage der Kreiselwirkung von symmetrischen und unsymmetrischen Rotoren, in: H. Irretier, R. Nordmann (Eds.), Schwingungen in rotierenden Maschinen, vol. IV, Springer, Vieweg, Braunschweig, Germany, 1997, pp. 97–104. [15] I. Santos, C. Saracho, J. Smith, J. Eiland, Contribution to experimental validation of linear and non-linear dynamic models for representing rotor-blade parametric coupled vibrations, Journal of Sound and Vibration 271 (3–5) (2004) 883–904. [16] G.W. Blankenship, R. Singh, Comparative study of selected gear mesh interface dynamic models, Advancing Power Transmission into the 21st Century and American Society of Mechanical Engineers, Design Engineering Division (Publication) DE 43 pt 1, 1992, pp. 137–146. [17] P. Velex, M. Ajmi, On the modelling of excitations in geared systems by transmission errors, Journal of Sound and Vibration 290 (3) (2006) 882–909. [18] A. Kahraman, R. Singh, Non-linear dynamics of a spur gear pair, Journal of Sound and Vibration 142 (1) (1990) 49–75. [19] A. Kahraman, R. Singh, Non-linear dynamics of a geared rotor-bearing system with multiple clearances, Journal of Sound and Vibration 144 (3) (1991) 469–506. [20] A. Kahraman, R. Singh, Interactions between time-varying mesh stiffness and clearance non-linearities in a geared system, Journal of Sound and Vibration 146 (1) (1991) 135–156. [21] A. Kahraman, R. Singh, Dynamics of an oscillator with both clearance and continuous non-linearities, Journal of Sound and Vibration 153 (1) (1992) 180–185. [22] H. Peeken, C. Troeder, G. Diekhans, Parametererregte Getriebeschwingungen, Teil 1 VDI-Z 122 (20) (1980) 869–877. [23] H. Peeken, C. Troeder, G. Diekhans, Parametererregte Getriebeschwingungen, Teil 2 VDI-Z 122 (21) (1980) 967–977. [24] H. Peeken, C. Troeder, G. Diekhans, Parametererregte Getriebeschwingungen, Teil 3 VDI-Z 122 (22) (1980) 1029–1043. [25] H. Peeken, C. Troeder, G. Diekhans, Parametererregte Getriebeschwingungen, Teil 4 VDI-Z 122 (23/24) (1980) 1101–1113. [26] P. Velex, D. Berthe, Dynamic tooth loads on geared train, in: Proceedings of the 1989 International Power Transmission and Gearing Conference: New Technologies for Power Transmissions of the 1990s, 1989, pp. 447–454. [27] P. Velex, M. Maatar, A mathematical model for analyzing the influence of shape deviations and mounting errors on gear dynamic behaviour, Journal of Sound and Vibration 191 (5) (1996) 629–660. [28] M. Ajmi, P. Velex, A model for simulating the quasi-static and dynamic behaviour of solid wide-faced spur and helical gears, Mechanism and Machine Theory 40 (2) (2005) 173–190. [29] R. Padieth, Exakte Ermittlung der Zahnform, Antriebtechnik 17 (10) (1978) 434–436. ¨ ¨ ¨ ¨ [30] C. Weber, K. Banaschek, Formanderung und Profilrucknahme bei gerad- und schragverzahnten Radern, Schriftenreihe Antriebstechnik (Heft 11). ¨ ¨ ¨ [31] H. Ziegler, Verzahnungssteifigkeit und Lastverteilung schragverzahnter Stirnrader, Ph.D. Thesis, Rheinisch-Westfalische Technische Hochschule Aachen, Germany, 1971. ¨ ¨ ¨ [32] G. Schmidt, Berechnung der W¨alzpressung schragverzahnter Stirnrader unter Berucksichtigung der Lastverteilung, Ph.D. Thesis, Technische ¨ Munchen, ¨ Universitat Germany, 1973. [33] L. Meirovitch, Fundamentals of Vibrations, McGraw-Hill, 2001.