Sage Journals: Discover world-class research

Abstract

Based on the planetary gear transmission system considering the coupling effects of friction and elastohydrodynamic lubrication, a torsional dynamic model considering friction, oil film, time-varying meshing stiffness, meshing damping, and gear backlash is established. The Runge–Kutta numerical method is used to solve the vibration equation of the system. The bifurcation diagram and largest Lyapunov exponent are used to analyze the dynamic characteristics of the system under different bifurcation parameters such as the excitation frequency, lubricant viscosity, sun–planet backlash, and planet–ring backlash. The numerical results demonstrate that with the increase of excitation frequency, the system exhibits rich nonlinear dynamic characteristics such as short-period motion, long-period motion, and chaotic motion. With the increase of lubricant viscosity, the chaotic motion of the system is suppressed at low excitation frequency and the periodic motion of the system increases at high excitation frequency. With the increase of sun–planet backlash, the chaotic motion of the system increases at low excitation frequency, and the bifurcation characteristics become complicated at high excitation frequency and enters chaotic motion in advance. With the increase of ring–planet backlash, the system delays into chaotic motion at low excitation frequency and bifurcates from single-period motion to multi-period motion in advance at high excitation frequency.

Keywords

Planetary gear transmission system mixed elastohydrodynamic lubrication friction oil film nonlinear dynamic characteristic bifurcation chaos Runge–Kutta method

Introduction

The gear system is one of the most widely used power transmission devices. In engineering practice, there will be a squeeze oil film between the meshing gear pairs. The coupling effect of oil film and friction constitutes a strong nonlinear excitation factor, which affects the vibration characteristics of the gear system. For planetary gear system with complicated structure, its vibration characteristics will be more complicated. Therefore, it is necessary to study the nonlinear dynamic characteristics of the gear system under the coupled effects of friction and lubrication.

For the study of the nonlinear characteristics of the gear system, in order to simplify the model and calculation, the lubrication factors are often ignored, and main factors considered are friction, time-varying meshing stiffness, backlash, and comprehensive error. For the single pair gear model, Kahraman and Singh¹ first coupled the nonlinear factors of clearance and stiffness to analyze the dynamic characteristics of the single pair gear system. Chen and Tang² further refined the backlash factors, and used constant backlash, time-varying backlash, and random backlash as analysis parameters. However, only the amplitude–frequency response and time–domain curve of the system were obtained, and the bifurcation characteristics of the system were not studied. With the development of dynamics theory, the friction was considered in the model. Liu et al.³ concluded that sliding friction can reduce the dynamic transmission error of the system based on the spur gear model considering sliding friction and illustrated by the bifurcation diagram that reducing the pinion speed can enhance the periodic motion of the system. Fang et al.⁴ analyzed the influence of random load and friction on the transient characteristics by Monte Carlo method. Han and Qi⁵ introduced the calculation method of friction excitation and time-varying meshing stiffness based on the 8 DOF helical gear model and obtained the dynamic response under steady state through numerical calculation. Wang et al.⁶ established a bending-torsion coupled spur gear system model considering factors such as stiffness fluctuation, transmission error, and backlash through an improved stiffness calculation model. The nonlinear analysis shows that the system exhibits different periods, subharmonic and chaotic motion under the coupling of bending and torsional vibrations. Xu and Qin⁷ established the vibration analysis model of the flexible spur ring gears by using the two-dimensional Timoshenko beam theory, and obtained that the vibration characteristics of the ring gear are different under different connection modes, providing reference for the dynamic characteristics analysis of the planetary gear. Aiming at the study of the dynamic characteristics of the planetary gear system, Li and Zhu⁸ used Poincaré-like cell-to-cell mapping method to obtain that the planetary gear system with certain parameters has four coexisting periodic orbits, while the shooting method does not have the ability to search for coexisting periodic tracks globally. Gao and Zhang⁹ analyzed the effects of layout parameters on the tooth surface friction and time-varying meshing stiffness based on the planetary gear torsion model. It was found that due to the existence of friction, the system jumps ahead and the gear’s impact motion increases. Sheng et al.¹⁰ further improved the model to increase the geometric eccentricity error of the gear and the bearing backlash and qualitatively studied the motion state of the system under different bifurcation parameters. Xiang et al.¹¹ established a nonlinear torsional model of a multistage gear transmission system consisting of planetary gears and two parallel gear stages, and concluded that the system has rich nonlinear dynamic characteristics such as periodic motion and nonperiodic motion and chaotic motion. Hou and Cao¹² added the torsional effect of the rotor on the basis of the planetary gear system and used the bifurcation diagram and the largest Lyapunov exponent (LLE) to explain the torsional nonlinear response of the system. Zhu et al.¹³ studied the effect of backlash on an encased differential planetary gear system. The results show that the smaller backlash causes bilateral shock of meshing teeth, which leads to chaotic motion of the system. The system is in stable quasiperiodic motion with the bigger backlash. The backlash change of any transmission stage has significant influence on gear meshing force and motion state of its own stage.

Lubrication factors are closely related to the dynamic characteristics of gear. Liu et al.¹⁴ reviewed the effect of lubricant oil on the efficiency, contact fatigue, and dynamics of gear systems, and concluded that lubricant oil with a relatively high viscosity can reduce the vibration and noise. Starting from the design of the gear body structure, the gear meshing vibration can also be reduced. Ramadani et al.^15,16 designed the gear body as a cellular lattice structure and reduced the stress concentration area through topology optimization. Finally, it was verified by experiments that the lattice structure can effectively reduce gear vibration, and the gaps in the lattice structure can also be filled polymer, which can further reduce vibration. The research of gear lubrication dynamics mainly focuses on single pair of meshing gears. Houser et al.¹⁷ experimentally studied the influence of friction on the vibration and noise of the gearbox, and concluded that the friction force dominates the load transmitted to the casing, and the oil viscosity and surface roughness have a great influence on the system vibration performance. Brancati et al.¹⁸ added the oil film damping force to the gear dynamic model by considering the damping lubricant effect in the meshing gear backlash and concluded that the existence of lubricant reduced the high-frequency vibration of gear. Theodossiades et al.¹⁹ studied the dynamic response of gears under the combined action of gear oil film force and bearing oil film force considering the effect of lubricant film and obtained that the hydrodynamic lubricant film behaves as a time-varying nonlinear spring-damper element. Zhang et al.²⁰ calculated the stiffness and damping of the oil film based on the line contact vibration model of the elastohydrodynamic lubrication (EHL). The results show that the viscosity effect of the lubricant can indeed improve the vibration performance compared with dry contact, and the contribution of the oil film stiffness to the total stiffness is small.

It can be seen from the above studies that the research model for gear lubrication is mainly the single pair of gear pairs, but there are few planetary gear models that consider the influence of lubrication factors. At the same time, the coupling effect of friction dynamics and lubrication has not attracted enough attention. It is rare to analyze the influence of lubrication factors on bifurcation and chaos characteristics in planetary gear system. Based on the above analysis, in order to make the planetary gear system model close to reality and improve the accuracy of the analysis results, both lubrication and friction should be considered in the dynamic model. Because the viscosity of the oil film affects the stiffness and damping of the system, the dynamic characteristics of the gear system will be more complicated. In this paper, the dynamic model of planetary gear system considering friction, oil film, time-varying meshing stiffness, meshing damping, and backlash is established and the system equation is solved by Runge–Kutta numerical method. Finally, the bifurcation and chaos characteristics of the system are analyzed by bifurcation diagram and LLE. At the same time, the motion sate of system is explained by combining the Poincaré map, phase diagram, time history, and fast fourier transform (FFT) spectrum.

The dynamic model of the system

The torsional dynamics model of the planetary gear system is established by using the concentrated mass method, as shown in Figure 1. In this model, the meshing gears are involute spur gears, and the meshing tooth surfaces are connected by spring-damping system. The mass and moment of inertia of each planet gear are the same, and parameters such as backlash, comprehensive meshing error, and damping ratio of the same type of members are also the same. The model consists basic components such as planet carrier ( $c$ ), ring gear ( $r$ ), sun gear ( $s$ ), and planet gear ( $n$ ), and the torsional displacements of the corresponding components are represented by $u_{c}$ , $u_{r}$ , $u_{s}$ , and $u_{n}$ , respectively, ( $n = 1, 2, \dots, N$ ).

Figure 1.

Dynamic model of planetary gear transmission system.

Dynamic equations of the system

Based on the literature,²¹ the coupling effects of friction and lubrication are further considered. The torsional dynamic vibration equation of the planetary gear system can be established by Newton’s law as follows

{\begin{cases} (\frac{I_{c}}{r_{c}^{2}} + N m_{n}) {\ddot{u}}_{c} - \sum_{n = 1}^{N} F_{rn} \cos (α_{rn}) - \sum_{n = 1}^{N} F_{sn} \cos (α_{sn}) = \frac{T_{out}}{r_{c}} \\ \frac{I_{r}}{r_{r}^{2}} {\ddot{u}}_{r} + \sum_{n = 1}^{N} F_{rn} - \sum_{n = 1}^{N} f_{rn} \frac{L_{rn}}{r_{r}} = 0 \\ \frac{I_{s}}{r_{s}^{2}} {\ddot{u}}_{s} + \sum_{n = 1}^{N} F_{sn} - \sum_{n = 1}^{N} f_{sn} \frac{L_{sn}}{r_{s}} = \frac{T_{in}}{r_{s}} \\ \frac{I_{n}}{r_{n}^{2}} {\ddot{u}}_{n} + F_{sn} - f_{sn} \frac{L_{sn}}{r_{s}} - F_{rn} + f_{rn} \frac{L_{rn}}{r_{r}} = 0 \end{cases}

(1)

Let $i = s n, r n$ , $j = s, c, r, n$ , then $s n$ is the internal meshing pair of sun–planet, $r n$ is the external meshing pair of ring–planet, $k_{i} (t)$ is equivalent stiffness, $e_{i} (t)$ is comprehensive meshing error, $c_{i}$ is equivalent damping, $F_{i}$ is dynamic engagement force, $f_{i}$ is friction force, $L_{i}$ is friction arm, $m_{j}$ , $I_{j}$ , and $r_{j}$ are mass, moment of inertia, and base circle radius, respectively, $T_{in}$ is input torque, $T_{out}$ is load torque.

The dynamic engagement force of the gear pair is related to the system damping and stiffness, and the equation of the time-varying dynamic engagement force can be expressed as

F_{i} = k_{i} (t) f (δ_{i}) + c_{i} f ({\dot{δ}}_{i})

(2)

where

f (δ_{i})

is the nonlinear function of backlash.

Under the premise of not affecting the results, the eccentricity error of the gear is ignored in the model, so the comprehensive meshing error is the static meshing error. The static transmission error of the gear mainly comes from the base joint error and tooth shape error of the gear pair, which can be expressed as the sine function with the meshing frequency as follows²²

e_{i} (t) = E_{i} \sin (w_{m} t + γ_{i})

(3)

where

E_{i}

is the magnitude of the meshing error,

γ_{i}

is the initial phase of the meshing error,

w_{m}

is the meshing frequency.

Considering that there may be three meshing states of each gear pair in the planetary gear transmission, namely normal meshing, empty meshing, and tooth back meshing, the segmental function $f (δ_{i})$ is used to describe the gear backlash,²³ as shown in Figure 2.

Figure 2.

Nonlinear function of backlash.

The expression of the relative displacement of the gear on the meshing line can be given by

f (δ_{i}) = {\begin{cases} δ_{i} - b_{i} δ_{i} > b_{i} \\ 0 | δ_{i} | \leq b_{i} \\ δ_{i} + b_{i} δ_{i} < - b_{i} \end{cases}

(4)

where

δ_{i}

is the relative displacement of the meshing point of the gear pair,

b_{i}

is half of the gear backlash.

Equivalent stiffness of meshing gear pairs considering squeeze oil film

According to the definition, the load applied by one or more pairs of gear teeth engaged on the meshing line to produce a unit deflection on the unit tooth width is called meshing stiffness.

Under the condition of mixed EHL, the deflection of spur gears under external force is the sum of gear elastic deformation and oil film deformation, so the equivalent meshing stiffness of spur gears is given by²⁴

k_{i} (t) = \frac{W}{δ_{e}} = \frac{W}{δ_{m} + δ_{l}} = \frac{W}{W / k_{m} (t) + W / k_{J Y} (t)} = \frac{1}{1 / k_{m} (t) + 1 / k_{J Y} (t)}

(5)

where

W

is the unit tooth width load,

δ_{e}

is the sum of the gear elastic deformation and the oil film deformation,

δ_{m}

is the gear elastic deformation,

δ_{l}

is oil film deformation,

k_{m} (t)

is the meshing stiffness of the gear pair,

k_{J Y} (t)

is the stiffness of the squeeze oil film.

In actual working conditions, when the gear tooth vibrate out of the stable meshing state, there is a squeezing effect due to the presence of the oil film in the backlash to reduce vibration. For the gear lubrication model, the meshing gear pair can be assumed to be in contact with two cylinders, and then simplified to the contact between cylinder with an equivalent radius of curvature and rigid plane, as shown in Figure 3.

Figure 3.

Elastohydrodynamic lubrication (EHL) model of gear pair.

Due to the entrainment effect, there is a squeeze oil film between the contact surfaces of the meshing gear pair. Assuming that the lubricating oil is an incompressible fluid, the basic equation of the squeeze oil film of the gear can be expressed as^25,26

t = \frac{A_{c} η L_{T}^{4}}{F_{J Y}} (\frac{1}{h^{2}} - \frac{1}{h_{\min}^{2}})

(6)

where

h

is the thickness of the squeeze oil film,

t

is the time required to reduce the film thickness from

h

h_{\min}

A_{c}

is a constant,

η

is the lubricant viscosity,

L_{T}

is the length of the area at the contact point and equal

{\begin{cases} L_{T} = \sqrt[4]{2 a^{2} \cdot B^{4} / (16 a^{2} + B^{4})} \\ a = \sqrt{8 F_{N} R (1 - ρ^{2}) / π B E} \end{cases}

(7)

where

F_{N}

is the normal force of the tooth profile along the direction of the meshing line,

B

is the tooth width,

R

is the comprehensive curvature radius,

ρ

is the Poisson’s ratio, and

E

is the elastic modulus.

The thickness of squeeze oil film can be expressed as

h = b_{0} / 2 - (r_{b 1} θ_{1} - r_{b 2} θ_{2})

(8)

where

b_{0}

is the average meshing backlash,

r_{b 1}

and

r_{b 2}

are the radius of the base circle of the equivalent gear,

θ_{1}

and

θ_{2}

are the angular displacements of the equivalent gear.

The basic equation (6) of squeeze oil film is derived and the squeezing force of the oil film can be expressed as

F_{J Y} = - 2 A_{c} η L_{T}^{4} v / h^{3}

(9)

then the oil film stiffness can be obtained by

k_{J Y} = 6 A_{c} η L_{T}^{4} v h^{- 4}

(10)

where

v

is the extrusion speed.

The meshing stiffness of the gear pair is obtained by using the Ishikawa formula. The time-varying meshing stiffness is developed by Fourier series with the meshing frequency as the fundamental frequency. In order to facilitate coupling with the oil film stiffness, the meshing stiffness only takes the harmonic fundamental frequency part, which can be expressed as

k_{m} (t) = k_{mav} (1 + k \cos (w_{m} t + ϕ))

(11)

where

k_{mav}

is the average meshing stiffness,

w_{m}

is the meshing frequency,

k

is the stiffness fluctuation coefficient,

φ

is the initial phase of meshing stiffness.

Since the oil film stiffness is also time-varying, the equivalent stiffness of the system can be obtained by the gear meshing stiffness and oil film stiffness in series when considering the coupling effect of lubrication. In order to illustrate the effect of the squeeze oil film on the stiffness, the stiffness comparison curves of the meshing pairs with and without the squeeze oil film are shown in Figure 4.

Figure 4.

Stiffness comparison curves.

It can be seen from Figure 4 that due to the continuous change of the single and double teeth of the gear, the stiffness of the gear changes periodically, which forms an internal excitation. The stiffness of the double-toothed region is greater than that of the single-toothed region. Compared with the contact time-varying meshing stiffness of the gear pair without oil film, the equivalent stiffness of the internal meshing pair with oil film and the external meshing pair with oil film both decrease, which is consistent with the conclusion of literature.²⁷

Equivalent damping of meshing gear pairs considering squeeze oil film

In the gear system with coupling effect of lubrication, the system has oil film damping due to the viscous effect of the oil film. When considering the lubrication effect, the damping is treated in the same way as the stiffness and the equivalent damping of the system can be obtained from oil film damping and meshing damping in series,²⁷ as shown in Figure 5.

Figure 5.

Equivalent damping and stiffness models of gear system under lubrication effect.

The equivalent damping of the gear system considering the squeeze oil film can be expressed as

c_{i} = \frac{1}{1 / c_{o} + 1 / c_{m}} = \frac{c_{0} c_{m}}{c_{0} + c_{m}}

(12)

where

c_{m}

is the meshing damping,

c_{0}

is the oil film damping

For the nonsteady state problem of pure squeeze lubrication, the Reynold²⁸ differential equation under isothermal conditions can be simplified as

\frac{\partial}{\partial \bar{ξ}} (h^{3} \frac{\partial p}{\partial \bar{ξ}}) = 12 \bar{η} \frac{\partial h}{\partial \bar{t}}

(13)

where

p

is the oil film pressure,

\bar{ξ}

is the relative rotation direction.

With $p = 0$ , $\bar{ξ} = \pm a$ , $\partial p / \partial \bar{ξ} = 0$ , and $\bar{ξ} = 0$ as the limits, the oil film pressure $p (\bar{ξ})$ is obtained by the integral transformation of equation (13). The oil film damping force can be obtained by performing a second integral transformation on the oil film pressure $p (\bar{ξ})$ in the high-pressure contact area $(- a, a)$ , which can be expressed as

F_{0} = c_{0} \cdot \partial h / \partial \bar{t} = - 8 a^{3} \bar{η} l / h^{3} \cdot \partial h / \partial \bar{t}

(14)

where

l

is the width of the oil film in the squeeze area,

\partial h / \partial \bar{t}

is the change in oil film thickness per unit time.

Finally, the expression of oil film damping is obtained by

c_{0} = - 8 a^{3} \bar{η} l / h^{3}

(15)

The meshing damping of the gear pairs can be expressed as¹⁰

c_{m} = 2 ζ \sqrt{k_{mav} / (1 / M_{p} + 1 / M_{q})}

(16)

where

ζ

is the meshing damping ratio,

k_{mav}

is the average meshing stiffness,

M

is the equivalent mass and

M = I / r^{2}

The friction force based on mixed EHL

During the gear transmission, the friction coefficient of the tooth surface changes significantly with the speed, load distribution, and tooth surface profile. Martin²⁹ found that the state of gear lubrication is constantly changing between full oil film EHL and boundary lubrication. In fact, the lubrication form of the gear meshing pairs is mixed EHL, which is composed of full oil film EHL and boundary lubrication.

The mixed EHL friction force of the meshing pair can be introduced as

f_{i} = λ μ F_{i}

(17)

where

μ

is the friction coefficient of mixed EHL,

F_{i}

is the dynamic meshing force of the gear pair,

λ

is the direction coefficient of the friction force^30,31 and equal

λ = sign (x) = {\begin{cases} 1 (x > 0) \\ 0 (x = 0) \\ - 1 (x < 0) \end{cases}

(18)

let

x = L_{p} - R_{b p} \tan (α)

, where

L_{p}

is the time-varying friction arm of the driving gear

p

R_{b p}

is the radius of base circle of the driving gear p,

α

is the pressure angle.

In order to obtain the friction force of mixed EHL, the friction coefficient of EHL needs to be determined. The friction coefficient of mixed EHL is related to the load distribution percentage function $ξ (0 < ξ < 1)$ . $ξ$ is defined as the function of the full oil film EHL load percentage, then the boundary lubrication load percentage function is $(1 - ξ)$ . Based on a large number of experiments and numerical calculations,^32–34 the load distribution percentage function and oil film thickness in the mixed lubrication state are described as

{\begin{cases} ξ = 1.21 {(h_{min} / σ)}^{0.64} / [1 + 0.37 {(h_{min} / σ)}^{1.26}] \\ h_{min} = 1.6 α_{p}^{0.6} {(η_{0} v_{m})}^{0.7} R^{0.43} E^{0.43} / {\bar{w}}^{0.13} \end{cases}

(19)

where

h_{min}

is the minimum oil film thickness,

σ

is the root mean square of roughness,

α_{p}

is the viscosity coefficient,

v_{m}

is the average speed of tooth surface meshing,

\bar{w}

is the load coefficient.

In the mixed lubrication state, the tooth surface contact load $P_{total}$ is composed of full oil film EHL pressure $P_{EHL}$ and boundary lubrication pressure $P_{B}$ . Similarly, the total tangential force $F_{total}$ is also divided into the sum of these two parts.

{\begin{cases} P_{total} = P_{EHL} + P_{B} \\ F_{total} = F_{EHL} + F_{B} \end{cases}

(20)

The mixed friction coefficient can be expressed as

μ = F_{total} / P_{total}

(21)

According to the full oil film EHL friction coefficient is $μ_{EHL}$ and the boundary lubrication friction coefficient is $μ_{B}$ , the tangential force can be given by

{\begin{cases} F_{EHL} = μ_{EHL} P_{EHL} = μ_{EHL} P_{total} ξ \\ F_{B} = μ_{B} P_{B} = μ_{B} P_{total} (1 - ξ) \end{cases}

(22)

Finally, the friction coefficient of mixed EHL can be introduced as

μ = ξ μ_{EHL} + (1 - ξ) μ_{B}

(23)

Under boundary lubrication condition, the friction coefficient is basically constant. The boundary lubrication friction coefficient measured through experiments is generally 0.1–0.2.²⁹ According to the full oil film EHL model,³⁵ the full oil film EHL coefficient can be expressed as

{\begin{cases} μ_{EHL} = \exp [f (S R, p h, η_{0}, σ)] p h^{b_{2}} {| S R |}^{b_{3}} v_{r}^{b_{6}} η_{0}^{b_{7}} R^{b_{8}} ξ^{0.2} \\ f (S R, p h, η_{0}, σ) = b_{1} + b_{4} p h | S R | \log_{10}^{η_{0}} + b_{5} \cdot \exp [- p h | S R | \log_{10}^{η_{0}}] + b_{9} \cdot \exp [σ] \end{cases}

(24)

where

S R

is the slip ratio of the gear pair,

p h

is the normal pressure,

v_{r}

is the relative sliding speed,

b_{1}

= −8.9,

b_{2}

= 1.1,

b_{3}

= 1,

b_{4}

= −0.3,

b_{5}

= 2.8,

b_{6}

= −0.1,

b_{7}

= −0.7,

b_{8}

= −0.4,

b_{9}

= −0.6.

In order to explain the influence of tooth surface friction on the dynamic engagement force of the meshing gear pair, the dynamic engagement force comparison curves of meshing gear pairs with and without friction are shown in Figure 6.

Figure 6.

Comparison curves of dynamic engagement force.

It can be seen from Figure 6 that compared with the dynamic engagement force of gear pair without friction, the system’s dynamic engagement force of the external gear pair with friction and the internal gear pair with friction both decrease, which indicates that the tooth surface friction suppresses the engagement force of the system.

Friction arm of internal and external meshing gear pairs

As shown in Figure 7, in the schematic diagram of the external meshing gear pair, $N_{1} N_{2}$ is the theoretical mesh line length, $B_{1} B_{2}$ is the actual mesh line length, $r_{s}$ and $r_{n}$ are the base circle radius of the sun gear and planet gear respectively, $r_{as}$ and $r_{an}$ are the top circle radius of the sun gear and the planet gear, respectively. Ignoring the in-gear and out-of-gear deformations when the gears are in contact, the friction arm of the external meshing gear pair can be expressed as

{\begin{cases} L_{s n} (t) = \sqrt{r_{a s}^{2} - r_{s}^{2}} - ε_{s n} S_{p b} + w_{s} r_{s} t \\ L_{n s} (t) = \sqrt{r_{a n}^{2} - r_{n}^{2}} - ε_{s n} S_{p b} + w_{s} r_{s} t \end{cases}

(25)

where

L_{s n} (t)

is the friction arm relative to the sun gear,

L_{n s} (t)

is the friction arm relative to the planetary gear,

ε_{s n}

is the coincidence degree coefficient, and

S_{p b}

is the pitch of the base circle.

Figure 7.

Schematic diagram of the external meshing gear pair.

Similarly, the friction arm of the internal meshing gear pair can be expressed as

{\begin{cases} L_{n r} (t) = \sqrt{r_{a n}^{2} - r_{n}^{2}} - ε_{r n} S_{p b} + w_{n} r_{n} t \\ L_{r n} (t) = \sqrt{r_{a r}^{2} - r_{r}^{2}} - ε_{r n} S_{p b} + w_{n} r_{n} t \end{cases}

(26)

where

L_{r n} (t)

is the friction arm relative to the ring gear,

L_{n r} (t)

is the friction arm relative to the planetary gear,

r_{a r}

is the root circle radius of the ring gear, and

r_{r}

is the base circle radius of the ring gear.

Dimensionless equations of the system

Due to the existence of gear backlash, the system constraints are incomplete. Equation (1) is a semi-positive definite system with rigid body displacement and indefinite solutions. In order to eliminate rigid body displacements and reduce the number of system equations to achieve dimensionality reduction, the relative displacements of the internal and external meshing pairs were introduced by

{\begin{cases} δ_{s n} = u_{s} - u_{c} \cos (α_{s n}) + u_{n} + e_{s n} \\ δ_{r n} = u_{r} - u_{c} \cos (α_{r n}) - u_{n} + e_{r n} \end{cases}

(27)

It should be pointed out that the introduced new coordinates are essentially obtained through bicontinuous mapping, which means that they have the same bifurcation characteristics. The coordinate transformation not only makes the system deterministic but also has intuitive physical meaning, especially in the research field of gear vibration and shock. With the effect of coordinate transformation and equation recombination, equation (1) can be rewritten as

{\begin{cases} {\ddot{δ}}_{s n} - (- \frac{1}{M_{s}} + \frac{μ_{s n} λ_{s n} L_{s n}}{M_{s} r {}_{s}} - \frac{\cos^{2} (α_{s n})}{M_{c}}) \sum_{n = 1}^{N} k_{s n} f (δ_{s n}) - (\frac{μ_{s n} λ_{s n} L_{s n}}{M_{s} r {}_{s}} - \frac{1}{M_{s}} - \frac{\cos^{2} (α_{s n})}{M_{c}}) \sum_{n = 1}^{N} c_{s n} {\dot{δ}}_{s n} \\ + \frac{\cos (α_{s n}) \cos (α_{r n})}{M_{c}} \sum_{n = 1}^{N} k_{r n} f (δ_{r n}) + \frac{\cos (α_{s n}) \cos (α_{r n})}{M_{c}} \sum_{n = 1}^{N} c_{r n} {\dot{δ}}_{r n} - (\frac{μ_{s n} λ_{s n} L_{n s}}{M_{n} r {}_{n}} - \frac{1}{M_{n}}) k_{s n} f (δ_{s n}) \\ - (\frac{μ_{s n} λ_{s n} L_{n s}}{M_{n} r {}_{n}} - \frac{1}{M_{n}}) c_{s n} {\dot{δ}}_{s n} - (\frac{1}{M_{n}} - \frac{μ_{r n} λ_{r n} L_{n r}}{M_{n} r {}_{n}}) k_{r n} f (δ_{r n}) - (\frac{1}{M_{n}} - \frac{μ_{r n} λ_{r n} L_{n r}}{M_{n} r {}_{n}}) c_{r n} {\dot{δ}}_{r n} \\ = \frac{T_{i n}}{M_{s} r {}_{s}} + \frac{T_{out} \cos (α_{s n})}{M_{c} r {}_{c}} + {\ddot{e}}_{s n} � 28 � \\ {\ddot{δ}}_{r n} - (- \frac{1}{M_{r}} + \frac{μ_{r n} λ_{r n} L_{r n}}{M_{r} r {}_{r}} - \frac{\cos^{2} (α_{r n})}{M_{c}}) \sum_{n = 1}^{N} k_{r n} f (δ_{r n}) - (\frac{μ_{r n} λ_{r n} L_{r n}}{M_{r} r {}_{r}} - \frac{1}{M_{r}} - \frac{\cos^{2} (α_{r n})}{M_{c}}) \sum_{n = 1}^{N} c_{r n} {\dot{δ}}_{r n} \\ + \frac{\cos (α_{s n}) \cos (α_{r n})}{M_{c}} \sum_{n = 1}^{N} k_{s n} f (δ_{s n}) + \frac{\cos (α_{s n}) \cos (α_{r n})}{M_{c}} \sum_{n = 1}^{N} c_{s n} {\dot{δ}}_{s n} - (\frac{1}{M_{n}} - \frac{μ_{s n} λ_{s n} L_{n s}}{M_{n} r {}_{n}}) k_{s n} f (δ_{s n}) \\ - (\frac{1}{M_{n}} - \frac{μ_{s n} λ_{s n} L_{n s}}{M_{n} r {}_{n}}) c_{s n} {\dot{δ}}_{s n} - (\frac{μ_{r n} λ_{r n} L_{n r}}{M_{n} r {}_{n}} - \frac{1}{M_{n}}) k_{r n} f (δ_{r n}) - (\frac{μ_{r n} λ_{r n} L_{n r}}{M_{n} r {}_{n}} - \frac{1}{M_{n}}) c_{r n} {\dot{δ}}_{r n} \\ = \frac{T_{out} \cos (α_{r n})}{M_{c} r {}_{c}} + {\ddot{e}}_{r n} \end{cases}

The dimensionless comparison analysis can be performed on the results at different dimension scales by introducing the dimensionless time variable $τ$ , nominal time scale $w_{n}$ , dimensionless excitation frequency $Ω$ , and letting $τ = w_{n} t$ , here

{\begin{cases} M_{1} = I_{1} / r_{1}^{2} \\ M_{s} = I_{s} / r_{s}^{2} \\ w_{n} = \sqrt{k_{s1} (1 / M_{s} + 1 / M_{1})} \\ Ω = w_{m} / w_{n} \end{cases}

(29)

Introducing the nominal displacement scale $b_{c}$ , the dimensionless displacement, velocity, acceleration, backlash, and lubricant viscosity can be expressed respectively as

{\begin{cases} {\bar{δ}}_{i} = δ_{i} / b_{c} \\ {\dot{\bar{δ}}}_{i} = δ_{i} / b_{c} w_{n} \\ {\ddot{\bar{δ}}}_{i} = δ_{i} / b_{c} w_{n}^{2} \\ {\bar{b}}_{i} = b_{i} / b_{c} \\ \bar{η} = η / η_{0} \end{cases}

(30)

The dimensionless backlash function of the system can be expressed as

f ({\bar{δ}}_{i}) = {\begin{cases} {\bar{δ}}_{i} - {\bar{b}}_{i} ​ {\bar{δ}}_{i} > {\bar{b}}_{i} \\ 0 | {\bar{δ}}_{i} | \leq {\bar{b}}_{i} \\ {\bar{δ}}_{i} + {\bar{b}}_{i} {\bar{δ}}_{i} < - {\bar{b}}_{i} \end{cases}

(31)

In order to facilitate the use of the variable-step Runge–Kutta numerical method, equation (28) in the form of the matrix vector can be expressed as

[M] \ddot{δ} + [C] \dot{δ} + [K] δ = F

(32)

where

δ

is the generalized coordinate vector,

[M]

is the inertia matrix,

[C]

is the damping matrix,

[K]

is the stiffness matrix, and

F

is the excitation vector.

Let $z_{1} = {\bar{δ}}_{s n}$ , $z_{2} = {\dot{\bar{δ}}}_{s n}$ , $z_{3} = {\bar{δ}}_{r n}$ , $z_{4} = {\dot{\bar{δ}}}_{r n}$ , the state equations of the system can be given by

{\begin{cases} {\dot{z}}_{1} = z_{2} \\ {\dot{z}}_{2} = (- \frac{1}{M_{s} w_{n}^{2}} + \frac{μ_{s n} λ_{s n} L_{s n}}{M_{s} r {}_{s}w {_{n}}^{2}} - \frac{\cos^{2} (α_{s n})}{M_{c} w_{n}^{2}}) \sum_{n = 1}^{N} {\bar{k}}_{s n} f ({\bar{δ}}_{s n}) + (\frac{μ_{s n} λ_{s n} L_{s n}}{M_{s} r {}_{s}w_{n}} - \frac{1}{M_{s} w_{n}} - \frac{\cos^{2} (α_{s n})}{M_{c} w_{n}}) \sum_{n = 1}^{N} c_{s n} {\dot{\bar{δ}}}_{s n} \\ - \frac{\cos (α_{s n}) \cos (α_{r n})}{M_{c} w_{n}^{2}} \sum_{n = 1}^{N} {\bar{k}}_{r n} f ({\bar{δ}}_{r n}) - \frac{\cos (α_{s n}) \cos (α_{r n})}{M_{c} w_{n}} \sum_{n = 1}^{N} c_{r n} {\dot{\bar{δ}}}_{r n} + (\frac{μ_{s n} λ_{s n} L_{n s}}{M_{n} r {}_{n}w {_{n}}^{2}} - \frac{1}{M_{n} w_{n}^{2}}) {\bar{k}}_{s n} f ({\bar{δ}}_{s n}) \\ + (\frac{μ_{s n} λ_{s n} L_{n s}}{M_{n} r {}_{n}w_{n}} - \frac{1}{M_{n} w_{n}}) c_{s n} {\dot{\bar{δ}}}_{s n} + (\frac{1}{M_{n} w_{n}^{2}} - \frac{μ_{r n} λ_{r n} L_{n r}}{M_{n} r {}_{n}w {_{n}}^{2}}) {\bar{k}}_{r n} f ({\bar{δ}}_{r n}) + (\frac{1}{M_{n} w_{n}} - \frac{μ_{r n} λ_{r n} L_{n r}}{M_{n} r {}_{n}w_{n}}) c_{r n} {\dot{\bar{δ}}}_{r n} \\ + \frac{T_{i n}}{M_{s} r {}_{s}w {_{n}}^{2} b_{c}} + \frac{T_{out} \cos (α_{s n})}{M_{c} r {}_{c}w {_{n}}^{2} b_{c}} + {\ddot{\bar{e}}}_{s n} \\ {\dot{z}}_{3} = z_{4} \\ {\dot{z}}_{4} = (- \frac{1}{M_{r} w_{n}^{2}} + \frac{μ_{r n} λ_{r n} L_{r n}}{M_{r} r {}_{r}w {_{n}}^{2}} - \frac{\cos^{2} (α_{r n})}{M_{c} w_{n}^{2}}) \sum_{n = 1}^{N} {\bar{k}}_{r n} f ({\bar{δ}}_{r n}) + (\frac{μ_{r n} λ_{r n} L_{r n}}{M_{r} r {}_{r}w_{n}} - \frac{1}{M_{r} w_{n}} - \frac{\cos^{2} (α_{r n})}{M_{c} w_{n}}) \sum_{n = 1}^{N} c_{r n} {\dot{\bar{δ}}}_{r n} \\ - \frac{\cos (α_{s n}) \cos (α_{r n})}{M_{c} w_{n}^{2}} \sum_{n = 1}^{N} {\bar{k}}_{s n} f ({\bar{δ}}_{s n}) - \frac{\cos (α_{s n}) \cos (α_{r n})}{M_{c} w_{n}} \sum_{n = 1}^{N} c_{s n} {\dot{\bar{δ}}}_{s n} + (\frac{1}{M_{n} w_{n}^{2}} - \frac{μ_{s n} λ_{s n} L_{n s}}{M_{n} r {}_{n}w {_{n}}^{2}}) {\bar{k}}_{s n} f ({\bar{δ}}_{s n}) \\ + (\frac{1}{M_{n} w_{n}} - \frac{μ_{s n} λ_{s n} L_{n s}}{M_{n} r {}_{n}w_{n}}) c_{s n} {\dot{\bar{δ}}}_{s n} + (\frac{μ_{r n} λ_{r n} L_{n r}}{M_{n} r {}_{n}w {_{n}}^{2}} - \frac{1}{M_{n} w_{n}^{2}}) {\bar{k}}_{r n} f ({\bar{δ}}_{r n}) + (\frac{μ_{r n} λ_{r n} L_{n r}}{M_{n} r {}_{n}w_{n}} - \frac{1}{M_{n} w_{n}}) c_{r n} {\dot{\bar{δ}}}_{r n} \\ + \frac{T_{out} \cos (α_{r n})}{M_{c} r {}_{c}w {_{n}}^{2} b_{c}} + {\ddot{\bar{e}}}_{r n} (33) \end{cases}

Numerical simulation and result analysis

The fourth to fifth order Runge–Kutta numerical method is used to solve the dynamic differential equation. In order to eliminate the transient response of the system, the first hundreds of cycles are omitted. Because the vibration characteristics of the internal and external meshing pairs are consistent, the relative displacement of the external meshing pair $δ_{s n}$ is analyzed as an example. The bifurcation diagram and LLE are used to analyze the bifurcation and chaos characteristics of the system. At the same time, the Poincaré map, phase diagram, time history, and FFT spectrum are used to illustrate the system’s motion state. The basic parameters of the planetary gear transmission system are shown in Table 1.³⁶ The benchmark parameters of the system obtained by calculating the data in Table 1 are as follows

Table 1.

Basic parameters of planetary gear system.

Parameters	Sun gear	Planet gear	Ring gear	Planet carrier
Number of teeth, Z	20	29	79	–
Module, m/mm	2.25
Base circle diameter, d/mm	42.29	61.31	167.03	113
Moment of inertia, I/ (kg.mm²)	146	297	17440	11178
Pressure angle, $α$ /（°）	20
Tooth width, B/mm	25	28	25	25
Mass, M/kg	0.53	0.35	1.85	4.5
Modulus of elasticity, E/Gpa	206
Poisson’s ratio, P	0.3
RMS value of roughness, $σ$ /um	0.06
Ambient viscosity of lubricant, $η_{0}$ /pa·s	0.01
Viscosity pressure coefficient of lubricant, $α_{p}$ /pa⁻¹	2.19 × 10⁻⁸
Load coefficient, $\bar{w}$	0.3
Addendum coefficient, $h_{a}$	1
Top backlash coefficient, $c_{b}$	0.25
Number of planets, $N$	3

$k$ = 0.25, $φ$ = 0, $k_{s n} = 4.37 \times 10^{8}$ N/m, $k_{r n} = 5.65 \times 10^{8}$ N/m, $E_{i}$ = 3 µm, $b_{i}$ = 50 µm, $γ_{i}$ = 0, $ζ$ = 0.088, $T_{i n}$ = 100 $N \cdot m$ .

For the gear system that considering the squeeze oil film, it can be seen from the basic equations (6)–(8) of the squeeze oil film that the parameters, such as excitation frequency, viscosity, and gear backlash, are related to the oil film thickness, then the oil film thickness will affect the oil film stiffness and oil film damping, which will affect the system’s dynamic response. In order to study the effect of squeeze oil film on the system, it is necessary to analyze the dynamic characteristics of the system from the perspective of excitation frequency, lubricant viscosity, and backlash.

To further clarify the effect of squeeze oil film on the dynamic response of the system, the time–domain curves of dimensionless relative vibration displacement of internal and external meshing gear pairs are obtained without and with the effect of squeeze oil film, as shown in Figure 8.

Figure 8.

Time–domain comparison curve of dimensionless relative vibration displacement.

It can be seen from Figure 8 that when the oil film is not considered, the dimensionless relative vibration displacement of the external gear pair changes in the range of [0.747, 1.256], and the dimensionless relative vibration displacement of the internal gear pair changes in the range of [0.638, 1.176]. When the oil film is considered, the dimensionless relative vibration displacement of the external gear pair changes in the range of [0.697, 1.206], and the dimensionless relative vibration displacement of the internal gear pair changes in the range of [0.593, 1.102]. After considering the effect of the oil film, the dimensionless relative vibration displacement amplitude and fluctuation range of the internal and external meshing gear pairs decrease, which is consistent with the conclusion in the literatures.^37,38 Through the above analysis, the correctness of the model and simulation results can be verified, and it also provides support for the nonlinear analysis later.

Nonlinear dynamics characteristic of the system with excitation frequency

The dimensionless gear backlash ${\bar{b}}_{i}$ is set to 1. When the dimensionless excitation frequency $Ω$ is considered as the bifurcation parameter, the bifurcation diagram and LLE of the system are shown in Figure 9.

Figure 9.

Bifurcation diagram and LLE of the system with $Ω$ .

It can be seen from Figure 9 that when dimensionless excitation frequency $Ω$ changes in the range of [0, 4], the system exhibits rich bifurcation characteristics, in which multiple bifurcation points and jump points appear. In the low-frequency region, the system first undergoes periodic motion, then enters chaotic motion via Hopf bifurcation, and then moves from chaotic motion to periodic motion through inverse doubling bifurcation. In the high-frequency region, the system moves from periodic motion to chaotic motion again through jump shock and finally changes from chaotic motion to periodic motion through inverse doubling bifurcation. In the whole excitation frequency region, the motion state of the system is mainly composed of two chaotic regions (LLE value is positive) and several separated periodic regions (LLE value is negative). At the same time, the two chaotic regions have narrow periodic regions.

In order to track the process and evolution of the system’s motion with the dimensionless excitation frequency $Ω$ , Poincaré map, phase diagram, time history, and FFT spectrum are obtained to analyze the nonlinear characteristics of the system in detail.

With the increase of the dimensionless excitation frequency $Ω$ , the system’s starting state of motion is 1 T-periodic motion; when $Ω$ increases to 0.76, the system enters 2 T-periodic motion through the Hopf bifurcation, as shown in Figure 10. When $Ω$ = 0.76, there are two discrete points on the Poincaré map. The phase diagram shows a closed curve formed by the phase trajectory around two cycles. The time history shows a cycle curve of two cycles. The peaks of the FFT spectrum are distributed at discrete points of $A Ω_{B} / 2$ (where: $A$ is a positive integer; $Ω_{B}$ is the fundamental frequency). When $Ω$ reaches 0.92, the system performs 4 T-periodic motion, as shown in Figure 11. When $Ω$ = 0.92, there are four discrete points on the Poincaré map. The phase diagram is a closed pattern formed by the phase trajectory around four cycles. The time history is a cycle graph that appears every four cycles. The peaks of the FFT spectrum are distributed at the discrete points of $A Ω_{B} / 4$ .

Figure 10.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when $Ω$ = 0.76.

Figure 11.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when $Ω$ = 0.92.

When $Ω$ increases to 0.98, the system enters chaotic motion through jump shock; when $Ω$ continues to increase to 1.02, the system enters 2 T-periodic motion; when $Ω$ reaches 1.14, the Hopf loop ruptures and the system enters chaotic motion, as shown in Figure 12. For example, when $Ω$ = 1.57, there is a strange attractor with a fractal structure on the Poincaré map, which is distributed in a cloud-like discrete point set in the certain area. Phase trajectories are intertwined and not closed, tending to fill a portion of the phase space. The time history shows an irregular aperiodic curve. In addition to the spikes corresponding to the periodic motion, the FFT spectrum also has a clear noise background. When $Ω$ increases to 1.68, the system enters the 2 T-periodic motion again. For example, when $Ω$ = 1.71, The 2 T-periodic motion of the system is shown in Figure 13. When $Ω$ is 1.8, the system enters 1 T-periodic motion, as shown in Figure 14. For example, when $Ω$ = 1.89, there is a discrete point on the Poincaré map. The phase diagram is a standard ellipse composed of phase trajectories around one cycle. The time history is a period curve, which appears every other period. The peaks of the FFT spectrum are distributed at discrete points of $A Ω_{B}$ .

Figure 12.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when $Ω$ = 1.57.

Figure 13.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when $Ω$ = 1.71.

Figure 14.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when $Ω$ = 1.89.

When $Ω$ belongs to the range of [2, 2.18], the system experiences several jump shocks, but it always moves for 1 T-periodic motion; when $Ω$ continues to increase to 2.2, the system changes from 1 T-periodic motion to chaotic motion, accompanied by a narrow 2 T-periodic motion in this chaotic region; when $Ω$ reaches 2.66, the system moves from chaotic motion to 8 T-periodic motion through inverse doubling bifurcation, as shown in Figure 15. For example when $Ω$ = 2.665, there are eight discrete points on the Poincaré map. The phase diagram is a closed curved ring formed by the phase trajectory around eight cycles. The time history is a cycle curve that cycles every eight cycles. The peaks of FFT spectrum are distributed at the discrete points of $A Ω_{B} / 8$ . When $Ω$ is 2.68, the system moves in 4 T-periodic motion. For example when $Ω$ = 2.685, The 4 T-periodic motion is shown in Figure 16. When $Ω$ increases to 2.7, the system moves in 2 T-periodic motion; when $Ω$ continues to increase to 2.94, the system finally moves to 1 T-periodic motion by the stable attractor.

Figure 15.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when $Ω$ = 2.665.

Figure 16.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when $Ω$ = 2.685.

In engineering practice, the excitation frequency corresponds to the speed of the system. In order to make the gear system have better stability and reliability, and extend the operating cycle life of the system, the speed range corresponding to chaotic motion and the critical speed (the point of bifurcation) that changes the state of the system must be avoided.

Nonlinear dynamics characteristics of the system with lubricant viscosity

The increase of lubricant viscosity will cause changes in friction coefficient, oil film damping, and oil film stiffness, which will affect the dynamic characteristics of the system. It is necessary to analyze the effect of lubricant viscosity on the system. The dimensionless excitation frequency $Ω$ is set to 1, the dimensionless gear backlash ${\bar{b}}_{i}$ is set to 1, and other parameters are unchanged. When the lubricant viscosity $\bar{η}$ is considered as the bifurcation parameter, the bifurcation diagram and LLE of the system are shown in Figure 17.

Figure 17.

Bifurcation diagram and LLE of the system with $\bar{η}$ .

It can be seen from Figure 17 that when the lubricant viscosity $\bar{η}$ is in the range of [0.01, 0.09], the system is mainly chaotic, accompanied by a period window, and the corresponding LLE changes between positive value and negative value; when $\bar{η}$ increases to 0.1, the system undergoes the jump shock from chaotic motion to the 4 T-periodic motion. At this moment, the LLE changes from negative value to positive value; when $\bar{η}$ continues to increase to 0.42, the system enters 4 T-quasi-periodic motion. After a short time, the system returns to the 4 T-periodic motion; when $\bar{η}$ continues to increase to 0.45, the system changes from 4 T-periodic motion to 2 T-periodic motion through inverse doubling bifurcation, and the corresponding LLE becomes positive value at this time; when $\bar{η}$ is in the range of [0.45, 1], the system experiences 2 T-periodic motion.

As the previous chapter has described in detail the motion characteristics of the Poincaré map, phase diagram, time history, and FFT spectrum of the system under different motion states, it will not be repeated later. Response of the system: chaotic motion $\to$ 4 T-periodic motion $\to$ 2 T-periodic motion, are shown in Figures 18 to 20.

Figure 18.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when $\bar{η}$ = 0.05.

Figure 19.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when $\bar{η}$ = 0.41.

Figure 20.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when $\bar{η}$ = 0.8.

In order to further explain the influence of the lubricant viscosity $\bar{η}$ on the bifurcation and chaos characteristics of the system, when the lubricant viscosity $\bar{η}$ is set to 0.05 and 0.12, the bifurcation diagram and LLE of the system are obtained, as shown in Figures 21 and 22. Compared with Figure 21, it is analyzed from Figure 22 that at low excitation frequency, the increase of the lubricant viscosity makes the system delay into chaotic motion and reduces the chaotic motion area of the system, and the positive value area corresponding to LLE is also reduced; at high excitation frequency, the increase of the lubricant viscosity complicates the system motion and increases the system’s multiple periodic motion, corresponding to the increase in the negative area of the LLE diagram. In the whole excitation frequency region, as the viscosity of the lubricant increases, the degree of chaos in the system decreases slightly.

Figure 21.

Bifurcation diagram and LLE of the system with $Ω$ when $\bar{η}$ = 0.05.

Figure 22.

Bifurcation diagram and LLE of the system with $Ω$ when $\bar{η}$ = 0.12.

Nonlinear dynamics characteristic of the system with gear backlash

In order to further analyze the effect of backlash on the system in detail, the backlash is refined into the backlash of the internal meshing pair and the backlash of the external meshing pair. The external meshing gear backlash is defined as the sun–planet backlash $b_{s n}$ , and the internal meshing gear backlash is defined as the ring–planet backlash $b_{r n}$ . The dimensionless excitation frequency $Ω$ is set to 1. When the dimensionless sun–planet backlash ${\bar{b}}_{s n}$ is considered as the bifurcation parameter, the bifurcation diagram and LLE of the system are obtained, as shown in Figure 23.

Figure 23.

Bifurcation diagram and LLE of the system with ${\bar{b}}_{s n}$ .

It can be seen from Figure 23 that when the dimensionless sun–planet backlash ${\bar{b}}_{s n}$ is in the range of [0, 5], the motion state of the system consists of chaotic motion and 2 T-periodic motion. With the increase of ${\bar{b}}_{s n}$ , the periodic region of the system is gradually reduced and the chaotic region is gradually increased, and the region with the positive value in the corresponding LLE is also increased.

When ${\bar{b}}_{s n}$ changes in the range of [0, 0.44], the system is mainly in chaotic motion, accompanied by a narrow periodic window, and the corresponding LLE value changes between positive and negative; when ${\bar{b}}_{s n}$ increases to 0.46, the system changes from chaotic motion to 2 T-periodic motion through jump shock, and the corresponding LLE changes from negative value to positive value; when ${\bar{b}}_{s n}$ increases to 2.5, the system changes from 2 T-periodic motion to chaotic motion through jump shock again; when ${\bar{b}}_{s n}$ continues to increase to 3.4, the system jumps from chaotic motion to 2 T-periodic motion, but the region of periodic motion is slightly reduced compared to the previous periodic region, and the chaotic motion region is slightly increased; when ${\bar{b}}_{s n}$ reaches 4.72, the system is chaotic and the LLE value becomes slightly larger, which indicates that the degree of chaos has increased slightly.

The Poincaré map, phase diagram, time history, and FFT spectrum are used to illustrate the system’s motion state. When ${\bar{b}}_{s n}$ is 0.83 and 3.2, respectively, the corresponding system’s motion states are 2 T-periodic motion and chaotic motion, as shown in Figures 24 and 25.

Figure 24.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when ${\bar{b}}_{s n}$ = 0.83.

Figure 25.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when ${\bar{b}}_{s n}$ = 3.2.

In order to further explain the influence of the dimensionless sun–planet backlash ${\bar{b}}_{s n}$ on the bifurcation and chaos characteristics of the system, when ${\bar{b}}_{s n}$ is set to 0.5 and 1.5, the bifurcation diagram and LLE of the system are obtained, as shown in Figures 26 and 27. Compared with Figure 26, it is analyzed from Figure 27 that at low excitation frequency, the increase of ${\bar{b}}_{s n}$ reduces the period window of the system in the chaotic region and increases the chaotic motion, and the corresponding positive value region in the LLE diagram increases; at high excitation frequency, the increase of ${\bar{b}}_{s n}$ complicates the system’s bifurcation characteristics and make system enter into chaotic motion in advance. In the whole excitation frequency range, the increase of ${\bar{b}}_{s n}$ will increase the degree of system vibration.

Figure 26.

Bifurcation diagram and LLE of the system with $Ω$ when ${\bar{b}}_{s n}$ = 0.5.

Figure 27.

Bifurcation diagram and LLE of the system with $Ω$ when ${\bar{b}}_{s n}$ = 1.5.

The dimensionless excitation frequency $Ω$ is set to 1, and other parameters remain unchanged. When the dimensionless ring–planet backlash ${\bar{b}}_{r n}$ is considered as the bifurcation parameter, the bifurcation diagram and LLE of the system are obtained, as shown in Figure 28.

Figure 28.

Bifurcation diagram and LLE of the system with ${\bar{b}}_{r n}$ .

It can be seen from Figure 28 that when the dimensionless ring–planet backlash ${\bar{b}}_{r n}$ changes in the range of [0, 5], the bifurcation and chaotic motion of the system are complicated, and the state of motion changes multiple times through jump shock and doubling bifurcation. When ${\bar{b}}_{r n} \leq$ 2.66, the chaotic motion and periodic motion of the system coexist, and the corresponding LLE changes between positive and negative; when ${\bar{b}}_{r n}$ >2.66, the system is chaotic motion, the corresponding LLE value is positive. At the same time, the LLE gradually increases, indicating that the degree of chaos of the system is enhanced.

When ${\bar{b}}_{r n}$ is in the range of [0, 1.32], the system changes between chaotic motion and 2 T-periodic motion; when ${\bar{b}}_{r n}$ increases to 1.6, the system moves from 2 T-perodic motion to 4 T-perodic motion by doubling bifurcation; when ${\bar{b}}_{r n}$ increases to 1.76, the system changes to chaotic motion by jump shock; when ${\bar{b}}_{r n}$ reaches 1.86, the system experiences 4 T-periodic motion again; when ${\bar{b}}_{r n}$ is in the range of [2, 5], the system mainly experiences chaotic motion, accompanied by the narrow periodic window.

The motion state of the system: 2 T-periodic motion $\to$ 4 T-periodic motion $\to$ chaotic motion, is shown in Figures 29 to 31.

Figure 29.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when ${\bar{b}}_{r n}$ = 0.86.

Figure 30.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when ${\bar{b}}_{r n}$ = 1.64.

Figure 31.

Poincaré map, phase diagram, time history, and FFT spectrum of the system when ${\bar{b}}_{r n}$ = 2.14.

In order to further explain the influence of the dimensionless ring–planet backlash ${\bar{b}}_{r n}$ on the bifurcation and chaos characteristics of the system, when ${\bar{b}}_{r n}$ is set to 0.5 and 1, the bifurcation diagram and LLE of the system are obtained, as shown in Figures 32 and 33. Compared with Figure 32, it is analyzed from Figure 33 that at low excitation frequency, the increase of ${\bar{b}}_{r n}$ causes the system to delay entering chaotic motion and reduce the period window in the chaotic region, which results in two chaotic regions becoming one chaotic region; at high excitation frequency, the increase of ${\bar{b}}_{r n}$ causes the system to move from 1 T-periodic motion into 4 T-periodic motion in advance and have the periodic window in the chaotic region. In the whole excitation frequency region, the increase of ${\bar{b}}_{r n}$ reduces the chaotic region of the system from three chaotic regions to two chaotic regions.

Figure 32.

Bifurcation diagram and LLE of the system with $Ω$ when ${\bar{b}}_{r n}$ = 0.5.

Figure 33.

Bifurcation diagram and LLE of the system with $Ω$ when ${\bar{b}}_{r n}$ = 1.5.

Conclusion

For the planetary gear system considering the coupling effects of friction and lubrication, based on the EHL theory, the stiffness and damping of the squeeze oil film are obtained. They are coupled in series with the meshing stiffness and meshing damping of the gear pair to obtain the system’s equivalent stiffness and equivalent damping. A torsional dynamic model of the planetary gear transmission system considering friction, oil film, time-varying meshing stiffness, meshing damping, and backlash is established. The Runge–Kutta numerical method is used to solve the system vibration equation. The bifurcation diagram, LLE, Poincaré map, phase diagram, time history, and FFT spectrum are used to analyze the bifurcation and chaos characteristics of the system. The conclusions obtained are as follows:

With the increase of the excitation frequency, the planetary gear system exhibits rich bifurcation characteristics, during which multiple bifurcation points and jumping points appear. The system experiences motion states such as short-period motion, long-period motion, and chaotic motion, which shows that the coupling of nonlinear factors such as equivalent stiffness, equivalent damping, and mixed EHL friction makes the planetary gear system exhibit rich nonlinear dynamic behavior.

With the increase of lubricant viscosity, the system begins to experience chaotic motion accompanied by the periodic window, and then jumps from chaotic motion to 4 T-periodic motion through jump shock. After the system temporarily appeared 4 T-quasiperiodic motion, it returned to 4 T-periodic motion. Finally, the system bifurcates from 4 T-periodic motion to 2 T-periodic motion through inverse doubling bifurcation. By changing the lubricant viscosity to analyze the response of the system with excitation frequency, it is concluded that with the increase of the lubricant viscosity, at low excitation frequency, the system delays into chaotic motion, and the range of chaotic motion decreases; at high excitation frequency, the system motion becomes more complex and multi-period motion appear, which results in an increase in the periodic motion area of the system. The chaotic degree of the system decreases slightly throughout the excitation frequency range.

With the increase of sun–planet backlash, the system first changes from chaotic motion to two-periodic motion, and after several changes between chaotic motion and 2 T-periodic motion, the system finally returns to chaotic motion. By changing the sun–planet backlash to analyze the response of the system with excitation frequency, it is concluded that with the increase of the sun–planet backlash, at low excitation frequency, the periodic window in the chaotic region of the system decreases and chaotic motion of the system increases; at high excitation frequency, the bifurcation characteristics of the system become complicated and system enters chaotic motion in advance. In the whole excitation frequency region, the degree of system vibration becomes larger, and the degree of chaos is basically unchanged.

With the increase of the ring–planet backlash, the system first changes between chaotic motion and 2 T-periodic motion, then changes from 2 T-periodic motion to 4 T-periodic motion by doubling bifurcation. Finally system enters into chaotic motion after several jump shocks. By changing the ring–planet backlash to analyze the response of the system with the excitation frequency, it is concluded that with the increase of the ring–planet backlash, at low excitation frequency, the system delays entering chaotic motion and the period window in the chaotic region becomes smaller, which leads to the reduction of two chaotic regions to one chaotic region in the system; at high excitation frequency, the system moves from 1 T-periodic motion to 4 T-periodic motion in advance, and a periodic window appears in the chaotic region.

Highlights

1. Under the coupling effect of friction and mixed elastohydrodynamic lubrication (EHL), a torsional dynamic model of the planetary gear system with friction, oil film, time-varying meshing stiffness, meshing damping, and gear backlash is established.

2. The dynamic characteristics of the planetary gear system are studied from the state of motion, and bifurcation and chaos characteristics, and mainly analyzed from different excitation parameters such as excitation frequency, lubricant viscosity, sun–planet backlash, and ring–planet backlash.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Natural Science Foundation of Liaoning Province of China (2020-MS-216), Liaoning BaiQianWan Talents Program (2019A31), China postdoctoral Science Foundation (2017M610496), State Key Laboratory of Mechanical Transmissions (SKLMT-KFKT-201605), and National Natural Science Foundation of China (5187096).

ORCID iD

Jingyue Wang

References

Kahraman

Singh

Interactions between time-varying mesh stiffness and clearance non-linearities in a geared system. J Sound Vib 1991; 146: 135–156.

Chen

Tang

JY.

Effect of backlash on dynamics of spur gear pair system with friction and time-varying stiffness. J Mech Eng 2009; 45: 161–167.

Liu

Jiang

Liu

, et al. Dynamic behavior analysis of spur gears with constant & variable excitations considering sliding friction influence. J Mech Sci Technol 2016; 30: 5363–5370.

Fang

Liang

Zuo

MJ.

Effects of friction and stochastic load on transient characteristics of a spur gear pair. Nonlinear Dyn 2018; 93: 599–609.

Han

HJ.

Dynamic response analysis of helical gear pair considering the interaction between friction and mesh stiffness. Meccanica 2019; 54: 2325–2337.

Wang

Zhang

Yao

, et al. Nonlinear characteristics of a multi-degree-of-freedom spur gear system with bending-torsional coupling vibration. Mech Syst Signal Process 2019; 121: 810–827.

Qin

DT.

Vibration characteristics of flexible spur ring gears using different connection types. J Low Freq Noise Vib Active Control 2019; 39: 246–257.

Zhu

RP.

Global analysis of a planetary gear train. Shock Vib 2014; 2014: 1–6.

Gao

Zhang

YD.

Nonlinear dynamics analysis of convergent transmission gear containing friction. J Vib Meas Diagn 2014; 34: 737–743.

10.

Sheng

Zhu

Jin

, et al. Bifurcation and chaos study on transverse-torsional coupled 2K-H planetary gear train with multiple clearances. J Cent South Univ 2016; 23: 86–101.

11.

Xiang

Zhang

Gao

, et al. Nonlinear dynamics of a multistage gear transmission system with multi-clearance. Int J Bifurcat Chaos 2018; 28: 1–12.

12.

Hou

Cao

Nonlinear dynamic analysis on planetary gears-rotor system in geared turbofan engines. Int J Bifurcat Chaos 2019; 29: 1–20.

13.

Zhu

Cheng

, et al. Impacts of backlash on nonlinear dynamic characteristic of encased differential planetary gear train. Shock Vib 2019; 2019: 1–15.

14.

Liu

Zhu

, et al. Effects of lubrication on gear performance: a review. Mech Mach Theory 2020; 145: 1–14.

15.

Ramadani

Kegl

Predan

, et al. Influence of cellular lattice body structure on gear vibration induced by meshing. J Mech Eng 2018; 64: 611–620.

16.

Ramadani

Belšak

Kegl

, et al. Topology optimization based design of lightweight and low vibration gear bodies. Int J Simul Model 2018; 17: 92–104.

17.

Houser

Vaishya

Sorenson

JD.

Vibro-acoustic effects of friction in gears: an experimental investigation. SAE Technical Paper, 2001-01-1516, 2001.

18.

Brancati

Rocca

Russo

A gear rattle model accounting for oil squeeze between the meshing gear teeth. P I Mech Eng D-J Aut 2005; 219: 1075–1083.

19.

Theodossiades

Tangasawi

Rahnejat

Gear teeth impacts in hydrodynamic conjunctions promoting idle gear rattle. J Sound Vib 2007; 303: 632–658.

20.

Zhang

Liu

Zhu

, et al. Oil film stiffness and damping in an elastohydrodynamic lubrication line contact-vibration. J Mech Sci Technol 2016; 30: 3031–3039.

21.

Zhu

Wang

, et al. Study on nonlinear dynamic model of planetary gear train sets with friction force. J Vib Shock 2010; 29: 217–220.

22.

Parker

RG.

A physical explanation for the effectiveness of planet phasing to suppress planetary gear vibration. J Sound Vib 2000; 236: 561–573.

23.

Wang

JJ.

System of gear vibration theory. Beijing: Science Press, 1997.

24.

Yuan

Dong

, et al. Dynamic loading analysis of involute gears considering lubrication performance. J Mech Eng 2012; 48: 10–16.

25.

Lin

ZG.

Lubrication of gear transmission. Beijing: Machinery Industry Press, 1980.

26.

Zhu

Analysis of gear bearing capacity. Beijing: Higher Education Press, 1992.

27.

Zhou

Xiao

ZL.

Stiffness and damping models for the oil film in line contact elastohydrodynamic lubrication and applications in the gear drive. Appl Math Model 2018; 61: 634–649.

28.

Wen

SZ.

Elastohydrodynamic lubrication. Beijing: Tsinghua University Press, 1992.

29.

Martin

KF.

A review of friction prediction in gear teeth. Wear 1978; 49: 201–238.

30.

Gunda

Singh

Effect of sliding friction on the dynamics of spur gear pair with realistic time-varying stiffness. J Sound Vib 2007; 301: 927–949.

31.

Cho

Singh

Prediction of dynamic friction forces in spur gears using alternate sliding friction formulations. J Sound Vib 2008; 309: 843–851.

32.

Zhu

YZ.

A computer program package for the prediction of EHL and mixed lubrication characteristics friction subsurface stresses and flash temperatures based on measured 3-D surface roughness. J

Tribol T 2001; 44: 383–390.

33.

Zhu

A full numerical solution to the mixed lubrication in point contacts. ASME J Tribol 2000; 122: 1–9.

34.

Dowson

Higginson

GR.

Elastohydrodynamic lubrication: the fundamental of roller and gear lubrication. Oxford: Pergamon Press, 1966.

35.

Castro

Seabrea

Coefficient of friction in mixed film lubrication: gear versus twin-disc. P I Mech Eng J-J Eng 2007; 22: 399–411.

36.

Zhang

Liu

Jiao

, et al. Vibration analysis of planetary gear trains based on a discrete-continuum dynamic model. J Mech Eng 2014; 50: 104–112.

37.

Cheon

GJ.

Analysis of nonlinear behavior of gear pairs considering hydrodynamic lubrication and sliding friction. J Mech Sci Tech 2009; 23: 2125–2137.

38.

Zou

Chang

GD.

Coupling of dynamics and elastohydrodynamic lubrication for spur gear. J Aerosp Power 2016; 31: 2010–2020.

Nonlinear dynamic characteristics of planetary gear transmission system considering squeeze oil film

Abstract

Keywords

Introduction

The dynamic model of the system

Dynamic equations of the system

Equivalent stiffness of meshing gear pairs considering squeeze oil film

Equivalent damping of meshing gear pairs considering squeeze oil film

The friction force based on mixed EHL

Friction arm of internal and external meshing gear pairs

Dimensionless equations of the system

Numerical simulation and result analysis

Nonlinear dynamics characteristic of the system with excitation frequency

Nonlinear dynamics characteristics of the system with lubricant viscosity

Nonlinear dynamics characteristic of the system with gear backlash

Conclusion

Highlights

Footnotes

Declaration of conflicting interests

Funding

ORCID iD

References