**4. Mesh stiffness models – Parametric excitations**

#### **4.1 Classic thin-slice approaches**

84 Mechanical Engineering

Note that the determinant of the stiffness matrix is zero which indicates a rigid-body mode (the mass matrix being diagonal). After multiplying the first line in (20) by *Rb I*1 02 , the second line by *Rb I*2 01 , adding the two equations and dividing all the terms by 2 2

semi-definite system (20) is transformed into the differential equation:

, cos *t b*

**3.6.2 A simple torsional-flexural model for spur gears** 

 ,

1

1 *k*

Fig. 5. Simplified torsional-flexural spur gear model.

The general expression of the structural vector **V M** (9) reduces to:

Re-writing the degree of freedom vector as 1 11 2 22 \*

parametrically excited differential system is obtained for linear free vibrations:

, equivalent mass

With 11 22 *x Rb Rb* 

*Rb I Rb I*

constant.

*m*

02 01 2 2 1 02 2 01 *I I*

gear axes of rotation 1 2

, relative apparent displacement

(bearing/shaft equivalent stiffness elements for instance) must be added.

*k t* ,**q**

when the pinion speed 1 and the output torque *Cr* are supposed to be

The simplest model which accounts for torsion and bending in spur gears is shown in Figure 5. It comprises 4 degrees of freedom, namely: 2 translations in the direction of the line of action 1 2 *V V*, (at pinion and gear centres respectively) and 2 rotations about the pinion and

> <sup>2</sup> *k V*<sup>1</sup> *V*2

> > 1 2 1 1 *<sup>T</sup>* **V0**

 

*<sup>T</sup>* **q** *v Rb v Rb* 

*Rb Rb* (22)

**Mq \* K q\* 0** *t* (23-1)

, the following

2  

(21)

. Because of the introduction of bending DOFs, some supports

,

*Ltx <sup>d</sup> mx k t x x F k M e M dM NLTE*

10 2 20 1 *I Rb I Rb* , the

2 2

*dt*

 

> From the results in section 2-5, it can be observed that, in the context of gear dynamic simulations, the mesh stiffness function defined as , , *L t k t k M dM* **q q** plays a key role.

> This function stems from a 'thin-slice' approach whereby the contact lines between the mating teeth are divided in a number of independent stiffness elements (with the limiting case presented here of an infinite set of non-linear time-varying elemental stiffness elements) as schematically represented in Figure 6.

Fig. 6. 'Thin-slice' model for time-varying mesh stiffness.

Since the positions of the teeth (and consequently the contact lines) evolve with time (or angular positions), the profiles slide with respect to each other and the stiffness varies because of the contact length and the individual tooth stiffness evolutions. The definition of mesh stiffness has generated considerable interest but mostly with the objective of calculating accurate static tooth load distributions and stress distributions. It has been shown by Ajmi and Velex (2005) that a classic 'thin-slice' model is sufficient for dynamic calculations as long as local disturbances (especially near the tooth edges) can be ignored. In this context, Weber and Banascheck (1953) proposed a analytical method of calculating tooth deflections of spur gears by superimposing displacements which arise from i) the contact

On the Modelling of Spur and Helical Gear Dynamic Behaviour 87

Although the results above are based on simplified bi-dimensional approaches, they are still widely used in gear design. For example, the mesh stiffness formulae in the ISO standard 6336 stem from Weber's analytical formulae which were modified to bring the values in closer agreement with the experimental results. Another important simplification brought by the ISO formulae is that the mesh stiffness per unit of contact length *k*0 is considered as

*k M dM k dM k L t*

Considering involute profiles, the contact lines in the base plane are inclined by the base

multiples of the apparent base pitch *Pba* and, when the pinion and the gear rotate, they all

*Pba*

*<sup>T</sup>*<sup>1</sup> *' T2*

 *Pba* 

Fig. 8. Base plane and contact lines ( *b* : face width; **z** : axial direction (direction of the axes

It transpires from this geometrical representation that the total length of contact between the pinion and the gear is likely to vary with time and, based on the simple stiffness equation (24), that mesh stiffness is time-varying and, consequently, contributes to the system

The extent of action on the base plane is an important property measured by the contact

which, in simple terms, represents the 'average number' of tooth pairs in contact

 *Pba* 

of rotation); 1 2 *T T*, : points of tangency on pinion and gear base circles and ' '

0 0

*<sup>b</sup>* (Figure 8) which is nil for spur gears. All contact lines are spaced by integer

,

**q** (24)

*T*2

**X**

1 2 *T T*, : limits of

approximately constant so that the following approximation can be introduced:

, ,

**q q**

*L t L t*

 

where *L t* ,**q** is the time-varying (possibly non-linear) contact length.

**4.2 Contact length variations for external spur and helical gears** 

undergo a translation in the **X** direction at a speed equal to *Rb*1 1 .

*b*

*'T1*

Contact line

helix angle

*b*

ratio

 

**z**

the contact area on base plane).

excitation via parametric excitations.

(possibly non integer) and is defined by:

between the teeth, ii) the tooth itself considered as a beam and, iii) the gear body (or foundation) influence. An analytical expression of the contact compliance was obtained using the 2D Hertzian theory for cylinders in contact which is singular as far as the normal approach between the parts (contact deflection) is concerned. The other widely-used formulae for tooth contact deflection comprise the analytical formula of Lundberg (1939), the approximate Hertzian approach originally used at Hamilton Standard (Cornell, 1981) and the semi-empirical formula developed by Palmgren (1959) for rollers. The tooth bending radial and tangential displacements were derived by equating the work produced by one individual force acting on the tooth profile and the strain energy of the tooth assimilated to a cantilever of variable thickness. Extensions and variants of the methodology were introduced by Attia (1964), Cornell (1981) and O'Donnell (1960, 1963) with regard to the foundation effects. Gear body contributions were initially evaluated by approximating them as part of an elastic semi-infinite plane loaded by the reactions at the junction with the tooth. A more accurate expression for this base deflection has been proposed by Sainsot *et al.* (2004) where the gear body is simulated by an elastic annulus instead of a half-plane. Figure 7 shows two examples of mesh stiffness functions (no contact loss) calculated by combining Weber's and Lundberg's results for a spur and a helical gear example. It can be observed that the stiffness fluctuations are stronger in the case of conventional spur gears compared with helical gears for which the contact variations between the teeth are smoother.

Fig. 7. Examples of mesh stiffness functions for errorless gears.

between the teeth, ii) the tooth itself considered as a beam and, iii) the gear body (or foundation) influence. An analytical expression of the contact compliance was obtained using the 2D Hertzian theory for cylinders in contact which is singular as far as the normal approach between the parts (contact deflection) is concerned. The other widely-used formulae for tooth contact deflection comprise the analytical formula of Lundberg (1939), the approximate Hertzian approach originally used at Hamilton Standard (Cornell, 1981) and the semi-empirical formula developed by Palmgren (1959) for rollers. The tooth bending radial and tangential displacements were derived by equating the work produced by one individual force acting on the tooth profile and the strain energy of the tooth assimilated to a cantilever of variable thickness. Extensions and variants of the methodology were introduced by Attia (1964), Cornell (1981) and O'Donnell (1960, 1963) with regard to the foundation effects. Gear body contributions were initially evaluated by approximating them as part of an elastic semi-infinite plane loaded by the reactions at the junction with the tooth. A more accurate expression for this base deflection has been proposed by Sainsot *et al.* (2004) where the gear body is simulated by an elastic annulus instead of a half-plane. Figure 7 shows two examples of mesh stiffness functions (no contact loss) calculated by combining Weber's and Lundberg's results for a spur and a helical gear example. It can be observed that the stiffness fluctuations are stronger in the case of conventional spur gears compared

with helical gears for which the contact variations between the teeth are smoother.

a - Spur gear

b - Helical gear

Fig. 7. Examples of mesh stiffness functions for errorless gears.

30 Although the results above are based on simplified bi-dimensional approaches, they are still widely used in gear design. For example, the mesh stiffness formulae in the ISO standard 6336 stem from Weber's analytical formulae which were modified to bring the values in closer agreement with the experimental results. Another important simplification brought by the ISO formulae is that the mesh stiffness per unit of contact length *k*0 is considered as approximately constant so that the following approximation can be introduced:

$$\int\_{L(t,\mathbf{q})} k(M)dM \equiv k\_0 \int\_{L(t,\mathbf{q})} dM = k\_0 L\left(t, \mathbf{q}\right) \tag{24}$$

where *L t* ,**q** is the time-varying (possibly non-linear) contact length.

#### **4.2 Contact length variations for external spur and helical gears**

Considering involute profiles, the contact lines in the base plane are inclined by the base helix angle *<sup>b</sup>* (Figure 8) which is nil for spur gears. All contact lines are spaced by integer multiples of the apparent base pitch *Pba* and, when the pinion and the gear rotate, they all undergo a translation in the **X** direction at a speed equal to *Rb*1 1 .

Fig. 8. Base plane and contact lines ( *b* : face width; **z** : axial direction (direction of the axes of rotation); 1 2 *T T*, : points of tangency on pinion and gear base circles and ' ' 1 2 *T T*, : limits of the contact area on base plane).

It transpires from this geometrical representation that the total length of contact between the pinion and the gear is likely to vary with time and, based on the simple stiffness equation (24), that mesh stiffness is time-varying and, consequently, contributes to the system excitation via parametric excitations.

The extent of action on the base plane is an important property measured by the contact ratio which, in simple terms, represents the 'average number' of tooth pairs in contact (possibly non integer) and is defined by:

On the Modelling of Spur and Helical Gear Dynamic Behaviour 89

b. it can observed that the time-varying part of the contact length disappears when either

c. harmonic analysis is possible by setting *k* 1,2,... in (27) and it is possible to represent the contact length variations for all possible values of profile and overlap contact ratios on a unique diagram. Figure 10 represents the RMS of contact length variations for a

> 

 ( 1 

1 1.5 2 2.5

Mesh stiffness can be determined using the Finite Elements Method but it is interesting to have orders of magnitude or approximate values at the design stage. For solid gears made of steel, an order of magnitude of the mesh stiffness per unit of contact length 0 *<sup>k</sup>* is 1.3 10 10 N/m². More accurate expressions can be derived from the ISO 6336 standard which,

0.8 *<sup>k</sup>* cos

*q*

 

 

above 1, contact length variations are very limited.

is below 2 and

 

has to avoided for a continuous

0.05 0.1 0.15

for a range of profile and transverse contact

0.2

(27)

below 1

The following conclusions can be drawn:

 

is an integer


0

ratios.

for solid gears, gives:

Fig. 10. Contour plot of the R.M.S. of *L L* / *<sup>m</sup>*

**4.3 Approximate expressions – Orders of magnitude** 

0.5

1

1.5

 

2

2.5

and *Sinc k* 1

realistic range of contact and overlap ratios. It shows that:

 

0.3 0.25 0.2 0.15 0.1 0.05

0

0.05

0.05


 


motion transfer) and /or 1

 

a. for spur gears, 0

 or 

$$\varepsilon\_{\alpha} = \frac{T\_1^{\prime} T\_2^{\prime}}{P b\_{\varrho}} = \frac{\sqrt{R a\_1^{\prime} - R b\_1^{\prime} \,^2} + \sqrt{R a\_2^{\prime} - R b\_2^{\prime} \,^2} - E \sin \alpha\_{\varrho}}{\pi \, m \cos a\_{\varrho}}\tag{25-1}$$

with 1 2 *Ra Ra* , : external radius of pinion, of gear; 1 2 *Rb Rb* , : base radius of pinion, of gear; *E OO* 1 2 : centre distance

In the case of helical gears, the overlap due to the helix is taken into account by introducing the overlap ratio defined as:

$$
\varepsilon\_{\beta} = \frac{b \tan \beta\_b}{P b\_a} = \frac{1}{\pi} \frac{b}{m} \frac{\tan \beta\_b}{\cos \alpha\_p} \tag{25-2}
$$

and the sum is defined as the total contact ratio.

Introducing the dimensionless time *m t T* where 1 1 *a m Pb <sup>T</sup> Rb* is the mesh period i.e. the time needed for a contact line to move by a base pitch on the base plane, a closed form expression of the contact length *L* for ideal gears is obtained under the form (Maatar & Velex, 1996), (Velex et al., 2011):

$$\frac{L(\tau)}{L\_m} = 1 + 2 \sum\_{k=1}^{n} \operatorname{Sinc}(k \varepsilon\_a) \operatorname{Sinc}(k \varepsilon\_\beta) \cos \left( \pi k \left( \varepsilon\_a + \varepsilon\_\beta - 2\tau \right) \right) \tag{26}$$

with: cos *<sup>m</sup> b <sup>b</sup> <sup>L</sup>* , average contact length

 sin *<sup>x</sup> Sinc x x* is the classic sine cardinal function which is represented in Figure 9.

Fig. 9. Evolutions of sin *<sup>x</sup> Sinc x x* .

The following conclusions can be drawn:

88 Mechanical Engineering

cos

 

cos *b b a p*

> *m Pb <sup>T</sup>*

(25-2)

*Rb* is the mesh period i.e. the

1 1 *a*

> 

 

(26)

for ideal gears is obtained under the form (Maatar &

   

(25-1)

*p*

22 22 ' ' 11 22 1 2 sin

*T T Ra Rb Ra Rb E*

with 1 2 *Ra Ra* , : external radius of pinion, of gear; 1 2 *Rb Rb* , : base radius of pinion, of gear;

In the case of helical gears, the overlap due to the helix is taken into account by introducing

tan 1 tan

*b b Pb m*

*m t T* 

where

time needed for a contact line to move by a base pitch on the base plane, a closed form

1 2 cos 2

is the classic sine cardinal function which is represented in Figure 9.

 

is defined as the total contact ratio.

*a p*

*Pb m*

*<sup>L</sup> Sinc k Sinc k k*

1

*x* 

.

 

*m k*

, average contact length

defined as:

: centre distance

 

> 

Introducing the dimensionless time

expression of the contact length *L*

*L*

*b*

Fig. 9. Evolutions of sin *<sup>x</sup> Sinc x*

*x* 

Velex, 1996), (Velex et al., 2011):

cos *<sup>m</sup>*

*<sup>b</sup> <sup>L</sup>* 

sin *<sup>x</sup> Sinc x*

 

*E OO* 1 2

the overlap ratio

and the sum

with:


	- contact length variations are significant when is below 2 and below 1
	- contact length is constant when 2 ( 1 has to avoided for a continuous motion transfer) and /or 1


Fig. 10. Contour plot of the R.M.S. of *L L* / *<sup>m</sup>* for a range of profile and transverse contact ratios.

#### **4.3 Approximate expressions – Orders of magnitude**

Mesh stiffness can be determined using the Finite Elements Method but it is interesting to have orders of magnitude or approximate values at the design stage. For solid gears made of steel, an order of magnitude of the mesh stiffness per unit of contact length 0 *<sup>k</sup>* is 1.3 10 10 N/m². More accurate expressions can be derived from the ISO 6336 standard which, for solid gears, gives:

$$k\_o \equiv \cos \beta \frac{0.8}{q} \tag{27}$$

On the Modelling of Spur and Helical Gear Dynamic Behaviour 91

where **V M** is the extended structural vector: structural vector completed by zeros to the

<sup>0</sup> 0 1 *Tm m*

The separation of the average and time-varying contributions in the mesh stiffness function

<sup>0</sup> 1 0 0 , , 1,2 *<sup>T</sup>*

01 2 **XX X X** 

which, when re-injected into (33) and after identifying like order terms leads to the

 0 0 00 0 , , 1,2 *<sup>T</sup> <sup>m</sup> k t eM t*

0 00 00 1

Interestingly, the left-hand sides of all the differential systems are identical and the analysis of the eigenvalues and corresponding eigenvectors of the homogeneous systems will provide useful information on the dynamic behaviour of the geared systems under

The following system is considered (the influence of damping on critical speeds being

from which the eigenvalues and eigenvectors are determined. The technical problems associated with the solution of (36) are not examined here and the reader may refer to

orthogonal eigenvectors **Φ<sup>P</sup>** have been determined which, to a great extent, control the gear

specialised textbooks. It is further assumed that a set of real eigenvalues

*m m k kt*

 

2

*M dt*

<sup>0</sup> 00 0 0 0

*k M dM k t* **Kt K VV K VV** (32)

*t eM t*

**MX K VV X F F 0 1 F2** (35-1)

*T T*

**MX K VV X VV X** (35-2)

<sup>0</sup> *<sup>m</sup> <sup>k</sup>* **<sup>T</sup> MX K VV X 0 0 0** (36)

*<sup>p</sup>* and real

**MX K VV X F F 0 1 F2** (33)

*<sup>T</sup>* **V V** (31)

*T T m*

) and an asymptotic expansion of the

... (34)

total number of DOF of the model

(30) can be simplified as:

as *kt k t <sup>m</sup>* 1

For most gears,

Main order:

*th* order:

ignored):

set dynamic behaviour.

Using an averaged structural vector as in (18):

 

*<sup>m</sup> k t* 

solution can be sought as a straightforward expansion of the form:

following series of constant coefficient differential systems:

consideration (critical speeds, modeshapes).

is usually a small parameter ( 1

*L t*

leads to the following state equations:

with

: helix angle (on pitch cylinder)

$$q = \mathbf{C}\_1 + \frac{\mathbf{C}\_2}{Zn\_1} + \frac{\mathbf{C}\_3}{Zn\_2} + \mathbf{C}\_4\mathbf{x}\_1 + \mathbf{C}\_5\frac{\mathbf{x}\_1}{Zn\_1} + \mathbf{C}\_6\mathbf{x}\_2 + \mathbf{C}\_7\frac{\mathbf{x}\_2}{Zn\_2} + \mathbf{C}\_8\mathbf{x}\_1^2 + \mathbf{C}\_9\mathbf{x}\_2^2$$

coefficients 1 9 *C C* ,..., have been tabulated and are listed in Table 1 below

<sup>3</sup> cos *i i <sup>Z</sup> Zn* , *<sup>i</sup>* 1,2 are the number of teeth of the equivalent virtual spur pinion ( 1 *<sup>i</sup>* ) and gear ( 2 *i* ).

*<sup>i</sup> x* , *i* 1,2 , are the profile shift coefficients on pinion( 1 *i* ) and gear( 2 *i* )


Table 1. Tabulated coefficients for mesh stiffness calculations according to ISO 6336.

### **5. Equations of motion – Dynamic behaviour**

#### **5.1 Differential system**

The equations of motion for undamped systems are derived by assembling all the elemental matrices and forcing term vectors associated with the gears but also the supporting members (shafts, bearings, casing, etc.) leading to a parametrically excited non-linear differential system of the form:

$$\left[\mathbf{M}\right]\ddot{\mathbf{X}} + \left[\mathbf{K}(t,\mathbf{X})\right]\mathbf{X} = \mathbf{F}\_0 + \mathbf{F}\_1\left(t,\mathbf{X},\delta e\left(\mathbf{M}\right)\right) + \mathbf{F}\_2\left(t,\dot{\Omega}\_{1,2}\right) \tag{28}$$

where **X** is the total DOF vector, **M** and *t*, **K X** are the global mass and stiffness matrices. Note that, because of the contact conditions between the teeth, the stiffness matrix can be non-linear (partial or total contact losses may occur depending on shape deviations and speed regimes). **F0** comprises the constant nominal torques; **F X <sup>1</sup>** *t eM* , , includes the contributions of shape deviations (errors, shape modifications, etc.); **F2** *t*,1,2 represents the inertial effects due to unsteady rotational speeds

#### **5.2 Linear behaviour - Modal analysis**

Considering linear (or quasi-linear) behaviour, the differential system can be re-written as:

$$\mathbf{E}\left[\mathbf{M}\right]\ddot{\mathbf{X}} + \left[\mathbf{K}(t)\right]\mathbf{X} = \mathbf{F}\_0(t) + \mathbf{F}\_1\left(t, \delta e\left(M\right)\right) + \mathbf{F}\_2\left(t, \dot{\Omega}\_{1,2}\right) \tag{29}$$

The time variations in the stiffness matrix *t* **K** are caused by the meshing and, using the formulation based on structural vectors, the constant and time-varying components can be separated as:

$$\mathbb{E}\left[\mathbf{K}(\mathbf{t})\right] = \left[\mathbf{K}\_0\right] + \int\limits\_{l(\ast)} k\left(M\right)\overline{\mathbf{V}}\left(\mathbf{M}\right)\overline{\mathbf{V}}\left(\mathbf{M}\right)^T dM\tag{30}$$

where **V M** is the extended structural vector: structural vector completed by zeros to the total number of DOF of the model

Using an averaged structural vector as in (18):

$$\overline{\mathbf{V}}\_{0} = \frac{1}{T\_{m}} \int\_{0}^{T\_{n}} \overline{\mathbf{V}}(M) \, dt \tag{31}$$

(30) can be simplified as:

90 Mechanical Engineering

coefficients 1 9 *C C* ,..., have been tabulated and are listed in Table 1 below

*<sup>i</sup> x* , *i* 1,2 , are the profile shift coefficients on pinion( 1 *i* ) and gear( 2 *i* )

**5. Equations of motion – Dynamic behaviour** 

the inertial effects due to unsteady rotational speeds

**5.2 Linear behaviour - Modal analysis** 

separated as:

2 3 1 2 2 2 1 41 5 62 7 81 92 1 2 1 2 *CC x x q C Cx C Cx C Cx Cx Zn Zn Zn Zn* 

, *<sup>i</sup>* 1,2 are the number of teeth of the equivalent virtual spur pinion ( 1 *<sup>i</sup>* )

*C*<sup>1</sup> *C*<sup>2</sup> *C*<sup>3</sup> *C*<sup>4</sup> *C*<sup>5</sup> *C*<sup>6</sup> *C*<sup>7</sup> *C*<sup>8</sup> *C*<sup>9</sup> 0,04723 0,15551 0,25791 -0,00635 -0,11654 -0,00193 -0,24188 0,00529 0,00182

The equations of motion for undamped systems are derived by assembling all the elemental matrices and forcing term vectors associated with the gears but also the supporting members (shafts, bearings, casing, etc.) leading to a parametrically excited non-linear

where **X** is the total DOF vector, **M** and *t*, **K X** are the global mass and stiffness matrices. Note that, because of the contact conditions between the teeth, the stiffness matrix can be non-linear (partial or total contact losses may occur depending on shape deviations

contributions of shape deviations (errors, shape modifications, etc.); **F2** *t*,1,2 represents

Considering linear (or quasi-linear) behaviour, the differential system can be re-written as:

The time variations in the stiffness matrix *t* **K** are caused by the meshing and, using the formulation based on structural vectors, the constant and time-varying components can be

> 

0

*L*

1,2 **MX K X F F 0 1 F2** (29)

*k M dM*

*T*

**Kt K VMVM** (30)

*t t t eM t* , ,

, 1,2 **MX K X X F F X F 0 1 <sup>2</sup>** (28)

includes the

*t*, *t eM t* , ,

and speed regimes). **F0** comprises the constant nominal torques; **F X <sup>1</sup>** *t eM* , ,

Table 1. Tabulated coefficients for mesh stiffness calculations according to ISO 6336.

with 

> <sup>3</sup> cos *i*

and gear ( 2 *i* ).

**5.1 Differential system** 

differential system of the form:

*i <sup>Z</sup> Zn*

: helix angle (on pitch cylinder)

$$\mathbb{E}\left[\mathbf{K}\left(\mathbf{t}\right)\right] = \left[\mathbf{K}\_0\right] + \int\limits\_{\boldsymbol{L}\left(t\right)} k\left(\boldsymbol{M}\right) d\boldsymbol{M} \,\overline{\mathbf{V}}\_0 \,\overline{\mathbf{V}}\_0^{\boldsymbol{T}} = \left[\mathbf{K}\_0\right] + k\_m \left(t\right) \overline{\mathbf{V}}\_0 \,\overline{\mathbf{V}}\_0^{\boldsymbol{T}}\tag{32}$$

The separation of the average and time-varying contributions in the mesh stiffness function as *kt k t <sup>m</sup>* 1 leads to the following state equations:

$$\mathbf{E}\left[\mathbf{M}\right]\ddot{\mathbf{X}} + \left[\mathbf{K}\_0\right] + k\_m \left(1 + \alpha \wp(t)\right) \mathbf{\overline{V}}\_0 \cdot \mathbf{\overline{V}}\_0^\top \begin{bmatrix} \mathbf{X} = \mathbf{F}\_0 + \mathbf{F}\_1 \left(t, \delta e\left(M\right)\right) + \mathbf{F}\_2 \left(t, \dot{\Omega}\_{1,2}\right) \end{bmatrix} \tag{33}$$

For most gears, is usually a small parameter ( 1 ) and an asymptotic expansion of the solution can be sought as a straightforward expansion of the form:

$$\mathbf{X} = \mathbf{X}\_0 + a\mathbf{X}\_1 + a^2\mathbf{X}\_2 + \dots \tag{34}$$

which, when re-injected into (33) and after identifying like order terms leads to the following series of constant coefficient differential systems:

Main order:

$$\mathbf{E}\left[\mathbf{M}\right]\ddot{\mathbf{X}}\_0 + \left[\left[\mathbf{K}\_0\right] + k\_m \overline{\mathbf{V}}\_0 \overline{\mathbf{V}}\_0^T\right] \mathbf{X}\_0 = \mathbf{F}\_0 + \mathbf{F}\_1\left(t, \delta e\left(M\right)\right) + \mathbf{F}\_2\left(t, \dot{\Omega}\_{1,2}\right) \tag{35-1}$$

*th* order:

$$\mathbb{E}\left[\mathbf{M}\right]\ddot{\mathbf{X}}\_{\ell} + \left[\left[\mathbf{K}\_{0}\right] + k\_{m}\overline{\mathbf{V}}\_{0}\,\overline{\mathbf{V}}\_{0}^{T}\right]\mathbf{X}\_{\ell} = -k\_{m}\boldsymbol{\rho}(\mathbf{t})\,\overline{\mathbf{V}}\_{0}\,\overline{\mathbf{V}}\_{0}^{T}\mathbf{X}\_{\ell-1}\tag{35-2}$$

Interestingly, the left-hand sides of all the differential systems are identical and the analysis of the eigenvalues and corresponding eigenvectors of the homogeneous systems will provide useful information on the dynamic behaviour of the geared systems under consideration (critical speeds, modeshapes).

The following system is considered (the influence of damping on critical speeds being ignored):

$$\mathbb{E}\left[\mathbf{M}\right]\ddot{\mathbf{X}}\_{\ell} + \left[\mathbb{E}\mathbf{K}\_{0}\right] + k\_{w}\overline{\mathbf{V}}\_{\mathbf{0}}\,\overline{\mathbf{V}}\_{\mathbf{0}}^{\mathrm{T}}\Big]\mathbf{X}\_{\ell} = \mathbf{0} \tag{36}$$

from which the eigenvalues and eigenvectors are determined. The technical problems associated with the solution of (36) are not examined here and the reader may refer to specialised textbooks. It is further assumed that a set of real eigenvalues*<sup>p</sup>* and real orthogonal eigenvectors **Φ<sup>P</sup>** have been determined which, to a great extent, control the gear set dynamic behaviour.

On the Modelling of Spur and Helical Gear Dynamic Behaviour 93

a. the assumption of proportional damping (Rayleigh's damping) which, in this case,

0 00 *<sup>T</sup>*

The damping matrix is supposed to be orthogonal with respect to the mode-shapes of the

Following Graig (1981), the damping matrix can be deduced by a truncated summation on a

Regardless of the technique employed, it should be stressed that both (41) and (42) depend

with significant percentages of strain energy in the meshing teeth. The methods also rely on the assumption of orthogonal mode shapes which is realistic when the modal density (number of modes per frequency range) is moderate so that inter-modal couplings can be

Based on the previous developments, the linear response of gears to mesh parametric excitations can be qualitatively assessed. Response peaks are to be expected at all tooth critical speeds and every sub-harmonic of these critical speeds because mesh stiffness time

*Nr* 2 *<sup>T</sup> p p*

*p p*

**C M <sup>Φ</sup> <sup>M</sup> <sup>Φ</sup>** (42)

corresponds to the range of variation for modes

 2 *<sup>T</sup> <sup>P</sup> p p* **Φp p C Φ** 

with: *a b*, , two constants to be adjusted from experimental results

b. the use of (a limited number of) modal damping factors *<sup>p</sup>*

undamped system with the averaged stiffness matrix such that:

 : modal damping factor associated with mode *p* , *p p k m* : modal stiffness and mass associated with mode *p*

> **C<sup>Φ</sup>** *diag k m* 2

> > 1

*<sup>p</sup> m <sup>p</sup>* 

or introducing the modal damping matrix **C<sup>Φ</sup>** :

limited number of modes *Nr* leading to the formula:

on estimated or measured modal damping factors *<sup>P</sup>*

rather sparse. It seems that 0.02 0.1 *<sup>P</sup>*

*<sup>m</sup> ab k* **C M K VV** (40)

:

*k m* (41-1)

0 *<sup>T</sup>* **Φp q C Φ** (41-2)

*<sup>P</sup> p p* , *p N* 1, mod (41-3)

for which the data in the literature is

leads to:

with: *<sup>p</sup>* 

with: *<sup>p</sup> p*

neglected.

**5.3.2 Linear response** 

*p*

*k m*

Focusing on dynamic tooth loads, it is interesting to introduce the percentage of modal strain energy stored in the gear mesh which, for a given pair (*<sup>p</sup>* , **Φ<sup>P</sup>** ), is defined as:

$$\rho\_p = k\_m \frac{\boldsymbol{\Phi}\_\mathbf{p}^T \overline{\mathbf{V}}\_0 \overline{\mathbf{V}}\_0^T \boldsymbol{\Phi}\_\mathbf{p}}{\boldsymbol{\Phi}\_\mathbf{p}^T \Big[\!\!\left[\!\!K\_0\!\right] + k\_m \overline{\mathbf{V}}\_0 \overline{\mathbf{V}}\_0^T \right] \boldsymbol{\Phi}\_\mathbf{p}} = \boldsymbol{\upsilon}\_{\boldsymbol{\Phi}\_\mathbf{p}^p}^2 \frac{k\_m}{k\_{\boldsymbol{\Phi}\_\mathbf{p}^p}}\tag{37}$$

with *T T <sup>p</sup> v* **Φp p V V Φ**

, *<sup>m</sup> <sup>p</sup> k k* : average mesh stiffness and modal stiffness associated with (*<sup>p</sup>* , **Φ<sup>p</sup>** )

It as been shown (Velex & Berthe., 1989) that *<sup>p</sup>* is a reliable indicator of the severity of one frequency with regard to the pinion-gear mesh and it can be used to identify the potentially critical speeds *<sup>p</sup>* for tooth loading which are those with the largest percentages of modal strain energy in the tooth mesh. If the only excitations are those generated by the meshing (the mesh frequency is *Z*1 1 ), the tooth critical speeds can be expressed in terms of pinion speed as:

$$\Omega\_1 = \alpha\_p \nmid kZ\_1 \qquad k = 1, 2, \dots \tag{38}$$

Based on the contact length variations and on the transmission error spectrum, the relative severity of the excitations can be anticipated.

*Remark:* The critical frequencies are supposed to be constant over the speed range (gyroscopic effects are neglected). Note that some variations can appear with the evolution of the torque versus speed (a change in the torque or load can modify the average mesh stiffness especially for modified teeth).

For the one DOF tosional model in Figure 4, there is a single critical frequency *<sup>m</sup> k m* whose expression can be developed for solid gears of identical face width leading to:

$$
\Omega\_1 \cong \frac{\Lambda}{k} \frac{\cos a\_p}{M \mathbb{Z}\_1^2} \sqrt{\frac{b}{B}} \sqrt{\cos \beta\_b} \sqrt{\varepsilon\_a} \sqrt{1 + u^2} \tag{39}
$$

where 1,2,... *<sup>k</sup>* represents the harmonic order; 0 <sup>8</sup>*<sup>k</sup>* ( is the density), for steel gears

<sup>3</sup> 210 <sup>1</sup> *ms* ; *M* is the module (in meter); *B* is the pinion or gear thickness (supposed identical); *b* is the effective contact width (which can be shorter than *B* because of chamfers for example); 1 2 *Z u <sup>Z</sup>* , speed ratio.

#### **5.3 Dynamic response**

#### **5.3.1 The problem of damping**

Energy dissipation is present in all geared systems and the amount of damping largely controls the amplification at critical speeds. Unfortunately, the prediction of damping is still a challenge and, most of the time; it is adjusted in order to fit with experimental evidence. Two classical procedures are frequently employed:

a. the assumption of proportional damping (Rayleigh's damping) which, in this case, leads to:

$$\mathbf{a}\begin{bmatrix} \mathbf{C} \end{bmatrix} = a\begin{bmatrix} \mathbf{M} \end{bmatrix} + b\begin{bmatrix} \begin{bmatrix} \mathbf{K}\_0 \end{bmatrix} + k\_m \overline{\mathbf{V}}\_0 \ \overline{\mathbf{V}}\_0^{\top} \end{bmatrix} \tag{40}$$

with: *a b*, , two constants to be adjusted from experimental results

b. the use of (a limited number of) modal damping factors *<sup>p</sup>* :

The damping matrix is supposed to be orthogonal with respect to the mode-shapes of the undamped system with the averaged stiffness matrix such that:

$$\mathbf{D}\_{\mathbf{p}}^{\top} \begin{bmatrix} \mathbf{C} \end{bmatrix} \boldsymbol{\Phi}\_{\mathbf{p}} = 2 \boldsymbol{\zeta}\_{\mathcal{P}} \sqrt{k\_{\boldsymbol{\Phi} \boldsymbol{\eta}} m\_{\boldsymbol{\Phi} \boldsymbol{\eta}}} \tag{41.1}$$

$$\mathbf{O}\_p^\top \begin{bmatrix} \mathbf{C} \end{bmatrix} \boldsymbol{\Phi}\_q = \mathbf{0} \tag{41.2}$$

with: *<sup>p</sup>* : modal damping factor associated with mode *p*

, *p p k m* : modal stiffness and mass associated with mode *p*

or introducing the modal damping matrix **C<sup>Φ</sup>** :

$$\left[\mathbf{C}\_{\Phi}\right] = \text{diag}\left(2\,\varphi\_{p}\sqrt{k\_{\Phi p}m\_{\Phi p}}\right), \; p = 1, N \bmod \tag{41.3}$$

Following Graig (1981), the damping matrix can be deduced by a truncated summation on a limited number of modes *Nr* leading to the formula:

$$\begin{bmatrix} \mathbf{C} \end{bmatrix} = \sum\_{p=1}^{\text{N}} \frac{2\boldsymbol{\varepsilon}\_{p}\,\boldsymbol{\alpha}\_{p}}{m\_{\boldsymbol{\Phi}\boldsymbol{\eta}}} \left( \begin{bmatrix} \mathbf{M} \end{bmatrix} \boldsymbol{\Phi}\_{p} \right) \left( \begin{bmatrix} \mathbf{M} \end{bmatrix} \boldsymbol{\Phi}\_{p} \right)^{T} \tag{42}$$

with: *<sup>p</sup> p p k m* 

92 Mechanical Engineering

Focusing on dynamic tooth loads, it is interesting to introduce the percentage of modal

0 00

frequency with regard to the pinion-gear mesh and it can be used to identify the potentially

strain energy in the tooth mesh. If the only excitations are those generated by the meshing (the mesh frequency is *Z*1 1 ), the tooth critical speeds can be expressed in terms of pinion

1 1 / 1,2,... *<sup>p</sup>*

Based on the contact length variations and on the transmission error spectrum, the relative

*Remark:* The critical frequencies are supposed to be constant over the speed range (gyroscopic effects are neglected). Note that some variations can appear with the evolution of the torque versus speed (a change in the torque or load can modify the average mesh

For the one DOF tosional model in Figure 4, there is a single critical frequency *<sup>m</sup>*

*b*

<sup>3</sup> 210 <sup>1</sup> *ms* ; *M* is the module (in meter); *B* is the pinion or gear thickness (supposed identical); *b* is the effective contact width (which can be shorter than *B* because of chamfers

Energy dissipation is present in all geared systems and the amount of damping largely controls the amplification at critical speeds. Unfortunately, the prediction of damping is still a challenge and, most of the time; it is adjusted in order to fit with experimental evidence.

cos 1 *<sup>p</sup> b*

 

(

whose expression can be developed for solid gears of identical face width leading to:

1

*k MZ B*

1 2

where 1,2,... *<sup>k</sup>* represents the harmonic order; 0 <sup>8</sup>*<sup>k</sup>*

cos

*<sup>k</sup> k v*

*T T*

**Φ V V Φ**

*p m T T p*

 **p p p p**

0 0 2

*m p*

*<sup>p</sup>* for tooth loading which are those with the largest percentages of modal

*k k*

*m*

*kZ k* (38)

2

*u*

(39)

*k m*

is the density), for steel gears

*<sup>p</sup>* , **Φ<sup>p</sup>** )

*<sup>p</sup>* is a reliable indicator of the severity of one

**<sup>Φ</sup> K VV <sup>Φ</sup>** (37)

*<sup>p</sup>* , **Φ<sup>P</sup>** ), is defined as:

strain energy stored in the gear mesh which, for a given pair (

, *<sup>m</sup> <sup>p</sup> k k* : average mesh stiffness and modal stiffness associated with (

with *T T*

critical speeds

speed as:

*<sup>p</sup> v* **Φp p V V Φ**

It as been shown (Velex & Berthe., 1989) that

severity of the excitations can be anticipated.

stiffness especially for modified teeth).

2 *Z u*

*<sup>Z</sup>* , speed ratio.

Two classical procedures are frequently employed:

for example); 1

**5.3 Dynamic response** 

**5.3.1 The problem of damping** 

Regardless of the technique employed, it should be stressed that both (41) and (42) depend on estimated or measured modal damping factors *<sup>P</sup>* for which the data in the literature is rather sparse. It seems that 0.02 0.1 *<sup>P</sup>* corresponds to the range of variation for modes with significant percentages of strain energy in the meshing teeth. The methods also rely on the assumption of orthogonal mode shapes which is realistic when the modal density (number of modes per frequency range) is moderate so that inter-modal couplings can be neglected.

#### **5.3.2 Linear response**

Based on the previous developments, the linear response of gears to mesh parametric excitations can be qualitatively assessed. Response peaks are to be expected at all tooth critical speeds and every sub-harmonic of these critical speeds because mesh stiffness time

On the Modelling of Spur and Helical Gear Dynamic Behaviour 95

introducing the unit Heaviside step function *H x* such that *H x* 1 if 0 *x* and

amplitudes of tooth modifications reducing the actual contact patterns, to spalls on the

b. the amplitude of the dynamic displacement **X** is sufficiently large so that the teeth can

Momentary contact losses can therefore occur when vibration amplitudes are sufficiently large; they are followed by a sequence of free flight within the tooth clearance until the teeth collide either on the driving flanks or on the back of the teeth (back strike). Such shocks are particularly noisy (rattle noise) and should be avoided whenever possible. Analytical investigations are possible using harmonic balance methods and approximations of *H*(*x*) (Singh et al., 1989), (Comparin & Singh, 1989), (Kahraman & Singh, 1990), (Kahraman & Singh, 1991), and numerical integrations can be performed by time-step schemes (Runge-

a. contact losses move the tooth critical frequencies towards the lower speeds (softening effect) which means that predictions based on a purely linear approach might be irrelevant. The phenomenon can be observed in Fig. 11 where the experimental peaks

Fig. 12. Dynamic response curves by numerical simulations – Amplitude jumps – Influence of the initial conditions (speed up vs speed down), *2% of the critical damping, Spur gears* 

It can observed that contact losses are related to the sign of <sup>0</sup> *M <sup>T</sup>* **V X**

separate ( **X** is periodic and can become negative in some part of the cycle).

which, it can be deduced that two cases have to be considered:

*e M* is larger than the normal approach 0

flanks (pitting) where contact can be lost, etc.

Kutta, Newmark, etc.). The most important conclusions are:

<sup>1</sup> *Z* 30 *,* <sup>2</sup> *Z* 45 *, M* 2*mm , standard tooth proportions.* 

are at lower speed than those predicted by the linear theory.

**dF M 1/2** *k M M H M dM* **n1** (44)

*<sup>T</sup>* **V X** which, typically, corresponds to large

*e M* , from

*H x* 0 otherwise. Finally, one obtains:

a.

variations may exhibit several harmonics with significant amplitudes. Figure 11, taken from Cai and Hayashi (1994), is a clear example of such typical dynamic response curves when the gear dynamic behaviour is dominated by one major tooth frequency*<sup>n</sup>* (and can be simulated by using the classic one DOF model). The amplifications associated with each peak depends on i) the excitation amplitude (Eq. (27) can provide some information on the amplitude associated with each mesh frequency harmonic) and ii) the level of damping for this frequency. For more complex gear sets, interactions between several frequencies can happen but, as far as the author is aware, the number of frequencies exhibiting a significant percentage of modal strain energy in the tooth mesh seems very limited (frequently less than 5) thus making it possible to anticipate the potential dangerous frequency coincidences for tooth durability.

Fig. 11. Examples of dynamic response curves (Cay & Hayashi, 1994).

#### **5.3.3 Contact condition – Contact losses and shocks**

Only compressive contact forces can exist on tooth flanks and using (11), this imposes the following unilateral condition in case of contact at point *M*:

$$\mathbf{dF\_{1/2}}\left(\mathbf{M}\right)\cdot\mathbf{n\_1} = k\left(M\right)\Delta\left(M\right)dM > 0\tag{43}$$

or, more simply, a positive mesh deflection *M* .

If *M* 0 , the contact at M is lost (permanently or temporarily) and the associated contact force is nil. These constraints can be incorporated in the contact force expression by

variations may exhibit several harmonics with significant amplitudes. Figure 11, taken from Cai and Hayashi (1994), is a clear example of such typical dynamic response curves when

simulated by using the classic one DOF model). The amplifications associated with each peak depends on i) the excitation amplitude (Eq. (27) can provide some information on the amplitude associated with each mesh frequency harmonic) and ii) the level of damping for this frequency. For more complex gear sets, interactions between several frequencies can happen but, as far as the author is aware, the number of frequencies exhibiting a significant percentage of modal strain energy in the tooth mesh seems very limited (frequently less than 5) thus making it possible to anticipate the potential dangerous frequency coincidences

*<sup>n</sup>* (and can be

the gear dynamic behaviour is dominated by one major tooth frequency

Fig. 11. Examples of dynamic response curves (Cay & Hayashi, 1994).

Only compressive contact forces can exist on tooth flanks and using (11), this imposes the

If *M* 0 , the contact at M is lost (permanently or temporarily) and the associated contact force is nil. These constraints can be incorporated in the contact force expression by

*k M M dM* 0 **dF M n 1/2 <sup>1</sup>** (43)

**5.3.3 Contact condition – Contact losses and shocks** 

following unilateral condition in case of contact at point *M*:

or, more simply, a positive mesh deflection *M* .

for tooth durability.

introducing the unit Heaviside step function *H x* such that *H x* 1 if 0 *x* and *H x* 0 otherwise. Finally, one obtains:

$$\mathbf{dF\_{1/2}}(\mathbf{M}) = k(M)\Delta(M)H(\Delta(M))dM\mathbf{n\_1} \tag{44}$$

It can observed that contact losses are related to the sign of <sup>0</sup> *M <sup>T</sup>* **V X** *e M* , from which, it can be deduced that two cases have to be considered:


Momentary contact losses can therefore occur when vibration amplitudes are sufficiently large; they are followed by a sequence of free flight within the tooth clearance until the teeth collide either on the driving flanks or on the back of the teeth (back strike). Such shocks are particularly noisy (rattle noise) and should be avoided whenever possible. Analytical investigations are possible using harmonic balance methods and approximations of *H*(*x*) (Singh et al., 1989), (Comparin & Singh, 1989), (Kahraman & Singh, 1990), (Kahraman & Singh, 1991), and numerical integrations can be performed by time-step schemes (Runge-Kutta, Newmark, etc.). The most important conclusions are:

a. contact losses move the tooth critical frequencies towards the lower speeds (softening effect) which means that predictions based on a purely linear approach might be irrelevant. The phenomenon can be observed in Fig. 11 where the experimental peaks are at lower speed than those predicted by the linear theory.

Fig. 12. Dynamic response curves by numerical simulations – Amplitude jumps – Influence of the initial conditions (speed up vs speed down), *2% of the critical damping, Spur gears*  <sup>1</sup> *Z* 30 *,* <sup>2</sup> *Z* 45 *, M* 2*mm , standard tooth proportions.* 

On the Modelling of Spur and Helical Gear Dynamic Behaviour 97

Fig. 14. Examples of quasi-static T.E. measurements and simulations – Spur gear

Fig. 15. T.E. measurements at various loads – Helical gear example.

**6.2 No-load transmission error (NLTE)** 

**6.3 Transmission errors under load** 

NASA measurements from www.grc.nasa.gov/WWW/RT2001/5000/5950oswald1.html.

No-load *T.E*. (*NLTE*) has already been introduced in (2); it can be linked to the results of gear testing equipment (single flank gear tester) and is representative of geometrical deviations. From a mathematical point of view, *NLTE* is derived by integrating (2) and is expressed as:

The concept of transmission error under load (TE) is clear when using the classic single degree of freedom torsional model (as Harris did) since it directly relies on the angles of

*NLTE*

(45)

cos *MAX b*

*E t*

(Velex and Maatar, 1996).


These phenomena are illustrated in the response curves in Figure 12.
