Study on bifurcation characteristics of multi-clearance bending torsional coupling gear transmission based on EHL

System model description

The 6-DOF bending torsional coupled vibration model of system considering friction, backlash, and support clearance is shown in Figure 1.

OPEN IN VIEWER

Figure 1.

6-DOF bending–torsional coupling vibration model.

In Figure 1, the angular displacement of driving gear and driven gear is represented by θ_s and θ_p, respectively. The lateral displacement is represented by X_s and X_p, respectively. The longitudinal displacement is represented by Y_s and Y_p, respectively. The base circle radius is represented by r_bs and r_bp, respectively. The stiffness, half support clearance, and damping in the x and y directions of the driving gear are represented by K_sx, K_sy, b_sx, b_sy, C_sx, and C_sy, respectively. The stiffness, half support clearance, and damping in the x and y directions of the driven gear are represented by K_px, K_py, b_px, b_py, C_px, and C_py, respectively. Here, i is x and y. The comprehensive meshing error, time-varying mesh stiffness, half backlash, meshing damping, and friction coefficient of the gear pair are expressed by e(t), K_sp(t), b_h, C_sp, and μ, respectively. The friction arm of driving gear and driven gear is represented by L_s and L_p, respectively. The eccentricity errors of driving gear and driven gear are expressed by e_s and e_p, respectively. The initial phase angles of eccentricity errors are α_s and α_p, respectively. The phase angle between the driven gear and the driving gear is δ. The linear displacements x_s and x_p produced by the angular displacements of the driving gear and the driven gear on the meshing line can be shown in equation (1)

{ x s = θ s r bs x p = θ p r bp

(1)

The comprehensive meshing error includes the static transmission error and eccentricity error of driving gear and driven gear. The harmonic period expansion of meshing transmission error by Fourier function can be expressed as follows

e ( t ) = E + ∑ b = 1 l e s sin ( b ω s t + α s ) + ∑ b = 1 l e p sin ( b ω p t + α p )

(2)

In equation (2), E is the average value of static transmission error. α_s and α_p are the initial phase angles of the bth harmonic component, respectively. l is the number of harmonic terms expanded by Fourier. e_s and e_p are the eccentricity error. ω_s and ω_p are rotation frequency of driving gear and driven gear (rad/s).

The non-linear functions of backlash and support backlash can be expressed in equation (3)

f ( X , b h ) = { X − b h X > b h 0 | X | ≤ b h X + b h X < − b h , f sx ( X s , b sx ) = { X s − b sx X s > b sx 0 | X s | ≤ b sx X s + b sx X s < − b sx , f sy ( Y s , b sy ) = { Y s − b sy Y s > b sy 0 | Y s | ≤ b sy Y s + b sy Y s < − b sy f px ( X p , b px ) = { X p − b px X p > b px 0 | X p | ≤ b px X p + b px X p < − b px , f py ( Y p , b py ) = { Y p − b py Y p > b py 0 | Y p | ≤ b py Y p + b py Y p < − b py

(3)

The dynamic meshing force of gear pair is defined as the sum of elastic restoring force and damping force, which can be expressed as follows

F N = K sp ( t ) f ( X , b h ) + C sp X ·

(4)

In equation (4), X is the relative displacement of the gear pair’s meshing point.

The damping coefficient is expressed as follows

C sp = 2 ξ K m / ( 1 / m s + 1 / m p )

(5)

In equation (5), ξ is the relative damping ratio of the meshing pair, m_s is the mass of the driving gear, m_p is the mass of the driven gear, and the relative damping ratio of the gear transmission system is generally 0.03 to 0.17.

Establishment of dynamic backlash b_h(t) based on fractal theory

The W-M model of gear pair meshing based on the fractal theory is shown in Figure 2. A is the highest point of surface profile of gear 1. B is the highest point of surface profile of gear 2. The research in Yin et al. ²⁷ shows that the machined surface has fractal characteristics, and the tooth surface is also machined surface. The fractal formula of backlash can be defined by W-M function as follows

b h ( t ) = ∑ k = 0 + ∞ λ b ( D − 2 ) k sin ( λ b k t ) ; λ b > 1 , 1 < D < 2

(6)

In equation (6), b_h(t) is the backlash varying with the time. λ_b is the characteristic scale coefficient. D is the fractal dimension, and only two-dimensional fractal is considered. k is the contact coefficient. t is time.

OPEN IN VIEWER

Figure 2.

W-M model of gear pair meshing.

Using equation (6), the simulation of backlash under different fractal dimensions is shown in Figure 3. Here, λ is 1.5 and D is 1.1–1.9, respectively. The larger the fractal dimension D, the more complex the micro morphology of W-M simulation. The curve changes more dramatically with the increase in fractal dimension D. The larger the fractal dimension D, the smoother the tooth surface and the smaller the backlash. The average backlash b_hm with different fractal dimension D is shown in Table 1.

OPEN IN VIEWER

Figure 3.

The simulation of backlash (b_h) with different fractal dimensions (D): (a) D = 1.1, (b) D = 1.2, (c) D = 1.3, (d) D = 1.4, (e) D = 1.5, (f) D = 1.6, (g) D = 1.7, (h) D = 1.8, and (i) D = 1.9.

OPEN IN VIEWER

Table 1.

The average backlash b_hm in different fractal dimension D.

D	1.1	1.2	1.3	1.4	1.5	1.6	1.7	1.8	1.9
b hm	1.5472	0.9744	0.6140	0.3872	0.2444	0.1544	0.0978	0.0620	0.0394

Calculation of time-varying tooth surface friction based on EHL

The schematic diagram of gear force analysis is shown in Figure 4.

OPEN IN VIEWER

Figure 4.

Friction analysis of gear pair.

In the process of gear transmission, the tooth surface will be affected by the dynamic meshing force F_N along the direction of the meshing line and the friction force F_f perpendicular to the meshing line. XY is the theoretical meshing line length. AD is the actual meshing line length. r_as and r_ap are the top circle radius of driving gear and driven gear. r_s and r_p are the pitch circle radius of driving gear and driven gear. ω_s and ω_p are the rotation frequency of driving gear and driven gear, rad/s. α_n is the meshing angle.

From the geometric relationship of gear transmission, we can get the length XY of theoretical meshing line, the length AD of actual meshing line, the distance XA between the start end of actual and theoretical meshing line, and the distance YD between the end of actual and theoretical meshing line. The calculation expression is as follows

{ XY = ( r s + r p ) sin α n AD = r as 2 − r bs 2 + r ap 2 − r bp 2 − ( r s + r p ) sin α n XA = ( r s + r p ) sin α n − r ap 2 − r bp 2 YD = ( r s + r p ) sin α n − r as 2 − r bs 2 P b = 2 π r bs Z s = 2 π r bp Z p

(7)

In equation (7), P_b is the pitch of gear. Z_s and Z_p are the number of teeth of driving gear and driven gear, respectively.

Therefore, the expression of the friction arm can be derived as follows

{ L s ( t ) = XA + P b + mod ( ω s r bs t , P b ) L p ( t ) = YD + mod ( ω p r bp t , P b )

(8)

In equation (8), L_s and L_p are the time-varying friction arms of the driving gear and driven gear, respectively, mod (·) is the complementary function.

The relative velocity u_e of two meshing contact tooth surfaces at any position is shown in equation (9)

u e ( t ) = L s ( t ) ω s − L p ( t ) ω p

(9)

Since the direction of friction force will change at the node, the direction coefficient λ can be used to represent the direction change. The direction coefficient is defined as

λ ( t ) = sgn ( u e ( t ) ) = { 1 u e ( t ) < 0 0 u e ( t ) = 0 − 1 u e ( t ) > 0

(10)

In equation (10), sgn(·) is a symbolic function.

The calculation formula of Coulomb’s law friction model is as follows

μ ( t ) = μ 0 sign ( AP − mod ( r bs ω s t , P b ) ) = { μ 0 AP < mod ( r bs ω s t , P b ) 0 AP = mod ( r bs ω s t , P b ) − μ 0 AP > mod ( r bs ω s t , P b )

(11)

In equation (11), μ₀ is the amplitude of friction coefficient.

The time-varying sliding friction force of tooth surface according to Coulomb’s law can be expressed as follows

F f ( t ) = λ ( t ) μ ( t ) F N

(12)

In equation (12), μ is the friction coefficient.

In aviation standard, the empirical formula of friction coefficient at any meshing point K synthesized according to the actual test results is as follows

μ = 0 . 12 ( W t R a 2 ν 0 u sp R ) 0 . 25

(13)

In equation (13), ν₀ is the initial viscosity of the lubricating oil (Pa/s). W_t is the tangential load per unit tooth width (N/mm). u_sp is the entrainment velocity (m/s). R is the comprehensive curvature radius of the gear pair (mm). R_a is the mean value of tooth surface roughness (μm). Ignoring the effect of roughness, R_a is equal to zero.

The friction coefficient is related to the geometric shape of the teeth, the surface roughness of the teeth, the relative sliding velocity of the teeth, the contact stress between the teeth, and the lubrication of the system, so it is difficult to accurately predict the friction coefficient. When the lubrication effect of tooth surface is considered, the gear contact problem becomes the linear contact EHL problem. Therefore, the friction coefficient is calculated by EHL theory.

By solving Reynolds equation, viscous pressure equation, dense pressure equation, energy balance equation, and load equation, the oil film pressure distribution and film thickness distribution are obtained. Combined with Ree–Eyring fluid equation, the oil film shear stress distribution is calculated, and then, the friction force and friction coefficient between two contact surfaces are obtained. It is assumed that the tooth surface is completely smooth, that is, the gear is in the state of complete EHL. In the lubrication analysis, the influence of temperature effect is ignored, and it is considered that the oil film thickness has not changed.

The Reynolds equation is

{ ∂ ∂ x [ ( ρ ν ) e h 3 ∂ p ∂ x ] = 12 u sp ∂ ( ρ * h ) dx u sp = ( u s + u p ) / 2 p x in = p x out = p x = 0 ∂ p ∂ x | x out = 0

(14)

In equation (14), x is the tangent direction of tooth surface. ρ is the density of lubricating oil (kg/m³). h is the oil film thickness (m). p is oil film pressure (Pa). ν is the viscosity of Ree-Eyring fluid (Pa·s). u_sp is the entrainment velocity (m/s), u_s, u_p are the tangential velocity of lubricating oil along gear axis (m/s). x_in and x_out are the entrance and exit positions of the calculation area in the x-direction. ρ* and (ρ/ν) _e are the equivalent parameters. The equivalent parameters are defined in equation (15)

{ ( ρ ν ) e = 12 ( ν e ρ ′ e ν ′ e − ρ ″ e ) ρ * = [ ρ ′ e ν e ( u s − u p ) + ρ e u s ] / u sp ρ e = 1 h ∫ 0 h ρ d z ρ ′ e = 1 h 2 ∫ 0 h ρ ∫ 0 z dz ′ ν * d z ρ ″ e = 1 h 3 ∫ 0 h ρ ∫ 0 z z ′ d z ′ ν * d z 1 ν e = 1 h ∫ 0 h d z ν * 1 ν ′ e = 1 h 2 ∫ 0 h z d z ν *

(15)

The boundary conditions for the differential equations are p(x_in) = 0, p(x_out) = 0, p(x_in) > 0, and (x_in < x < x_out). ν* is the equivalent viscosity of lubricating oil (Pa·s).

The film thickness equation is as follows

{ h ( x ) = h 0 + x 2 2 R + δ ( x ) δ ( x ) = − 2 π E ′ ∫ x in x out p ( Sr ) ln ( x − Sr ) 2 d Sr + C d

(16)

In equation (16), h₀ is the central film thickness (m) of the two gears without meshing deformation, which can be determined by the load balance condition. R is the comprehensive curvature radius (m). δ(x) is the elastic deformation displacement (m) due to the contact pressure of the tooth surface. E′ is the comprehensive elastic modulus (Pa) of the two gears. Sr is the sliding roll ratio when the gear is engaged. C_d is the viscosity pressure coefficient of Cameron’s viscosity pressure formula (m²/N), C_d = χ/15, and χ is the viscosity pressure coefficient of Barus’ viscosity pressure formula (m²/N).

The viscous pressure equation and the dense pressure equation can be expressed as follows

{ ν = ν 0 exp { ( ln ν 0 + 9 . 67 ) [ − 1 + ( 1 + 5 . 1 × 10 − 9 p ) z 0 ] } ρ = ρ 0 [ 1 + 0 . 6 × 10 − 9 p 1 + 1 . 7 × 10 − 9 p ] z 0 = χ / [ 5 . 1 × 10 − 9 × ( ln ν 0 + 9 . 67 ) ]

(17)

In equation (17), ν₀ is the initial viscosity of the lubricating oil (Pa·s). ρ₀ is the initial density of the lubricant (kg/m³).

The load balance equation can be expressed as follows

{ c w ( t ) W = ∫ x in x out p ( x ) dx c w ( t ) = { 1 3 ( 1 + Sr − x in x 0 ) x in ≤ Sr ≤ x in + x 0 1 x in + x 0 ≤ Sr ≤ x out − x 0 1 3 ( 1 + Sr − x out x 0 ) x out − x 0 ≤ Sr ≤ x out x 0 = π m n cos α n ( ε − 1 )

(18)

In equation (18), W is the equivalent cylinder external load (N). c_w(t) is the load factor of gear. Taking the meshing line of gear transmission as the coordinate axis and the node as the origin, x_in and x_out are the coordinates of entering and leaving the meshing point, respectively, m_n is the modulus and ε is the gear coincidence degree.

The constitutive equation of Ree-Eyring fluid can be expressed as follows

{ τ = ν * ∂ u ∂ z ν * = ν ( τ / τ 0 ) / sin h ( τ / τ 0 )

(19)

In equation (19), τ is the shear stress of oil film in x-direction (Pa). u is the velocity of lubricating oil in x-direction (m/s). τ₀ is the shear stress of Ree–Eyring fluid (Pa).

The formula of velocity and velocity gradient derived from the motion equation of fluid particle can be expressed in equation (20)

{ u = ( ∫ 0 z z ′ v * d z ′ − v e v ′ e h ∫ 0 z 1 v * d z ′ ) ∂ p ∂ x + v e h ( u p − u s ) ∫ 0 z 1 v * d z ′ + u s ∂ u ∂ z = ( z v * − v e v ′ e h v * ) ∂ p ∂ x + v e h ( u p − u s ) 1 v *

(20)

The friction force on the contact tooth surface of the gear and the velocity distribution of the lubricant are shown in Figure 5.

OPEN IN VIEWER

Figure 5.

Schematic diagram of velocity distribution and friction force.

Here, the sliding friction is mainly generated in the contact area (−b, b). F₁ and F₂ are the friction force on the surface of the two tooth surfaces, respectively. The sliding friction coefficient of the gear based on the EHL model can be expressed as follows

μ = F 1 + F 2 2 W = ∫ − b b τ | z = 0 d x + ∫ − b b τ | z = h d x 2 W

(21)

In equation (21), F 1 = ∫ − b b τ | z = 0 d x and F 2 = ∫ − b b τ | z = h d x .

The above mathematical model is dimensionless. Then, the multi-grid method is used to solve the pressure, and the numerical method is used to solve the friction coefficient and friction.

In order to verify the correctness of the proposed sliding friction coefficient, this article gives the friction coefficient calculated by EHL empirical formula in Xu et al. ³⁴ and makes a comparison. Calculation model expression of friction coefficient based on EHL empirical formula proposed by Xu et al. ³⁴ is as follows

μ = e f ( sr , P h , ν m , R rms ) P h b 2 | Sr | b 3 u sp b 6 ν m b 7 R b 8

(22)

In equation (22), P_h is the maximum contact Hertz stress (GPa). R_rms is the root mean square (RMS) value of roughness.

f (Sr, P_h, ν_m, R_rms) is calculated in equation (23) as follows

f ( Sr , P h , ν m , R rms ) = b 1 + b 4 | Sr | P h log 10 ( ν m ) + b 5 e − | sr | P h log 10 ( ν m ) + b 9 e R rms

(23)

The values of b₁, b₂, …, b₉ are shown in Table 2.

OPEN IN VIEWER

Table 2.

Parameters of EHL friction coefficient calculation model.

Parameters	b 1	b 2	b 3	b 4	b 5	b 6	b 7	b 8	b 9
Value	−8.9165	1.03303	1.03608	−0.3541	2.81208	−0.1006	0.75276	−0.3910	0.62031

EHL: elasto-hydrodynamic lubrication.

Maximum Hertz contact stress can be expressed as

P h = W ′ E ′ 2 π R

(24)

In equation (24), E′ is the comprehensive elasticity modulus, GPa, and calculated as follows

E ′ = 2 ( 1 − μ s 2 ) E s + ( 1 − μ p 2 ) E p

(25)

In equation (25), μ_s and μ_p are Poisson’s ratio of driving gear and driven gear, respectively. E_s and E_p are the elastic modulus of driving gear and driven gear, respectively.

The instantaneous velocity u_s and u_p at any meshing point is calculated as follows

{ u s = ω s ρ s = ω s ( s + XA ) u p = ω p ρ p = ω p ( XY − s − XA )

(26)

In equation (26), ρ_s and ρ_p are the curvature radius. s is the distance between the instantaneous meshing point and the actual meshing starting point, m.

The sliding roll ratio Sr, the tooth surface relative sliding velocity u_e, the entrainment velocity u_sp, and the rolling velocity u_r are calculated in equation (27)

{ Sr = u e u sp = 2 u e u r = 2 | u s − u p | u s + u p u e = | u s − u p | u r = u s + u p

(27)

Here, regardless of the influence of roughness, RMS value of roughness R_rms = 0.0 μm. Mobil dte24 anti-wear lubricating oil is selected, and its relevant parameters are shown in Table 3.

OPEN IN VIEWER

Table 3.

Lubricating oil parameters.

Parameters	Values	Unit
Lubricating oil initial density ρ0	871	kg/m3
Lubricating oil initial viscosity v0	0.0527	Pa/s
Lubricating oil kinematic viscosity v	31.5	cSt
Ree–Eyring fluid viscosity v*	0.824	Pa/s
Lubricating oil viscosity pressure coefficient χ	2.389 × 108	m2/N
Ree–Eyring fluid shear stress τ0	1.0 × 106	Pa
Elastic modulus Es, Ep	227 × 109	Pa
Poisson ratio μs, μp	0.3

In order to compare the calculation results of different friction coefficient models, the power is 100 kW and the input speed n_p is 5000 r/min. The change curve of tooth surface friction coefficient with slip roll ratio sr is shown in Figure 6. Figure 6(a) and (b) shows the curve of sliding friction coefficient with the entrainment velocity at 1.0, 2.0, and 3.0 m/s. EHL model comprehensively considers the load distribution, relative sliding velocity of tooth surface, gear surface morphology, and lubrication conditions. In this article, the friction coefficient calculation model based on EHL is used. From Figure 6, due to the calculation of EHL model without considering surface roughness, and assuming that the lubrication state is complete EHL, the value of sliding friction coefficient obtained by EHL model is smaller than that of Coulomb model. The coefficient of sliding friction calculated by EHL model increases first and then decreases with the decrease in sliding roll ratio Sr. The coefficient of sliding friction decreases with the increase in entrainment velocity u_sp.

OPEN IN VIEWER

Figure 6.

Change curve of friction coefficient with slip roll ratio Sr: (a) u_sp = 1.0 m/s, (b) u_sp = 2.0 m/s, and (c) u_sp = 3.0 m/s.

Calculation of time-varying mesh stiffness based on LTCA

LTCA method is used to calculate the time-varying mesh stiffness, and the results are compared with those of Ishikawa method. LTCA model of driving gear and driven gear is shown in Figure 7.

OPEN IN VIEWER

Figure 7.

LTCA model.

Under the load P, the driving gear goes through an approach Z. Due to the tooth deformation, contact load will become distributed. After contact deformation, the state can be described by the following equation for the tooth pair k

{ ∑ j = 1 n P j I + ∑ j = 1 n P j Π = P ( P j ( j = 1 , 2 , … , n ) ) [ F ] k [ P j ] k + [ w ] k = [ Z ] + [ d j ] k

(28)

In equation (28), j_k is a point along the relative principal direction, k = I, II, [F] _k is the normal flexibility coefficient matrix of the tooth pair k, [d_j] _k = [d_1k, d_2k, …, d_ik, …, djk_jk, d_nk]^T is the tooth surface clearance to contact the discrete point j of the ellipse long axis of the tooth pair k. [Z] = [1, 1, …, 1, …, 1]^T is the Normal displacement of gear teeth under load. [P_j] _k = [P_1k, P_2k, …, P_ik, …, Pjk_jk, P_nk]^T is the normal load to contact the discrete point j of the ellipse long axis of the tooth pair k. If p_jk > 0, then d_jk = 0. if p_jk = 0, then d_jk > 0. The tooth surface initial clearance w of gears I and II is [w] _k = [δ_T] _k +[b] _k (k = I, II). [w] _k = [w_1k, w_2k, …, w_ik, …, wjk_jk, w_nk]^T, [b] _k =[b₁, b₂, …, b_i, …, bn_n]^T, [δ_T] _k =[1, 1, …, 1, …, 1]^T. n is the number of discrete points. b_j is the tooth surface normal clearance (j = 1, …, n). δ_T is the geometric transmission error.

LTCA calculation is carried out for all the meshing positions in a meshing cycle, and Z, [P], and [d] are obtained. Z is the normal linear displacement transmission error of the gear after deformation under the normal load P of the current contact position, which is transformed into the angular displacement error Δφ_sp.

The elastic deformation of the gear teeth under normal load, including the geometric bending deformation δ_B of the gear teeth, the elastic deformation δ_M at the root of the gear teeth, and the contact deformation δ_C at the meshing point, is obtained. The relationship between the comprehensive elastic deformation Δφ_sp[T_sp(k)] and the transmission torque T_sp(k) of the gear pair is as follows

Δ φ sp [ T sp ( k ) ] = T sp ( k ) | T sp ( k ) | ( a + b T sp ( k ) + c T sp ( k ) 2 3 )

(29)

In equation (29), a, b, and c are the constant terms in the formula, which are determined by the geometric parameters of the gear and the coordinates of the meshing position, T_sp(k) (k = 1,2, …,5) is the torque of the kth meshing position in a meshing cycle. δ_B = a, δ_M = bT_sp(k), and δ_C = cT_sp(k)^2/3.

Through the finite element calculation and the gear geometric contact simulation analysis program, when the gear pair is in the meshing position k, solving equation (27), the load transmission error Δδ_T (0.1T_sp(k)), Δδ_T (0.5T_sp(k)), and Δδ_T (0.9T_sp(k)) are obtained. The results are brought into equation (28), and the coefficients a, b, and c are obtained, and the function relation between the loaded deformation of the driving gear and driven gear and the nominal load is obtained. The time-varying mesh stiffness of driving gear and driven gear is obtained. The following equation (30) is satisfied with the integrated time-varying mesh stiffness K_sp of driving gear and driven gear with the corresponding angular error Δφ_sp(T_sp(k))

K sp = T sp ( k ) r bs 2 Δ φ sp [ T sp ( k ) ]

(30)

The discrete value of the time-varying mesh stiffness calculated by LTCA is expanded into a harmonic periodic function with ω as the fundamental frequency by Fourier series. Its mathematical form can be described as follows

K sp ( t ) = K m + ∑ a = 1 l K a sin ( a ω t + ψ α )

(31)

In equation (31), K_m is the average value of the fluctuation of the meshing stiffness. ω is the meshing frequency of the gear pair, and the meshing period is T=2π/ω. l is the order of the harmonic term expanded by the Fourier series. K_a is the stiffness amplitude of the ath (a = 1, 2, 3, …, l) harmonic component. ψ_α is the initial phase angle of meshing of the ath (a = 1, 2, 3, …, l) harmonic component.

The calculation formula of Ishikawa method is that the contact gear is equivalent to the cantilever beam composed of rectangle and trapezoid to calculate the deformation, and the dangerous section of the gear is defined as the length of rectangle, the specific calculation will not be described here. Figure 8(a) shows the comparison curve between LTCA calculation results and Ishikawa calculation results. The time-varying mesh stiffness including the friction factor of tooth surface is calculated, as shown in Figure 8(b). The time-varying mesh stiffness of the gear considering the friction of the tooth surface is larger than that without considering the friction, and the variation range is about 8%. This is because when the friction factor of tooth surface is considered, the friction makes the meshing stiffness of tooth pair increase in the meshing stage and decrease in the meshing out stage. The larger the friction coefficient is, the greater the variation of stiffness will be.

OPEN IN VIEWER

Figure 8.

Time-varying meshing stiffness curve: (a) without tooth surface friction and (b) with tooth surface friction.