Sage Journals: Discover world-class research

Abstract

Graphical Abstract

Planetary gear set contains many kinds of nonlinear factors, such as time-varying meshing stiffness, bearing support, backlash, and comprehensive transmission error, which will deeply affect the gearbox response behaviors and vibration characteristics. During the operation of planetary gear set, there are several gear pairs meshing simultaneously, which will lead to multiple vibration sources and complicate vibration transmission. The revolution of the planet gears makes the meshing force of the gear pairs vary periodically, and the continuous motion of the meshing point also changes the length of vibration transmission path. This study aims to establish a hybrid model that can simulate the vibration signals of gearbox measuring points through fusing the dynamic model of planetary gear set and the transmission path model. Through the simulation analysis, the vibration characteristics of the measuring point signal of the planetary gearbox are revealed, and the crack fault of the planet gear tooth is effectively isolated. Finally, the availability and validity of the presented model are experimentally verified, and the magnitude match is up to 78%. The model can not only provide a theoretical basis for data-driven fault diagnosis methods, but also solve the pain point of diagnostic samples scale and balance by simulation means.

Keywords

Planetary gearbox dynamic modeling transmission path effect fault diagnosis spectrum analysis

Introduction

Compared with fixed-shaft gear transmission, planetary gear set is widely applied in various fields of modern industry because of its many unique merits. Planetary gear set can accomplish large transmission ratio output in a limited narrow space, and achieve the function of transmitting high power and speed with small volume, light weight, smooth transmission, and high efficiency. Therefore, planetary gear set is used in machine tools, electric power, construction cranes, metallurgy industry, especially in the transmission system of dual-use vehicles, ships, helicopters because of its remarkable power dividing characteristics.

In actual scenarios, the planetary gear set usually works under severe conditions of dynamic heavy load and variable operating conditions. In addition, the structure of planetary gear set is complex, including a large number of vital components such as sun gear, planet gear, ring gear, shafts, and bearings, etc., so the failure rate of planetary gear set remains high. As the critical part of power drive and transmission, once the planetary gearbox fails and cannot be eliminated in time, it will cause economic loss in terms of downtime and mission failure, or serious safety accidents. Therefore, the prognostics and health management (PHM) for planetary gearbox is very urgent and has important practical significance.

So far, there are two main fault diagnosis methods for gear transmission¹: data-driven method and physics-based method. The data-driven method purely relies on the analysis of the previous data from gearbox health status to implement diagnosis or prediction. Data-driven method has attracted the attention of scholars and developed rapidly with the rise of artificial intelligence in recent years. Physics-based method is divided into dynamic modeling and signal modeling.² Based on the in-depth analysis of the structural information and physical characteristics of the gear system, dynamic modeling uses the classical mechanical theory to analyze the dynamic characteristics and internal generation mechanism of the system. The gear system modeling based on dynamics can deepen the understanding of gear damage mechanism and reveal the impact of typical damage modes on the dynamic characteristics of the system in essence. It is of great significance for the fault feature extraction, degradation state identification and remaining useful life prediction of gear box. In addition, the dynamic modeling research can provide simulation data support for the subsequent development of fault diagnosis and prediction algorithms, such as providing a large number of training samples for various artificial intelligence algorithms to solve the problem of samples sparse and unbalance commonly faced by fault diagnosis. Therefore, it is necessary to establish a dynamic model which can reflect the fault and vibration mechanism of planetary gearbox system as much as possible.

In the past decades, some research achievements have been made by domestic and foreign researchers on the physical modeling of planetary gearbox system. Kahraman³ developed a nonlinear dynamic model of planetary gear set including manufacturing error, assembly error, planet gear spacing arrangement, tooth separations, and time-varying stiffness to study the influence of vital parameters on dynamic load sharing factor. Lin and Parker^4,5 set up a classical translational-torsional model of planetary gear set based on the follow-up coordinate system, the natural modes and parametric instability for contact ratios and mesh phasing of gear pairs were also investigated. Parker et al.⁶ established the semi-analytical formulation of the helicopter planetary transmission by using the finite element method, studied the dynamic response of the system under variable working conditions, and derived a detailed calculation process of the meshing phase relationships of the planetary gear set.⁷ Yuksel and Kahraman⁸ quantified the effect of tooth surface wear on the dynamic meshing force impact, and studied its influence mechanism on the nonlinear characteristics of the system. Chaari et al.^9–11 considered the influence of the gyroscopic effect and established the energetic Lagrange formulation of the planetary gear set, focusing on the effects of eccentricity, profile and tooth faults (crack and pitting) on the gear ratio and dynamic characteristics of the system. The lumped-parameter model and the FEM of the planetary gear train were set up by Ambarisha and Parker¹² respectively, and the qualitative and quantitative comparison of simulation results were made, and the law of the gear meshing phase under different bifurcation behaviors was also studied. Cheng et al.^13,14 developed a pure torsional dynamics model for a helicopter planetary gear set, and quantitatively analyzed and evaluated the gear damage level by using statistical index and gray relational analysis algorithm. Sondkar and Kahraman¹⁵ proposed a three-dimensional dynamic model of a double helical planetary gear set, and emphatically demonstrated the influence of key design parameters such as right-to-left stagger angle, support conditions on the dynamic response. Yi¹⁶ proposed the nonlinear dynamic model of wind power gearbox drive train, and evaluated the effect of average error, amplitude error and other factors on gear load sharing characteristics. Fan et al.¹⁷ established the mapping relationship between the length and position of the carrier crack and the torsional vibration response of the planetary gear set by using the FEM. Hammami et al.¹⁸ built a lumped-parameter torsional model of the gearbox with power cycle, simulated the response characteristics of the system under static load and variable load. Zhang et al.¹⁹ put forward a translational-torsional dynamic model of a two-stage planetary gear set, and studied the influence of system parameters on the natural mode by comparing with the FEM. Wei et al.²⁰ integrated the lumped-parameter and finite element method, constructed the dynamic model of multistage planetary gear set based on the element equivalence principle for shaft segment, ring gear and carrier, and studied its quasi-static and dynamic characteristics to guide the reliability design and noise control of the gear system. Liu et al.²¹ obtained the time-varying meshing stiffness of the gear cracked state using the FEM, and substituted it as a key parameter into the lumped-parameter model to simulate early crack response. Tatar et al.²² established a six degrees of freedom lumped-parameter dynamic model of a planetary gearbox, and used modal energy analysis to quantify the coupling level of the rotor system. The results show that the gyroscopic effect has a significant influence on the gearbox dominant mode. Mbarek et al.²³ discussed the influence of variable load and time-varying stiffness on the modal strain and energy distribution of planetary gear set through the hammer test. Sun et al.²⁴ proposed a signal model of planetary gear set considering the transfer path, and studied the effectiveness and advancement of cyclostationary analysis for fault diagnosis under non-stationary rotational speed excitation. Chung et al.²⁵ optimized the load sharing characteristics of the gear system using only floating members by simulating the dynamic model of the planetary gearbox in the off-road vehicle, so as to reduce the weight and miniaturize the transmission system. Zhang et al.²⁶ declared the tidal periods phenomenon of the fault meshing position of the sun gear inside the planetary gearbox, revealing the new formation mechanism of the fault vibration response and defining the minimum sample length required for such diagnosis. Han et al.²⁷ coupled the torsional model of planetary gear set with the MEC model of induction motor, and realized the fault isolation of gear system by identifying the fault characteristics of stator current spectra.

It can be seen from the above literature that dynamic modeling is conducive to in-depth understanding of the natural vibration characteristics and typical fault characteristics, and it is necessary to establish a nonlinear dynamic model that could approximate the real gear system as much as possible, including time-varying stiffness, backlash, transmission error, support stiffness, etc. In addition, the signal model of gearbox measuring point should be developed for verification by real industrial test. Therefore, it is necessary to further investigate the vibration transfer effects between the planetary gear components and the measuring point of the housing. Inalpolat²⁷ proposed a mathematical model to simulate the vibration transfer effect from the meshing forces of planet-ring pairs to the measuring point by using a weighted Hanning function. Later, the influence factors of the meshing forces from the planet-sun pairs were added to improve the transfer effect model.²⁸ Feng²⁹ drew on the above achievements and combined them with the planetary gear set signal model to develop a gearbox measuring point model for exploiting the generation mechanism of fault modulation sideband. Liang et al.^30,31 proposed Hamming and Modified Hamming successively to simulate the vibration transfer path, so as to solve the zero-setting defect of Hanning function when the planet gear rotates furthest away from the sensor. Zghal³² proclaimed a new lumped parameter model by using an energetic approach based on Lagrange formulation, in which the planet gear displacements were implicitly expressed in sensor positions to simulate vibration propagation effects. With the efforts of many researchers, some achievements have been made in the field of transfer path modeling, but there are still some defects that are tough to be overcome. Firstly, when Hanning function is used to represent the periodic modulation generated by the rotation of the planet gear, the end point will be set zero, and its continuous action range is only the minor arc close to the sensor. With the increase of the number of planet gears, this assumption will lead to the risk of error accumulation and even failure of crack tooth meshing simulation. Secondly, whether Hanning or Hamming function is applied as the weighting coefficient of vibration response, its essence is only an empirical expression to approximate the periodic variation of transmission path, and there is no accurate mathematical expression of transmission path for the specific gearbox physical object. Therefore, it is necessary to develop a simulation model that is as close to a realistic physical system as possible, reflecting essentially the typical damage mechanism of a planetary gearbox.

At present, a lot of efforts have been made in the dynamic modeling of planetary gear train, but the comprehensive model containing several nonlinear factors is rarely employed in the simulation. In addition, the accuracy and effectiveness of transmission path modeling for gearbox needs to be improved. In this paper, a translational-torsional model of the fault planetary gear set is established on the basis of classical mechanics theory, which contains several nonlinear factors such as time-varying meshing stiffness, variable bearing support, backlash, and comprehensive transmission error, and the dynamic meshing forces between gear pairs in the system are fully obtained. Then, the FEM is used to analyze the frequency response of the housing to acquire the exact mathematical expression of the vibration transmission effect from the meshing excitation to the sensor measuring point. The resultant vibration signal of the measuring point on the planetary gearbox is obtained by the joint simulation of the two models. Finally, the effectiveness and correctness of the hybrid model for vibration response simulation of fault gearbox are verified by bench test. The proposed planetary gearbox model can effectively deepen the understanding of system fault mechanism and response characteristics, provide simulation big data support for the advanced artificial intelligence algorithms in gearbox fault diagnosis to solve the pain point of samples scarcity and imbalance.

Dynamics modeling and simulation analysis of faulty planetary gear train

In this section, a translational-torsional model of the faulty planetary gear set using lumped parameter method is established, which has the following characteristics compared with the previous studies: (1) several nonlinear factors are incorporated to the model, including time-varying meshing stiffness, bearing support stiffness, system damping, gyroscopic, and centripetal forces, comprehensive transmission error and backlash; (2) some key dynamic parameters are evaluated using more accurate models, such as meshing stiffness calculated by the mainstream potential energy method, bearing support stiffness estimated by the time-varying model. The inclusion of these large number of nonlinear factors allows the presented model to simulate the gearbox physical system at a more refined level, and obtain the three-direction vibration response signals of all components in the system.

Dynamic equations

The planetary gear set usually consists of key components such as the sun gear (s), planet gear (p), ring gear (r), carrier (c), and support bearings, etc. While the planetary gear rotates, it also revolves around the sun gear with the carrier. The translational-torsional model takes into account the vibration of each key component in two translational directions and one torsional direction, whose mechanical sketch is shown in Figure 1. The coordinates $x_{i}$ , $y_{i}$ are transverse displacement and $u_{i}$ is torsional displacement, $i = c, r, s, p_{1}, . . ., p_{N}$ . N represents the number of planetary gears, $N = 4$ in this paper. The Oxy are follow-up rotating coordinate system fixed to the carrier with constant angular speed $ω_{c}$ relative to global coordinate system OXY.

Figure 1.

Translational-torsional model of the planetary gear set.

Based on the above mechanical model, the relative displacement between the sun-planet pair can be expressed by the following formula:

δ_{sp n} = (x_{p n} - x_{s}) \sin φ_{sp n} + (y_{p n} - y_{s}) \cos φ_{sp n} + u_{s} + u_{p n}

(1)

Similarly, the relative displacement of the ring-planet pair is as follows:

δ_{rp n} = (x_{r} - x_{p n}) \sin φ_{rp n} + (y_{r} - y_{p n}) \cos φ_{rp n} + u_{r} - u_{p n}

(2)

where $φ_{sp n} = φ_{p n} - α_{s}$ , $φ_{rp n} = φ_{p n} + α_{s}$ , $φ_{p n}$ represents the circumferential position of the $n - th$ planet gear, $α_{s}$ is the pressure angle of sun gear.

The relative displacement of carrier-planet pair is as follows:

\begin{matrix} δ_{c nx} = x_{c} - x_{p n} - u_{c} \sin φ_{n} \\ δ_{c ny} = y_{c} - y_{p n} + u_{c} \cos φ_{n} \\ δ_{c nu} = - \sin φ_{n} (x_{c} - x_{p n}) + \cos φ_{n} (y_{c} - y_{p n}) + u_{c} \end{matrix}

(3)

Therefore, the dynamic equations under the follow-up coordinate system of all components can be established by using Newton’s second law separately as follows:

For the carrier:

{\begin{matrix} m_{c} ({\overset{\cdot\cdot}{x}}_{c} - 2 ω_{c} {\overset{\cdot}{y}}_{c} - ω_{c}^{2} x_{c}) + \sum_{n = 1}^{N} c_{pn} {\overset{\cdot}{δ}}_{c nx} + \sum_{n = 1}^{N} k_{pn} δ_{c nx} + c_{c} {\overset{\cdot}{x}}_{c} + k_{c} x_{c} = 0 \\ m_{c} ({\overset{\cdot\cdot}{y}}_{c} + 2 ω_{c} {\overset{\cdot}{x}}_{c} - ω_{c}^{2} y_{c}) + \sum_{n = 1}^{N} c_{pn} {\overset{\cdot}{δ}}_{cny} + \sum_{n = 1}^{N} k_{pn} δ_{c ny} + c_{c} {\overset{\cdot}{y}}_{c} + k_{c} y_{c} = 0 \\ (I_{c} / r_{c}^{2}) {\overset{\cdot\cdot}{u}}_{c} + \sum_{n = 1}^{N} c_{pn} {\overset{\cdot}{δ}}_{c nu} + \sum_{n = 1}^{N} k_{pn} δ_{c nu} + c_{ct} {\overset{\cdot}{u}}_{c} / r_{c}^{2} + k_{ct} u_{c} / r_{c}^{2} = - T_{c} / r_{c} \end{matrix}

(4)

For the ring gear:

{\begin{matrix} m_{r} ({\overset{\cdot\cdot}{x}}_{r} - 2 ω_{c} {\overset{\cdot}{y}}_{r} - ω_{c}^{2} x_{r}) - \sum_{n = 1}^{N} F_{rp n} \sin φ_{rp n} + c_{r} {\overset{\cdot}{x}}_{r} + k_{r} x_{r} = 0 \\ m_{r} ({\overset{\cdot\cdot}{y}}_{r} + 2 ω_{c} {\overset{\cdot}{x}}_{r} - ω_{c}^{2} y_{r}) + \sum_{n = 1}^{N} F_{rp n} \cos φ_{rp n} + c_{r} {\overset{\cdot}{y}}_{r} + k_{r} y_{r} = 0 \\ (I_{r} / r_{r}^{2}) {\overset{\cdot\cdot}{u}}_{r} + \sum_{n = 1}^{N} F_{rp n} + c_{rt} {\overset{\cdot}{u}}_{r} / r_{r}^{2} + k_{rt} u_{r} / r_{r}^{2} = 0 \end{matrix}

(5)

For the sun gear:

{\begin{matrix} m_{s} ({\overset{\cdot\cdot}{x}}_{s} - 2 ω_{c} {\overset{\cdot}{y}}_{s} - ω_{c}^{2} x_{s}) + \sum_{n = 1}^{N} F_{sp n} \sin φ_{sp n} + c_{s} {\overset{\cdot}{x}}_{s} + k_{s} x_{s} = 0 \\ m_{s} ({\overset{\cdot\cdot}{y}}_{s} + 2 ω_{c} {\overset{\cdot}{x}}_{s} - ω_{c}^{2} y_{s}) + \sum_{n = 1}^{N} F_{sp n} \cos φ_{sp n} + c_{s} {\overset{\cdot}{y}}_{s} + k_{s} y_{s} = 0 \\ (I_{s} / r_{s}^{2}) {\overset{\cdot\cdot}{u}}_{s} + \sum_{n = 1}^{N} F_{sp n} + c_{st} {\overset{\cdot}{u}}_{s} / r_{s}^{2} + k_{st} u_{s} / r_{s}^{2} = \frac{T_{s}}{r_{s}} \end{matrix}

(6)

For the planet gears:

{\begin{matrix} m_{p n} ({\overset{\cdot\cdot}{x}}_{p n} - 2 ω_{c} {\overset{\cdot}{y}}_{p n} - ω_{c}^{2} x_{p n}) - F_{sp n} \sin φ_{sp n} \\ + F_{rp n} \sin φ_{rp n} - c_{p n} {\overset{\cdot}{δ}}_{c nx} - k_{p n} δ_{c nx} = 0 \\ m_{p n} ({\overset{\cdot\cdot}{y}}_{p n} + 2 ω_{c} {\overset{\cdot}{x}}_{p n} - ω_{c}^{2} y_{p n}) - F_{sp n} \cos φ_{sp n} \\ - F_{rp n} \sin φ_{rp n} - c_{p n} {\overset{\cdot}{δ}}_{c ny} - k_{p n} δ_{c ny} = 0 \\ (I_{p n} / r_{p n}^{2}) {\overset{\cdot\cdot}{u}}_{p n} + F_{sp n} - F_{rp n} = 0 \end{matrix}

(7)

where $F_{sp n}$ , $F_{rp n}$ are the meshing force between the sun-planet, ring-planet pair, respectively.

F_{j} = k_{j} δ_{j} + c_{j} {\overset{\cdot}{δ}}_{j}, j = sp n, rp n

(8)

Considering the tooth profile error in the real gearbox system due to the manufacturing and processing of gear teeth, $e_{sp n}$ , $e_{rp n}$ are adopted to characterize the comprehensive transmission error in sun-planet pair and ring-planet pair, respectively, which could be expression as follows⁹:

e_{j} (t) = E_{j} \sin (2 π f_{m} t + ϕ_{j})

(9)

where $E_{j}$ are the error amplitudes, $f_{m}$ is system meshing frequency, and $ϕ_{j}$ are the error initial phase.

Substitute equation (9) into equations (4)–(7) and write it in matrix form as follows:

\begin{matrix} M \overset{\cdot\cdot}{Q} + (C_{b} + C_{m} + 2 ω_{c} G_{g}) \overset{\cdot}{Q} \\ + (K_{b} + K_{m} - ω_{c}^{2} K_{Ω}) Q = F_{T} + F_{e} \end{matrix}

(10)

where the symbols meaning in the above equations can be found in Appendix.

In order to form a lubricating oil film within the meshing tooth profiles, appropriate backlash must be left to avoid gear meshing transmission jamming caused by tooth friction heat expansion. The backlash formula can be given by the following equation (11),^34,35

f (x, b) = {\begin{matrix} x - b, (x > b) \\ 0, | x | \leq b \\ x + b, (x < - b) \end{matrix}

(11)

where b is half of the backlash constant.

Newmark-β is a method that generalizes the linear acceleration and is widely used to solve dynamics differential equations of mechanical system because of unconditional convergence under set parameters and efficient solution. It is necessary to integrate the backlash formula with the dynamic equations into matrix form for using this efficient numerical simulation method.

Take the sun-planet1 as an example, the dynamic transmission error of the gear pair considering backlash along the action line can be expressed as³⁶:

ξ_{sp 1} = V_{sp 1} Q - e_{sp 1} (t)

(12)

where $V_{sp 1}$ is the corresponding coordinate transformation matrix, and:

\begin{matrix} V_{sp 1} = [0 0 0 0 0 0 \sin (φ_{sp 1}) \cos (φ_{sp 1}) 1 - \sin (φ_{sp 1}) \\ - \cos (φ_{sp 1}) 1 0 0 0 0 0 0 0 0 0] \end{matrix}

(13)

Thus, the dynamic meshing force within sun-planet1 pair can be calculated as follows:

F_{sp 1} = k_{sp 1} (γ_{0} ξ_{sp 1} + γ_{1} b) + c_{sp 1} γ_{0} {\overset{\cdot}{ξ}}_{sp 1}

(14)

where $γ_{0}$ , $γ_{1}$ are the backlash function, which could be obtained by the following equation:

γ_{1} = {\begin{matrix} - 1 & ξ_{sp 1} > b \\ 0 & else \\ 1 & ξ_{sp 1} < - b \end{matrix}, γ_{0} = {\begin{matrix} 1 & | ξ_{sp 1} | > b \\ 0 & else \end{matrix}

(15)

Substituting equations (12)–(15) into equations (4)–(7), a 21 degree-of-freedom planetary gear set dynamics matrix equations considering time-varying meshing stiffness, bearing support, comprehensive transmission error, the effect of gear backlash have been constructed.

Time-varying meshing stiffness evaluation of planetary gear train

Time-varying meshing stiffness is a key parameter in the study of the vibration characteristics of gear systems. At present, there are many kinds of evaluation algorithms for mesh stiffness, such as finite element method, RMAX software method, potential energy method, etc. In addition, ISO6336-1-2006 also provides some reference for mesh stiffness calculation. Mohammed³⁷ studied the mesh stiffness calculation method under large crack size. Chen and Shao^38–40 conducted an intensive study on the meshing stiffness of fixed-axle gear set and planetary gear set, and combined the finite element model with potential energy method to establish the non-penetrating crack stiffness models at different levels. Wan⁴¹ proposed an accurate analytical model to calculate time-varying stiffness which takes into account the fillet-foundation deflections and the flexibility between root circle and the base circle. Ma⁴² considered the transition curve between the tooth root circle and the base circle, and concluded that parabolic curves for crack path and limiting line were optimal by comparing the calculation results of different line types. Then it is improved the disadvantage of repeatedly considering the fillet-foundation deformation in the stiffness model.⁴³ Liang³³ proposed the analytical model of external mesh stiffness using a modified cantilever beam, and investigated the phase relationship of each meshing pair of the planetary gear set. Then, a mapping model of mesh stiffness with severe pitting was developed.⁴⁴ In general, the potential energy method is still the mainstream to evaluate the meshing stiffness. Reference to our previous research fruit,^45,46 the improved potential energy method is applied in this paper.

The improved potential energy method is to equivalent the gear tooth into a cantilever beam, and calculate the Hertz energy $U_{h}$ , bending potential energy $U_{b}$ , shear potential energy $U_{s}$ , compression potential energy $U_{a}$ , and fillet-foundation deformation potential energy $U_{f}$ stored in the meshing gear pair according to material and elasticity mechanics theories. The equivalent stiffness of the meshing gear pairs is obtained through synthesizing the stiffness corresponding to potential energy of each part. The equivalent stiffness of uniform crack is calculated as follows^38,45,47:

\begin{matrix} \frac{1}{k_{e}} = \frac{1}{k_{h}} + \frac{1}{k_{b 1, crack}} + \frac{1}{k_{s 1, crack}} + \frac{1}{k_{a 1}} + \frac{1}{k_{f 1}} + \frac{1}{k_{b 2}} \\ + \frac{1}{k_{s 2}} + \frac{1}{k_{a 2}} + \frac{1}{k_{f 2}} \end{matrix}

(16)

where the subscripts 1 and 2 represent gear and pinion, respectively.

k_{h} = \frac{π EL}{4 (1 - ν^{2})}

(17)

\frac{1}{k_{b, crack}} = \int_{0}^{d} \frac{{(x \cos α_{1} - h \sin α_{1})}^{2}}{E I_{x, crack}} d x

(18)

\frac{1}{k_{s, crack}} = \int_{0}^{d} \frac{1.2 \cos^{2} α_{1}}{G A_{x, crack}} d x

(19)

\frac{1}{k_{a}} = \int_{0}^{d} \frac{\sin^{2} α_{1}}{E A_{x}} d x

(20)

\frac{1}{k_{f}} = \frac{\cos^{2} α_{m}}{EL} {L^{*} {(\frac{u_{f}}{S_{f}})}^{2} + M^{*} (\frac{u_{f}}{S_{f}}) + P^{*} (1 + Q^{*} \tan^{2} α_{m})}

(21)

where the meanings and calculation methods of above formulas can be referred to Shen et al. and Mohammed et al.^47,48

As shown in Figure 2, the gear tooth containing non-uniform crack can be cut into a finite number of slices, and the crack depth of each slice is considered as a constant. The bending stiffness and shear stiffness of the gear teeth are calculated with reference to the uniform crack. The overall bending stiffness and shear stiffness are obtained by numerically integrating the tooth width, as shown in equations (22) and (23).

k_{b, total} = \int_{0}^{L} k_{b, slice} (z)

(22)

k_{s, total} = \int_{0}^{L} k_{s, slice} (z)

(23)

Figure 2.

Schematic diagram of stiffness calculation with non-uniform crack.⁴⁶

Therefore, equation (16) is transformed as:

\begin{matrix} \frac{1}{k_{e}} = \frac{1}{k_{h}} + \frac{1}{k_{b 1, total}} + \frac{1}{k_{s 1, total}} + \frac{1}{k_{a 1}} + \frac{1}{k_{f 1}} + \frac{1}{k_{b 2}} \\ + \frac{1}{k_{s 2}} + \frac{1}{k_{a 2}} + \frac{1}{k_{f 2}} \end{matrix}

(24)

Since the modeling of gear meshing stiffness is not the focus, only the case of uniform crack is investigated for simplicity. The structural parameters of a single-stage planetary gear set in a main reducer are shown in Table 1. The meshing stiffness of planet-ring pair with different crack levels when $β = 70^{°}$ (perfect, 0.5 mm crack, 1.0 mm crack, 1.5 mm crack, 2.0 mm crack, 2.5 mm crack) are evaluated using the improved potential method. The results are shown in Figure 3.

Table 1.

Structural parameters of a planetary gear set.

Parameters	Sun gear	Planet gear	Ring gear	Planet carrier
Number of teeth	28	34	96	–
Inertia (kg·mm²)	219	2387	34,657	15,300
Mass (kg)	0.3929	0.4639	2.614	7.0405
Base radius (mm)	26.3	31.9	90.1	–
Pressure angle (°)		α = 20
Module (mm)		m = 2
Teeth width (mm)		L = 18
Young’s modulus (Pa)		E = 2.068 × 10¹¹
Poisson’s ratio		u = 0.3
Density (kg/m³)		μ = 7800

Figure 3.

Mesh stiffness of different crack levels on a planet gear tooth.

In fact, there are simultaneous meshing actions of eight gear pairs during the operation of the planetary gear set. When the sun gear meshes with all four planet gears, there may be phase difference between their meshing forces. Similarly, there may be phase difference between the four ring-planet pairs. In addition, each planet gear is engaged with both the sun gear and ring gear, and the phase relationship between the two meshing forces is also of concern.

Parker and Lin⁷ conducted in-depth research on the meshing phase relationship of planetary gear set, and found that it is mainly related to the position angle of the planet gear and the number of gear tooth.

Without loss of generality, it is assumed that there is a crack q = 1.5 mm at tooth root (ring gear side) of planet gear p₁. According to the structural parameters listed in Table 1, the amplitude and phase results of eight groups meshing stiffness of the planetary gear set are shown in Figure 4 and Table 2 referring to.⁷ $θ_{m}$ is the rotation angular displacement of the planet gear in one mesh period, whose expression can be referred to.³³ $γ_{rs}$ is the relative phase between the nth ring–planet pair and the $n - th$ sun–planet pair. $γ_{sn}$ / $γ_{rn}$ represent the relative phase between the nth sun–planet/ring–planet pair concerning the first sun–planet pair, $n = 1, 2, . . ., N$ .

Figure 4.

Meshing stiffness curves of planetary gear set with a planet gear crack: (a) $k_{sp 1} / k_{sp 2} / k_{sp 3} / k_{sp 4}$ , (b) $k_{rp 1}$ , and (c) $k_{rp 2} / k_{rp 3} / k_{rp 4}$ .

Table 2.

Relative phases of the planetary gear set.

$γ_{s 1}$	$γ_{s 2}$	$γ_{s 3}$	$γ_{s 4}$	$γ_{r 1}$	$γ_{r 2}$	$γ_{r 3}$	$γ_{r 4}$	$γ_{rs}$
0	0	0	0	0	0	0	0	0.5

It can be seen from Figure 4(a) that the meshing stiffness amplitudes of four groups of sun-planet pairs keep the same and the phases remain synchronous as the perfect condition because the crack on the ring gear side is deemed that the sun-planet pairs are perfect. Figure 4(b) shows a significant amplitude decrease due to the presence of tooth crack, as shown by the yellow circle, which is distinguishing from the meshing stiffness of other three planet-ring pairs in Figure 4(c). In addition, the four phases of the planet-ring meshing forces remain the same. The $γ_{rs} = 0.5$ shown by the red arrow indicates that the corresponding internal and external meshing phases for each planet gear are just opposite.

Bearing stiffness evaluation

As a precision mechanical part, rolling bearing has the remarkable advantages of low friction coefficient and power transmission efficiency. It is widely applied in various fields such as construction machinery and national defense equipment. As an important supporting component of rotating machinery, rolling bearing facilitates the smooth rotation at high speed, reduction friction, and bearing heavy load. It plays an important role in mechanical operation and determines the performance and reliability of mechanical equipment. Rolling bearings are widely used in the planetary gear set of gearbox because of the several components with complex motion relationships. In the planetary gear set listed in Table 1, the planet bearing is HK3020, whose parameters are shown in Table 3 below:

Table 3.

Parameters of planetary gear bearing.

Parameters	Value
Diameter of inner race (mm)	$r_{i}$ = 30
Diameter of outer race (mm)	$r_{o}$ = 37
Needle diameter (mm)	D = 2.5
Needle number	Z = 25
Needle length (mm)	l = 13.8
Static rated load (N)	C₀ = 39,500
Reference speed (rpm)	$n_{r}$ = 6300

As an important part of the gearbox, the mechanical properties of the bearing will also affect the dynamic behavior and vibration response characteristics of the system. From the dynamic model of planetary gear set established in the previous section, bearing stiffness is one of the most critical dynamic parameter. Therefore, it is necessary to evaluate the bearing stiffness.

It is assumed that the bearing only bears radial load $F_{r}$ during operation. when the roller passing through the load area in turn, it will produce elastic contact deformation with the inner and outer ring under the action of external force. In this process, the bearing stiffness will vary periodically with the number of rollers in the bearing area. There will be two extreme conditions during the operation of rolling bearing: (1) the action line of radial force just locates in the middle of the two rollers, as shown in Figure 5(a), namely even pressure; (2) the action line of radial force coincides with the roller center, as shown in Figure 5(b), namely odd pressure.

Figure 5.

Variation of the rollers in spatial position: (a) even pressure and (b) odd pressure.

Without loss of generality, the bearing stiffness under odd pressure state can be evaluated and the bearing radial deformation is obtained by the following empirical formula⁴⁹:

δ_{r} = 7.68 \times 1 0^{- 5} \frac{Q_{max}^{0.9}}{l^{0.8} \cos α_{z}}

(25)

where $Q_{max}$ is the nominal maximum load of the roller on the action line of $F_{r}$ , which could be calculated by the approximate formula as follow:

Q_{max} = \frac{5 F_{r}}{Z \cos α_{z}}

(26)

Substituting equation (26) into equation (25), the bearing stiffness under odd pressure is obtained:

k_{bo} = \frac{F_{r}}{δ_{r}}

(27)

The bearing stiffness under even pressure can also be calculated by the above formula. It should be noted that the maximum nominal load on the roller is the decomposition force on both sides of the action line in this case.

Assuming that the bearing stiffness is isotropic, $k_{h} (t)$ , $k_{v} (t)$ represent bearing support stiffness in the horizontal and vertical directions, respectively, which can be described as:

k_{h, v} (t) = k_{stc} + k_{flct} \sin (2 π f_{bp} t + β_{0})

(28)

k_{stc} = \frac{(k_{bo} + k_{be})}{2}

(29)

k_{flct} = \frac{(k_{be} - k_{bo})}{2}

(30)

where $f_{bp}$ is the roller passing frequency, equal to the product of the rotation frequency of the bearing cage and the number of rollers; $β_{0}$ is the initial phase of roller; $k_{stc}$ , $k_{flct}$ are the static value and time-varying fluctuation component of bearing stiffness, respectively.

Numerical simulation and analysis

In this section, the dynamic model of faulty planetary gear set will be numerically simulated to investigate the vibration response characteristics at different crack levels. The main physical parameters have been listed in Table 1, and other dynamic parameters are shown in Table 4. It is assumed that the sun gear and carrier as input and output components are subjected to the applied torques $T_{s} = 10 Nm$ and $T_{c} = 44.3 Nm$ respectively; The input speed of sun gear is 133 r/min; And the amplitude of comprehensive transmission error and backlash constant b are 50 and 30 μm, respectively. The simulation duration is 3.2 s. Three health conditions of the gearbox are selected for numerical simulation: perfect (case 1), 1.5 mm crack (case 2), 2.5 mm crack (case 3) in planet gear tooth. The operating parameters of the planetary gear set under these conditions are shown in Table 5.

Table 4.

Dynamic parameters of planetary gear set.

Parameters	Value
Bearing damping (N·s/m)	$c_{c} = c_{r} = c_{s} = c_{p} = 3000$
Torsional damping (N·s/m)	$c_{ct} = 0.3357; c_{st} = 0.0352; c_{rt} = 39500;$
Torsional stiffness (N/m)	$k_{rt} = 1.0 \times 10^{10}$ ; $k_{st} = 5 \times 10^{5}$ , $k_{ct} = 1.7 \times 10^{6}$
Meshing damping ratio	$ζ = 0.07$

Table 5.

Operation parameters of the planetary gear set.

Parameters	Value
Meshing frequency f_m (Hz)	48.17
Local fault frequency f_p (Hz)	1.417
Output frequency f_c (Hz)	0.5
Bear passing frequency f_bp (Hz)	16.39

Assumed the crack occurs at planet gear p₁, the dynamic transmission errors $δ_{rp 1}$ between the p₁ and ring gear for the three operating conditions are shown in Figure 6 below.

Figure 6.

Dynamic transmission error $δ_{rp 1}$ in three operating conditions: (a) perfect, (b) 1.5 mm crack, and (c) 2.5 mm crack.

As can be seen in Figure 6(a), the response signal in case 1 is stable, and the basic vibration response comes from the periodic variation of the gear system meshing stiffness. When there is a small-scale crack (case 2) on the planet gear and the fault tooth enters the meshing duration, there are a series of sudden spikes in the time-domain response signal, but the spikes are still very weak at this case, as shown by the red arrows in Figure 6(b). And during the time when the faulty tooth is not involved in the engagement, the response amplitude is basically no different from that in the perfect case. With the extension of the crack length, the peak abrupt changes of the vibration response become more and more obvious, and present vivid periodic characteristics, as shown in Figure 6(c). The spikes period is T_p, which corresponds to the local fault frequency f_p of the planet gear in Table 5. The above features are distinct fault symptoms, which imply that there must be a tooth fault on the planet gear p₁.

Figure 7 shows the frequency-domain fragments of dynamic transmission error $δ_{rp 1}$ in case 3. From the full frequency band shown in Figure 7(a), the dominant component of the response signal is still the system meshing frequency $f_{m}$ and its high-order harmonics $n f_{m}$ (n is integer), namely $f_{main}$ , but the amplitude of other frequency components are very weak relative to the meshing frequency and the spectral details need to be observed at a larger scale. The frequency bands near the first two meshing frequencies are extracted as shown in Figure 7(b). There are many asymmetric modulation sidebands spaced by f_p on both sides of the first two meshing frequencies such as $f_{m} \pm n f_{p}$ , $2 f_{m} \pm n f_{p}$ . This phenomenon indicates that the planet gear is most likely to be a fault. In the low frequency band (0–20 Hz) amplified as shown in Figure 7(c), the planet gear fault characteristic f_p and its high order harmonics are clearly visible, while the passing frequency of the support bearing $f_{bp}$ is still difficult to identify after amplification. Based on the analysis of Figure 7(b) and (c), there is sufficient evidence that a certain fault has occurred in the planet gear tooth, while the stiffness change of the normal support bearing has little effect on the vibration response.

Figure 7.

Spectrum segments of $δ_{rp 1}$ in case 3: (a) 0–300 Hz, (b) 40–105 Hz, and (c) 0–20 Hz.

Figure 8 shows the dynamic transmission errors of the four planet-ring pairs in case 3. According to the comparison of the four subgraphs, the peak mutation of the response signal $δ_{rp 1}$ is the most obvious, followed by that $δ_{rp 3}$ associated with p₃ which is in the symmetrical position of p₁, as shown in Figure 8(c). The abrupt changes in $δ_{rp 2}$ , $δ_{rp 4}$ are already extremely weak and difficult to observe. The simulation results show that although only the tooth of p₁ is cracked, the fault response of p₁ will be transmitted to the response signals of other planet gears in the complex meshing process of planetary gear set, especially for the planet gear in the symmetrical position with the cracked gear. The above conclusions are very important for the discrimination of the planetary gear set regarding the positioning and isolation for the faulty planet gear.

Figure 8.

Dynamic transmission errors of four planet-ring pairs in case 3: (a) $δ_{rp 1}$ , (b) $δ_{rp 2}$ , (c) $δ_{rp 3}$ , and (d) $δ_{rp 4}$ .

The dynamic meshing forces of four groups planet-ring pairs in case 3 can be calculated by equation (8), as shown in Figure 9 below:

Figure 9.

The dynamic meshing forces of four groups planet-ring pairs in case 3: (a) $F_{rp 1}$ , (b) $F_{rp 2}$ , (c) $F_{rp 3}$ , and (d) $F_{rp 4}$ .

The time-domain characteristics of the four sets meshing forces presented in Figure 9 are agreement with the dynamic transmission errors in general. The dynamic meshing force $F_{rp 1}$ between the p1 and gear ring is also selected for frequency-domain analysis, and the results are shown in Figure 10. Comparing Figures 7 with 10, it can be seen that although the frequency amplitudes of the two signals are quite different, the frequency components and distribution of each frequency band are relatively consistent. The time-frequency analysis also indicates the presence of a gear tooth fault in the planet gear, which is not repeated here.

Figure 10.

Spectrum segments of $F_{rp 1}$ in case 3: (a) 0–300 Hz, (b) 40–105 Hz, and (c) 0–20 Hz.

Transmission path modeling of planetary gearbox

The fault dynamics model established in the previous chapter can simulate the three-direction responses of each component of the planetary gear set. However, in the engineering practice of prognostics and health management (PHM) for gearbox system, the sensor devices are usually arranged on the surface of the gearbox housing to monitor its vibration level to avoid destroying the structural integrity. Therefore, in order to simulate the vibration signal of the sensor measuring point, it is necessary to establish a mathematical model of the transmission path from the components of planetary gear set to the sensing device. In this section, the transmission paths of the measuring point signal are analyzed at first. Considering the periodic direction change and time variability of excitation source induced by components meshing, the frequency response functions from meshing forces of gear pairs to the measuring point are obtained by the finite element method and fitted to a multi-level Fourier series. Thus, the accurate mathematical expression of the transmission path model is acquired.

Transmission path effect analysis of measuring point signal

During the transmission process of the planetary gear set, the planet gears engage with both the sun gear and ring gear. There are three main transmission paths of the planet-ring pair meshing vibration, as shown in Figure 11. In the path 1, the vibration is transmitted directly from the gear meshing pair through the ring gear to the sensor (the transmission path of the superior arc is ignored due to the long length), as shown in the white dashed line. The path 2 is from the meshing pair through the planet gear body, sun gear, the planet gear closest to the sensor position, and then to the sensor in turn, as shown by the black dashed line. The path 3 is more complex and transmitted from the meshing pair through the planet gear body, planet gear bearing, support shaft, the carrier in turn, and then from the carrier to the sensor through the closest planet gear, as shown by the red dashed line.²⁹

Figure 11.

Vibration transmission paths of the planet-ring pair.

Compared with the path 1, the path 2 and path 3 are longer and involve more parts. Considering the effects of backlash and damping, the vibration attenuation is more severe, and the vibration characteristics are seriously submerged. Therefore, the path 1 is mainly discussed for the vibration transmission path of the planet-ring pair.

There are three main vibration transmission paths from the planet-sun pair to the sensor, as shown in Figure 12. In the path 1, the vibration signal is transmitted from the meshing pair through planet gear, the ring gear, and then to the sensor in turn, as shown by the white dashed line. In the path 2, the vibration signal is transmitted from the meshing pair through the sun gear, the sun gear shaft, support bearing of sun gear shaft, the housing, and sensor in turn, as shown by the black dashed line. In the path 3, the vibration signal is transmitted from the meshing pair through the planet gear, the planet bearing and shaft, the carrier, the carrier bearing, the housing, and the sensor in turn, as shown by the red dashed line. The length of the last two transmission paths does not change with the rotation of the carrier.

Figure 12.

Vibration transmission paths of planet-sun pair.

Likewise, the characteristics of the vibration signal may change or useful characteristic information may be lost due to the long route of the last two transmission paths, especially through the complex decay process of bearings. Therefore, the transmission path 1 will be considered as a priority.

As shown above, the shortest and most efficient transmission path for both planet-ring pairs and planet-sun pairs ends up being through the engagement of the planet-ring pairs. Therefore, this part will be modeled in detail.

Mathematical modeling of transmission path

For the planetary gear set with the inner ring gear fixed on the housing, the transmission path from the meshing point of the planet-ring pair to the sensor arrangement point changes periodically during the revolution of planet gear. During the process of the engagement point close to the sensor position, the transmission path gradually becomes shorter and the vibration signal perceived by the sensor becomes stronger and stronger; while when the engagement point gradually moves away from the sensor position, the transmission path becomes longer and the vibration signal obtained by the sensor continuously becomes weaker. The Hanning window and Hamming window functions are often used to describe the effects of transfer paths,^27–31 but these functions have certain drawbacks. The aforementioned window functions are only a simplified mathematical formulation for the transfer path effect. Firstly, for specific gearbox housing objects, different physical properties such as size, structure, and material will reflect different transfer characteristics, which cannot be accurately described by the aforementioned window function. Secondly, the Hanning window only considers the meshing vibration between the part of the inner ring gear circumference (1/N) near the sensor and the planet gear, and the simulated transmission path error will be great when the number of planet gears is large.²⁷ If the failure occurs within the neglected meshing arc, the Hanning window function will fail the simulation. Therefore, this paper will apply finite element method to simulate the housing in modal test, and fit an accurate transmission path model of the planetary gearbox vibration response by calculating the frequency response function between the ring gear and the sensor.

Figure 13(a) and (b) show the 3D geometric model and finite element model of a planetary gearbox housing, respectively. The hexahedral element with fast convergence and high precision is adopted as the element type. The element side length is set 10 mm, and the housing material is steel. The force hammer experiment is simulated through the software to calculate the frequency response function of the housing. The arrow in Figure 13(b) shows the force hammer loading point, and the ring gear is marked successively in clockwise direction at equal intervals. In engineering practice, the sensors are usually arranged in the Z-direction perpendicular to the horizontal plane. Therefore, the vertical acceleration signals of the nodes 1–7 are extracted, and the amplitude frequency functions of all response point to the excitation position are calculated as shown in Figure 14 below. As can be seen from the figure, the farther the nodes of the ring gear are from the excitation point, the smaller the amplitudes of frequency response functions. Due to the structural symmetry of the housing, the frequency response functions of nodes on the left half ring gear to the excitation point can be obtained by mirror.

Figure 13.

Planetary gearbox model: (a) 3D geometric model and (b) FEM.

Figure 14.

Amplitude frequency function of housing.

According to the simulation condition of the dynamic model in the previous chapter, the system meshing frequency f_m is 48.17 Hz, as shown in Table 5. The extracted frequency response amplitudes of the ring gear circumference under this case are shown in Table 6 below.

Table 6.

Frequency response amplitude of the nodes 1–7 under meshing frequency.

Nodes	1	2	3	4	5	6	7
Amplitude (m/s²/N)	0.0163	0.0152	0.0118	0.0101	0.0093	0.0077	0.0070

By fitting the frequency response amplitude of each discrete response point with Fourier series, the frequency response function from the ring gear circumference to the excitation position can be obtained as follows equation (31):

h (ψ, f_{m}) = a_{0} + \sum_{p = 1}^{P} a_{l} \cos (p ψ)

(31)

where $ψ$ is the angle variable, P is the fitting series, $a_{0}$ , $a_{l}$ are the Fourier expansion coefficient, whose value are related to the housing structure, material, and constraints.

Considering many kinds of error indicators, the frequency response amplitudes in Table 6 are fitted by fourth-order Fourier series, and the results are shown in Figure 15.

Figure 15.

Fourier series fitting curve.

To facilitate the joint simulation with the dynamics model of the planetary gear set, the independent variable $ψ$ in equation (31) is converted to the time domain, $ψ = 2 π f_{c} t + φ_{p n}$ . Then, the amplitude frequency function of the vibration signal meshing between the $n - th$ planet gear and the entire ring circumference can be expressed as follows:

h_{n} (t, f_{m}) = {a_{n}}_{0} + \sum_{p = 1}^{P} a_{nl} \cos (2 p π f_{c} t + p φ_{p n})

(32)

Characteristic analysis of resultant signal

In this section, the dynamic model of planetary gear set and transmission path model are fused for a joint numerical simulation, and finally the virtual measuring point signal of gearbox, namely resultant signal, could be obtained. Through the time domain and frequency domain analysis of resultant signal, the transmission characteristics and dynamic characteristics of the gearbox vibration signal are deepened, and the gear tooth fault mechanism of the planetary gear set is revealed. These analyses are conducive to the development of fault diagnosis algorithms.

The mathematical model of transmission path is established under the global coordinate system OXY, while the dynamic model of planetary gear set is developed based on the follow-up coordinate system Oxy. In order to achieve a joint simulation, the coordinate systems need to be unified. Firstly, the dynamic meshing forces $F_{rpn}$ of the planet-ring pairs are projected to the follow-up coordinate system:

{\begin{matrix} F_{rp nx} = F_{rp n} \sin φ_{rpn} \\ F_{rp ny} = - F_{rp n} \cos φ_{rpn} \end{matrix}

(33)

where $F_{rp nx}$ , $F_{rp ny}$ are the projections of the meshing forces in x, y direction respectively.

According to the transformation relationship between the follow-up coordinate system Oxy and the global coordinate system OXY,⁴ the projections of $F_{rp nx}$ , $F_{rp ny}$ to the global coordinate Y-direction lead to:

F_{rp nY} = F_{rp nx} \sin (2 π f_{c} t) + F_{rp ny} \cos (2 π f_{c} t)

(34)

where $F_{rp nY}$ is the dynamic meshing component force of the $n - th$ planet-ring pair in the Y-direction of the global coordinate system.

Therefore, the dynamic meshing forces $F_{rpn}$ of the planet-ring pairs obtained from the dynamics model simulation of the planetary gear set in case 1 and case 3 are projected in the global coordinate system, as shown in Figures 16 and 17 below.

Figure 16.

$F_{rp nY}$ in case 1: (a) $F_{rp 1 Y}$ , (b) $F_{rp 2 Y}$ , (c) $F_{rp 3 Y}$ , and (d) $F_{rp 4 Y}$ .

Figure 17.

$F_{rp nY}$ in case 3: (a) $F_{rp 1 Y}$ , (b) $F_{rp 2 Y}$ , (c) $F_{rp 3 Y}$ , and (d) $F_{rp 4 Y}$ .

Figure 16 shows that the four groups dynamic meshing forces of planet-ring pairs in case 1 present obvious amplitude-frequency modulation after two projections, and the phase difference between the adjacent projection forces is $π / 2$ . Comparing Figure 17 with Figure 16, there are visible force shocks with period $T_{p}$ for $F_{rp 1 Y}$ and $F_{rp 3 Y}$ in case 3. The shock amplitudes of $F_{rp 1 Y}$ are greater than $F_{rp 3 Y}$ , while the shock amplitudes of $F_{rp 2 Y}$ , $F_{rp 4 Y}$ are too small to distinguish. These time-domain fault characteristics presented are completely consistent with Figure 9.

The mathematical model of the resultant acceleration signal at the sensor measuring point is given as follows, and the time-domain simulation results are shown in Figure 18.

ac c_{sensor} = \sum_{n = 1}^{N} F_{rp nY} \cdot h_{n} (t, f_{m})

(35)

Figure 18.

Time-domain resultant accelerations: (a) case 1 and (b) case 3.

Observing the time-domain resultant accelerations under two cases in Figure 18, the obvious features are that both signals reflect the high-frequency vibration caused by time-varying meshing stiffness and the low-frequency modulation due to carrier rotation. Comparing the two subgraphs, it can be seen that there are weak periodic impacts on the acceleration signal due to tooth fault in case 3, as shown by the red arrow in Figure 18(b), which requires detailed observation to be known. However, at the time of fault tooth meshing out, the acceleration signals are basically the same in the two cases. Therefore, gear fault diagnosis from the perspective of time domain is not easy, and the signal could be converted to frequency domain for analysis.

Figure 19 shows the frequency domain comparison of acceleration signals under the two conditions. The main frequency components $f_{main}$ of the two signals have basically not changed. The obvious difference lies in the presence of very small modulation bands of bearing pass frequency near the main frequency components in case 1, as shown in Figure 19(a). But these weak bearing pass frequency bands are obliterated by the large bands brought by the faulty gear in case 3. In order to explore the detailed features of the signals in the frequency domain, the frequency band near the second-order meshing frequency is analyzed with different amplification levels as shown in Figure 20.

Figure 19.

Frequency-domain resultant accelerations: (a) case 1 and (b) case 3.

Figure 20.

Resultant acceleration frequency-domain sections: (a) 75–115 Hz, (b) 75–85 Hz, and (c) 92–102 Hz.

Figure 20 shows the frequency band characteristics around the second order meshing frequency (75–115 Hz). Common ground of both cases is that there are several modulation sidebands with strong energy interval of $4 f_{c}$ near $2 f_{m}$ .The main difference is the presence of weak sidebands with spacing $f_{bp}$ on both sides of $2 f_{m}$ in case 1. The frequency band near the left sideband $2 f_{m} - f_{bp}$ is enlarged as shown in Figure 20(b), which shows that this sideband also continues to be modulated. The right sideband is also similar. In case 3, all sidebands generated by the bearing through the frequency have been completely covered. Zooming in on the adjacent frequency band of $2 f_{m}$ , as shown in Figure 20(c), only the sideband peaks of $2 f_{m} \pm 4 f_{c}$ are clear except for extremely weak $2 f_{m} \pm f_{c}$ in case 1, while there are more complex modulation sidebands in case 3. Sizable sidebands appear in the following locations: $2 f_{m} \pm k f_{c} \pm h f_{p}$ (k, h are integers). The above frequency-domain analysis clearly implies that there must be a planet gear tooth fault in case 3, which is also consistent with the previous simulation assumptions and the conclusions of the literatures.^28–30,32

From the above time-frequency domain amplitude analysis, it is clear that the time-varying stiffness variation of the perfect bearing has a very limited effect on the resultant acceleration signals. In addition, $4 f_{c}$ is the most considerable symmetric modulation sideband in both cases, because the symmetric distribution of the planet gear makes the other frequency components almost cancel each other except for $4 f_{c}$ . In case 3, the planet gear symmetry is relatively weakened due to a planet gear tooth fault, so that the modulation bands at intervals of $\pm f_{c}$ , $\pm 2 f_{c}$ , $\pm 3 f_{c}$ are significantly enhanced, compared to case 1 as shown in Figure 20(c).

Comparison and verification

In this chapter, the resultant simulation signals of presented path model and the classical window function path model are compared in the time domain. Then, the consistency of the resultant signals deduced from the two models with the physical test signals are compared by building a fault diagnosis test platform of transmission system in the laboratory environment. Finally, it is verified that the resultant signal of presented model is accurate and effective to simulate the experimental signals of the faulty planetary gear set through the spectral comparative analysis.

Simulation comparison

To illustrate the advancement of transfer path model proposed in this paper, it is compared with the typical window functions developed in,^28–33 as shown in Figure 21. Liang et al.³¹ considered that the measured point signal could be regarded as the summation of window weighted vibration of each planet gear because the vibration of planet gears contains both the information of sun-planet mesh and ring-planet mesh, as shown in the equation (36):

a (t) = \sum_{n = 1}^{N} w_{n} (t) a_{n} (t)

(36)

Figure 21.

Typical window functions for transmission path effect.³²

where $w_{n} (t)$ is the weighting function which counts for the effect of the transmission path; $a_{n} (t)$ represents the acceleration of the n-th planet.

Through the comparative analysis by Liang et al.³² it was concluded that the resultant vibration signal has relatively good consistency with the experimental signal when the modified Hamming window with $α = 0.1$ was applied to simulate the transmission path effect, that is, $w_{n} (t) = e^{α {(mod (ω_{c} t + ψ_{n}, 2 π) - π)}^{2}} H_{n} (t)$ . Therefore, this modified window function is focused as the comparison object for the presented path effect model.

Figure 22 shows the resultant acceleration signals of faulty gearbox measuring points under two different transmission path models. As can be seen from the figure that the amplitude of the simulation acceleration signal weighted by the modified Hamming window is significantly larger than the simulation signal using the presented fitting model.

Figure 22.

Comparison of resultant acceleration using two transfer path models.

Experimental verification

To verify the accuracy and effectiveness of the planetary gearbox hybrid model with a tooth fault, a test platform of gear transmission system for fault diagnosis and prediction was built. The platform consists of drive motor, bevel gear pair, two-stage planetary gear set (main reducer), load motor, lubricating oil pump, control system, and data acquisition system, as shown in Figure 23.

Figure 23.

Gearbox fault diagnosis and prediction test platform: (a) test rig and (b) data acquisition system.

The rated power of the driving motor is 40 hp, and the speed is adjustable in the range of 0–4000 r/min. The rated power of load motor is 100 hp, and the load torque is adjustable in the range of 0–142 Nm. The control system and data acquisition system of the experimental platform are integrated. The control system is used to set the working mode and bench parameters, while the data acquisition system could collect eight vibration signals, one rotating signal, and one torque signal synchronously, and the sampling frequency range is 2.56–51.2 kHz. The accelerometer is model 608A11 of IMI company, whose sensitivity is 100 mV/g and frequency range is 0.5–10 kHz.

The main body of the test platform is a three-stage gear reducer, which has a total reduction ratio of 39.85:1, and other specific parameters are configured in Table 7.

Table 7.

Parameters of gear transmission system.

Gear system	Bevel stage		First planetary stage			Second planetary stage
Gear	Input	Output	Sun	Planet	Ring	Sun	Planet	Ring
Teeth number	18	36	32	40 (3)	112	28	34 (4)	96
Transmission ratio	2		4.5			4.4286

The test platform is capable of simulating multiple fault modes and damage evolution processes of key components in the transmission system, providing a validation environment for damage dynamics analysis, degradation state identification and remaining useful life prediction.

To eliminate the influence of environment noise and the original errors caused by installation and manufacturing, all gears are considered perfect before manual fault implantation, while the test under this condition was defined as a benchmark test. Sensors were arranged on the side and perpendicular to the horizontal plane of the main reducer housing, as shown in Figure 24. Spare parts for the second-stage planetary gear set p₁ was processed with the wire cutting for root cracks of virous sizes to simulate different crack levels, as shown in Figure 25. The input speed and load torque of the transmission system were controlled and adjusted by the drive motor and the load motor respectively, which could make the test settings of the second-stage planetary gear set consistent with the settings of the simulation conditions.

Figure 24.

Sensor layout in main reducer.

Figure 25.

Multiscale crack implanted in planet gear.

Several sets of experiments were conducted under different operating conditions to analyze the effects of planet gear cracks on the vibration characteristics of gearbox. Figure 26 shows the time-domain vertical acceleration signals of the benchmark test and fault test with q = 2.5 mm.

Figure 26.

Time-domain experimental signals: (a) benchmark test and (b) fault test.

Visually comparing the simulated signals shown in Figure 22 with the experimental signals shown in Figure 26(b), the overall waveform and amplitude of the resultant signal using the proposed path model relatively match the experimental signal about 78%, while the amplitude of resultant signal weighted the modified Hamming window obviously has a large deviation. This phenomenon occurs because the modified Hamming window is only an extremely simplified empirical representation model class for path effect without rigorous mechanical analysis and the value $α$ in equation (40) also is a rough preset value, which actually appears in Liang et al.³² as well.

Actually, the amplitude of resultant signal using present path model is also smaller than that of the experimental signal evidently, because the resultant signal ignores multiple other minor vibration transmission paths. In addition, the experimental signal will inevitably be mixed with a variety of noises effect, which is also an important factor in the amplitude difference. It is worth noting that there is no significant difference in the time domain curves between the two experimental conditions, except that the overall signal amplitude in the crack condition is slightly larger than the benchmark test. The periodic impact symptoms are not clearly observed, which are overwhelmed by the noise mixed in the measuring process. Therefore, the frequency domain view may be a better perspective to analyze the fault symptoms.

Figure 27 depicts the comparison of the spectrum segments of the test signals in the two cases, covering the frequency band near the first four order meshing frequencies. As can be seen from the figure that the spectral structure and components of the signals in the two cases are basically the identical, and the difference lies in the magnitudes of the frequency components. In the vicinity of the first four order meshing frequencies, except for $f_{main}$ with significant amplitudes, the sidebands $f_{main} \pm 4 f_{c}$ are no longer dominant compared to other frequency components such as $f_{main} \pm f_{c}$ , $f_{main} \pm 2 f_{c}$ , $f_{main} \pm 3 f_{c}$ , even the amplitude is significantly smaller than the others, which is different from the simulation cases shown in Figure 20 especially for case 1. That is because there are inevitable system mounting and gear manufacturing errors for a real physical gearbox, which could significantly weaken the symmetry of the planetary gear set and reduce the phases extinction effect due to the homogeneous distribution of the four planet gears. Overall, sizable sideband amplitudes emerge in the following locations in two cases: $f_{main} \pm k f_{c} \pm h f_{p}$ . However, the sideband amplitudes with tooth crack are almost greatly amplified compared to the benchmark test, except for the sideband $3 f_{m} - 2 f_{c}$ . The occurrence of fault characteristic sidebands in the benchmark test is acceptable in practice, because there may be certain invisible minor faults in the benchmark condition due to material defects or manufacturing errors, which will introduce sidebands into the test frequency spectrum. Then, the remarkable growth of the fault sidebands in the crack case implies that the planet gear in the gearbox system must occur an apparent fault,^27–29 which is also consistent with the analysis conclusions of the resultant simulation signal.

Figure 27.

Spectrum segments of two conditions: (a) 46–50 Hz, (b) 94–98 Hz, (c) 142.5–146.5 Hz, and (d) 191–195 Hz.

In order to implement a more intuitive comparison, a quantitative analysis of the sidebands ( $f_{main} \pm f_{p}$ ) on both sides of the first four meshing frequencies are performed and the results are shown in Figure 28 below.

Figure 28.

Comparison of fault modulation sidebands near the first four order meshing frequencies.

As can be seen from the above figure, there is a significant increase in the amplitudes of the sidebands $f_{main} \pm f_{p}$ near the first four order meshing frequencies in the crack case relative to the benchmark case, ranging from 144% to 6015%, which further confirms the inference that the planet gear tooth have failed.

In addition, the low frequency band of two signals are amplified in detail as shown in Figure 29 below. The response amplitudes are extremely weak for both the benchmark case and crack case, including the fault characteristic frequency $f_{p}$ and the bearing passing frequency $f_{bp}$ . The phenomenon gives a hint that the analysis of the fault modulation sidebands may be a more effective way for fault diagnosis. Besides, the stiffness fluctuation of healthy bearing in the gearbox is exceedingly faint or even negligible for the vibration response in fault diagnosis engineering applications.

Figure 29.

The low-frequency segments: (a) 0–5 Hz and (b) 12–18 Hz.

In summary, the vibration responses of different crack levels are analyzed in the time and frequency domains, and the accuracy and validity of the hybrid model of the sensor measuring point signal on the gearbox are verified. The test signals of the faulty transmission system could be simulated and restored by the proposed model to a large extent.

Conclusion

In this paper, a 21 degree-of-freedom translational-torsional model of the faulty planetary gear set is first developed to simulate the three-direction vibration responses of each component, incorporating several nonlinear factors such as time-varying mesh stiffness and damping, time-varying bearing support, backlash, and comprehensive transmission error. Finite element technology is introduced into virtual frequency response analysis to obtain an accurate transfer path model from the gear components to gearbox sensor measuring point considering the influence of simultaneous multi-point meshing vibration sources, which is different from the existing empirical representation window function. The resultant signal obtained by fusing the two aforementioned models can simulate the gearbox measuring point signal about 78%. In the time domain, the resultant signal has significant modulation phenomenon of carrier rotation frequency, but the fault symptoms and their periodicity are attenuate and not easily detectable; in the frequency domain, a large number of modulation sidebands are present when the planet gear occurs a crack, and the identification and localization of these sidebands can serve better for fault diagnosis. Several benchmark tests and fault implantation tests are conducted to confirm the accuracy and validity of the hybrid model for fault diagnosis. In addition, both simulation and experimental data support the interesting conclusion that the time-varying nature of healthy bearing stiffness has a very limited effect on the response characteristics of the gearbox system. The measuring point model of gearbox established in this paper can not only effectively reveal the fault mechanism and behavior characteristics, but also profoundly reveal the conduction rule of fault response inside the gearbox system, and provide front-end simulation big data support for promoting fault diagnosis and prediction using artificial intelligence algorithms.

Footnotes

Appendix

Handling Editor: Chenhui Liang

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: This work is financial supported by the National Key Research and Development Program of China (Grant No. 2018YFB1702401), the National Natural Science Foundation of China (Grant No. 51975576, 52105133), and the Defense Industrial Technology Development Program (Grant No. WDZC20205500301, WDZC20205250305).

ORCID iD

Jiao Hu

References

Liang

Zuo

Feng

Dynamic modeling of gearbox faults: a review. Mech Syst Signal Process 2018; 98: 852–876.

Parra

Vicuña

CM.

Two methods for modeling vibrations of planetary gearboxes including faults: comparison and validation. Mech Syst Signal Process 2017; 92: 213–225.

Kahraman

Load sharing characteristics of planetary transmissions. Mech Mach Theory 1994; 29: 1151–1165.

Lin

Parker

RG.

Analytical characterization of the unique properties of planetary gear free vibration. J Vib Acoust 1999; 121: 316–321.

Lin

Parker

RG.

Planetary gear parametric instability caused by mesh stiffness variation. J Sound Vib 2002; 249: 129–145.

Parker

Agashe

Vijayakar

SM.

Dynamic response of a planetary gear system using a finite element/contact mechanics model. Asme J Mech Des 2000; 122: 304–310.

Parker

Lin

Mesh phasing relationships in planetary and epicyclic gears. J Mech Des 2004; 126: 365–370.

Yuksel

Kahraman

Dynamic tooth loads of planetary gear sets having tooth profile wear. Mech Mach Theory 2004; 39: 695–715.

Chaari

Fakhfakh

Hbaieb

, et al. Influence of manufacturing errors on the dynamic behavior of planetary gears. Int J Adv Manuf Technol 2006; 27: 738–746.

10.

Chaari

Fakhfakh

Haddar

Dynamic analysis of a planetary gear failure caused by tooth pitting and cracking. J Fail Anal Prev 2006; 6: 73–78.

11.

Chaari

Fakhfakh

Haddar

Analytical investigation on the effect of gear teeth faults on the dynamic response of a planetary gear set. Noise Vib Worldw 2006; 37: 9–17.

12.

Ambarisha

Parker

RG.

Nonlinear dynamics of planetary gears using analytical and finite element models. J Sound Vib 2007; 302: 577–595.

13.

Cheng

Quantitative damage detection for planetary gear sets based on physical models. Chin J Mech Eng 2012; 25: 190–196.

14.

Cheng

Zhang

Crack level estimation approach for planetary gearbox based on simulation signal and GRA. J Sound Vib 2012; 331: 5853–5863.

15.

Sondkar

Kahraman

A dynamic model of a double-helical planetary gear set. Mech Mach Theory 2013; 70: 157–174.

16.

Zhang

Guo

, et al. Dynamic modeling and analysis of load sharing characteristics of wind turbine gearbox. Adv Mech Eng 2015; 7: 177–190.

17.

Fan

Wang

, et al. Nonlinear dynamic modeling of a helicopter planetary gear train for carrier plate crack fault diagnosis. Chin J Aeronaut 2016; 29: 675–687.

18.

Hammami

Fernandez Del Rincon

Chaari

, et al. Effects of variable loading conditions on the dynamic behaviour of planetary gear with power recirculation. Measurement 2016; 94: 306–315.

19.

Zhang

Wang

, et al. Dynamic modeling and vibration characteristics of a two-stage closed-form planetary gear train. Mech Mach Theory 2016; 97: 12–28.

20.

Wei

Zhang

Qin

, et al. A coupling dynamics analysis method for a multistage planetary gear system. Mech Mach Theory 2017; 110: 27–49.

21.

Liu

Yang

Zhang

Resultant vibration signal model based fault diagnosis of a single stage planetary gear train with an incipient tooth crack on the sun gear. Renew Energy 2018; 122: 65–79.

22.

Tatar

Schwingshackl

Friswell

MI.

Dynamic behaviour of three-dimensional planetary geared rotor systems. Mech Mach Theory 2019; 134: 39–56.

23.

Mbarek

Hammami

Fernandez

, et al. Influence of the non-linear hertzian stiffness on the dynamic behavior of planetary gear during run up condition. Adv Cond Monit Mach Non-Stat Oper 2019; 30: 30–38.

24.

Sun

Yang

Gryllias

, et al. Cyclostationary modeling for local fault diagnosis of planetary gear vibration signals. J Sound Vib 2020; 471: 165–175.

25.

Chung

Han

, et al. Analytical study of floating effects on load sharing characteristics of planetary gearbox for off-road vehicle. Adv Mech Eng 2020; 12: 168–178.

26.

Zhang

Zuo

Wei

, et al. Motion periods of sun gear dynamic fault meshing positions in planetary gear systems. Measurement 2020; 162: 107897.

27.

Han

Wang

Ding

, et al. Magnetic equivalent modeling of stator currents for localized fault detection of planetary gearboxes coupled to electric motors. IEEE Trans Ind Electron 2021; 68: 2575–2586.

28.

Inalpolat

Kahraman

A theoretical and experimental investigation of modulation sidebands of planetary gear sets. J Sound Vib 2009; 323: 677–696.

29.

Inalpolat

Kahraman

A dynamic model to predict modulation sidebands of a planetary gear set having manufacturing errors. J Sound Vib 2010; 329: 371–393.

30.

Feng

Zuo

MJ.

Vibration signal models for fault diagnosis of planetary gearboxes. J Sound Vib 2012; 331: 4919–4939.

31.

Liang

Zuo

Hoseini

MR.

Understanding vibration properties of a planetary gear set for fault detection. IEEE Conf Progn Health Manag 2014; 30: 1–6.

32.

Liang

Zuo

Hoseini

MR.

Vibration signal modeling of a planetary gear set for tooth crack detection. Eng Fail Anal 2015; 48: 185–200.

33.

Zghal

Graja

Dziedziech

, et al. A new modeling of planetary gear set to predict modulation phenomenon. Mech Syst Signal Process 2019; 127: 234–261.

34.

Liang

Zuo

Pandey

Analytically evaluating the influence of crack on the mesh stiffness of a planetary gear set. Mech Mach Theory 2014; 76: 20–38.

35.

Kahraman

Singh

Non-linear dynamics of a geared rotor-bearing system with multiple clearances. J Sound Vib 1991; 144: 469–506.

36.

Kahraman

Singh

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

37.

Yang

Tang

, et al. Dynamic analysis for a spur geared rotor system with tooth tip chipping based on an improved time-varying mesh stiffness model. Mech Mach Theory 2021; 165: 104435.

38.

Mohammed

Rantatalo

Aidanpää

JO.

Improving mesh stiffness calculation of cracked gears for the purpose of vibration-based fault analysis. Eng Fail Anal 2013; 34: 235–251.

39.

Chen

Shao

Dynamic simulation of spur gear with tooth root crack propagating along tooth width and crack depth. Eng Fail Anal 2011; 18: 2149–2164.

40.

Chen

Shao

Mesh stiffness calculation of a spur gear pair with tooth profile modification and tooth root crack. Mech Mach Theory 2013; 62: 63–74.

41.

Chen

Zhai

Shao

, et al. Analytical model for mesh stiffness calculation of spur gear pair with non-uniformly distributed tooth root crack. Eng Fail Anal 2016; 66: 502–514.

42.

Wan

Cao

, et al. An improved time-varying mesh stiffness algorithm and dynamic modeling of gear-rotor system with tooth root crack. Eng Fail Anal 2014; 42: 157–177.

43.

Song

Pang

, et al. Time-varying mesh stiffness calculation of cracked spur gears. Eng Fail Anal 2014; 44: 179–194.

44.

Pang

Feng

, et al. Improved time-varying mesh stiffness model of cracked spur gears. Eng Fail Anal 2015; 55: 271–287.

45.

Liang

Zhang

Liu

, et al. The influence of tooth pitting on the mesh stiffness of a pair of external spur gears. Mech Mach Theory 2016; 106: 1–15.

46.

Yang

Xia

Han

, et al. Vibration analysis for tooth crack detection in a spur gear system with clearance nonlinearity. Int J Mech Sci 2019; 157–158: 648–661.

47.

Shen

Zhang

, et al. Dynamic analysis of planetary gear with Root Crack in Sun gear based on improved time-varying mesh stiffness. Appl Sci 2020; 10: 8379.

48.

Mohammed

Rantatalo

Aidanpää

, et al. Vibration signal analysis for gear fault diagnosis with various crack progression scenarios. Mech Syst Signal Process 2013; 41: 176–195.

49.

Sainsot

Velex

Duverger

. Contribution of gear body to tooth deflections—a new bidimensional analytical formula. J Mech Des 2004; 126: 748–752.

Vibration signal modeling of planetary gearbox based on dynamics and finite element method for tooth crack diagnosis

Abstract

Keywords

Introduction

Dynamics modeling and simulation analysis of faulty planetary gear train

Dynamic equations

Time-varying meshing stiffness evaluation of planetary gear train

Bearing stiffness evaluation

Numerical simulation and analysis

Transmission path modeling of planetary gearbox

Transmission path effect analysis of measuring point signal

Mathematical modeling of transmission path

Characteristic analysis of resultant signal

Comparison and verification

Simulation comparison

Experimental verification

Conclusion

Footnotes

Appendix

Declaration of conflicting interests

Funding

ORCID iD

References