# Root Crack Identification of Sun Gear in Planetary Gear System Combining Fault Dynamics with VMD Algorithm

Authors: Hongwei Fan; Yiqing Yang; Hongwei Ma; Xuhui Zhang; Xiang Wan; Xiangang Cao; Qinghua Mao; Chao Zhang; Qi Liu
Journal: Shock and Vibration
 (2021)

---

## Abstract

Planetary gearbox is widely used in various low-speed machines due to its large transmission ratio. However, the gears in a planetary gearbox are prone to the mechanical faults due to the complex dynamic heavy load. Vibration frequencies caused by an early tooth root crack of sun gear are usually difficult to accurately extract, so its fault diagnosis is one of the main challenges of planetary gearbox reliability. In this paper, a simplified tooth root crack model of sun gear is proposed, and then a rigid-flexible coupled dynamics model of the whole planetary gear system is constructed. By the numerical simulation, the fault frequencies caused by a tooth root crack of sun gear are obtained. A Variational Mode Decomposition (VMD) algorithm for the vibration frequency extraction is proposed. The measured vibration signals are decomposed into the sparse Intrinsic Mode Functions (IMFs) by the VMD, and then the IMFs are further analyzed by the spectral method to accurately extract the crack-induced frequency components. The experimental results show that the proposed dynamics model and VMD method are feasible; an error between the characteristic frequencies from the tested signal analysis and the theoretical calculation is less than 1%.


---

## Body

## 1. Introduction

Because of large transmission ratio and compact structure, the planetary gearbox is widely used in various low-speed heavy machines. However, in the actual operation process, due to the complex dynamic load, the planetary gearbox is prone to various kinds of mechanical failure [1]. When the gearbox is subjected to a heavy or impact load, it is very easy to cause the tooth root crack. Because the early crack has a small size, the vibration energy caused by the crack is weak, so the crack is difficult to identify in time, and further the crack can develop into a serious fault, such as tooth breaking. Therefore, it is significant to study the feature analysis technology of gear crack in a planetary gearbox for the crack diagnosis.At present, there are two kinds of tooth root crack modeling methods [2, 3], one is the parabolic crack model and the other is the penetrating crack model with simultaneous propagation along both tooth width and thickness directions. Chaari et al. [4, 5] put forward a penetrating tooth root crack model. Chen and Shao [6, 7] proposed a parabolic crack model with simultaneous propagation along both tooth width and thickness directions. Ma et al. [8] established a crack model with propagation only along the thickness direction, in which the crack depth is simulated by a parabola whose starting point is on the tooth transition circle and symmetric about the center line of tooth. Based on the tooth crack model, a whole dynamics model of planetary gear system can be built. There are three kinds of dynamics model of planetary gear system: lumped parameter model [9, 10], finite element model [11], and rigid-flexible coupled model [12,13]. The lumped parameter model of a planetary gearbox was established, and the influence of the installation error of planet carrier on the vibration response was studied by Zhai et al. [14]. Jiang et al. [15] built a lumped parameter model of planetary gear system to solve the excitation forces that generate vibrations. Shen et al. [16] were devoted to analyzing fault mechanism and dynamic characteristics of gear wear and proposed a torsional dynamics model of planetary gear system with tooth wear. Liu et al. [17] established a 3D finite element model of planetary gearbox and found that the effect of crack on the vibration response is mainly reflected in the frequency domain. Liu et al. [18] proposed a dynamics finite element model of planetary gear system, in which there is a localized fault to be considered for the outer race of the planet bearing, and Liu et al. [19] also presented a rigid-flexible coupled dynamics model for planetary gear system, both the flexible ring gear, flexible supports of the ring gear and sun gear, and a rectangular local fault in the planet bearing are formulated in the model. Duan et al. [20] carried out a system transmission error investigation of gearbox by proposing a rigid-flexible coupled dynamics model, and the modal test was launched to verify its effectiveness. Fan et al. [21] developed a rigid-flexible coupled dynamics model of the planetary gearbox considering the flexibility of ring gear and sun gear shaft based on the shell theory and Timoshenko beam theory, respectively. The lumped parameter model cannot take into account the influence of gear tooth deformation on the system, and the whole finite element model is inefficient. Therefore, the rigid-flexible coupled model is adopted in this paper.There are two kinds of common vibration signal analysis methods, including frequency-domain and time-frequency analysis. The frequency-domain signal analysis based on the Fourier transform is suitable for the periodic signals [22]. Time-frequency analysis is appropriate for the nonstationary signals. Typical time-frequency analysis methods include the short-time Fourier transform (STFT) [23], wavelet transform (WT) [24, 25], and Hilbert–Huang transform (HHT) [26, 27]. Huang et al. [28] proposed an Empirical Mode Decomposition (EMD) method based on the whole decomposition of signals, which avoid the selection of wavelet functions. Wang et al. [29] studied a signal analysis method based on EMD and extracted the effective signal components from the background noise. Yan and Liu [30] constructed an improved weak signal detection method to reduce the boundary effect of EMD. Ensemble Empirical Mode Decomposition (EEMD) method is proposed by Wu and Huang [31] to solve the mode confusion problem of EMD. Ma et al. [32] proposed a weak signal detection method by combining EMD with singular spectrum analysis. Due to the low computational efficiency of EEMD, Dragomiretskiy et al. [33] proposed a new signal decomposition method in 2014, i.e., Variational Mode Decomposition (VMD); it can avoid the mode confusion and has high efficiency, antinoise ability, and high-frequency resolution. Li and Zhu [34] applied the VMD for the mechanical fault diagnosis; the results show that the VMD is better than EEMD in signal decomposition effect. Wang et al. [35] used the VMD to realize the feature extraction of gear tooth crack in a single-stage gearbox, which verified the effectiveness of the VMD method for gear fault feature extraction.In order to reduce the modeling error and the machining difficulty of the cracked sun gear, a simplified gear tooth crack model is proposed in this paper. Combined with the crack model, a rigid-flexible coupled model of planetary gear system with tooth crack of sun gear is constructed, and the dynamics simulation is designed to find the crack-induced signal frequency characteristics. In order to identify the crack from the tested vibration signals, a VMD-based signal analysis method is proposed to extract the featured frequency of tooth crack. The structure of this paper is as follows. Section1 introduces the fault diagnosis principle based on vibration analysis and the theoretical characteristic frequency calculation. In Section 2, the crack modeling is introduced and the dynamics simulation is completed. Section 3 introduces the VMD-based fault feature extraction method, and in Section 4, the proposed model and analysis method are verified by experiments.


## 2. Frequency-Domain Diagnosis Basis for Planetary Gear System

In this paper, the crack mechanism of a planetary gear system is analyzed by dynamics simulation, and the crack type and its degree are judged by the characteristic frequencies [36].


### 2.1. Working Frequency of Planetary Gear System

A planetary gearbox consists of the sun gear, planet gear, ring gear, and planet carrier. If its basic components include two central gearsZ and a planet carrier X, the planetary gear system type is 2Z-X, as shown in Figure 1.Figure 1
Schematic diagram of 2Z-X planetary gear system.The dynamics response of a 2Z-X planetary gear system is more complex than that of a fixed-axis gear system because of the relative motion among all of the gears. Therefore, the calculation of transmission ratio and characteristic frequencies cannot directly refer to that of a fixed-axis gear system.The transmission ratio of a planetary gear system is(1)i=fsfc=ZrZs+1,where fc is the rotational frequency of planet carrier, Zr is the tooth number of ring gear, and Zs is the tooth number of sun gear.Since the rotational speed of sun gear is equal to the input speed of transmission shaft, the rotational frequency of sun gear is(2)fs=n60,where fs is the rotational frequency of sun gear and n is the input rotational speed of transmission shaft.The fundamental difference between the planetary and fixed-axis gear systems is that the planetary gear system has a planet carrier rotating around the main axis, and the planet gear installed on the planet carrier has both rotation and revolution. According to the relative motion principle, a public frequency can be added to the whole planetary gear system, which is equal to the rotational frequency of planet carrier and has an opposite direction, while the relative motion among the components of planetary gear system remains constant.(3)fs−fcfr−fc=−ZrZs,where fr is the rotational frequency of ring gear.As can be seen from equation (3), the rotational frequency of planet carrier is(4)fc=ZsZS+Zr⋅fs.Similarly, the rotational frequency of planet gear is(5)fp=fs−fc⋅ZsZP−fc,where fp is the rotational frequency of planet gear and Zp is the tooth number of planet gear.In the transmission process of a gear pair, the meshing frequencies of two gears engaged with each other are equal. Therefore, in a planetary gear system, the meshing frequency of each gear pair is equal, i.e.,(6)fm=fc⋅Zr=Zs⋅ZrZs+Zr⋅fs,where fm is the meshing frequency of planetary gear system.


### 2.2. Fault Frequency of Planetary Gear System

If there is a crack in one of the teeth of sun gear in a planetary gear system, the cracked gear meshes with all the gears in a rotational period relative to the planet carrier. Therefore, the fault frequency of sun gear is(7)fsg=fmZs⋅n,where fsg is the fault frequency of sun gear.Considering the modulation of fault frequency and fundamental frequency, the vibration signal model of a planetary gear system can be simply expressed as [37](8)xt=htc1+Acos2πfgt+ϕ⋅cos2πfmt+Bsin2πfgt+φ+θ,where A and B are the amplitude and frequency modulation intensity caused by gear fault, A, B > 0; Φ, φ, and θ are the initial phase; fm is the meshing frequency; fsg is the fault frequency; c is a dimensionless constant, usually c = 1; and h(t) describes an amplitude modulation of vibration transfer path.If the vibration sensor is installed at the top of ring gear, the distance from the fault position to the sensor changes with the rotation of sun gear when a local fault occurs.(9)ht=1−cos2πfst.The amplitude modulation of vibration signal can be written as [37](10)htat=1−cos2πfst⋅1+Acos2πfsgt+ϕ=1−cos2πfst+Acos2πfsgt+ϕ−12Acos2πfs+fsgt+ϕ+cos2πfs−fsgt−ϕ.The Fourier transform is applied to equation (10):(11)Hf∗Af=δf−δf−fs+Aδf−fsexpjϕ−12Aδf−fsg+fsexpjϕ+δf−fs−fsgexp−jϕ.As can be seen from equation (11) [37], when a crack occurs in the sun gear, the peak value in the amplitude spectrum appears in the rotational frequency fs and its multiple frequencies mfs, the fault frequency fsg and its multiple frequencies nfsg, and the combined frequencies fm ± fs and nfsg ± mfm. If mfs and nfsg are considered as the modulation frequencies, the other peaks mainly appear at 1/4fs and n/4fs, fm and mfm, and the combined frequencies n/4fsg ± mfm(m, n = 1, 2, ...).


## 2.1. Working Frequency of Planetary Gear System

A planetary gearbox consists of the sun gear, planet gear, ring gear, and planet carrier. If its basic components include two central gearsZ and a planet carrier X, the planetary gear system type is 2Z-X, as shown in Figure 1.Figure 1
Schematic diagram of 2Z-X planetary gear system.The dynamics response of a 2Z-X planetary gear system is more complex than that of a fixed-axis gear system because of the relative motion among all of the gears. Therefore, the calculation of transmission ratio and characteristic frequencies cannot directly refer to that of a fixed-axis gear system.The transmission ratio of a planetary gear system is(1)i=fsfc=ZrZs+1,where fc is the rotational frequency of planet carrier, Zr is the tooth number of ring gear, and Zs is the tooth number of sun gear.Since the rotational speed of sun gear is equal to the input speed of transmission shaft, the rotational frequency of sun gear is(2)fs=n60,where fs is the rotational frequency of sun gear and n is the input rotational speed of transmission shaft.The fundamental difference between the planetary and fixed-axis gear systems is that the planetary gear system has a planet carrier rotating around the main axis, and the planet gear installed on the planet carrier has both rotation and revolution. According to the relative motion principle, a public frequency can be added to the whole planetary gear system, which is equal to the rotational frequency of planet carrier and has an opposite direction, while the relative motion among the components of planetary gear system remains constant.(3)fs−fcfr−fc=−ZrZs,where fr is the rotational frequency of ring gear.As can be seen from equation (3), the rotational frequency of planet carrier is(4)fc=ZsZS+Zr⋅fs.Similarly, the rotational frequency of planet gear is(5)fp=fs−fc⋅ZsZP−fc,where fp is the rotational frequency of planet gear and Zp is the tooth number of planet gear.In the transmission process of a gear pair, the meshing frequencies of two gears engaged with each other are equal. Therefore, in a planetary gear system, the meshing frequency of each gear pair is equal, i.e.,(6)fm=fc⋅Zr=Zs⋅ZrZs+Zr⋅fs,where fm is the meshing frequency of planetary gear system.


## 2.2. Fault Frequency of Planetary Gear System

If there is a crack in one of the teeth of sun gear in a planetary gear system, the cracked gear meshes with all the gears in a rotational period relative to the planet carrier. Therefore, the fault frequency of sun gear is(7)fsg=fmZs⋅n,where fsg is the fault frequency of sun gear.Considering the modulation of fault frequency and fundamental frequency, the vibration signal model of a planetary gear system can be simply expressed as [37](8)xt=htc1+Acos2πfgt+ϕ⋅cos2πfmt+Bsin2πfgt+φ+θ,where A and B are the amplitude and frequency modulation intensity caused by gear fault, A, B > 0; Φ, φ, and θ are the initial phase; fm is the meshing frequency; fsg is the fault frequency; c is a dimensionless constant, usually c = 1; and h(t) describes an amplitude modulation of vibration transfer path.If the vibration sensor is installed at the top of ring gear, the distance from the fault position to the sensor changes with the rotation of sun gear when a local fault occurs.(9)ht=1−cos2πfst.The amplitude modulation of vibration signal can be written as [37](10)htat=1−cos2πfst⋅1+Acos2πfsgt+ϕ=1−cos2πfst+Acos2πfsgt+ϕ−12Acos2πfs+fsgt+ϕ+cos2πfs−fsgt−ϕ.The Fourier transform is applied to equation (10):(11)Hf∗Af=δf−δf−fs+Aδf−fsexpjϕ−12Aδf−fsg+fsexpjϕ+δf−fs−fsgexp−jϕ.As can be seen from equation (11) [37], when a crack occurs in the sun gear, the peak value in the amplitude spectrum appears in the rotational frequency fs and its multiple frequencies mfs, the fault frequency fsg and its multiple frequencies nfsg, and the combined frequencies fm ± fs and nfsg ± mfm. If mfs and nfsg are considered as the modulation frequencies, the other peaks mainly appear at 1/4fs and n/4fs, fm and mfm, and the combined frequencies n/4fsg ± mfm(m, n = 1, 2, ...).


## 3. Crack Modeling and Dynamics Analysis

### 3.1. Geometric Model of Planetary Gear System

The 2Z-X planetary gear system used is composed of a sun gear, a planet carrier, a ring gear, and four planet gears. A 3D model of the planetary gear system is shown in Figure 2; the design parameters are shown in Table 1. On the premise of not affecting the analysis results, the auxiliary parts are simplified; for example, both the bearings and bolts used to connect the planet gears and the planet carrier are treated as cylinders.Figure 2
3D model of planetary gear system.Table 1
Design parameters of planetary gear system.
Types of gearsRing gearPlanet gearSun gearNumber of teeth1003628Module (mm)1Pressure angle (°)20Tooth width (mm)10MaterialSteelYoung’s modulus (Pa×1011)2.07Siméon Denis Poisson0.29Input speed (rpm)1800


### 3.2. Rigid-Flexible Coupled Model with Crack

In the above planetary gear system, the sun gear contacts with four planet gears to transmit motion and torque. When the sun gear is used as an input, it is simultaneously subjected to forces and torques from four planet gears. Therefore, the sun gear is prone to a mechanical fault. During gear transmission, the bending stress is concentrated at the tooth root, and the crack begins and extends from the tooth root [38], as shown in Figure 3.Figure 3
Tooth root crack model of sun gear.In Figure3, φ is the included angle between the tooth center line and the tangent line of root fillet curve, and the area between φ = 30° and φ = 34° is the stress concentration area [6]. The included angle between the crack extension line and tooth center line is defined as the crack extension angle α. In order to approach the crack growth law under the real working condition, φ = 32° and α = 74° are taken in the following dynamics simulation. The chord of the approximate arc of the tooth root crack growth path is taken as the crack width growth line, i.e., w, shown in Figure 3, and the tooth width direction is called the crack length extension line, i.e., l, shown in Figure 3.In order to study the frequency characteristics of different-scale cracks, three kinds of crack were designed, as shown in Table2. It is known that the tooth width of sun gear is l = 10 mm, and the maximum crack width is w = 1 mm. The crack width and length were expanded simultaneously in a proportion of 20%, 40%, and 80% relative to the chord of tooth root arc and tooth width, respectively.Table 2
Crack parameters of sun gear.
Crack caseNo crack1# crack2# crack3# crackCrack length,l (mm)0248Crack width,w (mm)00.20.40.8In order to reduce the difficulty of making crack, the crack width propagation path is treated approximately as an arc chord, as shown in Figure4.Figure 4
3D crack model of sun gear. (a) 1# crack. (b) 2# crack. (c) 3# crack.
(a)(b)(c)In order to obtain the crack-induced vibration frequency characteristics of sun gear, the finite element model is used to treat the sun gear with flexibility. The flexibilization steps are shown in Figure5.Figure 5
Flexibilization steps of sun gear.As shown in Figure5, for a finite element model, the generation quality of mesh is closely related to the element type and quantity. In this paper, ANSYS is used to make sun gear flexible. Considering the crack shape, a tetrahedron element is suitable for the gear with crack. For the tetrahedral element, the SOLID187 uses an optimized algorithm compared to the SOLID92 to achieve a high-quality model. When generating the mesh, the part without crack is firstly divided, and then the crack area is refined separately. The meshed sun gear is shown in Figure 6. Figure 6(a) is the entire mesh model; Figure 6(b) is the crack mesh model.Figure 6
Mesh model of sun gear. (a) Entire mesh. (b) Crack mesh.
(a)(b)When a finite element model is introduced into a rigid-body dynamics system and connected to other parts, there is a region that does not deform, i.e., the rigid region. In this paper, ADAMS is used for the dynamics simulation; the rigid region is composed of the nodes of two interfaces and the surrounding nodes. Select those nodes as the master nodes and the nodes of end circles of sun gear as the slave nodes; the rigid region is created, as shown in Figure7(a); the rigid-flexible coupled model of planetary gear system is obtained by connecting the flexible sun gear with other rigid bodies through rigid regions and interface nodes, as shown in Figure 7(b).Figure 7
Dynamics model of planetary gear system. (a) Rigid region model. (b) Rigid-flexible coupled model.
(a)(b)


### 3.3. Dynamics Analysis

In dynamics simulation, the input shaft of torque and speed is the shaft where the sun gear is located, and the output shaft is that connected by the planet carrier. The parameters are set as shown in Table3.Table 3
Parameter setup used for dynamics simulation.
ParameterSetupInput frequency (Hz)30Input torque (N·m)5Function typeSTEPThe input rotational frequency is 30 Hz, the torque is 5 N·m, and the direction of torque is opposite to the direction of rotational frequency. In order to avoid the sudden loading, the STEP function was adopted, which is expressed as STEP (time,0,0,0.1,10800D), the initial speed is 0, the speed is 10800°/s at 0.1 s, i.e., the frequency is 30 Hz, and the torque gradually increased to 5 N·m in the form of STEP (time,0,0,0.1,5000). In the model settings, the rigid-rigid contact was used between four planet gears and the ring gear, the rigid-flexible contact was used between four planet gears and the sun gear, and the other parameters were considered in terms of [39].In the simulation, the ring gear was fixed and the contact force between the gears was added to simulate the meshing. The contact force between the gears changed the angular acceleration with the meshing process. Therefore, the characteristic frequencies of the tooth crack were analyzed by the angular acceleration signal of sun gear.The dynamics simulation was carried out and the results of angular acceleration in the time and frequency domains were obtained under different crack cases, as shown in Figure8. In Figures 8(a) and 8(b), four subfigures from top to bottom are the angular acceleration results in the time and frequency domains under no crack, 1# crack, 2# crack, and 3# crack.Figure 8
Angular acceleration of sun gear for different crack cases. (a) Time-domain waveform. (b) Frequency-domain spectrum.
(a)(b)As can be seen from Figure8(a), the amplitude of the angular acceleration is obviously smaller when the sun gear has no crack than a crack fault, and the amplitude increases with the increase of the crack size. As can be seen from Figure 8(b), the amplitude peaks in the frequency spectrum mainly appear in the meshing frequency fm and its multiple frequencies mfm, the fault frequencies fsg of sun gear and its multiple frequencies nfsg, and the combined frequencies mfm + nfsg.


## 3.1. Geometric Model of Planetary Gear System

The 2Z-X planetary gear system used is composed of a sun gear, a planet carrier, a ring gear, and four planet gears. A 3D model of the planetary gear system is shown in Figure 2; the design parameters are shown in Table 1. On the premise of not affecting the analysis results, the auxiliary parts are simplified; for example, both the bearings and bolts used to connect the planet gears and the planet carrier are treated as cylinders.Figure 2
3D model of planetary gear system.Table 1
Design parameters of planetary gear system.
Types of gearsRing gearPlanet gearSun gearNumber of teeth1003628Module (mm)1Pressure angle (°)20Tooth width (mm)10MaterialSteelYoung’s modulus (Pa×1011)2.07Siméon Denis Poisson0.29Input speed (rpm)1800


## 3.2. Rigid-Flexible Coupled Model with Crack

In the above planetary gear system, the sun gear contacts with four planet gears to transmit motion and torque. When the sun gear is used as an input, it is simultaneously subjected to forces and torques from four planet gears. Therefore, the sun gear is prone to a mechanical fault. During gear transmission, the bending stress is concentrated at the tooth root, and the crack begins and extends from the tooth root [38], as shown in Figure 3.Figure 3
Tooth root crack model of sun gear.In Figure3, φ is the included angle between the tooth center line and the tangent line of root fillet curve, and the area between φ = 30° and φ = 34° is the stress concentration area [6]. The included angle between the crack extension line and tooth center line is defined as the crack extension angle α. In order to approach the crack growth law under the real working condition, φ = 32° and α = 74° are taken in the following dynamics simulation. The chord of the approximate arc of the tooth root crack growth path is taken as the crack width growth line, i.e., w, shown in Figure 3, and the tooth width direction is called the crack length extension line, i.e., l, shown in Figure 3.In order to study the frequency characteristics of different-scale cracks, three kinds of crack were designed, as shown in Table2. It is known that the tooth width of sun gear is l = 10 mm, and the maximum crack width is w = 1 mm. The crack width and length were expanded simultaneously in a proportion of 20%, 40%, and 80% relative to the chord of tooth root arc and tooth width, respectively.Table 2
Crack parameters of sun gear.
Crack caseNo crack1# crack2# crack3# crackCrack length,l (mm)0248Crack width,w (mm)00.20.40.8In order to reduce the difficulty of making crack, the crack width propagation path is treated approximately as an arc chord, as shown in Figure4.Figure 4
3D crack model of sun gear. (a) 1# crack. (b) 2# crack. (c) 3# crack.
(a)(b)(c)In order to obtain the crack-induced vibration frequency characteristics of sun gear, the finite element model is used to treat the sun gear with flexibility. The flexibilization steps are shown in Figure5.Figure 5
Flexibilization steps of sun gear.As shown in Figure5, for a finite element model, the generation quality of mesh is closely related to the element type and quantity. In this paper, ANSYS is used to make sun gear flexible. Considering the crack shape, a tetrahedron element is suitable for the gear with crack. For the tetrahedral element, the SOLID187 uses an optimized algorithm compared to the SOLID92 to achieve a high-quality model. When generating the mesh, the part without crack is firstly divided, and then the crack area is refined separately. The meshed sun gear is shown in Figure 6. Figure 6(a) is the entire mesh model; Figure 6(b) is the crack mesh model.Figure 6
Mesh model of sun gear. (a) Entire mesh. (b) Crack mesh.
(a)(b)When a finite element model is introduced into a rigid-body dynamics system and connected to other parts, there is a region that does not deform, i.e., the rigid region. In this paper, ADAMS is used for the dynamics simulation; the rigid region is composed of the nodes of two interfaces and the surrounding nodes. Select those nodes as the master nodes and the nodes of end circles of sun gear as the slave nodes; the rigid region is created, as shown in Figure7(a); the rigid-flexible coupled model of planetary gear system is obtained by connecting the flexible sun gear with other rigid bodies through rigid regions and interface nodes, as shown in Figure 7(b).Figure 7
Dynamics model of planetary gear system. (a) Rigid region model. (b) Rigid-flexible coupled model.
(a)(b)


## 3.3. Dynamics Analysis

In dynamics simulation, the input shaft of torque and speed is the shaft where the sun gear is located, and the output shaft is that connected by the planet carrier. The parameters are set as shown in Table3.Table 3
Parameter setup used for dynamics simulation.
ParameterSetupInput frequency (Hz)30Input torque (N·m)5Function typeSTEPThe input rotational frequency is 30 Hz, the torque is 5 N·m, and the direction of torque is opposite to the direction of rotational frequency. In order to avoid the sudden loading, the STEP function was adopted, which is expressed as STEP (time,0,0,0.1,10800D), the initial speed is 0, the speed is 10800°/s at 0.1 s, i.e., the frequency is 30 Hz, and the torque gradually increased to 5 N·m in the form of STEP (time,0,0,0.1,5000). In the model settings, the rigid-rigid contact was used between four planet gears and the ring gear, the rigid-flexible contact was used between four planet gears and the sun gear, and the other parameters were considered in terms of [39].In the simulation, the ring gear was fixed and the contact force between the gears was added to simulate the meshing. The contact force between the gears changed the angular acceleration with the meshing process. Therefore, the characteristic frequencies of the tooth crack were analyzed by the angular acceleration signal of sun gear.The dynamics simulation was carried out and the results of angular acceleration in the time and frequency domains were obtained under different crack cases, as shown in Figure8. In Figures 8(a) and 8(b), four subfigures from top to bottom are the angular acceleration results in the time and frequency domains under no crack, 1# crack, 2# crack, and 3# crack.Figure 8
Angular acceleration of sun gear for different crack cases. (a) Time-domain waveform. (b) Frequency-domain spectrum.
(a)(b)As can be seen from Figure8(a), the amplitude of the angular acceleration is obviously smaller when the sun gear has no crack than a crack fault, and the amplitude increases with the increase of the crack size. As can be seen from Figure 8(b), the amplitude peaks in the frequency spectrum mainly appear in the meshing frequency fm and its multiple frequencies mfm, the fault frequencies fsg of sun gear and its multiple frequencies nfsg, and the combined frequencies mfm + nfsg.


## 4. VMD Time-Frequency Analysis Method of Vibration Signals

The frequency resolution of the spectrum in Figure8(b) is low, so it is difficult to accurately find the fault frequency values. Therefore, the VMD algorithm [40] is proposed for the time-frequency analysis of complex vibration signals.


### 4.1. Principle of VMD Algorithm

VMD is the decomposition of signalx(t) into a number of sparse IMFs mk(t), each of which fluctuates around the central frequency ωk, i.e., the bandwidth is limited around the central frequency. The bandwidth can be obtained by estimating the square norm of the gradient of the frequency shift signals.First, the HHT transforms the real modemk into an analytic signal mk+ with a unilateral spectrum with nonnegative frequencies:(12)mk+t=δt+jπt∗mkt.Then, by multiplying with an exponential harmonic wave whose frequencies are adjusted to their respective central frequencies, the spectrum is shifted to the base band:(13)mkBt=δt+jπt∗mkte−jωkt.Finally, the bandwidth is estimated by the square norm of the gradient:(14)△ω=∂tδt+jπt∗mkte−jωkt22.Thus, all the estimated mode component bandwidths can be expressed as a constrained variational model:(15)minmkt,ωk=∑k=1K∂tδt+jπt∗mkte−jωkt22,s.t.∑k=1Kmkt=xt,where δ(·) is the Dirac Delta function, ∗ is the convolution operation, and K is the number of the IMFs.In order to transform the optimization problem expressed by the above formula into an unconstrained form, a penalty factorα and a Lagrange multiplier are added to accelerate the convergence and strengthen the constraint. The objective function that needs to be minimized is transformed into an augmented Lagrange function:(16)Lmkt,ωk,λt=α∑k=1K∂tδt+jπt∗mkte−jωkt22+xt−∑k=1Kmkt22+λt,xt−∑k=1Kmkt,where λ(t) is the Lagrange multiplier and <·,·> is the inner product operation.The minimization problem represented by equation (15) is solved by an Alternate Direction Multiplier (ADM) method to find the saddle point of an augmented Lagrange problem represented by equation (16). The equivalent minimization problem of equation (16) is solved, and each IMF is updated by the optimal solution mk(t):(17)mkt=argminmkLmkt,ωk,λt=argminmkα∑k=1K∂tδt+jπt∗mkte−jωkt22+xt−∑k=1Kmkt+λt222.In the frequency domain, the solution of equation (17) is(18)m^kω=x^ω−∑i≠km^iω+1/2λ^ω1+2αω−ωk2.The IMF in the time domain can be obtained by Fourier inverse transform of equation (18).The equivalent minimization problem of equation (16) is solved, and the central frequency ωk of each IMF mk(t) is updated with the optimal solution:(19)ωk=argminωkLmkt,ωk,λt=argminωk∑k=1K∂tδt+jπt∗mkte−jωkt22.The central frequencyωk is the centroid of the IMF power spectrum m^kω , i.e.,(20)ωk=∫0∞ωm^kω2dω∫0∞m^kω2dω.According to the solution process of the above constrained variational problem, the VMD can divide the frequency band according to the frequency characteristic of the signal and realize the adaptive decomposition of the signal by updating each mode and its corresponding central frequency.The flowchart of the VMD algorithm is shown in Figure9.Step 1. Set m^k0tω^k0λ^0t, and n = 0; define the convergence threshold and number K of IMFs to be separated.Step 2. k=1,2,...,K and ω≥0; update each IMF mkt and its central frequency ωk:(21)m^kn+1ω=x^ω−∑i<km^in+1ω−∑i>km^inω+1/2λ^nω1+2αω−ωkn2,ωkn+1=∫0∞ωm^kn+1ω2dω∫0∞m^kn+1ω2dω.Step 3. ω≥0; update the Lagrange multiplier:(22)λ^n+1ω=λ^nω+τx^ω−∑km^kn+1ω,where τ is the updating parameter of Lagrange multiplier.Figure 9
Flowchart of VMD algorithm.Check the convergence condition:(23)∑km^kn+1t−m^knt22m^knt22<ε.If the formula is true, letmkt=m^kn+1t, ωkt=ωkn+1, terminate the decomposition, and get K IMFs with the finite bandwidth. Otherwise, let n=n+1; return to Step 2.


### 4.2. Frequency-Domain Feature Extraction Using VMD

In the VMD algorithm, the value ofK determines the number of IMFs in the decomposed signal. Taking 3# crack as an example, when K takes different values, the signal is decomposed by the VMD, and the frequency of each IMF is shown in Table 4. From Table 4, the central frequency of IMF6 is 4439 Hz, which is close to the maximum value of 4499 Hz in the case of K = 7, and the central frequency of IMF4 is close to that of IMF5 in the case of K = 7. In addition, since the planetary gear system consists of six moving parts, K = 6 is reasonable.Table 4
Central frequency of each IMF for 3# crack.
KCentral frequency of each IMF (Hz)IMF1IMF2IMF3IMF4IMF5IMF6IMF7211822735—————3118219543457————4390.8112227253767———5390.81182195428464439——6390.811821954273534574439—7390.8118219542475284637674499WhenK = 6, the angular acceleration signal for 3# crack is decomposed by the VMD, as shown in Figure 10, the VMD decomposition results for other crack cases are shown in Figure 11.Figure 10
Each IMF of angular acceleration signal for 3# crack. (a) Time-domain waveform. (b) Frequency-domain spectrum.
(a)(b)Figure 11
Each IMF of angular acceleration signal for no crack and 1# to 3# crack. (a) No crack. (b) 1# crack. (c) 2# crack. (d) 3# crack.
(a)(b)(c)(d)Figure10 shows that the IMF frequencies correspond one to one with the central frequencies in Table 4 and are all the characteristic frequencies for sun gear crack. The central frequencies and amplitudes of the IMFs in Figure 11 are listed in Table 5. According to Table 5, Figure 12 is obtained.Table 5
Central frequency and amplitude of each IMF.
IMFNo crack1# crack2# crack3# crackf (Hz)A (105)f (Hz)A (105)f (Hz)A (105)f (Hz)A (105)IMF1360.71.8665117.536400.87.502390.89.014IMF211822.42415035.21611926.678118210.280IMF317942.78520146.29816834.23419547.243IMF425553.73125854.34824655.81727356.026IMF531262.82936773.44734874.37534576.187IMF641683.52946793.16645294.18444395.868Total—17.164—30.011—32.79—44.618Figure 12
Central frequency and amplitude of each IMF.As can be seen from Figure12 and Table 5 above, the amplitude shown by the black line marked by rhombic dot is lower than others. When the cracks appear in the sun gear, the amplitudes corresponding to the central frequencies of each IMF are higher than those with no crack, and the errors of the central frequencies of each IMF are small.In the case of 3# crack, the errors between the IMF central frequencies and their theoretical values are shown in Table 6.Table 6
Error between the simulated VMD IMF central frequency and theoretical value for 3# crack.
IMFVMD frequency (Hz)Theoretical frequency (Hz)Corresponding featureError (%)IMF1390.8398.2fm−11/4fsg1.86IMF211821194.82fm−5fsg1.07IMF319541944.63fm−1/4fsg0.48IMF427352741.24fm + 5/4fsg0.23IMF534573467.84fm + 9fsg0.31IMF644394427.97fm−7/4fsg0.25As shown in Figure13 and Table 6, the frequencies of each IMF correspond to the rotational frequency fs of sun gear, the meshing frequency fm, the fault frequencies fsg and n/4fsg of sun gear, and the combined results of the above frequencies, respectively. The errors between the characteristic frequencies from dynamics simulation and theoretical calculation are less than 2%. It can be seen that the central frequencies of IMF components obtained by VMD can be used as the featured parameters identifying cracks in a planetary gear system.Figure 13
Error between the simulated VMD IMF central frequency and theoretical value for 3# crack.


## 4.1. Principle of VMD Algorithm

VMD is the decomposition of signalx(t) into a number of sparse IMFs mk(t), each of which fluctuates around the central frequency ωk, i.e., the bandwidth is limited around the central frequency. The bandwidth can be obtained by estimating the square norm of the gradient of the frequency shift signals.First, the HHT transforms the real modemk into an analytic signal mk+ with a unilateral spectrum with nonnegative frequencies:(12)mk+t=δt+jπt∗mkt.Then, by multiplying with an exponential harmonic wave whose frequencies are adjusted to their respective central frequencies, the spectrum is shifted to the base band:(13)mkBt=δt+jπt∗mkte−jωkt.Finally, the bandwidth is estimated by the square norm of the gradient:(14)△ω=∂tδt+jπt∗mkte−jωkt22.Thus, all the estimated mode component bandwidths can be expressed as a constrained variational model:(15)minmkt,ωk=∑k=1K∂tδt+jπt∗mkte−jωkt22,s.t.∑k=1Kmkt=xt,where δ(·) is the Dirac Delta function, ∗ is the convolution operation, and K is the number of the IMFs.In order to transform the optimization problem expressed by the above formula into an unconstrained form, a penalty factorα and a Lagrange multiplier are added to accelerate the convergence and strengthen the constraint. The objective function that needs to be minimized is transformed into an augmented Lagrange function:(16)Lmkt,ωk,λt=α∑k=1K∂tδt+jπt∗mkte−jωkt22+xt−∑k=1Kmkt22+λt,xt−∑k=1Kmkt,where λ(t) is the Lagrange multiplier and <·,·> is the inner product operation.The minimization problem represented by equation (15) is solved by an Alternate Direction Multiplier (ADM) method to find the saddle point of an augmented Lagrange problem represented by equation (16). The equivalent minimization problem of equation (16) is solved, and each IMF is updated by the optimal solution mk(t):(17)mkt=argminmkLmkt,ωk,λt=argminmkα∑k=1K∂tδt+jπt∗mkte−jωkt22+xt−∑k=1Kmkt+λt222.In the frequency domain, the solution of equation (17) is(18)m^kω=x^ω−∑i≠km^iω+1/2λ^ω1+2αω−ωk2.The IMF in the time domain can be obtained by Fourier inverse transform of equation (18).The equivalent minimization problem of equation (16) is solved, and the central frequency ωk of each IMF mk(t) is updated with the optimal solution:(19)ωk=argminωkLmkt,ωk,λt=argminωk∑k=1K∂tδt+jπt∗mkte−jωkt22.The central frequencyωk is the centroid of the IMF power spectrum m^kω , i.e.,(20)ωk=∫0∞ωm^kω2dω∫0∞m^kω2dω.According to the solution process of the above constrained variational problem, the VMD can divide the frequency band according to the frequency characteristic of the signal and realize the adaptive decomposition of the signal by updating each mode and its corresponding central frequency.The flowchart of the VMD algorithm is shown in Figure9.Step 1. Set m^k0tω^k0λ^0t, and n = 0; define the convergence threshold and number K of IMFs to be separated.Step 2. k=1,2,...,K and ω≥0; update each IMF mkt and its central frequency ωk:(21)m^kn+1ω=x^ω−∑i<km^in+1ω−∑i>km^inω+1/2λ^nω1+2αω−ωkn2,ωkn+1=∫0∞ωm^kn+1ω2dω∫0∞m^kn+1ω2dω.Step 3. ω≥0; update the Lagrange multiplier:(22)λ^n+1ω=λ^nω+τx^ω−∑km^kn+1ω,where τ is the updating parameter of Lagrange multiplier.Figure 9
Flowchart of VMD algorithm.Check the convergence condition:(23)∑km^kn+1t−m^knt22m^knt22<ε.If the formula is true, letmkt=m^kn+1t, ωkt=ωkn+1, terminate the decomposition, and get K IMFs with the finite bandwidth. Otherwise, let n=n+1; return to Step 2.


## 4.2. Frequency-Domain Feature Extraction Using VMD

In the VMD algorithm, the value ofK determines the number of IMFs in the decomposed signal. Taking 3# crack as an example, when K takes different values, the signal is decomposed by the VMD, and the frequency of each IMF is shown in Table 4. From Table 4, the central frequency of IMF6 is 4439 Hz, which is close to the maximum value of 4499 Hz in the case of K = 7, and the central frequency of IMF4 is close to that of IMF5 in the case of K = 7. In addition, since the planetary gear system consists of six moving parts, K = 6 is reasonable.Table 4
Central frequency of each IMF for 3# crack.
KCentral frequency of each IMF (Hz)IMF1IMF2IMF3IMF4IMF5IMF6IMF7211822735—————3118219543457————4390.8112227253767———5390.81182195428464439——6390.811821954273534574439—7390.8118219542475284637674499WhenK = 6, the angular acceleration signal for 3# crack is decomposed by the VMD, as shown in Figure 10, the VMD decomposition results for other crack cases are shown in Figure 11.Figure 10
Each IMF of angular acceleration signal for 3# crack. (a) Time-domain waveform. (b) Frequency-domain spectrum.
(a)(b)Figure 11
Each IMF of angular acceleration signal for no crack and 1# to 3# crack. (a) No crack. (b) 1# crack. (c) 2# crack. (d) 3# crack.
(a)(b)(c)(d)Figure10 shows that the IMF frequencies correspond one to one with the central frequencies in Table 4 and are all the characteristic frequencies for sun gear crack. The central frequencies and amplitudes of the IMFs in Figure 11 are listed in Table 5. According to Table 5, Figure 12 is obtained.Table 5
Central frequency and amplitude of each IMF.
IMFNo crack1# crack2# crack3# crackf (Hz)A (105)f (Hz)A (105)f (Hz)A (105)f (Hz)A (105)IMF1360.71.8665117.536400.87.502390.89.014IMF211822.42415035.21611926.678118210.280IMF317942.78520146.29816834.23419547.243IMF425553.73125854.34824655.81727356.026IMF531262.82936773.44734874.37534576.187IMF641683.52946793.16645294.18444395.868Total—17.164—30.011—32.79—44.618Figure 12
Central frequency and amplitude of each IMF.As can be seen from Figure12 and Table 5 above, the amplitude shown by the black line marked by rhombic dot is lower than others. When the cracks appear in the sun gear, the amplitudes corresponding to the central frequencies of each IMF are higher than those with no crack, and the errors of the central frequencies of each IMF are small.In the case of 3# crack, the errors between the IMF central frequencies and their theoretical values are shown in Table 6.Table 6
Error between the simulated VMD IMF central frequency and theoretical value for 3# crack.
IMFVMD frequency (Hz)Theoretical frequency (Hz)Corresponding featureError (%)IMF1390.8398.2fm−11/4fsg1.86IMF211821194.82fm−5fsg1.07IMF319541944.63fm−1/4fsg0.48IMF427352741.24fm + 5/4fsg0.23IMF534573467.84fm + 9fsg0.31IMF644394427.97fm−7/4fsg0.25As shown in Figure13 and Table 6, the frequencies of each IMF correspond to the rotational frequency fs of sun gear, the meshing frequency fm, the fault frequencies fsg and n/4fsg of sun gear, and the combined results of the above frequencies, respectively. The errors between the characteristic frequencies from dynamics simulation and theoretical calculation are less than 2%. It can be seen that the central frequencies of IMF components obtained by VMD can be used as the featured parameters identifying cracks in a planetary gear system.Figure 13
Error between the simulated VMD IMF central frequency and theoretical value for 3# crack.


## 5. Experimental Study

### 5.1. Test Setup

The test platform of a planetary gear system was established, as shown in Figure14.Figure 14
Test platform.The test platform is composed of motor, frequency converter, planetary gearbox (one-stage transmission, deceleration ratio is 4.57 : 1), parallel-axle spur gearbox (two-stage transmission, deceleration ratio is 8.63 : 1), magnetic powder brake, power supply (adjustable current source), vibration sensor, signal collector, and upper computer.For easy machining, three kinds of sun gear crack were made by Electrical Discharge Machining (EDM), and the actual “rectangle” is used instead of the ideal “trilateral” cracks shown in Figure4 due to the limit of fabrication process, which does not affect the characteristic frequency extraction. Crack parameters are shown in Table 7.Table 7
Actual crack parameters of tested gear.
Crack number1# crack2# crack3# crackCrack length (mm)1.433.928.47Crack width (mm)0.540.691.34Crack height (mm)0.50Set the motor speed to 1800 rpm and the torque to 5 N·m. The 608A11-type accelerometer with a sensitivity of 0.103 V/g was used to measure the vertical vibration. The Spectral Quest data acquisition instrument was used, its sampling frequency is 10.24 kHz, and the sampled signal length is 6 s to 7 s.


## 5.1. Test Setup

The test platform of a planetary gear system was established, as shown in Figure14.Figure 14
Test platform.The test platform is composed of motor, frequency converter, planetary gearbox (one-stage transmission, deceleration ratio is 4.57 : 1), parallel-axle spur gearbox (two-stage transmission, deceleration ratio is 8.63 : 1), magnetic powder brake, power supply (adjustable current source), vibration sensor, signal collector, and upper computer.For easy machining, three kinds of sun gear crack were made by Electrical Discharge Machining (EDM), and the actual “rectangle” is used instead of the ideal “trilateral” cracks shown in Figure4 due to the limit of fabrication process, which does not affect the characteristic frequency extraction. Crack parameters are shown in Table 7.Table 7
Actual crack parameters of tested gear.
Crack number1# crack2# crack3# crackCrack length (mm)1.433.928.47Crack width (mm)0.540.691.34Crack height (mm)0.50Set the motor speed to 1800 rpm and the torque to 5 N·m. The 608A11-type accelerometer with a sensitivity of 0.103 V/g was used to measure the vertical vibration. The Spectral Quest data acquisition instrument was used, its sampling frequency is 10.24 kHz, and the sampled signal length is 6 s to 7 s.


## 6. Results and Discussion

The vibration signals of sun gear under three kinds of actual crack were measured. The time-domain waveform and frequency-domain spectrum are shown in Figure15.Figure 15
Time-domain waveform and frequency-domain spectrum for no crack and three crack cases.The results of no crack and 1# to 3# cracks are listed from top to bottom. It can be seen that the amplitude of time-domain waveform increases with the increase of crack size. The peak value appears in the meshing frequency fm and its multiple frequencies mfm, local fault frequency fsg of sun gear and its multiple frequencies nfsg, 1/4fsg and its multiple frequencies n/4fsg, and combined frequencies mfm + nfsg and mfm + n/fsg.The VMD results of vibration signals under 1# to 3# cracks at K = 6 are shown in Figure 16; the central frequencies and amplitudes of each IMF are shown in Table 8.Figure 16
VMD results for no crack and three crack cases. (a) No crack. (b) 1# crack. (c) 2# crack. (d) 3# crack.
(a)(b)(c)(d)Table 8
Central frequency and amplitude of each IMF for 3# crack.
No crack1# crack2# crack3# crackf (Hz)A (10−3)f (Hz)A (10−3)f (Hz)A (10−3)f (Hz)A (10−3)IMF130.062.0730.062.5130.062.3530.062.96IMF2651.32.91641.32.32661.31.70651.32.54IMF319641.0219440.8519740.9719541.85IMF428261.0025851.4726152.0726152.22IMF532861.0832162.6132365.1732365.42IMF645090.8543691.1044991.3945592.06Total—8.93—10.86—13.65—17.05From Figure17, it can be seen that with the increase of crack size, the total vibration energy increases, which proves that this experimental design is correct. However, due to the sun gear change and the complex gap distribution among the gear elements in an actual planetary gear system, these nonlinear factors inevitably affect the value of each IMF. The amplitudes corresponding to all the IMF central frequencies except the IMF2 of 1# to 3# cracks and the IMF3 of 1# to 2# cracks are larger than those with no crack.Figure 17
Central frequency and amplitude of each IMF.In the case of 3# crack, the errors between the VMD and the theoretical frequency values of each IMF are shown in Table 9 and Figure 18. The frequencies of each IMF correspond to the rotational frequency fs of sun gear, the meshing frequency fm, the fault frequency of sun gear fsg and n/4fsg, and the combined results of the above frequencies, respectively. It can be seen that the errors between the VMD and the theoretical values are less than 1%.Table 9
Error between the tested IMF central frequency and theoretical value for 3# crack.
IMFTheoretical frequency (Hz)VMD frequency (Hz)Corresponding featureError (%)IMF13030.06fs0.20IMF2656651.3fm0.72IMF31944.619543fm−1/4fsg0.48IMF42600.626154fm−1/4fsg0.55IMF53233.132365fm−2/4fsg0.09IMF64568.645597fm−1/4fsg0.21Figure 18
Error between the tested VMD IMF central frequency and theoretical value for 3# crack.


## 7. Conclusions

This paper proposes a tooth root crack identification method of sun gear in a planetary gear system combining the fault dynamics with VMD algorithm, which provides a positive contribution to the fault diagnosis of a complex mechanical transmission system.(1)
In order to overcome the shortcomings of the existing tooth crack models, a gear crack model with the simultaneous propagation in both length and width directions is proposed. And a rigid-flexible coupled model of planetary gear system is established by treating the sun gear as a flexible body and connecting it with other rigid components. The proposed model is appropriate for the fault frequency determination considering the computational accuracy and efficiency.(2)
In order to overcome the shortcomings of the existing EMD algorithms, the VMD method for the nonstationary weak signal is used. Through the analysis of simulated and tested signals, the total amplitude corresponding to each IMF central frequency obtained by VMD can reflect the level change of different crack scales, and the error between the tested VMD frequency and the theoretical value is less than 1%, which proves the proposed VMD algorithm effective.(3)
The test platform of a planetary gear system was established, and the tooth root crack of sun gear was made. The tested results show that the meshing frequency, fault frequency, and their combined frequency are prominent in the amplitude spectrum of vibration signal with crack, so the gear crack model and rigid-flexible coupled system model established in this paper are feasible and the VMD method is also reliable, which provides an effective approach for the crack diagnosis of sun gear in a planetary gear system.



---

*Source: 5561417-2021-05-17.xml*