Behavior of lubricant fluid film in gears under dynamic conditions

Tribology International 62 (2013) 37–48

Contents lists available at SciVerse ScienceDirect

Tribology International journal homepage: www.elsevier.com/locate/triboint

Behavior of lubricant fluid film in gears under dynamic conditions Marco Barbieri a,n, Antonius A. Lubrecht b, Francesco Pellicano a a b

Universita di Modena e Reggio Emilia, Dipartimento di Ingegneria Enzo Ferrari, ViaVignolese 905, 41125 Modena, Italy Laboratoire de Mecanique des Contacts et des Solides, UMR-CNRS 5514, INSA de Lyon, Lyon, France

a r t i c l e i n f o

abstract

Article history: Received 26 June 2012 Received in revised form 14 January 2013 Accepted 16 January 2013 Available online 1 February 2013

This paper presents a new method for modeling the fluid film lubrication in gears, considering the actual meshing conditions and gear dynamics. The model takes into account both the elastohydrodynamic lubrication (EHL) and the dynamic load between the mating tooth pair. The EHL film is described as a fully flooded elliptical contact. The present approach is validated by means of comparison with other methods found in the literature, in which dynamic effects were neglected. The effect of the gear dynamics on the fluid film is investigated. It is shown that pressure and film thickness are strongly modified by the dynamics of the gear pair. The dependence of the dynamic gear lubrication on dimensionless parameters is investigated: a new dimensionless inertia parameter is added to the standard Moes’ parameters. These parameters are useful to describe the lubrication conditions in gear pairs. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Transient EHL Gear dynamics Multigrid

1. Introduction Understanding the dynamic behavior of gear pairs is crucial in order to reduce vibration and noise, as well as to increase the lifetime of gear systems, and this is a fortiori true for spur gear pairs. For this reason, many approaches have been proposed in the past years to describe the dynamics of gear systems. In 1958, Harris [1] proposed a method to reduce vibrations in spur gear pairs by means of profile reliefs. Harris assumed that the coupling is stiffer when two pairs of teeth are in contact, more deformable when only one couple is mating. In 1990, Kahraman et al. [2] described the dynamics of a spur gear pair using a 1-dof lumped parameter model. Such model was validated by Kahraman and Blankenship [3], who also proved the effectiveness of profile modifications in reducing gear vibrations. In the last years, many models for gear dynamics were proposed, taking into account misalignments, profile errors or their complicating effects; however, few efforts were done to understand the dissipation mechanism. Since the contact between the teeth is lubricated, the presence of the lubricant can play an important role in the dynamics of the gear pair, as well as the dynamics of the gear pair will affect the lubrication condition: dynamics and lubrication are coupled problems. The lubrication regime for industrial gear systems is generally elastohydrodynamic, i.e. it is necessary to keep into account the

n

Corresponding author. Tel.: þ39 059 2056312; fax: þ39 059 2056126 E-mail addresses: [email protected] (M. Barbieri), [email protected] (A.A. Lubrecht), [email protected] (F. Pellicano). 0301-679X/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.triboint.2013.01.017

local deformation of the mating surfaces due to lubricant pressure, in order to obtain a proper estimation of the oil pressure and the minimum film thickness. Indeed, since the first results due to Grubin [4] in 1949, the coupling between the Reynolds equation and the elastic deformation was considered in gears. The minimum film thickness obtained from EHL computation was high enough to explain the low wear rate measured on tooth flanks. More recently, the problem of transient EHL lubrication in gear pairs has been faced by many authors. In fact, gear contacts are not steady from a lubrication point of view, even if gears are operating in steady conditions. In 1997, Larsson [5] presented a solution for the lubricated line contact problem under a square-wave load, to simulate the varying number of teeth in contact. In 2004, Wang et al. [6] used a trapezoidal varying load and took into account thermal effects and non-Newtonian lubricants. In 2009, Evans et al. developed a method to model the mixed lubrication conditions in gear pairs [7]. All these approaches neglect the inertial effects in the meshing gears, so they are good approximations of the real behavior only for low rotational speeds. To point out the role of gear dynamics on the contact, recently Osman and Velex [8] have developed a model for gear contact fatigue, considering dynamic tooth loads. In 2010, Li et al. [9] published a model for transient mixed lubrication in gear pairs. In the same year, some results on the coupling between EHL lubrication in gear pair and the gear dynamics were presented in [11]. Recently, Li et al. [12] presented some results on the effect of gear dynamics on gear lubrication, using two decoupled models for the gear dynamics and the EHL lubrication. In 2012, De La Cruz et al. [10] published an interesting tribo-dynamic model for a multimesh transmission: their model considers actual coupling between EHL and gear dynamics, but assumes constant stiffness in gear meshing.

38

M. Barbieri et al. / Tribology International 62 (2013) 37–48

Nomenclature

r a0 aB Z l x Rx

wph n

dimensionless density transversal pressure angle in working condition (rad) dimensionless pressure viscosity coefficient dimensionless dynamic viscosity dimensionless coefficient inside the Reynolds equation dimensionless coefficient inside the Reynolds equation dimensionless equivalent radius of curvature along the profile n parameter correcting the elliptic case pffiffiffiffiffiffiffiffiffiffiffiffiffiffi2ph for pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi n n2 n wpn ¼ ðpD Þ=ð4k Þ 1k =ðK ‘ ð 1kn2 ÞE‘ ð 1kn2 ÞÞ h

Z0 Z k wn w‘

ni Oi

r0 r t Wgi x n

a b Dn E0 E0i E‘ f Fg fm fn Fw H h H0

reference viscosity (Pa s) dynamic viscosity (Pa s) ellipticity ratio k ¼ a=b external load at the pitch point wn ¼ T g1 =Rb1 equivalent external line load (N/m) Poisson’s modulus for the ith gear angular velocity for the ith gear (rad/s) reference density (kg/m3) density (kg/m3) transmission ratio of the gear pair t ¼ O2 =O1 rotation angle of the ith gear (rad) contact position along the line of action (m) reference values at the pitch point in the stationary conditions minor half axis of the contact ellipse (m) major half axis of the contact ellipse (m) reference ratio of the radii at pitch point Dn ¼ Ry =Rnx  equivalent Young modulus E0 ¼ 2 ð1n21 Þ=E1 þ 1 ð1n22 Þ=E2 (Pa) Young modulus for the ith gear (Pa) complete elliptic integral of the second kind load sharing coefficient total contact force in gear model (N) mesh frequency of the gear pair f m ¼ O1 =Z 1 (Hz) purely natural frequency f n ¼ 1=ð2pÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi elastic kel =mge (Hz) face-width of the gear (m) dimensionless film thickness film thickness (m) dimensionless rigid body approach

The elastohydrodynamic lubrication problem can be solved by means of numerical integration. In the case of gear tooth contact, the conditions are by no means steady state, because the load, the relative curvature and the relative velocity change during gear meshing. To this extent, a transient EHL modeling is needed. Furthermore, usual load conditions in gear systems lead to contact pressures up to 2 GPa, and highly loaded contacts require a very accurate discretization of the EHL computational domain. For all these reasons, in order to speed up the computation and to have accurate results, the multilevel solver proposed by Venner and Lubrecht [13] is here adapted to solve the elliptical and transient problem. The purpose of the present work is to investigate the effect of higher gear speeds (close to the natural frequency of the gear pair) on the film thickness and the contact pressure distribution. As an extra result, the dynamic transmission error of the gear pair under lubricated condition can be computed. Some dimensionless parameters, useful to identify the dynamic lubrication regime in gear pairs, are proposed.

h0 In Igi K‘ kel Ln m Mn mge magc ni P p p0 pi Rx Ry Rai Rbi Rix Riy T t tn Tge Tgi ui Um um W X x xa ,xb xg Y y ya ,yb z Zi

rigid body approach (m) 2 dimensionless inertia parameter In ¼ mge =ðRnx t n pnh Þ inertia of the ith gear (kg m2) complete elliptic integral of the first kind elastic stiffness (N/m) Moes’ dimensionless parameter at the pitch point 1=4  Ln ¼ aB E0 ð2Z0 unm Þ=ðE0 Rnx Þ normal module (m) Moes’ dimensionless parameter at the pitch point 3=4 2  M n ¼ w=ðE0 Rnx Þ ð2Z0 unm Þ=ðE0 Rnx Þ equivalent mass of the gear pair (kg) magnitude of crowning (m) rotational speed for the ith gear (RPM) dimensionless pressure pressure (Pa) reference pressure (Pa) contact pressure on the ith contact (Pa) equivalent radius of curvature along the profile (m) equivalent radius of curvature along the facewidth (m) outer diameter for the ith gear (m) base radius of the ith gear (m) radius of curvature along the profile for the ith gear (m) radius of curvature along the facewidth for the ith gear (m) dimensionless time time (s) reference time (s) equivalent torque (N) torque applied to the ith gear (Nm) tangent component of velocity for the ith gear (m/s) dimensionless mean velocity mean velocity (m/s) transmitted power (W) dimensionless x coordinate coordinate of the contact ellipse along the profile limits for x (m) transmission error on the line of action (m) dimensionless y coordinate coordinate of the contact ellipse along the face width limits for y (m) pressure viscosity index z ¼ aB p0 =ðlnðZ0 Þ þ 9:67Þ number of teeth for the ith gear

2. A model for a gear pair In order to describe the dynamics of a lubricated gear pair, the dynamic equilibrium of each gear is imposed; the equation of motion is coupled to the transient Reynold’s equation and to the elastic deformation equation. Spur gear pairs in the presence of crowning are considered in this work: this means that the theoretical dry contact pattern is elliptical and the EHL point lubrication problem is to be solved.

2.1. Dynamic model of the gear pair Fig. 1 depicts the dynamic model of a gear pair with contact ratio between 1 and 2: the load is shared by one or two EHL contacts, provided that one or two pairs of teeth are in contact at a given time. If the total contact force is denoted with Fg, the

M. Barbieri et al. / Tribology International 62 (2013) 37–48

39

Fig. 1. Dynamic model of the gear pair.

equations of the gear system are the following: 8 < Ig1 W€ g1 ðtÞ ¼ T g1 Rb1 F g : I W€ ðtÞ ¼ T þ R F g g2

g2

g2

ð1Þ

b2

where the external torques are assumed constant, and Fg is an unknown function F g ¼ F g ðW_ g1 , W_ g2 , Wg1 , Wg2 ,tÞ. This set of two equations can be reduced to one by introducing the transmission error on the line of action xg(t), together with the equivalent mass mge and the equivalent torque Tge   R xg ðtÞ ¼ Rb1 Wg1 ðtÞ b2 Wg2 ðtÞ Rb1 !1   2 2 Rb1 T g1 Rb2 T g2 T g1 Rb1 Rb2 þ , T ge ¼ mge þ ð2Þ mge ¼ ¼ Ig1 Ig2 Ig1 Ig2 Rb1 The resulting equation of motion reads mge x€ g ðtÞ þF g ðx_ g ,xg ,tÞ ¼ T ge

ð3Þ

The total load F g ðx_ g ,xg ,tÞ is given by the two simultaneous contact pressure distributions p1 ðx,y,tÞ and p2 ðx,y,tÞ, the transmission error is expressed as the approach h0 of the contacting bodies, with h0 ¼ xg ðtÞ; Eq. (3) can be written as Z 1Z 1 Z 1Z 1 2 T g1 d h0 ¼ p1 ðx,y,tÞ dx dy þ p2 ðx,y,tÞ dx dy mge 2 Rb1 1 1 1 1 dt ð4Þ where x and y indicate the position with respect to the axes of the contact ellipse, as shown in Fig. 2. Note that Ry b Rx , see Fig. 2. The major and the minor half-axes are b and a. Eq. (4) is the new dynamic model proposed here, it is new and more realistic than those present in the literature, but two fluid fields must be numerically determined at each time step. This model is called ‘‘dynamic model’’. The dynamic model can be simplified by using a coefficient f(t) in order to describe how the total load Fg(t) is shared by the two simultaneous contacts Z 1Z 1 2 T g1 d h0 1 ¼ f ðtÞ pðx,y,tÞ dx dy ð5Þ mge 2 Rb1 1 1 dt Of course, imposing f(t) a priori is somewhat arbitrary, this function is generally imposed in a heuristic way. Let us call this model ‘‘simplified dynamic model’’. Furthermore, one can neglect the inertial effects in the contacting bodies, i.e. Ig1 ¼ Ig2 ¼ e, e-0; this limiting case of Eq. (4) becomes Z 1Z 1 T g1 1 pðx,y,tÞ dx dy ¼ ð6Þ f ðtÞ Rb1 1 1 This model can be called ‘‘quasi static’’, as inertia is neglected, but transient EHL is considered. It is to note that Eq. (6) is the same model proposed in [5,6,14]. The load sharing between tooth pairs was simulated by Larsson [5] through a step function, whereas Wang [6] used a

Fig. 2. Contact ellipse over a tooth flank.

trapezoidal function; here both formulations will be considered for comparison. The step function is 1 when only one pair of teeth is in contact, 1/2 when two pairs of teeth are meshing at the same time. The trapezoidal function reduces the step, in order to keep into account tooth bending, as depicted in Fig. 3. It is to note that the trapezoidal formulation is considered in order to draw comparisons with [6]. 2.2. Lubrication model The fluid film lubrication between the contacting teeth is modeled considering two fully flooded point contacts. The EHL solution is isothermal and does not take into account the effect of the surface roughness. The transient Reynolds equation reads ! ! 3 3 @ rh @p @ rh @p @ @ ð7Þ þ um ðrhÞ ðrhÞ ¼ 0 @x 12Z @x @y 12Z @y @x @t where p, h, r, and Z are function of the position (x,y) and time t. The mean velocity um ¼ um ðtÞ varies along the contact; density– pressure dependence is described by means of the Dowson– Higginson relationship [15]

rðpÞ ¼ r0

5:9  108 Pa þ1:34p 5:9  108 Pa þp

ð8Þ

where the pressure p is expressed in Pascal. Viscosity–pressure dependence is modeled using Roelands’ formulation [16]     a p p z ZðpÞ ¼ Z0 exp B 0 1 þ 1 ð9Þ p0 z The equation of the film thickness, that describes the deformation in the contacting surfaces, is Z 1Z 1 0 0 x2 y2 2 pðx0 ,y0 ,tÞ dx dy qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð10Þ þ þ h ¼ h0 þ 0 2Rx 2Ry pE 1 1 ðxx0 Þ2 þðyy0 Þ2

M. Barbieri et al. / Tribology International 62 (2013) 37–48

1

Load sharing factor - f

Load sharing factor - f

40

0.8 0.6 0.4 0.2 0 -5

-4 -3 -2 -1 0 1 2 3 4 Position on the contact line - [mm]

1 0.8 0.6 0.4 0.2 0 -5

5

-4 -3 -2 -1 0 1 2 3 4 Position on the contact line - [mm]

5

Fig. 3. Different load sharing functions: (a) step function; (b) trapezoidal function.

and the equivalent radius, Eq. (11), becomes

Table 1 Size of the computational domain. M

n

M n o 200 M n Z 200 n

2

xa

xb

ya ¼ yb

4:5an

1:5an 1:25an

3:0an =kn 1:1an =kn

2:083an

Values at the pitch point.

Rx ðxÞ ¼

Rx ¼

R1x R2x , R1x þ R2x

Ry ¼

R1y R2y R1y þ R2y

pðxa ,y,tÞ ¼ 0

8y,t

pðxb ,y,tÞ ¼ 0 pðx,ya ,tÞ ¼ 0

8y,t 8x,t

Ry ¼

pðx,yb ,tÞ ¼ 0

8x,t

ð12Þ

Rb1 tanða0 Þ , ð1 þ tÞ

R x ðxÞ ¼

Rx ðxÞ Rnx

ð15Þ

where R x ðxÞ is a dimensionless equivalent radius. The radius of curvature in the y direction is constant; it depends on the value of crowning magc and on the face-width of the gear Fw

ð11Þ

The radius R1x is measured along the profile of the first gear (pinion), R1y can be obtained from the measure of the pinion crowning; the same definitions apply to R2x and R2y. The effect of cavitation is taken into account by checking if the pressure is positive over all the computational domain. The boundary conditions are

ð14Þ

where t is the transmission ratio of the gear pair and a0 is the actual transversal pressure angle. The reference value for Rx is the value at the pitch point Rnx ¼ Rx ð0Þ Rnx ¼

where Rx and Ry are the equivalent radii of curvature along the profile and the facewidth of the gear, see Fig. 2. They are related to the actual radii along the tooth profiles

tx þ ð1tÞ Rb1 tanða0 Þx þ R2b1 tan2 ða0 Þ Rb1 ð1þ tÞtanða0 Þ

Fw 8mag c

ð16Þ

In order to define the rolling/sliding velocities, it is necessary to compute the component of the velocity tangent to the profiles of pinion and gear u1 ðxÞ ¼ O1 R1x ðxÞ u2 ðxÞ ¼ O2 R2x ðxÞ ¼ tO1 R2x ðxÞ The mean velocity at the contact point is   u1 ðxÞ þ u2 ðxÞ 1t ¼ O1 Rb1 tanða0 Þ þ um ðxÞ ¼ x 2 2

ð17Þ

ð18Þ

The reference value at pitch point unm is

with x A ½xa ,xb , yA ½ya ,yb . The domain is set considering the Moes’ parameter Mn [17], in Table 1. This parameter is defined in the nomenclature; it will be discussed in the following. Note that the problem is time dependent, so that it is necessary to choose a reference configuration during gear meshing and freeze the dimension of the computational domain for that position. Here the pitch point is chosen as reference: an, bn are the half-axes of the contact ellipse at the pitch point. It is to note that in a spur gear pair in static conditions, the contact area is maximum close to the pitch point, where the load is shared by the minimum number of tooth pairs.

unm ¼ O1 Rb1 tan ða0 Þ 2.4. Dimensionless equations and parameters

In order to reduce the total number of parameters needed to solve the problem, it is useful to introduce the following non dimensional quantities x y Rn , Y ¼ n , H ¼ x2 h n a a an p Z r P¼ n, Z¼ , r¼ ph Z0 r0

X¼

Um ¼

2.3. Time-dependent parameters In a spur gear pair, the radii of curvature and the velocities of the contact point are functions of the contact position along the line of action. Denoting with x the coordinate along the line of action, Fig. 4, one can easily calculate R1x and R2x

ð19Þ

tn um , an

T¼

t tn

ð20Þ

The reference values an and pnh refer to the hertzian solution at the pitch point; Z0 , r0 refer to the environmental values, tn is defined as follows: tn ¼

an unm

ð21Þ

0

R1x ¼ Rb1 tanða Þ þ x R2x ¼ Rb2 tanða0 Þx

ð13Þ

This choice allows to have the same coefficient of the wedge and of the squeeze terms inside the Reynolds equation at the pitch

M. Barbieri et al. / Tribology International 62 (2013) 37–48

41

[mm]

Fig. 4. Line of action in a spur gear pair.

point. The dimensionless EHL problem reads      @ @ @P @ @P @  þ U m x x r H  ðr HÞ ¼ 0 @X @X @Y @Y @X @T Z 1Z 1 0 0 1 2 Dn 2 2 PðX 0 ,Y 0 ,TÞ dX dY qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Y þ 2 wpn X þ H ¼ H0 þ h 2 p 2R x 1 1 ðXX 0 Þ2 þ ðYY 0 Þ2 ð22Þ n

n

Where D ¼ Ry =Rx , wpn is a correction term for the elliptical h problem (see the nomenclature for details) and x is defined as follows:

x¼

r H3 lZ

where l ¼

12unm Z0 R2x an3 pnh

ð23Þ

As usual for EHL problems, it is useful to introduce a minimum set of dimensionless parameters: here Moes’ parameters are defined with reference to the stationary EHL elliptical problem at the pitch point; details about the solution of such problem can be found in Nijenbanning et al. [18]. Moes’ parameters for

this problem are !1=4   T g1 2Z0 unm 3=4 ð1 þ tÞ5 M ¼ ¼ 3 2 E0 Rnx Rb1 tan2 ða0 Þ 8Z30 E0 O31 E0 Rnx   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2Z0 unm 1=4 4 Ln ¼ aB E0 ¼ 2Z0 a4B E03 O1 ð1 þ tÞ 0 n E Rx n

wn

ð24Þ

where wn ¼ T g1 =Rb1 is the external load. Using the aforementioned transformations, the equation of motion (4) can be reduced to its dimensionless form: Z 1Z 1 2 d H0 0 0 In ¼ P 1 ðX 0 ,Y 0 ,TÞ dX dY 2 1 1 dT Z 1Z 1 2p 0 0 þ P2 ðX 0 ,Y 0 ,TÞ dX dY  n ð25Þ 3 k 1 1 where In is a dimensionless inertia In ¼

2

mge n n2

Rx t ph n

¼

2pmge unm 1 3Rnx wn kn

ð26Þ

42

M. Barbieri et al. / Tribology International 62 (2013) 37–48

Fig. 5. Flowchart of the method for a gear pair having 1 o ea o 2.

Mn, Ln, Dn and In are the four fundamental parameters describing the problem of gear lubrication. The problem is also slightly dependent on the laws of variation of the equivalent curvature and of the mean velocity at the contact point, as well as the pressure viscosity index z, if Roeland’s viscosity pressure relationship (9) is used. Fig. 5 clarifies how the presence of two contacts is taken into account. The first step (k¼ 0) is a stationary EHL solution for contact 1 at the pitch point, while the pressure on contact 2 is set to zero. The active contact at the pitch point is denoted by A, whereas B is the inactive contact. Within the first mesh period (m¼1), the contact 2 enters at time TscB1 and the contact 1 exits at time TecA1. At the end of the mesh period, the indexes A and B are reversed. The main issue in the code is the initialization of the pressure for contact B: since the solution method is iterative, an initial guess pressure has to be chosen. In highly loaded contacts, the guess pressure distribution must be very close to the final one in order to achieve convergence. In the present paper, a continuation method has been developed to define initial guess pressure. To this extent, Fig. 6 shows how the guess pressure is initialized: the pressure on contact B is set equal to the pressure on contact A at the previous step, the time varying parameters R x and Um on

contact B are set equal to the parameters on contact A and the stationary solution for the fixed rigid approach H0k is found for contact B. Later on, R x and Um are homotopically deformed up to their prescribed values in 10 steps, so that the resulting pressure is a suitable guess value for PBijk . Further details about the model can be in found in Ref. [19].

3. Validation of the static solver In order to validate the present method, a comparison with results of Wang et al. [6] has been carried out. In that paper, a transient line EHL solution for the static case, see Eq. (6), is presented, both for the thermal and for the iso-thermal case. The main goal of Wang’s work was to estimate the friction coefficient, and the temperature increase along the line of contact, so that a thermal EHL solver was needed; in the present work the focus is on the dynamic in the normal direction with respect to the contact (squeezing film), thus the thermal effects are neglected at the first stage. Nonetheless, Wang’s results show that neglecting thermal effects leads to overestimating the film thickness in the access portion of the contact segment.

M. Barbieri et al. / Tribology International 62 (2013) 37–48

In order to correlate the iso-thermal line solution of Ref. [6] with the present point solution, the approach proposed by Nijenbanning [18] is considered: a point load w, equivalent to the line load w‘ , is

Fig. 6. Detail of the method used to define the initial guess pressure on the contact B. Table 2 Gear data from Ref. [6]. Data Number of teeth Normal module (mm) Pressure angle (deg) Outer radius (mm) Crowning magnitude (mm) Face width (mm) Young’s modulus (MPa) Poisson’s coefficient Center distance (mm)

Z m

a Ra mag c Fw E

n a0g

t

Transmission ratio Transmitted power (W) Rotational speed (RPM) Torque at pinion (Nmm)

P n1 Tg1

Pinion

Gear

35 2.00 20.00 36.50 0.020 20.00 226,000 0.30 174.50

140 2.00 20.00 141.04 0.020 20.00 226,000 0.30

43

defined by imposing that the width of the contact pattern of the Hertz line solution is equal to the minor axis of the contact ellipse in the point solution. Such equivalent load is expressed as follows: sffiffiffiffiffiffiffiffi 2 32Dn wn‘ 2Rnx wpnh n w ¼ ð27Þ 3p pE0 kn Data of the case study are collected in Table 2. In Ref. [6], the outer radii Ra of both pinion and gear were not defined: here they are chosen so that the change of the number of mating tooth pairs occurs as in Fig. 4(c) of Ref. [6]. The rotational speed is inferred from pressure graphs in [6]: the pitch point pressure matches if the rotational speed is set to 1200 RPM. The comparison in Fig. 7(a) shows a very good agreement in terms of central pressure pc ¼ pð0,0,tÞ. The solution of Ref. [6] is less sharp close to the load jumps; far from the pitch point (x ¼ 0) the agreement is no more excellent, this is because the equivalent load is defined at the pitch point. Note that the total contact spans between 7 4:154 mm and the contact length is given by the position of load jumps. In Ref. [6], the contact is between 7 5:2 mm, which is non realistic. Since the contact is highly loaded, the central pressure is very close to the hertzian value. Fig. 7(b) shows that central and minimum film thickness are both higher and sharper than Ref. [6]. The number of time steps is 180 in [6], and 1537 in the present elliptic solution: here this number is chosen in order to the have the same spatial and time resolution in terms of the dimensionless variables. It is well known that, for highly loaded contacts, a poor accuracy in the numerical solution may result in underestimating the film thickness [20]. In order to check this hypothesis, a further comparison has been drawn with stationary EHL interpolating formulae presented in [21]. This comparison proves that the present solution is almost stationary far from the load jumps: at the pitch point, the central film thickness is 0:592 mm, while the value obtained using the formula of Ref. [21] is 0:607 mm, the difference is 2.5%.

4. Results for a case study The main goal of the present work is to point out the role played by the gear dynamics on gear lubrication conditions; in order to achieve such result, a highly loaded high speed gear pair is considered (see Table 3). Note that the load is higher than in common gear pairs, so that the contact pressure is raised up to 3 GPa. The first step is the characterization of the dynamical system, to this extent it is useful to consider a simplified problem: the effect of the lubricant is neglected, so that Eq. (4) can be solved considering the Hertzian pressure. Eq. (4) is now linearized in the neighborhood of the static equilibrium position at the pitch

0.25 12,000 1200 95,490

1000

0.8

900

0.75 0.7

hc, hmin [ m]

pc [MPa]

800 700 600 500

0.65 0.6 0.55 0.5 0.45

400 300 -6

0.4 -4 -2 0 2 4 Position on the contact line - [mm]

Fig. 7. Comparisons with results in Ref. [6]: (a) present approach.

pc from Ref. [6],

6

0.35 -6

pc present approach; (b)

-4 -2 0 2 4 Position on the contact line - [mm] hc from Ref. [6],

hmin from Ref. [6],

6

hc present approach,

hmin

44

M. Barbieri et al. / Tribology International 62 (2013) 37–48

point; the resulting equivalent stiffness is the following: 3w 2 hn0

n

kel ¼

ð28Þ

Table 3 Data set of the case study. Data Number of teeth Normal module (mm) Pressure angle (deg) Outer radius (mm) Crowning magnitude (mm) Young’s modulus (MPa) Poisson’s coefficient Center distance (mm)

Z m

a Ra mag c E

n

Gear

28 3.00 20.00 46.55 0.131 206,000 0.30

43 3.00 20.00 69.85 0.131 206,000 0.30

a0g

111.00

n1 Tg1

15,000 470,000

3500

16

3000

14

2500

12 hc [ m]

pc [MPa]

Rotational speed (RPM) Torque at pinion (Nmm)

Pinion

Using the equivalent stiffness kel and the equivalent mass mge as defined in (2), the corresponding ‘‘natural frequency’’ fn can be computed. The gear pair under investigation presents a natural frequency f n ¼ 5142 Hz; the rotational speed of 15,000 RPM which corresponds to a mesh frequency f m ¼ 7000 Hz. Note that the mesh frequency is the excitation frequency of the system. The test case described in Table 3 is characterized by the following dimensionless parameters: M n ¼ 2401:12, Ln ¼ 27:40, In ¼ 31:13, Dn ¼ 0:06. This combination of parameters shows that, if steady-state conditions are assumed, the lubrication regime at the pitch point is elastic piezo-viscous, i.e. elastic deformation in the contacting gears must be taken into account, as well as the pressure–density dependence. The test case is analyzed by means of the full ‘‘dynamic model’’ described by Eq. (4). Eqs. (22) and (25) are solved starting from stationary EHL conditions at the pitch point for contact 2. Note that, at the pitch point contact 1 is inactive and gives no load contribution. The initial central pressure and central film thickness can be found in Fig. 8. Moving

2000 1500

8 6

1000

4

500 0

10

2 0

20 40 60 80 Position on the contact line -

100 [mm]

0

120

0

Fig. 8. Transient response for the first 12 Tm:

20 40 60 80 Position on the contact line contact 2,

100 [mm]

120

contact 1.

0.8 0.5 0.6 0.4

-0.5

0.2

V0(T)

H0,V0

0

0

-1

-0.2 -1.5

-0.4

-2 1175

1180 1185 1190 1195 Dimensionless time - T

1200

-0.6 -2

-1.5

-1

-0.5

H0(T) X: 7007 Y: 0.3154

0.3

|FFT(H0(f))|

0.25 0.2 0.15 0.1 X: 5217 X: 1790 Y: 0.02678 Y: 0.0106

0.05

X: 1.401e+04 Y: 0.03205

00

00 60 00 80 00 10 00 0 12 00 0 14 00 0 16 00 0

40

20

0

0

Frequency (Hz) Fig. 9. Steady-state solution (a) time history,

H0 ,

V0; (b) steady-state phase portrait; (c) spectrum of the dimensionless rigid approach.

M. Barbieri et al. / Tribology International 62 (2013) 37–48

the gear pair along the contact line, at x ¼ 3:5 mm contact 1 becomes active: its pressure is set assuming steady-state conditions on the entering contact. At x ¼ 5:76 mm, contact 2 expires, and contact 1 bears all the load. During the one tooth pair contact, the central pressure is minimum: this is because the system is forced over the resonance frequency f m 4 f n . At 12.10 mm, contact 2 becomes active again. This scheme is iterated up to convergence to a steady state solution in terms of approach H0 and velocity V0. The time integration is performed using a second-order implicit scheme, which is more stable for stiff problems. Nonetheless, since the multilevel method is iterative, implicit time integration requires to solve (iteratively) the problem for each time step. For the initial time step and for all those steps in which a contact becomes active or non-active, the integration scheme is switched to first-order implicit, and much more iterations are needed to get the converged solution. From Fig. 8(a), it possible to see that for x ¼ 16:27 the pressure suddenly drops to zero in contact 2, and the corresponding central film thickness is raised up to eight times the stationary value, see Fig. 8(b). This behavior corresponds to a loss of contact and is due to the dynamic effect in the gear pair: the initial conditions were arbitrarily assumed as stationary EHL on a single contact. It is to stress that the problem is continuously transient, due to the entering/exit of contacts, nonetheless it gives a periodic steady-state response in terms of the rigid approach H0. Fig. 9(a) shows the steady-state response of the system after 156 mesh periods. It is to note that the velocity V0 is not C0

1.4

1 0.8 0.6 0.4 0.2 0

0

5

25

30

Position on the contact line -

10

[mm]

15

20

Fig. 10. Dimensionless contact force:

contact 2,

35

40

contact 1.

everywhere: when the number of contacts changes, there is a jump in acceleration, due to the sudden stiffness change. This behavior is shown in Fig. 9(b) as well: the orbit in the phase portrait is non smooth, and the singular points are the maxima and the minima of V0. The spectrum of the rigid approach (Fig. 9(c)) presents a major peak at the forcing frequency, i.e. the mesh frequency f m ¼ 7000 Hz, and a lower peak at 2f m . A further peak at 5217 Hz is present: it is due to the transient response at the elastic natural frequency f n ¼ 5142 Hz, which can also be computed by means of Eq. (28). There is a third peak at low frequency 1790 Hz ¼7007  5217 Hz. The combination of the frequencies is expected because of the quadratic type nonlinearities. The contact force in each contact can be computed by numerical integration of the pressure. Its value, normalized with respect to wn , is shown in Fig. 10. At the present rotational velocity (n1 ¼ 15,000 RPM), the maximum dynamic load amplification is 1.29. The steady state response proves that the load is almost equally shared by the contacts, when they are both active. A difference is present close to the maxima, where entering and exit of contacts occur: this behavior is an effect of coupling between the two EHL contacts. Central and maximum pressure in the EHL contacts are shown in Fig. 11. The maximum pressure is 3.28 GPa, which exceeds commonly accepted values for gear contacts, as discussed before. As expected in highly loaded contacts, the maximum pressure is always located in the center of the contact for the most of the contact time. Nonetheless, close to the beginning of the contact (Fig. 11(b), x ¼ 44:5 mm), the maximum is moving through the contact area, and the two curves do not overlap. The steady-state film thickness is represented in Fig. 12(a). Let us follow contact 1: the minimum film thickness is found close to the starting contact point (point A), minðhmin Þ ¼ 1:9 mm; the corresponding central film thickness is 2:2 mm. As expected, the increment of the equivalent curvature radius Rx ðxÞ is responsible for an increase of both the central and the minimum film thickness along the contact, so that the most critical point is still the beginning of the contact; nonetheless, the oscillation of the film thickness due to the dynamic effect is prevalent. After point A of Fig. 12(a), the fluid film is developed and in point C a local maximum of the central film thickness occurs. The effect of the exit of contact 2 at x ¼ 5:76 mm is delayed and the following minimum of the central film thickness of contact 1 (point D, x ¼ 7:38 mm) occurs 1.3 mm after the contact exit. Note that the corresponding minimum in the central pressure (Fig. 11(a)) is further delayed of 1.7 mm with respect to point D. This delay is explained by transient effects in EHL film: in highly loaded contacts the viscosity in the center of contact is high, and thus the Poiseuille term in the Reynolds equation is negligible. For this

3500

3300

3000

3200

2500

pc,pmax [MPa]

pc,pmax [MPa]

Normalized contact force

1.2

2000 1500 1000

45

3100 3000 2900 2800

500 2700 0

0

2

4 6 8 10 12 14 16 Position on the contact line - [mm]

Fig. 11. Central and minimum pressure:

18

pc contact 2,

3

3.5 4 4.5 5 5.5 6 Position on the contact line - [mm]

pmax contact 2,

pc contact 1,

pmax contact 1.

46

M. Barbieri et al. / Tribology International 62 (2013) 37–48

E

3

3

D A

B

hchm [ m]

hchmin [ m]

F

C

2.5 2

3.2

1.5

2.8 2.6 2.4

1 2.2 0.5 2

A

B

C

D

E

F

hc contact 2,

hmin contact 2,

reason, the Reynolds equation can be reduced to a transport equation of rh. The transport of rh can be observed in Fig. 12(b), where a section of the film thickness along the plane X,h is shown. At the begin of contact 1 (point A), the film thickness is flat in the center of the contact and has a ripple at the end, as usual in stationary EHL. In point B, the central film thickness is unchanged, a ripple of h is moving from right to left and has a minimum in the center of the contact; this explains why the increment of central film thickness is delayed. In point C the ripple has moved away and a maximum in central film thickness occurs. At point D the minimum of h is due to the presence of another moving ripple, which was started by the exit of contact 2. At point E the conditions are again almost stationary and the central part of the film thickness is flat; the global maximum of central film thickness occurs. Note that the point E is close to the maximum value of the central pressure, which is located 1.78 mm after point E. In point F, where contact 1 exits, the central film thickness is reduced, because of the ripple due to entering of contact 2, which happened 2.5 mm before. In point F, the minimum of film thickness occurs on the left side of the contact, because the moving ripple is more relevant than the stationary ripple on the right side. Fig. 13 shows the contour plots of h, which correspond to the contact positions A–F of Fig. 12. The contours are linearly spaced. Points A and E present almost stationary conditions; in point B the motion of the ripple can be observed; in points C and F two different ripples are present, one on the left and one on the right side. At point D, the contact area is notably smaller and the global minimum of h occurs in the center of the contact. A comparison with stationary EHL value can be drawn: if the pitch point is considered (i.e. x ¼ 0 for contact 2 in Fig. 12(a)), the dynamic central and minimum film thickness are 2:495 mm and 2:127 mm, while the corresponding stationary EHL values are 2:415 mm and 2:163 mm. Therefore, the stationary EHL model is well representing the conditions at the pitch point, it can be used to estimate lubrication conditions in gear pairs. Point D presents the most dangerous non stationary conditions: the effective contact area is reduced, as well as the central contact pressure, and both the central and the minimum film thickness are smaller than the corresponding stationary EHL values. This effect cannot be modeled without considering coupling of EHL lubrication with the dynamics of the gear pair. A more refined comparison can be drawn with the transient static model: i.e. the transient EHL solver coupled with Eq. (6). Fig. 14 shows that in the dynamic model the shape of the contact pressure is given mainly by the dynamics of the mating teeth, and the variables are smoother thanks to the inertia. It is clear that the pressure spikes in the transient static model are non-realistic. The minimum pressure in the dynamic model is found close to the pitch point (x ¼ 8 mm), when only one pair of teeth is meshing: this is because the system behaves like an oscillator excited over

1

X

Position on the contact line - [mm] Fig. 12. Film thickness, steady-state conditions: (a)

0 0. 2 0. 4 0. 6 0. 8

-1 -0 .8 -0 .6 -0 .4 -0 .2

1

16

14

12

8

10

6

4

2

0

0

hc contact 1,

hmin contact 1; (b) film thickness for reference positions.

Fig. 13. Contour plots of the film thickness at different contact positions (see Fig. 12(a)). (a) Point A. (b) Point B. (c) Point C. (d) Point D. (e) Point E. (f) Point F.

the natural frequency. In terms of film thickness, the dynamic effects are responsible for a thinner film in the access part of the gear meshing with respect to the one computed using the static transient model. Conversely, the film in the recess part is thicker than in the transient static simulation.

5. Conclusions In this paper a fully dynamic, transient EHL model applied to gear dynamics is presented. The model is capable to describe the

47

3.2

3600 3400 3200 3000 2800 2600 2400 2200 2000 1800 1600

3 2.8

hc, hmin [ m]

pc, pmax [MPa]

M. Barbieri et al. / Tribology International 62 (2013) 37–48

2.6 2.4 2.2 2 1.8 1.6

2

4 6 8 10 12 14 Position on the contact line - [mm]

1

2

4 6 8 10 12 14 Position on the contact line - [mm]

16

Fig. 14. Comparison between static transient model and dynamic model: (a) dimensional pressure, hc static transient model, pmax static transient model, pc dynamic model, pmax dynamic model; (b) central and minimum film thickness, hc static transient model, hmin static transient model, hc dynamic model, hmin dynamic model.

effect of coupling between transient point EHL lubrication in spur gear pairs and the dynamic behavior of the gear pair itself. The model is applied to gears having a contact ratio lower than 2, so that two EHL contacts share the total dynamic load. The generalization to higher contact ratios is formally straightforward. The bending stiffness of the teeth and of the wheel has been neglected in order to focus the attention on the lubrication conditions. The model is validated by comparison with the literature. After transforming the governing equations into a dimensionless form, it is shown that three dimensionless parameters (Mn, Ln, and In) characterize the dynamic EHL conditions of the system. The dynamic model has been applied to a highly loaded spur gear pair; the problem is locally unsteady, due to the sudden entering/ exit of new contacts, nonetheless the system reaches a periodic steady-state response. The steady-state contact force is not the same in the two contacts: this is an indication that the two contacts interact with each other. Even though the test case is far from the elastic resonance, gear dynamics largely affects the fluid film lubrication between the teeth: the maxima and the minima of film thickness along the contact can be correlated to the minima and the maxima in the dynamic contact force, provided that delays due to the squeeze term inside the Reynolds equation are considered. The worst lubrication conditions can be found close to the pitch point in the access part of the contact, when the effect of the exit of the other contact is propagated into the fluid film, leading to a low film thickness spreading over a small contact area. A comparison with standard static transient EHL models (i.e. assuming a given static external load on a single lubricated contact) shows that neglecting the dynamic effects leads to the appearance of unrealistic pressure spikes along the contact segment, and to underestimation of the severe lubrication conditions which occur in the access portion of the contact segment. Moreover, if the dynamic model is considered, the minimum contact pressure occurs close to the pitch point, whereas static models assume that the contact pressure is maximum near the pitch point.

þ

x i,j1=2,k P i,j1,k ðx i,j1=2,k þ x i,j þ 1=2,k ÞP i,j,k þ x i,j þ 1=2,k P i,j þ 1,k DY 2

U m 

1:5r i,j,k Hi,j,k 2r i1,j,k Hi,j,k þ 0:5r i2,j,k Hi2,j,k

DX 1:5r i,j,k Hi,j,k 2r i,j,k Hi,j,k1 þ 0:5r i,j,k Hi,j,k2 DT

¼0

ðA:1Þ

where

x i 7 1=2,j,k ¼

x i,j,k þ x i 7 1,j,k

x i,j 7 1=2,k ¼

2

x i,j,k þ x i,j 7 1,k 2

ðA:2Þ

The steps are chosen so that: DX ¼ DY ¼ DT. The equation of motion (25) is discretized as follows: 8 PP 1:5V 0k 2V 0k1 þ 0:5V 0k2 DX DY PP 2p > > ¼ P 1,i,j,k þ DXInDY P2,i,j,k  n n > < DT In 3k I i j i j > 1:5H0k 2H0k1 þ 0:5H0k2 > > ¼ V 0k : DT

ðA:3Þ Whenever a new contact becomes active, the initial pressure at k1 is set as if the contact was a stationary EHL, and a first-order discretization scheme is used for the first step. The corresponding discretized squeeze term reads   r i,j,k Hi,j,k r i,j,k Hi,j,k1 @ ðr HÞ ðA:4Þ ¼ @T DT i,j,k while the equation of motion becomes 8 V 0k V 0k1 DX DY PP DX DY PP 2p > > ¼ P 1,i,j,k þ n P 2,i,j,k  n n > < DT In I 3 k I i j i j > H0k H0k1 > > ¼ V 0k : DT

ðA:5Þ

For the first steps after the new contact entering, the time step is reduced in order to avoid numerical instability. References

Appendix A In this section, a brief overview of the numerical scheme is given. If a new contact has not entered at the k1 step, Eq. (22) is discretized for the k step using a second-order finite difference scheme as follows:

x i1=2,j,k P i1,j,k ðx i1=2,j,k þ x i þ 1=2,j,k ÞPi,j,k þ x i þ 1=2,j,k Pi þ 1,j,k DX 2

[1] Harris SL. Dynamic loads on the teeth of spur gears. Proceedings of the Institution of Mechanical Engineers 1958;172(2):87–100. [2] Kahraman A, Singh R. Non-linear dynamics of a spur gear pair. Journal of Sound and Vibrations 1990;142(1):49–75. [3] Kahraman A, Blankenship GW. Experiments on nonlinear dynamic behavior of an oscillator with clearance and periodically time-varying parameters. Journal of Applied Mechanics 1997;64:217–26. [4] Grubin AN. Fundamentals of the hydrodynamic theory of lubrication of heavily loaded cylindrical surfaces; 1949. Moscow: TsNIITMASh. Book 30 (Engl. transl.).

48

M. Barbieri et al. / Tribology International 62 (2013) 37–48

[5] Larsson R. Transient non-newtonian elastohydrodynamic lubrication analysis of an involute spur gear. Wear 1997;207:67–73. [6] Wang T, Li H, Tong J, Yang P. Transient thermoelastohydrodynamic lubrication analysis of an involute spur gear. Tribology International; 2004 p. 37. [7] Evans HP, Snidle RW, Sharif KJ. Deterministic mixed lubrication modelling using roughness measurements in gear applications. Tribology International 2009;42(10):1406–17. [8] Osman T, Velex P. A model for the simulation of the interactions between dynamic tooth loads and contact fatigue in spur gears. Tribology International 2012;46(1):84–96. [9] Li S, Kahraman A. A transient mixed elastohydrodynamic lubrication model for spur gear pairs. Journal of Tribology 2010:132. [10] De la Cruz M, Chong WWF, Teodorescu M, Theodossiades S, Rahnejat H. Transient mixed thermo-elastohydrodynamic lubrication in multi-speed transmissions. Tribology International 2012;49:17–29. [11] Barbieri M, Pellicano F. Coupling of two EHL-lubricated contacts in gear dynamics In: STLE/ASME 2010 international joint tribology conference proceedings; 2010. [12] Li S, Kahraman A. Influence of dynamic behaviour on elastohydrodynamic lubrication of spur gears. Proceedings of the IMechE Part J: Engineering Tribology 2011:225.

[13] Venner CH, Lubrecht AA. Multilevel methods in lubrication. In: Dowson D, editor. Tribology series, vol. 37. Elsevier; 2000. [14] Kumar P, Saini PK, Tandon P. Transient elastohydrodynamic lubrication analysis of an involute spur gear using couple-stress fluid. Proceedings of the IMechE Part J: Engineering Tribology 2007:221. [15] Dowson D, Higginson GR. Elastohydrodynamic lubrication, the fundamentals of roller and gear lubrication. Oxford: Pergamon Press; 1966. [16] Roelands CJA, Vlugter JC, Waterman HI. The viscosity–temperature–pressure relationship of lubricating oils and its correlation with chemical constitution. Journal of Basic Engineering 1963;85(4):601–7. [17] Moes H, Bosma R. Design charts for optimum bearing configuration, I the Full Journal Bearing. ASME Journal of Tribology 1971;93:302–6. [18] Nijenbanning G, Venner CH, Moes H. Film thickness in elastohydrodynamically lubricated elliptic contacts. WEAR 1994;176:217–29. [19] Barbieri M. Gear dynamics under EHL conditions, Universita di Modena e Reggio Emilia, PhD dissertation; 2010. [20] Venner CH. EHL Film thickness computations at low speeds: risk of artificial trends as a result of poor accuracy and implications for mixed lubrication modelling. Journal of Engineering Tribology 2005;219:285–90. [21] Canzi A, Venner CH, Lubrecht AA. Film thickness prediction in elastohydrodynamically lubricated elliptical contacts. Proceedings of the IMechE Part J: Journal of Engineering Tribology 2010;224(9):917–23.